[Master Index] [Index for PublicToolbox/regutools]

# heb_new

## (PublicToolbox/regutools/heb_new.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    lambda = heb_new(lambda_0,alpha,s,beta,omega)

### Help Text

    HEB_NEW Newton iteration with Hebden model (utility routine for LSQI).

     lambda = heb_new(lambda_0,alpha,s,beta,omega)

     Uses Newton iteration with a Hebden (rational) model to find the
     solution lambda to the secular equation
        || L (x_lambda - x_0) || = alpha ,
     where x_lambda is the solution defined by Tikhonov regularization.

     The initial guess is lambda_0.

     The norm || L (x_lambda - x_0) || is computed via s, beta and omega.
     Here, s holds either the singular values of A, if L = I, or the
     c,s-pairs of the GSVD of (A,L), if L ~= I.  Moreover, beta = U'*b
     and omega is either V'*x_0 or the first p elements of inv(X)*x_0.

### Cross-Reference Information

This function is called by

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

### Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\heb_new.m

    function lambda = heb_new(lambda_0,alpha,s,beta,omega)
    %HEB_NEW Newton iteration with Hebden model (utility routine for LSQI).
    %
    % lambda = heb_new(lambda_0,alpha,s,beta,omega)
    %
    % Uses Newton iteration with a Hebden (rational) model to find the
    % solution lambda to the secular equation
    %    || L (x_lambda - x_0) || = alpha ,
    % where x_lambda is the solution defined by Tikhonov regularization.
    %
    % The initial guess is lambda_0.
    %
    % The norm || L (x_lambda - x_0) || is computed via s, beta and omega.
    % Here, s holds either the singular values of A, if L = I, or the
    % c,s-pairs of the GSVD of (A,L), if L ~= I.  Moreover, beta = U'*b
    % and omega is either V'*x_0 or the first p elements of inv(X)*x_0.

    % Reference: T. F. Chan, J. Olkin & D. W. Cooley, "Solving quadratically
    % constrained least squares using block box unconstrained solvers",
    % BIT 32 (1992), 481-495.
    % Extension to the case x_0 ~= 0 by Per Chr. Hansen, UNI-C, 11/20/91.

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

    % Set defaults.
    thr = eps;     % Relative stopping criterion.
    it_max = 50;   % Max number of iterations.

    % Initialization.
    if (lambda_0 < 0)
      error('Initial guess lambda_0 must be nonnegative')
    end
    [p,ps] = size(s);
    if (ps==2), mu = s(:,2); s = s(:,1)./s(:,2); end
    s2 = s.^2;

    % Iterate.
    lambda = lambda_0^2; step = 1; it = 0;
    while (abs(step) > thr*lambda & it < it_max), it = it+1;
      e = s./(s2 + lambda); f = 2.*e;
      if (ps==1)
        Lx = e.*beta - f.*omega;
      else
        Lx = e.*beta - f.*mu.*omega;
      end
      norm_Lx = norm(Lx);
      Lv = Lx./(s2 + lambda);
      step = (norm_Lx - alpha)*norm_Lx^2/((Lv'*Lx)*alpha);
      lambda = lambda + step;
    end

    % Terminate with an error if too many iterations.
    if (abs(step) > thr*lambda), error('Max. number of iterations reached'), end

    lambda = sqrt(lambda);

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