summaryrefslogtreecommitdiffstats
path: root/scilab_doc
diff options
context:
space:
mode:
authorMichaŽl Baudin <michael.baudin@scilab.org>2009-10-23 16:40:59 +0200
committerMichaŽl Baudin <michael.baudin@scilab.org>2009-10-23 16:40:59 +0200
commit8c698ae3174add8a8613bc5789309137eddd0baa (patch)
tree6ca45d679c2849bc45a8efdd7614c366b68c783b /scilab_doc
parentb43168874bbac8775b79d2dd99336f17ac6a523f (diff)
downloadscilab-8c698ae3174add8a8613bc5789309137eddd0baa.zip
scilab-8c698ae3174add8a8613bc5789309137eddd0baa.tar.gz
fminsearch: updated manual/ O'Neill test cases
Diffstat (limited to 'scilab_doc')
-rw-r--r--scilab_doc/neldermead/mckinnon-history-simplex.pngbin10992 -> 6524 bytes
-rw-r--r--scilab_doc/neldermead/method-neldermead.tex648
-rw-r--r--scilab_doc/neldermead/method-spendley.tex3
-rw-r--r--scilab_doc/neldermead/neldermead-dimension-nfevals.pngbin12140 -> 4581 bytes
-rw-r--r--scilab_doc/neldermead/neldermead-dimension-rho.pngbin11013 -> 4328 bytes
-rw-r--r--scilab_doc/neldermead/neldermead-neldermead-so.blg5
-rw-r--r--scilab_doc/neldermead/neldermead-neldermead-so.pdfbin558881 -> 591781 bytes
-rw-r--r--scilab_doc/neldermead/neldermead-neldermead-so.toc21
-rw-r--r--scilab_doc/neldermead/quad2bis-nm-history-logfopt.pngbin12858 -> 3672 bytes
-rw-r--r--scilab_doc/neldermead/quad2bis-nm-history-sigma.pngbin10614 -> 3714 bytes
-rw-r--r--scilab_doc/neldermead/quad2bis-nm-simplexcontours.pngbin32449 -> 7334 bytes
-rw-r--r--scilab_doc/neldermead/scripts/neldermead_ONEILL.sce126
12 files changed, 612 insertions, 191 deletions
diff --git a/scilab_doc/neldermead/mckinnon-history-simplex.png b/scilab_doc/neldermead/mckinnon-history-simplex.png
index a1f81b7..6b363b6 100644
--- a/scilab_doc/neldermead/mckinnon-history-simplex.png
+++ b/scilab_doc/neldermead/mckinnon-history-simplex.png
Binary files differ
diff --git a/scilab_doc/neldermead/method-neldermead.tex b/scilab_doc/neldermead/method-neldermead.tex
index 397f2f7..42186af 100644
--- a/scilab_doc/neldermead/method-neldermead.tex
+++ b/scilab_doc/neldermead/method-neldermead.tex
@@ -239,12 +239,12 @@ simplex in the Nelder-Mead algorithm.
239\label{fig-nm-moves} 239\label{fig-nm-moves}
240\end{figure} 240\end{figure}
241 241
242The figure \ref{fig-nm-moves-reflection} 242The figures \ref{fig-nm-moves-reflection}
243to \ref{fig-nm-moves-shrinkafterco} present the 243to \ref{fig-nm-moves-shrinkafterco} present the
244detailed situations when each type of step occur. 244detailed situations when each type of step occur.
245We emphasize that these figures are not the result of 245We emphasize that these figures are not the result of
246numerical experiments. These figures been created in order 246numerical experiments. These figures been created in order
247to illustrate specific points of the algorithm. 247to illustrate the following specific points of the algorithm.
248 248
249\begin{itemize} 249\begin{itemize}
250\item Obviously, the expansion step is performed when the 250\item Obviously, the expansion step is performed when the
@@ -373,7 +373,7 @@ We measure the rate of convergence defined by
373\begin{eqnarray} 373\begin{eqnarray}
374\label{rho-rate-convergence} 374\label{rho-rate-convergence}
375\rho(S_0,n) = \textrm{lim sup}_{k\rightarrow \infty} 375\rho(S_0,n) = \textrm{lim sup}_{k\rightarrow \infty}
376\left(\sum_{i=0,k-1} \frac{\sigma(S_{i+1})}{\sigma(S_i)}\right)^{1/k} 376\left(\sum_{i=0,k-1} \frac{\sigma(S_{i+1})}{\sigma(S_i)}\right)^{1/k}.
377\end{eqnarray} 377\end{eqnarray}
378That definition can be viewed as the geometric mean of the ratio of the 378That definition can be viewed as the geometric mean of the ratio of the
379oriented lengths between successive simplices and the minimizer 0. 379oriented lengths between successive simplices and the minimizer 0.
@@ -381,7 +381,7 @@ This definition implies
381\begin{eqnarray} 381\begin{eqnarray}
382\label{rho-rate-convergence2} 382\label{rho-rate-convergence2}
383\rho(S_0,n) = \textrm{lim sup}_{k\rightarrow \infty} 383\rho(S_0,n) = \textrm{lim sup}_{k\rightarrow \infty}
384\left( \frac{\sigma(S_{k+1})}{\sigma(S_0)}\right)^{1/k} 384\left( \frac{\sigma(S_{k+1})}{\sigma(S_0)}\right)^{1/k}.
385\end{eqnarray} 385\end{eqnarray}
386 386
387According to the definition, the algorithm is convergent if $\rho(S_0,n) < 1$. 387According to the definition, the algorithm is convergent if $\rho(S_0,n) < 1$.
@@ -888,19 +888,41 @@ presents the multi-directional direct search algorithm.
888 888
889The function we try to minimize is the following quadratic 889The function we try to minimize is the following quadratic
890in 2 dimensions 890in 2 dimensions
891
892\begin{eqnarray} 891\begin{eqnarray}
893f(x_1,x_2) = x_1^2 + x_2^2 - x_1 x_2 892f(x_1,x_2) = x_1^2 + x_2^2 - x_1 x_2.
894\end{eqnarray} 893\end{eqnarray}
895 894
896The stopping criteria is based on the relative size of the simplex 895The stopping criteria is based on the relative size of the simplex
897with respect to the size of the initial simplex 896with respect to the size of the initial simplex
898
899\begin{eqnarray} 897\begin{eqnarray}
900\sigma(S) < tol \times \sigma(S_0) 898\sigma_+(S) < tol \times \sigma_+(S_0),
901\end{eqnarray} 899\end{eqnarray}
900where the tolerance is set to $tol=10^{-8}$.
901
902The initial simplex is a regular simplex with unit length.
903
904The following Scilab script allows to perform the optimization.
905
906\lstset{language=scilabscript}
907\begin{lstlisting}
908function [ y , index ] = quadratic ( x , index )
909 y = x(1)^2 + x(2)^2 - x(1) * x(2);
910endfunction
911nm = neldermead_new ();
912nm = neldermead_configure(nm,"-numberofvariables",2);
913nm = neldermead_configure(nm,"-function",quadratic);
914nm = neldermead_configure(nm,"-x0",[2.0 2.0]');
915nm = neldermead_configure(nm,"-maxiter",100);
916nm = neldermead_configure(nm,"-maxfunevals",300);
917nm = neldermead_configure(nm,"-tolxmethod",%f);
918nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-8);
919nm = neldermead_configure(nm,"-simplex0method","spendley");
920nm = neldermead_configure(nm,"-method","variable");
921nm = neldermead_search(nm);
922neldermead_display(nm);
923nm = neldermead_destroy(nm);
924\end{lstlisting}
902 925
903The initial simplex is computed from the coordinate axis and the unit length.
904The numerical results are presented in table \ref{fig-nm-numexp1-table}. 926The numerical results are presented in table \ref{fig-nm-numexp1-table}.
905 927
906\begin{figure}[htbp] 928\begin{figure}[htbp]
@@ -909,12 +931,12 @@ The numerical results are presented in table \ref{fig-nm-numexp1-table}.
909\begin{tabular}{|l|l|} 931\begin{tabular}{|l|l|}
910\hline 932\hline
911Iterations & 65 \\ 933Iterations & 65 \\
912Function Evaluations & 127 \\ 934Function Evaluations & 130 \\
913$x_0$ & $(2.0,2.0)$ \\ 935$x_0$ & $(2.0,2.0)$ \\
914Relative tolerance on simplex size & $10^{-8}$ \\ 936Relative tolerance on simplex size & $10^{-8}$ \\
915Exact $x^\star$ & $(0.,0.)$\\ 937Exact $x^\star$ & $(0.,0.)$\\
916Computed $x^\star$ & $(7.3e-10 , -2.5e-9)$\\ 938Computed $x^\star$ & $(-2.519D-09 , 7.332D-10)$\\
917Computed $f(x^\star)$ & $8.7e-18$\\ 939Computed $f(x^\star)$ & $8.728930e-018$\\
918\hline 940\hline
919\end{tabular} 941\end{tabular}
920%\end{tiny} 942%\end{tiny}
@@ -924,11 +946,8 @@ $f(x_1,x_2) = x_1^2 + x_2^2 - x_1 x_2$}
924\label{fig-nm-numexp1-table} 946\label{fig-nm-numexp1-table}
925\end{figure} 947\end{figure}
926 948
927
928The various simplices generated during the iterations are 949The various simplices generated during the iterations are
929presented in figure \ref{fig-nm-numexp1-historysimplex}. 950presented in figure \ref{fig-nm-numexp1-historysimplex}.
930The method use reflections in the early iterations. Then there
931is no possible improvement using reflections and shrinking is necessary.
932 951
933\begin{figure} 952\begin{figure}
934\begin{center} 953\begin{center}
@@ -947,11 +966,11 @@ step-by-step evolution of the simplex with the Spendley et al. algorithm.
947\begin{center} 966\begin{center}
948\includegraphics[width=10cm]{quad2bis-nm-history-sigma.png} 967\includegraphics[width=10cm]{quad2bis-nm-history-sigma.png}
949\end{center} 968\end{center}
950\caption{Nelder-Mead numerical experiment -- history of length of simplex} 969\caption{Nelder-Mead numerical experiment -- History of logarithm of length of simplex}
951\label{fig-nm-numexp1-sigma} 970\label{fig-nm-numexp1-sigma}
952\end{figure} 971\end{figure}
953 972
954The convergence is quite fast in this case, since less than 60 iterations 973The convergence is quite fast in this case, since less than 70 iterations
955allow to get a function value lower than $10^{-15}$, as shown in 974allow to get a function value lower than $10^{-15}$, as shown in
956figure \ref{fig-nm-numexp1-logfopt}. 975figure \ref{fig-nm-numexp1-logfopt}.
957 976
@@ -959,7 +978,7 @@ figure \ref{fig-nm-numexp1-logfopt}.
959\begin{center} 978\begin{center}
960\includegraphics[width=10cm]{quad2bis-nm-history-logfopt.png} 979\includegraphics[width=10cm]{quad2bis-nm-history-logfopt.png}
961\end{center} 980\end{center}
962\caption{Nelder-Mead numerical experiment -- history of logarithm of function} 981\caption{Nelder-Mead numerical experiment -- History of logarithm of function}
963\label{fig-nm-numexp1-logfopt} 982\label{fig-nm-numexp1-logfopt}
964\end{figure} 983\end{figure}
965 984
@@ -976,13 +995,42 @@ The more $a$ is large, the more difficult the problem is
976to solve with the simplex algorithm. 995to solve with the simplex algorithm.
977 996
978We set the maximum number of function evaluations to 400. 997We set the maximum number of function evaluations to 400.
979The initial simplex is computed from the coordinate axis and the unit length. 998The initial simplex is a regular simplex with unit length.
999The stopping criteria is based on the relative size of the simplex
1000with respect to the size of the initial simplex
1001\begin{eqnarray}
1002\sigma_+(S) < tol \times \sigma_+(S_0),
1003\end{eqnarray}
1004where the tolerance is set to $tol=10^{-8}$.
1005
1006The following Scilab script allows to perform the optimization.
1007
1008\lstset{language=scilabscript}
1009\begin{lstlisting}
1010a = 100.0;
1011function [ y , index ] = quadratic ( x , index )
1012 y = a * x(1)^2 + x(2)^2;
1013endfunction
1014nm = neldermead_new ();
1015nm = neldermead_configure(nm,"-numberofvariables",2);
1016nm = neldermead_configure(nm,"-function",quadratic);
1017nm = neldermead_configure(nm,"-x0",[10.0 10.0]');
1018nm = neldermead_configure(nm,"-maxiter",400);
1019nm = neldermead_configure(nm,"-maxfunevals",400);
1020nm = neldermead_configure(nm,"-tolxmethod",%f);
1021nm = neldermead_configure(nm,"-tolsimplexizerelative",1.e-8);
1022nm = neldermead_configure(nm,"-simplex0method","spendley");
1023nm = neldermead_configure(nm,"-method","variable");
1024nm = neldermead_search(nm);
1025neldermead_display(nm);
1026nm = neldermead_destroy(nm);
1027\end{lstlisting}
980 1028
981The numerical results are presented in table \ref{fig-nm-numexp2-table}, 1029The numerical results are presented in table \ref{fig-nm-numexp2-table},
982where the experiment is presented for $a=100$. We can check that the 1030where the experiment is presented for $a=100$. We can check that the
983number of function evaluation (161 function evaluations) is much lower than the number 1031number of function evaluation (161 function evaluations) is much lower than the number
984for the fixed shape Spendley et al. method (400 function evaluations) 1032for the fixed shape Spendley et al. method (400 function evaluations)
985and that the function value at optimum is very accurate ($f(x^\star)\approx 1.e-17$ 1033and that the function value at optimum is very accurate ($f(x^\star)\approx 10^{-17}$
986compared to Spendley's et al. $f(x^\star) \approx 0.08$). 1034compared to Spendley's et al. $f(x^\star) \approx 0.08$).
987 1035
988\begin{figure}[h] 1036\begin{figure}[h]
@@ -992,14 +1040,16 @@ compared to Spendley's et al. $f(x^\star) \approx 0.08$).
992\hline 1040\hline
993& Nelder-Mead & Spendley et al.\\ 1041& Nelder-Mead & Spendley et al.\\
994\hline 1042\hline
995Iterations & 83 & 340 \\ 1043Iterations & 82 & 340 \\
996Function Evaluations & 161 & 400 \\ 1044Function Evaluations & 164 & Max=400 \\
997$a$ & $100.0$ & - \\ 1045$a$ & $100.0$ & $100.0$ \\
998$x_0$ & $(10.0,10.0)$ & - \\ 1046$x_0$ & $(10.0,10.0)$ & $(10.0,10.0)$ \\
999Relative tolerance on simplex size & - \\ 1047Initial simplex & regular & regular \\
1000Exact $x^\star$ & $(0.,0.)$ & -\\ 1048Initial simplex length & 1.0 & 1.0 \\
1001Computed $x^\star$ & $(2.e-10, -3.e-9)$& $(0.001,0.2)$\\ 1049Relative tolerance on simplex size & $10^{-8}$ & $10^{-8}$ \\
1002Computed $f(x^\star)$ & $1.e-17$ & $0.08$\\ 1050Exact $x^\star$ & $(0.,0.)$ & $(0.,0.)$ \\
1051Computed $x^\star$ & $(-2.D-10 -1.D-09)$ & $(0.001,0.2)$\\
1052Computed $f(x^\star)$ & $1.D-017$ & $0.08$\\
1003\hline 1053\hline
1004\end{tabular} 1054\end{tabular}
1005%\end{tiny} 1055%\end{tiny}
@@ -1012,7 +1062,7 @@ to the fixed shaped Spendley et al. method.}
1012 1062
1013In figure \ref{fig-nm-numexp2-scaling}, we analyze the 1063In figure \ref{fig-nm-numexp2-scaling}, we analyze the
1014behavior of the method with respect to scaling. 1064behavior of the method with respect to scaling.
1015We check that the method behave very smoothly, with a very 1065We check that the method behaves very smoothly, with a very
1016small number of additional function evaluations when the 1066small number of additional function evaluations when the
1017scaling deteriorates. This shows how much the Nelder-Mead algorithms 1067scaling deteriorates. This shows how much the Nelder-Mead algorithms
1018improves over Spendley's et al. method. 1068improves over Spendley's et al. method.
@@ -1022,17 +1072,19 @@ improves over Spendley's et al. method.
1022%\begin{tiny} 1072%\begin{tiny}
1023\begin{tabular}{|l|l|l|l|} 1073\begin{tabular}{|l|l|l|l|}
1024\hline 1074\hline
1025$a$ & Function evaluations & Computed $f(x^\star)$ & Computed $x^\star$\\ 1075$a$ & Function & Computed $f(x^\star)$ & Computed $x^\star$\\
1026$1.0$ & 139 & $8.0e-18$ & $(2.e-9 -1.e-9)$\\ 1076& Evaluations & & \\
1027$10.0$ & 151 & $7.0e-17$ & $(5.e-10 2.e-9)$\\ 1077\hline
1028$100.0$ & 161 & $1.0e-17$ & $(2.e-10 -3.e-9)$ \\ 1078$1.0$ & 147 & $1.856133e-017$ & $(1.920D-09 , -3.857D-09)$\\
1029$1000.0$ & 165 & $1.0e-17$ & $(-1.e-010 9.e-10)$\\ 1079$10.0$ & 156 & $6.299459e-017$ & $(2.482D-09 , 1.188D-09)$\\
1030$10000.0$ & 167 & $3.0e-17$ & $(5.0e-11,-1.0e-10)$ \\ 1080$100.0$ & 164 & $1.140383e-017$ & $(-2.859D-10 , -1.797D-09)$ \\
1081$1000.0$ & 173 & $2.189830e-018$ & $(-2.356D-12 , 1.478D-09)$\\
1082$10000.0$ & 189 & $1.128684e-017$ & $(2.409D-11 , -2.341D-09)$ \\
1031\hline 1083\hline
1032\end{tabular} 1084\end{tabular}
1033%\end{tiny} 1085%\end{tiny}
1034\end{center} 1086\end{center}
1035\caption{Numerical experiment with Spendley's et al. method on a badly scaled quadratic function} 1087\caption{Numerical experiment with Nelder-Mead method on a badly scaled quadratic function}
1036\label{fig-nm-numexp2-scaling} 1088\label{fig-nm-numexp2-scaling}
1037\end{figure} 1089\end{figure}
1038 1090
@@ -1051,103 +1103,83 @@ in n-dimensions
1051f(\bold{x}) = \sum_{i=1,n} x_i^2. 1103f(\bold{x}) = \sum_{i=1,n} x_i^2.
1052\end{eqnarray} 1104\end{eqnarray}
1053 1105
1054The initial simplex is computed from the coordinate axis and the unit length. 1106The initial simplex is given to the solver.
1055The initial guess is at 0 so that the first vertex is the origin ; 1107The first vertex is the origin ; this vertex is never updated during the iterations.
1056this vertex is never updated during the iterations. 1108The other vertices are based on uniform random numbers in the interval $[-1,1]$.
1109The vertices $i=2,n+1$ are computed from
1110\begin{eqnarray}
1111\bv_i^{(0)} = 2 rand(n,1) - 1,
1112\end{eqnarray}
1113as prescribed by \cite{HanNeumann2006}.
1114In Scilab, the \scifunction{rand} function returns a matrix of
1115uniform random numbers in the interval $[0,1)$.
1057 1116
1058The figure \ref{fig-nm-numexp3-dimension} presents the results of this 1117The stopping criteria is based on the absolute size of the simplex, i.e.
1059experiment for $n=1,19$. 1118the simulation is stopped when
1060 1119\begin{eqnarray}
1061During the iterations, no shrink steps are performed. The 1120\sigma_+(S) < tol,
1062algorithm performs reflections, inside and outside contractions. 1121\end{eqnarray}
1063The figure \ref{fig-nm-numexp3-steps} shows the detailed sequence of 1122where the tolerance is set to $tol=10^{-8}$.
1064iterations for $n=10$. We see that there is no general
1065pattern for the iterations. One can check, however, that there
1066are never no more than $n$ consecutive reflection steps, which is
1067as expected. After one or more contractions, the reflection
1068steps move the worst vertices toward better function values.
1069But there are only $n+1$ vertices so that the $n$ worst
1070vertices are moved in at most $n$ reflection steps.
1071 1123
1072\begin{figure}[htbp] 1124We perform the experiment for $n=1,\ldots,19$.
1073\begin{center} 1125For each experiment, we compute the convergence rate from
1074%\begin{tiny} 1126\begin{eqnarray}
1075\begin{verbatim} 1127\rho(S_0,n) = \left( \frac{\sigma(S_{k})}{\sigma(S_0)}\right)^{1/k},
1076I I I I I I I I I I I I I I I I I I I I R R R R R R R R R R R R R I O 1128\end{eqnarray}
1077R R R R R R R R R R R R R I R I I R I R O I I I I R I R I I R O I R R 1129where $k$ is the number of iterations.
1078R R I R I R I R R R R R R R I R R R R I R I I R I R I I I R R I I I R
1079R R I R R I R R R R R R I R I R R R R R I R R O R R I O I O R R R R I
1080I I O R I R R R R R R I I I R R I I R R R O R I I R R R I R I I O I R
1081I R R O I I R R R R I R R O I R R O R I R I R I R R I R I R R R I I I
1082I I O R R R R I I I R R R I I I R R R I I I I R R R R I I R R R R I R
1083R R I O I R R I I R R R R O I R I I R R R R R R O R R R O I R R I I I
1084I O R I I I R I I I I R R I I R R I R R R R R R I R R I I R R O R I I
1085O R I R O O R O I I R I I I R I I R R R R R R R R R R R R I R R O I R
1086I O I R I I I I R I I R I I R I R O R I O R I R I R R R O R I R R R I
1087I R I R R R I R I R R R R I I R R I R R R R I I R R R R I I R I R I I
1088O I R I I R R R R R R R R I O I R R I I I R I R I I I I R R R R I R R
1089I R I R R R R I I R R R I I I I R I I I I R I R R I R I O R R R I R I
1090O I R R I I R R I R R R R O R R R R I O R R R I R I I I I R I R R R R
1091R I I R I I R R R R R O R R R I R R R R I R I R R I R I I R R I I R R
1092I I I I R R R R R R R R R I R R O R R R R R O I I I I R I I R O I I R
1093R R R R I I I R R
1094\end{verbatim}
1095%\end{tiny}
1096\end{center}
1097\caption{Numerical experiment with Nelder-Mead method on a generalized
1098quadratic function - steps of the algorithm : I = inside contraction, O = outside contraction,
1099R = reflection, S = shrink}
1100\label{fig-nm-numexp3-steps}
1101\end{figure}
1102 1130
1103The figure \ref{fig-nm-numexp3-nbsteps} presents the number and 1131The following Scilab script allows to perform the optimization.
1104the kind of steps performed during the iterations for $n=1,19$.
1105It appears that the number of shrink steps and expansion steps is zero, as expected.
1106More interesting is that the number of reflection is
1107larger than the number of inside contraction when $n$
1108is large. The number of outside contraction is always
1109the smallest in this case.
1110 1132
1111\begin{figure}[htbp] 1133\lstset{language=scilabscript}
1112\begin{center} 1134\begin{lstlisting}
1113%\begin{tiny} 1135function [ f , index ] = quadracticn ( x , index )
1114\begin{tabular}{|l|l|l|l|l|l|} 1136 f = sum(x.^2);
1115\hline 1137endfunction
1116$n$ & \# Reflections & \# Expansion & \# Inside & \# Outside & \#Shrink\\ 1138//
1117 & & & Contractions & Contractions & \\ 1139// solvepb --
1118\hline 1140// Find the solution for the given number of dimensions
11191 & 0 & 0 & 27 & 0 & 0\\ 1141//
11202 & 0 & 0 & 5 & 49 & 0\\ 1142function [nbfevals , niter , rho] = solvepb ( n )
11213 & 54 & 0 & 45 & 36 & 0\\ 1143 rand("seed",0)
11224 & 93 & 0 & 74 & 34 & 0\\ 1144 nm = neldermead_new ();
11235 & 123 & 0 & 101 & 33 & 0\\ 1145 nm = neldermead_configure(nm,"-numberofvariables",n);
11246 & 170 & 0 & 122 & 41 & 0\\ 1146 nm = neldermead_configure(nm,"-function",quadracticn);
11257 & 202 & 0 & 155 & 35 & 0\\ 1147 nm = neldermead_configure(nm,"-x0",zeros(n,1));
11268 & 240 & 0 & 178 & 41 & 0\\ 1148 nm = neldermead_configure(nm,"-maxiter",2000);
11279 & 267 & 0 & 205 & 40 & 0\\ 1149 nm = neldermead_configure(nm,"-maxfunevals",2000);
112810 & 332 & 0 & 234 & 38 & 0\\ 1150 nm = neldermead_configure(nm,"-tolxmethod",%f);
112911 & 381 & 0 & 267 & 36 & 0\\ 1151 nm = neldermead_configure(nm,"-tolsimplexizerelative",0.0);
113012 & 476 & 0 & 299 & 32 & 0\\ 1152 nm = neldermead_configure(nm,"-tolsimplexizeabsolute",1.e-8);
113113 & 473 & 0 & 316 & 42 & 0\\ 1153 nm = neldermead_configure(nm,"-simplex0method","given");
113214 & 545 & 0 & 332 & 55 & 0\\ 1154 coords (1,1:n) = zeros(1,n);
113315 & 577 & 0 & 372 & 41 & 0\\ 1155 for i = 2:n+1
113416 & 635 & 0 & 396 & 46 & 0\\ 1156 coords (i,1:n) = 2.0 * rand(1,n) - 1.0;
113517 & 683 & 0 & 419 & 52 & 0\\ 1157 end
113618 & 756 & 0 & 445 & 55 & 0\\ 1158 nm = neldermead_configure(nm,"-coords0",coords);
113719 & 767 & 0 & 480 & 48 & 0\\ 1159 nm = neldermead_configure(nm,"-method","variable");
1138\hline 1160 nm = neldermead_search(nm);
1139\end{tabular} 1161 si0 = neldermead_get ( nm , "-simplex0" );
1140%\end{tiny} 1162 sigma0 = optimsimplex_size ( si0 , "sigmaplus" );
1141\end{center} 1163 siopt = neldermead_get ( nm , "-simplexopt" );
1142\caption{Numerical experiment with Nelder-Mead method on a generalized quadratic function -- number 1164 sigmaopt = optimsimplex_size ( siopt , "sigmaplus" );
1143and kinds of steps performed} 1165 niter = neldermead_get ( nm , "-iterations" );
1144\label{fig-nm-numexp3-nbsteps} 1166 rho = (sigmaopt/sigma0)^(1.0/niter);
1145\end{figure} 1167 nbfevals = neldermead_get ( nm , "-funevals" );
1168 mprintf ( "%d %d %d %f\n", n , nbfevals , niter , rho );
1169 nm = neldermead_destroy(nm);
1170endfunction
1171// Perform the 20 experiments
1172for n = 1:20
1173 [nbfevals niter rho] = solvepb ( n );
1174 array_rho(n) = rho;
1175 array_nbfevals(n) = nbfevals;
1176 array_niter(n) = niter;
1177end
1178\end{lstlisting}
1146 1179
1147We check that the number of function evaluations 1180The figure \ref{fig-nm-numexp3-dimension} presents the results of this
1148increases approximately linearly with the dimension of the problem in 1181experiment. The rate of convergence, as measured by $\rho(S_0,n)$
1149figure \ref{fig-nm-numexp3-fvn}. A rough rule of thumb is that, for $n=1,19$, 1182converges rapidly toward 1.
1150the number of function evaluations is equal to $100n$.
1151 1183
1152\begin{figure}[htbp] 1184\begin{figure}[htbp]
1153\begin{center} 1185\begin{center}
@@ -1156,25 +1188,26 @@ the number of function evaluations is equal to $100n$.
1156\hline 1188\hline
1157$n$ & Function evaluations & Iterations & $\rho(S_0,n)$\\ 1189$n$ & Function evaluations & Iterations & $\rho(S_0,n)$\\
1158\hline 1190\hline
11591 & 56 & 28 & 0.5125321059829373\\ 11911 & 56 & 27 & 0.513002 \\
11602 & 111 & 55 & 0.71491052830553248\\ 11922 & 113 & 55 & 0.712168 \\
11613 & 220 & 136 & 0.87286283470760984\\ 11933 & 224 & 139 & 0.874043 \\
11624 & 314 & 202 & 0.91247307800713973\\ 11944 & 300 & 187 & 0.904293 \\
11635 & 397 & 258 & 0.93107793607270162\\ 11955 & 388 & 249 & 0.927305 \\
11646 & 503 & 334 & 0.94628781077508028\\ 11966 & 484 & 314 & 0.941782 \\
11657 & 590 & 393 & 0.95404424343636474\\ 11977 & 583 & 383 & 0.951880 \\
11668 & 687 & 460 & 0.96063768057900478\\ 11988 & 657 & 430 & 0.956872 \\
11679 & 767 & 513 & 0.96471820169933631\\ 11999 & 716 & 462 & 0.959721 \\
116810 & 887 & 605 & 0.97000569588245511\\ 120010 & 853 & 565 & 0.966588 \\
116911 & 999 & 685 & 0.97343652480535203\\ 120111 & 910 & 596 & 0.968266 \\
117012 & 1151 & 808 & 0.97745310525741003\\ 120212 & 1033 & 685 & 0.972288 \\
117113 & 1203 & 832 & 0.97803465666405531\\ 120313 & 1025 & 653 & 0.970857 \\
117214 & 1334 & 933 & 0.98042500139065414\\ 120414 & 1216 & 806 & 0.976268 \\
117315 & 1419 & 991 & 0.98154526298964495\\ 120515 & 1303 & 864 & 0.977778 \\
117416 & 1536 & 1078 & 0.98305435726547608\\ 120616 & 1399 & 929 & 0.979316 \\
117517 & 1643 & 1155 & 0.98416149958157839\\ 120717 & 1440 & 943 & 0.979596 \\
117618 & 1775 & 1257 & 0.98544909490809807\\ 120818 & 1730 & 1193 & 0.983774 \\
117719 & 1843 & 1296 & 0.98584701106083183\\ 120919 & 1695 & 1131 & 0.982881 \\
121020 & 1775 & 1185 & 0.983603 \\
1178\hline 1211\hline
1179\end{tabular} 1212\end{tabular}
1180%\end{tiny} 1213%\end{tiny}
@@ -1183,11 +1216,16 @@ $n$ & Function evaluations & Iterations & $\rho(S_0,n)$\\
1183\label{fig-nm-numexp3-dimension} 1216\label{fig-nm-numexp3-dimension}
1184\end{figure} 1217\end{figure}
1185 1218
1219We check that the number of function evaluations
1220increases approximately linearly with the dimension of the problem in
1221figure \ref{fig-nm-numexp3-fvn}. A rough rule of thumb is that, for $n=1,19$,
1222the number of function evaluations is equal to $100n$.
1223
1186\begin{figure} 1224\begin{figure}
1187\begin{center} 1225\begin{center}
1188\includegraphics[width=10cm]{neldermead-dimension-nfevals.png} 1226\includegraphics[width=10cm]{neldermead-dimension-nfevals.png}
1189\end{center} 1227\end{center}
1190\caption{Nelder-Mead numerical experiment -- number of function evaluations 1228\caption{Nelder-Mead numerical experiment -- Number of function evaluations
1191depending on the number of variables} 1229depending on the number of variables}
1192\label{fig-nm-numexp3-fvn} 1230\label{fig-nm-numexp3-fvn}
1193\end{figure} 1231\end{figure}
@@ -1203,7 +1241,7 @@ explained by Han \& Neumann.
1203\begin{center} 1241\begin{center}
1204\includegraphics[width=10cm]{neldermead-dimension-rho.png} 1242\includegraphics[width=10cm]{neldermead-dimension-rho.png}
1205\end{center} 1243\end{center}
1206\caption{Nelder-Mead numerical experiment -- rate of convergence 1244\caption{Nelder-Mead numerical experiment -- Rate of convergence
1207depending on the number of variables} 1245depending on the number of variables}
1208\label{fig-nm-numexp3-rho} 1246\label{fig-nm-numexp3-rho}
1209\end{figure} 1247\end{figure}
@@ -1223,24 +1261,74 @@ particularities
1223\item the expansion is greedy, i.e. the expansion point is accepted if it is better than the lower point, 1261\item the expansion is greedy, i.e. the expansion point is accepted if it is better than the lower point,
1224\item an automatic restart is performed if a factorial test shows that the 1262\item an automatic restart is performed if a factorial test shows that the
1225computed optimum is greater than a local point computed with a relative 1263computed optimum is greater than a local point computed with a relative
1226epsilon equal to 1.e-3. 1264epsilon equal to 1.e-3 and a step equal to the length of the initial simplex.
1227\end{itemize} 1265\end{itemize}
1228 1266
1229The following tests are presented by O'Neill : 1267In order to get an accurate view on O'Neill's factorial test,
1268we must describe explicitely the algorithm.
1269This algorithm is given a vector of lengths, stored in the
1270$step$ variable. It is also given a small value $\epsilon$.
1271The algorithm is presented in figure \ref{algo-factorialtest}.
1272
1273\begin{figure}[htbp]
1274\begin{algorithmic}
1275\STATE $\bx \gets \bx^\star$
1276\STATE $istorestart = FALSE$
1277\FOR{$i = 1$ to $n$}
1278 \STATE $\delta = step ( i ) * \epsilon$
1279 \IF { $\delta == 0.0$ }
1280 \STATE $\delta = \epsilon$
1281 \ENDIF
1282 \STATE $x ( i ) = x ( i ) + \delta$
1283 \STATE $fv = f ( x ) $
1284 \IF { $ fv < fopt $}
1285 \STATE $istorestart = TRUE$
1286 \STATE break
1287 \ENDIF
1288 \STATE $x ( i ) = x ( i ) - \delta - \delta $
1289 \STATE $fv = f ( x ) $
1290 \IF { $ fv < fopt $}
1291 \STATE $istorestart = TRUE$
1292 \STATE break
1293 \ENDIF
1294 \STATE $x ( i ) = x ( i ) + \delta$
1295\ENDFOR
1296\end{algorithmic}
1297\caption{O'Neill's factorial test}
1298\label{algo-factorialtest}
1299\end{figure}
1300
1301O'Neill's factorial test requires a large number of function evaluations, namely
1302$2^n$ function evaluations. In O'Neill's implementation, the parameter
1303$\epsilon$ is set to the constant value $1.e-3$.
1304In Scilab's implementation, this parameter can be customized,
1305thanks to the \scivar{-restarteps} option. Its default value is \scivar{\%eps},
1306the machine epsilon. In O'Neill's implementation, the parameter \scivar{step}
1307is equal to the vector of length used in order to compute the initial
1308simplex. In Scilab's implementation, the two parameters are different,
1309and the \scivar{step} used in the factorial test can be customized
1310with the \scivar{-restartstep} option. Its default value is 1.0, which is
1311expanded into a vector with size $n$.
1312
1230 1313
1314The following tests are presented by O'Neill :
1231\begin{itemize} 1315\begin{itemize}
1232\item Rosenbrock's parabolic valley \cite{citeulike:1903787} 1316\item Rosenbrock's parabolic valley \cite{citeulike:1903787}
1233\begin{eqnarray} 1317\begin{eqnarray}
1234\label{nm-oneill-rosenbrock} 1318\label{nm-oneill-rosenbrock}
1235f(x_1,x_2) = 100(x_2 - x_1^2)^2 + (1-x_1)^2 1319f(x_1,x_2) = 100(x_2 - x_1^2)^2 + (1-x_1)^2
1236\end{eqnarray} 1320\end{eqnarray}
1237with starting point $(x_1,x_2) = (-1.2,1.0)$ 1321with starting point $\bx_0=(x_1,x_2) = (-1.2,1)^T$. The function value at initial guess
1322is $f(\bx_0)=24.2$. The solution is $\bx^\star=(1,1)^T$ where the function
1323value is $f(\bx^\star)=0$.
1238\item Powell's quartic function \cite{Powell08011962} 1324\item Powell's quartic function \cite{Powell08011962}
1239\begin{eqnarray} 1325\begin{eqnarray}
1240\label{nm-oneill-powell} 1326\label{nm-oneill-powell}
1241f(x_1,x_2,x_3,x_4) = (x_1 + 10x_2)^2 + 5 ( x_3 - x_4)^2 + (x_2 - 2x_3)^4 + 10 (x_1 - x_4)^4 1327f(x_1,x_2,x_3,x_4) = (x_1 + 10x_2)^2 + 5 ( x_3 - x_4)^2 + (x_2 - 2x_3)^4 + 10 (x_1 - x_4)^4
1242\end{eqnarray} 1328\end{eqnarray}
1243with starting point $(x_1,x_2,x_3,x_4) = (3,-1,0,1)$ 1329with starting point $\bx_0=(x_1,x_2,x_3,x_4) = (3,-1,0,1)^T$. The function value at initial guess
1330is $f(\bx_0)=215.$. The solution is $\bx^\star=(0,0,0,0)^T$ where the function
1331value is $f(\bx^\star)=0.$.
1244\item Fletcher and Powell's helical valley \cite{R.Fletcher08011963} 1332\item Fletcher and Powell's helical valley \cite{R.Fletcher08011963}
1245\begin{eqnarray} 1333\begin{eqnarray}
1246\label{nm-oneill-fletcherpowell} 1334\label{nm-oneill-fletcherpowell}
@@ -1250,10 +1338,17 @@ f(x_1,x_2,x_3) = 100\left(x_3 + 10\theta(x_1,x_2)\right)^2
1250where 1338where
1251\begin{eqnarray} 1339\begin{eqnarray}
1252\label{nm-oneill-fletcherpowelltheta} 1340\label{nm-oneill-fletcherpowelltheta}
12532\pi \theta(x_1,x_2) &=& \arctan(x_2,x_1), x_1>0\\ 13412\pi \theta(x_1,x_2) &=&
1254&=& \pi + \arctan(x_2,x_1), x_1<0\\ 1342\left\{
1343\begin{array}{ll}
1344\arctan(x_2,x_1), & \textrm{ if } x_1>0\\
1345\pi + \arctan(x_2,x_1), & \textrm{ if } x_1<0
1346\end{array}
1347\right.
1255\end{eqnarray} 1348\end{eqnarray}
1256with starting point $(x_1,x_2,x_3) = (-1,0,0)$. Note 1349with starting point $\bx_0 = (x_1,x_2,x_3) = (-1,0,0)$. The function value at initial guess
1350is $f(\bx_0)=2500$. The solution is $\bx^\star=(1,0,0)^T$ where the function
1351value is $f(\bx^\star)=0.$. Note
1257that since $\arctan(0/0)$ is not defined neither 1352that since $\arctan(0/0)$ is not defined neither
1258the function $f$ on the line $(0,0,x_3)$. This line is excluded 1353the function $f$ on the line $(0,0,x_3)$. This line is excluded
1259by assigning a very large value to the function. 1354by assigning a very large value to the function.
@@ -1262,26 +1357,158 @@ by assigning a very large value to the function.
1262\label{nm-oneill-powers} 1357\label{nm-oneill-powers}
1263f(x_1,\ldots,x_{10}) = \sum_{i=1,10} x_i^4 1358f(x_1,\ldots,x_{10}) = \sum_{i=1,10} x_i^4
1264\end{eqnarray} 1359\end{eqnarray}
1265with starting point $(x_1,\ldots,x_{10}) = (1,\ldots,1)$ 1360with starting point $\bx_0 = (x_1,\ldots,x_{10}) = (1,\ldots,1)$. The function value at initial guess
1361is $f(\bx_0)=10$. The solution is $\bx^\star=(0,\ldots,0)^T$ where the function
1362value is $f(\bx^\star)=0.$.
1266\end{itemize} 1363\end{itemize}
1267 1364
1268The parameters are set to 1365The parameters are set to (following O'Neill's notations)
1269
1270\begin{itemize} 1366\begin{itemize}
1271\item $REQMIN=10^{-16}$, the absolute tolerance on the variance of the function 1367\item $REQMIN=10^{-16}$, the absolute tolerance on the variance of the function
1272values in the simplex, 1368values in the simplex,
1273\item $STEP = 1.0$, the absolute side length of the initial simplex, 1369\item $STEP = 1.0$, the absolute side length of the initial simplex,
1274\item $ICOUNT$, the maximum number of function evaluations. 1370\item $ICOUNT=1000$, the maximum number of function evaluations.
1275\end{itemize} 1371\end{itemize}
1276 1372
1373The following Scilab script allows to define the objective functions.
1374
1375\lstset{language=scilabscript}
1376\begin{lstlisting}
1377// Rosenbrock's "banana" function
1378// initialguess [-1.2 1.0]
1379// xoptimum [1.0 1.0}
1380// foptimum 0.0
1381function [ y , index ] = rosenbrock ( x , index )
1382y = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
1383endfunction
1384// Powell's quartic valley
1385// initialguess [3.0 -1.0 0.0 1.0]
1386// xoptimum [0.0 0.0 0.0 0.0]
1387// foptimum 0.0
1388function [ f , index ] = powellquartic ( x , index )
1389 f = (x(1)+10.0*x(2))^2 + 5.0 * (x(3)-x(4))^2 + (x(2)-2.0*x(3))^4 + 10.0 * (x(1) - x(4))^4
1390endfunction
1391// Fletcher and Powell helical valley
1392// initialguess [-1.0 0.0 0.0]
1393// xoptimum [1.0 0.0 0.0]
1394// foptimum 0.0
1395function [ f , index ] = fletcherpowellhelical ( x , index )
1396 rho = sqrt(x(1) * x(1) + x(2) * x(2))
1397 twopi = 2 * %pi
1398 if ( x(1)==0.0 ) then
1399 f = 1.e154
1400 else
1401 if ( x(1)>0 ) then
1402 theta = atan(x(2)/x(1)) / twopi
1403 elseif ( x(1)<0 ) then
1404 theta = (%pi + atan(x(2)/x(1))) / twopi
1405 end
1406 f = 100.0 * (x(3)-10.0*theta)^2 + (rho - 1.0)^2 + x(3)*x(3)
1407 end
1408endfunction
1409// Sum of powers
1410// initialguess ones(10,1)
1411// xoptimum zeros(10,1)
1412// foptimum 0.0
1413function [ f , index ] = sumpowers ( x , index )
1414 f = sum(x(1:10).^4);
1415endfunction
1416\end{lstlisting}
1417
1418The following Scilab function solves an optimization problem,
1419given the number of parameters, the cost function and the
1420initial guess.
1421
1422\lstset{language=scilabscript}
1423\begin{lstlisting}
1424//
1425// solvepb --
1426// Find the solution for the given problem.
1427// Arguments
1428// n : number of variables
1429// cfun : cost function
1430// x0 : initial guess
1431//
1432function [nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( n , cfun , x0 )
1433 tic();
1434 nm = neldermead_new ();
1435 nm = neldermead_configure(nm,"-numberofvariables",n);
1436 nm = neldermead_configure(nm,"-function",cfun);
1437 nm = neldermead_configure(nm,"-x0",x0);
1438 nm = neldermead_configure(nm,"-maxiter",1000);
1439 nm = neldermead_configure(nm,"-maxfunevals",1000);
1440 nm = neldermead_configure(nm,"-tolxmethod",%f);
1441 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
1442 // Turn ON the tolerance on variance
1443 nm = neldermead_configure(nm,"-tolvarianceflag",%t);
1444 nm = neldermead_configure(nm,"-tolabsolutevariance",1.e-16);
1445 nm = neldermead_configure(nm,"-tolrelativevariance",0.0);
1446 // Turn ON automatic restart
1447 nm = neldermead_configure(nm,"-restartflag",%t);
1448 nm = neldermead_configure(nm,"-restarteps",1.e-3);
1449 nm = neldermead_configure(nm,"-restartstep",1.0);
1450 // Turn ON greedy expansion
1451 nm = neldermead_configure(nm,"-greedy",%t);
1452 // Set initial simplex to axis-by-axis (this is already the default anyway)
1453 nm = neldermead_configure(nm,"-simplex0method","axes");
1454 nm = neldermead_configure(nm,"-simplex0length",1.0);
1455 nm = neldermead_configure(nm,"-method","variable");
1456 //nm = neldermead_configure(nm,"-verbose",1);
1457 //nm = neldermead_configure(nm,"-verbosetermination",1);
1458 //
1459 // Perform optimization
1460 //
1461 nm = neldermead_search(nm);
1462 //neldermead_display(nm);
1463 niter = neldermead_get ( nm , "-iterations" );
1464 nbfevals = neldermead_get ( nm , "-funevals" );
1465 fopt = neldermead_get ( nm , "-fopt" );
1466 xopt = neldermead_get ( nm , "-xopt" );
1467 nbrestart = neldermead_get ( nm , "-restartnb" );
1468 status = neldermead_get ( nm , "-status" );
1469 nm = neldermead_destroy(nm);
1470 cputime = toc();
1471 mprintf ( "=============================\n")
1472 mprintf ( "status = %s\n" , status )
1473 mprintf ( "xopt = [%s]\n" , strcat(string(xopt)," ") )
1474 mprintf ( "fopt = %e\n" , fopt )
1475 mprintf ( "niter = %d\n" , niter )
1476 mprintf ( "nbfevals = %d\n" , nbfevals )
1477 mprintf ( "nbrestart = %d\n" , nbrestart )
1478 mprintf ( "cputime = %f\n" , cputime )
1479 //mprintf ( "%d %d %e %d %f\n", nbfevals , nbrestart , fopt , niter , cputime );
1480endfunction
1481\end{lstlisting}
1482
1483The following Scilab script solves the 4 cases.
1484
1485\lstset{language=scilabscript}
1486\begin{lstlisting}
1487// Solve Rosenbrock's
1488x0 = [-1.2 1.0].';
1489[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 2 , rosenbrock , x0 );
1490
1491// Solve Powell's quartic valley
1492x0 = [3.0 -1.0 0.0 1.0].';
1493[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 4 , powellquartic , x0 );
1494
1495// Solve Fletcher and Powell helical valley
1496x0 = [-1.0 0.0 0.0].';
1497[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 3 , fletcherpowellhelical , x0 );
1498
1499// Solve Sum of powers
1500x0 = ones(10,1);
1501[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 10 , sumpowers , x0 );
1502\end{lstlisting}
1503
1277The table \ref{fig-nm-oneill-table} presents the results which were 1504The table \ref{fig-nm-oneill-table} presents the results which were
1278computed by O'Neill compared with our software. 1505computed by O'Neill compared with Scilab's.
1279For most experiments, the results are very close in terms of 1506For most experiments, the results are very close in terms of
1280number of function evaluations. 1507number of function evaluations.
1281The problem \#4 exhibits a completely different behavior than the 1508The problem \#4 exhibits a different behavior than the
1282results presented by O'Neill. For us, the maximum number of function evaluations 1509results presented by O'Neill. For Scilab, the tolerance on variance
1283is reached (i.e. 1000 function evaluations), whereas for O'Neill, the algorithm 1510of function values is reach after 3 restarts, whereas for O'Neill, the algorithm
1284is restarted and gives the result with 474 function evaluations. 1511is restarted once and gives the result with 474 function evaluations.
1285We did not find any explanation for this behavior. A possible cause of 1512We did not find any explanation for this behavior. A possible cause of
1286difference may be the floating point system which are different and may 1513difference may be the floating point system which are different and may
1287generate different simplices in the algorithms. 1514generate different simplices in the algorithms.
@@ -1291,31 +1518,30 @@ that the numerical experiment were performed by O'Neill
1291on a ICL 4-50 where the two problem 1 and 2 were solved in 3.34 seconds 1518on a ICL 4-50 where the two problem 1 and 2 were solved in 3.34 seconds
1292and the problems 3 and 4 were solved in 22.25 seconds. 1519and the problems 3 and 4 were solved in 22.25 seconds.
1293 1520
1294
1295\begin{figure}[htbp] 1521\begin{figure}[htbp]
1296\begin{center} 1522\begin{center}
1297%\begin{tiny} 1523%\begin{tiny}
1298\begin{tabular}{|l|l|l|l|l|l|l|} 1524\begin{tabular}{|l|l|l|l|l|l|l|}
1299\hline 1525\hline
1300Author & Problem & Function & No. of Restarts & Function Value & Iterations & CPU\\ 1526Author & Problem & Function & Number Of & Function & Iterations & CPU\\
1301& & Evaluations & & & & Time \\ 1527 & & Evaluations & Restarts & Value & & Time \\
1302\hline 1528\hline
1303O'Neill & 1 & 148 & 0 & 3.19 e-9 & ? & ? \\ 1529O'Neill & 1 & 148 & 0 & 3.19e-9 & ? & ? \\
1304Baudin & 1 & 149 & 0 & 1.15 e-7 & 79 & 0.238579 \\ 1530Scilab & 1 & 155 & 0 & 1.158612e-007 & 80 & 0.625000 \\
1305\hline 1531\hline
1306O'Neill & 2 & 209 & 0 & 7.35 e-8 & ? & ? \\ 1532O'Neill & 2 & 209 & 0 & 7.35e-8 & ? & ? \\
1307Baudin & 2 & 224 & 0 & 1.07 e-8 & 126 & 0.447958 \\ 1533Scilab & 2 & 234 & 0 & 1.072588e-008 & 126 & 0.938000 \\
1308\hline 1534\hline
1309O'Neill & 3 & 250 & 0 & 5.29 e-9 & ? & ? \\ 1535O'Neill & 3 & 250 & 0 & 5.29e-9 & ? & ? \\
1310Baudin & 3 & 255 & 0 & 4.56 e-8 & 137 & 0.627493 \\ 1536Scilab & 3 & 263 & 0 & 4.560288e-008 & 137 & 1.037000 \\
1311\hline 1537\hline
1312O'Neill & 4 & 474 & 1 & 3.80 e-7 & ? & ? \\ 1538O'Neill & 4 & 474 & 1 & 3.80e-7 & ? & ? \\
1313Baudin & 4 & 999 & 0 & 5.91 e-9 & 676 & - \\ 1539Scilab & 4 & 616 & 3 & 3.370756e-008 & 402 & 2.949000 \\
1314\hline 1540\hline
1315\end{tabular} 1541\end{tabular}
1316%\end{tiny} 1542%\end{tiny}
1317\end{center} 1543\end{center}
1318\caption{Numerical experiment with Nelder-Mead method on O'Neill test cases - O'Neill results and our results} 1544\caption{Numerical experiment with Nelder-Mead method on O'Neill test cases - O'Neill results and Scilab's results}
1319\label{fig-nm-oneill-table} 1545\label{fig-nm-oneill-table}
1320\end{figure} 1546\end{figure}
1321 1547
@@ -1330,7 +1556,6 @@ method to converge to a non stationnary point.
1330 1556
1331Consider a simplex in two dimensions with vertices at 0 (i.e. the origin), 1557Consider a simplex in two dimensions with vertices at 0 (i.e. the origin),
1332$\bv^{(n+1)}$ and $\bv^{(n)}$. Assume that 1558$\bv^{(n+1)}$ and $\bv^{(n)}$. Assume that
1333
1334\begin{eqnarray} 1559\begin{eqnarray}
1335\label{mckinnon-sortedfv} 1560\label{mckinnon-sortedfv}
1336f(0) < f(\bv^{(n+1)}) < f(\bv^{(n)}). 1561f(0) < f(\bv^{(n+1)}) < f(\bv^{(n)}).
@@ -1339,7 +1564,6 @@ f(0) < f(\bv^{(n+1)}) < f(\bv^{(n)}).
1339The centroid of the simplex is $\overline{\bv} = \bv^{(n+1)}/2$, the midpoint 1564The centroid of the simplex is $\overline{\bv} = \bv^{(n+1)}/2$, the midpoint
1340of the line joining the best and second vertex. The reflected 1565of the line joining the best and second vertex. The reflected
1341point is then computed as 1566point is then computed as
1342
1343\begin{eqnarray} 1567\begin{eqnarray}
1344\label{mckinnon-reflection} 1568\label{mckinnon-reflection}
1345\br^{(n)} = \overline{\bv} + \rho ( \overline{\bv} - \bv^{(n)} ) 1569\br^{(n)} = \overline{\bv} + \rho ( \overline{\bv} - \bv^{(n)} )
@@ -1350,7 +1574,6 @@ Assume that the reflection point $\br^{(n)}$ is rejected, i.e. that
1350$f(\bv^{(n)}) < f(\br^{(n)})$. In this case, the inside contraction 1574$f(\bv^{(n)}) < f(\br^{(n)})$. In this case, the inside contraction
1351step is taken and the point $\bv^{(n+2)}$ is computed using the 1575step is taken and the point $\bv^{(n+2)}$ is computed using the
1352reflection factor $-\gamma = -1/2$ so that 1576reflection factor $-\gamma = -1/2$ so that
1353
1354\begin{eqnarray} 1577\begin{eqnarray}
1355\label{mckinnon-insidecontraction} 1578\label{mckinnon-insidecontraction}
1356\bv^{(n+2)} = \overline{\bv} - 1579\bv^{(n+2)} = \overline{\bv} -
@@ -1361,14 +1584,12 @@ reflection factor $-\gamma = -1/2$ so that
1361Assume then that the inside contraction point is accepted, i.e. $f(\bv^{(n+2)}) < f(\bv^{(n+1)})$. 1584Assume then that the inside contraction point is accepted, i.e. $f(\bv^{(n+2)}) < f(\bv^{(n+1)})$.
1362If this sequence of steps repeats, the simplices are subject to the 1585If this sequence of steps repeats, the simplices are subject to the
1363following linear recurrence formula 1586following linear recurrence formula
1364
1365\begin{eqnarray} 1587\begin{eqnarray}
1366\label{mckinnon-reccurence} 1588\label{mckinnon-reccurence}
13674 \bv^{(n+2)} - \bv^{(n+1)} + 2 \bv^{(n)} = 0 15894 \bv^{(n+2)} - \bv^{(n+1)} + 2 \bv^{(n)} = 0
1368\end{eqnarray} 1590\end{eqnarray}
1369 1591
1370Their general solutions are of the form 1592Their general solutions are of the form
1371
1372\begin{eqnarray} 1593\begin{eqnarray}
1373\bv^{(n)} = \lambda_1^k a_1 + \lambda_2^k a_2 1594\bv^{(n)} = \lambda_1^k a_1 + \lambda_2^k a_2
1374\end{eqnarray} 1595\end{eqnarray}
@@ -1417,9 +1638,58 @@ We consider here the more regular case $\tau=3$, $\theta=6$
1417and $\phi = 400$, i.e. the function is defined by 1638and $\phi = 400$, i.e. the function is defined by
1418\begin{eqnarray} 1639\begin{eqnarray}
1419\label{mckinnon-function3} 1640\label{mckinnon-function3}
1420f(x_1,x_2) &=& - 2400 x_1^3 + x_2 + x_2^2, \qquad x_1\leq 0, \\ 1641f(x_1,x_2) &=&
1421&=& 6 x_1^3 + x_2 + x_2^2, \qquad x_1\geq 0. 1642\left\{
1643\begin{array}{ll}
1644- 2400 x_1^3 + x_2 + x_2^2, & \textrm{ if } x_1\leq 0, \\
16456 x_1^3 + x_2 + x_2^2, & \textrm{ if } x_1\geq 0.
1646\end{array}
1647\right.
1422\end{eqnarray} 1648\end{eqnarray}
1649The solution is $\bx^\star = (0 , -0.5 )^T$.
1650
1651The following Scilab script solves the optimization problem.
1652
1653\lstset{language=scilabscript}
1654\begin{lstlisting}
1655function [ f , index ] = mckinnon3 ( x , index )
1656 if ( length ( x ) ~= 2 )
1657 error ( 'Error: function expects a two dimensional input\n' );
1658 end
1659 tau = 3.0;
1660 theta = 6.0;
1661 phi = 400.0;
1662 if ( x(1) <= 0.0 )
1663 f = theta * phi * abs ( x(1) ).^tau + x(2) * ( 1.0 + x(2) );
1664 else
1665 f = theta * x(1).^tau + x(2) * ( 1.0 + x(2) );
1666 end
1667endfunction
1668lambda1 = (1.0 + sqrt(33.0))/8.0;
1669lambda2 = (1.0 - sqrt(33.0))/8.0;
1670coords0 = [
16711.0 1.0
16720.0 0.0
1673lambda1 lambda2
1674];
1675x0 = [1.0 1.0]';
1676nm = nmplot_new ();
1677nm = nmplot_configure(nm,"-numberofvariables",2);
1678nm = nmplot_configure(nm,"-function",mckinnon3);
1679nm = nmplot_configure(nm,"-x0",x0);
1680nm = nmplot_configure(nm,"-maxiter",200);
1681nm = nmplot_configure(nm,"-maxfunevals",300);
1682nm = nmplot_configure(nm,"-tolfunrelative",10*%eps);
1683nm = nmplot_configure(nm,"-tolxrelative",10*%eps);
1684nm = nmplot_configure(nm,"-simplex0method","given");
1685nm = nmplot_configure(nm,"-coords0",coords0);
1686nm = nmplot_configure(nm,"-simplex0length",1.0);
1687nm = nmplot_configure(nm,"-method","variable");
1688nm = nmplot_search(nm);
1689nmplot_display(nm);
1690nm = nmplot_destroy(nm);
1691\end{lstlisting}
1692
1423 1693
1424The figure \ref{fig-nm-numexp-mckinnon} shows the contour plot of this function and the first 1694The figure \ref{fig-nm-numexp-mckinnon} shows the contour plot of this function and the first
1425steps of the Nelder-Mead method. 1695steps of the Nelder-Mead method.
@@ -1612,19 +1882,19 @@ Author & Step & $f(\bv_1^\star)$ & Function & Angle (\degre)\\
1612& Tolerance & & Evaluations & \\ 1882& Tolerance & & Evaluations & \\
1613\hline 1883\hline
1614Torzcon & 1.e-1 & 7.0355e-5 & 1605 & 89.396677792198 \\ 1884Torzcon & 1.e-1 & 7.0355e-5 & 1605 & 89.396677792198 \\
1615Baudin & 1.e-1 & 8.2272e-5 & 530 & 87.7654 \\ 1885Scilab & 1.e-1 & 8.2272e-5 & 530 & 87.7654 \\
1616\hline 1886\hline
1617Torzcon & 1.e-2 & 6.2912e-5 & 1605 & 89.935373548613 \\ 1887Torzcon & 1.e-2 & 6.2912e-5 & 1605 & 89.935373548613 \\
1618Baudin & 1.e-2 & 7.4854e-5 & 1873 & 89.9253 \\ 1888Scilab & 1.e-2 & 7.4854e-5 & 1873 & 89.9253 \\
1619\hline 1889\hline
1620Torzcon & 1.e-3 & 6.2912e-5 & 3600 & 89.994626919197 \\ 1890Torzcon & 1.e-3 & 6.2912e-5 & 3600 & 89.994626919197 \\
1621Baudin & 1.e-3 & 7.4815e-5 & 2135 & 90.0001 \\ 1891Scilab & 1.e-3 & 7.4815e-5 & 2135 & 90.0001 \\
1622\hline 1892\hline
1623Torzcon & 1.e-4 & 6.2912e-5 & 3670 & 89.999288284747 \\ 1893Torzcon & 1.e-4 & 6.2912e-5 & 3670 & 89.999288284747 \\
1624Baudin & 1.e-4 & 7.481546e-5 & 2196 & 89.9991 \\ 1894Scilab & 1.e-4 & 7.481546e-5 & 2196 & 89.9991 \\
1625\hline 1895\hline
1626Torzcon & 1.e-5 & 6.2912e-5 & 3750 & 89.999931862232 \\ 1896Torzcon & 1.e-5 & 6.2912e-5 & 3750 & 89.999931862232 \\
1627Baudin & 1.e-5 & 7.427212e-5 & 4626 & 89.999990 \\ 1897Scilab & 1.e-5 & 7.427212e-5 & 4626 & 89.999990 \\
1628\hline 1898\hline
1629\end{tabular} 1899\end{tabular}
1630%\end{tiny} 1900%\end{tiny}
diff --git a/scilab_doc/neldermead/method-spendley.tex b/scilab_doc/neldermead/method-spendley.tex
index b6a4724..0c2ae6c 100644
--- a/scilab_doc/neldermead/method-spendley.tex
+++ b/scilab_doc/neldermead/method-spendley.tex
@@ -401,8 +401,7 @@ on simplex size.
401We set the maximum number of function evaluations to 400. 401We set the maximum number of function evaluations to 400.
402The initial simplex is a regular simplex with length unity. 402The initial simplex is a regular simplex with length unity.
403 403
404The following Scilab script uses the neldermead algorithm to perform the 404The following Scilab script allows to perform the optimization.
405optimization.
406 405
407\lstset{language=scilabscript} 406\lstset{language=scilabscript}
408\begin{lstlisting} 407\begin{lstlisting}
diff --git a/scilab_doc/neldermead/neldermead-dimension-nfevals.png b/scilab_doc/neldermead/neldermead-dimension-nfevals.png
index 0e32aae..b5dc68e 100644
--- a/scilab_doc/neldermead/neldermead-dimension-nfevals.png
+++ b/scilab_doc/neldermead/neldermead-dimension-nfevals.png
Binary files differ
diff --git a/scilab_doc/neldermead/neldermead-dimension-rho.png b/scilab_doc/neldermead/neldermead-dimension-rho.png
index 0ab98ff..459ac43 100644
--- a/scilab_doc/neldermead/neldermead-dimension-rho.png
+++ b/scilab_doc/neldermead/neldermead-dimension-rho.png
Binary files differ
diff --git a/scilab_doc/neldermead/neldermead-neldermead-so.blg b/scilab_doc/neldermead/neldermead-neldermead-so.blg
new file mode 100644
index 0000000..57c6cd3
--- /dev/null
+++ b/scilab_doc/neldermead/neldermead-neldermead-so.blg
@@ -0,0 +1,5 @@
1This is BibTeX, Version 0.99cThe top-level auxiliary file: neldermead-neldermead-so.aux
2A level-1 auxiliary file: chapter-notations.aux
3A level-1 auxiliary file: method-neldermead.aux
4The style file: plain.bst
5Database file #1: neldermead.bib
diff --git a/scilab_doc/neldermead/neldermead-neldermead-so.pdf b/scilab_doc/neldermead/neldermead-neldermead-so.pdf
index 55c1403..b41a233 100644
--- a/scilab_doc/neldermead/neldermead-neldermead-so.pdf
+++ b/scilab_doc/neldermead/neldermead-neldermead-so.pdf
Binary files differ
diff --git a/scilab_doc/neldermead/neldermead-neldermead-so.toc b/scilab_doc/neldermead/neldermead-neldermead-so.toc
new file mode 100644
index 0000000..493f697
--- /dev/null
+++ b/scilab_doc/neldermead/neldermead-neldermead-so.toc
@@ -0,0 +1,21 @@
1\contentsline {chapter}{\numberline {1}Nelder-Mead method}{4}{chapter.1}
2\contentsline {section}{\numberline {1.1}Introduction}{4}{section.1.1}
3\contentsline {subsection}{\numberline {1.1.1}Overview}{4}{subsection.1.1.1}
4\contentsline {subsection}{\numberline {1.1.2}Algorithm}{5}{subsection.1.1.2}
5\contentsline {section}{\numberline {1.2}Geometric analysis}{9}{section.1.2}
6\contentsline {section}{\numberline {1.3}Convergence properties on a quadratic}{9}{section.1.3}
7\contentsline {subsection}{\numberline {1.3.1}With default parameters}{13}{subsection.1.3.1}
8\contentsline {subsection}{\numberline {1.3.2}With variable parameters}{17}{subsection.1.3.2}
9\contentsline {section}{\numberline {1.4}Numerical experiments}{20}{section.1.4}
10\contentsline {subsection}{\numberline {1.4.1}Quadratic function}{21}{subsection.1.4.1}
11\contentsline {subsubsection}{Badly scaled quadratic function}{22}{section*.3}
12\contentsline {subsection}{\numberline {1.4.2}Sensitivity to dimension}{24}{subsection.1.4.2}
13\contentsline {subsection}{\numberline {1.4.3}O'Neill test cases}{27}{subsection.1.4.3}
14\contentsline {subsection}{\numberline {1.4.4}Convergence to a non stationnary point}{32}{subsection.1.4.4}
15\contentsline {subsection}{\numberline {1.4.5}Han counter examples}{36}{subsection.1.4.5}
16\contentsline {subsubsection}{First counter example}{36}{section*.4}
17\contentsline {subsubsection}{Second counter example}{36}{section*.5}
18\contentsline {subsection}{\numberline {1.4.6}Torczon's numerical experiments}{37}{subsection.1.4.6}
19\contentsline {subsubsection}{Penalty \#1}{38}{section*.6}
20\contentsline {section}{\numberline {1.5}Conclusion}{40}{section.1.5}
21\contentsline {chapter}{Bibliography}{41}{section.1.5}
diff --git a/scilab_doc/neldermead/quad2bis-nm-history-logfopt.png b/scilab_doc/neldermead/quad2bis-nm-history-logfopt.png
index 83f1bd5..c5d310d 100644
--- a/scilab_doc/neldermead/quad2bis-nm-history-logfopt.png
+++ b/scilab_doc/neldermead/quad2bis-nm-history-logfopt.png
Binary files differ
diff --git a/scilab_doc/neldermead/quad2bis-nm-history-sigma.png b/scilab_doc/neldermead/quad2bis-nm-history-sigma.png
index 92207c4..d1ae6bb 100644
--- a/scilab_doc/neldermead/quad2bis-nm-history-sigma.png
+++ b/scilab_doc/neldermead/quad2bis-nm-history-sigma.png
Binary files differ
diff --git a/scilab_doc/neldermead/quad2bis-nm-simplexcontours.png b/scilab_doc/neldermead/quad2bis-nm-simplexcontours.png
index 8a33ec8..232fa86 100644
--- a/scilab_doc/neldermead/quad2bis-nm-simplexcontours.png
+++ b/scilab_doc/neldermead/quad2bis-nm-simplexcontours.png
Binary files differ
diff --git a/scilab_doc/neldermead/scripts/neldermead_ONEILL.sce b/scilab_doc/neldermead/scripts/neldermead_ONEILL.sce
new file mode 100644
index 0000000..a9be6eb
--- /dev/null
+++ b/scilab_doc/neldermead/scripts/neldermead_ONEILL.sce
@@ -0,0 +1,126 @@
1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2// Copyright (C) 2008-2009 - INRIA - Michael Baudin
3//
4// This file must be used under the terms of the CeCILL.
5// This source file is licensed as described in the file COPYING, which
6// you should have received as part of this distribution. The terms
7// are also available at
8// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10//
11// Performs O'Neill test cases
12//
13// Rosenbrock's "banana" function
14// initialguess [-1.2 1.0]
15// xoptimum [1.0 1.0}
16// foptimum 0.0
17function [ y , index ] = rosenbrock ( x , index )
18y = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
19endfunction
20// Powell's quartic valley
21// initialguess [3.0 -1.0 0.0 1.0]
22// xoptimum [0.0 0.0 0.0 0.0]
23// foptimum 0.0
24function [ f , index ] = powellquartic ( x , index )
25 f = (x(1)+10.0*x(2))^2 + 5.0 * (x(3)-x(4))^2 + (x(2)-2.0*x(3))^4 + 10.0 * (x(1) - x(4))^4
26endfunction
27// Fletcher and Powell helical valley
28// initialguess [-1.0 0.0 0.0]
29// xoptimum [1.0 0.0 0.0]
30// foptimum 0.0
31function [ f , index ] = fletcherpowellhelical ( x , index )
32 rho = sqrt(x(1) * x(1) + x(2) * x(2))
33 twopi = 2 * %pi
34 if ( x(1)==0.0 ) then
35 f = 1.e154
36 else
37 if ( x(1)>0 ) then
38 theta = atan(x(2)/x(1)) / twopi
39 elseif ( x(1)<0 ) then
40 theta = (%pi + atan(x(2)/x(1))) / twopi
41 end
42 f = 100.0 * (x(3)-10.0*theta)^2 + (rho - 1.0)^2 + x(3)*x(3)
43 end
44endfunction
45// Sum of powers
46// initialguess ones(10,1)
47// xoptimum zeros(10,1)
48// foptimum 0.0
49function [ f , index ] = sumpowers ( x , index )
50 f = sum(x(1:10).^4);
51endfunction
52
53//
54// solvepb --
55// Find the solution for the given problem.
56// Arguments
57// n : number of variables
58// cfun : cost function
59// x0 : initial guess
60//
61function [nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( n , cfun , x0 )
62 tic();
63 nm = neldermead_new ();
64 nm = neldermead_configure(nm,"-numberofvariables",n);
65 nm = neldermead_configure(nm,"-function",cfun);
66 nm = neldermead_configure(nm,"-x0",x0);
67 nm = neldermead_configure(nm,"-maxiter",1000);
68 nm = neldermead_configure(nm,"-maxfunevals",1000);
69 nm = neldermead_configure(nm,"-tolxmethod",%f);
70 nm = neldermead_configure(nm,"-tolsimplexizemethod",%f);
71 // Turn ON the tolerance on variance
72 nm = neldermead_configure(nm,"-tolvarianceflag",%t);
73 nm = neldermead_configure(nm,"-tolabsolutevariance",1.e-16);
74 nm = neldermead_configure(nm,"-tolrelativevariance",0.0);
75 // Turn ON automatic restart
76 nm = neldermead_configure(nm,"-restartflag",%t);
77 nm = neldermead_configure(nm,"-restarteps",1.e-3);
78 nm = neldermead_configure(nm,"-restartstep",1.0);
79 // Turn ON greedy expansion
80 nm = neldermead_configure(nm,"-greedy",%t);
81 // Set initial simplex to axis-by-axis (this is already the default anyway)
82 nm = neldermead_configure(nm,"-simplex0method","axes");
83 nm = neldermead_configure(nm,"-simplex0length",1.0);
84 nm = neldermead_configure(nm,"-method","variable");
85 //nm = neldermead_configure(nm,"-verbose",1);
86 //nm = neldermead_configure(nm,"-verbosetermination",1);
87 //
88 // Perform optimization
89 //
90 nm = neldermead_search(nm);
91 //neldermead_display(nm);
92 niter = neldermead_get ( nm , "-iterations" );
93 nbfevals = neldermead_get ( nm , "-funevals" );
94 fopt = neldermead_get ( nm , "-fopt" );
95 xopt = neldermead_get ( nm , "-xopt" );
96 nbrestart = neldermead_get ( nm , "-restartnb" );
97 status = neldermead_get ( nm , "-status" );
98 nm = neldermead_destroy(nm);
99 cputime = toc();
100 mprintf ( "=============================\n")
101 mprintf ( "status = %s\n" , status )
102 mprintf ( "xopt = [%s]\n" , strcat(string(xopt)," ") )
103 mprintf ( "fopt = %e\n" , fopt )
104 mprintf ( "niter = %d\n" , niter )
105 mprintf ( "nbfevals = %d\n" , nbfevals )
106 mprintf ( "nbrestart = %d\n" , nbrestart )
107 mprintf ( "cputime = %f\n" , cputime )
108 //mprintf ( "%d %d %e %d %f\n", nbfevals , nbrestart , fopt , niter , cputime );
109endfunction
110
111// Solve Rosenbrock's
112x0 = [-1.2 1.0].';
113[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 2 , rosenbrock , x0 );
114
115// Solve Powell's quartic valley
116x0 = [3.0 -1.0 0.0 1.0].';
117[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 4 , powellquartic , x0 );
118
119// Solve Fletcher and Powell helical valley
120x0 = [-1.0 0.0 0.0].';
121[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 3 , fletcherpowellhelical , x0 );
122
123// Solve Sum of powers
124x0 = ones(10,1);
125[nbfevals , niter , nbrestart , fopt , cputime ] = solvepb ( 10 , sumpowers , x0 );
126