VB Least Squares Example

Imports System

Imports CenterSpace.NMath.Core


Namespace CenterSpace.NMath.Examples.VisualBasic

  A .NET example in Visual Basic demonstrating the features of the classes for solving ordinary least
  squares problems.
  Module LeastSquaresExample

    Sub Main()

      A general m x n system with random entries.
      Dim Rng As New RandGenUniform(-1, 1)
      Rng.Reset(&H124)
      Dim Rows As Integer = 20
      Dim Cols As Integer = 4
      Dim A As New DoubleMatrix(Rows, Cols, Rng)

      Construct a class to solve the least squares problem using a
      Cholesky factorization of the normal equations.
      Dim CholLsq As 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 = " & CholLsq.IsGood)

      Create a random right-hand side to solve for.
      Dim B As New DoubleVector(A.Rows, Rng)

      Use Cholesky least squares method to solve.
      Dim X As DoubleVector = CholLsq.Solve(B)

      Look at the norm squared of the residual Ax - b.
      Dim ResidualNorm As Double = CholLsq.ResidualNormSqr(B)
      Console.WriteLine("Cholesky: residual = " & ResidualNorm)    5.89829729588437
      Console.WriteLine()

      Solve the same system with a QR decomposition and compare the results.
      Dim QRLsq As New DoubleQRLeastSq(A)
      Console.WriteLine("Rank of A from QR = " & QRLsq.Rank)
      X = QRLsq.Solve(B)
      ResidualNorm = QRLsq.ResidualNormSqr(B)
      Console.WriteLine("QR: residual = " & ResidualNorm)   5.89829729588437
      Console.WriteLine()

      How about SVD?
      Dim SVDLsq As New DoubleSVDLeastSq(A)
      Console.WriteLine("Rank of A from SVD = " & SVDLsq.Rank)
      X = SVDLsq.Solve(B)
      ResidualNorm = SVDLsq.ResidualNormSqr(B)
      Console.WriteLine("SVD: residual = " & 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.
      Dim NormalRNG As New RandGenNormal(0, 3.0E-23)
      NormalRNG.Reset(&H131)
      Dim SmallNormalDeviates As New DoubleVector(A.Rows, NormalRNG)
      A.Col(1)(Slice.All) = DoubleVector.Add(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 (Not CholLsq.IsGood) Then
        Console.WriteLine("Can not use Cholesky method on A. A is rank deficient.")
      Else
        X = CholLsq.Solve(B)
        ResidualNorm = CholLsq.ResidualNormSqr(B)
        Console.Write("Cholesky ill conditioned: residual = ")
        Console.WriteLine(ResidualNorm)
      End If
      Console.WriteLine()

      Solve the same system with a QR decomposition and compare the results.
      QRLsq.Factor(A)
      Console.WriteLine("Rank of A from QR = " & QRLsq.Rank)
      X = QRLsq.Solve(B)
      ResidualNorm = QRLsq.ResidualNormSqr(B)
      Console.WriteLine("QR  ill conditioned: residual = " & ResidualNorm)   6.03604367999183
      Console.WriteLine()

      How about SVD?
      SVDLsq.Factor(A)
      Console.WriteLine("Rank of A from SVD = " & SVDLsq.Rank)
      X = SVDLsq.Solve(B)
      ResidualNorm = SVDLsq.ResidualNormSqr(B)
      Console.Write("SVD  ill conditioned: residual = ")
      Console.WriteLine(ResidualNorm)   5.21216647666608

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

    End Sub
  End Module
End Namespace