www.pudn.com > gcv.rar > gcv.m

function [rpar, G] = gcv(U, s, g, method) 
% GCV  Generalized cross-validation.广义交叉校验 
% 
%    Given a matrix of left singular vectors U, a vector of singular 
%    values s, a data vector g, and a regularization method, 
% 
%        [rpar, G] = gcv(U, s, g, method); 
% 
%    returns a regularization parameter rpar chosen by generalized 
%    cross-validation. Also returned is the value G of the generalized 
%    cross-validation function 
% 
%                 RSS 
%            G =  --- 
%                 T^2 
% 
%    where T = n - sum(f_i) is an effective number of degrees of 
%    freedom and f_i are the filter factors of the regularization 
%    method. 
% 
%    The regularization method must be one of the following: 
% 
%     method = 'ridge' or 'Tikhonov': ridge regression/Tikhonov reg. 
%     method = 'tsvd' :               truncated SVD 
% 
%    If no method is specified, ridge regression is performed. The 
%    returned regularization parameter rpar is the ridge parameter of 
%    ridge regression or the truncation parameter of the truncated 
%    singular value decomposition. 
 
% References: 
%     Hansen, P. C., 1998: Rank-Deficient and Discrete Ill-Posed 
%        Problems. SIAM Monogr. on Mathematical Modeling and 
%        Computation, SIAM. 
% 
%     Wahba, G., 1990: Spline Models for Observational Data. CBMS-NSF 
%        Regional Conference Series in Applied Mathematics, Vol. 59, 
%        SIAM. 
% 
%     Adapted from various routines in Per-Christian Hansen's 
%     regularization toolbox. 
 
  % Check inputs 
  error(nargchk(3, 4, nargin)) 
 
  if nargin == 3 
    method = 'ridge';  % default regularization method 
  end 
 
  [n, q]      = size(U); 
  s           = s(:); 
  if length(s) ~= q 
    error('Input s must be vector of singular values.'); 
  end 
 
  % Coefficients in expansion of solution in terms of right singular 
  % vectors 
  fc          = U(:, 1:q)'*g; 
  s2          = s.^2; 
 
  % Least squares residual 
  rss0        = 0; 
  if n > q 
    rss0      = sum((g - U(:, 1:q)*fc).^2); 
  end 
 
  switch lower(method) 
   case {'ridge', 'tikhonov'} 
    % Accuracy of regularization parameter 
    h_tol     = ((q^2 + q + 1)*eps)^(1/2); 
 
    % Heuristic upper bound on regularization parameter 
    h_max     = max(s); 
 
    % Heuristic lower bound 
    h_min     = min(s) * h_tol; 
 
    % Find minimizer of GCV function 
    minopt    = optimset('TolX', h_tol, 'Display', 'off'); 
    [rpar, G] = fminbnd('gcvfctn', h_min, h_max, minopt, s2, fc, ... 
			rss0, n-q); 
   case {'tsvd'} 
    % compute residual sum of squares for truncation parameters 1,...,q 
    rss       = flipud(cumsum(flipud([fc(2:q).^2; rss0]))); 
    if n > q 
      [G, rpar] = min(rss./(n-[1:q])'.^2); 
    else 
      [G, rpar] = min(rss(1:n-1)./(n-[1:(n-1)])'.^2); 
    end 
   otherwise 
    error('Unknown regularization method.') 
  end