gcv

(PublicToolbox/regutools/gcv.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[reg_min,G,reg_param] = gcv(U,s,b,method)

Help Text

GCV Plot the GCV function and find its minimum.

 [reg_min,G,reg_param] = gcv(U,s,b,method)
 [reg_min,G,reg_param] = gcv(U,sm,b,method)  ,  sm = [sigma,mu]

 Plots the GCV-function
          || A*x - b ||^2
    G = -------------------
        (trace(I - A*A_I)^2
 as a function of the regularization parameter reg_param.
 Here, A_I is a matrix which produces the regularized solution.

 The following methods are allowed:
    method = 'Tikh' : Tikhonov regularization   (solid line )
    method = 'tsvd' : truncated SVD or GSVD     (o markers  )
    method = 'dsvd' : damped SVD or GSVD        (dotted line)
 If method is not specified, 'Tikh' is default.

 If any output arguments are specified, then the minimum of G is
 identified and the corresponding reg. parameter reg_min is returned.

Cross-Reference Information

dsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\dsvd.m
tsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\tsvd.m

regudemo                           C:\BrainStorm_2001\PublicToolbox\regutools\regudemo.m

Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\gcv.m

function [reg_min,G,reg_param] = gcv(U,s,b,method)
%GCV Plot the GCV function and find its minimum.
%
% [reg_min,G,reg_param] = gcv(U,s,b,method)
% [reg_min,G,reg_param] = gcv(U,sm,b,method)  ,  sm = [sigma,mu]
%
% Plots the GCV-function
%          || A*x - b ||^2
%    G = -------------------
%        (trace(I - A*A_I)^2
% as a function of the regularization parameter reg_param.
% Here, A_I is a matrix which produces the regularized solution.
%
% The following methods are allowed:
%    method = 'Tikh' : Tikhonov regularization   (solid line )
%    method = 'tsvd' : truncated SVD or GSVD     (o markers  )
%    method = 'dsvd' : damped SVD or GSVD        (dotted line)
% If method is not specified, 'Tikh' is default.
%
% If any output arguments are specified, then the minimum of G is
% identified and the corresponding reg. parameter reg_min is returned.

% Per Christian Hansen, UNI-C, 03/16/93.

% Reference: G. Wahba, "Spline Models for Observational Data",
% SIAM, 1990.

% Set defaults.
if (nargin==3), method='Tikh'; end  % Default method.
npoints = 100;                      % Number of points on the curve.
smin_ratio = 16*eps;                % Smallest regularization parameter.

% Initialization.
[m,n] = size(U); [p,ps] = size(s);
beta = U'*b; beta2 = b'*b - beta'*beta;
if (ps==2)
  s = s(p:-1:1,1)./s(p:-1:1,2); beta = beta(p:-1:1);
end
if (nargout > 0), find_min = 1; else find_min = 0; end

if (method(1:4)=='Tikh' | method(1:4)=='tikh')

  reg_param = zeros(npoints,1); G = reg_param; s2 = s.^2;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  ratio = 1.2*(s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
  delta0 = 0;
  if (m > n & beta2 > 0), delta0 = beta2; end
  for i=1:npoints
    f1 = (reg_param(i)^2)./(s2 + reg_param(i)^2);
    fb = f1.*beta(1:p); rho2 = fb'*fb + delta0;
    G(i) = rho2/(m - n + sum(f1))^2;
  end 
  loglog(reg_param,G,'-'), xlabel('lambda'), ylabel('G(lambda)')
  title('GCV function')
  if (find_min)
    [minG,minGi] = min(G); reg_min = reg_param(minGi);
    HoldState = ishold; hold on;
    loglog(reg_min,minG,'*',[reg_min,reg_min],[minG/1000,minG],':')
    title(['GCV function, minimum at ',num2str(reg_min)])
    if (~HoldState), hold off; end
  end

elseif (method(1:4)=='tsvd' | method(1:4)=='tgsv')

  rho2(p-1) = beta(p)^2;
  if (m > n & beta2 > 0), rho2(p-1) = rho2(p-1) + beta2; end
  for k=p-2:-1:1, rho2(k) = rho2(k+1) + beta(k+1)^2; end
  for k=1:p-1
    G(k) = rho2(k)/(m - k + (n - p))^2;
  end
  reg_param = [1:p-1]';
  semilogy(reg_param,G,'o'), xlabel('k'), ylabel('G(k)')
  title('GCV function')
  if (find_min)
    [minG,reg_min] = min(G);
    HoldState = ishold; hold on;
    semilogy(reg_min,minG,'*',[reg_min,reg_min],[minG/1000,minG],'--')
    title(['GCV function, minimum at ',num2str(reg_min)])
    if (~HoldState), hold off; end
  end

elseif (method(1:4)=='dsvd' | method(1:4)=='dgsv')

  reg_param = zeros(npoints,1); G = reg_param;
  reg_param(npoints) = max([s(p),s(1)*smin_ratio]);
  ratio = (s(1)/reg_param(npoints))^(1/(npoints-1));
  for i=npoints-1:-1:1, reg_param(i) = ratio*reg_param(i+1); end
  delta0 = 0;
  if (m > n & beta2 > 0), delta0 = beta2; end
  for i=1:npoints
    f1 = reg_param(i)./(s + reg_param(i));
    fb = f1.*beta(1:p); rho2 = fb'*fb + delta0;
    G(i) = rho2/(m - n + sum(f1))^2;
  end 
  loglog(reg_param,G,':'), xlabel('lambda'), ylabel('G(lambda)')
  title('GCV function')
  if (find_min)
    [minG,minGi] = min(G); reg_min = reg_param(minGi);
    HoldState = ishold; hold on;
    loglog(reg_min,minG,'*',[reg_min,reg_min],[minG/1000,minG],'--')
    tiel(['GCV function, minimum at ',num2str(reg_min)])
    if (~HoldState), hold off; end
  end

elseif (method(1:4)=='mtsv')

  error('The MTSVD method is not supported')

else, error('Illegal method'), end