The following program calculates the minimum point of a multi-variable function using Newton's enhanced method. This method is implemented some enhancements that make the Newton's method more efficient.
The program prompts you to enter for each variable:
1. The maximum number of iterations
2. The values for the initial set of variables
3. The values for the tolerances for each variable.
The program displays intermediate results to show the iteration progress and also the following final results:
1. The coordinates of the minimum value.
2. The fine step sizes for each variable.
3. The minimum function value.
4. 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 ! Enhanced Newton's Method for Optimization
20 !---------------------------------------------------------------------------
30 !
40 ! This program implements the Modified Newton Method that uses a line search
50 ! step to improve the refined guess. As a result, this method can solve the
60 ! Rosenbrook function in contrast with the pure Newton Method which cannot.
70 !
80 !---------------------------------------------------------------------------
90 !
100 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
110 ! DECISION-MAKING CONSTRUCTS
120 !
130 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
140 !
150 N = 2
160 !
170 DIM X(N), G(N)
180 DIM T0(N)
190 DIM D(N), J(N, N)
200 !
210 DISP "NEWTON'S OPTIMIZATIOM METHOD (version 2)"
220 INPUT "Enter maximum number iterations? "; I9
230 For I = 1 To N
240 DISP "Enter guess for X(";I;")";
250 INPUT X(I)
260 DISP "Enter tolerance for X(";I;")";
270 INPUT T0(I)
280 Next I
290 !
300 I0 = 0
310 REPEAT
320 I0 = I0 + 1
330 If I0 > I9 Then
340 DISP "Reached maximum iterations limit"
350 LEAVE
360 End If
370 !
380 CALL GetDrv1(N, X, G)
390 !
400 ! test if gradient is shallow enough
410 N0 = 0
420 For I = 1 To N
430 N0 = N0 + G(I)^2
440 Next I
450 N0 = Sqr(N0)
460 ! If N0 < E Then LEAVE
470 !
480 CALL GetDrv2(N, X, J)
490 !
500 MAT J = INV(J)
510 MAT G = J * G
520 !
530 MAT D = G
540 !
550 L8 = 0.1
560 CALL LSrch(X, N, L8, D, K0)
570 If K0 = 0 Then
580 DISP "Linear Search failed"
590 LEAVE
600 End If
610 !
620 For I = 1 To N
630 D(I) = L8 * D(I)
640 X(I) = X(I) + D(I)
650 Next I
660 !
670 S0 = 1
680 For I = 1 To N
690 If Abs(D(I)) > T0(I) Then
700 S0 = 0
710 I = N
720 End If
730 Next I
740 REM
750 CALL MyFx(X, N, F)
760 DISP "F = ";F;" ";
770 For I = 1 To N
780 DISP "X=(";I;")=";X(I);" ";
790 Next I
800 DISP
810 !
820 UNTIL S0 = 1
830 !
840 CALL MyFx(X, N, F)
850 DISP "**********FINAL RESULTS************"
860 DISP "Optimum at:"
870 For I = 1 To N
880 DISP "X(";I;")=";X(I) @ PAUSE
890 Next I
900 For I = 1 To N
910 DISP "Delta X(";I;")=";D(I) @ PAUSE
920 Next I
930 DISP "Function value ="; F @ PAUSE
940 DISP "Number of iterations = ";I0
950 !
960 SUB MyFx(X(), N, R)
970 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
980 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
990 END SUB
1000 !
1010 SUB FrstDrv(N, X(), N0, R)
1020 X0 = X(N0)
1030 h = 0.01 * (1 + Abs(X0))
1040 X(N0) = X0 + h
1050 CALL MyFx(X, N, F1)
1060 X(N0) = X0 - h
1070 CALL MyFx(X, N, F2)
1080 X(N0) = X0
1090 R = (F1 - F2) / 2 / h
1100 END SUB
1110 !
1120 SUB SecDerv(N, X(), N0, J0, R)
1130 ! calculate second derivative?
1140 If N0 = J0 Then
1150 CALL MyFx(X, N, F0)
1160 X0 = X(N0)
1170 H1 = 0.01 * (1 + Abs(X0))
1180 X(N0) = X0 + H1
1190 CALL MyFx(X, N, F1)
1200 X(N0) = X0 - H1
1210 CALL MyFx(X, N, F2)
1220 X(N0) = X0
1230 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1240 Else
1250 X0 = X(N0)
1260 Y0 = X(J0)
1270 H1 = 0.01 * (1 + Abs(X0))
1280 H2 = 0.01 * (1 + Abs(Y0))
1290 ! calculate F3
1300 X(N0) = X0 + H1
1310 X(J0) = Y0 + H2
1320 CALL MyFx(X, N, F3)
1330 ! calculate F4
1340 X(N0) = X0 - H1
1350 X(J0) = Y0 - H2
1360 CALL MyFx(X, N, F4)
1370 ! calculate F5
1380 X(N0) = X0 + H1
1390 X(J0) = Y0 - H2
1400 CALL MyFx(X, N, F5)
1410 ! calculate F6
1420 X(N0) = X0 - H1
1430 X(J0) = Y0 + H2
1440 CALL MyFx(X, N, F6)
1450 X(N0) = X0
1460 X(J0) = Y0
1470 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
1480 End If
1490 END SUB
1500 !
1510 SUB GetDrv1(N, X(), D1())
1520 For I = 1 To N
1530 CALL FrstDrv(N, X, I, D1(I))
1540 Next I
1550 END SUB
1560 !
1570 SUB GetDrv2(N, X(), D2(,))
1580 For I = 1 To N
1590 For J = 1 To N
1600 CALL SecDerv(N, X, I, J, D2(I, J))
1610 Next J
1620 Next I
1630 END SUB
1640 !
1650 SUB MyFxEx(N,X(),D(), L, R)
1660 DIM X2(N)
1670 For I = 1 To N
1680 X2(I) = X(I) + L * D(I)
1690 Next I
1700 CALL MyFx(X2, N, R)
1710 END SUB
1720 !
1730 SUB LSrch(X(), N, L, D(), R)
1740 M8 = 100
1750 T0 = 0.000001
1760 R = 1
1770 REM
1780 I0 = 0
1790 REPEAT
1800 I0 = I0 + 1
1810 If I0 > M8 Then
1820 R = 0
1830 LEAVE
1840 End If
1850 h = 0.01 * (1 + Abs(L))
1860 CALL MyFxEx(N, X, D, L, F0)
1870 CALL MyFxEx(N, X, D, L + h, F1)
1880 CALL MyFxEx(N, X, D, L - h, F2)
1890 D1 = (F1 - F2) / 2 / h
1900 D2 = (F1 - 2 * F0 + F2) / h ^ 2
1910 If D2 = 0 Then LEAVE
1920 D0 = D1 / D2
1930 L = L - D0
1940 UNTIL Abs(D0) < T0
1950 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..
! Enhanced Newton's Method for Optimization
!---------------------------------------------------------------------------
!
! This program implements the Modified Newton Method that uses a line search
! step to improve the refined guess. As a result, this method can solve the
! Rosenbrook function in contrast with the pure Newton Method which cannot.
!
!---------------------------------------------------------------------------
!
! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
! DECISION-MAKING CONSTRUCTS
!
! USES MATH ROM TO SUPPOT MAT STATEMENTS
!
N = 2
!
DIM X(N), G(N)
DIM T0(N)
DIM D(N), J(N, N)
!
DISP "NEWTON'S OPTIMIZATIOM METHOD (version 2)"
INPUT "Enter maximum number iterations? "; I9
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
!
I0 = 0
REPEAT
I0 = I0 + 1
If I0 > I9 Then
DISP "Reached maximum iterations limit"
LEAVE
End If
!
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
!
CALL GetDrv2(N, X, J)
!
MAT J = INV(J)
MAT G = J * G
!
MAT D = G
!
L8 = 0.1
CALL LSrch(X, N, L8, D, K0)
If K0 = 0 Then
DISP "Linear Search failed"
LEAVE
End If
!
For I = 1 To N
D(I) = L8 * D(I)
X(I) = X(I) + D(I)
Next I
!
S0 = 1
For I = 1 To N
If Abs(D(I)) > T0(I) Then
S0 = 0
I = N
End If
Next I
CALL MyFx(X, N, F)
DISP "F = ";F;" ";
For I = 1 To N
DISP "X=(";I;")=";X(I);" ";
Next I
DISP
!
UNTIL S0 = 1
!
CALL MyFx(X, N, F)
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 ="; F @ PAUSE
DISP "Number of iterations = ";I0
!
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 FrstDrv(N, X(), N0, R)
X0 = X(N0)
h = 0.01 * (1 + Abs(X0))
X(N0) = X0 + h
CALL MyFx(X, N, F1)
X(N0) = X0 - h
CALL MyFx(X, N, F2)
X(N0) = X0
R = (F1 - F2) / 2 / h
END SUB
!
SUB SecDerv(N, X(), N0, J0, R)
! calculate second derivative?
If N0 = J0 Then
CALL MyFx(X, N, F0)
X0 = X(N0)
H1 = 0.01 * (1 + Abs(X0))
X(N0) = X0 + H1
CALL MyFx(X, N, F1)
X(N0) = X0 - H1
CALL MyFx(X, N, F2)
X(N0) = X0
R = (F1 - 2 * F0 + F2) / H1 ^ 2
Else
X0 = X(N0)
Y0 = X(J0)
H1 = 0.01 * (1 + Abs(X0))
H2 = 0.01 * (1 + Abs(Y0))
! calculate F3
X(N0) = X0 + H1
X(J0) = Y0 + H2
CALL MyFx(X, N, F3)
! calculate F4
X(N0) = X0 - H1
X(J0) = Y0 - H2
CALL MyFx(X, N, F4)
! calculate F5
X(N0) = X0 + H1
X(J0) = Y0 - H2
CALL MyFx(X, N, F5)
! calculate F6
X(N0) = X0 - H1
X(J0) = Y0 + H2
CALL MyFx(X, N, F6)
X(N0) = 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 FrstDrv(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 SecDerv(N, X, I, J, D2(I, J))
Next J
Next I
END SUB
!
SUB MyFxEx(N,X(),D(), L, R)
DIM X2(N)
For I = 1 To N
X2(I) = X(I) + L * D(I)
Next I
CALL MyFx(X2, N, R)
END SUB
!
SUB LSrch(X(), N, L, D(), R)
M8 = 100
T0 = 0.000001
R = 1
I0 = 0
REPEAT
I0 = I0 + 1
If I0 > M8 Then
R = 0
LEAVE
End If
h = 0.01 * (1 + Abs(L))
CALL MyFxEx(N, X, D, L, F0)
CALL MyFxEx(N, X, D, L + h, F1)
CALL MyFxEx(N, X, D, L - h, F2)
D1 = (F1 - F2) / 2 / h
D2 = (F1 - 2 * F0 + F2) / h ^ 2
If D2 = 0 Then LEAVE
D0 = D1 / D2
L = L - D0
UNTIL Abs(D0) < T0
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 ! Enhanced Newton's Method for Optimization
20 !---------------------------------------------------------------------------
30 !
40 ! This program implements the Modified Newton Method that uses a line search
50 ! step to improve the refined guess. As a result, this method can solve the
60 ! Rosenbrook function in contrast with the pure Newton Method which cannot.
70 !
80 !---------------------------------------------------------------------------
90 !
130 ! USES MATH ROM TO SUPPOT MAT STATEMENTS
140 !
150 N = 2
160 !
170 DIM X(N), G(N)
180 DIM T0(N)
190 DIM D(N), J(N, N)
200 !
210 DISP "NEWTON'S OPTIMIZATIOM METHOD (version 2)"
220 INPUT "Enter maximum number iterations? "; I9
230 For I = 1 To N
240 DISP "Enter guess for X(";I;")";
250 INPUT X(I)
260 DISP "Enter tolerance for X(";I;")";
270 INPUT T0(I)
280 Next I
290 !
300 I0 = 0
310 !REPEAT
320 I0 = I0 + 1
330 If I0 <= I9 Then 380
340 DISP "Reached maximum iterations limit"
350 GOTO 840
360 !End If
370 !
380 CALL GetDrv1(N, X, G)
390 !
400 ! test if gradient is shallow enough
410 N0 = 0
420 For I = 1 To N
430 N0 = N0 + G(I)^2
440 Next I
450 N0 = Sqr(N0)
460 ! If N0 < E Then GOTO 840
470 !
480 CALL GetDrv2(N, X, J)
490 !
500 MAT J = INV(J)
510 MAT G = J * G
520 !
530 MAT D = G
540 !
550 L8 = 0.1
560 CALL LSrch(X, N, L8, D, K0)
570 If K0 <> 0 Then 600
580 DISP "Linear Search failed"
590 GOTO 840
600 !End If
610 !
620 For I = 1 To N
630 D(I) = L8 * D(I)
640 X(I) = X(I) + D(I)
650 Next I
660 !
670 S0 = 1
680 For I = 1 To N
690 If Abs(D(I)) <= T0(I) Then 720
700 S0 = 0
710 I = N
720 !End If
730 Next I
740 REM
750 CALL MyFx(X, N, F)
760 DISP "F = ";F;" ";
770 For I = 1 To N
780 DISP "X=(";I;")=";X(I);" ";
790 Next I
800 DISP
810 !
820 IF S0 = 1 THEN 310
830 !
840 CALL MyFx(X, N, F)
850 DISP "**********FINAL RESULTS************"
860 DISP "Optimum at:"
870 For I = 1 To N
880 DISP "X(";I;")=";X(I) @ PAUSE
890 Next I
900 For I = 1 To N
910 DISP "Delta X(";I;")=";D(I) @ PAUSE
920 Next I
930 DISP "Function value ="; F @ PAUSE
940 DISP "Number of iterations = ";I0
950 !
960 SUB MyFx(X(), N, R)
970 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
980 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
990 END SUB
1000 !
1010 SUB FrstDrv(N, X(), N0, R)
1020 X0 = X(N0)
1030 h = 0.01 * (1 + Abs(X0))
1040 X(N0) = X0 + h
1050 CALL MyFx(X, N, F1)
1060 X(N0) = X0 - h
1070 CALL MyFx(X, N, F2)
1080 X(N0) = X0
1090 R = (F1 - F2) / 2 / h
1100 END SUB
1110 !
1120 SUB SecDerv(N, X(), N0, J0, R)
1130 ! calculate second derivative?
1140 If N0 <> J0 Then 1250
1150 CALL MyFx(X, N, F0)
1160 X0 = X(N0)
1170 H1 = 0.01 * (1 + Abs(X0))
1180 X(N0) = X0 + H1
1190 CALL MyFx(X, N, F1)
1200 X(N0) = X0 - H1
1210 CALL MyFx(X, N, F2)
1220 X(N0) = X0
1230 R = (F1 - 2 * F0 + F2) / H1 ^ 2
1235 GOTO 1480
1240 !Else
1250 X0 = X(N0)
1260 Y0 = X(J0)
1270 H1 = 0.01 * (1 + Abs(X0))
1280 H2 = 0.01 * (1 + Abs(Y0))
1290 ! calculate F3
1300 X(N0) = X0 + H1
1310 X(J0) = Y0 + H2
1320 CALL MyFx(X, N, F3)
1330 ! calculate F4
1340 X(N0) = X0 - H1
1350 X(J0) = Y0 - H2
1360 CALL MyFx(X, N, F4)
1370 ! calculate F5
1380 X(N0) = X0 + H1
1390 X(J0) = Y0 - H2
1400 CALL MyFx(X, N, F5)
1410 ! calculate F6
1420 X(N0) = X0 - H1
1430 X(J0) = Y0 + H2
1440 CALL MyFx(X, N, F6)
1450 X(N0) = X0
1460 X(J0) = Y0
1470 R = (F3 - F6 - F5 + F4) / (4 * H1 * H2)
1480 !End If
1490 END SUB
1500 !
1510 SUB GetDrv1(N, X(), D1())
1520 For I = 1 To N
1530 CALL FrstDrv(N, X, I, D1(I))
1540 Next I
1550 END SUB
1560 !
1570 SUB GetDrv2(N, X(), D2(,))
1580 For I = 1 To N
1590 For J = 1 To N
1600 CALL SecDerv(N, X, I, J, D2(I, J))
1610 Next J
1620 Next I
1630 END SUB
1640 !
1650 SUB MyFxEx(N,X(),D(), L, R)
1660 DIM X2(N)
1670 For I = 1 To N
1680 X2(I) = X(I) + L * D(I)
1690 Next I
1700 CALL MyFx(X2, N, R)
1710 END SUB
1720 !
1730 SUB LSrch(X(), N, L, D(), R)
1740 M8 = 100
1750 T0 = 0.000001
1760 R = 1
1770 REM
1780 I0 = 0
1790 !REPEAT
1800 I0 = I0 + 1
1810 If I0 <= M8 Then 1850
1820 R = 0
1830 GOTO 1950
1840 !End If
1850 h = 0.01 * (1 + Abs(L))
1860 CALL MyFxEx(N, X, D, L, F0)
1870 CALL MyFxEx(N, X, D, L + h, F1)
1880 CALL MyFxEx(N, X, D, L - h, F2)
1890 D1 = (F1 - F2) / 2 / h
1900 D2 = (F1 - 2 * F0 + F2) / h ^ 2
1910 If D2 = 0 Then 1950
1920 D0 = D1 / D2
1930 L = L - D0
1940 IF Abs(D0) < T0 THEN 1790
1950 END SUB
Copyright (c) Namir Shammas. All rights reserved.