l_corner

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

Function Synopsis

[reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M)

Help Text

L_CORNER Locate the "corner" of the L-curve.

 [reg_c,rho_c,eta_c] =
        l_corner(rho,eta,reg_param)
        l_corner(rho,eta,reg_param,U,s,b,method,M)
        l_corner(rho,eta,reg_param,U,sm,b,method,M) ,  sm = [sigma,mu]

 Locates the "corner" of the L-curve in log-log scale.

 It is assumed that corresponding values of || A x - b ||, || L x ||,
 and the regularization parameter are stored in the arrays rho, eta,
 and reg_param, respectively (such as the output from routine l_curve).

 If nargin = 3, then no particular method is assumed, and if
 nargin = 2 then it is issumed that reg_param = 1:length(rho).

 If nargin >= 6, then the following methods are allowed:
    method = 'Tikh'  : Tikhonov regularization
    method = 'tsvd'  : truncated SVD or GSVD
    method = 'dsvd'  : damped SVD or GSVD
    method = 'mtsvd' : modified TSVD,
 and if no method is specified, 'Tikh' is default.  If the Spline Toolbox
 is not available, then only 'Tikh' and 'dsvd' can be used.

 An eighth argument M specifies an upper bound for eta, below which
 the corner should be found.

Cross-Reference Information

dsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\dsvd.m
fnder                              C:\BrainStorm_2001\PublicToolbox\regutools\fnder.m
ppcut                              C:\BrainStorm_2001\PublicToolbox\regutools\ppcut.m
sp2pp                              C:\BrainStorm_2001\PublicToolbox\regutools\sp2pp.m
spleval                            C:\BrainStorm_2001\PublicToolbox\regutools\spleval.m
spmak                              C:\BrainStorm_2001\PublicToolbox\regutools\spmak.m
tsvd                               C:\BrainStorm_2001\PublicToolbox\regutools\tsvd.m

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

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

function [reg_c,rho_c,eta_c] = l_corner(rho,eta,reg_param,U,s,b,method,M)
%L_CORNER Locate the "corner" of the L-curve.
%
% [reg_c,rho_c,eta_c] =
%        l_corner(rho,eta,reg_param)
%        l_corner(rho,eta,reg_param,U,s,b,method,M)
%        l_corner(rho,eta,reg_param,U,sm,b,method,M) ,  sm = [sigma,mu]
%
% Locates the "corner" of the L-curve in log-log scale.
%
% It is assumed that corresponding values of || A x - b ||, || L x ||,
% and the regularization parameter are stored in the arrays rho, eta,
% and reg_param, respectively (such as the output from routine l_curve).
%
% If nargin = 3, then no particular method is assumed, and if
% nargin = 2 then it is issumed that reg_param = 1:length(rho).
%
% If nargin >= 6, then the following methods are allowed:
%    method = 'Tikh'  : Tikhonov regularization
%    method = 'tsvd'  : truncated SVD or GSVD
%    method = 'dsvd'  : damped SVD or GSVD
%    method = 'mtsvd' : modified TSVD,
% and if no method is specified, 'Tikh' is default.  If the Spline Toolbox
% is not available, then only 'Tikh' and 'dsvd' can be used.
%
% An eighth argument M specifies an upper bound for eta, below which
% the corner should be found.

% The following routines from the Spline Toolbox are needed if
% method differs from 'Tikh' or 'dsvd':
% fnder, ppbrk, ppcut, ppmak, sp2pp, spbrk, spmak.

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

% Set default regularization method.
if (nargin <= 3)
  method = 'none';
  if (nargin==2), reg_param = [1:length(rho)]'; end
else
  if (nargin==6), method = 'Tikh'; end
end

% Set threshold for skipping very small singular values in the
% L-curve analysis.
s_thr = eps;  % Neglect singular values less than s_thr.

% Set default parameters for treatment of discrete L-curve.
deg   = 2;  % Degree of local smooting polynomial.
q     = 2;  % Half-width of local smoothing interval.
order = 4;  % Order of fitting 2-D spline curve.

% Initialization.
if (length(rho) < order)
  error('Too few data points for L-curve analysis')
end
if (nargin > 3)
  [p,ps] = size(s); [m,n] = size(U);
  if (ps==2), s = s(p:-1:1,1)./s(p:-1:1,2); U = U(:,p:-1:1); end
  beta = U'*b; xi = beta./s;
end

% Restrict the analysis of the L-curve according to M (if specified)
% and s_thr.
if (nargin==8)
  index = find(eta < M);
  rho = rho(index); eta = eta(index); reg_param = reg_param(index);
  s = s(index); beta = beta(index); xi = xi(index);
end

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

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Initialization.
  [reg_m,reg_n] = size(reg_param);
  phi = zeros(reg_m,reg_n); dphi = phi; psi = phi; dpsi = phi;
  s2 = s.^2; beta2 = beta.^2; xi2 = xi.^2;

  % Compute some intermediate quantities.
  for i = 1:length(reg_param)
    f  = s2./(s2 + reg_param(i)^2); cf = 1 - f;
    f1 = -2*f.*cf/reg_param(i);
    f2 = -f1.*(3-4*f)/reg_param(i);
    phi(i)  = sum(f.*f1.*xi2);
    psi(i)  = sum(cf.*f1.*beta2);
    dphi(i) = sum((f1.^2 + f.*f2).*xi2);
    dpsi(i) = sum((-f1.^2 + cf.*f2).*beta2);
  end

  % Now compute the first and second derivatives of eta and rho
  % with respect to lambda;
  deta  =  phi./eta;
  drho  = -psi./rho;
  ddeta =  dphi./eta - deta.*(deta./eta);
  ddrho = -dpsi./rho - drho.*(drho./rho);

  % Convert to derivatives of log(eta) and log(rho).
  dlogeta  = deta./eta;
  dlogrho  = drho./rho;
  ddlogeta = ddeta./eta - (dlogeta).^2;
  ddlogrho = ddrho./rho - (dlogrho).^2;

  % Let g = curvature.
  g = (dlogrho.*ddlogeta - ddlogrho.*dlogeta)./...
      (dlogrho.^2 + dlogeta.^2).^(1.5);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmax,gi] = max(g);
  if (gmax < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    rho_c = rho(gi); eta_c = eta(gi); reg_c = reg_param(gi);
  end

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

  % The L-curve is discrete and may include unwanted fine-grained
  % corners.  Use local smoothing, followed by fitting a 2-D spline
  % curve to the smoothed discrete L-curve.

  % Check if the Spline Toolbox exists, otherwise return.
  if (exist('spdemos')~=2)
    error('The Spline Toolbox in not available so l_corner cannot be used')
  end

  % For TSVD, TGSVD, and MTSVD, restrict the analysis of the L-curve
  % according to s_thr.
  if (nargin > 3)
    index = find(s > s_thr);
    rho = rho(index); eta = eta(index); reg_param = reg_param(index);
    s = s(index); beta = beta(index); xi = xi(index);
  end

  % Convert to logarithms.
  lr = length(rho);
  lrho = log(rho); leta = log(eta); slrho = lrho; sleta = leta;

  % For all interior points k = q+1:length(rho)-q-1 on the discrete
  % L-curve, perform local smoothing with a polynomial of degree deg
  % to the points k-q:k+q.
  v = [-q:q]'; A = zeros(2*q+1,deg+1); A(:,1) = ones(length(v),1);
  for j = 2:deg+1, A(:,j) = A(:,j-1).*v; end
  for k = q+1:lr-q-1
    cr = A\lrho(k+v); slrho(k) = cr(1);
    ce = A\leta(k+v); sleta(k) = ce(1);
  end

  % Fit a 2-D spline curve to the smoothed discrete L-curve.
  sp = spmak([1:lr+order],[slrho';sleta']);
  pp = ppcut(sp2pp(sp),[4,lr+1]);

  % Extract abscissa and ordinate splines and differentiate them.
  P     = spleval(pp);  dpp   = fnder(pp);
  D     = spleval(dpp); ddpp  = fnder(pp,2);
  DD    = spleval(ddpp);
  ppx   = P(1,:);       ppy   = P(2,:);
  dppx  = D(1,:);       dppy  = D(2,:);
  ddppx = DD(1,:);      ddppy = DD(2,:);

  % Compute the corner of the spline curve via max. curvature.
  % Define curvature = 0 where both dppx and dppy are zero.
  k1    = dppx.*ddppy - ddppx.*dppy;
  k2    = (dppx.^2 + dppy.^2).^(1.5);
  I_nz  = find(k2 ~= 0);
  kappa = zeros(1,length(dppx));
  kappa(I_nz) = -k1(I_nz)./k2(I_nz);
  [kmax,ikmax] = max(kappa);
  x_corner = ppx(ikmax); y_corner = ppy(ikmax);

  % Locate the point on the discrete L-curve which is closest to the
  % corner of the spline curve.  Prefer a point below and to the
  % left of the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  if (kmax < 0)
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    index = find(lrho < x_corner & leta < y_corner);
    if (length(index) > 0)
      [dummy,rpi] = min((lrho(index)-x_corner).^2 + (leta(index)-y_corner).^2);
      rpi = index(rpi);
    else
      [dummy,rpi] = min((lrho-x_corner).^2 + (leta-y_corner).^2);
    end
    reg_c = reg_param(rpi); rho_c = rho(rpi); eta_c = eta(rpi);
  end

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

  % The L-curve is differentiable; computation of curvature in
  % log-log scale is easy.

  % Initialization.
  [reg_m,reg_n] = size(reg_param);
  phi = zeros(reg_m,reg_n); dphi = phi; psi = phi; dpsi = phi;
  beta2 = beta.^2; xi2 = xi.^2;

  % Compute some intermediate quantities.
  for i = 1:length(reg_param)
    f  = s./(s + reg_param(i)); cf = 1 - f;
    f1 = -f.*cf/reg_param(i);
    f2 = -2*f1.*cf/reg_param(i);
    phi(i)  = sum(f.*f1.*xi2);
    psi(i)  = sum(cf.*f1.*beta2);
    dphi(i) = sum((f1.^2 + f.*f2).*xi2);
    dpsi(i) = sum((-f1.^2 + cf.*f2).*beta2);
  end

  % Now compute the first and second derivatives of eta and rho
  % with respect to lambda;
  deta  =  phi./eta;
  drho  = -psi./rho;
  ddeta =  dphi./eta - deta.*(deta./eta);
  ddrho = -dpsi./rho - drho.*(drho./rho);

  % Convert to derivatives of log(eta) and log(rho).
  dlogeta  = deta./eta;
  dlogrho  = drho./rho;
  ddlogeta = ddeta./eta - (dlogeta).^2;
  ddlogrho = ddrho./rho - (dlogrho).^2;

  % Let g = curvature.
  g = (dlogrho.*ddlogeta - ddlogrho.*dlogeta)./...
      (dlogrho.^2 + dlogeta.^2).^(1.5);

  % Locate the corner.  If the curvature is negative everywhere,
  % then define the leftmost point of the L-curve as the corner.
  [gmax,gi] = max(g);
  if (gmax < 0)
    lr = length(rho);
    reg_c = reg_param(lr); rho_c = rho(lr); eta_c = eta(lr);
  else
    rho_c = rho(gi); eta_c = eta(gi); reg_c = reg_param(gi);
  end

else, error('Illegal method'), end