Deleted vectorized computation feature. Deleted neldermead_contour. Fixed the demos.
[scilab.git] / scilab / modules / optimization / help / en_US / leastsq.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - INRIA
5  * 
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at    
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  -->
13 <refentry version="5.0-subset Scilab" xml:id="leastsq" xml:lang="en"
14           xmlns="http://docbook.org/ns/docbook"
15           xmlns:xlink="http://www.w3.org/1999/xlink"
16           xmlns:svg="http://www.w3.org/2000/svg"
17           xmlns:ns4="http://www.w3.org/1999/xhtml"
18           xmlns:mml="http://www.w3.org/1998/Math/MathML"
19           xmlns:db="http://docbook.org/ns/docbook">
20   <info>
21     <pubdate>$LastChangedDate$</pubdate>
22   </info>
23
24   <refnamediv>
25     <refname>leastsq</refname>
26
27     <refpurpose>Solves non-linear least squares problems</refpurpose>
28   </refnamediv>
29
30   <refsynopsisdiv>
31     <title>Calling Sequence</title>
32
33     <synopsis>[fopt,[xopt,[grdopt]]]=leastsq(fun, x0)
34 [fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, x0)
35 [fopt,[xopt,[grdopt]]]=leastsq(fun, cstr, x0)
36 [fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, cstr, x0)
37 [fopt,[xopt,[grdopt]]]=leastsq(fun, dfun, cstr, x0, algo)
38 [fopt,[xopt,[grdopt]]]=leastsq([imp], fun [,dfun] [,cstr],x0 [,algo],[df0,[mem]],[stop])</synopsis>
39   </refsynopsisdiv>
40
41   <refsection>
42     <title>Parameters</title>
43
44     <variablelist>
45       <varlistentry>
46         <term>fopt</term>
47
48         <listitem>
49           <para>value of the function <emphasis>f(x)=||fun(x)||^2</emphasis>
50           at <literal>xopt</literal></para>
51         </listitem>
52       </varlistentry>
53
54       <varlistentry>
55         <term>xopt</term>
56
57         <listitem>
58           <para>best value of <literal>x</literal> found to minimize
59           <emphasis>||fun(x)||^2</emphasis></para>
60         </listitem>
61       </varlistentry>
62
63       <varlistentry>
64         <term>grdopt</term>
65
66         <listitem>
67           <para>gradient of <emphasis>f</emphasis> at
68           <literal>xopt</literal></para>
69         </listitem>
70       </varlistentry>
71
72       <varlistentry>
73         <term>fun</term>
74
75         <listitem>
76           <para>a scilab function or a list defining a function from
77           <emphasis>R^n</emphasis> to <emphasis>R^m</emphasis> (see more
78           details in DESCRIPTION).</para>
79         </listitem>
80       </varlistentry>
81
82       <varlistentry>
83         <term>x0</term>
84
85         <listitem>
86           <para>real vector (initial guess of the variable to be
87           minimized).</para>
88         </listitem>
89       </varlistentry>
90
91       <varlistentry>
92         <term>dfun</term>
93
94         <listitem>
95           <para>a scilab function or a string defining the Jacobian matrix of
96           <literal>fun</literal> (see more details in DESCRIPTION).</para>
97         </listitem>
98       </varlistentry>
99
100       <varlistentry>
101         <term>cstr</term>
102
103         <listitem>
104           <para>bound constraints on <literal>x</literal>. They must be
105           introduced by the string keyword <literal>'b'</literal> followed by
106           the lower bound <literal>binf</literal> then by the upper bound
107           <literal>bsup</literal> (so <literal>cstr</literal> appears as
108           <literal>'b',binf,bsup</literal> in the calling sequence). Those
109           bounds are real vectors with same dimension than
110           <literal>x0</literal> (-%inf and +%inf may be used for dimension
111           which are unrestricted).</para>
112         </listitem>
113       </varlistentry>
114
115       <varlistentry>
116         <term>algo</term>
117
118         <listitem>
119           <para>a string with possible values: <literal>'qn'</literal> or
120           <literal>'gc'</literal> or <literal>'nd'</literal>. These strings
121           stand for quasi-Newton (default), conjugate gradient or
122           non-differentiable respectively. Note that <literal>'nd'</literal>
123           does not accept bounds on <literal>x</literal>.</para>
124         </listitem>
125       </varlistentry>
126
127       <varlistentry>
128         <term>imp</term>
129
130         <listitem>
131           <para>scalar argument used to set the trace mode.
132           <literal>imp=0</literal> nothing (except errors) is reported,
133           <literal>imp=1</literal> initial and final reports,
134           <literal>imp=2</literal> adds a report per iteration,
135           <literal>imp&gt;2</literal> add reports on linear search. Warning,
136           most of these reports are written on the Scilab standard
137           output.</para>
138         </listitem>
139       </varlistentry>
140
141       <varlistentry>
142         <term>df0</term>
143
144         <listitem>
145           <para>real scalar. Guessed decreasing of
146           <emphasis>||fun||^2</emphasis> at first iteration.
147           (<literal>df0=1</literal> is the default value).</para>
148         </listitem>
149       </varlistentry>
150
151       <varlistentry>
152         <term>mem</term>
153
154         <listitem>
155           <para>integer, number of variables used to approximate the Hessean
156           (second derivatives) of <emphasis>f</emphasis> when
157           <literal>algo</literal><literal>='qn'</literal>. Default value is
158           around 6.</para>
159         </listitem>
160       </varlistentry>
161
162       <varlistentry>
163         <term>stop</term>
164
165         <listitem>
166           <para>sequence of optional parameters controlling the convergence of
167           the algorithm. They are introduced by the keyword
168           <literal>'ar'</literal>, the sequence being of the form
169           <literal>'ar',nap, [iter [,epsg [,epsf [,epsx]]]]</literal></para>
170
171           <variablelist>
172             <varlistentry>
173               <term>nap</term>
174
175               <listitem>
176                 <para>maximum number of calls to <literal>fun</literal>
177                 allowed.</para>
178               </listitem>
179             </varlistentry>
180
181             <varlistentry>
182               <term>iter</term>
183
184               <listitem>
185                 <para>maximum number of iterations allowed.</para>
186               </listitem>
187             </varlistentry>
188
189             <varlistentry>
190               <term>epsg</term>
191
192               <listitem>
193                 <para>threshold on gradient norm.</para>
194               </listitem>
195             </varlistentry>
196
197             <varlistentry>
198               <term>epsf</term>
199
200               <listitem>
201                 <para>threshold controlling decreasing of
202                 <literal>f</literal></para>
203               </listitem>
204             </varlistentry>
205
206             <varlistentry>
207               <term>epsx</term>
208
209               <listitem>
210                 <para>threshold controlling variation of <literal>x</literal>.
211                 This vector (possibly matrix) of same size as
212                 <literal>x0</literal> can be used to scale
213                 <literal>x</literal>.</para>
214               </listitem>
215             </varlistentry>
216           </variablelist>
217         </listitem>
218       </varlistentry>
219     </variablelist>
220   </refsection>
221
222   <refsection>
223     <title>Description</title>
224
225     <para><emphasis>fun</emphasis> being a function from
226     <emphasis>R^n</emphasis> to <emphasis>R^m</emphasis> this routine tries to
227     minimize w.r.t. <emphasis>x</emphasis>, the function:</para>
228
229     <informalequation>
230       <mediaobject>
231         <imageobject>
232           <imagedata align="center" fileref="../mml/leastsq_equation_1.mml" />
233         </imageobject>
234       </mediaobject>
235     </informalequation>
236
237     <para>which is the sum of the squares of the components of
238     <emphasis>fun</emphasis>. Bound constraints may be imposed on
239     <literal>x</literal>.</para>
240   </refsection>
241
242   <refsection>
243     <title>How to provide fun and dfun</title>
244
245     <para><literal>fun</literal> can be either a usual scilab function (case
246     1) or a fortran or a C routine linked to scilab (case 2). For most
247     problems the definition of <emphasis>fun</emphasis> will need
248     supplementary parameters and this can be done in both cases.</para>
249
250     <variablelist>
251       <varlistentry>
252         <term>case 1:</term>
253
254         <listitem>
255           <para>when <literal>fun</literal> is a Scilab function, its calling
256           sequence must be: <literal>y=fun(x
257           [,opt_par1,opt_par2,...])</literal>. When <literal>fun</literal>
258           needs optional parameters it must appear as
259           <literal>list(fun,opt_par1,opt_par2,...)</literal> in the calling
260           sequence of <literal>leastsq</literal>.</para>
261         </listitem>
262       </varlistentry>
263
264       <varlistentry>
265         <term>case 2:</term>
266
267         <listitem>
268           <para>when <literal>fun</literal> is defined by a Fortran or C
269           routine it must appear as <literal>list(fun_name,m
270           [,opt_par1,opt_par2,...])</literal> in the calling sequence of
271           <literal>leastsq</literal>, <literal>fun_name</literal> (a string)
272           being the name of the routine which must be linked to Scilab (see
273           <link linkend="link">link</link>). The generic calling sequences for
274           this routine are:</para>
275
276           <programlisting role = ""><![CDATA[ 
277 In Fortran:    subroutine fun(m, n, x, params, y)
278                integer m,n
279                double precision x(n), params(*), y(m)
280
281 In C:          void fun(int *m, int *n, double *x, double *params, double *y)
282  ]]></programlisting>
283
284           <para>where <literal>n</literal> is the dimension of vector
285           <literal>x</literal>, <literal>m</literal> the dimension of vector
286           <literal>y</literal> (which must store the evaluation of
287           <emphasis>fun</emphasis> at <emphasis>x</emphasis>) and
288           <literal>params</literal> is a vector which contains the optional
289           parameters <emphasis>opt_par1, opt_par2, ...</emphasis> (each
290           parameter may be a vector, for instance if
291           <emphasis>opt_par1</emphasis> has 3 components, the description of
292           <emphasis>opt_par2</emphasis> begin from
293           <literal>params(4)</literal> (fortran case), and from
294           <literal>params[3]</literal> (C case), etc... Note that even if
295           <literal>fun</literal> doesn't need supplementary parameters you
296           must anyway write the fortran code with a <literal>params</literal>
297           argument (which is then unused in the subroutine core).</para>
298         </listitem>
299       </varlistentry>
300     </variablelist>
301
302     <para>In many cases it is adviced to provide the Jacobian matrix
303     <literal>dfun</literal> (<emphasis>dfun(i,j)=dfi/dxj</emphasis>) to the
304     optimizer (which uses a finite difference approximation otherwise) and as
305     for <literal>fun</literal> it may be given as a usual scilab function or
306     as a fortran or a C routine linked to scilab.</para>
307
308     <variablelist>
309       <varlistentry>
310         <term>case 1:</term>
311
312         <listitem>
313           <para>when <literal>dfun</literal> is a scilab function, its calling
314           sequence must be: <literal>y=dfun(x [, optional
315           parameters])</literal> (notes that even if <literal>dfun</literal>
316           needs optional parameters it must appear simply as
317           <literal>dfun</literal> in the calling sequence of
318           <literal>leastsq</literal>).</para>
319         </listitem>
320       </varlistentry>
321
322       <varlistentry>
323         <term>case 2:</term>
324
325         <listitem>
326           <para>when <literal>dfun</literal> is defined by a Fortran or C
327           routine it must appear as <literal>dfun_name</literal> (a string) in
328           the calling sequence of <literal>leastsq</literal>
329           (<literal>dfun_name</literal> being the name of the routine which
330           must be linked to Scilab). The calling sequences for this routine
331           are nearly the same than for <literal>fun</literal>:</para>
332
333           <programlisting role = ""><![CDATA[ 
334 In Fortran:    subroutine dfun(m, n, x, params, y)
335                integer m,n
336                double precision x(n), params(*), y(m,n)
337
338 In C:          void fun(int *m, int *n, double *x, double *params, double *y)
339  ]]></programlisting>
340
341           <para>in the C case <emphasis>dfun(i,j)=dfi/dxj</emphasis> must be
342           stored in <literal>y[m*(j-1)+i-1]</literal>.</para>
343         </listitem>
344       </varlistentry>
345     </variablelist>
346   </refsection>
347
348   <refsection>
349     <title>Remarks</title>
350
351     <para>Like <link linkend="datafit">datafit</link>,
352     <literal>leastsq</literal> is a front end onto the <link
353     linkend="optim">optim</link> function. If you want to try the
354     Levenberg-Marquard method instead, use <link
355     linkend="lsqrsolve">lsqrsolve</link>.</para>
356
357     <para>A least squares problem may be solved directly with the <link
358     linkend="optim">optim</link> function ; in this case the function <link
359     linkend="NDcost">NDcost</link> may be useful to compute the derivatives
360     (see the <link linkend="NDcost">NDcost</link> help page which provides a
361     simple example for parameters identification of a differential
362     equation).</para>
363   </refsection>
364
365   <refsection>
366     <title>Examples</title>
367
368     <programlisting role="example"><![CDATA[ 
369 // We will show different calling possibilities of leastsq on one (trivial) example
370 // which is non linear but doesn't really need to be solved with leastsq (applying
371 // log linearizes the model and the problem may be solved with linear algebra). 
372 // In this example we look for the 2 parameters x(1) and x(2) of a simple
373 // exponential decay model (x(1) being the unknow initial value and x(2) the
374 // decay constant): 
375
376 function y = yth(t, x)
377    y  = x(1)*exp(-x(2)*t) 
378 endfunction  
379
380 // we have the m measures (ti, yi):
381 m = 10;
382 tm = [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5]';
383 ym = [0.79, 0.59, 0.47, 0.36, 0.29, 0.23, 0.17, 0.15, 0.12, 0.08]';
384 wm = ones(m,1); // measure weights (here all equal to 1...)
385
386 // and we want to find the parameters x such that the model fits the given 
387 // datas in the least square sense:
388 // 
389 //  minimize  f(x) = sum_i  wm(i)^2 ( yth(tm(i),x) - ym(i) )^2   
390
391 // initial parameters guess
392 x0 = [1.5 ; 0.8];
393
394 // in the first examples, we define the function fun and dfun 
395 // in scilab language
396 function e = myfun(x, tm, ym, wm)
397    e = wm.*( yth(tm, x) - ym )
398 endfunction
399
400 function g = mydfun(x, tm, ym, wm)
401    v = wm.*exp(-x(2)*tm)
402    g = [v , -x(1)*tm.*v]
403 endfunction
404
405 // now we could call leastsq:
406
407 // 1- the simplest call
408 [f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),x0)
409
410 // 2- we provide the Jacobian
411 [f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),mydfun,x0)
412
413 // a small graphic (before showing other calling features)
414 tt = linspace(0,1.1*max(tm),100)';
415 yy = yth(tt, xopt);
416 xbasc()
417 plot2d(tm, ym, style=-2)
418 plot2d(tt, yy, style = 2)
419 legend(["measure points", "fitted curve"]);
420 xtitle("a simple fit with leastsq")
421
422 // 3- how to get some information (we use imp=1)
423 [f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0)
424
425 // 4- using the conjugate gradient (instead of quasi Newton)
426 [f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0,"gc")
427
428 // 5- how to provide bound constraints (not useful here !)
429 xinf = [-%inf,-%inf]; xsup = [%inf, %inf];
430 [f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),"b",xinf,xsup,x0) // without Jacobian
431 [f,xopt, gropt] = leastsq(list(myfun,tm,ym,wm),mydfun,"b",xinf,xsup,x0) // with Jacobian 
432
433 // 6- playing with some stopping parameters of the algorithm
434 //    (allows only 40 function calls, 8 iterations and set epsg=0.01, epsf=0.1)
435 [f,xopt, gropt] = leastsq(1,list(myfun,tm,ym,wm),mydfun,x0,"ar",40,8,0.01,0.1)
436
437
438 // 7 and 8: now we want to define fun and dfun in fortran then in C
439 //          Note that the "compile and link to scilab" method used here
440 //          is believed to be OS independant (but there are some requirements, 
441 //          in particular you need a C and a fortran compiler, and they must 
442 //          be compatible with the ones used to build your scilab binary).
443  
444 // 7- fun and dfun in fortran
445
446 // 7-1/ Let 's Scilab write the fortran code (in the TMPDIR directory):
447 f_code = ["      subroutine myfun(m,n,x,param,f)"
448           "*     param(i) = tm(i), param(m+i) = ym(i), param(2m+i) = wm(i)"
449           "      implicit none"
450           "      integer n,m"
451           "      double precision x(n), param(*), f(m)"
452           "      integer i"
453           "      do i = 1,m"
454           "         f(i) = param(2*m+i)*( x(1)*exp(-x(2)*param(i)) - param(m+i) )"
455           "      enddo"
456           "      end ! subroutine fun"
457           ""
458           "      subroutine mydfun(m,n,x,param,df)"
459           "*     param(i) = tm(i), param(m+i) = ym(i), param(2m+i) = wm(i)"
460           "      implicit none"
461           "      integer n,m"
462           "      double precision x(n), param(*), df(m,n)"
463           "      integer i"
464           "      do i = 1,m"
465           "         df(i,1) =  param(2*m+i)*exp(-x(2)*param(i))"
466           "         df(i,2) = -x(1)*param(i)*df(i,1)"
467           "      enddo"
468           "      end ! subroutine dfun"];
469
470 mputl(f_code,TMPDIR+'/myfun.f')
471
472 // 7-2/ compiles it. You need a fortran compiler !
473 names = ["myfun" "mydfun"]
474 flibname = ilib_for_link(names,"myfun.o",[],"f",TMPDIR+"/Makefile");
475
476 // 7-3/ link it to scilab (see link help page)
477 link(flibname,names,"f") 
478
479 // 7-4/ ready for the leastsq call: be carreful don't forget to
480 //      give the dimension m after the routine name !
481 [f,xopt, gropt] = leastsq(list("myfun",m,tm,ym,wm),x0)  // without Jacobian
482 [f,xopt, gropt] = leastsq(list("myfun",m,tm,ym,wm),"mydfun",x0) // with Jacobian
483
484 // 8- last example: fun and dfun in C
485
486 // 8-1/ Let 's Scilab write the C code (in the TMPDIR directory):
487 c_code = ["#include <math.h>"
488           "void myfunc(int *m,int *n, double *x, double *param, double *f)"
489           "{"
490           "  /*  param[i] = tm[i], param[m+i] = ym[i], param[2m+i] = wm[i] */"
491           "  int i;"
492           "  for ( i = 0 ; i &lt; *m ; i++ )"
493           "    f[i] = param[2*(*m)+i]*( x[0]*exp(-x[1]*param[i]) - param[(*m)+i] );"
494           "  return;"
495           "}"
496           ""
497           "void mydfunc(int *m,int *n, double *x, double *param, double *df)"
498           "{"
499           "  /*  param[i] = tm[i], param[m+i] = ym[i], param[2m+i] = wm[i] */"
500           "  int i;"
501           "  for ( i = 0 ; i &lt; *m ; i++ )"
502           "    {"
503           "      df[i] = param[2*(*m)+i]*exp(-x[1]*param[i]);"
504           "      df[i+(*m)] = -x[0]*param[i]*df[i];"
505           "    }"
506           "  return;"
507           "}"];
508
509 mputl(c_code,TMPDIR+'/myfunc.c')
510
511 // 8-2/ compiles it. You need a C compiler !
512 names = ["myfunc" "mydfunc"]
513 clibname = ilib_for_link(names,"myfunc.o",[],"c",TMPDIR+"/Makefile");
514
515 // 8-3/ link it to scilab (see link help page)
516 link(clibname,names,"c") 
517
518 // 8-4/ ready for the leastsq call
519 [f,xopt, gropt] = leastsq(list("myfunc",m,tm,ym,wm),"mydfunc",x0)
520  ]]></programlisting>
521   </refsection>
522
523   <refsection>
524     <title>See Also</title>
525
526     <simplelist type="inline">
527       <member><link linkend="lsqrsolve">lsqrsolve</link></member>
528
529       <member><link linkend="optim">optim</link></member>
530
531       <member><link linkend="NDcost">NDcost</link></member>
532
533       <member><link linkend="datafit">datafit</link></member>
534
535       <member><link linkend="external">external</link></member>
536
537       <member><link linkend="qpsolve">qpsolve</link></member>
538     </simplelist>
539   </refsection>
540 </refentry>