VB.NET QR Decomp Example

Imports System

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

Namespace CenterSpace.NMath.Matrix.Examples.VisualBasic

  ' A .NET example in VB.NET demonstrating the features of the QR decomposition classes.
  Module QRDecompExample

    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 = 6
      Dim Cols As Integer = 3
      Dim A As New DoubleMatrix(Rows, Cols, Rng)

      ' Construct a QR decomposition of A.
      Dim Decomp As New DoubleQRDecomp(A)

      ' Look at the components of the factorization AP = QR.
      Console.WriteLine()
      Console.Write("Q = ")
      Console.WriteLine(Decomp.Q.ToString("F4"))

      Console.WriteLine()
      Console.Write("R = ")
      Console.WriteLine(Decomp.R.ToString("F4"))

      Console.WriteLine()
      Console.Write("P = ")
      Console.WriteLine(Decomp.P)

      ' Q = 6x3 [ -0.1871  0.5159  0.0534  
      '            0.0654 -0.3618  0.0361 
      '            0.1424  0.3315 -0.4280
      '            0.2198 -0.5969 -0.5066
      '            0.7840  0.0122  0.5601  
      '           -0.5267 -0.3695  0.4922 ]
      '
      ' R = 3x3 [ -1.1782  0.2637 0.1833 
      '            0.0000 -0.8685 -0.0966
      '            0.0000  0.0000 0.7454 ]
      '
      ' P = 3x3 [ 0 1 0 
      '           0 0 1 
      '           1 0 0 ]

      ' Usually, you do not need to access the components of the 
      ' decomposition explicitly. There are member functions for
      ' multiplying by the decomposition components without 
      ' explicitly constructing them. As an example, let's solve
      ' the least squares problem Ax = b using a QR decomposition.
      ' To do this we write A = QR, compute the vector QTb 
      ' (QT is the transpose of the matrix Q), and solve the upper-
      ' triagular system Rx = QTb for x. (NOTE: the NMath Matrix library
      ' contains class DoubleQRLeastSq for solving least squares problems
      ' using a QR decomposition).
      Dim B As New DoubleVector(A.Rows, Rng)
      Dim W As DoubleVector = Decomp.QTx(B)
      Dim X As DoubleVector = Decomp.Px(Decomp.RInvx(W))

      Console.WriteLine()
      Console.Write("Least squares solution = ")
      Console.WriteLine(X.ToString())

      ' Note that we had to multiply by the permutation matrix to 
      ' get the components of the solution back in the right order
      ' since pivoting was done.
      '
      ' Suppose that you wanted to compute a QR decomposition without pivoting.
      ' The class DoubleQRDecompServer provides control over the pivoting 
      ' stategy used by the QR decomposition algorithm.
      Dim DecompServer As New DoubleQRDecompServer()
      DecompServer.Pivoting = False

      ' Use the server to get a decomposition that was computed without pivoting.
      Decomp = DecompServer.GetDecomp(A)

      ' The permutation matrix should be I, the identity matrix.
      Console.WriteLine()
      Console.Write("P with no pivoting: ")
      Console.WriteLine(Decomp.P.ToString())

      ' Use the decomposition to solve the least squares problem again. This
      ' time we do not need to multiply by P.
      W = Decomp.QTx(B)
      X = Decomp.RInvx(W)
      Console.WriteLine()
      Console.WriteLine("Solution to the least squares problem with no pivoting:")
      Console.WriteLine(X.ToString())

      ' You can also set the columns of A to be 'initial' columns. If the ith column
      ' of A is set to be an intial column, it is moved to the beginning of AP 
      ' before the computation, and fixed in place during the computation. If a 
      ' column of A is set to be a 'free' column, it may be interchanged during the
      ' computation with any other free column.
      DecompServer.SetInitialColumn(0)   ' Do not move the first column of A.
      Decomp = DecompServer.GetDecomp(A)
      DecompServer.SetFreeColumn(0)   ' Set the first column back to a free column.

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

    End Sub
  End Module
End Namespace