[Master Index] [Index for PublicToolbox/regutools]

# 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

This function calls

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

This function is called by

    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

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON