C# QR Decomp Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.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.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      int rows = 6;
      int cols = 3;
      DoubleMatrix A = new DoubleMatrix(rows, cols, rng);

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

      // Look at the components of the factorization AP = QR.
      Console.WriteLine();
      Console.WriteLine("Q = {0}", decomp.Q.ToString("F4"));

      Console.WriteLine();
      Console.WriteLine("R = {0}", decomp.R.ToString("F4"));

      Console.WriteLine();
      Console.WriteLine("P = {0}", decomp.P);

      // Q = 6x3 [ -0.1871  0.5159  0.0534  
      //            0.0654 -0.3618  0.0361 
      //            0.1424  0.3315 -0.4280
      //            0.2198 -0.5969 -0.5066
      //            0.7840  0.0122  0.5601  
      //           -0.5267 -0.3695  0.4922 ]
      //
      // R = 3x3 [ -1.1782  0.2637 0.1833 
      //            0.0000 -0.8685 -0.0966
      //            0.0000  0.0000 0.7454 ]
      //
      // P = 3x3 [ 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, let's 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-
      // triagular system Rx = QTb for x. (NOTE: the NMath Matrix library
      // contains class DoubleQRLeastSq for solving least squares problems
      // using a QR decomposition).
      DoubleVector b = new DoubleVector(A.Rows, rng);
      DoubleVector w = decomp.QTx(b);
      DoubleVector x = decomp.Px(decomp.RInvx(w));

      Console.WriteLine();
      Console.WriteLine("Least squares solution = {0}", 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 
      // stategy used by the QR decomposition algorithm.
      DoubleQRDecompServer 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: {0}", decomp.P.ToString());

      // 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();
      Console.WriteLine("Solution to the least squares problem with no pivoting:");
      Console.WriteLine(x.ToString());

      // You can also set the columns of A to be 'initial' columns. If the ith column
      // of A is set to be an intial 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 'free' column, 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();
    }
  }
}