std_form.m

这是在网上下的一个东东

function [A_s,b_s,L_p,K,M] = std_form(A,L,b,W) 
%STD_FORM Transform a general-form reg. problem into one in standard form. 
% 
% [A_s,b_s,L_p,K,M] = std_form(A,L,b)      (method 1) 
% [A_s,b_s,L_p,x_0] = std_form(A,L,b,W)    (method 2) 
% 
% Transforms a regularization problem in general form 
%    min { || A x - b ||^2 + lambda^2 || L x ||^2 } 
% into one in standard form 
%    min { || A_s x_s - b_s ||^2 + lambda^2 || x_s ||^2 } . 
% 
% Two methods are available.  In both methods, the regularized 
% solution to the original problem can be written as 
%    x = L_p*x_s + d 
% where L_p and d depend on the method as follows: 
%    method = 1: L_p = pseudoinverse of L, d  = K*M*(b - A*L_p*x_s) 
%    method = 2: L_p = A-weighted pseudoinverse of L, d = x_0. 
% 
% The transformation from x_s back to x can be carried out by means 
% of the subroutine gen_form. 
 
% Reference: P. C. Hansen, "Rank-Deficient and Discrete Ill-PosedProblems. 
% Numerical Aspects of Linear Inversion", SIAM, Philadelphia, 1997. 
 
% Per Christian Hansen, IMM, April 8, 2001. 
 
% Nargin determines which method. 
if (nargin==3) 
 
  % Initialization for method 1. 
  [m,n] = size(A); [p,np] = size(L); 
  if (np~=n), error('A and L must have the same number of columns'), end 
 
  % Special treatment of the case where L is square. 
  if (p==n) 
    L_p = inv(L); K = []; M = []; A_s = A/L; b_s = b; 
    return 
  end 
 
  % Compute a QR factorization of L'. 
  [K,R] = qr(full(L')); R = R(1:p,:); 
 
  % Compute a QR factorization of A*K(:,p+1:n)). 
  [H,T] = qr(A*K(:,p+1:n)); T = T(1:n-p,:); 
 
  % Compute the transformed quantities. 
  L_p = K(:,1:p)/R';   %(R\(K(:,1:p)'))'; 
  K   = K(:,p+1:n); 
  M   = T\(H(:,1:n-p)'); 
  A_s = H(:,n-p+1:m)'*A*L_p; 
  b_s = H(:,n-p+1:m)'*b; 
 
else 
 
  % Initialization for method 2. 
  [m,n] = size(A); [p,nl] = size(L); nu = n-p; 
  if (nl~=n), error('A and L must have the same number of columns'), end 
 
  % Special treatment of the case where L is square. 
  if (p==n) 
    L_p = inv(L); A_s = A/L; b_s = b; 
    x_0 = zeros(n,1); K = x_0; % Fix output name. 
    return 
  end 
 
  % Compute T and x_0; 
  [T,x_0] = pinit(W,A,b); 
  b_s = b - A*x_0; 
 
  % Compute the transformed quantities. 
  L1 = inv(L(:,1:p)); 
  L_p = [L1;zeros(nu,p)] - W*(T(:,1:p)*L1); 
  A_s = A*L_p; 
 
  % Fix output name. 
  K = x_0; 
 
end