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