quasiopt.m

这是在网上下的一个东东

function [reg_min,Q,reg_param] = quasiopt(U,s,b,method) 
%QUASIOPT Quasi-optimality criterion for choosing the regularization parameter. 
% 
% [reg_min,Q,reg_param] = quasiopt(U,s,b,method) 
% [reg_min,Q,reg_param] = quasiopt(U,sm,b,method)  ,  sm = [sigma,mu] 
% 
% Plots the quasi-optimality function Q for the following methods: 
%    method = 'Tikh' : Tikhonov regularization   (solid line ) 
%    method = 'tsvd' : truncated SVD or GSVD     (o markers  ) 
%    method = 'dsvd' : damped SVD or GSVD        (dotted line) 
% If no method is specified, 'Tikh' is default. 
% 
% If any output arguments are specified, then the minimum of Q is 
% identified and the corresponding reg. parameter reg_min is returned. 
 
% Per Christian Hansen, IMM, Feb. 21, 2001. 
 
% Set defaults. 
npoints = 200;  % Number of points for 'Tikh' and 'dsvd'. 
if (nargin==3), method = 'Tikh'; end   % Default method. 
 
% Initialization. 
[m,n] = size(U); [p,ps] = size(s); 
if (ps==2), s = s(p:-1:1,1)./s(p:-1:1,2); U = U(:,p:-1:1); end 
xi = (U'*b)./s; 
if (nargout > 0), find_min = 1; else find_min = 0; end 
 
% Compute the quasioptimality function Q. 
if (strncmp(method,'Tikh',4) | strncmp(method,'tikh',4)) 
    
  % Compute a vector of Q-values. 
  Q = zeros(npoints,1); reg_param = Q; 
  reg_param(npoints) = s(p); 
  ratio = (s(1)/s(p))^(1/(npoints-1)); 
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end 
  for i=1:npoints 
    Q(i) = quasifun(reg_param(i),s,xi); 
  end 
  
 % Find the minimum, if requested. 
 if (find_min) 
   [minQ,minQi] = min(Q); % Initial guess. 
   reg_min = fminbnd('quasifun',... 
     reg_param(min(minQi+1,npoints)),reg_param(max(minQi-1,1)),... 
     optimset('Display','off'),s,xi); % Minimizer. 
   minQ = quasifun(reg_min,s,xi); % Minimum of function. 
 end 
 
elseif (strncmp(method,'dsvd',4) | strncmp(method,'dgsv',4)) 
    
  % Compute a vector of Q-values. 
  Q = zeros(npoints,1); reg_param = Q; 
  reg_param(npoints,1) = s(p); 
  ratio = (s(1)/s(p))^(1/(npoints-1)); 
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end 
  for i=1:npoints 
    Q(i) = quasifun(reg_param(i),s,xi,1); 
  end 
  
 % Find the minimum, if requested. 
 if (find_min) 
   [minQ,minQi] = min(Q); % Initial guess.  
   reg_min = fminbnd('quasifun',... 
     reg_param(min(minQi+1,npoints)),reg_param(max(minQi-1,1)),... 
     optimset('Display','off'),s,xi,1); % Minimizer. 
   minQ = quasifun(reg_min,s,xi,1); % Minimum of function. 
 end 
 
elseif (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4)) 
    
  % Compute the quasi-optimality function. 
  Q = abs(xi); reg_param = [1:p]'; 
   
  % Find the minimum, if requested. 
  if (find_min) 
    [minQ,minQi] = min(Q); reg_min = reg_param(minQi); 
  end 
 
else, error('Illegal method'), end 
 
% Plot the function. 
if (strncmp(method,'tsvd',4) | strncmp(method,'tgsv',4)) 
  semilogy(reg_param,Q,'o'), xlabel('k'), ylabel('Q(k)') 
  title('Quasi-optimality function') 
  if (find_min) 
    ax = axis; 
    HoldState = ishold; hold on; 
    semilogy([reg_min,reg_min],[minQ,minQ/1000],'--') 
    axis(ax); 
    if (~HoldState), hold off; end 
    title(['Quasi-optimality function, minimum at ',num2str(reg_min)]) 
  end 
else 
  if (method(1:4)=='Tikh' | method(1:4)=='tikh' | ... 
      method(1:4)=='dsvd' | method(1:4)=='dgsv'      ) 
    loglog(reg_param,Q), xlabel('\lambda'), ylabel('Q(\lambda)') 
  else 
    loglog(reg_param,Q,':'), xlabel('\lambda'), ylabel('Q(\lambda)') 
  end 
  title('Quasi-optimality function') 
  if (find_min) 
    ax = axis; 
    HoldState = ishold; hold on; 
    loglog([reg_min,reg_min],[minQ,minQ/1000],'--') 
    axis(ax) 
    if (~HoldState), hold off; end 
    title(['Quasi-optimality function, minimum at ',num2str(reg_min)]) 
  end 
end