lagrange

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

Function Synopsis

[La,dLa,lambda0] = lagrange(U,s,b,more)

Help Text

LAGRANGE Plot the Lagrange function for Tikhonov regularization.

 [La,dLa,lambda0] = lagrange(U,s,b,more)
 [La,dLa,lambda0] = lagrange(U,sm,b,more)  ,  sm = [sigma,mu]

 Plots the Lagrange function
    La(lambda) = || A x - b ||^2 + lambda^2*|| L x ||^2
 and its first derivative dLa = dLa/dlambda versus lambda.
 Here, x is the Tikhonov regularized solution.

 If nargin = 4, || A x - b || and || L x || are also plotted.

 Returns La, dLa, and the value lambda0 of lambda for which
 dLa has its minimum.

Cross-Reference Information

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

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

function [La,dLa,lambda0] = lagrange(U,s,b,more)
%LAGRANGE Plot the Lagrange function for Tikhonov regularization.
%
% [La,dLa,lambda0] = lagrange(U,s,b,more)
% [La,dLa,lambda0] = lagrange(U,sm,b,more)  ,  sm = [sigma,mu]
%
% Plots the Lagrange function
%    La(lambda) = || A x - b ||^2 + lambda^2*|| L x ||^2
% and its first derivative dLa = dLa/dlambda versus lambda.
% Here, x is the Tikhonov regularized solution.
%
% If nargin = 4, || A x - b || and || L x || are also plotted.
%
% Returns La, dLa, and the value lambda0 of lambda for which
% dLa has its minimum.

% Per Christian Hansen, UNI-C, 07/14/92.

% Initialization.
[m,n] = size(U); [p,ps] = size(s); npoints = 40;
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
xi = beta(1:p)./s;

% Compute the L-curve.
eta = zeros(npoints,1); rho = eta;
lambda(npoints,1) = s(p);
ratio = (s(1)/s(p))^(1/(npoints-1));
for i=npoints-1:-1:1, lambda(i) = ratio*lambda(i+1); end
for i=1:npoints
  f = fil_fac(s,lambda(i));
  eta(i) = norm(f.*xi);
  rho(i) = norm((1-f).*beta(1:p));
end
if (m > n & beta2 > 0), rho = sqrt(rho.^2 + beta2); end

% Compute the Lagrange function and its derivative.
La = rho.^2 + (lambda.^2).*(eta.^2);
dLa = 2*lambda.*(eta.^2);
[mindLa,mindLi] = min(dLa); lambda0 = lambda(mindLi);

% Plot the functions.
if (nargin==3)
  loglog(lambda,La,'-',lambda,dLa,'--',lambda0,mindLa,'o')
  title('--- La   - - - dLa/dlambda')
else
  loglog(lambda,La,'-',lambda,dLa,'--',lambda,eta,':',lambda,rho,'-.',...
         lambda0,mindLa,'o')
  title('--- La   - - - dLa/dlambda   ... || Lx ||   -.-. || Ax-b ||')
end
xlabel('lambda')