heb_new.m

求解离散病态问题的正则化方法matlab 工具箱

function lambda = heb_new(lambda_0,alpha,s,beta,omega)
%HEB_NEW Newton iteration with Hebden model (utility routine for LSQI).
%
% lambda = heb_new(lambda_0,alpha,s,beta,omega)
%
% Uses Newton iteration with a Hebden (rational) model to find the
% solution lambda to the secular equation
%    || L (x_lambda - x_0) || = alpha ,
% where x_lambda is the solution defined by Tikhonov regularization.
%
% The initial guess is lambda_0.
%
% The norm || L (x_lambda - x_0) || is computed via s, beta and omega.
% 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.

% Reference: T. F. Chan, J. Olkin & D. W. Cooley, "Solving quadratically
% constrained least squares using block box unconstrained solvers",
% BIT 32 (1992), 481-495.
% Extension to the case x_0 ~= 0 by Per Chr. Hansen, IMM, 11/20/91.

% Per Christian Hansen, IMM, 12/29/97.

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

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

% Iterate, using Hebden-Newton iteration, i.e., solve the nonlinear
% problem || L x ||^(-2) - alpha^(-2) = 0.  This version was found
% experimentally to work slight鎦 better than Newton's method for
% alpha-values near || L x^exact ||.
lambda = lambda_0; step = 1; it = 0;
while (abs(step) > thr*lambda & it < it_max), it = it+1;
  e = s./(s2 + lambda^2); f = s.*e;
  if (ps==1)
    Lx = e.*beta - f.*omega;
  else
    Lx = e.*beta - f.*mu.*omega;
  end
  norm_Lx = norm(Lx);
  Lv = lambda^2*Lx./(s2 + lambda^2);
  step = (lambda/4)*(norm_Lx^2 - alpha^2)/(Lv'*Lx); % Newton step.
  step = (norm_Lx^2/alpha^2)*step; % Hebden step.
  lambda = lambda + step;
  if (lambda < 0), lambda = 2*lambda_0; lambda_0 = 2*lambda_0; end
end

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