C# Least Squares 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 classes for solving ordinary least
  /// squares problems.
  /// </summary>
  class LeastSquaresExample
  {

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

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

      // 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.
      DoubleVector 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

      // Solve the same system with a QR decomposition and compare the results.
      DoubleQRLeastSq 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

      // How about SVD?
      DoubleSVDLeastSq 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

      // 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.
      RandGenNormal normalRNG = new RandGenNormal(0, 3e-23);
      normalRNG.Reset(0x131);
      DoubleVector 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);
      }

      // 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);

      // 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();
    }
  }
}