discrep.m 5.52 KB
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
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