Skip to content
Snippets Groups Projects
Commit bd4f1101 authored by tsal's avatar tsal
Browse files

Adding v2.2 as from homepage

parents
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
Showing
with 2407 additions and 0 deletions
Manual.pdf 0 → 100644
+ 0
− 0
View file @ bd4f1101
File added
readme.html 0 → 100644
+ 72
− 0
View file @ bd4f1101
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML>
<HEAD>
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso-8859-1">
<TITLE>immoptibox</TITLE>
</HEAD>
<BODY BACKGROUND="bar.gif" TEXT="#000000" BGCOLOR="#ffffff" LINK="#0000EE" VLINK="#551A8B" ALINK="#FF0000">
<center>
<table BORDER="0" CELLSPACING="0" CELLPADDING="0">
<tr> <td valign="top" width="600">
<IMG SRC="pixel.gif" ALT="" WIDTH="600" HEIGHT="1" VSPACE="1" HSPACE="0"> </td></tr>
<tr> <td valign="top" width="600">
<table BORDER="0" CELLSPACING="0" CELLPADDING="0">
<tr><td valign="top" width="50">
<IMG SRC="pixel.gif" ALT="" WIDTH="50" HEIGHT="1" VSPACE="20" HSPACE="0"><br>
<center>
<IMG SRC="https://archive.compute.dtu.dk/swarm/img/dtu.jpg" ALT="DTU"
WIDTH="35" border="0"> </center> </td>
<td valign="top" WIDTH="5">
<IMG SRC="pixel.gif" ALT="" WIDTH="5" HEIGHT="1" HSPACE="0">
</td>
<td valign="top" WIDTH="545">
<IMG SRC="pixel.gif" ALT="" WIDTH="545" HEIGHT="1" HSPACE="0">
<center>
<h1><tt>immoptibox</tt></h1>
<h2>A MATLAB TOOLBOX FOR OPTIMIZATION AND DATA FITTING</h2>
<p>
<h3> Hans Bruun Nielsen <br>
<a href="http://www.imm.dtu.dk">Informatics and Mathematical Modelling</a><br>
<a href="http://www.dtu.dk">Technical University of Denmark</a> </h3>
</center>
</table> </td></tr>
<p><p>
<tr><td width="600" align="center">
<font size=+2 color="red">Version 2.2. &nbsp; November 2010 </font>
<!-- <br> &nbsp; <br>
<blink> <img src="greenbal.gif"> </blink> &nbsp;
<font size=+3 color="green">
The site is under construction </font> </td> -->
</table>
<p> &nbsp; <p>
<a name="TOP"></a>
<table BORDER="0" CELLSPACING="0" CELLPADDING="0">
<tr> <td valign="top" width="600">
<IMG SRC="pixel.gif" ALT="" WIDTH="595" HEIGHT="1" VSPACE="1" HSPACE="5"> </td>
<tr> <td width="600">
The toolbox
contains a number of functions for optimization an data fitting. The selection
of algorithms was guided by courses
in optimization,
but the M<font size=-1>ATLAB</font> functions in the toolbox are expected also
to have wider interest.
<p>
To get the toolbox, download <tt>
<a href="immoptibox.zip">immoptibox.zip</a></tt>
to the directory, where you save your M<font size=-1>ATLAB</font> files, use
<b>unzip</b> to unpack it, and update your M<font size=-1>ATLAB</font> path.
<p>
A pdf version of the manual can be downloaded from
<a href="Manual.pdf">here</a>.
</table>
</BODY>
</HTML>
src/Contents.m 0 → 100644
+ 62
− 0
View file @ bd4f1101
% IMMOPTIBOX Optimization and Data Fitting Toolbox
% Version 2.2, November 19, 2010
% Copyright (c) 2010 by Hans Bruun Nielsen and IMM.
%
% General unconstrained optimization
% dampnewton - Damped Neton method. Demands analytical expreesions
% for the gradient and Hessian
% linesearch - Choice between soft and exact line search.
% ucminf - Unconstrained minimization of a scalar function of a
% vector variable. Demands analytical expression for
% the gradient. Based on BFGS updating of the inverse
% Hessian and soft line search.
%
% Unconstrained, nonlinear least squares problems
% dogleg - Powell's dog-leg method. Demands analytical
% expression for the Jacobian.
% marquardt - Levenberg-Marquardt method. Demands analytical
% expression for the Jacobian.
% smarquardt - Levenberg-Marquardt method with successive updating
% of approximations to the Jacobian.
%
% Data fitting with cubic splines
% splinefit - Weighted least squares fit of a cubic spline to
% given data points. Possibility of assigning boundary
% conditions.
% splineval - Evaluate a cubic spline s as computed by SPLINEFIT.
% splinedif - Evaluate s', s'' or s'''.
%
% Robust estimation
% huberobj - Value and gradient of Huber estimator.
% Allows one-sided Huber function.
% linhuber - Minimizer of an extended linear Huber estimation
% problem. Allows one-sided Huber function.
% nonlinhuber - Minimizer of an extended nonlinear Huber estimation
% problem. Allows one-sided Huber function.
%
% Multiexponential fitting
% mexpfit - Weighted least squares fit of a multiexponential
% model to given data points. Algorithm based on
% separability.
%
% Nonlinear systems of equations
% nonlinsys - Solve nonlinear system of equations. Dog Leg
% method with updating of approximate Jacobian.
%
% Auxiliary programs
% checkgrad - Check implementation of gradient (or Jacobian)
% by means of finite differences.
%
% Test problems
% uctpget - Define test problem for unconstrained minimization.
% uctpval - Evaluate test problem.
%
% Data files
% optic.dat - Optic fibre data.
% osl.dat - Data from an optically stimulated luminiscence
% experiment.
% peaks.dat - Data with peaks and "shoulders".
% wild.dat - Data with "wild points".
% efit1.dat - Data for exponential fitting.
% efit2.dat - As efit1.dat, except that there are "wild points".
\ No newline at end of file
src/checkgrad.m 0 → 100644
+ 96
− 0
View file @ bd4f1101
function [maxJ, err, index] = checkgrad(fun, x, h, varargin)
%CHECKGRAD : Check the user's implementation of a nonlinear
% (vector) function and its gradient (Jacobian).
% The implementation must be a MATLAB function with declaration
% function [r, J] = fun(x,p1,p2,...)
% p1,p2,... are parameters of the function. In connection with
% nonlinear data fitting they may be arrays with coordinates of
% the data points. r and J should be the vector function and
% its Jacobian evaluated at x .
%
% Call [maxJ, err, index] = checkgrad(fun, x, h, p1,p2,...)
%
% Input parameters
% fun : Handle to the function.
% x : The point where we wish to check the Jacobian.
% h : Steplength used in difference approximation to J .
% May be a vector with h(j) used in the j'th coordinate
% direction.
% p1,p2,... are passed dirctly to the function FUN .
%
% Output parameters
% maxJ : The maximal absolute value of the elements in J .
% err : Vector with three elements. The maximal absolute deviation
% between an element in J and the corresponding difference
% approximation, computed by
% err(1): Forward difference with steplength h
% err(2): Backward difference with steplength h/2
% err(3): "Extrapolated difference", see below.
% index : 3*2 array, with index(k,:) giving the position in J
% where err(k) occurs.
%
% If the function is "smooth" and h is sufficiently small (but not
% so small that rounding errors dominate), then we can expect
% error(1) "=" A*h, error(2) "=" -A*(h/2)
% where "=" means "approximately equal to" and A is a constant, that
% depends on the function and x , but not on h . The "extrapolated
% approximation" is
% (forward approximation + 2(backward approximation))/3
% and its error is of the order of magnitude error(2)^2 .
%
% The algorithm is further discussed in "Checking Gradients", which
% can be downloaded from http://www.imm.dtu.dk/~hbn/Software/#AUX
% Version 10.10.19. hbn(a)imm.dtu.dk
% Initialize
sx = size(x); n = max(sx);
if min(sx) ~= 1, error('x must be a vector'), end
sh = size(h); nh = max(sh);
if nh > 1 & (min(sh) ~= 1 | nh ~= n)
error('h must be a scalar or a vector of the same length as x')
end
if nh == 1, h = repmat(h,1,n); end
ieq = find(x(:)+h(:) == x(:) | x(:)-h(:)/2 == x(:));
if ~isempty(ieq)
error(['h was too small in direction(s) ' int2str(ieq')])
end
% Check call
[f J] = feval(fun, x, varargin{:}); % offset point
sf = size(f); m = max(sf);
if min(sf) ~= 1, error('f must be a vector'), end
sJ = size(J);
if m == 1 % Scalar case. Allow both row and column
ok = (min(sJ) == 1) & (max(sJ) == n);
J = J(:).';
else % Vector case. Insist on m*n
ok = norm([m n] - sJ) == 0;
end
if ~ok
error('The sizes of x , f and J do not match'), end
% Check results
err = zeros(3,1); index = zeros(3,2);
maxJ = max(max(abs(J))); df = zeros(m,3);
% loop through the coordinate directions.
for j = 1 : n
% Make differences for all function components
xx = x; xx(j) = x(j) + h(j); % forward point
dx = xx(j) - x(j); % true difference in j'th direction
ff = feval(fun, xx, varargin{:});
df(:,1) = (ff(:) - f(:))/dx; % forward approximation
xx(j) = x(j) - h(j)/2; dx = x(j) - xx(j); % backward point
fb = feval(fun, xx, varargin{:});
df(:,2) = (f(:) - fb(:))/dx; % backward approximation
df(:,3) = (2*df(:,2) + df(:,1))/3; % extrapolated approximation
% Store errors
for k = 1 : 3
dif = df(:,k) - J(:,j);
[md, i] = max(abs(dif)); % largest difference and its index
if md > abs(err(k)) % New maximum
err(k) = dif(i); index(k,:) = [i j];
end
end;
end % j-loop
src/dampnewton.m 0 → 100644
+ 130
− 0
View file @ bd4f1101
function [X, info, perf] = dampnewton(fun, x0, opts, varargin)
%DAMPNEWTON Levenberg-Marquardt damped Newton method for unconstrained
% optimization: Find xm = argmin{f(x)} , where x is an n-vector and the
% scalar function f with gradient g (with elements g_i = Df/Dx_i )
% and Hessian H (with elements H_{j,k} = D^2f/(Dx_j Dx_k) ) must be
% given by a Matlab function with with declaration
% function [f, g, H] = fun(x, p1,p2,...)
% p1,p2,... are parameters of the function.
%
% Call [X, info] = dampnewton(fun, x0)
% [X, info] = dampnewton(fun, x0, opts, p1,p2,...)
% [X, info, perf] = dampnewton(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for xm .
% opts : Either a struct with fields 'tau', 'tolg', 'tolx' and 'maxeval'
% or a vector with the values of these options,
% opts = [tau tolg tolx maxeval].
% tau is used in starting value for damping parameter:
% mu = tau) * ||H0||_inf, where H0 is the Hessian of f at x0.
% The other options are used in stopping criteria:
% ||g(x)||_inf <= tolg or
% ||dx||2 <= tolx*(tolx + ||x||_2) or
% no. of function evaluations exceeds maxeval .
% Default tau = 1e-3, tolg = 1e-4, tolx = 1e-8, maxeval = 100.
% If the input opts has less than 4 elements, it is augmented by
% the default values. Also, zeros and negative elements are
% replaced by the default values.
% p1,p2,.. are passed directly to the function FUN .
%
% Output parameters
% X : If perf is present, then array, holding the iterates xk
% columnwise. Otherwise, computed solution vector.
% info : Performance information, vector with 7 elements:
% info(1:4) = final values of
% [f(x) ||g||_inf ||dx||_2 mu/||H(x)||_inf] .
% info(5:6) = no. of iteration steps and function evaluations.
% info(7) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : No. of function evaluations exceeded
% -1 : x is not a real valued vector
% -2 : f is not a real valued scalar
% -3 : g is not a real valued vector or
% H is not a real valued matrix
% -4 : Dimension mismatch in x, g, H
% -5 : H is not symmetric
% perf : Struct with fields
% f : values of f(xk) ,
% ng : values of || g(xk) ||_inf ,
% mu : values of the damping mu(xk) .
% Version 10.10.08. hbn(a)imm.dtu.dk
% Check parameters and function call
if nargin < 2, stop = -1;
else
[stop x n] = checkx(x0);
if ~stop
[stop f g H] = checkfgH(fun,x0,varargin{:});
end
end
if stop
X = x0; perf = []; info = [repmat(NaN,1,4) 0 1 stop];
return
end
if nargin < 3, opts = []; end
opts = checkopts('dampnewton', opts); % use default options where required
tau = opts(1); tolg = opts(2); tolx = opts(3); maxeval = opts(4);
% Finish initialization
ng = norm(g,inf); mu = tau * norm(H,inf);
if mu == 0 % zero initial Hessian. Steepest descent direction
mu = ng / max(norm(x), sqrt(eps));
end
Trace = nargout > 2;
if Trace
o = ones(1, maxeval);
X = x * o; perf = [f; ng; mu] * o;
end
nu = 2; nh = 0; stop = 0;
info(5:7) = [0 1 0]; kit = 0;
% Iterate
while ~stop
if ng <= tolg, stop = 1;
else
[h mu] = geth(H,g,mu);
nh = norm(h); nx = tolx + norm(x);
if nh <= tolx*nx, stop = 2; end
end
if ~stop
xnew = x + h; h = xnew - x;
[stop fn gn Hn] = checkfgH(fun,xnew,varargin{:});
info(6) = info(6) + 1;
if ~stop
df = f - fn; accept = 0;
if df > 0
accept = 1; dL = (h'*(mu*h - g))/2;
elseif fn <= f + abs(f)*(1 + 100*eps) % Try gradient
df = (g + gn)'*(g - gn);
if df > 0, accept = 2; end
end
if accept % Update x and modify mu
kit = kit + 1;
if accept == 1 & dL > 0
mu = mu * max(1/3, 1 - (2*df/dL - 1)^3); nu = 2;
else
mu = mu*nu; nu = 2*nu;
end
x = xnew; f = fn; g = gn; H = Hn; ng = norm(g,inf);
if Trace
jtr = kit + 1;
X(:,jtr) = xnew; perf(:,jtr) = [f ng mu]'; end
else % Same x, increase mu
mu = mu*nu; nu = 2*nu;
end
if info(6) == maxeval, stop = 3; end
end
end
end
% Set return values
if Trace
jj = 1:jtr; X = X(:,jj);
perf = struct('f',perf(1,jj), 'ng',perf(2,jj), 'mu',perf(3,jj));
else, X = x; end
if stop < 0, f = NaN; ng = NaN; end
info = [f ng nh mu / norm(H,inf) kit info(6) stop];
\ No newline at end of file
src/dogleg.m 0 → 100644
+ 164
− 0
View file @ bd4f1101
function [X, info, perf] = dogleg(fun, x0, opts, varargin)
%DOGLEG Powell's dog-leg method for least squares.
% Find xm = argmin{f(x)} , where x is an n-vector and
% f(x) = 0.5 * sum(r_i(x)^2) .
% The functions r_i(x) (i=1,...,m) and the Jacobian matrix J(x)
% (with elements J(i,j) = Dr_i/Dx_j ) must be given by a MATLAB
% function with declaration
% function [r, J] = fun(x, p1,p2,...)
% The gradient of f is g = J' * r.
% p1,p2,... are parameters of the function. In connection with nonlinear
% data fitting they may be arrays with coordinates of the data points.
%
% Call [X, info] = dogleg(fun, x0)
% [X, info] = dogleg(fun, x0, opts, p1,p2,...)
% [X, info, perf] = dogleg(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for xm .
% opts : Either a struct with fields 'Delta', 'tolg', 'tolx', 'tolr' and
% 'maxeval', or a vector with the values of these options,
% opts = [Delta tolg tolx tolr maxeval].
% Delta initial trust region radius.
% The other options are used in stopping criteria:
% ||g(x)||_inf <= tolg or
% ||dx||_2 <= tolx*(tolx + ||x||_2) or
% ||r(x)||_inf <= tolr or
% no. of function evaluations exceeds maxeval .
% Default Delta = [0.1(1+||x0||), tolg = 1e-4, tolx = 1e-8,
% tolr = 1e-6, maxeval = 100.
% If the input opts has less than 5 elements, it is augmented by
% the default values. Also, zeros and negative elements are
% replaced by the default values.
% p1,p2,.. are passed dirctly to the function FUN .
%
% Output parameters
% X : If perf is present, then array, holding the iterates xk
% columnwise. Otherwise, computed solution vector.
% info : Performance information, vector with 7 elements:
% info(1:4) = final values of [f(x) ||g||_inf ||dx||_2 Delta]
% info(5:6) = no. of iteration steps and function evaluations
% info(7) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : Stopped by small r-vector
% 4 : No. of iteration steps exceeded
% -1 : x is not a real valued vector
% -2 : r is not a real valued column vector
% -3 : J is not a real valued matrix
% -4 : Dimension mismatch in x, r, J
% -5 : Overflow during computation
% perf : Struct with fields
% f : values of f(xk) ,
% ng : values of || g(xk) ||_inf ,
% Delta : values of the trust region radius.
% Version 10.10.08. hbn(a)imm.dtu.dk
% Check parameters and function call
if nargin < 2, stop = -1;
else
[stop x n] = checkx(x0);
if ~stop, [stop f r J] = checkrJ(fun,x0,varargin{:}); k = 1; end
end
if ~stop
g = -(J'*r); ngi = norm(g,inf); % negative gradient and its norm
if isinf(ngi), stop = -5; end
else
f = NaN; ngi = NaN;
end
if stop
X = x0; perf = []; info = [f ngi 0 opts(1) 0 1 stop];
return
end
% Finish initialization
if nargin < 3, opts = []; end
opts = checkopts('dogleg', opts); % use default options where required
if opts(1) == 0, opts(1) = 0.1 * (1 + norm(x)); end
Delta = opts(1); tolg = opts(2); tolx = opts(3); tolr = opts(4);
maxeval = opts(5);
Trace = nargout > 2;
if Trace
o = ones(1, maxeval);
X = x * o; perf = [f; ngi; Delta] * o;
end
nu = 2; nx = norm(x); fact = 1; stop = 0;
reduce = 0; nstep = 0; kit = 0;
% Iterate
while ~stop
if isinf(ngi), stop = -5;
elseif ngi <= tolg, stop = 1;
elseif Delta <= tolx*(tolx + nx), stop = 2;
elseif norm(r,inf) <= tolr, stop = 3;
elseif k >= maxeval, stop = 4;
else
if fact % Factorize and compute hGN
[Q R] = qr(J,0); [U S V] = svd(R);
s = diag(S); i = find(s > 100*eps*s(1));
Qr = -(r'*Q)'; UQr = U(:,i)'*Qr;
hGN = V(:,i)*(UQr ./ s(i));
nhGN = norm(hGN); fact = 0;
end
if nhGN > Delta % Include gradient
nh = Delta;
ng = norm(g); alpha = (ng / norm(J*g))^2;
gn = alpha*g; ngn = alpha*ng;
if ngn >= Delta
h = (Delta/ngn) * gn;
dLpre = Delta*(2*ngn - Delta)/(2*alpha);
else % Dog leg
b = hGN - gn; bb = b'*b; gb = gn'*b;
c = (Delta + ngn)*(Delta - ngn);
if gb > 0
beta = c / (gb + sqrt(gb^2 + c * bb));
else
beta = (sqrt(gb^2 + c * bb) - gb)/bb;
end
h = gn + beta*b;
dLpre = .5*alpha*(1 - beta)^2*ng^2 + beta*(2-beta)*f;
end
else
h = hGN; nh = nhGN;
dLpre = f;
if nh <= tolx*(tolx + nx), stop = 2; end
end
end
if ~stop
xnew = x + h; h = xnew - x; nstep = norm(h);
dL = f - .5*norm(r + J*h)^2;
[stop fn rn Jn] = checkrJ(fun,xnew,varargin{:}); k = k + 1;
if ~stop
df = f - fn;
if df > 0 & dL > 0 % Update x
kit = kit + 1; updx = 1;
x = xnew; f = fn; J = Jn; r = rn; fact = 1;
g = -(J'*r); ngi = norm(g,inf);
rho = dL / df;
else, rho = -1; updx = 0; end
% Update Delta
if abs(rho-1) < 0.2 & nh > Delta/3 & reduce <= 0
Delta = 3*Delta; nu = 2; reduce = 0;
elseif rho < 0.25
Delta = Delta/nu;
nu = 2*nu; reduce = 2;
else
reduce = reduce - 1;
end
if Trace && updx % store iterate
X(:,kit+1) = xnew; perf(:,kit+1) = [fn ngi Delta]';
end
if k > maxeval, stop = 3; end
end
end
end
% Set return values
if Trace
ii = 1 : kit+1; X = X(:,ii);
perf = struct('f',perf(1,ii), 'ng',perf(2,ii), 'Delta',perf(3,ii));
else, X = x; end
if stop < 0, f = NaN; ngi = NaN; end
info = [f ngi nstep Delta kit k stop];
\ No newline at end of file
src/efit1.dat 0 → 100644
+ 48
− 0
View file @ bd4f1101
% Exponential fit Data
% Hans Bruun Nielsen, IMM, DTU. 99.03.10
0.020 0.090542
0.040 0.124569
0.060 0.179367
0.080 0.195654
0.100 0.269707
0.120 0.286027
0.140 0.289892
0.160 0.317475
0.180 0.308191
0.200 0.336995
0.220 0.348371
0.240 0.321337
0.260 0.299423
0.280 0.338972
0.300 0.304763
0.320 0.288903
0.340 0.300820
0.360 0.303974
0.380 0.283987
0.400 0.262078
0.420 0.281593
0.440 0.267531
0.460 0.218926
0.480 0.225572
0.500 0.200594
0.520 0.197375
0.540 0.182440
0.560 0.183892
0.580 0.152285
0.600 0.174028
0.620 0.150874
0.640 0.126220
0.660 0.126266
0.680 0.106384
0.700 0.118923
0.720 0.091868
0.740 0.128926
0.760 0.119273
0.780 0.115997
0.800 0.105831
0.820 0.075261
0.840 0.068387
0.860 0.090823
0.880 0.085205
0.900 0.067203
src/efit2.dat 0 → 100644
+ 40
− 0
View file @ bd4f1101
0.02 0.2919
0.04 0.3408
0.06 0.3920
0.08 0.4214
0.10 0.4259
0.12 0.4516
0.14 0.5250
0.16 0.4464
0.18 0.4451
0.20 0.4329
0.22 0.4209
0.24 0.4093
0.26 0.3855
0.28 0.3824
0.30 0.3651
0.32 0.3709
0.34 0.3603
0.36 0.3372
0.38 0.3334
0.40 0.3102
0.42 0.3113
0.44 0.1700
0.46 0.2914
0.48 0.2797
0.50 0.2863
0.52 0.2570
0.54 0.2561
0.56 0.2553
0.58 0.2570
0.60 0.2585
0.62 0.2440
0.64 0.2497
0.66 0.2357
0.68 0.2297
0.70 0.2391
0.72 0.2288
0.74 0.1750
0.76 0.2309
0.78 0.2182
0.80 0.2056
src/huberobj.m 0 → 100644
+ 75
− 0
View file @ bd4f1101
function [f, S, r, J, g] = huberobj(fun, x, gamma, Htype, varargin)
%HUBEROBJ Value and gradient of the Huber objective function
% for the vector function given by
% function [r, J] = fun(x,p1,p2,...)
% p1,p2,... are parameters of the function and J is the Jacobean.
% If nargout < 4, then the function only needs to return r, and the
% gradient is not computed.
%
% Call [f, S, r] = huberobj(fun, x, gamma)
% [f, S, r] = huberobj(fun, x, gamma, Htype, p1,p2,...)
% [f, S, r, J, g] = huberobj(.....)
%
% Input parameters
% fun : Handle to the function.
% x : n-vector, argument of fun
% gamma : Huber threshold.
% Htype : Choice of Huber function,
% 1 : one-sided, r_i > 0 are neglected,
% 2 : one-sided, all r_i <= gamma are active,
% otherwise, all abs(r_i) <= gamma are active (default).
% p1,p2,.. are passed directly to the function FUN .
%
% Output parameters
% f : Huber objective function.
% S : Struct with the Huber active set. Fields
% s : Huber sign vector,
% A : active set (indices of small elements in r),
% N : inactive set (indices of large elements in r),
% L : vector [length(S.A), length(S.N)].
% r,J : Output from the evaluation of fun .
% If nargout < 4 , then fun only needs to return r(x) ,
% and the gradient is not computed.
% g : Gradient of f .
% Version 10.10.05. hbn(a)imm.dtu.dk
% Evaluate the function and set rounding error threshold
if nargout < 4
r = feval(fun,x, varargin{:});
tau = 10*eps*norm(r,inf);
else
[r J] = feval(fun,x, varargin{:});
tau = 10*eps*max(norm(r,inf), norm(J,inf)*norm(x,inf));
end
if gamma < tau
error(sprintf('gamma must be at least %9.2e',tau))
end
thr = gamma + tau;
% Huber sign vector
s = sign(r);
if nargin < 4 | abs(Htype - 1.5) ~= 0.5, Htype = 3; end
if Htype == 1 % Neglect positive contributions
w = find(-thr <= r & r <= tau); s(w) = 0;
p = find(s > 0); s(p) = 0;
elseif Htype == 2 % Negative contributions are active
w = find(r <= thr); s(w) = 0;
else % simple Huber function
w = find(abs(r) <= thr); s(w) = 0;
end
% Active set
N = find(s);
S = struct('s',s', 'A',w', 'N',N', 'L',[length(w) length(N)]);
% Compute function
rA = r(S.A); rN = r(S.N);
f = norm(rA)^2/(2*gamma) + norm(rN,1) - .5*gamma*S.L(2);
if nargout > 3 % gradient with check of rounding errors
g = J(S.A,:)'*rA/gamma + J(S.N,:)'*s(S.N);
thr = 10*eps*(norm(J(S.A,:),1)*norm(rA,inf)/gamma + norm(J(S.N,:),1));
if norm(g,inf) <= thr, g = zeros(size(g)); end
end % gradient
\ No newline at end of file
src/immoptset.m 0 → 100644
+ 85
− 0
View file @ bd4f1101
function opts = immoptset(mlfun, varargin)
%IMMOPTSET Set options for the immoptibox function MLFUN
% Call opts = immoptset(mlfun, p1,v1, p2,v2, ...)
%
% Input parameters
% mlfun : String with the name of the function or a handle to it.
% p1,p2,... : Strings with option names.
% v1,v2,... : Real numbers with the values of the options.
%
% Output parameter
% opts : Struct with the option names and their values. Options
% that do not appear as input are assigned their default value.
% Version 10.10.19. cv(a)imm.dtu.dk and hbn(a)imm.dtu.dk
% Define default options of respective functions:
dampnewton = struct('tau',1e-3, 'tolg',1e-4, 'tolx',1e-8, 'maxeval',100);
linesearch = struct('choice',1, 'cp1',1e-3, 'cp2',0.99, 'maxeval',10, ...
'amax',10);
ucminf = struct('Delta',1, 'tolg',1e-4, 'tolx',1e-8, 'maxeval',100);
dogleg = struct('Delta',0, 'tolg',1e-4, 'tolx',1e-8, 'tolr',1e-6, ...
'maxeval',100); % NB Default Delta depends on x
marquardt = struct('tau',1e-3, 'tolg',1e-4, 'tolx',1e-8, 'maxeval',100);
smarquardt = struct('tau',1e-3, 'tolg',1e-4, 'tolx',1e-8, ...
'maxeval',0, 'relstep',1e-7); % NB Default maxeval depends on n
nonlinhuber = struct('gamma',0, 'Htype',0, 'tau',1e-3, 'tolg',1e-4, ...
'tolx',1e-8, 'maxeval',100);
mexpfit = struct('const',0, 'tau',1e-3, 'tolg',1e-6, 'tolx',1e-10, ...
'maxeval',100);
nonlinsys = struct('Delta',0, 'tolg',1e-6, 'tolx',1e-8, ...
'maxeval',100, 'relstep',1e-6);
% List of functions:
funlist = struct('dampnewton',dampnewton, 'linesearch',linesearch, ...
'ucminf',ucminf, 'dogleg',dogleg, 'marquardt',marquardt, ...
'smarquardt',smarquardt, 'nonlinhuber',nonlinhuber, ...
'mexpfit',mexpfit, 'nonlinsys',nonlinsys);
funnames = fieldnames(funlist);
% Check mlfun
if isa(mlfun,'function_handle'), mlfun = lower(func2str(mlfun));
elseif ischar(mlfun), mlfun = lower(mlfun);
else
error('mlfun must be a string or a function handle')
end
foundname = 0; kfun = 0;
while ~foundname && (kfun < length(funnames))
kfun = kfun + 1;
foundname = strcmp(funnames(kfun), mlfun);
end
if ~foundname
error(['IMMOPTSET is not prepared for ', mlfun])
end
% Check existence of options and consistency of varargin
opts = funlist.(mlfun); optnames = fieldnames(opts);
lst = length(varargin);
for i = 1 : 2 : lst-1
if ~ischar(varargin{i})
error(['The value ' num2str(varargin{i}) ...
' is not assigned to an option'])
end
foundname = 0; j = 0;
while ~foundname && j < length(optnames)
j = j + 1; foundname = strcmpi(optnames{j},varargin{i});
end
if ~foundname
error(['Option ' num2str(varargin{i}) ...
' does not exist in function ',mlfun])
end
if ischar(varargin{i+1})
error(['No value specified for option ',num2str(varargin{i})])
end
% Update option
opts.(optnames{j}) = varargin{i+1};
end
if mod(lst,2) % odd number of elements in varargin
if ischar(varargin{lst})
error('No value assigned to the last option')
else
error(['The last value is not assigned to an option'])
end
end
\ No newline at end of file
src/linesearch.m 0 → 100644
+ 199
− 0
View file @ bd4f1101
function [xn,fn,gn,info,perf] = ...
linesearch(fun, x,f,g, h, opts, varargin)
%LINESEARCH Find am = argmin_{a > 0}{ P(a) = f(x+a*h) } , where x and
% h are given n-vectors and the scalar function f and its gradient g
% (with elements g(i) = Df/Dx_i ) must be given by a MATLAB function with
% declaration
% function [f, g] = fun(x,p1,p2,...)
% p1,p2,... are parameters of the function.
%
% Call [xn,fn,gn,info] = linesearch(fun,x,f,g,h)
% [xn,fn,gn,info] = linesearch(fun,x,f,g,h,opts,p1,p2,...)
% [xn,fn,gn,info,perf] = linesearch(......)
%
% Input parameters
% fun : Handle to the function.
% x : Current x.
% f,g : f(x) and g(x).
% h : Step vector.
% opts : Either a struct with fields 'choice', 'cp1', 'cp2', 'maxeval',
% 'amax' or a vector with the values of these options,
% opts = [choice cp1 cp2 maxeval amax].
% choice = 0 : exact line search.
% Otherwise: soft line search (Default).
% cp1, cp2 : options for stopping criteria.
% choice = 0: |P'(a)| <= cp1*|P'(0)| or c-b <= cp2*c ,
% where [b,c] is the current interval for a.
% Default cp1 = cp2 = 1e-3 .
% Otherwise: P(a) <= P(0) + a*cp1*P'(0) and
% P'(a) >= cp2*P'(0).
% Default cp1 = 1e-3, cp2 = 0.99 .
% maxeval : Maximum number of function evaluations.
% Default: maxeval = 10.
% amax : Maximal allowable step. Default amax = 10.
% p1,p2,.. are passed dirctly to the function FUN .
%
% Output parameters
% xn : x + am*h
% fn,gn : f(xn) and g(xn).
% info : Performance information, vector with 3 elements
% info(1) > 0 : am. Successfull call
% = 0 : h is not downhill or it is so large and maxeval
% so small, that a better point was not found.
% = -1 : x is not a real valued vector.
% = -2 : f is not a real valued scalar.
% = -3 : g or h is not a real valued vector.
% = -4 : g or h has different length from x.
% info(2) = slope ratio at the solution, P'(am)/P'(0) .
% info(3) = number of function evaluations used.
% perf : Struct with fields
% alpha : values of a ,
% phi : values of P(a) = f(x+a*h) ,
% slope : values of P'(a) .
% Version 10.09.28. hbn(a)imm.dtu.dk
% Initial check
if nargin < 5, error('Too few input parameters'), end
% Check OPTS
if nargin < 6 | isempty(opts), opts = 1; end
opts = checkopts('linesearch', opts); % use default options where required
% if opts(1) == 0, opts = checkopts(opts, [0 1e-3 1e-3 10 10]);
% else, opts = checkopts(opts, [1 1e-3 0.99 10 10]); end
choice = opts(1) ~= 0;
cp1 = opts(2); cp2 = opts(3); maxeval = opts(4); amax = opts(5);
% Default return values and simple checks
xn = x; fn = f; gn = g; info = [0 1 0];
[info(1) n] = check(x,f,g,h);
if info(1), return, else, stop = 0; end
x = x(:); h = h(:); % both are treated as column vectors
% Check descent condition
f0 = fn; df0 = dot(h,gn);
if df0 >= -10*eps*norm(h)*norm(gn) % not significantly downhill
info(1) = 0; return
end
Trace = nargout > 4;
if Trace
o = ones(1, maxeval);
X = x * o; perf = [0; f0; df0] * o;
end
% Finish initialization
if choice % soft line search
slope0 = cp1*df0; slopethr = cp2*df0;
else % exact line search
slope0 = 0; slopethr = cp1*abs(df0);
end
% Get an initial interval for am
a = 0; fa = fn; dfa = df0; stop = 0;
b = min(1, amax);
while ~stop
[stop fb g] = checkfgH(fun,x+b*h,varargin{:}); info(3) = info(3)+1;
if stop, info(1) = stop;
else
dfb = dot(g,h);
if Trace, perf(:,info(3)) = [b; fb; dfb]; end
if fb < f0 + slope0*b % new lower bound
info(1:2) = [b dfb/df0];
if choice, a = b; fa = fb; dfa = dfb; end
xn = x + b*h; fn = fb; gn = g;
if (dfb < min(slopethr,0)) && (info(3) < maxeval) && (b < amax)
% Augment right hand end
if ~choice, a = b; fa = fb; dfa = dfb; end
if 2.5*b >= amax, b = amax; else, b = 2*b; end
else, stop = 1; end
else, stop = 1; end
end
end % phase 1: expand interval
if stop >= 0 % OK so far. Check stopping criteria
stop = (info(3) >= maxeval) | (b >= amax & dfb < slopethr)... % Cannot improve
| (choice & (a > 0 & dfb >= slopethr)); % OK
end
if stop
if Trace
ii = 1 : info(3);
perf = struct('alpha',perf(1,1:k), 'phi',perf(2,1:k), 'slope',perf(3,1:k));
end
return
end
% Refine interval. Use auxiliary array xfd
xfd = [a b b; fa fb fb; dfa dfb dfb];
while ~stop
c = interpolate(xfd,n);
[stop fc g] = checkfgH(fun,x+c*h,varargin{:}); info(3) = info(3)+1;
if stop, info(1) = stop;
else
xfd(:,3) = [c; fc; dot(g,h)];
if Trace, perf(:,info(3)) = xfd(:,3); end
if choice % soft line search
if fc < f0 + slope0*c % new lower bound
info(1:2) = [c xfd(3,3)/df0];
xn = x + c*h; fn = fc; gn = g;
xfd(:,1) = xfd(:,3);
stop = xfd(3,3) > slopethr;
else % new upper bound
xfd(:,2) = xfd(:,3);
end
else % exact line search
if fc < fn % better approximant
info(1:2) = [c xfd(3,3)/df0];
xn = x + c*h; fn = fc; gn = g;
end
if xfd(3,3) < 0, xfd(:,1) = xfd(:,3); % new lower bound
else, xfd(:,2) = xfd(:,3); end % new upper bound
stop = abs(xfd(3,3)) <= slopethr...
| diff(xfd(1,1:2)) < cp2*xfd(1,2);
end
end
stop = stop | info(3) >= maxeval;
end % refine
% Return values
if Trace
ii = 1 : info(3);
perf = struct('alpha',perf(1,ii), 'phi',perf(2,ii), 'slope',perf(3,ii));
end
%============ Auxiliary functions ========================
function t = interpolate(xfd,n);
% Minimizer of parabola given by xfd = [a b; f(a) f(b); f'(a) dummy]
a = xfd(1,1); b = xfd(1,2); d = b - a; df = xfd(3,1);
C = diff(xfd(2,1:2)) - d*df;
if C >= 5*n*eps*b % Minimizer exists
A = a - .5*df*(d^2/C); d = 0.1*d;
t = min(max(a+d, A), b-d); % Ensure significant resuction
else
t = (a+b)/2;
end
function [err, n] = check(x,f,g,h)
% Check x
err = 0; sx = size(x); n = max(sx);
if (min(sx) ~= 1) | ~isreal(x) | any(isnan(x(:))) | isinf(norm(x(:)))
err = -1;
else
% Check f
sf = size(f);
if any(sf ~= 1) | ~isreal(f) | any(isnan(f(:))) | any(isinf(f(:)))
err = -2;
else
err = checkvec(g, n);
if ~err, err = checkvec(h, n); end
end
end
function err = checkvec(v,n)
sv = size(v);
if (min(sv) ~= 1) | ~isreal(v) | any(isnan(v(:))) | isinf(norm(v(:)))
err = -3;
elseif max(sv) ~= n, err = -4;
else, err = 0; end
\ No newline at end of file
src/linhuber.m 0 → 100644
+ 362
− 0
View file @ bd4f1101
function [X, info, perf] = linhuber(A,b, c, L, par, init)
%LINHUBER Minimizer of the piecewise quadratic
% f(x) = sum(phi(ri(x))) + c'*x + 0.5mu*||L*x||2^2 ,
% where ri(x) is the i'th element in r(x) = b - A*x and phi is the
% (possibly one-sided) Huber function with threshold gamma.
%
% Call: [X, info] = linhuber(A,b, c, L, par)
% [X, info] = linhuber(A,b, c, L, par, init)
% [X, info, perf] = linhuber(......)
%
% Input parameters
% A,b : m*n matrix and m-vector, respectively.
% c : n-vector or empty, in which case the term c'*x is omitted.
% L : Matrix with n columns or empty. In the case mu = par(2) > 0
% an empty L is treated as L = I, the identity matrix.
% par : Vector with one, two or three elements.
% par(1) = gamma, Huber threshold.
% par(2) = mu. Default: mu = 0.
% par(3) : Choice of Huber function,
% 1 : one-sided, rho(r) = 0 for r > 0,
% 2 : one-sided, all ri <= gamma are active
% Otherwise: standard Huber (default), ie
% equations with abs(ri)<=gamma are active.
% init: Choice of starting point,
% init is a vector: given x0,
% init is a struct like info below: given active set and
% factorization,
% If nargin < 6, x0 is the least squares solution to Ax "=" b.
%
% Output parameters
% X : If perf is present, then array, holding the iterates
% columnwise. Otherwise, computed minimizer of f.
% info : Struct with information on the performance and computed solution.
% Fields
% pf : Vector [f(x) ||g(x)||inf #iterations #factorizations] ,
% where x is the computed minimizer.
% pf(1) = -Inf indicates an unbounded problem.
% S : Struct with the Huber active set at the solution.
% Fields
% s : Huber sign vector,
% A : active set (indices of small residuals),
% N : inactive set (indices of large residuals),
% L : vector [length(S.A), length(S.N)].
% R : n*n array with Cholesky factor of active set.
% p : Permutation vector used in the factorization.
% perf : Iteration history. Struct with fields
% f : values of f(x) ,
% ng : values of || g(x) ||inf
% nA : number of active equations.
%
% Method
% Essentially (without downdating of the factorization) as described
% in H.B. Nielsen, "Computing a minimizer of a piecewise quadratic -
% Implementation", IMM-REP-1998-14, IMM, DTU.
% Version 10.10.05. hbn(a)imm.dtu.dk
% Initialize
X = []; S = []; R = []; p = [];
pf = [-Inf NaN 0 0]; perf = [];
[aux Htype] = lhaux(A,b,c,L,par);
[m n] = size(A); kmax = 2*m; fact = 0;
kcomp = 1;
if nargin > 5 % given info on starting point
if isstruct(init) % given active set and factorization
kcomp = 2;
S = init.S;
[f g] = lhobj(zeros(n,1),b,S,c,L,A,aux);
if aux(3) == 0 % no L-contribution. Reuse factorization
[x u cf R p fact] = lhfacsolv(S,A,L,g,aux,S,init.R,init.p);
else
[x u cf R p fact] = lhfacsolv(S,A,L,g,aux);
end
else, x = init(:); end % given x0
else % Start with least sqares solution to Ax "=" b
x = pinv(A)*b; fact = 1;
end
Trace = nargout > 2;
if Trace
X = zeros(n,kmax); perf = zeros(3,kmax);
end
% Iterate
k = 0; stop = 0;
while ~stop & k < kmax
Sp = S;
k = k + 1;
[r S tau] = lhsign(A,b,x,aux,Htype);
[f g] = lhobj(x,r,S,c,L,A,aux); ng = norm(g,inf);
if Trace, X(:,k) = x; perf(:,k) = [f; ng; S.L(1)]; end
if ng == 0 | ( k > kcomp & isequal(Sp,S) )
stop = 1; t = 0;
else
[h u lc R p rf] = lhfacsolv(S,A,L,g,aux,Sp,R,p); fact = fact + rf;
t = lhlines(S,r,u,lc,aux,Htype);
if t > 0, x = x + t*h;
else, stop = 1; end
end
end
if t*ng == 0, pf(1) = f; end
pf(2:4) = [ng k-1 fact];
info = struct('pf',pf, 'S',S, 'R',R, 'p',p);
if Trace
kk = 1 : k; X = X(:,kk);
perf = struct('f',perf(1,kk), 'ng',perf(2,kk), 'nA',perf(3,kk));
else, X = x; end
%%%%%%%%%%%%%%%%%%%% Auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [aux, Htype] = lhaux(A,b,c,L,par)
% Simple checks and set thresholds. 04.08.24
aux = [];
e = [A(:); b(:); c(:); L(:); par(:)];
v = isreal(e) & all(~isinf(e)) & all(~isnan(e));
if ~v
error('Indata contains non-real elements')
end
[m n] = size(A); sb = size(b);
v = n > 0 & (min(sb) == 1 & max(sb) == m);
if ~v
disp('A must be a matrix and ')
error(sprintf('b must be a vector of length %g',m))
end
sc = size(c); sL = size(L); tn = sprintf(' %g',n);
v = max(sc) == 0 | (min(sc) == 1 & max(sc) == n);
if ~v
error(['c must be empty or a vector of length' tn])
end
v = max(sL) == 0 | sL(2) == n;
if ~v
error(['L must be empty or a matrix with' tn ' columns'])
end
sp = size(par);
if min(sp) ~= 1
error('par must be a scalar or a vector')
end
if max(sp) < 2, par = [par 0]; end
if ~all(par(1:2) >= 0)
error('par(1:2) = [gamma mu] cannot be negative')
end
gamma = par(1); mu = par(2); gm = gamma*mu;
% Rounding error threshold
tau = 100*eps*norm(b,inf);
if par(1) < tau
error(sprintf('gamma must be at least %9.2e',tau))
end
% Diverse norms and thresholds
nAi = norm(A,inf); nA1 = norm(A,1);
sigma0 = 10*sqrt(eps * nA1*nAi); sigma = sigma0;
if gm > 0
sgm = sqrt(gm);
if isempty(L) % L = I
nL1 = 1; nLi = 1;
if sgm >= sigma0, sigma = sgm; end
else
nL1 = norm(L,1); nLi = norm(L,inf);
if gm*nL1*nLi >= sigma0^2
nBB = (nA1 + sgm*nL1) * max(nAi, sgm*nLi);
sigma = 10*sqrt(eps * nBB);
sigma0 = sigma;
end
end
else
nL1 = 0; nLi = 0;
end % mu-contribution
% aux = [gamma tau gamma*mu sigma sigma0 ...
% ||A||inf gamma*(||A||1 + mu*||L||1*||L||inf) gamma*||c||inf]
aux = [gamma tau gm sigma sigma0 ...
nAi gamma*nA1 + gm*nL1*nLi gamma*norm(c,inf)];
if length(par) < 3 | abs(par(3) - 1.5) ~= 0.5, Htype = 3;
else, Htype = par(3); end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [h,u,lc,R,p,refac] = lhfacsolv(S,A,L,g,aux,Sp,R,p)
% Factorize Hessian and find Newton step. 04.08.25
n = length(g);
if nargin < 6 | isempty(Sp), refac = 1;
else % possible update from previos factorization
w = zeros(size(A,1),1); w(S.A) = 1; w(Sp.A) = w(Sp.A) - 1;
if any(w == -1), refac = 1; % downdates
else
refac = 0; E = find(w); % augmented active set
if ~isempty(E) % update factorization
[Q R p1] = qr([R; A(E,p)],0); p = p(p1);
end
end
end
if refac % refactorize
d = aux(4)*ones(size(g));
if aux(3) > 0 & ~isempty(L) % mu-contribution with L ~= I
C = [A(S.A,:); sqrt(aux(3))*L; diag(d)];
else
C = [A(S.A,:); diag(d)];
end
[Q R p] = qr(C,0);
end
% Check for almost rank deficient
rd = any(abs(diag(R)) < 10*aux(5));
% Get solution to Hx = -g
g = -g(:); x(p) = R \ (g(p)' / R)'; x = x(:);
if rd % Check consistency
v(p) = R \ (aux(4)^2*x(p)' / R)'; v = v(:); h = x + v;
if S.L(1) > 0, gR = A(S.A,:)' * (A(S.A,:)*h);
else, gR = zeros(n,1); end
if aux(3) > 0 % mu-contribution
if isempty(L), gR = gR + aux(3)*h;
else, gR = gR + aux(3)*(L' * (L*h)); end
end
if norm(g - gR) > sqrt(eps)*norm(g) % inconsistent (?)
h = v / norm(v);
end
else % full rank. iterative refinement step
if S.L(1) > 0, Hx = A(S.A,:)' * (A(S.A,:)*x);
else, Hx = zeros(n,1); end
if aux(3) > 0 % mu-contribution
if isempty(L), Hx = Hx + aux(3)*x;
else, Hx = Hx + aux(3)*(L' * (L*x)); end
end
v(p) = R \ ((g(p) - Hx(p))' / R)'; h = x + v(:);
end
% u = A*h, modified for 'zero' elements
u = A*h;
z = find(abs(u) <= eps * aux(4) * norm(h,inf) );
u(z) = 0;
% Offset values for line search coefficients
if aux(3) > 0 % mu-contribution
if isempty(L), c2 = aux(3)*norm(h)^2;
else, c2 = aux(3)*norm(L*h)^2; end
else
c2 = 0;
end
lc = [dot(g,h) c2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [kv, lc] = lhkinks(S,r,u,lc,aux,Htype)
% Get kink values and finish initialization of line search coeffs. 04.08.23
gam = aux(1); tau = aux(2);
if nargin < 6 | abs(Htype - 1.5) ~= 0.5, Htype = 3; end
thr = [-1 0; -inf 1; -1 1]; thr = gam*thr(Htype,:); % Huber thresholds
kv = zeros(2*length(r),2); k = 0;
if ~isempty(S.A)
% Remove leaving 'maybes'
act = S.A; w = ones(size(act));
o = find((u(act) > 0 & abs(r(act) - thr(1)) <= tau) | ...
(u(act) < 0 & abs(r(act) - thr(2)) <= tau));
w(o) = 0; act = act(find(w));
if ~isempty(act) % remaining actives
lc(2) = lc(2) + norm(u(act))^2; % update initial line search coeffs.
n = find(u(act) < 0); p = find(u(act) > 0);
ii = act([p; n])'; li = length(ii);
if li % get leave values
d = [repmat(thr(1),length(p),1); repmat(thr(2),length(n),1)];
kv(k+[1:li],:) = [(r(ii) - d)./u(ii) -ii]; k = k + li;
end
end
end % actives
% Enter-leave kinks
i = find(u(S.N) .* r(S.N) > 0); E = S.N(i);
if Htype == 1 % add neglected components
E = [E(:); find(r > tau & u > 0)]';
end
if ~isempty(E)
i = E(find(u(E) > 0)); j = E(find(u(E) < 0));
ii = [i j]'; li = length(ii);
d = [repmat(thr([2 1]),length(i),1); repmat(thr,length(j),1)];
kv(k+[1:li],:) = [(r(ii) - d(:,1))./u(ii) ii]; k = k + li;
kv(k+[1:li],:) = [(r(ii) - d(:,2))./u(ii) -ii]; k = k + li;
end
% sort by kink values
[sk s] = sort(kv(1:k,1));
% if k > 2 & kv(s(k),2) < 0, s = s(1:k-1); end % remove last leave
kv = kv(s,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function t = lhlines(S,r,u,lc,aux,Htype)
% Piecewise linear line search. 04.08.23
[kv lc] = lhkinks(S,r,u,lc,aux,Htype);
c1 = lc(1); c2 = lc(2);
tj = 0; j = 0; nk = size(kv,1); fd = 0;
while ~fd & j < nk
j = j+1; t1 = tj; tj = kv(j,1);
dt = tj - t1; cj = c2*dt - c1;
if c2 & (cj >= 0) % passed minimum
t = t1 + c1/c2; fd = 1;
else % update coefficients
c1 = -cj; ii = kv(j,2);
c2 = c2 + sign(ii)*u(abs(ii))^2;
end
end
if fd, return
elseif nk > 1, t = mean(kv(nk-1:nk,1));
elseif nk == 1, t = kv(1,1) + aux(2)/abs(u(abs(kv(1,2))));
elseif c2, t = c1/c2;
else, t = -1; end % error return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f,g] = lhobj(x,r,S,c,L,A,aux)
% Value and gradient of objective function. 04.08.24
% Compute function
fA = r(S.A); fN = r(S.N); gamma = aux(1);
f = norm(fA)^2/(2*gamma) + norm(fN,1) - .5*gamma*S.L(2);
if ~isempty(c), f = f + dot(c,x); end
if aux(3) > 0
if isempty(L), Lx = x; else, Lx = L*x; end
f = f + aux(3)/(2*gamma)*norm(Lx)^2;
end
if nargout > 1 % gamma*gradient
fv = -gamma*S.s'; fv(S.A) = -fA; g = A'*fv;
if ~isempty(c), g = g + gamma*c; end
if aux(3) > 0
if isempty(L), g = g + aux(3)*x;
else, g = g + L'*(L*(aux(3)*x)); end
end
% Threshold for rounding errors
thr = 100*eps*(aux(7)*norm(fv,inf) + aux(7)*norm(x,inf) + aux(8));
if norm(g,inf) <= thr, g = zeros(size(g)); end
end % gradient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r,S,tau] = lhsign(A,b,x,aux,Htype)
% Residual and Huber sign vector etc
% Rounding error and effective thresholds
tau = max(aux(2), 10*eps*aux(6)*norm(x,inf));
thr = aux(1) + tau;
% Residual and sign vector
r = b - A*x; s = sign(r);
if nargin < 5 | abs(Htype - 1.5) ~= 0.5 % simple Huber function
w = find(abs(r) <= thr); s(w) = 0;
elseif Htype == 1 % Neglect positive contributions
w = find(-thr <= r & r <= tau); s(w) = 0;
p = find(s > 0); s(p) = 0;
else % Negative contributions are active
w = find(r <= thr); s(w) = 0;
end
N = find(s);
% Return result
S = struct('s',s', 'A',w', 'N',N', 'L',[length(w) length(N)]);
\ No newline at end of file
src/marquardt.m 0 → 100644
+ 126
− 0
View file @ bd4f1101
function [X, info, perf] = marquardt(fun, x0, opts, varargin)
%MARQUARDT Levenberg-Marquardt's method for least squares.
% Find xm = argmin{f(x)} , where x is an n-vector and
% f(x) = 0.5 * sum(r_i(x)^2) .
% The functions r_i(x) (i=1,...,m) and the Jacobian matrix J(x)
% (with elements J(i,j) = Dr_i/Dx_j ) must be given by a MATLAB
% function with declaration
% function [r, J] = fun(x,p1,p2,...)
% p1,p2,... are parameters of the function. In connection with
% nonlinear data fitting they may be arrays with coordinates of
% the data points.
%
% Call [X, info] = marquardt(fun, x0)
% [X, info] = marquardt(fun, x0, opts, p1,p2,...)
% [X, info, perf] = marquardt(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for xm .
% opts : Either a struct with fields 'tau', 'tolg', 'tolx' and 'maxeval',
% or a vector with the values of these options,
% opts = [tau tolg tolx maxeval].
% tau used in starting value for Marquardt parameter:
% mu = tau * max( [J(x0)'*J(x0)]_(ii) ) .
% The other options are used in stopping criteria:
% ||g||_inf <= tolg or
% ||dx||_2 <= tolx*(tolx + ||x||_2) or
% no. of function evaluations exceeds maxeval .
% Default tau = 1e-3, tolg = 1e-4, tolx = 1e-8, maxeval = 100.
% If the input opts has less than 4 elements, it is augmented by
% the default values. Also, zeros and negative elements are
% replaced by the default values.
% p1,p2,.. are passed directly to the function FUN .
%
% Output parameters
% X : If perf is present, then array, holding the iterates
% columnwise. Otherwise, computed solution vector.
% info : Performance information, vector with 7 elements:
% info(1:4) = final values of
% [f(x) ||g(x)||inf ||dx||2 mu/max( [J(x)'*J(x)]_(ii) )] .
% info(5:6) = no. of iteration steps and function evaluations.
% info(7) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : No. of function evaluations exceeded
% -1 : x is not a real valued vector
% -2 : r is not a real valued column vector
% -3 : J is not a real valued matrix
% -4 : Dimension mismatch in x, r, J
% -5 : Overflow during computation
% perf : Struct with fields
% f : values of f(xk) ,
% ng : values of || g(xk) ||_inf ,
% mu : values of Marquardt damping parameter.
% Version 10.11.03. hbn(a)imm.dtu.dk
% Check parameters and function call
if nargin < 2, stop = -1;
else
[stop x n] = checkx(x0);
if ~stop, [stop f r J] = checkrJ(fun,x0,varargin{:}); k = 1; end
end
if ~stop
g = J'*r; ng = norm(g,inf); A = J'*J;
if isinf(ng) | isinf(norm(A(:),inf)), stop = -5; end
else
f = NaN; ng = NaN;
end
if stop
X = x0; perf = []; info = [f ng 0 opts(1) 0 1 stop];
return
end
if nargin < 3, opts = []; end
opts = checkopts('marquardt', opts); % use default options where required
tau = opts(1); tolg = opts(2); tolx = opts(3); maxeval = opts(4);
% Finish initialization
mu = tau * max(diag(A));
Trace = nargout > 2;
if Trace
o = ones(1, maxeval);
X = x * o; perf = [f; ng; mu] * o;
end
nu = 2; nh = 0; stop = 0; kit = 0;
% Iterate
while ~stop
if ng <= tolg, stop = 1;
else
[h mu] = geth(A,g,mu);
nh = norm(h); nx = tolx + norm(x);
if nh <= tolx*nx, stop = 2; end
end
if ~stop
xnew = x + h; h = xnew - x; dL = (h'*(mu*h - g))/2;
[stop fn rn Jn] = checkrJ(fun, xnew, varargin{:}); k = k + 1;
if ~stop
if length(rn) ~= length(r)
df = f - fn;
else % more accurate
df = ( (r - rn)' * (r + rn) )/2;
end
if (dL > 0) & (df > 0) % Update x and modify mu
kit = kit + 1;
x = xnew; f = fn; J = Jn; r = rn;
A = J'*J; g = J'*r; ng = norm(g,inf);
mu = mu * max(1/3, 1 - (2*df/dL - 1)^3); nu = 2;
if Trace
X(:,kit+1) = xnew; perf(:,kit+1) = [fn ng mu]'; end
else % Same x, increase mu
mu = mu*nu; nu = 2*nu;
end
if k > maxeval, stop = 3;
elseif isinf(ng) | isinf(norm(A(:),inf)), stop = -5; end
end
end
end
% Set return values
if Trace
ii = 1 : kit+1; X = X(:,ii);
perf = struct('f',perf(1,ii), 'ng',perf(2,ii), 'mu',perf(3,ii));
else, X = x; end
if stop < 0, f = NaN; ng = NaN; end
info = [f ng nh mu/max(diag(A)) kit k stop];
\ No newline at end of file
src/mexpfit.m 0 → 100644
+ 204
− 0
View file @ bd4f1101
function [Z, c, info, perf] = mexpfit(tyw, z0, opts)
%MEXPFIT Least squares fit with a sum of decaying exponentials: Find
% zm = argmin{f(z)} , where z is a p-vector with positive elements,
% f(z) = ||W*r(z)||^2 , r(z) = y - F(z)*c(z) ,
% where F(i,:) = [exp(-z1*ti) ... exp(-zp*ti)] or
% F(i,:) = [exp(-z1*ti) ... exp(-zp*ti) 1] ,
% depending on the option 'const', and c(z) is the least squares
% solution to W*F*c "=" W*y .
%
% Call [Z, c, info] = mexpfit(tyw, z0)
% [Z, c, info] = mexpfit(tyw, z0, opts)
% [Z, c, info, perf] = mexpfit(...)
%
% Input parameters
% tyw : Data points and weights. Array with 2 or 3 columns,
% tyw(:,1:2) : Abscissas and ordinates of data points.
% tyw(:,3) : Weights. If tyw has less than 3 columns, then
% all weights are set to 1.
% z0 : Starting guess for z .
% opts : Either a struct with fields 'const', 'tau', 'tolg', 'tolx',
% and 'maxeval', or a vector with the values of these options,
% opts = [const tau tolg tolx maxeval].
% const : If positive then there is a constant term.
% tau : Used in starting value for Marquardt parameter for the
% optimization, mu = tau * max( [J(z0)'*J(z0)]_(ii) ) .
% The other options are used in stopping criteria:
% ||g||_inf <= tolg or
% ||dz||_2 <= tolx*(tolx + ||z||_2) or
% no. of function evaluations exceeds maxeval .
% Default const = 0, tau = 1e-3, tolg = 1e-6, tolx = 1e-10,
% maxeval = 100. If the input opts has less than 5 elements, it
% is augmented by the default values. Also, zeros and negative
% elements are replaced by the default values.
%
% Output parameters
% Z : If perf is present, then an array, holding the z-iterates
% columnwise. Otherwise, computed solution vector zm.
% c : If info.fail = 0, then c holds coefficients [c1, ..., cq], where
% opts.const > 0 : q = p+1; cq is the constant term
% otherwise : q = p .
% info : Struct with information on performance. Fields
% fail: Successful run if fail = 0. Otherwise info.msg tells
% what went wrong.
% msg: Text message about performance.
% vals: Vector with final values of f(z), ||g(z)||inf, ||dz||2
% its: Vector with number of iterations and function evaluations.
% perf : Struct with fields
% f : values of f(zk) ,
% ng : values of || g(zk) ||_inf ,
% mu : values of Marquardt damping parameter.
% Version 10.10.29. hbn(a)imm.dtu.dk
% Initialise with check of input
Z = []; c = []; perf = [];
if nargin < 3, opts = []; end
[z info opts nc t y W] = init(tyw, z0, opts);
if info.fail, return, end
p = length(z);
if p == 1 & nc == 1 % simple case
[Z c f] = simple(t, y, W);
if ischar(c), info.fail = 1; info.msg = c;
else, info.vals(1) = f; info.its = [0 1]; end
return
end
% General case. Remaining options and initial values
tau = opts(2); tolg = opts(3); tolx = opts(4); maxeval = opts(5);
[r J c] = func(z, t,y,W, nc); nv1 = 1;
A = J'*J; g = J'*r; f = (r'*r)/2; ng = norm(g,inf);
mu = tau * max(diag(A));
Trace = nargout > 3;
if Trace
Z = z*ones(1,maxeval); perf = [f; ng; mu]*ones(1,maxeval);
end
nu = 2; mok = 0; nh = 0; n = size(A,1); stop = 0;
while ~stop
if ng <= tolg, stop = 1;
else
hh = (A + mu*eye(n))\(-g);
nh = norm(hh); nz = tolx + norm(z);
if nh <= tolx*nz, stop = 2;
elseif nh >= nz/eps, stop = 4; end % Almost singular ?
end
if ~stop
% Control h to guarantee znew > 0
if p == 1
hneg = -0.75*z; hpos = z;
else
hneg = ([z(2:p); -z(p)] - z)/3;
hpos = ([2*z(1); z(1:p-1)] - z)/3;
end
hmod = max( min(hh, hpos), hneg);
znew = z + hmod; h = znew - z;
dr = J*h; dL = -(dr' * (2*r + dr)); % twice linear gain
[rn Jn c] = func(znew, t,y,W, nc); nv1 = nv1 + 1;
if ischar(c), info.fail = 1; info.msg = c; return, end
fn = (rn'*rn)/2; df = (r - rn)' * (r + rn); % twice actual gain
if (dL > 0) & (df > 0)
mok = mok + 1;
z = znew; f = fn; J = Jn; r = rn;
A = J'*J; g = J'*r; ng = norm(g,inf);
if norm(hh - hmod) % step was reduced. Increase mu
mu = mu*nu; nu = 2*nu;
else
mu = mu * max(1/3, 1 - (2*df/dL - 1)^3); nu = 2;
end
if Trace, Z(:,mok+1) = z; perf(:,mok+1) = [f; ng; mu]; end
else % Marquardt fail
mu = mu*nu; nu = 2*nu;
if mok > n & nu > 8, stop = 5; end
end
if nv1 >= maxeval, stop = 3; end
end
end
% Set return values
if Trace
jj = 1 : (mok+1); Z = Z(:,jj);
perf = struct('f',perf(1,jj), 'ng',perf(2,jj), 'mu',perf(3,jj));
else, Z = z; end
if stop == 1, info.msg = 'Iterations stopped by a small gradient';
elseif stop == 2, info.msg = 'Iterations stopped by a small z-step';
elseif stop == 3, info.msg = 'Iterations stopped by maxeval';
elseif stop == 4, info.msg = 'Iterations stopped by extreme step';
elseif stop == 5, info.msg = 'Iterations stopped by stalling'; end
info.vals = [f ng nh];
info.its = [mok nv1];
% ========== auxiliary functions =================================
function [z, info, opts, nc, t,y,W] = init(tyw, z0, opts)
Z = NaN; t = []; y = []; W = [];
info = struct('fail',0, 'msg','', 'vals',[NaN NaN NaN], 'its',[0 0 0]);
sz = size(z0);
if ~isreal(z0) | min(sz) ~= 1
info.fail = 1;
info.msg = 'z must be a vector with real elements'; return
end
if any(z0 <= 0)
info.fail = 1;
info.msg = 'z must must have strictly positive elements'; return
end
opts = checkopts('mexpfit', opts);
[m q] = size(tyw); p = length(z0);
if opts(1), nc = p+1; else, nc = p; end
if q > 2 % remove points with non-positive weight
i = find(tyw(:,3) > 0); ma = length(i);
if ma < m, tyw = tyw(i,:); m = ma; end
else, tyw = [tyw ones(m,1)]; end
if p + nc > m
info.fail = 1; info.msg = 'Underdetermined problem'; return
end
t = tyw(:,1); y = tyw(:,2); W = tyw(:,3)*ones(1,p);
z = sort(z0(:),1,'descend');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r, J, c] = func(z, t,y,W, nc)
% Function for exponential fit
p = length(z); m = length(t);
E = W .* exp(t * (-z')); y = W(:,1).*y;
if nc > p, F = [E W(:,1)]; else, F = E; end
[Q R] = qr(F, 0); rc = rcond(R); pp = size(R,2);
if rc < max(pp,10) * eps
r = []; J = [];
c = 'Rank deficient system'; return
elseif rc < sqrt(eps) % use iterative refinement
md = max(abs(diag(R))) * sqrt(eps);
RR = qr([R; md*eye(pp)]); RR = triu(RR(1:pp,:));
c = RR \ (y' * Q)';
r = y - F*c;
c = c + RR \ (r' * Q)';
else
md = 0;
c = R \ (y' * Q)';
end
r = y - F*c; % residual. Now get Jacobian - without iterative refinement
tE = (t * ones(1,p)) .* E; rho = (r' * tE)';
tEc = tE .* (ones(m,1) * c(1:p)');
if nc > p, B = [diag(rho); zeros(1,p)]; else, B = diag(rho); end
J = tEc - Q * (R' \ (F'*tEc - B));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [z, c, f] = simple(t, y, w)
% One exponential only. Use logaritmic transformation
% Ignore negative ordinates
m = length(y); a = find(y > 0); ma = length(a);
if ma < 2
c = 'Simple case. Less than 2 active points';
f = NaN; return
end
if ma < m, t = t(a); y = y(a); w = w(a); end
wt = w .* y; % transformed weights
x = [wt wt.*t] \ (wt .* log(y));
z = -x(2); c = exp(x(1));
r = w .* (y - c * exp(-z*t));
f = 0.5 * norm(r)^2;
\ No newline at end of file
src/nonlinhuber.m 0 → 100644
+ 116
− 0
View file @ bd4f1101
function [X, S, info, perf] = nonlinhuber(fun, x0, opts, varargin)
%NONLINHUER Levenberg-Marquardt type method for Huber estimation.
% Find xm = argmin{f_gamma(x)} , where x is an n-vector and
% f_gamma(x) = 0.5/gamma * sum(phi(r_i(x))) , the Huber objective function.
% The functions r_i(x) (i=1,...,m) and the Jacobian matrix J(x)
% (with elements J(i,j) = Df_i/Dx_j ) must be given by a MATLAB function
% with declaration
% function [r, J] = fun(x,p1,p2,...)
% p1,p2,... are parameters of the function.
%
% Call [X, S, info] = nonlinhuber(fun, x0)
% [X, S, info] = nonlinhuber(fun, x0, opts, p1,p2,...)
% [X, S, info, perf] = nonlinhuber(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for Huber estimator.
% opts : Either a struct with fields 'gamma', 'Htype', 'tau', 'tolg',
% 'tolx' and 'maxeval', or a vector with the values of these
% options, opts = [gamma Htype tau tolg tolx maxeval].
% gamma: Huber threshold.
% Htype: Choice of Huber function,
% Htype = 1: one-sided, r_i > 0 are neglected,
% 2: one-sided, all r_i <= gamma are active,
% otherwise, all abs(r_i) <= gamma are active.
% tau: used to get starting value for Marquardt parameter:
% mu = tau * max( [J(x0)'*J(x0)]_(ii) ) .
% The remaining options are used in stopping criteria:
% ||g||_inf <= tolg or
% ||dx||_2 <= tolx*(tolx + ||x||_2) or
% no. of function evaluations exceeds maxeval .
% Default gamma = ||r(x0)||inf, Htype = 0,
% tau = 1e-3, tolg = 1e-4, tolx = 1e-8, maxeval = 100.
% If the input opts has less than 6 elements, it is augmented by
% the default values. Also, zeros and negative elements are
% replaced by the default values.
% p1,p2,.. are passed directly to the function FUN .
%
% Output parameters
% X : If perf is present, then array, holding the iterates columnwise,
% otherwise computed solution.
% S : Struct with the Huber active set for the last col. in X. Fields
% s: Huber sign vector,
% A: active set (indices of small components in r),
% N: indices of large components,
% L: vector [length(S.A), length(S.N)].
% info: Performance information, vector with 7 elements:
% info(1:4) = final values of
% [f(x) ||g(x)||inf ||dx||2 mu/max( [J(x)'*J(x)]_(ii) )] ,
% info(5:6) = no. of iteration steps and function evaluations.
% info(7) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : No. of function evaluations exceeded.
% perf : Struct with fields
% f : values of f(xk) ,
% ng : values of || g(xk) ||_inf ,
% mu : values of Marquardt damping parameter.
% Version 10.11.16. hbn(a)imm.dtu.dk
% Check options
if nargin < 3, opts = []; end
opts = checkopts('nonlinhuber', opts); % use default options where required
gamma = opts(1); Htype = opts(2);
tau = opts(3); tolg = opts(4); tolx = opts(5); maxeval = opts(6);
% Initialize
x = x0(:);
if gamma == 0 % get gamma such that all eqs are active
gamma = norm(feval(fun,x, varargin{:}), inf);
neval = 1;
else, neval = 0; end
[f S r J g] = huberobj(fun,x,gamma,Htype,varargin{:}); neval = neval + 1;
A = J'*J; mu = tau * max(diag(A));
ng = norm(g,inf);
nit = 1; % 1 + no. of iterations
nu = 2; nh = 0; stop = 0;
Trace = nargout > 3;
if Trace
o = ones(1,maxeval);
X = x*o; perf = [f; ng; mu]*o;
end
info = zeros(size(x));
% Iterate
while ~stop
if ng <= tolg, stop = 1;
else
[h info] = linhuber(-J,r,[],[],[gamma mu Htype],info);
nh = norm(h); nx = tolx + norm(x);
if nh <= tolx*nx, stop = 2; end
end
if ~stop
xnew = x + h; h = xnew - x; dL = f - info.pf(1);
[fn Sn rn Jn gn] = huberobj(fun,xnew,gamma,Htype,varargin{:});
neval = neval + 1; df = f - fn;
if (dL > 0) & (df > 0) % Update x and modify mu
nit = nit + 1;
x = xnew; f = fn; J = Jn; r = rn; S = Sn; g = gn;
ng = norm(g,inf);
mu = mu * max(1/3, 1 - (2*df/dL - 1)^3); nu = 2;
if Trace
X(:,nit) = xnew; perf(:,nit) = [fn ng mu]'; end
else % Same x, increase mu
mu = mu*nu; nu = 2*nu;
end
if neval == maxeval, stop = 3; end
end
end
% Set return values
if Trace
ii = 1 : nit; X = X(:,ii);
perf = struct('f',perf(1,ii), 'ng',perf(2,ii), 'mu',perf(3,ii));
else, X = x; end
info = [f ng nh mu/max(diag(A)) nit-1 neval stop];
\ No newline at end of file
src/nonlinsys.m 0 → 100644
+ 212
− 0
View file @ bd4f1101
function [X, info, perf] = nonlinsys(fun, x0, opts, B0, varargin)
%NONLINSYS Powell's dog-leg method for non-linear system of equations.
% Find an n-vector xm, such that f_i(xm) = 0 , i=1,...,n.
% The vector function f(x) (with components f_i, i=1,...,n) must be
% given by a MATLAB function with declaration
% function f = fun(x, p1,p2,...)
% % % p1,p2,... are parameters of the function.
%
% Call
% [X, info] = nonlinsys(fun, x0)
% [X, info] = nonlinsys(fun, x0, opts)
% [X, info] = nonlinsys(fun, x0, opts, B0, p1,p2,...)
% [X, info, perf] = nonlinsys(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for xm .
% opts : Vector with five elements:
% opts(1) initial trust region radius Delta
% opts(2:4) used in stopping criteria:
% ||f||inf <= opts(2) or
% ||dx||2 <= opts(3)*(opts(3) + ||x||2) or
% no. of iteration steps exceeds opts(4) .
% opts(5) "relative" step length for difference approximations.
% Default opts = [0.1(1+||x0||) 1e-6 1e-8 100 1e-6]
% If the input opts has less than 5 elements, it is
% augmented by the default values.
% B0 : Initial approximation to the Jacobian of f.
% If B0 is not given, a forward difference approximation
% to it is used.
% p1,p2,.. are passed dirctly to the function FUN .
%
% Output parameters
% X : If perf is present, then array, holding the iterates
% columnwise. Otherwise, computed solution vector.
% info : Performance information, vector with 6 elements:
% info(1:3) = final values of
% [||f(x)||inf ||dx||2 Delta]
% info(4:5) = no. of iteration steps and function avaluations
% info(6) = 1 : Stopped by small f-vector
% 2 : Stopped by small x-step
% 3 : No. of iteration steps exceeded
% -1 : x is not a real valued vector
% -2 : f is not a real valued column vector
% -3 : Dimension mismatch in x, f, B0
% -4 : Maybe started at a saddle point
% -5 : Overflow during computation
% perf : Array, holding
% perf(1,:) = values of || f(x) ||inf
% perf(2,:) = Delta-values.
% Version 04.04.10. hbn(a)imm.dtu.dk
% Check parameters and function call
if nargin < 2
[stop,x,f,opts,B,D,info,n] = checkinput(fun, [], [], []);
elseif nargin < 3
[stop,x,f,opts,B,D,info,n] = checkinput(fun, x0, [], []);
elseif nargin < 4
[stop,x,f,opts,B,D,info,n] = checkinput(fun, x0, opts, []);
else
[stop,x,f,opts,B,D,info,n] = checkinput(fun, x0, opts, B0, varargin{:});
end
if stop
perf = []; X = x; return
end
% Finish initialization
Delta = opts(1); kmax = opts(4);
Trace = nargout > 2;
if Trace
X = repmat(x,1,kmax+1);
perf = repmat([norm(f,inf); Delta],1,kmax+1);
end
k = 1; nx = norm(x); nf = norm(f,inf); F = .5*norm(f)^2;
ku = 0; refac = 0; % For "extra" updates
newD = 0; % Signifies whether D should be recomputed
nstep = 0;
% Iterate
while ~stop
if isinf(nf), stop = -5;
elseif nf <= opts(2), stop = 1;
elseif Delta <= opts(3)*(opts(3) + nx), stop = 2;
elseif k >= kmax, stop = 3;
else
% Newton step
if newD, [stop D hN] = getDandh(B,f);
else, hN = D*(-f); end
nhN = norm(hN);
if nhN <= opts(3)*(opts(3) + nx), stop = 2;
elseif ~stop
if nhN > Delta % Include gradient
g = B'*(-f); ng = norm(g); alpha = (ng / norm(B*g))^2;
gn = alpha*g; ngn = alpha*ng;
if ngn >= Delta
h = (Delta/ngn) * gn;
else % Dog leg
b = hN - gn; bb = b'*b; gb = gn'*b;
c = (Delta + ngn)*(Delta - ngn);
if gb > 0
beta = c / (gb + sqrt(gb^2 + c * bb));
else
beta = (sqrt(gb^2 + c * bb) - gb)/bb;
end
h = gn + beta*b;
end
else
h = hN;
end
end
end
if ~stop
ku = mod(ku,n) + 1;
if abs(h(ku)) < .8*norm(h) % extra step for updating B and D
xu = x;
if x(ku) == 0, xu(ku) = opts(5)^2;
else, xu(ku) = x(ku) + opts(5)*abs(x(ku)); end
[stop,B,D,Fu,fu,hu,nhu,newD] = updBD(fun,x,xu,f,B,D,0,varargin{:});
info(5) = info(5)+1;
end
if ~stop
xnew = x + h;
[stop,B,D,Fn,fn,h,nstep,newD] = updBD(fun,x,xnew,f,B,D,newD,varargin{:});
info(5) = info(5)+1;
if ~stop
% Check trust region
dF = F - Fn;
if dF > 0 % Update x
dL = F - .5*norm(f + B*h)^2;
x = xnew; F = Fn; f = fn;
nf = norm(f,inf); nx = norm(x);
if dL > 0, rho = dF / dL; else, rho = 0; end
else, rho = 0; end
if rho > .75, Delta = max(Delta, 3*nstep);
elseif rho < .25, Delta = Delta/2; end
k = k + 1;
if Trace, X(:,k) = x; perf(:,k) = [nf Delta]'; end
end
end
end
end
% Set return values
if Trace
X = X(:,1:k); perf = perf(:,1:k);
else, X = x; end
if stop < 0, nf = NaN; end
info = [nf nstep Delta k-1 info(5) stop];
%%%%%%%%%%%%%%%%%%%% Auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [stop,x,f,opts,B,D,info,n] = checkinput(fun, x0, opts, B0, varargin)
% Check indata
info = [NaN zeros(1,5)]; f = NaN; B = NaN; D = NaN; n = NaN;
stop = 0;
if isempty(x0), stop = -1; x = [];
else
[stop x n] = checkx(x0);
if ~stop
[stop F f] = checkfJ(fun,x0,varargin{:});
info([1 5]) = [norm(f,inf) 1];
nf = norm(f, inf);
if ~stop & length(f) ~= n, stop = -3; end
if ~stop
% Finish initialization
opts = checkopts(opts, [.1*(1+norm(x,inf)) 1e-6 1e-8 100 1e-6]);
% Jacobian
sB = size(B0);
if sum(sB) == 0 % placeholder
[stop B] = Dapprox(fun,x,opts(5),f,varargin{:});
info(5) = info(5) + n;
elseif any(sB ~= n), stop = -3;
else, B = B0; end
% Check gradient
if ~stop
g = B'*(-f); ng = norm(g,inf);
if isinf(ng), stop = -5; end
end
% Get initial inverse Jacobian and check D*f
if ~stop
[stop D hN] = getDandh(B,f);
end
end
end
end
if stop, info(6) = stop; end
function [stop, D, hN] = getDandh(B,f)
% Get inverse of approximate Jacobian, and check for overflow
[U S V] = svd(B);
s = diag(S); i = find(s > 100*eps*s(1));
if isempty(i), stop = -4; D = NaN; hN = NaN;
else
D = V(:,i) * diag(1./s(i)) * U(:,i)';
hN = D*(-f);
if isinf(norm(hN)), stop = -5; else, stop = 0; end
end
function [stop,B,D,Fn,fn,h,nh,newD] = updBD(fun,x,xn,f,B,D,newD,varargin)
% Evaluate at new point and update B and D
[stop Fn fn] = checkfJ(fun,xn,varargin{:});
if ~stop
h = xn - x; nh = norm(h); y = fn - f;
B = B + ((y - B*h)/nh^2) * h';
if ~newD
Dy = D*y; hDy = dot(h,Dy);
if abs(hDy) < sqrt(eps)*nh, newD = 1;
else
D = D + ((h - Dy)/hDy)*(h'*D);
end
end
end
\ No newline at end of file
src/optic.dat 0 → 100644
+ 200
− 0
View file @ bd4f1101
0.00 17.4529
0.16 16.1061
0.32 14.2163
0.48 14.1534
0.64 12.8829
0.80 12.3122
0.96 12.0819
1.12 11.5606
1.28 11.0695
1.44 10.5746
1.60 9.94791
1.76 9.60335
1.92 9.31752
2.08 8.66679
2.24 8.74171
2.40 8.20382
2.56 7.95962
2.72 7.85510
2.88 7.52395
3.04 7.24599
3.20 7.00317
3.36 6.95045
3.52 6.41164
3.68 6.37972
3.84 6.21461
4.00 6.15310
4.16 5.79189
4.32 5.64620
4.48 5.31505
4.64 5.30164
4.80 4.83174
4.96 5.03015
5.12 4.50429
5.28 4.45064
5.44 4.16157
5.60 4.30957
5.76 4.32854
5.92 3.64635
6.08 3.83736
6.24 3.82534
6.40 3.80036
6.56 3.40770
6.72 3.41325
6.88 3.40724
7.04 3.10846
7.20 3.50621
7.36 2.83327
7.52 3.20235
7.68 2.95075
7.84 2.45310
8.00 2.70146
8.16 2.47252
8.32 2.64226
8.48 2.16450
8.64 1.99754
8.80 2.09466
8.96 1.91567
9.12 2.05951
9.28 2.12426
9.44 2.19687
9.60 2.17699
9.76 1.49295
9.92 1.36021
10.08 1.40184
10.24 1.95915
10.40 1.16781
10.56 2.03454
10.72 1.33061
10.88 1.38842
11.04 1.36576
11.20 1.29500
11.36 1.07069
11.52 0.93425
11.68 1.50590
11.84 1.13312
12.00 0.71918
12.16 1.09057
12.32 0.96107
12.48 0.72843
12.64 1.04109
12.80 0.62160
12.96 0.43058
13.12 0.58460
13.28 0.84591
13.44 0.41301
13.60 1.02259
13.76 0.46250
13.92 1.03739
14.08 0.37508
14.24 0.44261
14.40 0.15308
14.56 0.29276
14.72 0.24420
14.88 0.18037
15.04 0.08417
15.20 -0.20766
15.36 -0.29230
15.52 0.52910
15.68 -0.07353
15.84 -0.03561
16.00 -0.09620
16.16 0.03885
16.32 -0.22107
16.48 0.20720
16.64 0.54575
16.80 -0.54713
16.96 -0.07585
17.12 0.05642
17.28 0.16973
17.44 0.22153
17.60 0.25021
17.76 -0.37508
17.92 -0.25668
18.08 -0.42503
18.24 -0.85608
18.40 -0.30293
18.56 -0.40283
18.72 -0.23680
18.88 -0.60772
19.04 -0.42272
19.20 -0.26501
19.36 -0.48423
19.52 -0.76451
19.68 -0.78856
19.84 -0.51568
20.00 -0.33855
20.16 -0.61050
20.32 -1.04386
20.48 -0.16048
20.64 -1.03415
20.80 -0.58645
20.96 -0.79920
21.12 -0.52262
21.28 -0.76960
21.44 -0.49903
21.60 -0.96200
21.76 -0.59385
21.92 -1.20157
22.08 -0.73861
22.24 -0.48007
22.40 -0.77746
22.56 -0.73861
22.72 -1.10676
22.88 -0.59015
23.04 -0.93101
23.20 -0.39636
23.36 -1.01842
23.52 -0.38110
23.68 -0.72566
23.84 -0.34918
24.00 -1.25337
24.16 -1.02906
24.32 -0.91667
24.48 -1.12526
24.64 -0.50505
24.80 -0.81816
24.96 -1.20712
25.12 -0.81816
25.28 -0.99345
25.44 -1.15856
25.60 -1.22054
25.76 -0.67895
25.92 -1.02490
26.08 -1.06837
26.24 -1.24274
26.40 -1.64187
26.56 -1.08780
26.72 -0.74601
26.88 -1.20620
27.04 -0.95228
27.20 -0.96986
27.36 -1.01010
27.52 -0.42735
27.68 -1.16550
27.84 -1.25754
28.00 -1.63586
28.16 -1.11694
28.32 -0.98096
28.48 -1.46381
28.64 -1.19047
28.80 -1.64881
28.96 -1.39351
29.12 -0.89308
29.28 -1.44947
29.44 -1.12526
29.60 -1.10260
29.76 -1.56834
29.92 -1.41062
30.08 -1.53365
30.24 -1.26817
30.40 -0.51800
30.56 -0.97125
30.72 -1.52209
30.88 -1.50081
31.04 -0.91528
31.20 -0.57858
31.36 -0.94257
31.52 -1.27511
31.68 -1.23117
31.84 -0.90095
src/osl.dat 0 → 100644
+ 0
− 0
View file @ bd4f1101
This diff is collapsed.
src/peaks.dat 0 → 100644
+ 200
− 0
View file @ bd4f1101
0.025 1.014940
0.050 1.108103
0.075 1.142839
0.100 1.246959
0.125 1.138845
0.150 1.074924
0.175 1.280720
0.200 1.209563
0.225 1.158840
0.250 1.204884
0.275 1.100952
0.300 1.234133
0.325 1.186674
0.350 1.227846
0.375 1.081482
0.400 1.331683
0.425 1.240807
0.450 1.035223
0.475 1.271694
0.500 1.363072
0.525 1.322924
0.550 1.413546
0.575 1.432737
0.600 1.254249
0.625 1.157762
0.650 1.577044
0.675 1.492731
0.700 1.512275
0.725 1.537095
0.750 1.509848
0.775 1.767650
0.800 1.735955
0.825 1.931329
0.850 2.174533
0.875 2.501963
0.900 2.778346
0.925 3.034479
0.950 3.862237
0.975 4.504688
1.000 5.744424
1.025 6.884948
1.050 8.648960
1.075 10.262895
1.100 11.024679
1.125 10.249237
1.150 8.643818
1.175 6.793879
1.200 5.516481
1.225 4.555589
1.250 3.743961
1.275 3.187197
1.300 2.782900
1.325 2.505164
1.350 2.213014
1.375 2.017841
1.400 1.993876
1.425 1.800502
1.450 1.776591
1.475 1.599553
1.500 1.622683
1.525 1.664446
1.550 1.481009
1.575 1.553224
1.600 1.501848
1.625 1.506422
1.650 1.475327
1.675 1.315966
1.700 1.501002
1.725 1.653921
1.750 1.599891
1.775 1.524022
1.800 1.522410
1.825 1.637173
1.850 1.614317
1.875 1.643041
1.900 1.746422
1.925 1.914642
1.950 1.629529
1.975 1.993535
2.000 1.932339
2.025 2.213690
2.050 2.015819
2.075 2.353167
2.100 2.533129
2.125 2.666678
2.150 2.894890
2.175 3.015484
2.200 3.335135
2.225 3.697025
2.250 4.198146
2.275 4.569405
2.300 5.305350
2.325 5.801749
2.350 6.333429
2.375 7.266943
2.400 8.110484
2.425 9.049076
2.450 9.810732
2.475 10.419047
2.500 10.724160
2.525 10.829967
2.550 10.473012
2.575 10.245331
2.600 9.552201
2.625 9.083635
2.650 8.632697
2.675 8.223566
2.700 7.805174
2.725 7.638578
2.750 7.337614
2.775 6.991945
2.800 6.719296
2.825 6.437725
2.850 6.293143
2.875 6.047208
2.900 6.001650
2.925 5.866614
2.950 5.961151
2.975 5.769354
3.000 5.617959
3.025 5.507371
3.050 5.473228
3.075 5.473841
3.100 5.551242
3.125 5.315186
3.150 5.303839
3.175 5.219292
3.200 5.130185
3.225 5.080369
3.250 4.995783
3.275 4.754496
3.300 4.571512
3.325 4.450684
3.350 4.331275
3.375 4.143868
3.400 3.760158
3.425 3.667662
3.450 3.337121
3.475 3.095104
3.500 2.957190
3.525 2.784806
3.550 2.449690
3.575 2.182533
3.600 2.151452
3.625 2.057031
3.650 1.877548
3.675 1.894320
3.700 1.607671
3.725 1.614608
3.750 1.380580
3.775 1.258307
3.800 1.368161
3.825 1.114716
3.850 1.067723
3.875 1.082565
3.900 0.955738
3.925 1.105710
3.950 0.991701
3.975 0.864384
4.000 0.872222
4.025 0.835979
4.050 0.896251
4.075 0.655884
4.100 0.794968
4.125 0.711081
4.150 0.853375
4.175 0.862044
4.200 0.804956
4.225 0.767763
4.250 0.519182
4.275 0.754353
4.300 0.658779
4.325 0.704697
4.350 0.688704
4.375 0.601854
4.400 0.754242
4.425 0.713936
4.450 0.650790
4.475 0.758422
4.500 0.814667
4.525 0.557820
4.550 0.651485
4.575 0.544053
4.600 0.710648
4.625 0.681199
4.650 0.580708
4.675 0.571060
4.700 0.547449
4.725 0.698090
4.750 0.583716
4.775 0.587784
4.800 0.479798
4.825 0.620946
4.850 0.837167
4.875 0.609467
4.900 0.445811
4.925 0.596537
4.950 0.435846
4.975 0.578822
\ No newline at end of file
src/private/Dapprox.m 0 → 100644
+ 16
− 0
View file @ bd4f1101
function [err, J] = Dapprox(fun,x,d,f,varargin)
% Approximate Jacobian by forward differences
J = zeros(length(f),length(x));
xx = x;
for j = 1 : length(x)
if x(j) == 0, xp = d^2;
else, xp = x(j) + d*abs(x(j)); end
xx(j) = xp; fp = feval(fun,xx,varargin{:});
J(:,j) = (fp - f)/(xp - x(j));
xx(j) = x(j);
end
% Check J
if ~isreal(J) | any(isnan(J(:))) | any(isinf(J(:)))
err = -6;
else, err = 0; end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment