Nelder Mead Optimization

[rawi]

Edit: Improved version
Improvements:
Maximum number of iterations is set
Better output

2nd Edit: Error correction in FUNCX
3rd edit: Typo in program corrected, program shortened

For me one of the best things of the HP71B was the curve fit module with the possbility to optimize functions with up to 5 dimensions.

This program uses Nelder Mead method to optimize functions with the number of dimensions only limited by memory. Nelder Mead optimization is a type of simplex routine that uses only values of the function. It is very robust and no derivates are needed.

Preparation:
Define function FUNCX() that is to optimize. FUNCX() has to take values of parameters from list L1 and return the function value.
The program determines parameters so that function is minimized. If you want to maximize use the negative.
Put in starting values for the optimization process in List L1. See that the list has the same dimension as the dimension of the problem.
Program asks for the initial step-width S0 for the computation of the starting simplex (e.g. 0.5), for the termination criterion E0 (small number, e.g. 1E-5) and for the maximum number of iterations (e.g. 200). Type them and press OK. Screen shows number of iterations, minimum reached so far and mean difference between simplex points and stops after the mean difference between points is less than termination criterion E0.

Result:
fmin: Function value at optimum in Y


Parameters:
D0 = Dimension (min: 2, is taken from size of L1)
S0 = Starting step width for simplex (e.g. 0.5)
Y = Value of function
L1 = List of x values (at the beginning starting values, at the end result)
L2 = Reflection point coordinates
L3 = List of function values of simplex
E0 = Termination criterion (if average side length of simplex is smaller, process stops)
M1 = Matrix of simplex

Example:
Minimize the following function for the parameters x1-x5::
f(x1,x2,x3,x4,x5) = (x1-2)²+(x2-1)²+(x3-1.5)²+x4²+(x5+1)²+2*((x1-2)*(x2-1))^2+0.5*((x2-1)*(x5+1))^2
Exact solution is x1 = 2, X2 = 1, X3 = 1.5, X4=0 and X5 = -1 with f =0.

Type in the following program:
Function FUNCX():
EXPORT FUNCX()
BEGIN
LOCAL Q1;
(L1(1)-2.)^2+(L1(2)-1)^2+(L1(3)-1.5)^2▶Q1;
Q1+L1(4)^2+(L1(5)+1)^2▶Q1;
Q1+2*((L1(1)-2)*(L1(2)-1))^2▶Q1;
Q1+0.5*((L1(2)-1)*(L1(5)+1))^2▶Q1;
RETURN Q1;
END;

We take {0,0,0,0,0} as starting value. Put this in list L1. If there are more than 5 cells in L1 delete them. Execute NMO. You are asked for S0 and E0 and maximum number of iterations ITM. We take S0=0.5 and E0=1.E-5 and ITM=200. Press OK. Program shows minimum reached during computation and stops after a few seconds with the message “DONE” and number of Iterations (159), Termination criterion (6.23E-6), Minimum Function Fmin (1.44E-11) and difference between maximum and minimum in simplex (1.2E-10). So it has converged.
Result: List 1 shows: 2.00000, 1.00000, 1.50000, -1.8984E-6, -1.00000

Code:

EXPORT NMO(S0,E0,ITM)
BEGIN
// Function to minimize FUNCX()
// S0 Starting value of simplex size
// E0 Termination criterion of simplex
// ITM Maximum number of iterations
// Parameters of function: L1
LOCAL I1,I2,I3,I4,IP,J1,J2,D0;
LOCAL fmin,fmax,fdiff;
LOCAL imin,imax,flag,ZAK;
Size(L1)▶D0;
1▶ZAK;
0▶J1;
0▶I4;
MAKEMAT(0,D0,D0+1)▶M1;
MAKELIST(0.,J1,1,D0)▶L2;
MAKELIST(0.,J1,1,D0+1)▶L3;
For I1 FROM 1 TO D0 DO //Simplex
 J1+I1▶J1;
 1▶J2;
 FOR I2 FROM I1+1 TO D0+1 DO
  √(J1)/I1/J2*S0▶M1(I1,I2);
  I1+1▶J2;
 END;
 FOR I2 FROM 1 TO D0+1 DO
  M1(I1,I2)+L1(I1)▶M1(I1,I2);
 END;
END;
0▶flag;
REPEAT // Main loop
 If flag=0 THEN 
  FOR I1 FROM 1 TO D0+1 DO
   FOR I2 FROM 1 TO D0 DO
    M1(I2,I1)▶L1(I2);
   END;
   FUNCX()▶L3(I1);
  END;
 END;
 0▶flag;
 1ᴇ99▶fmin;
 −1ᴇ99▶fmax; 
 FOR I1 FROM 1 TO D0+1 DO
  IF L3(I1)<fmin Then
   I1▶imin;
   L3(I1)▶fmin;
  END;
  IF L3(I1)>fmax Then
   I1▶imax;
   L3(I1)▶fmax;
  END;
 END;
 FOR I1 FROM 1 TO D0 DO // Reflection center
  0▶L2(I1);
  FOR I2 FROM 1 TO D0+1 DO
   L2(I1)+M1(I1,I2)▶L2(I1);
  END;
  (L2(I1)-M1(I1,imax))/D0▶L2(I1);
  2.*L2(I1)-M1(I1,imax)▶L1(I1); // New point
 END;
 FUNCX()▶Y;
 IF Y < fmin THEN // New point is min 
 FOR I1 FROM 1 TO D0 DO
   3.*L2(I1)-2.*M1(I1,imax)▶L1(I1);
  END;
  Y▶fmin;
  FUNCX()▶Y;
  IF Y > fmin THEN
   FOR I1 FROM 1 TO D0 DO
    L1(I1)-L2(I1)+M1(I1,imax)▶L1(I1);
   END;
  ELSE
   Y▶fmin;
  END;
 END;
 IF Y > fmax THEN // New point is max
  FOR I1 FROM 1 TO D0 DO
   1.5*L3(I1)-0.5*M1(I1,imax)▶L1(I1);
  END;
  FUNCX()▶Y;
  If Y > fmax Then
   FOR I1 FROM 1 TO D0 DO
    FOR I2 FROM 1 TO D0+1 DO
     (M1(I1,I2)+M1(I1,imin))*0.5▶M1(I1,I2);
    END;
   END;
   1▶flag; 
  END; 
 END;
 IF Y < fmax THEN // Updating simplex
  FOR I1 FROM 1 TO D0 DO
   L1(I1)▶M1(I1,imax);
  END;
  Y▶L3(imax);
 END;
 0▶ZAK; // Termination criterion
 FOR I1 FROM 1 TO D0 DO
  FOR I2 FROM 1 TO D0+1 DO
   ZAK+(M1(I1,I2)-M1(I1,imin))^2▶ZAK;
  END;
 END;
 √(ZAK/D0)▶ZAK;
 I4+1▶I4;
 PRINT("--- Iter,fmin,fmax,EPS");
 PRINT(I4);
 PRINT(fmin);
 PRINT(fmax);
 PRINT(ZAK);
 (I4>=ITM) + (ZAK<E0)▶IP; 
UNTIL IP>0;
FOR I1 FROM 1 TO D0 DO // Final work
 M1(I1,imin)▶L1(I1);
END;
fmin▶Y;
fmax-fmin▶fdiff;
Print("--- DONE -----------");
Print("Iter: "+I4);
Print("Eps:  "+ZAK);
Print("Fmin: "+fmin);
Print("Fdiff:"+fdiff);
Print("Paramater: L1");
END;