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
Copyright (c) Namir Shammas. All rights reserved.