Using a Gradient Method
>function map %fv (t; f, x, v) := f(x+t*v);
>function vmin (f,x,v,tstart,tmax=1,tmin=epsilon) ...
## Minimize f(x+t*v)
y=f(x);
t=tstart;
repeat
ynew=f(x+t*v);
while ynew<y;
t=t*2;
until t>tmax;
end;
t=fmin("%fv",0,t;f,x,v);
return {x+t*v,t};
endfunction
>function f([x,y]) &= x^2+2*y^2
2 2
2 y + x
>function df([x,y]) &= gradient(f(x,y),[x,y])
[2 x, 4 y]

>x0=[1,1]; {x1,t}=vmin("f",x0,-df(x0),0.1); x1,
[ 0.444444445171 -0.111111109657 ]
>{x2,t}=vmin("f",x1,-df(x1),t); x2,
[ 0.0740740740087 0.0740740731982 ]
>X=(x0_x1_x2)'; ...
> plot2d("f",r=2,niveau=[f(x0),f(x1),f(x2)]); ...
> plot2d(X[1],X[2],color=red,>add); insimg;

>X
1 0.444444445171 0.0740740740087
1 -0.111111109657 0.0740740731982

>function newtonvf([x,y]) &= -df(x,y).invert(hessian(f(x,y),[x,y]))
[ - x - y ]
>x0=[1,1]; x1=vmin("f",x0,newtonvf(x0),0.1)
[ 1.03672626039e-012 1.03672626039e-012 ]
>x0+newtonvf(x0)
[ 0 0 ]
>function f([x,y]) &= x^2+y^4
4 2
y + x
>function df([x,y]) &= gradient(f(x,y),[x,y])
3
[2 x, 4 y ]
>x=[1,1]; t=0.1; X=x; ...
> loop 1 to 100; {x,t}=vmin("f",x,-df(x),t); X=X_x; end; ...
> X=X';
>plot2d("f",r=1.2,niveau=fmap(X[1,1:3],X[2,1:3])); ...
> plot2d(X[1],X[2],>add,color=red); insimg;

>X[:,-10:-1]
Column 1 to 3:
0.000120628156681 -5.90109250917e-005 0.000116746853007
-0.0348745189166 -0.0347481888263 -0.0344982644339
Column 4 to 6:
-5.7144026496e-005 0.000113084545099 -5.53661693139e-005
-0.0343759568583 -0.0341339342155 -0.034015450507
Column 7 to 9:
0.000109580698827 -5.3693072045e-005 0.000106283681283
-0.0337809421765 -0.0336660666202 -0.033438690026
Column 10:
-5.20882874773e-005
-0.0333272632694

>k=5:cols(X); plot2d(log(fmap(X[1,k],X[2,k]))/k); insimg;

>function newtonvf([x,y]) &= -df(x,y).invert(hessian(f(x,y),[x,y]))
[ y ]
[ - x - - ]
[ 3 ]
>x=[1,1]; t=0.1; X=x; ...
> loop 1 to 20; {x,t}=vmin("f",x,newtonvf(x),t); X=X_x; end; ...
> X=X';
>plot2d("f",r=1.2,niveau=fmap(X[1,1:3],X[2,1:3])); ...
> plot2d(X[1],X[2],>add,color=red); insimg;

>X[:,-1]
5.66579793109e-011
9.612814048e-006
>k=3:cols(X); plot2d(log(fmap(X[1,k],X[2,k]))/k); insimg;

>x=[1,1]; X=x; ...
> loop 1 to 20; x=x+newtonvf(x); X=X_x; end; ...
> X=X';
>plot2d("f",r=1.2,niveau=fmap(X[1,1:3],X[2,1:3])); ...
> plot2d(X[1],X[2],>add,color=red); insimg;

>X[:,-1]
0
0.000300728659822
>k=3:cols(X); plot2d(log(fmap(X[1,k],X[2,k]))/k); insimg;

>function f([x,y]) &= 100*(y-x^2)^2+(1-x)^2;

>plot2d("f",r=3,>contour); plot2d(1,1,>points,>add); insimg;

>function df([x,y]) &= gradient(f(x,y),[x,y])
2 2
[- 400 x (y - x ) - 2 (1 - x), 200 (y - x )]
>x=[0,0]; t=0.1; X=x; ...
> loop 1 to 100; {x,t}=vmin("f",x,-df(x),t); X=X_x; end; ...
> X=X';
>plot2d("f",r=1.2,>contour); plot2d(1,1,>points,>add); ...
> plot2d(X[1],X[2],>add,color=red); insimg;

>function newtonvf([x,y]) &= -df(x,y).invert(hessian(f(x,y),[x,y]));
>x=[0,0]; t=0.1; X=x; ...
> loop 1 to 20; {x,t}=vmin("f",x,newtonvf(x),t); X=X_x; end; ...
> X=X';
>plot2d("f",r=1.2,>contour); plot2d(1,1,>points,>add); ...
> plot2d(X[1],X[2],>add,color=red); insimg;

>X[:,-1]
1
1
>load "Nelder Mead Demo"
>function f([x,y]) &= 100*(y-x^2)^2+(1-x)^2
2 2 2
100 (y - x ) + (1 - x)
>&solve(gradient(f(x,y),[x,y]))
[[y = 1, x = 1]]
>plot3d("f",r=3,>contour,angle=160°); insimg;

>plot2d("f",r=3,>contour); plot2d(1,1,>add,>points); insimg;

>s=[0,1;1,0;0,0];
>nelderdemo(s,"f"), insimg;
0.985415224543 0.969003877292
0.966540627635 0.934359321519
0.987483912325 0.97603469899

>nelder("f",[0,0])
[ 0.999999942023 0.999999814945 ]
>type vmin
function vmin (f, x, v, tstart, tmax, tmin)
## Default for tmax : 1
## Default for tmin : 2.22044604925e-012
y=f(x);
t=tstart;
repeat
ynew=f(x+t*v);
while ynew<y;
t=t*2;
until t>tmax;
end;
t=fmin("%fv",0,t;f,x,v);
return {x+t*v,t};
endfunction
>type %fv
function map %fv (t, f, x, v)
useglobal; return f(x+t*v);
endfunction
>function f([x,y]) := x+y^2
>f([4,5]), f(4,5)
29
29

>plot2d("%fv",0,2;"f",[-1,-1],[1,1],>insimg);

>integrate("%fv",0,2;"f",[-1,-1],[1,1])
0.666666666667
>fmin("%fv",0,2;"f",[-1,-1],[1,1])
0.499999994731
Examples Homepage