VB.Net Program to Find a Function Minimum

Using Marquardt's Method

by Namir Shammas

The following program uses the Marquardt method to find the minimum of a function.

Click here to download a ZIP file containing the project files for this program.

The program prompts you to either use the predefined default input values or to enter the following:

1. The initial guesses for the optimum point for each variable

2. The gradient norm tolerance.

3. The function tolerance

4. The maximum number of iterations.

5. The Alpha factor (a large value like 1E+4 is suitable).

6. The value for C1 and C2 factors (values for C1 = 0.1 and C2 = 10 are generally adequate)

In case you choose the default input values, the program displays these values and proceeds to find the optimum point. In the case you select being prompted, the program displays the name of each input variable along with its default value. You can then either enter a new value or simply press Enter to use the default value. This approach allows you to quickly and efficiently change only a few input values if you so desire.

The program displays the following final results:

1. The coordinates of the minimum point.

2. The minimum function value.

3. The number of iterations

The current code finds the minimum for the following function:

f(x1,x2) = x1 - x2 + 2 * x1 ^ 2 + 2 * x1 * x2 + x2 ^ 2

Using, for each variable, an initial value of 0,  initial step size of 0.1, minimum step size of 1e-7, and using a function tolerance of 1e-7. Here is the sample console screen:

Here is the listing for the main module.  The module contains several test functions: 

Module Module1

  Sub Main()
    Dim nNumVars As Integer = 2
    Dim fX() As Double = {0, 0}
    Dim fParam() As Double = {0.5, 0.5}
    Dim fToler() As Double = {0.000001, 0.000001}
    Dim nIter As Integer = 0
    Dim nMaxIter As Integer = 1000
    Dim fEpsFx As Double = 0.0000001
    Dim fAlpha As Double = 10000.0
    Dim C1 As Double = 0.1
    Dim C2 As Double = 10
    Dim I As Integer
    Dim fBestF
    Dim sAnswer As String, sErrorMsg As String = ""
    Dim oOpt As CMarquardt1
    Dim MyFx As MyFxDelegate = AddressOf Fx3
    Dim SayFx As SayFxDelegate = AddressOf SayFx3

    oOpt = New CMarquardt1

    Console.WriteLine("Marquardt's Method for Optimization")
    Console.WriteLine("Finding the minimum of function:")
    Console.WriteLine(SayFx())
    Console.Write("Use default input values? (Y/N) ")
    sAnswer = Console.ReadLine()
    If sAnswer.ToUpper() = "Y" Then
      For I = 0 To nNumVars - 1
        Console.WriteLine("X({0}) = {1}", I + 1, fX(I))
        Console.WriteLine("Tolerance({0}) = {1}", I + 1, fToler(I))
      Next
      Console.WriteLine("Function tolerance = {0}", fEpsFx)
      Console.WriteLine("Maxumum cycles = {0}", nMaxIter)
      Console.WriteLine("Alpha = {0}", fAlpha)
      Console.WriteLine("C1 = {0}", C1)
      Console.WriteLine("C2 = {0}", C2)
    Else
      For I = 0 To nNumVars - 1
        fX(I) = GetIndexedDblInput("X", I + 1, fX(I))
        fToler(I) = GetIndexedDblInput("Tolerance", I + 1, fToler(I))
      Next
      fEpsFx = GetDblInput("Function tolerance", fEpsFx)
      nMaxIter = GetIntInput("Maxumum cycles", nMaxIter)
      fAlpha = GetDblInput("Alpha", fAlpha)
      C1 = GetDblInput("C1", C1)
      C2 = GetDblInput("C1", C2)
    End If

    Console.WriteLine("******** FINAL RESULTS *************")
    fBestF = oOpt.CalcOptim(nNumVars, fX, fParam, fToler, fEpsFx, nMaxIter, fAlpha, C1, C2, nIter, sErrorMsg, MyFx)
    If sErrorMsg.Length > 0 Then
      Console.WriteLine("** NOTE: {0} ***", sErrorMsg)
    End If
    Console.WriteLine("Optimum at")
    For I = 0 To nNumVars - 1
      Console.WriteLine("X({0}) = {1}", I + 1, fX(I))
    Next
    Console.WriteLine("Function value = {0}", fBestF)
    Console.WriteLine("Number of iterations = {0}", nIter)
    Console.WriteLine()
    Console.Write("Press Enter to end the program ...")
    Console.ReadLine()
  End Sub

  Function GetDblInput(ByVal sPrompt As String, ByVal fDefInput As Double) As Double
    Dim sInput As String

    Console.Write("{0}? ({1}): ", sPrompt, fDefInput)
    sInput = Console.ReadLine()
    If sInput.Trim().Length > 0 Then
      Return Double.Parse(sInput)
    Else
      Return fDefInput
    End If
  End Function

  Function GetIntInput(ByVal sPrompt As String, ByVal nDefInput As Integer) As Integer
    Dim sInput As String

    Console.Write("{0}? ({1}): ", sPrompt, nDefInput)
    sInput = Console.ReadLine()
    If sInput.Trim().Length > 0 Then
      Return Integer.Parse(sInput)
    Else
      Return nDefInput
    End If
  End Function

  Function GetIndexedDblInput(ByVal sPrompt As String, ByVal nIndex As Integer, ByVal fDefInput As Double) As Double
    Dim sInput As String

    Console.Write("{0}({1})? ({2}): ", sPrompt, nIndex, fDefInput)
    sInput = Console.ReadLine()
    If sInput.Trim().Length > 0 Then
      Return Double.Parse(sInput)
    Else
      Return fDefInput
    End If
  End Function

  Function GetIndexedIntInput(ByVal sPrompt As String, ByVal nIndex As Integer, ByVal nDefInput As Integer) As Integer
    Dim sInput As String

    Console.Write("{0}({1})? ({2}): ", sPrompt, nIndex, nDefInput)
    sInput = Console.ReadLine()
    If sInput.Trim().Length > 0 Then
      Return Integer.Parse(sInput)
    Else
      Return nDefInput
    End If
  End Function

  Function SayFx1() As String
    Return "F(X) = 10 + (X(1) - 2) ^ 2 + (X(2) + 5) ^ 2"
  End Function

  Function Fx1(ByVal N As Integer, ByRef X() As Double, ByRef fParam() As Double) As Double
    Return 10 + (X(0) - 2) ^ 2 + (X(1) + 5) ^ 2
  End Function

  Function SayFx2() As String
    Return "F(X) = 100 * (X(1) - X(2) ^ 2) ^ 2 + (X(2) - 1) ^ 2"
  End Function

  Function Fx2(ByVal N As Integer, ByRef X() As Double, ByRef fParam() As Double) As Double
    Return 100 * (X(0) - X(1) ^ 2) ^ 2 + (X(1) - 1) ^ 2
  End Function

  Function SayFx3() As String
    Return "F(X) = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2"
  End Function

  Function Fx3(ByVal N As Integer, ByRef X() As Double, ByRef fParam() As Double) As Double
    Return X(0) - X(1) + 2 * X(0) ^ 2 + 2 * X(0) * X(1) + X(1) ^ 2
  End Function

End Module

Notice that the user-defined functions have accompanying helper functions to display the mathematical expression of the function being optimized. For example, function Fx1 has the helper function SayFx1 to list the function optimized in Fx1. Please observe the following rules::

The program uses the following class to optimize the objective function:

Public Delegate Function MyFxDelegate(ByVal nNumVars As Integer, ByRef fX() As Double, ByRef fParam() As Double) As Double
Public Delegate Function SayFxDelegate() As String

Public Class CMarquardt1
  Dim m_MyFx As MyFxDelegate

  Protected Function MyFxEx(ByVal nNumVars As Integer, _
              ByRef fX() As Double, ByRef fParam() As Double, _
              ByRef fDeltaX() As Double, ByVal fLambda As Double) As Double
    Dim I As Integer
    Dim fXX(nNumVars) As Double

    For I = 0 To nNumVars - 1
      fXX(I) = fX(I) + fLambda * fDeltaX(I)
    Next I

    MyFxEx = m_MyFx(nNumVars, fXX, fParam)
  End Function

  Protected Function FirstDeriv(ByVal nNumVars As Integer, _
               ByRef fX() As Double, ByRef fParam() As Double, _
               ByVal nIdxI As Integer) As Double
    Dim fXt, h, Fp, Fm As Double

    fXt = fX(nIdxI)
    h = 0.01 * (1 + Math.Abs(fXt))
    fX(nIdxI) = fXt + h
    Fp = m_MyFx(nNumVars, fX, fParam)
    fX(nIdxI) = fXt - h
    Fm = m_MyFx(nNumVars, fX, fParam)
    fX(nIdxI) = fXt
    FirstDeriv = (Fp - Fm) / 2 / h
  End Function

  Protected Function SecondDeriv(ByVal nNumVars As Integer, _
                 ByRef fX() As Double, ByRef fParam() As Double, _
                 ByVal nIdxI As Integer, ByVal nIdxJ As Integer) As Double
    Dim fXt, fYt, fHX, fHY, F0, Fp, Fm As Double
    Dim Fpp, Fmm, Fpm, Fmp, fResult As Double

    ' calculate second derivative?
    If nIdxI = nIdxJ Then
      F0 = m_MyFx(nNumVars, fX, fParam)
      fXt = fX(nIdxI)
      fHX = 0.01 * (1 + Math.Abs(fXt))
      fX(nIdxI) = fXt + fHX
      Fp = m_MyFx(nNumVars, fX, fParam)
      fX(nIdxI) = fXt - fHX
      Fm = m_MyFx(nNumVars, fX, fParam)
      fX(nIdxI) = fXt
      fResult = (Fp - 2 * F0 + Fm) / fHX ^ 2
    Else
      fXt = fX(nIdxI)
      fYt = fX(nIdxJ)
      fHX = 0.01 * (1 + Math.Abs(fXt))
      fHY = 0.01 * (1 + Math.Abs(fYt))
      ' calculate Fpp
      fX(nIdxI) = fXt + fHX
      fX(nIdxJ) = fYt + fHY
      Fpp = m_MyFx(nNumVars, fX, fParam)
      ' calculate Fmm
      fX(nIdxI) = fXt - fHX
      fX(nIdxJ) = fYt - fHY
      Fmm = m_MyFx(nNumVars, fX, fParam)
      ' calculate Fpm
      fX(nIdxI) = fXt + fHX
      fX(nIdxJ) = fYt - fHY
      Fpm = m_MyFx(nNumVars, fX, fParam)
      ' calculate Fmp
      fX(nIdxI) = fXt - fHX
      fX(nIdxJ) = fYt + fHY
      Fmp = m_MyFx(nNumVars, fX, fParam)
      fX(nIdxI) = fXt
      fX(nIdxJ) = fYt
      fResult = (Fpp - Fmp - Fpm + Fmm) / (4 * fHX * fHY)
    End If

    Return fResult
  End Function


  Protected Sub GetFirstDerives(ByVal nNumVars As Integer, _
                  ByRef fX() As Double, ByRef fParam() As Double, _
                  ByRef fFirstDerivX() As Double)

    Dim I As Integer

    For I = 0 To nNumVars - 1
      fFirstDerivX(I) = FirstDeriv(nNumVars, fX, fParam, I)
    Next I
  End Sub

  Protected Sub GetSecondDerives(ByVal nNumVars As Integer, _
                  ByRef fX() As Double, ByRef fParam() As Double, _
                  ByRef fSecondDerivX(,) As Double)

    Dim I, J As Integer

    For I = 0 To nNumVars - 1
      For J = 0 To nNumVars - 1
        fSecondDerivX(I, J) = SecondDeriv(nNumVars, fX, fParam, I, J)
      Next J
    Next I
  End Sub

  Protected Function LinSearch_DirectSearch(ByVal nNumVars As Integer, ByRef fX() As Double, ByRef fParam() As Double, _
            ByRef fLambda As Double, ByRef fDeltaX() As Double, ByVal fInitStep As Double, ByVal fMinStep As Double) As Boolean
    Dim F1, F2 As Double

    F1 = MyFxEx(nNumVars, fX, fParam, fDeltaX, fLambda)

    Do
      F2 = MyFxEx(nNumVars, fX, fParam, fDeltaX, fLambda + fInitStep)
      If F2 < F1 Then
        F1 = F2
        fLambda += fInitStep
      Else
        F2 = MyFxEx(nNumVars, fX, fParam, fDeltaX, fLambda - fInitStep)
        If F2 < F1 Then
          F1 = F2
          fLambda -= fInitStep
        Else
          ' reduce search step size
          fInitStep /= 10
        End If
      End If
    Loop Until fInitStep < fMinStep

    Return True

  End Function

  Public Function CalcOptim(ByVal nNumVars As Integer, ByRef fX() As Double, ByRef fParam() As Double, _
                  ByRef fToler() As Double, ByVal fEpsFx As Double, ByVal nMaxIter As Integer, _
                  ByRef fAlpha As Double, ByVal C1 As Double, ByVal C2 As Double, _
                  ByRef nIter As Integer, ByRef sErrorMsg As String, _
                  ByVal MyFx As MyFxDelegate) As Double

    Const ALPHA_MAX = 1.0E+100

    Dim I As Integer
    Dim F1, F2, fNorm As Double
    Dim nInnerLoopIter As Integer
    Dim bStop As Boolean
    Dim fX2(nNumVars) As Double
    Dim g(nNumVars) As Double
    Dim G2(nNumVars) As Double
    Dim DeltaX(nNumVars) As Double
    Dim J(nNumVars, nNumVars) As Double
    Dim J2(nNumVars, nNumVars) As Double
    Dim fIdenMat(nNumVars, nNumVars) As Double
    Dim fIdenMat2(nNumVars, nNumVars) As Double
    Dim nIndex(nNumVars) As Integer

    On Error GoTo HandleErr

    m_MyFx = MyFx

    ' validate iterations parameters
    If C1 <= 0 Or C1 >= 1 Then C1 = 0.1
    If C2 <= 1 Then C2 = 10

    F2 = MyFx(nNumVars, fX, fParam)

    nIter = 1
    Do
      nIter += 1
      If nIter > nMaxIter Then
        sErrorMsg = "Reached maximum iterations limit"
        Exit Do
      End If

      ' copy last Fx value
      F1 = F2

      GetFirstDerives(nNumVars, fX, fParam, g)
      MatrixLibVb.DuplicateVector(g, G2)

      ' test if gradient is shallow enough
      fNorm = 0
      For I = 0 To nNumVars - 1
        fNorm += g(I) ^ 2
      Next I
      fNorm = Math.Sqrt(fNorm)
      If fNorm < fEpsFx Then Exit Do

      MatrixLibVb.DuplicateVector(fX, fX2) ' copy fX
      ' get matrix J
      GetSecondDerives(nNumVars, fX, fParam, J)
      ' create identitty matrix
      MatrixLibVb.MatIdentity(fIdenMat)
      ' initialize inner-loop counter
      nInnerLoopIter = 0

      Do
        nInnerLoopIter += 1
        If nInnerLoopIter > nMaxIter Then Exit Function

        MatrixLibVb.DuplicateVector(G2, g)
        MatrixLibVb.MatMultVal(fIdenMat, fAlpha, fIdenMat2)
        MatrixLibVb.MatAddMat(J, fIdenMat2, J2)
        MatrixLibVb.MV_LUDecomp(J2, nIndex, nNumVars)
        MatrixLibVb.MV_LUBackSubst(J2, nIndex, nNumVars, g)
        For I = 0 To nNumVars - 1
          DeltaX(I) = g(I)
          fX(I) -= DeltaX(I)
        Next I

        F2 = MyFx(nNumVars, fX, fParam)

        If F2 < F1 Then
          fAlpha = C1 * fAlpha
          bStop = True
        Else
          fAlpha = C2 * fAlpha
          MatrixLibVb.DuplicateVector(fX2, fX) ' restore array fX
          bStop = False
          If fAlpha >= ALPHA_MAX Then Return F2
        End If
      Loop Until bStop

    Loop

ExitProc:
    Return F2

HandleErr:
    sErrorMsg = "Calculation error: " & Err.Description
    Resume ExitProc
  End Function

End Class

BACK

Copyright (c) Namir Shammas. All rights reserved.