newton

(PublicToolbox/regutools/newton.m in BrainStorm 2.0 (Alpha))

Function Synopsis

lambda = newton(lambda_0,delta,s,beta,omega,delta_0)

Help Text

NEWTON Newton iteration (utility routine for DISCREP).

 lambda = newton(lambda_0,delta,s,beta,omega,delta_0)

 Uses Newton iteration to find the solution lambda to the equation
    || A x_lambda - b || = delta ,
 where x_lambda is the solution defined by Tikhonov regularization.

 The initial guess is lambda_0.

 The norm || A x_lambda - b || is computed via s, beta, omega and
 delta_0.  Here, s holds either the singular values of A, if L = I,
 or the c,s-pairs of the GSVD of (A,L), if L ~= I.  Moreover,
 beta = U'*b and omega is either V'*x_0 or the first p elements of
 inv(X)*x_0.  Finally, delta_0 is the incompatibility measure.

Cross-Reference Information

discrep                            C:\BrainStorm_2001\PublicToolbox\regutools\discrep.m

Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\newton.m

function lambda = newton(lambda_0,delta,s,beta,omega,delta_0)
%NEWTON Newton iteration (utility routine for DISCREP).
%
% lambda = newton(lambda_0,delta,s,beta,omega,delta_0)
%
% Uses Newton iteration to find the solution lambda to the equation
%    || A x_lambda - b || = delta ,
% where x_lambda is the solution defined by Tikhonov regularization.
%
% The initial guess is lambda_0.
%
% The norm || A x_lambda - b || is computed via s, beta, omega and
% delta_0.  Here, s holds either the singular values of A, if L = I,
% or the c,s-pairs of the GSVD of (A,L), if L ~= I.  Moreover,
% beta = U'*b and omega is either V'*x_0 or the first p elements of
% inv(X)*x_0.  Finally, delta_0 is the incompatibility measure.

% Reference: V. A. Morozov, "Methods for Solving Incorrectly Posed
% Problems", Springer, 1984; Chapter 26.

% Per Christian Hansen, UNI-C, 02/21/92.

% Set defaults.
thr = 100*sqrt(eps);  % Relative stopping criterion.
it_max = 50;      % Max number of iterations.
% lambda_0 = eps;

% Initialization.
if (lambda_0 < 0)
  error('Initial guess lambda_0 must be nonnegative')
end
[p,ps] = size(s);
if (ps==2), sigma = s(:,1); mu = s(:,2); s = s(:,1)./s(:,2); end
s2 = s.^2;

% Iterate; avoid negative values of lambda.
lambda = lambda_0^2; step = 1; it = 0;
while (abs(step) > thr*lambda & it < it_max), it = it+1;
  f = s2./(s2 + lambda);
  if (ps==1)
    r  = (1-f).*(beta - s.*omega);
    dr = f.*(beta - omega)./(s2 + lambda);
  else
    r  = (1-f).*(beta - sigma.*omega);
    dr = f.*(beta - sigma.*omega)./(s2 + lambda);
  end
  res = sqrt(r'*r + delta_0^2);
  step = -(res - delta)*res/(dr'*r);
  lambda = lambda + step;
  if (lambda <= 0), lambda = eps*max(s2); end
end

% Terminate with an error if too many iterations.
if (abs(step) > thr*lambda)
  error('Max. number of iterations reached')
end

lambda = sqrt(lambda);