Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
function [X,info,perf,B] = smarquardt(fun, x0, opts, B0, varargin)
%SMARQUARDT Secant version of Levenberg-Marquardt's method for least
% squares: Find xm = argmin{f(x)} , where x = [x_1, ..., x_n] and
% f(x) = 0.5 * sum(r_i(x)^2) .
% The functions r_i(x) (i=1,...,m) must be given by a MATLAB
% function with declaration
% function r = 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] = smarquardt(fun, x0)
% [X, info] = smarquardt(fun, x0, opts)
% [X, info] = smarquardt(fun, x0, opts, B0, p1,p2,...)
% [X, info, perf] = smarquardt(.....)
% [X, info, perf, B] = smarquardt(.....)
%
% Input parameters
% fun : Handle to the function.
% x0 : Starting guess for xm .
% opts : Either a struct with fields 'tau', 'tolg', 'tolx', 'maxeval'
% and 'relstep', or a vector with the values of these options,
% opts = [tau tolg tolx maxeval relstep].
% tau used in starting value for Marquardt parameter:
% mu = tau * max( [B0'*B0]_(ii) ) ,
% where B0 is an approximate Jacobian at x0 .
% relstep "relative" step length for difference approximations.
% The other options are used in stopping criteria:
% ||B(x)'*r(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 + 10*n, relstep = 1e-7.
% 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.
% B0 : Initial approximation to J. 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 7 elements:
% info(1:4) = final values of
% [f(x) ||B(x)'*r(x)||inf ||dx||2 mu/max( [B(x)'*B(x)]_(ii) )] .
% info(5:6) = no. of iteration steps and function evaluations
% info(6) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : No. of iteration steps exceeded
% -1 : x is not a real valued vector
% -2 : r is not a real valued column vector
% -4 : Dimension mismatch in x, r, B0
% -5 : Overflow during computation
% -6 : Error in approximate Jacobian
% info(7) = no. of iterations.
% perf : Struct with fields
% f : values of f(xk) ,
% ng : values of || B'*r(x) ||_inf ,
% mu : values of Marquardt damping parameter.
% B : Computed approximation to Jacobian at the solution.
% Version 10.11.08. hbn(a)imm.dtu.dk
% Check parameters and function call
f = NaN; ng = NaN; perf = []; B = [];
info = zeros(1,7);
if nargin < 2, stop = -1;
else
[stop x n] = checkx(x0);
if ~stop
[stop f r] = checkrJ(fun,x0,varargin{:}); info(6) = 1;
if ~stop
% Finish initialization
if nargin < 3, opts = []; end
opts = checkopts('smarquardt', opts); % use default options where required
tau = opts(1); tolg = opts(2); tolx = opts(3); relstep = opts(5);
if opts(4) > 0, maxeval = opts(4); else, maxeval = 100 + 10*n; end
% Jacobian
if nargin > 3 % B0 is given
sB = size(B0); m = length(r);
if sum(sB) == 0 % placeholder
[stop B] = Dapprox(fun,x,relstep,r,varargin{:});
info(6) = info(6) + n;
elseif any(sB ~= [m n]), stop = -4;
else, B = B0; end
else
[stop B] = Dapprox(fun,x,relstep,r,varargin{:});
info(6) = info(6) + n;
end
% Check gradient and J'*J
if ~stop
g = B'*r; ng = norm(g,inf); A = B'*B;
if isinf(ng) | isinf(norm(A(:),inf)), stop = -5; end
end
end
end
end
if stop
X = x0; info([1:5 7]) = [f ng 0 tau 0 stop];
return
end
% 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
% Iterate
k = 1; nu = 2; nh = 0;
ng0 = ng;
ku = 0; % direction of last update
kit = 0; % no. of iteration steps
while ~stop
if ng <= opts(2), 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;
[stop fn rn] = checkrJ(fun,xnew,varargin{:}); info(6) = info(6)+1;
if ~stop
% Update B
ku = mod(ku,n) + 1;
if abs(h(ku)) < .8*norm(h) % extra step
xu = x;
if x(ku) == 0, xu(ku) = opts(5)^2;
else, xu(ku) = x(ku) + opts(5)*abs(x(ku)); end
[stop fu ru] = checkrJ(fun,xu,varargin{:}); info(6) = info(6)+1;
if ~stop
hu = xu - x;
B = B + ((ru - r - B*hu)/norm(hu)^2) * hu';
end
end
B = B + ((rn - r - B*h)/norm(h)^2) * h';
k = k + 1;
dL = (h'*(mu*h - g))/2;
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; r = rn;
mu = mu * max(1/3, 1 - (2*df/dL - 1)^3); nu = 2;
if Trace
X(:,kit+1) = x; perf(:,kit+1) = [fn norm(B'*rn,inf) mu]'; end
else % Same x, increase mu
mu = mu*nu; nu = 2*nu;
end
if info(5) > maxeval, stop = 3;
else
g = B'*r; ng = norm(g,inf); A = B'*B;
if isinf(ng) | isinf(norm(A(:),inf)), stop = -5; end
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, tau = NaN; else, tau = mu/max(diag(A)); end
info([1:5 7]) = [f ng nh tau kit stop];