discrep.m

A comparison of methods for inverting helioseismic data

function [x_delta,lambda] = discrep(U,s,V,b,delta,x_0)
%DISCREP Discrepancy principle criterion for choosing the reg. parameter.
%
% [x_delta,lambda] = discrep(U,s,V,b,delta,x_0)
% [x_delta,lambda] = discrep(U,sm,X,b,delta,x_0)  ,  sm = [sigma,mu]
%
% Least squares minimization with a quadratic inequality constraint:
%    min || x - x_0 ||       subject to   || A x - b || <= delta
%    min || L (x - x_0) ||   subject to   || A x - b || <= delta
% where x_0 is an initial guess of the solution, and delta is a
% positive constant.  Requires either the compact SVD of A saved as
% U, s, and V, or part of the GSVD of (A,L) saved as U, sm, and X.
% The regularization parameter lambda is also returned.
%
% If delta is a vector, then x_delta is a matrix such that
%    x_delta = [ x_delta(1), x_delta(2), ... ] .
%
% If x_0 is not specified, x_0 = 0 is used.

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

% Per Christian Hansen, IMM, August 6, 2007.

% Initialization.
m = size(U,1);          n = size(V,1);
[p,ps] = size(s);       ld  = length(delta);
x_delta = zeros(n,ld);  lambda = zeros(ld,1);  rho = zeros(p,1);
if (min(delta)<0)
  error('Illegal inequality constraint delta')
end
if (nargin==5), x_0 = zeros(n,1); end
if (ps == 1), omega = V'*x_0; else omega = V\x_0; end

% Compute residual norms corresponding to TSVD/TGSVD.
beta = U'*b;
if (ps == 1)
  delta_0 = norm(b - U*beta);
  rho(p) = delta_0^2;
  for i=p:-1:2
    rho(i-1) = rho(i) + (beta(i) - s(i)*omega(i))^2;
  end
else
  delta_0 = norm(b - U*beta);
  rho(1) = delta_0^2;
  for i=1:p-1
    rho(i+1) = rho(i) + (beta(i) - s(i,1)*omega(i))^2;
  end
end

% Check input.
if (min(delta) < delta_0)
  error('Irrelevant delta < || (I - U*U'')*b ||')
end

% Determine the initial guess via rho-vector, then solve the nonlinear
% equation || b - A x ||^2 - delta_0^2 = 0 via Newton's method.
if (ps == 1)
    
  % The standard-form case.
  s2 = s.^2;
  for k=1:ld
    if (delta(k)^2 >= norm(beta - s.*omega)^2 + delta_0^2)
      x_delta(:,k) = x_0;
    else
      [dummy,kmin] = min(abs(rho - delta(k)^2));
      lambda_0 = s(kmin);
      lambda(k) = newton(lambda_0,delta(k),s,beta,omega,delta_0);
      e = s./(s2 + lambda(k)^2); f = s.*e;
      x_delta(:,k) = V(:,1:p)*(e.*beta + (1-f).*omega);
    end
  end
  
elseif (m>=n)
    
  % The overdetermined or square genera-form case.
  omega = omega(1:p); gamma = s(:,1)./s(:,2);
  x_u   = V(:,p+1:n)*beta(p+1:n);
  for k=1:ld
    if (delta(k)^2 >= norm(beta(1:p) - s(:,1).*omega)^2 + delta_0^2)
      x_delta(:,k) = V*[omega;U(:,p+1:n)'*b];
    else
      [dummy,kmin] = min(abs(rho - delta(k)^2));
      lambda_0 = gamma(kmin);
      lambda(k) = newton(lambda_0,delta(k),s,beta(1:p),omega,delta_0);
      e = gamma./(gamma.^2 + lambda(k)^2); f = gamma.*e;
      x_delta(:,k) = V(:,1:p)*(e.*beta(1:p)./s(:,2) + ...
                               (1-f).*s(:,2).*omega) + x_u;
    end
  end
  
else

  % The underdetermined general-form case.
  omega = omega(1:p); gamma = s(:,1)./s(:,2);
  x_u   = V(:,p+1:m)*beta(p+1:m);
  for k=1:ld
    if (delta(k)^2 >= norm(beta(1:p) - s(:,1).*omega)^2 + delta_0^2)
      x_delta(:,k) = V*[omega;U(:,p+1:m)'*b];
    else
      [dummy,kmin] = min(abs(rho - delta(k)^2));
      lambda_0 = gamma(kmin);
      lambda(k) = newton(lambda_0,delta(k),s,beta(1:p),omega,delta_0);
      e = gamma./(gamma.^2 + lambda(k)^2); f = gamma.*e;
      x_delta(:,k) = V(:,1:p)*(e.*beta(1:p)./s(:,2) + ...
                               (1-f).*s(:,2).*omega) + x_u;
    end
  end
end

%-------------------------------------------------------------------

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, 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), sigma = s(:,1); s = s(:,1)./s(:,2); end
s2 = s.^2;

% Use Newton's method to solve || b - A x ||^2 - delta^2 = 0.
% It was found experimentally, that this formulation is superior
% to the formulation || b - A x ||^(-2) - delta^(-2) = 0.
lambda = lambda_0; step = 1; it = 0;
while (abs(step) > thr*lambda & abs(step) > thr & it < it_max), it = it+1;
  f = s2./(s2 + lambda^2);
  if (ps==1)
    r = (1-f).*(beta - s.*omega);
    z = f.*r;
  else
    r = (1-f).*(beta - sigma.*omega);
    z = f.*r;
  end
  step = (lambda/4)*(r'*r + (delta_0+delta)*(delta_0-delta))/(z'*r);
  lambda = lambda - step;
  % If lambda < 0 then restart with smaller initial guess.
  if (lambda < 0), lambda = 0.5*lambda_0; lambda_0 = 0.5*lambda_0; end
end

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