VB Least Squares Example

← All NMath Code Examples

 

Imports System

Imports CenterSpace.NMath.Core
Imports CenterSpace.NMath.Matrix

Namespace CenterSpace.NMath.Matrix.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

← All NMath Code Examples
Top