HP-71B 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.

The program prompts you to enter for each variable (i.e. dimension):

1. The gradient norm tolerance.

2. The maximum number of iterations.

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

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

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

6. The tolerance values for each variable.

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

Here is a sample session. The first figure shows the input.:

The second figure shows the output:

Here is the BASIC listing:  

10 ! Marquardt's Method for Optimization
20 !
30 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
40 ! DECISION-MAKING CONSTRUCTS
50 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 N = 2
80 A9 = 1E+100 ! Alpha max
90 !
100 DIM X(N), X2(N), G(N)
110 DIM T0(N), D(N), J(N, N)
120 DIM O1(N,N), O2(N,N)
130 !
140 DISP "MARQUARDT OPTIMIZATION METHOD"
150 INPUT "Enter gradient norm tolerance value? "; E
160 INPUT "Enter maximum number iterations? "; I9
170 INPUT "Enter value for Alpha factor? "; A1
180 INPUT "Enter value for C1 factor? "; C1
190 INPUT "Enter value for C2 factor? "; C2
200 For I = 1 To N
210 DISP "Enter guess for X(";I;")";
220 INPUT X(I)
230 DISP "Enter tolerance for X(";I;")";
240 INPUT T0(I)
250 Next I
260 !
270 CALL MyFx(X, N, F2)
280 !
290 I0 = 0
300 REPEAT
310 I0 = I0 + 1
320 If I0 > I9 Then
330 DISP "Reached maximum iterations limit"
340 LEAVE
350 End If
360 !
370 ! copy last Fx value
380 F1 = F2
390 !
400 CALL GetDrv1(N, X, G)
410 !
420 ! test if gradient is shallow enough
430 N0 = 0
440 For I = 1 To N
450 N0 = N0 + G(I)^2
460 Next I
470 N0 = Sqr(N0)
480 If N0 < E Then LEAVE
490 !
500 ! copy vector X
510 MAT X2 = X
520 ! get matrix J
530 CALL GetDrv2(N, X, J)
540 ! create identity matrix
550 MAT O1 = IDN(N,N)
560 ! Initialize inner loop
570 I1 = 0
580 !
590 REPEAT
600 I1 = I1 + 1
610 If I1 > I9 Then LEAVE
620 !
630 MAT O2 = (A1) * O1
640 MAT J = J + O2
650 MAT J = INV(J)
660 MAT G = J * G
670 !
680 MAT D = G
690 MAT X = X - D
700 !
710 CALL MyFx(X, N, F2)
720 !
730 DISP "F = ";F2;" ";
740 For I = 1 To N
750   DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
760 Next I
770 DISP
780 !
790 If F2 < F1 Then
800 A1 = C1 * A1
810 S0 = 1
820 Else
830 A1 = C2 * A1
840 ! restore array X
850 MAT X = X2
860 S0 = 0
870 If A1 >= A9 Then
880 DISP "**********SOLUTION FAILS************"
890 End If
900 End If
910 !
920 S2 = 1
930 For I = 1 To N
940 If Abs(D(I)) > T0(I) Then S2 = 0
950 Next I
960 !
970 UNTIL (S0 = 1) OR (S2 = 1)
980 !
990 If I1 > I9 Then LEAVE
1000 If A1 >= A9 Then LEAVE
1010 !
1020 UNTIL S2 = 1
1030 !
1040 CALL MyFx(X, N, F2)
1050 DISP "**********FINAL RESULTS************"
1060 DISP "Optimum at:"
1070 For I = 1 To N
1080 DISP "X(";I;")=";X(I) @ PAUSE
1090 Next I
1100 For I = 1 To N
1110 DISP "Delta X(";I;")=";D(I) @ PAUSE
1120 Next I
1130 DISP "Function value ="; F2 @ PAUSE
1140 DISP "Number of iterations = ";I0
1150 !
1160 ! END
1170 !
1180 SUB MyFx(X(), N, R)
1190 !R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1200 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1210 END SUB
1220 !
1230 SUB MyFxEx(N,X(),D(), l8, R)
1240 DIM X0(N)
1250 For I = 1 To N
1260 X0(I) = X(I) + l8 * D(I)
1270 Next I
1280 CALL MyFx(X0, N, R)
1290 END SUB
1300 !
1310 SUB LSrch(X(), N, l8, D(), R)
1320 M8 = 100
1330 T8 = 0.000001
1340 R = 1
1350 K = 0
1360 REPEAT
1370 K = K + 1
1380 If K > M8 Then
1390 R = 0
1400 LEAVE
1410 End If
1420 h = 0.01 * (1 + Abs(l8))
1430 CALL MyFxEx(N, X, D, l8, F0)
1440 CALL MyFxEx(N, X, D, l8 + h, F1)
1450 CALL MyFxEx(N, X, D, l8 - h, F2)
1460 D1 = (F1 - F2) / 2 / h
1470 D2 = (F1 - 2 * F0 + F2) / h ^ 2
1480 If D2 = 0 Then LEAVE
1490 D0 = D1 / D2
1500 l8 = l8 - D0
1510 UNTIL Abs(D0) < T8
1520 END SUB
1530 !
1540 SUB FDrv(N, X(), I0, R)
1550 X0 = X(I0)
1560 h = 0.01 * (1 + Abs(X0))
1570 X(I0) = X0 + h
1580 CALL MyFx(X, N, F1)
1590 X(I0) = X0 - h
1600 CALL MyFx(X, N, F2)
1610 X(I0) = X0
1620 R = (F1 - F2) / 2 / h
1630 END SUB
1640 !
1650 SUB SDrv(N, X(), I0, J0, R)
1660 ! calculate second derivative?
1670 If I0 = J0 Then
1680 CALL MyFx(X, N, F0)
1690 X0 = X(I0)
1700 H1 = 0.01 * (1 + Abs(X0))
1710 X(I0) = X0 + H1
1720 CALL MyFx(X, N, F1)
1730 X(I0) = X0 - H1
1740 CALL MyFx(X, N, F2)
1750 X(I0) = X0
1760 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1770 Else
1780 X0 = X(I0)
1790 Y0 = X(J0)
1800 H1 = 0.01 * (1 + Abs(X0))
1810 H2 = 0.01 * (1 + Abs(Y0))
1820 ! calculate F3
1830 X(I0) = X0 + H1
1840 X(J0) = Y0 + H2
1850 CALL MyFx(X, N, F3)
1860 ! calculate F4
1870 X(I0) = X0 - H1
1880 X(J0) = Y0 - H2
1890 CALL MyFx(X, N, F4)
1900 ! calculate F5
1910 X(I0) = X0 + H1
1920 X(J0) = Y0 - H2
1930 CALL MyFx(X, N, F5)
1940 ! calculate F6
1950 X(I0) = X0 - H1
1960 X(J0) = Y0 + H2
1970 CALL MyFx(X, N, F6)
1980 X(I0) = X0
1990 X(J0) = Y0
2000 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
2010 End If
2020 END SUB
2030 !
2040 SUB GetDrv1(N, X(), D1())
2050 For I = 1 To N
2060 CALL FDrv(N, X, I, D1(I))
2070 Next I
2080 END SUB
2090 !
2100 SUB GetDrv2(N, X(), D2(,))
2110 For I = 1 To N
2120 For J = 1 To N
2130 CALL SDrv(N, X, I, J, D2(I, J))
2140 Next J
2150 Next I
2160 END SUB
The above listing uses the JPCROM to employ more structured constructs for loops and decision-making. The above listing uses the REPEAT-UNTIL and LOOP-ENDLOOP loops and also the IF-THEN-ELSE-ENDIF enhanced version of the basic IF statement. These program flow control statements limits or eliminates using the GOTO statements and make the programs easier to read and update. Here is the version of the listing that is indented and void of line numbers. I am including it here since it is easier to read than the above listing although it will not run on either the HP-71B handheld computer or the EMU71 emulator..   
! Marquardt's Method for Optimization
!
! USES JPCROM TO SUPPORT ADVANCED LOOPS AND 
! DECISION-MAKING CONSTRUCTS
!
! USES MATH ROM TO SUPPOT MAT STATEMENTS
N = 2
A9 = 1E+100 ! Alpha max
!
DIM X(N), X2(N), G(N)
DIM T0(N), D(N), J(N, N)
DIM O1(N,N), O2(N,N)
!
DISP "MARQUARDT OPTIMIZATION METHOD"
INPUT "Enter gradient norm tolerance value? "; E
INPUT "Enter maximum number iterations? "; I9
INPUT "Enter value for Alpha factor? "; A1
INPUT "Enter value for C1 factor? "; C1
INPUT "Enter value for C2 factor? "; C2
For I = 1 To N
  DISP "Enter guess for X(";I;")";
  INPUT X(I)
  DISP "Enter tolerance for X(";I;")";
  INPUT T0(I)  
Next I
!
CALL MyFx(X, N, F2)
!    
I0 = 0
REPEAT
  I0 = I0 + 1
  If I0 > I9 Then
    DISP "Reached maximum iterations limit"
    LEAVE
  End If
  !
  ! copy last Fx value
  F1 = F2
  !
  CALL GetDrv1(N, X, G)
  !
  ! test if gradient is shallow enough
  N0 = 0
  For I = 1 To N
    N0 = N0 + G(I)^2
  Next I  
  N0 = Sqr(N0)
  If N0 < E Then LEAVE
  !
  ! copy vector X
  MAT X2 = X
  ! get matrix J
  CALL GetDrv2(N, X, J)
  ! create identity matrix
  MAT O1 = IDN(N,N)
  ! Initialize inner loop
  I1 = 0
  !
  REPEAT
    I1 = I1 + 1
    If I1 > I9 Then LEAVE
    !
    MAT O2 = (A1) * O1  
    MAT J = J + O2 
    MAT J = INV(J)
    MAT G = J * G
    !
    MAT D = G
    MAT X = X - D
    !
    CALL MyFx(X, N, F2)
    !
		DISP "F = ";F2;" ";
		For I = 1 To N
		  DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
		Next I
		DISP
    !
    If F2 < F1 Then
      A1 = C1 * A1
      S0 = 1
    Else
      A1 = C2 * A1
      ! restore array X
      MAT X = X2 
      S0 = 0
      If A1 >= A9 Then 
        DISP "**********SOLUTION FAILS************"
      End If  
    End If
    !
    S2 = 1
    For I = 1 To N
      If Abs(D(I)) > T0(I) Then S2 = 0
    Next I
    !
  UNTIL (S0 = 1) OR (S2 = 1)
  !
  If I1 > I9 Then LEAVE
  If A1 >= A9 Then LEAVE
  !
UNTIL S2 = 1
!
CALL MyFx(X, N, F2)
DISP "**********FINAL RESULTS************"
DISP "Optimum at:"
For I = 1 To N
  DISP "X(";I;")=";X(I) @ PAUSE
Next I
For I = 1 To N
  DISP "Delta X(";I;")=";D(I) @ PAUSE
Next I
DISP "Function value ="; F2 @ PAUSE
DISP "Number of iterations = ";I0
!  
! END
!
SUB MyFx(X(), N, R)
  !R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
  R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
END SUB
!
SUB MyFxEx(N,X(),D(), l8, R)
  DIM X0(N)  
  For I = 1 To N
    X0(I) = X(I) + l8 * D(I)
  Next I  
  CALL MyFx(X0, N, R)
END SUB
!
SUB LSrch(X(), N, l8, D(), R) 
  M8 = 100
  T8 = 0.000001
  R = 1  
  K = 0
  REPEAT
    K = K + 1
    If K > M8 Then
      R = 0
      LEAVE
    End If
    h = 0.01 * (1 + Abs(l8))
    CALL MyFxEx(N, X, D, l8, F0)
    CALL MyFxEx(N, X, D, l8 + h, F1)
    CALL MyFxEx(N, X, D, l8 - h, F2)
    D1 = (F1 - F2) / 2 / h
    D2 = (F1 - 2 * F0 + F2) / h ^ 2
    If D2 = 0 Then LEAVE
    D0 = D1 / D2
    l8 = l8 - D0
  UNTIL Abs(D0) < T8
END SUB
!
SUB FDrv(N, X(), I0, R) 
  X0 = X(I0)
  h = 0.01 * (1 + Abs(X0))
  X(I0) = X0 + h
  CALL MyFx(X, N, F1)
  X(I0) = X0 - h
  CALL MyFx(X, N, F2)
  X(I0) = X0
  R = (F1 - F2) / 2 / h
END SUB
!
SUB SDrv(N, X(), I0, J0, R) 
  ! calculate second derivative?
  If I0 = J0 Then
    CALL MyFx(X, N, F0)
    X0 = X(I0)
    H1 = 0.01 * (1 + Abs(X0))
    X(I0) = X0 + H1
    CALL MyFx(X, N, F1)
    X(I0) = X0 - H1
    CALL MyFx(X, N, F2)
    X(I0) = X0
    R = (F1 - 2 * F0 + F2) / H1 ^ 2
  Else
    X0 = X(I0)
    Y0 = X(J0)
    H1 = 0.01 * (1 + Abs(X0))
    H2 = 0.01 * (1 + Abs(Y0))
    ! calculate F3
    X(I0) = X0 + H1
    X(J0) = Y0 + H2
    CALL MyFx(X, N, F3)
    ! calculate F4
    X(I0) = X0 - H1
    X(J0) = Y0 - H2
    CALL MyFx(X, N, F4)
    ! calculate F5
    X(I0) = X0 + H1
    X(J0) = Y0 - H2
    CALL MyFx(X, N, F5)
    ! calculate F6
    X(I0) = X0 - H1
    X(J0) = Y0 + H2
    CALL MyFx(X, N, F6)
    X(I0) = X0
    X(J0) = Y0
    R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
  End If
END SUB
!
SUB GetDrv1(N, X(), D1())  
  For I = 1 To N
    CALL FDrv(N, X, I, D1(I))
  Next I
END SUB
!
SUB GetDrv2(N, X(), D2(,))
  For I = 1 To N 
    For J = 1 To N 
      CALL SDrv(N, X, I, J, D2(I, J))
    Next J
  Next I
END SUB

Here is a version of the first listing that DOES NOT rely on the JPCROM. This version should run on the basic HP-71B with the Math module.

10 ! Marquardt's Method for Optimization
20 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 N = 2
80 A9 = 1E+100 ! Alpha max
90 !
100 DIM X(N), X2(N), G(N)
110 DIM T0(N), D(N), J(N, N)
120 DIM O1(N,N), O2(N,N)
130 !
140 DISP "MARQUARDT OPTIMIZATION METHOD"
150 INPUT "Enter gradient norm tolerance value? "; E
160 INPUT "Enter maximum number iterations? "; I9
170 INPUT "Enter value for Alpha factor? "; A1
180 INPUT "Enter value for C1 factor? "; C1
190 INPUT "Enter value for C2 factor? "; C2
200 For I = 1 To N
210 DISP "Enter guess for X(";I;")";
220 INPUT X(I)
230 DISP "Enter tolerance for X(";I;")";
240 INPUT T0(I)
250 Next I
260 !
270 CALL MyFx(X, N, F2)
280 !
290 I0 = 0
300 !REPEAT
310 I0 = I0 + 1
320 If I0 <= I9 Then 350
330 DISP "Reached maximum iterations limit"
340 GOTO 1040
350 !End If
360 !
370 ! copy last Fx value
380 F1 = F2
390 !
400 CALL GetDrv1(N, X, G)
410 !
420 ! test if gradient is shallow enough
430 N0 = 0
440 For I = 1 To N
450 N0 = N0 + G(I)^2
460 Next I
470 N0 = Sqr(N0)
480 If N0 < E Then 1040
490 !
500 ! copy vector X
510 MAT X2 = X
520 ! get matrix J
530 CALL GetDrv2(N, X, J)
540 ! create identity matrix
550 MAT O1 = IDN(N,N)
560 ! Initialize inner loop
570 I1 = 0
580 !
590 !REPEAT
600 I1 = I1 + 1
610 If I1 > I9 Then 990
620 !
630 MAT O2 = (A1) * O1
640 MAT J = J + O2
650 MAT J = INV(J)
660 MAT G = J * G
670 !
680 MAT D = G
690 MAT X = X - D
700 !
710 CALL MyFx(X, N, F2)
720 !
730 DISP "F = ";F2;" ";
740 For I = 1 To N
750   DISP "X=(";I;")=";X(I);" Delta(";I;")=";D(I);" ";
760 Next I
770 DISP
780 !
790 If F2 >= F1 Then 820
800 A1 = C1 * A1
810 S0 = 1
815 GOTO 900
820 !Else
830 A1 = C2 * A1
840 ! restore array X
850 MAT X = X2
860 S0 = 0
870 If A1 >= A9 Then DISP "**********SOLUTION FAILS************"
900 !End If
910 !
920 S2 = 1
930 For I = 1 To N
940 If Abs(D(I)) > T0(I) Then S2 = 0
950 Next I
960 !
970 If (S0 = 1) OR (S2 = 1) Then 590
980 !
990 If I1 > I9 Then 1040
1000 If A1 >= A9 Then 1040
1010 !
1020 If S2 = 1 THEN 300
1030 !
1040 CALL MyFx(X, N, F2)
1050 DISP "**********FINAL RESULTS************"
1060 DISP "Optimum at:"
1070 For I = 1 To N
1080 DISP "X(";I;")=";X(I) @ PAUSE
1090 Next I
1100 For I = 1 To N
1110 DISP "Delta X(";I;")=";D(I) @ PAUSE
1120 Next I
1130 DISP "Function value ="; F2 @ PAUSE
1140 DISP "Number of iterations = ";I0
1150 !
1160 ! END
1170 !
1180 SUB MyFx(X(), N, R)
1190 !R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1200 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1210 END SUB
1220 !
1230 SUB MyFxEx(N,X(),D(), l8, R)
1240 DIM X0(N)
1250 For I = 1 To N
1260 X0(I) = X(I) + l8 * D(I)
1270 Next I
1280 CALL MyFx(X0, N, R)
1290 END SUB
1300 !
1310 SUB LSrch(X(), N, l8, D(), R)
1320 M8 = 100
1330 T8 = 0.000001
1340 R = 1
1350 K = 0
1360 !REPEAT
1370 K = K + 1
1380 If K > M8 Then
1390 R = 0
1400 GOTO 1520
1410 End If
1420 h = 0.01 * (1 + Abs(l8))
1430 CALL MyFxEx(N, X, D, l8, F0)
1440 CALL MyFxEx(N, X, D, l8 + h, F1)
1450 CALL MyFxEx(N, X, D, l8 - h, F2)
1460 D1 = (F1 - F2) / 2 / h
1470 D2 = (F1 - 2 * F0 + F2) / h ^ 2
1480 If D2 = 0 Then 1520
1490 D0 = D1 / D2
1500 l8 = l8 - D0
1510 If Abs(D0) < T8 Then 1360
1520 END SUB
1530 !
1540 SUB FDrv(N, X(), I0, R)
1550 X0 = X(I0)
1560 h = 0.01 * (1 + Abs(X0))
1570 X(I0) = X0 + h
1580 CALL MyFx(X, N, F1)
1590 X(I0) = X0 - h
1600 CALL MyFx(X, N, F2)
1610 X(I0) = X0
1620 R = (F1 - F2) / 2 / h
1630 END SUB
1640 !
1650 SUB SDrv(N, X(), I0, J0, R)
1660 ! calculate second derivative?
1670 If I0 <> J0 Then 1770
1680 CALL MyFx(X, N, F0)
1690 X0 = X(I0)
1700 H1 = 0.01 * (1 + Abs(X0))
1710 X(I0) = X0 + H1
1720 CALL MyFx(X, N, F1)
1730 X(I0) = X0 - H1
1740 CALL MyFx(X, N, F2)
1750 X(I0) = X0
1760 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1765 GOTO 2010
1770 !Else
1780 X0 = X(I0)
1790 Y0 = X(J0)
1800 H1 = 0.01 * (1 + Abs(X0))
1810 H2 = 0.01 * (1 + Abs(Y0))
1820 ! calculate F3
1830 X(I0) = X0 + H1
1840 X(J0) = Y0 + H2
1850 CALL MyFx(X, N, F3)
1860 ! calculate F4
1870 X(I0) = X0 - H1
1880 X(J0) = Y0 - H2
1890 CALL MyFx(X, N, F4)
1900 ! calculate F5
1910 X(I0) = X0 + H1
1920 X(J0) = Y0 - H2
1930 CALL MyFx(X, N, F5)
1940 ! calculate F6
1950 X(I0) = X0 - H1
1960 X(J0) = Y0 + H2
1970 CALL MyFx(X, N, F6)
1980 X(I0) = X0
1990 X(J0) = Y0
2000 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
2010 !End If
2020 END SUB
2030 !
2040 SUB GetDrv1(N, X(), D1())
2050 For I = 1 To N
2060 CALL FDrv(N, X, I, D1(I))
2070 Next I
2080 END SUB
2090 !
2100 SUB GetDrv2(N, X(), D2(,))
2110 For I = 1 To N
2120 For J = 1 To N
2130 CALL SDrv(N, X, I, J, D2(I, J))
2140 Next J
2150 Next I
2160 END SUB

BACK

Copyright (c) Namir Shammas. All rights reserved.