summaryrefslogtreecommitdiffstats log msg author committer range
blob: 9209dff992b7b2180466c8b56b237305f3e5854d (plain)
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610  \chapter{The \scifunction{fminsearch} function} In this chapter, we analyze the implementation of the \scifunction{fminsearch} which is provided in Scilab. In the first part, we describe the specific choices of this implementation with respect to the Nelder-Mead algorithm. In the second part, we present some numerical experiments which allows to check that the feature is behaving as expected, by comparison to Matlab's \scifunction{fminsearch}. \section{\scifunction{fminsearch}'s algorithm} \subsection{The algorithm} The algorithm used is the Nelder-Mead algorithm. This corresponds to the "variable" value of the "-method" option of the \scifunction{neldermead}. The "non greedy" version is used, that is, the expansion point is accepted only if it improves over the reflection point. \subsection{The initial simplex} The fminsearch algorithm uses a special initial simplex, which is an heuristic depending on the initial guess. The strategy chosen by fminsearch corresponds to the -simplex0method flag of the neldermead component, with the "pfeffer" method. It is associated with the - simplex0deltausual = 0.05 and -simplex0deltazero = 0.0075 parameters. Pfeffer's method is an heuristic which is presented in "Global Optimization Of Lennard-Jones Atomic Clusters" by Ellen Fan \cite{Fan2002}. It is due to L. Pfeffer at Stanford. See in the help of optimsimplex for more details. \subsection{The number of iterations} In this section, we present the default values for the number of iterations in fminsearch. The options input argument is an optionnal data structure which can contain the options.MaxIter field. It stores the maximum number of iterations. The default value is 200n, where n is the number of variables. The factor 200 has not been chosen by chance, but is the result of experiments performed against quadratic functions with increasing space dimension. This result is presented in "Effect of dimensionality on the nelder-mead simplex method" by Lixing Han and Michael Neumann \cite{HanNeumann2006}. This paper is based on Lixing Han's PhD, "Algorithms in Unconstrained Optimization" \cite{Han2000}. The study is based on numerical experiment with a quadratic function where the number of terms depends on the dimension of the space (i.e. the number of variables). Their study shows that the number of iterations required to reach the tolerance criteria is roughly 100n. Most iterations are based on inside contractions. Since each step of the Nelder-Mead algorithm only require one or two function evaluations, the number of required function evaluations in this experiment is also roughly 100n. \subsection{The termination criteria} The algorithm used by \scifunction{fminsearch} uses a particular termination criteria, based both on the absolute size of the simplex and the difference of the function values in the simplex. This termination criteria corresponds to the "-tolssizedeltafvmethod" termination criteria of the \scifunction{neldermead} component. The size of the simplex is computed with the $\sigma-+$ method, which corresponds to the "sigmaplus" method of the \scifunction{optimsimplex} component. The tolerance associated with this criteria is given by the "TolX" parameter of the \scifunction{options} data structure. Its default value is 1.e-4. The function value difference is the difference between the highest and the lowest function value in the simplex. The tolerance associated with this criteria is given by the "TolFun" parameter of the \scifunction{options} data structure. Its default value is 1.e-4. \section{Numerical experiments} In this section, we analyse the behaviour of Scilab's \scifunction{fminsearch} function, by comparison of Matlab's \scifunction{fminsearch}. We especially analyse the results of the optimization, so that we can check that the algorithm is indeed behaving the same way, even if the implementation is completely different. We consider the unconstrained optimization problem \cite{citeulike:1903787} \begin{eqnarray} \min f(\bx) \end{eqnarray} where $\bx\in\RR^2$ and the objective function $f$ is defined by \begin{eqnarray} f(\bx) = 100*(x_2-x_1^2)^2+(1-x_1)^2. \end{eqnarray} The initial guess is \begin{eqnarray} \bx^0 = ( -1.2 , 1.)^T, \end{eqnarray} where the function value is \begin{eqnarray} f(\bx^0) = 24.2. \end{eqnarray} The global solution of this problem is \begin{eqnarray} \bx^\star = ( 1 , 1.)^T \end{eqnarray} where the function value is \begin{eqnarray} f(\bx^\star) = 0. \end{eqnarray} \subsection{Algorithm and numerical precision} In this section, we are concerned by the comparison of the behavior of the two algorithms. We are going to check that the algorithms produces the same intermediate and final results. We also analyze the numerical precision of the results, by detailing the number of significant digits. To make a more living presentation of this topic, we will include small scripts which allow to produce the output that we are going to analyze. Because of the similarity of the languages, in order to avoid confusion, we will specify, for each script, the language we use by a small comment. Scripts and outputs written in Matlab's language will begin with \lstset{language=matlabscript} \begin{lstlisting} % Matlab % ... \end{lstlisting} while script written in Scilab's language will begin with \lstset{language=scilabscript} \begin{lstlisting} // Scilab // ... \end{lstlisting} The following Matlab script allows to see the behaviour of Matlab's \scifunction{fminsearch} function on Rosenbrock's test case. \lstset{language=matlabscript} \begin{lstlisting} % Matlab format long banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2; [x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1]) output.message \end{lstlisting} When this script is launched in Matlab, the following output is produced. \lstset{language=matlabscript} \begin{lstlisting} >> % Matlab >> format long >> banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2; >> [x,fval] = fminsearch(banana,[-1.2, 1]) >> [x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1]) x = 1.000022021783570 1.000042219751772 fval = 8.177661197416674e-10 exitflag = 1 output = iterations: 85 funcCount: 159 algorithm: 'Nelder-Mead simplex direct search' message: [1x194 char] >> output.message ans = Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 \end{lstlisting} The following Scilab script allows to solve the problem with Scilab's \scifunction{fminsearch}. \lstset{language=scilabscript} \begin{lstlisting} // Scilab format(25) function y = banana (x) y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2; endfunction [x , fval , exitflag , output] = fminsearch ( banana , [-1.2 1] ) output.message \end{lstlisting} The output associated with this Scilab script is the following. \lstset{language=scilabscript} \begin{lstlisting} -->// Scilab -->format(25) -->function y = banana (x) --> y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2; -->endfunction -->[x , fval , exitflag , output] = fminsearch ( banana , [-1.2 1] ) output = algorithm: "Nelder-Mead simplex direct search" funcCount: 159 iterations: 85 message: [3x1 string] exitflag = 1. fval = 0.0000000008177661099387 x = 1.0000220217835567027009 1.0000422197517710998227 -->output.message ans = !Optimization terminated: ! ! ! !the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 ! ! ! !and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 ! \end{lstlisting} Because the two softwares do not use the same formatting rules to produce their outputs, we must perform additionnal checking in order to check our results. The following Scilab script displays the results with 16 significant digits. \lstset{language=scilabscript} \begin{lstlisting} // Scilab // Print the result with 15 significant digits mprintf ( "%.15e" , fval ); mprintf ( "%.15e %.15e" , x(1) , x(2) ); \end{lstlisting} The previous script produces the following output. \lstset{language=scilabscript} \begin{lstlisting} -->// Scilab -->mprintf ( "%.15e" , fval ); 8.177661099387146e-010 -->mprintf ( "%.15e %.15e" , x(1) , x(2) ); 1.000022021783557e+000 1.000042219751771e+000 \end{lstlisting} These results are reproduced verbatim in the table \ref{fig-fminsearch-comparison}. \begin{figure}[htbp] \begin{center} \begin{tabular}{|l|l|l|} \hline Matlab Iterations & 85 &\\ Scilab Iterations & 85 &\\ \hline Matlab Function Evaluations & 159 &\\ Scilab Function Evaluations & 159 &\\ \hline Matlab $\bx^\star$ & 1.000022021783570 & 1.000042219751772 \\ Scilab $\bx^\star$ & 1.000022021783557e+000 & 1.000042219751771e+000 \\ \hline Matlab $f(\bx^\star)$ & 8.177661197416674e-10 &\\ Scilab $f(\bx^\star)$ & 8.177661099387146e-010 &\\ \hline \end{tabular} \end{center} \caption{Numerical experiment with Rosenbrock's function -- Comparison of results produced by Matlab and Scilab.} \label{fig-fminsearch-comparison} \end{figure} We must compute the common number of significant digits in order to check the consistency of the results. The following Scilab script computes the relative error between Scilab and Matlab results. \lstset{language=scilabscript} \begin{lstlisting} // Scilab // Compare the result xmb = [1.000022021783570 1.000042219751772 ]; err = norm(x - xmb) / norm(xmb); mprintf ( "Relative Error on x : %e\n", err ); fmb = 8.177661197416674e-10; err = abs(fval - fmb) / abs(fmb); mprintf ( "Relative Error on f : %e\n", err ); \end{lstlisting} The previous script produces the following output. \lstset{language=scilabscript} \begin{lstlisting} // Scilab Relative Error on x : 9.441163e-015 Relative Error on f : 1.198748e-008 \end{lstlisting} We must take into account for the floating point implementations of both Matlab and Scilab. In both these numerical softwares, double precision floating point numbers are used, i.e. the relative precision is both these softwares is $\epsilon \approx 10^{-16}$. That implies that there are approximately 16 significant digits. Therefore, the relative error on $x$, which is equivalent to 15 significant digits, is acceptable. If we now consider the relative error on $f$, which is equivalent to only 8 significant digits, that may sound as a problem. This corresponds to the square root of the relative precision, because $\sqrt{\epsilon} = \approx 10^{-8}$. In fact, this is the best that we can expect from an optimization algorithm (\cite{Brent73algorithmsfor,Gill81MurrayWright}). Therefore, the result is as close as possible to the result produced by Matlab. More specifically : \begin{itemize} \item the optimum $x$ is the same up to 15 significant digits, \item the function value at optimum is the same up to 8 significant digits, \item the number of iterations is the same, \item the number of function evaluations is the same, \item the exit flag is the same, \item the content of the output is the same (but the string is not display the same way). \end{itemize} The output of the two functions is the same. We must now check that the algorithms performs the same way, that is, produces the same intermediate steps. The following Matlab script allows to get deeper information by printing a message at each iteration with the "Display" option. \lstset{language=matlabscript} \begin{lstlisting} % Matlab opt = optimset('Display','iter'); [x,fval,exitflag,output] = fminsearch(banana,[-1.2, 1] , opt ); \end{lstlisting} The previous script produces the following output. \lstset{language=matlabscript} \begin{lstlisting} % Matlab Iteration Func-count min f(x) Procedure 0 1 24.2 1 3 20.05 initial simplex 2 5 5.1618 expand 3 7 4.4978 reflect 4 9 4.4978 contract outside 5 11 4.38136 contract inside 6 13 4.24527 contract inside 7 15 4.21762 reflect 8 17 4.21129 contract inside 9 19 4.13556 expand 10 21 4.13556 contract inside 11 23 4.01273 expand 12 25 3.93738 expand 13 27 3.60261 expand 14 28 3.60261 reflect 15 30 3.46622 reflect 16 32 3.21605 expand 17 34 3.16491 reflect 18 36 2.70687 expand 19 37 2.70687 reflect 20 39 2.00218 expand 21 41 2.00218 contract inside 22 43 2.00218 contract inside 23 45 1.81543 expand 24 47 1.73481 contract outside 25 49 1.31697 expand 26 50 1.31697 reflect 27 51 1.31697 reflect 28 53 1.1595 reflect 29 55 1.07674 contract inside 30 57 0.883492 reflect 31 59 0.883492 contract inside 32 61 0.669165 expand 33 63 0.669165 contract inside 34 64 0.669165 reflect 35 66 0.536729 reflect 36 68 0.536729 contract inside 37 70 0.423294 expand 38 72 0.423294 contract outside 39 74 0.398527 reflect 40 76 0.31447 expand 41 77 0.31447 reflect 42 79 0.190317 expand 43 81 0.190317 contract inside 44 82 0.190317 reflect 45 84 0.13696 reflect 46 86 0.13696 contract outside 47 88 0.113128 contract outside 48 90 0.11053 contract inside 49 92 0.10234 reflect 50 94 0.101184 contract inside 51 96 0.0794969 expand 52 97 0.0794969 reflect 53 98 0.0794969 reflect 54 100 0.0569294 expand 55 102 0.0569294 contract inside 56 104 0.0344855 expand 57 106 0.0179534 expand 58 108 0.0169469 contract outside 59 110 0.00401463 reflect 60 112 0.00401463 contract inside 61 113 0.00401463 reflect 62 115 0.000369954 reflect 63 117 0.000369954 contract inside 64 118 0.000369954 reflect 65 120 0.000369954 contract inside 66 122 5.90111e-005 contract outside 67 124 3.36682e-005 contract inside 68 126 3.36682e-005 contract outside 69 128 1.89159e-005 contract outside 70 130 8.46083e-006 contract inside 71 132 2.88255e-006 contract inside 72 133 2.88255e-006 reflect 73 135 7.48997e-007 contract inside 74 137 7.48997e-007 contract inside 75 139 6.20365e-007 contract inside 76 141 2.16919e-007 contract outside 77 143 1.00244e-007 contract inside 78 145 5.23487e-008 contract inside 79 147 5.03503e-008 contract inside 80 149 2.0043e-008 contract inside 81 151 1.12293e-009 contract inside 82 153 1.12293e-009 contract outside 83 155 1.12293e-009 contract inside 84 157 1.10755e-009 contract outside 85 159 8.17766e-010 contract inside Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 \end{lstlisting} The following Scilab script set the "Display" option to "iter" and run the \scifunction{fminsearch} function. \lstset{language=scilabscript} \begin{lstlisting} // Scilab opt = optimset ( "Display" , "iter" ); [x , fval , exitflag , output] = fminsearch ( banana , [-1.2 1] , opt ); \end{lstlisting} \lstset{language=scilabscript} \begin{lstlisting} // Scilab Iteration Func-count min f(x) Procedure 0 3 24.2 1 3 20.05 initial simplex 2 5 5.161796 expand 3 7 4.497796 reflect 4 9 4.497796 contract outside 5 11 4.3813601 contract inside 6 13 4.2452728 contract inside 7 15 4.2176247 reflect 8 17 4.2112906 contract inside 9 19 4.1355598 expand 10 21 4.1355598 contract inside 11 23 4.0127268 expand 12 25 3.9373812 expand 13 27 3.602606 expand 14 28 3.602606 reflect 15 30 3.4662211 reflect 16 32 3.2160547 expand 17 34 3.1649126 reflect 18 36 2.7068692 expand 19 37 2.7068692 reflect 20 39 2.0021824 expand 21 41 2.0021824 contract inside 22 43 2.0021824 contract inside 23 45 1.8154337 expand 24 47 1.7348144 contract outside 25 49 1.3169723 expand 26 50 1.3169723 reflect 27 51 1.3169723 reflect 28 53 1.1595038 reflect 29 55 1.0767387 contract inside 30 57 0.8834921 reflect 31 59 0.8834921 contract inside 32 61 0.6691654 expand 33 63 0.6691654 contract inside 34 64 0.6691654 reflect 35 66 0.5367289 reflect 36 68 0.5367289 contract inside 37 70 0.4232940 expand 38 72 0.4232940 contract outside 39 74 0.3985272 reflect 40 76 0.3144704 expand 41 77 0.3144704 reflect 42 79 0.1903167 expand 43 81 0.1903167 contract inside 44 82 0.1903167 reflect 45 84 0.1369602 reflect 46 86 0.1369602 contract outside 47 88 0.1131281 contract outside 48 90 0.1105304 contract inside 49 92 0.1023402 reflect 50 94 0.1011837 contract inside 51 96 0.0794969 expand 52 97 0.0794969 reflect 53 98 0.0794969 reflect 54 100 0.0569294 expand 55 102 0.0569294 contract inside 56 104 0.0344855 expand 57 106 0.0179534 expand 58 108 0.0169469 contract outside 59 110 0.0040146 reflect 60 112 0.0040146 contract inside 61 113 0.0040146 reflect 62 115 0.0003700 reflect 63 117 0.0003700 contract inside 64 118 0.0003700 reflect 65 120 0.0003700 contract inside 66 122 0.0000590 contract outside 67 124 0.0000337 contract inside 68 126 0.0000337 contract outside 69 128 0.0000189 contract outside 70 130 0.0000085 contract inside 71 132 0.0000029 contract inside 72 133 0.0000029 reflect 73 135 0.0000007 contract inside 74 137 0.0000007 contract inside 75 139 0.0000006 contract inside 76 141 0.0000002 contract outside 77 143 0.0000001 contract inside 78 145 5.235D-08 contract inside 79 147 5.035D-08 contract inside 80 149 2.004D-08 contract inside 81 151 1.123D-09 contract inside 82 153 1.123D-09 contract outside 83 155 1.123D-09 contract inside 84 157 1.108D-09 contract outside 85 159 8.178D-10 contract inside Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 \end{lstlisting} A close inspection at the data reveals that the two softwares produces indeed the same intermediate results. \subsection{Plot features} In this section, we check that the plotting features of the \scifunction{fminsearch} function are the same. The following output function plots in the current graphic window the value of the current parameter $\bx$. To let Matlab load that script, save the content in a .m file, in a directory known by Matlab. \lstset{language=matlabscript} \begin{lstlisting} % Matlab function stop = outfun(x, optimValues, state) stop = false; hold on; plot(x(1),x(2),'.'); drawnow \end{lstlisting} The following Matlab script allows to perform the optimization so that the output function is called back at each iteration. \lstset{language=matlabscript} \begin{lstlisting} % Matlab options = optimset('OutputFcn', @outfun); [x fval] = fminsearch(banana, [-1.2, 1], options) \end{lstlisting} This produces the plot which is presented in figure \ref{fig-fminsearch-matlab-outputfun}. \begin{figure} \begin{center} \includegraphics[width=10cm]{testFminsearchPlotMatlab.png} \end{center} \caption{Plot produced by Matlab's fminsearch, with customized output function.} \label{fig-fminsearch-matlab-outputfun} \end{figure} The following Scilab script sets the "OutputFcn" option and then calls the \scifunction{fminsearch} in order to perform the optimization. \lstset{language=scilabscript} \begin{lstlisting} // Scilab function outfun ( x , optimValues , state ) plot( x(1),x(2),'.'); endfunction opt = optimset ( "OutputFcn" , outfun); [x fval] = fminsearch ( banana , [-1.2 1] , opt ); \end{lstlisting} The previous script produces the plot which is presented in figure \ref{fig-fminsearch-scilab-outputfun}. \begin{figure} \begin{center} \includegraphics[width=10cm]{testFminsearchPlotScilab.png} \end{center} \caption{Plot produced by Scilab's fminsearch, with customized output function.} \label{fig-fminsearch-scilab-outputfun} \end{figure} Except for the size of the dots (which can be configured in both softwares), the graphics are exactly the same.