# C# Least Squares Example

← All NMath Code Examples

```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();

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

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