[Master Index] [Index for PublicToolbox/regutools]

# 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

This function calls

    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

This function is called by

    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

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON