summaryrefslogtreecommitdiffstats
path: root/scilab_doc
diff options
context:
space:
mode:
authorMichaŽl Baudin <michael.baudin@scilab.org>2009-10-06 16:57:41 +0200
committerMichaŽl Baudin <michael.baudin@scilab.org>2009-10-06 16:57:41 +0200
commitfa95791eb8e4c6c04e8213dd68b41c15cb349d99 (patch)
treedb3df7085355154333da6aa4241723ad07725c62 /scilab_doc
parente81014e45c050cddfd7bc0bc7748d8da45513b48 (diff)
downloadscilab-fa95791eb8e4c6c04e8213dd68b41c15cb349d99.zip
scilab-fa95791eb8e4c6c04e8213dd68b41c15cb349d99.tar.gz
Created demos for -other- fminsearch
Diffstat (limited to 'scilab_doc')
-rw-r--r--scilab_doc/neldermead/method-spendley.tex111
-rw-r--r--scilab_doc/neldermead/neldermead-spendley-so.pdfbin301166 -> 301107 bytes
-rw-r--r--scilab_doc/neldermead/scripts/fminsearch_demo.m78
-rw-r--r--scilab_doc/neldermead/scripts/onehump.m (renamed from scilab_doc/neldermead/scripts/fminsearch_banana.m)9
-rw-r--r--scilab_doc/neldermead/scripts/outfun.m16
-rw-r--r--scilab_doc/neldermead/testFminsearchPlotMatlab.pngbin0 -> 2928 bytes
-rw-r--r--scilab_doc/neldermead/testFminsearchPlotScilab.pngbin0 -> 3183 bytes
7 files changed, 169 insertions, 45 deletions
diff --git a/scilab_doc/neldermead/method-spendley.tex b/scilab_doc/neldermead/method-spendley.tex
index 62a9df1..adb1990 100644
--- a/scilab_doc/neldermead/method-spendley.tex
+++ b/scilab_doc/neldermead/method-spendley.tex
@@ -194,7 +194,7 @@ optimum is searched.
194 194
195\begin{figure} 195\begin{figure}
196\begin{center} 196\begin{center}
197\includegraphics[width=10cm]{spendley-steps-shrink.pdf} 197\includegraphics[width=8cm]{spendley-steps-shrink.pdf}
198\end{center} 198\end{center}
199\caption{Spendley et al. simplex moves -- The shrink step is the 199\caption{Spendley et al. simplex moves -- The shrink step is the
200only possible move.} 200only possible move.}
@@ -235,9 +235,6 @@ In all cases, the difficulty is that a sequence of simplices produced
235by the Nelder-Mead simplex method can come arbitrarily close to 235by the Nelder-Mead simplex method can come arbitrarily close to
236degeneracy. 236degeneracy.
237 237
238%% TODO...
239
240
241\section{Numerical experiments} 238\section{Numerical experiments}
242 239
243In this section, we present some numerical experiments 240In this section, we present some numerical experiments
@@ -269,6 +266,8 @@ where $\|.\|_2$ is the euclidian norm defined by
269\begin{eqnarray} 266\begin{eqnarray}
270\|\bx\|_2 = \sum_{i=1,n}x_i^2. 267\|\bx\|_2 = \sum_{i=1,n}x_i^2.
271\end{eqnarray} 268\end{eqnarray}
269In this experiment, we use $tol=10^{-8}$ as the relative tolerance
270on simplex size.
272 271
273The initial simplex is a regular simplex with length unity. 272The initial simplex is a regular simplex with length unity.
274 273
@@ -306,11 +305,12 @@ The numerical results are presented in table \ref{fig-spendley-numexp1-table}.
306\hline 305\hline
307Iterations & 49 \\ 306Iterations & 49 \\
308Function Evaluations & 132 \\ 307Function Evaluations & 132 \\
309$x_0$ & $(2.0,2.0)$ \\ 308$\bx_0$ & $(2.0,2.0)$ \\
310Relative tolerance on simplex size & $10^{-8}$ \\ 309Relative tolerance on simplex size & $10^{-8}$ \\
311Exact $x^\star$ & $(0.,0.)$\\ 310Exact $\bx^\star$ & $(0.,0.)$\\
312Computed $x^\star$ & $(2.169e-10, 2.169e-10)$\\ 311Computed $\bx^\star$ & $(2.169e-10, 2.169e-10)$\\
313Computed $f(x^\star)$ & $4.706e-20$\\ 312Exact $f(\b\bx^\star)$ & $0.$\\
313Computed $f(\bx^\star)$ & $4.706e-20$\\
314\hline 314\hline
315\end{tabular} 315\end{tabular}
316%\end{tiny} 316%\end{tiny}
@@ -370,6 +370,30 @@ f(x_1,x_2) = a x_1^2 + x_2^2,
370where $a>0$ is a chosen scaling parameter. 370where $a>0$ is a chosen scaling parameter.
371The more $a$ is large, the more difficult the problem is 371The more $a$ is large, the more difficult the problem is
372to solve with the simplex algorithm. 372to solve with the simplex algorithm.
373Indeed, let us compute the Hessian matrix associated with the
374cost function. We have
375\begin{eqnarray}
376\bH = \left(
377\begin{array}{cc}
3782 a & 0 \\
3790 & 2
380\end{array}
381\right).
382\end{eqnarray}
383Therefore, the eigenvalues of the Hessian matrix
384are $2a$ and $2$, which are stricly
385positive if $a>0$. Hence, the cost function is stricly convex and
386there is only one global solution, that is $\bx^\star = (0,0,\ldots,0)^T$.
387The ratio between these two eigenvalues is $a$. This leads
388to an elongated valley, which is extremely narrow when $a$ is large.
389
390The stopping criteria is based on the relative size of the simplex
391with respect to the size of the initial simplex
392\begin{eqnarray}
393\sigma_+(S) < tol \times \sigma_+(S_0).
394\end{eqnarray}
395In this experiment, we use $tol=10^{-8}$ as the relative tolerance
396on simplex size.
373 397
374We set the maximum number of function evaluations to 400. 398We set the maximum number of function evaluations to 400.
375The initial simplex is a regular simplex with length unity. 399The initial simplex is a regular simplex with length unity.
@@ -379,8 +403,9 @@ optimization.
379 403
380\lstset{language=scilabscript} 404\lstset{language=scilabscript}
381\begin{lstlisting} 405\begin{lstlisting}
406a = 100;
382function y = quadratic (x) 407function y = quadratic (x)
383 y = 100 * x(1)^2 + x(2)^2; 408 y = a * x(1)^2 + x(2)^2;
384endfunction 409endfunction
385nm = nmplot_new (); 410nm = nmplot_new ();
386nm = nmplot_configure(nm,"-numberofvariables",2); 411nm = nmplot_configure(nm,"-numberofvariables",2);
@@ -405,9 +430,9 @@ nm = nmplot_destroy(nm);
405 430
406 431
407The numerical results are presented in table \ref{fig-spendley-numexp1-table}, 432The numerical results are presented in table \ref{fig-spendley-numexp1-table},
408where the experiment is presented for $a=100$. One can check that the 433where the experiment is presented for $a=100$. We can check that the
409number of function evaluation is equal to its maximum limit, even if the value of the 434number of function evaluations is equal to its maximum limit, even if the value of the
410function at optimum is very inaccurate ($f(x^\star) \approx 0.08$). 435function at optimum is very inaccurate ($f(\bx^\star) \approx 0.08$).
411 436
412\begin{figure}[h] 437\begin{figure}[h]
413\begin{center} 438\begin{center}
@@ -417,11 +442,11 @@ function at optimum is very inaccurate ($f(x^\star) \approx 0.08$).
417Iterations & 340 \\ 442Iterations & 340 \\
418Function Evaluations & 400 \\ 443Function Evaluations & 400 \\
419$a$ & $100.0$ \\ 444$a$ & $100.0$ \\
420$x_0$ & $(10.0,10.0)$ \\ 445$\bx_0$ & $(10.0,10.0)$ \\
421Relative tolerance on simplex size & $10^{-8}$ \\ 446Relative tolerance on simplex size & $10^{-8}$ \\
422Exact $x^\star$ & $(0.,0.)$\\ 447Exact $\bx^\star$ & $(0.,0.)$\\
423Computed $x^\star$ & $(0.001,0.2)$\\ 448Computed $\bx^\star$ & $(0.001,0.2)$\\
424Computed $f(x^\star)$ & $0.08$\\ 449Computed $f(\bx^\star)$ & $0.08$\\
425\hline 450\hline
426\end{tabular} 451\end{tabular}
427%\end{tiny} 452%\end{tiny}
@@ -430,14 +455,12 @@ Computed $f(x^\star)$ & $0.08$\\
430\label{fig-spendley-numexp2-table} 455\label{fig-spendley-numexp2-table}
431\end{figure} 456\end{figure}
432 457
433
434The various simplices generated during the iterations are 458The various simplices generated during the iterations are
435presented in figure \ref{fig-spendley-numexp2-historysimplex}. 459presented in figure \ref{fig-spendley-numexp2-historysimplex}.
436The method use reflections in the early iterations. Then there 460The method use reflections in the early iterations. Then there
437is no possible improvement using reflections and shrinking is necessary. 461is no possible improvement using reflections, so that shrinking is necessary.
438But the shrinking makes the simplex very small so that a large number of 462But the repeated shrink steps makes the simplex very small, leading to a large number of
439iterations are necessary to improve the function value. 463iterations. This is a limitation of the method, which is based on a simplex
440This is a limitation of the method, which is based on a simplex
441which can vary its size, but not its shape. 464which can vary its size, but not its shape.
442 465
443\begin{figure} 466\begin{figure}
@@ -459,7 +482,7 @@ when $a>10$
459%\begin{tiny} 482%\begin{tiny}
460\begin{tabular}{|l|l|l|} 483\begin{tabular}{|l|l|l|}
461\hline 484\hline
462$a$ & Function evaluations & Computed $f(x^\star)$ \\ 485$a$ & Function evaluations & Computed $f(\bx^\star)$ \\
463$1.0$ & 160 & $2.35e-18$ \\ 486$1.0$ & 160 & $2.35e-18$ \\
464$10.0$ & 222 & $1.2e-17$ \\ 487$10.0$ & 222 & $1.2e-17$ \\
465$100.0$ & 400 & $0.083$ \\ 488$100.0$ & 400 & $0.083$ \\
@@ -476,7 +499,11 @@ $10000.0$ & 400 & $56.08$ \\
476\subsection{Sensitivity to dimension} 499\subsection{Sensitivity to dimension}
477 500
478In this section, we try to study the convergence of the 501In this section, we try to study the convergence of the
479Spendley et al. algorithm with respect to the number of variables. 502Spendley et al. algorithm with respect to the number of variables,
503as presented by Han \& Neumann in \cite{HanNeumann2006}.
504We emphasize, though, that Han \& Neumann present their numerical
505experiment with the Nelder-Mead algorithm, while we present
506in this section the Spendley et al. algorithm.
480The function we try to minimize is the following quadratic function 507The function we try to minimize is the following quadratic function
481in $n$-dimensions 508in $n$-dimensions
482\begin{eqnarray} 509\begin{eqnarray}
@@ -484,14 +511,16 @@ in $n$-dimensions
484f(\bx) = \sum_{i=1,n} x_i^2. 511f(\bx) = \sum_{i=1,n} x_i^2.
485\end{eqnarray} 512\end{eqnarray}
486 513
487The initial guess is at 0 so that this vertex is never updated 514The initial guess is the origin $\bx_0 = (0,0,\ldots,0)^T$,
515which is also the global solution of the problem.
516We have $f(\bx_0)=0$ so that this vertex is never updated
488during the iterations. 517during the iterations.
489The initial simplex is computed with a random number generator. 518The initial simplex is computed with a random number generator.
490The first vertex of the initial simplex is the origin. 519The first vertex of the initial simplex is the origin.
491The other vertices are uniform in the $[-1,1]$ interval. 520The other vertices are uniform in the $[-1,1]$ interval.
492for $i = 2, 2, \ldots , n+1$. 521An absolute termination criteria on the size of the simplex is used,
493An absolute termination criteria on the size of the simplex is used. 522that is, the algorithm is stopped when the inequality
494More precisely, the method was stopped when $\sigma(Sk) \leq 10^{-8}$ is satisfied. 523$\sigma_+(S_k) \leq 10^{-8}$ is satisfied.
495 524
496For this test, we compute the rate of convergence as presented 525For this test, we compute the rate of convergence as presented
497in Han \& Neuman \cite{HanNeumann2006}. This rate is defined as 526in Han \& Neuman \cite{HanNeumann2006}. This rate is defined as
@@ -502,7 +531,7 @@ in Han \& Neuman \cite{HanNeumann2006}. This rate is defined as
502\end{eqnarray} 531\end{eqnarray}
503where $k$ is the number of iterations. 532where $k$ is the number of iterations.
504That definition can be viewed as the geometric mean of the ratio of the 533That definition can be viewed as the geometric mean of the ratio of the
505oriented lengths between successive simplices and the minimizer 0. 534oriented lengths between successive simplices.
506This definition implies 535This definition implies
507\begin{eqnarray} 536\begin{eqnarray}
508\label{rho-sp-rate-convergence2} 537\label{rho-sp-rate-convergence2}
@@ -513,9 +542,10 @@ This definition implies
513If $k$ is the number of iterations required to obtain convergence, as 542If $k$ is the number of iterations required to obtain convergence, as
514indicated by the termination criteria, the rate of convergence is practically computed as 543indicated by the termination criteria, the rate of convergence is practically computed as
515\begin{eqnarray} 544\begin{eqnarray}
516\rho(S_0,n,k) = \left( \frac{\sigma(S_{k})}{\sigma(S_0)}\right)^{1/k} 545\rho(S_0,n,k) = \left( \frac{\sigma(S_{k})}{\sigma(S_0)}\right)^{1/k}.
517\end{eqnarray} 546\end{eqnarray}
518 547
548The following Scilab script allows to perform the optimization.
519 549
520\lstset{language=scilabscript} 550\lstset{language=scilabscript}
521\begin{lstlisting} 551\begin{lstlisting}
@@ -635,20 +665,20 @@ $n$ & \#Iterations & \# Reflections & \# Reflection & \#Shrink\\
635\end{tabular} 665\end{tabular}
636%\end{tiny} 666%\end{tiny}
637\end{center} 667\end{center}
638\caption{Numerical experiment with Spendley et al method on a generalized quadratic function -- number 668\caption{Numerical experiment with Spendley et al method on a
639and kinds of steps performed} 669generalized quadratic function -- Number of iterations
670and types of steps performed}
640\label{fig-sp-numexp3-nbsteps} 671\label{fig-sp-numexp3-nbsteps}
641\end{figure} 672\end{figure}
642 673
643The figure \ref{fig-sp-numexp3-dimension} presents the number of function 674The figure \ref{fig-sp-numexp3-dimension} presents the number of function
644evaluations depending on the number of variables. 675evaluations depending on the number of variables.
645 676We can see that the number of function evaluations
646One can see that the number of function evaluations
647increases approximately linearly with the dimension of the problem in 677increases approximately linearly with the dimension of the problem in
648figure \ref{fig-sp-numexp3-fvn}. A rough rule of thumb is that, for $n=1,20$, 678figure \ref{fig-sp-numexp3-fvn}. A rough rule of thumb is that, for $n=1,20$,
649the number of function evaluations is equal to $30n$. 679the number of function evaluations is equal to $30n$:
650This test is in fact the best that we can expect from this algorithm : since 680most iterations are shrink steps and approximately 30 iterations are required,
651most iterations are shrink steps, most iterations improves the function value. 681almost independently of $n$.
652 682
653\begin{figure}[htbp] 683\begin{figure}[htbp]
654\begin{center} 684\begin{center}
@@ -697,9 +727,8 @@ depending on the number of variables}
697 727
698The table \ref{fig-sp-numexp3-dimension} also shows the interesting 728The table \ref{fig-sp-numexp3-dimension} also shows the interesting
699fact that the convergence rate is almost constant and 729fact that the convergence rate is almost constant and
700very close to 1. This is a consequence of the shrink steps, 730very close to $1/2$. This is a consequence of the shrink steps,
701which are dividing the size of the simplex at each iterations by 2. 731which are dividing the size of the simplex at each iteration by 2.
702Therefore, the convergence rate is close to $1/2$.
703 732
704\section{Conclusion} 733\section{Conclusion}
705 734
@@ -713,7 +742,7 @@ functions and generates a very slow method in these
713cases. 742cases.
714 743
715In the last experiment, we have explored what happens when the number of 744In the last experiment, we have explored what happens when the number of
716iterations is increasing. The rate of convergence in this case is close to $1/2$ 745iterations is increasing. In this experiment, the rate of convergence is close to $1/2$
717and the number of function evaluations is a linear function of the number 746and the number of function evaluations is a linear function of the number
718of variables (approximately $30n$ in this case). 747of variables (approximately $30n$).
719 748
diff --git a/scilab_doc/neldermead/neldermead-spendley-so.pdf b/scilab_doc/neldermead/neldermead-spendley-so.pdf
index 19a77bc..01e42c8 100644
--- a/scilab_doc/neldermead/neldermead-spendley-so.pdf
+++ b/scilab_doc/neldermead/neldermead-spendley-so.pdf
Binary files differ
diff --git a/scilab_doc/neldermead/scripts/fminsearch_demo.m b/scilab_doc/neldermead/scripts/fminsearch_demo.m
new file mode 100644
index 0000000..8a8c200
--- /dev/null
+++ b/scilab_doc/neldermead/scripts/fminsearch_demo.m
@@ -0,0 +1,78 @@
1% Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2% Copyright (C) 2009 - Digiteo - 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% Basic use
11format long
12banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2;
13[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1])
14output.message
15
16% Use tolerance on X
17opt = optimset('TolX',1.e-2);
18[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
19
20% Use tolerance on F
21opt = optimset('TolFun',1.e-10);
22[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
23
24% Use MaxIter
25opt = optimset('MaxIter',10);
26[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
27
28% Use MaxFunEvals
29opt = optimset('MaxFunEvals',10);
30[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
31
32% Set Display to 'final'
33opt = optimset('Display','final');
34[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
35
36% Set Display to 'iter'
37opt = optimset('Display','iter');
38[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
39
40% Set Display to 'off'
41opt = optimset('Display','off');
42[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
43
44% Set Display to 'notify'
45opt = optimset('Display','notify');
46[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
47
48% Set Display to 'off', when there is a maximum number of iterations reached
49opt = optimset('Display','off' , 'MaxIter' , 10 );
50[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
51
52% Set Display to 'notify', when there is a maximum number of iterations reached
53opt = optimset('Display','notify' , 'MaxIter' , 10 );
54[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
55
56% Set Display to 'iter', when there is a maximum number of iterations reached
57opt = optimset('Display','iter' , 'MaxIter' , 10 );
58[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt );
59
60% Sample use of an output function
61options = optimset('OutputFcn', @outfun);
62hold on
63objfun=@(x) exp(x(1))*(4*x(1)^2+2*x(2)^2+x(1)*x(2)+2*x(2));
64[x fval] = fminsearch(objfun, [-1 1], options)
65hold off
66
67% Sample use of a plot function : optimplotfval
68options = optimset('PlotFcns',@optimplotfval);
69[x ffinal] = fminsearch(@onehump,[2,1],options)
70
71% Sample use of a plot function : optimplotfval
72options = optimset('PlotFcns',@optimplotx);
73[x ffinal] = fminsearch(@onehump,[2,1],options)
74
75% Sample use of a plot function : optimplotfunccount
76options = optimset('PlotFcns',@optimplotfunccount);
77[x ffinal] = fminsearch(@onehump,[2,1],options)
78
diff --git a/scilab_doc/neldermead/scripts/fminsearch_banana.m b/scilab_doc/neldermead/scripts/onehump.m
index d4dfeef..376cbe3 100644
--- a/scilab_doc/neldermead/scripts/fminsearch_banana.m
+++ b/scilab_doc/neldermead/scripts/onehump.m
@@ -7,8 +7,9 @@
7% are also available at 7% are also available at
8% http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt 8% http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9 9
10format long 10function f = onehump(x)
11banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2; 11r = x(1)^2 + x(2)^2;
12[x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1]) 12s = exp(-r);
13output.message 13f = x(1)*s+r/20;
14
14 15
diff --git a/scilab_doc/neldermead/scripts/outfun.m b/scilab_doc/neldermead/scripts/outfun.m
new file mode 100644
index 0000000..02dcd36
--- /dev/null
+++ b/scilab_doc/neldermead/scripts/outfun.m
@@ -0,0 +1,16 @@
1% Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2% Copyright (C) 2009 - Digiteo - 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% Use output function
11function stop = outfun(x, optimValues, state)
12stop = false;
13hold on;
14plot(x(1),x(2),'.');
15drawnow
16
diff --git a/scilab_doc/neldermead/testFminsearchPlotMatlab.png b/scilab_doc/neldermead/testFminsearchPlotMatlab.png
new file mode 100644
index 0000000..e78702c
--- /dev/null
+++ b/scilab_doc/neldermead/testFminsearchPlotMatlab.png
Binary files differ
diff --git a/scilab_doc/neldermead/testFminsearchPlotScilab.png b/scilab_doc/neldermead/testFminsearchPlotScilab.png
new file mode 100644
index 0000000..164d426
--- /dev/null
+++ b/scilab_doc/neldermead/testFminsearchPlotScilab.png
Binary files differ