C# Least Squares Example

using System;

using CenterSpace.NMath.Core;


namespace CenterSpace.NMath.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the classes for solving ordinary least
  /// squares problems.
  /// </summary>
  class LeastSquaresExample
  {

    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 = 20;
      int cols = 4;
      var A = new DoubleMatrix( rows, cols, rng );

      // Construct a class to solve the least squares problem using a
      // Cholesky factorization of the normal equations.
      var cholLsq = new DoubleCholeskyLeastSq( A );

      Console.WriteLine();
      
      // Make sure the Cholesky factorization is good (for example, that matrix A
      // has full rank).
      Console.WriteLine( "Cholesky least squares is good = {0}", cholLsq.IsGood );

      // Create a random right-hand side to solve for.
      var b = new DoubleVector( A.Rows, rng );

      // Use Cholesky least squares method to solve.
      DoubleVector x = cholLsq.Solve( b );

      // Look at the norm squared of the residual Ax - b.
      double residualNorm = cholLsq.ResidualNormSqr( b );
      Console.WriteLine( "Cholesky: residual = {0}", residualNorm ); // 5.89829729588437
      Console.WriteLine();

      // Solve the same system with a QR decomposition and compare the results.
      var qrLsq = new DoubleQRLeastSq( A );
      Console.WriteLine( "Rank of A from QR = {0}", qrLsq.Rank );
      x = qrLsq.Solve( b );
      residualNorm = qrLsq.ResidualNormSqr( b );
      Console.WriteLine( "QR: residual = {0}", residualNorm ); // 5.89829729588437
      Console.WriteLine();

      // How about SVD?
      var svdLsq = new DoubleSVDLeastSq( A );
      Console.WriteLine( "Rank of A from SVD = {0}", svdLsq.Rank );
      x = svdLsq.Solve( b );
      residualNorm = svdLsq.ResidualNormSqr( b );
      Console.WriteLine( "SVD: residual = {0}", residualNorm ); // 5.89829729588438
      Console.WriteLine();

      // All methods give about the same residual. But what if A is nearly
      // rank deficient? Set the 2nd column of A to be almost equal to the 3rd.
      var normalRNG = new RandGenNormal( 0, 3e-23 );
      normalRNG.Reset( 0x131 );
      var smallNormalDeviates = new DoubleVector( A.Rows, normalRNG );
      A.Col( 1 )[Slice.All] = A.Col( 2 ) + smallNormalDeviates;

      // Use Cholesky least squares method to solve ill-conditioned.
      cholLsq.Factor( A );
      // Depending on your hardware, A could be so close to being rank
      // deficient that the Cholesky factorization fails. 
      if ( !cholLsq.IsGood )
      {
        Console.WriteLine( "Can not use Cholesky method on A. A is rank deficient." );
      }
      else
      {
        x = cholLsq.Solve( b );
        residualNorm = cholLsq.ResidualNormSqr( b );
        Console.WriteLine( "Cholesky ill conditioned: residual = {0}", residualNorm );
      }
      Console.WriteLine();

      // Solve the same system with a QR decomposition and compare the results.
      qrLsq.Factor( A );
      Console.WriteLine( "Rank of A from QR = {0}", qrLsq.Rank );
      x = qrLsq.Solve( b );
      residualNorm = qrLsq.ResidualNormSqr( b );
      Console.WriteLine( "QR  ill conditioned: residual = {0}", residualNorm );
      Console.WriteLine();

      // How about SVD?
      svdLsq.Factor( A );
      Console.WriteLine( "Rank of A from SVD = {0}", svdLsq.Rank );
      x = svdLsq.Solve( b );
      residualNorm = svdLsq.ResidualNormSqr( b );
      Console.WriteLine( "SVD  ill conditioned: residual = {0}", residualNorm );

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