Eigen/MinPack: non linear optimization#

Multidimensional Rootfinding: Powell Hybrid#

ctxboost.PowellHybrid(f, fjac, guess)#

This is a modified version of Powell’s Hybrid method as implemented in the hybrj algorithm in minpack. The Hybrid algorithm retains the fast convergence of Newton’s method but will also reduce the residual when Newton’s method is unreliable. The algorithm uses a generalized trust region to keep each step under control.

See also: Moré et al. [439], Eigen [223], Wikipedia [1558], Powell [498].

Parameters:

f:: a callback matrix function defining the system.
fjac:: a callback matrix function defining the Jacobian of the system system.
guess:: a matrix containing the initial guess for the solution

Results:

x1:: a matrix containing the solution.

Example:

This is an example from the original manual of MINPACK (Moré et al. [439]).

The routine in Python:

def demo_PowellHybridCtx(ctx):

    # The function to optimize
    def XmatHybrd(x, fvec):
        n = x.size
        for k in range(n):
            temp = (3.0 - 2.0 * x[k]) * x[k]
            temp1 = ctx.zero()
            if (k!=0): temp1 = x[k-1]
            temp2 = ctx.zero()
            if (k != n-1): temp2 = x[k+1]
            fvec[k] = temp - temp1 - 2.0*temp2 + 1.0

    # The Jacobi function
    def XmatHybrdJ(x, fjac):
        n = x.size
        for k in range(n):
            for j in range(n):
                fjac[k, j] = ctx.zero()
            fjac[k, k] = 3.0 - 4.0 * x[k]
            if (k != 0): fjac[k, k - 1] = -1
            if (k != n - 1): fjac[k, k + 1] = -2

    # This defines the start vector
    n = 9
    matInput = ctx.matZeros(n,1)
    matInput[0] = 1
    matInput[1] = 2  # entries 2 .. 8 are 0.
    print('matInput: \n', matInput)

    # This executes the PowellHybrd solver
    matRes = ctx.powellHybrd(XmatHybrd, XmatHybrdJ, matInput)

    # Check the result
    Y = ctx.matZeros(n,1)
    XmatHybrd(matRes, Y)
    print('matRes: \n', matRes)
    print('Y: \n', Y)

The same routine in Visual Basic:

 Public Sub XmatHybrd(x As xrealMatrix, fvec As xrealMatrix)
     Dim n As Int32 = x.size
     For k As Int32 = 0 To n - 1
         Dim temp = (3.0 - 2.0 * x(k)) * x(k)
         Dim temp1 = xreal.Zero()
         If (k<>0) Then temp1 = x(k-1)
         Dim temp2 = xreal.Zero()
         If (k <> n-1) Then temp2 = x(k+1)
         fvec(k) = temp - temp1 - 2.0*temp2 + 1.0
     Next
 End Sub

 Public Sub XmatHybrdJ(x As xrealMatrix, fjac As xrealMatrix)
     Dim n As Int32 = x.size
     For k As Int32 = 0 To n - 1
         For j As Int32 = 0 To n - 1
             fjac(k, j) = 0
         Next
         fjac(k, k) = 3.0 - 4.0 * x(k)
         if (k <> 0) Then fjac(k, k - 1) = -1
         if (k <> n - 1) Then fjac(k, k + 1) = -2
     Next
 End Sub

 Sub DemoPowellHybrdClassxreal()
     Dim n As Int32 = 9
     Dim matInput = xreal.MatZeros(n,1)
     matInput(0) = 1
     matInput(1) = 2  ' entries 2 .. 8 are 0.
     Dim matX = xreal.PowellHybrd(AddressOf XmatHybrd, AddressOf XmatHybrdJ, matInput)
     Console.WriteLine("")
     matX.print("X (solution):", 10)
     Dim matEval = xreal.MatZeros(n,1)
     XmatHybrd(matX, matEval)
     matEval.print("matEval =  F(X=solution):", 10)
End Sub

The output of these routines:

Hello DemoPowellHybrdClassDbl()

X (solution):
-0.570654511600659,
-0.681628342291231,
-0.701732452563471,
-0.704212940083752,
-0.701369047627289,
-0.691865643379914,
-0.665792012154689,
-0.596034201280817,
-0.416412062998472,

matEval =  F(X=solution):
6.56011067690088E-09,
-4.17547307840493E-09,
-5.19316567526573E-09,
-2.39601338769546E-09,
2.02249372804886E-09,
4.81791939677123E-09,
2.57950016901987E-09,
-3.88373844195655E-09,
-1.35886191188206E-10,

Nonlinear LeastSquares: Levenberg-Marquardt#

ctxboost.LevenbergMarquardt(f, fjac, matInput)#

This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. The algorithm uses a generalized trust region to keep each step under control.

See also: Moré et al. [439], Eigen [223], Wikipedia [1552].

Parameters:

f:: a callback matrix function defining the system.
fjac:: a callback matrix function defining the Jacobian of the system system.
guess:: a matrix containing the initial guess for the solution

Results:

x1:: a matrix containing the solution.

Example:

This is an example from the original manual of MINPACK (\(Moré1980\)).

The routine in Python:

def demo_LevenbergCtx(ctx):

    # The function to optimize
    def XmatLM(x, fvec):
        m = 15
        y = ctx.matZeros(m,1)
        y[0] = 1.4e-1
        y[1] = 1.8e-1
        y[2] = 2.2e-1
        y[3] = 2.5e-1
        y[4] = 2.9e-1
        y[5] = 3.2e-1
        y[6] = 3.5e-1
        y[7] = 3.9e-1
        y[8] = 3.7e-1
        y[9] = 5.8e-1
        y[10] = 7.3e-1
        y[11] = 9.6e-1
        y[12] = 1.34e0
        y[13] = 2.1e0
        y[14] = 4.39e0
        for i in range(m):
            tmp1 = i + 1
            tmp2 = m - i
            tmp3 = tmp1
            if (i >= 8): tmp3 = tmp2
            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3))

    # The Jacobi function
    def XmatLMJ(x, fjac):
        m = 15
        for i in range(m):
            tmp1 = i + 1
            tmp2 = m - i
            tmp3 = tmp1
            if (i >= 8): tmp3 = tmp2
            tmp4 = (x[1] * tmp2 + x[2] * tmp3)
            tmp4 = tmp4 * tmp4
            fjac[i, 0] = -1
            fjac[i, 1] = tmp1 * tmp2 / tmp4
            fjac[i, 2] = tmp1 * tmp3 / tmp4

    # This defines the start vector
    n = 3; m = 15
    matInput = ctx.matZeros(n,1)
    matInput[0] = 1
    matInput[1] = 2  # entries 2 .. 8 are 0.
    print('matInput: \n', matInput)
    # This executes the Levenberg solver
    matRes = ctx.levenberg(XmatLM, XmatLMJ, matInput)
    # Check the result
    Y = ctx.matZeros(m,1)
    XmatLM(matRes, Y)
    print('matRes: \n', matRes)
    print('Y: \n', Y)

The same routine in Visual Basic:

Public Sub XmatLM(x As xrealMatrix, fvec As xrealMatrix)
    Dim m As Int32 = 15
    Dim y = xreal.MatZeros(m,1)
    y(0) = 1.4e-1
    y(1) = 1.8e-1
    y(2) = 2.2e-1
    y(3) = 2.5e-1
    y(4) = 2.9e-1
    y(5) = 3.2e-1
    y(6) = 3.5e-1
    y(7) = 3.9e-1
    y(8) = 3.7e-1
    y(9) = 5.8e-1
    y(10) = 7.3e-1
    y(11) = 9.6e-1
    y(12) = 1.34e0
    y(13) = 2.1e0
    y(14) = 4.39e0
    For i As Int32 = 0 To m - 1
        Dim tmp1 = i + 1
        Dim tmp2 = m - i
        Dim tmp3 = tmp1
        If (i >= 8) Then tmp3 = tmp2
        fvec(i) = y(i) - (x(0) + tmp1/(x(1)*tmp2 + x(2)*tmp3))
    Next
End Sub

Public Sub XmatLMJ(x As xrealMatrix, fjac As xrealMatrix)
    Dim m As Int32 = 15
    For i As Int32 = 0 To m - 1
        Dim tmp1 = i + 1
        Dim tmp2 = m - i
        Dim tmp3 = tmp1
        If (i >= 8) Then tmp3 = tmp2
        Dim tmp4 = (x(1) * tmp2 + x(2) * tmp3)
        tmp4 = tmp4 * tmp4
        fjac(i, 0) = -1
        fjac(i, 1) = tmp1 * tmp2 / tmp4
        fjac(i, 2) = tmp1 * tmp3 / tmp4
    Next
End Sub

Sub DemoLevenbergClassxreal()
    Dim n As Int32 = 3
    Dim m As Int32 = 15
    Dim matInput = xreal.MatZeros(n,1)
    matInput(0) = 1
    matInput(1) = 2  ' entries 2 .. 8 are 0.
    Dim matRes = xreal.Levenberg(AddressOf XmatLM, AddressOf XmatLMJ, matInput)
    Console.WriteLine("")
    matRes.print("X (solution):", 10)
    Dim Y = xreal.MatZeros(m,1)
    XmatLM(matRes, Y)
    Y.print("Y =  F(X=solution):", 10)
End Sub

This produces the following output:

Hello DemoLevenbergClassDbl()

X (solution):
0.0824105765758334,
1.1330366534715,
2.34369463894115,

matEval =  F(X=solution):
0.00588109515673704,
0.000265360346254795,
-0.000274673051589042,
-0.00654152299741256,
0.000823003778696318,
0.00129950005693674,
0.00446310734534455,
0.0199629386905548,
-0.0822160569476201,
0.0182119488681469,
0.0148111570102206,
0.0147099692233311,
0.0112079895785155,
0.00420403028888439,
-0.00680784758001085,

Eigen/MinPack: non linear optimization

Contents

Eigen/MinPack: non linear optimization#

Multidimensional Rootfinding: Powell Hybrid#

Nonlinear LeastSquares: Levenberg-Marquardt#