Code covered by the BSD License  

Highlights from
OpenOpt

image thumbnail

OpenOpt

by

 

25 Nov 2006 (Updated )

nonSmoothSolve (similar to fsolve), non-smooth & noisy local + some global solvers; works in Octave

nonSmoothSolveEx.m
clc
disp('----------------------------------')
%opts = optimset(@fsolve);
opts.TolFun = 1e-3;
opts.MaxIter = 1e4;
opts.MaxFunEvals = 1e5;
opts.iterPrint = -1;%no output for nonSmoothSolve

N = 40;% fsolve() fails even with N=2 for MATLAB & 4 for Octave
noise = 1e-5;

% warning! if noise is big modifying of TolFun, diffInt is recomended
disp(['noisy: N = ' int2str(N) '; noise = ' num2str(noise)])
A = sin(pascal(N)) + N*speye(N);
b = A*ones(N,1);

ff = @(x) A*x-b;
Noise = @(x) noise*rand(length(x),1);
FF = @(x) ff(x)+Noise(x);
 
x0 = 15*cos(1:N).';
disp('fsolve:')
tic
[x fval exitflag output] = fsolve(FF, x0, opts);
toc
disp(['fsolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval without noise = '])
disp(norm(ff(x), inf)) % max(abs(A*x-b))
disp('nonSmoothSolve:')
tic
[x fval exitflag output] = nonSmoothSolve(ff, x0, opts);
toc
% other way for nonSmoothSolve call
% prob = ooAssign(FF, x0, 'iterPrint=-1');
% tic
% r = ooRun(prob, @nonSmoothSolve);
% toc
% x = r.xf;

disp(['nonSmoothSolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval without noise = '])
disp(norm(ff(x), inf)) % max(abs(A*x-b))


disp('----------------------------------')
% example 2: non-smooth
N = 15; 
A = 1.25.^(1./hilb(N));

b = A*ones(N,1);

ff = @(x) A*(200*abs(x-0.005))-b;
%non-smooth, convex

x0 = N*cos(1:N).';

tic
[x fval exitflag output] = fsolve(ff, x0, opts);
toc
disp(['fsolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval without noise = '])
disp(norm(ff(x), inf)) % max(abs(A*x-b))


opts.iterPrint = 100;
tic
[x fval exitflag output] = nonSmoothSolve(ff, x0, opts);
toc
disp(['nonSmoothSolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval without noise = '])
disp(norm(ff(x), inf)) % max(abs(A*x-b))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%other ways (more flexible) for nonSmoothSolve calling
disp('----------------------------------')
N = 10; 
A = 1.25.^(1./hilb(N));

b = A*ones(N,1);

ff = @(x) A*(200*abs(x-0.005))-b;
%non-smooth, convex

x0 = N*cos(1:N).';

prob = nonSmoothAssign(ff, x0, 'iterPrint=100;MaxIter=1e4;MaxFunEvals=1e5');
prob.nonSmoothSolve.solver = @ralg;
tic
r = ooRun(prob, @nonSmoothSolve);
toc

%suppose we somehow know radius of sphere with center x0
% that certainly containes x* 
% then we can try using ShorEllipsoid instead of default ralg
prob.r0 = norm(N*cos(1:N).' - ones(N,1));
prob.nonSmoothSolve.solver = @ShorEllipsoid;
% but ShorEllipsoid is good only for small-scale problems, 
% with nVars<=10
disp('----------ShorEllipsoid---------')
tic
r = ooRun(prob, @nonSmoothSolve);
toc




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% ill-conditioned
% (uncomment this one if you want)
% N = 100;
% A = hilb(N);
% disp(['ill-conditioned: N = ' int2str(N) ' cond(A) = ' sprintf('%.2g',cond(A))])
% b = A*ones(N,1);
% ff = @(x) A*x-b;
% x0 = sqrt(N)*sin(1:N).';
% % % % % % fsolve also handles ill-conditioned well
% % % % % % disp('fsolve:')
% % % % % % tic
% % % % % % [x fval exitflag output] = fsolve(ff, x0, opts);
% % % % % % toc
% % % % % % disp(['fsolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval = '])
% % % % % % disp(norm(A*x-b, inf)) % max(abs(A*x-b))
% disp('nonSmoothSolve:')
% tic
% [x fval exitflag output] = nonSmoothSolve(ff, x0, opts);
% toc
% disp(['nonSmoothSolve: funccount= ' int2str(output.funcCount) '; exitflag = ' int2str(exitflag) '; fval = '])
% disp(norm(A*x-b, inf)) % max(abs(A*x-b))


% example 1 (noisy): 
% objfun(x) = (Ax-b) + random noise
% noise is often appear from numerical integration, solving ODE systems,
% other long difficult calculations where numerical errors are growing up

Contact us