C# QR Decomp Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the QR decomposition classes.
  /// </summary>
  class QRDecompExample
  {

    static void Main( string[] args )
    {
      // A general m x n system with random entries.
      var rng = new RandGenUniform( -1, 1 );
      rng.Reset( 0x124 );
      int rows = 6;
      int cols = 3;
      var A = new DoubleMatrix( rows, cols, rng );

      // Construct a QR decomposition of A.
      var decomp = new DoubleQRDecomp( A );

      Console.WriteLine();

      // Look at the components of the factorization AP = QR.
      Console.WriteLine( "Q =" );
      Console.WriteLine( decomp.Q.ToTabDelimited( "G3" ) );

      Console.WriteLine( "R =" );
      Console.WriteLine( decomp.R.ToTabDelimited( "G3" ) );

      Console.WriteLine( "P =" );
      Console.WriteLine( decomp.P.ToTabDelimited() );

      // Q =
      // -0.187  0.516   0.0534
      // 0.0654  -0.362  0.0361
      // 0.142   0.331   -0.428
      // 0.22    -0.597  -0.507
      // 0.784   0.0122  0.56
      // -0.527  -0.37   0.492

      // R =
      // -1.18   0.264   0.183
      // 0       -0.869  -0.0966
      // 0       0       0.745

      // P =
      // 0       1       0
      // 0       0       1
      // 1       0       0

      // Usually, you do not need to access the components of the 
      // decomposition explicitly. There are member functions for
      // multiplying by the decomposition components without 
      // explicitly constructing them. As an example, lets solve
      // the least squares problem Ax = b using a QR decomposition.
      // To do this we write A = QR, compute the vector QTb 
      // (QT is the transpose of the matrix Q), and solve the upper-
      // triangular system Rx = QTb for x. (NOTE: the NMath Matrix library
      // contains class DoubleQRLeastSq for solving least squares problems
      // using a QR decomposition).
      var b = new DoubleVector( A.Rows, rng );
      DoubleVector w = decomp.QTx( b );
      DoubleVector x = decomp.Px( decomp.RInvx( w ) );

      Console.WriteLine( "Least squares solution =" );
      Console.WriteLine( x.ToString() );

      // Note that we had to multiply by the permutation matrix to 
      // get the components of the solution back in the right order
      // since pivoting was done.
      //
      // Suppose that you wanted to compute a QR decomposition without pivoting.
      // The class DoubleQRDecompServer provides control over the pivoting 
      // strategy used by the QR decomposition algorithm.
      var decompServer = new DoubleQRDecompServer();
      decompServer.Pivoting = false;

      // Use the server to get a decomposition that was computed without pivoting.
      decomp = decompServer.GetDecomp( A );

      // The permutation matrix should be I, the identity matrix.
      Console.WriteLine();
      Console.WriteLine( "P with no pivoting =" );
      Console.WriteLine(decomp.P.ToTabDelimited() );

      // Use the decomposition to solve the least squares problem again. This
      // time we do not need to multiply by P.
      w = decomp.QTx( b );
      x = decomp.RInvx( w );

      Console.WriteLine( "Least squares solution, no pivoting =" );
      Console.WriteLine( x.ToString() );

      // You can also set the columns of A to be initialcolumns. If the ith column
      // of A is set to be an initial column, it is moved to the beginning of AP 
      // before the computation, and fixed in place during the computation. If a 
      // column of A is set to be a freecolumn, it may be interchanged during the
      // computation with any other free column.
      decompServer.SetInitialColumn( 0 ); // Do not move the first column of A.
      decomp = decompServer.GetDecomp( A );
      decompServer.SetFreeColumn( 0 ); // Set the first column back to a free column.

      Console.WriteLine();
      Console.WriteLine( "Press Enter Key" );
      Console.Read();
    }
  }
}