HP-71B Program to Find a Function Minimum

Using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) Method

by Namir Shammas

The following program uses the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method to find the minimum of a function. This method is a quasi-Newton method. That is, the BFGS method is based on Newton's method but performs different calculations to update the guess refinements.

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

1. The maximum number of iterations.

2. The tolerance for the minimized function,

3. The tolerance for the gradient.

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

5. The tolerance for the guess refinement for each variable.

The program displays intermediate values for the function and the variables. 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 ! Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)
20 !
30 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
40 ! DECISION-MAKING CONSTRUCTS
50 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 !
110 N = 2
120 !
130 DIM T0(N), X(N), G1(N,1), G2(N,1), G(N,1)
140 DIM D(N,1), S(N,1), B(N, N), M(N, N)
150 DIM N1(N, N), M1(N, N), M2(N, N), M3(N, N)
160 DIM M4(N, N), M5(N, N), M6(N, N), M7(N, N)
170 DIM M8(N, N), M9(N, N), Z1(N, N), Z2(N, N)
180 !
190 DISP "Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)"
200 INPUT "Enter maximum number of iterations? "; I9
210 INPUT "Enter function tolerance? "; E1
220 INPUT "Enter gradient tolerance? "; E2
230 N = N
240 For I = 1 To N
250 DISP "Enter value for X(";I;")";
260 INPUT X(I)
270 DISP "Enter tolerance value for X(";I;")";
280 INPUT T0(I)
290 Next I
300 !
310 ! set matrix B as an indentity matrix
320 MAT B = IDN(N,N)
330 !
340 I0 = 0
350 ! calculate initial gradient
360 CALL GetDrv1(N, X, G1)
370 !
380 ! start main loop
390 LOOP
400 !
410 I0 = I0 + 1
420 If I0 > I9 Then
430 DISP "Reached iteration limits"
440 LEAVE
450 End If
460 !
470 ! calculate vector S() and reverse its sign
480 MAT S = B * G1
490 MAT S = (-1) * S
500 !
510 CALL CalcNorm(S, N, N0)
520 MAT S = (1/N0) * S
530 !
540 L8 = 0.1
550 CALL LSrch(X, N, L8, S, 0.1, 0.00001, R)
560 ! calculate optimum X()
570 MAT D = (L8) * S
580 For I = 1 To N
590 X(I) = X(I) + D(I,1)
600 Next I
610 !
620 ! get new gradient
630 CALL GetDrv1(N, X, G2)
640 MAT G = G2 - G1
650 MAT G1 = G2
660 !
670 ! test for convergence
680 S0 = 1
690 For I = 1 To N
700 If Abs(D(I,1)) > T0(I) Then
710 S0 = 0
720 I=N
730 End If
740 Next I
750 !
760 If S0 = 1 Then
770 DISP "Exit due to values of L8 * S()"
780 DISP "Lamda="; L8
790 DISP "Array S is:"
800 For I = 1 To N
810 DISP "S(";I;")=";S(I,1);" ";
820 Next I
830 DISP
840 LEAVE
850 End If
860 !
870 ! exit if gradient is low
880 CALL CalcNorm(G, N, N0)
890 If N0 < E2 Then
900 DISP "Exit due to G Norm"
910 DISP "Norm=";N0
920 LEAVE
930 End If
940 !
950 !-------------------------------------------------
960 ! Start elaborate process to upgare matrix B
970 !
980 ! M1 = G as column matrix
990 MAT M1 = G
1000 ! M2 = D as column matrix
1010 MAT M2 = D
1020 ! M3 = G^T
1030 MAT M3 = TRN(G)
1040 ! M4 = D^T
1050 MAT M4 = TRN(D)
1060 ! M5 = D D^T
1070 MAT M5 = M2 * M4
1080 ! M6 = D^T G
1090 MAT M6 = M4 * M1
1100 ! M7 = D G^T
1110 MAT M7 = M2 * M3
1120 ! M8 = G D^T
1130 MAT M8 = M1 * M4
1140 !objML.MatMultMat M3, B, M9
1150 MAT M9 = M3 * B
1160 ! M9 = G^T [B] G
1170 MAT M9 = M1 * M9
1180 ! Z1 = D G^T [B]
1190 MAT Z1 = M7 * B
1200 ! Z2 = [B] G D^T
1210 MAT Z2 = B * M8
1220 !
1230 T1 = M6(1, 1) ! D^T G
1240 T2 = (1 + M9(1, 1) / T1) / T1
1250 MAT M5 = (T2) * M5
1260 MAT Z1 = (1/T1) * Z1
1270 MAT Z2 = (1/T2) * Z2
1280 MAT B = B + M5
1290 MAT B = B - Z1
1300 MAT B = B - Z2
1310 !
1320 CALL MyFx(X, N, F)
1330 DISP "F = ";F;" ";
1340 For I = 1 To N
1350 DISP "X(";I;")=";X(I);" ";
1360 DISP "G(";I;")=";G(I,1);" ";
1370 Next I
1380 DISP
1390 !
1400 END LOOP
1410 !
1420 DISP "********FINAL RESULTS**********"
1430 For I = 1 To N
1440 DISP "X(";I;")=";X(I) @ PAUSE
1450 Next I
1460 For I = 1 To N
1470 DISP "Gradient(";I;")=";G(I,1) @ PAUSE
1480 Next I
1490 CALL MyFx(X, N, F)
1500 DISP "Function value = ";F @ PAUSE
1510 DISP "Iterations = ";I0
1520 !
1530 ! END
1540 !
1550 SUB CalcNorm(X(,), N, N0)
1560 N0 = 0
1570 For I = 1 To N
1580 N0 = N0 + X(I,1)^2
1590 Next I
1600 N0 = Sqr(N0)
1610 END SUB
1620 !
1630 SUB MyFx(X(), N, R)
1640 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1650 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1660 End SUB
1670 !
1680 SUB FDrv(N, X(), I, R)
1690 X0 = X(I)
1700 h = 0.01 * (1 + Abs(X0))
1710 X(I) = X0 + h
1720 CALL MyFx(X, N, F1)
1730 X(I) = X0 - h
1740 CALL MyFx(X, N, F2)
1750 X(I) = X0
1760 R = (F1 - F2) / 2 / h
1770 End SUB
1780 !
1790 Sub GetDrv1(N, X(), D1(,))
1800 For I = 1 To N
1810 CALL FDrv(N, X, I, D1(I,1))
1820 Next I
1830 End Sub
1840 !
1850 SUB MyFxEx(N, X(), D(,), L8, R)
1860 DIM X2(N)
1870 !
1880 For I = 1 To N
1890 X2(I) = X(I) + L8 * D(I,1)
1900 Next I
1910 CALL MyFx(X2, N, R)
1920 End SUB
1930 !
1940 SUB LSrch(X(), N, L8, D(,), I9, I0, R)
1950 CALL MyFxEx(N, X, D, L8, F1)
1960 REPEAT
1970 CALL MyFxEx(N, X, D, L8 + I9, F2)
1980 If F2 < F1 Then
1990 F1 = F2
2000 L8 = L8 + I9
2010 Else
2020 CALL MyFxEx(N, X, D, L8 - I9, F2)
2030 If F2 < F1 Then
2040 F1 = F2
2050 L8 = L8 - I9
2060 Else
2070 ! reduce search step size
2080 I9 = I9 / 10
2090 End If
2100 End If
2110 UNTIL I9 < I0
2120 R = 1
2130 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..   
! Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)
!
! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
! DECISION-MAKING CONSTRUCTS
!
! USES MATH ROM TO SUPPOT MAT STATEMENTS
!
N = 2
!
DIM T0(N), X(N), G1(N,1), G2(N,1), G(N,1)
DIM D(N,1), S(N,1), B(N, N), M(N, N)
DIM N1(N, N), M1(N, N), M2(N, N), M3(N, N)
DIM M4(N, N), M5(N, N), M6(N, N), M7(N, N)
DIM M8(N, N), M9(N, N), Z1(N, N), Z2(N, N)
!
DISP "Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)"
INPUT "Enter maximum number of iterations? "; I9
INPUT "Enter function tolerance? "; E1
INPUT "Enter gradient tolerance? "; E2
N = N
For I = 1 To N
  DISP "Enter value for X(";I;")";
  INPUT X(I)
  DISP "Enter tolerance value for X(";I;")";
  INPUT T0(I)
Next I
!
! set matrix B as an indentity matrix
MAT B = IDN(N)
!
I0 = 0
! calculate initial gradient
CALL GetDrv1(N, X, G1)
!
! start main loop
LOOP
  !
  I0 = I0 + 1
  If I0 > I9 Then
    DISP "Reached iteration limits"
    LEAVE
  End If
  !
  ! calculate vector S() and reverse its sign
  MAT S = B * G1
  MAT S = (-1) * S
  !
  CALL CalcNorm(S, N, N0)
  MAT S = (1/N0) * S
  !
  L8 = 0.1
  CALL LSrch(X, N, L8, S, 0.1, 0.00001, R)
  ! calculate optimum X()
  MAT D = (L8) * S
  For I = 1 To N
    X(I) = X(I) + D(I,1)
  Next I
  !
  ! get new gradient
  CALL GetDrv1(N, X, G2)
  MAT G = G2 - G1
  MAT G1 = G2
  !
  ! test for convergence
  S0 = 1
  For I = 1 To N
    If Abs(D(I,1)) > T0(I) Then
      S0 = 0
      I=N
    End If
  Next I
  !
  If S0 = 1 Then
    DISP "Exit due to values of L8 * S()"
    DISP "Lamda="; L8
    DISP "Array S is:"
    For I = 1 To N
      DISP "S(";I;")=";S(I,1);" ";
    Next I
    DISP
    LEAVE
  End If
  !
  ! exit if gradient is low
  CALL CalcNorm(G, N, N0)
  If N0 < E2 Then
    DISP "Exit due to G Norm"
    DISP "Norm=";N0
    LEAVE
  End If
  !
  !-------------------------------------------------
  ! Start elaborate process to upgare matrix B
  !
  ! M1 = G as column matrix
  MAT M1 = G
  ! M2 = D as column matrix
  MAT M2 = D
  ! M3 = G
  MAT M3 = TRN(G)
  ! M4 = D
  MAT M4 = TRN(D)
  ! M5 = D D
  MAT M5 = M2 * M4
  ! M6 = D G
  MAT M6 = M4 * M1
  ! M7 = D G
  MAT M7 = M2 * M3
  ! M8 = G D
  MAT M8 = M1 * M4
  !objML.MatMultMat M3, B, M9
  MAT M9 = M3 * B
  ! M9 = G [B] G
  MAT M9 = M1 * M9
  ! Z1 = D G [B]
  MAT Z1 = M7 * B
  ! Z2 = [B] G D
  MAT Z2 = B * M8
  !
  T1 = M6(1, 1) ! D G
  T2 = (1 + M9(1, 1) / T1) / T1
  MAT M5 = (T2) * M5
  MAT Z1 = (1/T1) * Z1
  MAT Z2 = (1/T2) * Z2
  MAT B = B + M5
  MAT B = B - Z1
  MAT B = B - Z2
  !
  CALL MyFx(X, N, F)
  DISP "F = ";F;" ";
  For I = 1 To N
    DISP "X(";I;")=";X(I);" ";
    DISP "G(";I;")=";G(I,1);" ";
  Next I
  DISP
  !
END LOOP
!
DISP "********FINAL RESULTS**********"
For I = 1 To N
  DISP "X(";I;")=";X(I) @ PAUSE
Next I
For I = 1 To N
  DISP "Gradient(";I;")=";G(I,1) @ PAUSE
Next I
CALL MyFx(X, N, F)
DISP "Function value = ";F @ PAUSE
DISP "Iterations = ";I0
!
! END
!
SUB CalcNorm(X(,), N, N0)
  N0 = 0
  For I = 1 To N
    N0 = N0 + X(I,1)^2
  Next I
  N0 = Sqr(N0)
END SUB
!
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 FDrv(N, X(), I, R)
  X0 = X(I)
  h = 0.01 * (1 + Abs(X0))
  X(I) = X0 + h
  CALL MyFx(X, N, F1)
  X(I) = X0 - h
  CALL MyFx(X, N, F2)
  X(I) = X0
  R = (F1 - F2) / 2 / h
End SUB
!
Sub GetDrv1(N, X(), D1(,))
  For I = 1 To N
    CALL FDrv(N, X, I, D1(I,1))
  Next I
End Sub
!
SUB MyFxEx(N, X(), D(,), L8, R)
  DIM X2(N)
  !
  For I = 1 To N
    X2(I) = X(I) + L8 * D(I,1)
  Next I
  CALL MyFx(X2, N, R)
End SUB
!
SUB LSrch(X(), N, L8, D(,), I9, I0, R)
  CALL MyFxEx(N, X, D, L8, F1)
  REPEAT
    CALL MyFxEx(N, X, D, L8 + I9, F2)
    If F2 < F1 Then
      F1 = F2
      L8 = L8 + I9
    Else
      CALL MyFxEx(N, X, D, L8 - I9, F2)
      If F2 < F1 Then
        F1 = F2
        L8 = L8 - I9
      Else
        ! reduce search step size
        I9 = I9 / 10
      End If
    End If
  UNTIL I9 < I0
  R = 1
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 ! Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)
20 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 !
110 N = 2
120 !
130 DIM T0(N), X(N), G1(N,1), G2(N,1), G(N,1)
140 DIM D(N,1), S(N,1), B(N, N), M(N, N)
150 DIM N1(N, N), M1(N, N), M2(N, N), M3(N, N)
160 DIM M4(N, N), M5(N, N), M6(N, N), M7(N, N)
170 DIM M8(N, N), M9(N, N), Z1(N, N), Z2(N, N)
180 !
190 DISP "Broyden-Fletcher-Goldfarb-Shanno Method (BFGS)"
200 INPUT "Enter maximum number of iterations? "; I9
210 INPUT "Enter function tolerance? "; E1
220 INPUT "Enter gradient tolerance? "; E2
230 N = N
240 For I = 1 To N
250 DISP "Enter value for X(";I;")";
260 INPUT X(I)
270 DISP "Enter tolerance value for X(";I;")";
280 INPUT T0(I)
290 Next I
300 !
310 ! set matrix B as an indentity matrix
320 MAT B = IDN(N,N)
330 !
340 I0 = 0
350 ! calculate initial gradient
360 CALL GetDrv1(N, X, G1)
370 !
380 ! start main loop
390 ! LOOP *******************
400 !
410 I0 = I0 + 1
420 If I0 <= I9 Then 450
430 DISP "Reached iteration limits"
440 GOTO 1420
450 !End If
460 !
470 ! calculate vector S() and reverse its sign
480 MAT S = B * G1
490 MAT S = (-1) * S
500 !
510 CALL CalcNorm(S, N, N0)
520 MAT S = (1/N0) * S
530 !
540 L8 = 0.1
550 CALL LSrch(X, N, L8, S, 0.1, 0.00001, R)
560 ! calculate optimum X()
570 MAT D = (L8) * S
580 For I = 1 To N
590 X(I) = X(I) + D(I,1)
600 Next I
610 !
620 ! get new gradient
630 CALL GetDrv1(N, X, G2)
640 MAT G = G2 - G1
650 MAT G1 = G2
660 !
670 ! test for convergence
680 S0 = 1
690 For I = 1 To N
700 If Abs(D(I,1)) <= T0(I) Then 730
710 S0 = 0
720 I=N
730 !End If
740 Next I
750 !
760 If S0 = 0 Then 850
770 DISP "Exit due to values of L8 * S()"
780 DISP "Lamda="; L8
790 DISP "Array S is:"
800 For I = 1 To N
810 DISP "S(";I;")=";S(I,1);" ";
820 Next I
830 DISP
840 GOTO 1420
850 !End If
860 !
870 ! exit if gradient is low
880 CALL CalcNorm(G, N, N0)
890 If N0 >= E2 Then 930
900 DISP "Exit due to G Norm"
910 DISP "Norm=";N0
920 GOTO 1420
930 !End If
940 !
950 !-------------------------------------------------
960 ! Start elaborate process to upgare matrix B
970 !
980 ! M1 = G as column matrix
990 MAT M1 = G
1000 ! M2 = D as column matrix
1010 MAT M2 = D
1020 ! M3 = G
1030 MAT M3 = TRN(G)
1040 ! M4 = D
1050 MAT M4 = TRN(D)
1060 ! M5 = D D
1070 MAT M5 = M2 * M4
1080 ! M6 = D G
1090 MAT M6 = M4 * M1
1100 ! M7 = D G
1110 MAT M7 = M2 * M3
1120 ! M8 = G D
1130 MAT M8 = M1 * M4
1140 !objML.MatMultMat M3, B, M9
1150 MAT M9 = M3 * B
1160 ! M9 = G [B] G
1170 MAT M9 = M1 * M9
1180 ! Z1 = D G [B]
1190 MAT Z1 = M7 * B
1200 ! Z2 = [B] G D
1210 MAT Z2 = B * M8
1220 !
1230 T1 = M6(1, 1) ! D G
1240 T2 = (1 + M9(1, 1) / T1) / T1
1250 MAT M5 = (T2) * M5
1260 MAT Z1 = (1/T1) * Z1
1270 MAT Z2 = (1/T2) * Z2
1280 MAT B = B + M5
1290 MAT B = B - Z1
1300 MAT B = B - Z2
1310 !
1320 CALL MyFx(X, N, F)
1330 DISP "F = ";F;" ";
1340 For I = 1 To N
1350 DISP "X(";I;")=";X(I);" ";
1360 DISP "G(";I;")=";G(I,1);" ";
1370 Next I
1380 DISP
1390 !
1400 GOTO 390 ! END LOOP *******************
1410 !
1420 DISP "********FINAL RESULTS**********"
1430 For I = 1 To N
1440 DISP "X(";I;")=";X(I) @ PAUSE
1450 Next I
1460 For I = 1 To N
1470 DISP "Gradient(";I;")=";G(I,1) @ PAUSE
1480 Next I
1490 CALL MyFx(X, N, F)
1500 DISP "Function value = ";F @ PAUSE
1510 DISP "Iterations = ";I0
1520 !
1530 ! END
1540 !
1550 SUB CalcNorm(X(,), N, N0)
1560 N0 = 0
1570 For I = 1 To N
1580 N0 = N0 + X(I,1)^2
1590 Next I
1600 N0 = Sqr(N0)
1610 END SUB
1620 !
1630 SUB MyFx(X(), N, R)
1640 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1650 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1660 End SUB
1670 !
1680 SUB FDrv(N, X(), I, R)
1690 X0 = X(I)
1700 h = 0.01 * (1 + Abs(X0))
1710 X(I) = X0 + h
1720 CALL MyFx(X, N, F1)
1730 X(I) = X0 - h
1740 CALL MyFx(X, N, F2)
1750 X(I) = X0
1760 R = (F1 - F2) / 2 / h
1770 End SUB
1780 !
1790 Sub GetDrv1(N, X(), D1(,))
1800 For I = 1 To N
1810 CALL FDrv(N, X, I, D1(I,1))
1820 Next I
1830 End Sub
1840 !
1850 SUB MyFxEx(N, X(), D(,), L8, R)
1860 DIM X2(N)
1870 !
1880 For I = 1 To N
1890 X2(I) = X(I) + L8 * D(I,1)
1900 Next I
1910 CALL MyFx(X2, N, R)
1920 End SUB
1930 !
1940 SUB LSrch(X(), N, L8, D(,), I9, I0, R)
1950 CALL MyFxEx(N, X, D, L8, F1)
1960 !REPEAT
1970 CALL MyFxEx(N, X, D, L8 + I9, F2)
1980 If F2 >= F1 Then 2010
1990 F1 = F2
2000 L8 = L8 + I9
2005 GOTO 2100
2010 !Else
2020 CALL MyFxEx(N, X, D, L8 - I9, F2)
2030 If F2 >= F1 Then 2060
2040 F1 = F2
2050 L8 = L8 - I9
2055 GOTO 2090
2060 !Else
2070 ! reduce search step size
2080 I9 = I9 / 10
2090 !End If
2100 !End If
2110 If I9 >= I0 Then 1960
2120 R = 1
2130 End SUB

BACK

Copyright (c) Namir Shammas. All rights reserved.