tikhcstr.m

这是在网上下的一个东东

function [x_lambda,rho,eta] = tikhcstr(A,b,G,d,L,lambda,x_0)
%TIKHCSTR Tikhonov regularization with linear inequality constraints.
% 
% [x_lambda,rho,eta] = tikhcstr(A,b,G,d,L,lambda,x_0)
% 
% Computes the constrained Tikhonov regularized solution x_lambda
% that solves the problem
%    min { || A x - b ||^2 + lambda^2 || L (x - x_0) ||^2 }
% subject to the linear inequality constraints
%    G x >= d . 
%
% If x_0 is not specified, then x_0 = 0 is used.
% If L is specified as the empty matrix [], then L = I is used.
%
% If lambda is a vector, then x_lambda is a matrix such that
%   x_lambda = [ x_lambda(1), x_lambda(2), ... ] .
%
% The solution and residual norms are returned in eta and rho.
%
% If G and d define an empty set of constraints, then NaNs are returned.

% Ann-Charlotte Berglund, IMM & UNI-C, June 28, 1999.

% Initialization.
[m,n] = size(A);
if nargin < 6
  error('Too few input arguments'),
end
if nargin < 7
  x_0 = zeros(n,1);
end
if min(lambda) < 0
  error('Illegal regularization parameter lambda'),
end
if isempty(L)
  L = speye(n);
end

% If no constraints then compute the ordinary Tikhonov solution.
if isempty(G) | isempty(d)
  if norm(L-speye(n),'fro')==0
    [U,s,V] = csvd(A);
    [x_lambda,rho,eta] = tikhonov(U,s,V,b,lambda,x_0);
  else
    [U,sm,X] = cgsvd(A,L);
    [x_lambda,rho,eta] = tikhonov(U,sm,X,b,lambda,x_0);
  end
  return;
end

% Otherwise prepare for computing the constrainted solution.
b = b(:);
d = d(:);
x_0 = x_0(:);
[mG,nG] = size(G);

% First try to find a feasible point.
[p,status] = fpoint([],G,d,x_0);
if ~strcmp(status,'feasible')
  % There exists no feasible point for the given constraints.
  x_lambda = NaN*ones(n,length(lambda));
  rho = NaN*ones(1,length(lambda));
  eta = NaN*ones(1,length(lambda));
  return
end

% Constrained Tikhonov regularization.
for i = 1:length(lambda)
  
  x_lambda(:,i) = tcsub([A;lambda(i)*L],[b;lambda(i)*L*x_0],...
                        G,d,x_0,0,'tikhcstr',mG,nG);

  % Residual and solution norms.
  if nargout > 1, eta(i,1) = norm(L*x_lambda(:,i)); end
  if nargout == 3, rho(i,1) = norm(A*x_lambda(:,i)-b); end
  
end % Loop over lambda-values.