quasiopt

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

Function Synopsis

[reg_min,Q,reg_param] = quasiopt(U,s,b,method)

Help Text

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.

Cross-Reference Information

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

Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\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.

% Reference: T. Kitagawa, "A deterministic approach to optimal
% regularization - the finite dimensional case", Japan J. Appl.
% Math. 4 (1987), 371-391.

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

% Set defaults.
npoints = 80;  % 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 (method(1:4)=='Tikh' | method(1:4)=='tikh')

  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
    f = (s.^2)./(s.^2 + reg_param(i)^2);
    Q(i) = norm((1 - f).*f.*xi);
  end

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

  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
    f = s./(s + reg_param(i));
    Q(i) = norm((1 - f).*f.*xi);
  end

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

  Q = abs(xi); reg_param = [1:p]';

else, error('Illegal method'), end

% Plot the function.
if (method=='tsvd' | method=='tgsv')
  semilogy(reg_param,Q,'o'), xlabel('k')
  title('Quasi-optimality function')
  if (find_min)
    [minQ,minQi] = min(Q); reg_min = reg_param(minQi);
    HoldState = ishold; hold on;
    semilogy([reg_min,reg_min],[minQ,minQ/1000],'--')
    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')
  else
    loglog(reg_param,Q,':'), xlabel('lambda')
  end
  title('Quasi-optimality function')
  if (find_min)
    [minQ,minQi] = min(Q); reg_min = reg_param(minQi);
    HoldState = ishold; hold on;
    loglog([reg_min,reg_min],[minQ,minQ/1000],'--')
    if (~HoldState), hold off; end
    title(['Quasi-optimality function, minimum at ',num2str(reg_min)])
  end
end