VB.NET Iteratively Reweighted Least Sq 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 classes for solving iteratively reweighted
  ' least squares problems.
  Module IterativelyReweightedLeastSqExample

    Sub Main()

      ' Set up a least squares problem, Ax = b, with random data.
      Dim RNG As New RandGenUniform(-2, 2, 124)
      Dim Rows As Integer = 10
      Dim Cols As Integer = 2
      Dim A As New DoubleMatrix(Rows, Cols, RNG)
      Dim X As New DoubleVector(Cols, RNG)
      ' Fix up the right hand side b so that x
      ' is the exact solution, then throw in some outliers.
      Dim B As DoubleVector = NMathFunctions.Product(A, X)
      ' Throw in a few outliers...
      B(1) = 23
      B(4) = -10

      ' Create an iteratively reweighted least squares instance
      ' and use it to solve the problem using the default settings.
      ' The default weighting is DoubleBisquareWeightingFunction which
      ' uses the bisquare weighting algorithm.
      Dim Irls As New DoubleIterativelyReweightedLeastSq()

      ' Solve. The third parameter below specifies prepending a column of ones to the
      ' data in A representing a constant term in the model (which should come out
      ' to be zero in the solution from the way we cooked the data). Note that our
      ' input matrix A will not actually be changed.
      Dim Solution As DoubleVector = Irls.Solve(A, B, True)
      Console.WriteLine("Solution with bisquare weighting")
      Console.WriteLine(Solution.ToString())
      Console.WriteLine()
      Console.WriteLine("||residuals|| = " & Irls.Residuals.TwoNorm())
      Console.WriteLine()
      If (Irls.Iterations >= Irls.MaxIterations) Then
        Console.WriteLine("The algorithm did not converge in {0} iterations.", Irls.MaxIterations)
      Else
        Console.WriteLine("Algorithm converged in {0} iterations.", Irls.Iterations)
      End If

      ' Change some of the settings that control the iteration.
      Irls.MaxIterations = 300
      Irls.Tolerance = 0.0000001

      ' The convergence function is a delegate that may specified by the user for
      ' determining if the algorithm has converged and iteration terminated. The
      ' delegate takes as arguments the previous and current solutions and residuals
      ' and the tolerance and returns a bool. See the ResidualsChanged function 
      ' below.
      Dim ResidualsUnchanged As New DoubleIterativelyReweightedLeastSq.ToleranceMetFunction(AddressOf ResidualsUnchangedFunction)
      Irls.ConvergenceFunction = ResidualsUnchanged

      ' Change the weighting function used from the default bisquare weighting to the
      ' fair weighting function. See the class DoubleFairWeightingFunction for 
      ' particulars.
      Irls.WeightsFunction = New DoubleFairWeightingFunction()

      ' Solve the problem with the new settings.
      Solution = Irls.Solve(A, B, True)

      Console.WriteLine()
      Console.WriteLine("Solution with fair weighting is ")
      Console.WriteLine(Solution.ToString())
      Console.WriteLine()
      Console.WriteLine("||residuals|| = " & Irls.Residuals.TwoNorm())
      Console.WriteLine()
      If (Irls.Iterations >= Irls.MaxIterations) Then
        Console.WriteLine("The algorithm did not converge in {0} iterations.", Irls.MaxIterations)
      Else
        Console.WriteLine("Algorithm converged in {0} iterations.", Irls.Iterations)
      End If

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

    End Sub

    ' Convergence function for use in the iteratively reweighted least squares
    ' algorithm. This particular convergence function returns true when the residuals
    ' from the current iterations are relatively the same as the residuals in the
    ' previous iteration.
    Private Function ResidualsUnchangedFunction(ByVal Tolerance As Double, ByVal LastSolution As DoubleVector, _
      ByVal CurrenSolution As DoubleVector, ByVal LastResiduals As DoubleVector, ByVal CurrentResiduals As DoubleVector) As Boolean

      Dim MaxAbsDiff As Double = NMathFunctions.MaxAbsValue(CurrentResiduals - LastResiduals)
      Return (MaxAbsDiff / NMathFunctions.MaxAbsValue(CurrentResiduals)) < Tolerance

    End Function

  End Module
End Namespace