The following program uses the Davidson-Fletcher-Powell (DFP) method to find the minimum of a function. This method is a quasi-Newton method. That is, the DFP algorithm 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 ! Davidson-Fletcher-Powell Method (DFP)
20 !
30 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
40 ! DECISION-MAKING CONSTRUCTS
50 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 !
80 N=2
90 !
100 DIM T0(N), X(N), G1(N,1), G2(N,1)
110 DIM G(N,1), S(N,1), D(N,1), B(N, N)
120 DIM M(N, N),M9(N, N)
130 DIM M1(N, 1), M2(N, N), M3(N, N)
140 DIM M4(N, N), V1(N,1)
150 REM
160 DISP "Davidson-Fletcher-Powell Method (DFP)"
170 INPUT "Enter maximum number of iterations? "; I9
180 INPUT "Enter function tolerance? "; E1
190 INPUT "Enter gradient tolerance? "; E2
200 For I = 1 To N
210 DISP "Enter value for X(";I;")";
220 INPUT X(I)
230 DISP "Enter tolerance value for X(";I;")";
240 INPUT T0(I)
250 Next I
260 !
270 ! set matrix B as an indentity matrix
280 MAT B = IDN(N,N)
290 !
300 I0 = 0
310 ! calculate initial gradient
320 CALL GetDrv1(N, X, G1)
330 !
340 ! start main loop
350 LOOP
360 !
370 I0 = I0 + 1
380 If I0 > I9 Then
390 DISP "Reached iteration limits"
400 LEAVE
410 End If
420 !
430 ! calcuaLte vector S() and reverse its sign
440 MAT S = B * G1
450 MAT S = (-1) * S
460 ! normailize vector S
470 CALL CalcNorm(S, N, N0)
480 MAT S = (1/N0) * S
490 !
500 L8 = 0.1
510 CALL LSrch(X, N, L8, S, 0.1, 0.00001, R)
520 ! calculate optimum X() with the given L8
530 MAT D = (L8) * S
540 For I = 1 To N
550 X(I) = X(I) + D(I,1)
560 Next I
570 !
580 ! get new gradient
590 CALL GetDrv1(N, X, G2)
600 ! calculate vector G() and shift G2() into G1()
610 MAT G = G2 - G1
620 MAT G1 = G2
630 !
640 ! test for convergence
650 S0 = 1
660 For I = 1 To N
670 If Abs(D(I,1)) > T0(I) Then
680 S0 = 0
690 I = N
700 End If
710 Next I
720 !
730 If S0 = 1 Then
740 DISP "Exit due to values of L8 * S()"
750 DISP "Lamda="; L8
760 DISP "Array S is:"
770 For I = 1 To N
780 DISP "S(";I;")=";S(I,1);" ";
790 Next I
800 DISP
810 LEAVE
820 End If
830 !
840 ! exit if gradient is low
850 CALL CalcNorm(G, N, N0)
860 If N0 < E2 Then
870 DISP "Exit due to G Norm"
880 DISP "Norm=";N0
890 LEAVE
900 End If
910 !
920 !-------------------------------------------------
930 ! Start elaborate process to update matrix B
940 !
950 ! First calculate matrix M (stored as M)
960 ! M1 = S in matrix form
970 MAT M1 = S
980 ! M2 = S^T in matrix form
990 MAT M2 = TRN(S)
1000 ! M3 = G in matrix form
1010 MAT M3 = G
1020 ! M = S * S^T
1030 MAT M = M1 * M2
1040 ! M4 = S^T * G
1050 MAT M4 = M2 * M3
1060 ! calculate natrix M
1070 MAT M = (L8 / M4(1,1)) * M
1080 ! Calculate matrix N (stored as M9)
1090 ! V1 = {B] G
1100 MAT V1 = B * G
1110 ! M1 = V1 in column matrix form ([B] G)
1120 MAT M1 = V1
1130 ! M2 = V1 in I0 matrix form ([B] G)^T
1140 MAT M2 = TRN(V1)
1150 ! M9 = ([B] G) * ([B] G)^T
1160 MAT M9 = M1 * M2
1170 ! M1 = G^T
1180 MAT M1 = TRN(G)
1190 ! M2 = G
1200 MAT M2 = G
1210 ! M3 = G^T [B]
1220 MAT M3 = M1 * B
1230 ! M4 = G^T [B] G
1240 MAT M4 = M3 * M2
1250 ! calculate matrix N
1260 MAT M9 = (-1/M4(1,1)) * M9
1270 !
1280 ! B = B + M
1290 MAT B = B + M
1300 ! B = B + N
1310 MAT B = B + M9
1320 !
1330 CALL MyFx(X, N, F)
1340 DISP "F = ";F;" ";
1350 For I = 1 To N
1360 DISP "X(";I;")=";X(I);" ";
1370 DISP "G(";I;")=";G(I,1);" ";
1380 Next I
1390 DISP
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 For I = 1 To N
1880 X2(I) = X(I) + L8 * D(I,1)
1890 Next I
1900 CALL MyFx(X2, N, R)
1910 End SUB
1920 !
1930 SUB LSrch(X(), N, L8, D(,), I9, I0, R)
1940 CALL MyFxEx(N, X, D, L8, F1)
1950 REPEAT
1960 CALL MyFxEx(N, X, D, L8 + I9, F2)
1970 If F2 < F1 Then
1980 F1 = F2
1990 L8 = L8 + I9
2000 Else
2010 CALL MyFxEx(N, X, D, L8 - I9, F2)
2020 If F2 < F1 Then
2030 F1 = F2
2040 L8 = L8 - I9
2050 Else
2060 ! reduce search step size
2070 I9 = I9 / 10
2080 End If
2090 End If
2100 UNTIL I9 < I0
2110 R = 1
2120 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..
! Davidson-Fletcher-Powell Method (DFP)
!
! 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)
DIM G(N,1), S(N,1), D(N,1), B(N, N)
DIM M(N, N),M9(N, N)
DIM M1(N, 1), M2(N, N), M3(N, N)
DIM M4(N, N), V1(N,1)
DISP "Davidson-Fletcher-Powell Method (DFP)"
INPUT "Enter maximum number of iterations? "; I9
INPUT "Enter function tolerance? "; E1
INPUT "Enter gradient tolerance? "; E2
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,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
! normailize vector 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() with the given L8
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)
! calculate vector G() and shift G2() into G1()
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 update matrix B
!
! First calculate matrix M (stored as M)
! M1 = S in matrix form
MAT M1 = S
! M2 = S in matrix form
MAT M2 = TRN(S)
! M3 = G in matrix form
MAT M3 = G
! M = S * S
MAT M = M1 * M2
! M4 = S * G
MAT M4 = M2 * M3
! calculate natrix M
MAT M = (L8 / M4(1,1)) * M
! Calculate matrix N (stored as M9)
! V1 = {B] G
MAT V1 = B * G
! M1 = V1 in column matrix form ([B] G)
MAT M1 = V1
! M2 = V1 in I0 matrix form ([B] G)
MAT M2 = TRN(V1)
! M9 = ([B] G) * ([B] G)
MAT M9 = M1 * M2
! M1 = G
MAT M1 = TRN(G)
! M2 = G
MAT M2 = G
! M3 = G [B]
MAT M3 = M1 * B
! M4 = G [B] G
MAT M4 = M3 * M2
! calculate matrix N
MAT M9 = (-1/M4(1,1)) * M9
!
! B = B + M
MAT B = B + M
! B = B + N
MAT B = B + M9
!
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 ! Davidson-Fletcher-Powell Method (DFP)
20 !
60 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
70 !
80 N=2
90 !
100 DIM T0(N), X(N), G1(N,1), G2(N,1)
110 DIM G(N,1), S(N,1), D(N,1), B(N, N)
120 DIM M(N, N),M9(N, N)
130 DIM M1(N, 1), M2(N, N), M3(N, N)
140 DIM M4(N, N), V1(N,1)
150 REM
160 DISP "Davidson-Fletcher-Powell Method (DFP)"
170 INPUT "Enter maximum number of iterations? "; I9
180 INPUT "Enter function tolerance? "; E1
190 INPUT "Enter gradient tolerance? "; E2
200 For I = 1 To N
210 DISP "Enter value for X(";I;")";
220 INPUT X(I)
230 DISP "Enter tolerance value for X(";I;")";
240 INPUT T0(I)
250 Next I
260 !
270 ! set matrix B as an indentity matrix
280 MAT B = IDN(N,N)
290 !
300 I0 = 0
310 ! calculate initial gradient
320 CALL GetDrv1(N, X, G1)
330 !
340 ! start main loop
350 !LOOP
360 !
370 I0 = I0 + 1
380 If I0 <= I9 Then 410
390 DISP "Reached iteration limits"
400 GOTO 1420
410 !End If
420 !
430 ! calcuLate vector S() and reverse its sign
440 MAT S = B * G1
450 MAT S = (-1) * S
460 ! normailize vector S
470 CALL CalcNorm(S, N, N0)
480 MAT S = (1/N0) * S
490 !
500 L8 = 0.1
510 CALL LSrch(X, N, L8, S, 0.1, 0.00001, R)
520 ! calculate optimum X() with the given L8
530 MAT D = (L8) * S
540 For I = 1 To N
550 X(I) = X(I) + D(I,1)
560 Next I
570 !
580 ! get new gradient
590 CALL GetDrv1(N, X, G2)
600 ! calculate vector G() and shift G2() into G1()
610 MAT G = G2 - G1
620 MAT G1 = G2
630 !
640 ! test for convergence
650 S0 = 1
660 For I = 1 To N
670 If Abs(D(I,1)) <= T0(I) Then 700
680 S0 = 0
690 I = N
700 !End If
710 Next I
720 !
730 If S0 <> 1 Then 820
740 DISP "Exit due to values of L8 * S()"
750 DISP "Lamda="; L8
760 DISP "Array S is:"
770 For I = 1 To N
780 DISP "S(";I;")=";S(I,1);" ";
790 Next I
800 DISP
810 GOTO 1420
820 !End If
830 !
840 ! exit if gradient is low
850 CALL CalcNorm(G, N, N0)
860 If N0 >= E2 Then 900
870 DISP "Exit due to G Norm"
880 DISP "Norm=";N0
890 GOTO 1420
900 !End If
910 !
920 !-------------------------------------------------
930 ! Start elaborate process to update matrix B
940 !
950 ! First calculate matrix M (stored as M)
960 ! M1 = S in matrix form
970 MAT M1 = S
980 ! M2 = S in matrix form
990 MAT M2 = TRN(S)
1000 ! M3 = G in matrix form
1010 MAT M3 = G
1020 ! M = S * S
1030 MAT M = M1 * M2
1040 ! M4 = S * G
1050 MAT M4 = M2 * M3
1060 ! calculate natrix M
1070 MAT M = (L8 / M4(1,1)) * M
1080 ! Calculate matrix N (stored as M9)
1090 ! V1 = {B] G
1100 MAT V1 = B * G
1110 ! M1 = V1 in column matrix form ([B] G)
1120 MAT M1 = V1
1130 ! M2 = V1 in I0 matrix form ([B] G)
1140 MAT M2 = TRN(V1)
1150 ! M9 = ([B] G) * ([B] G)
1160 MAT M9 = M1 * M2
1170 ! M1 = G
1180 MAT M1 = TRN(G)
1190 ! M2 = G
1200 MAT M2 = G
1210 ! M3 = G [B]
1220 MAT M3 = M1 * B
1230 ! M4 = G [B] G
1240 MAT M4 = M3 * M2
1250 ! calculate matrix N
1260 MAT M9 = (-1/M4(1,1)) * M9
1270 !
1280 ! B = B + M
1290 MAT B = B + M
1300 ! B = B + N
1310 MAT B = B + M9
1320 !
1330 CALL MyFx(X, N, F)
1340 DISP "F = ";F;" ";
1350 For I = 1 To N
1360 DISP "X(";I;")=";X(I);" ";
1370 DISP "G(";I;")=";G(I,1);" ";
1380 Next I
1390 DISP
1400 GOTO 350 ! 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 For I = 1 To N
1880 X2(I) = X(I) + L8 * D(I,1)
1890 Next I
1900 CALL MyFx(X2, N, R)
1910 End SUB
1920 !
1930 SUB LSrch(X(), N, L8, D(,), I9, I0, R)
1940 CALL MyFxEx(N, X, D, L8, F1)
1950 !REPEAT
1960 CALL MyFxEx(N, X, D, L8 + I9, F2)
1970 If F2 >= F1 Then 2000
1980 F1 = F2
1990 L8 = L8 + I9
1995 GOTO 2090
2000 !Else
2010 CALL MyFxEx(N, X, D, L8 - I9, F2)
2020 If F2 >= F1 Then 2050
2030 F1 = F2
2040 L8 = L8 - I9
2045 GOTO 2080
2050 !Else
2060 ! reduce search step size
2070 I9 = I9 / 10
2080 End If
2090 End If
2100 If I9 >= I0 Then 1950
2110 R = 1
2120 End SUB
Copyright (c) Namir Shammas. All rights reserved.