[Master Index] [Index for PublicToolbox/regutools]

# newton

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

### Function Synopsis

    lambda = newton(lambda_0,delta,s,beta,omega,delta_0)

### Help Text

    NEWTON Newton iteration (utility routine for DISCREP).

     lambda = newton(lambda_0,delta,s,beta,omega,delta_0)

     Uses Newton iteration to find the solution lambda to the equation
        || A x_lambda - b || = delta ,
     where x_lambda is the solution defined by Tikhonov regularization.

     The initial guess is lambda_0.

     The norm || A x_lambda - b || is computed via s, beta, omega and
     delta_0.  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.  Finally, delta_0 is the incompatibility measure.

### Cross-Reference Information

This function is called by

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

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

    function lambda = newton(lambda_0,delta,s,beta,omega,delta_0)
    %NEWTON Newton iteration (utility routine for DISCREP).
    %
    % lambda = newton(lambda_0,delta,s,beta,omega,delta_0)
    %
    % Uses Newton iteration to find the solution lambda to the equation
    %    || A x_lambda - b || = delta ,
    % where x_lambda is the solution defined by Tikhonov regularization.
    %
    % The initial guess is lambda_0.
    %
    % The norm || A x_lambda - b || is computed via s, beta, omega and
    % delta_0.  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.  Finally, delta_0 is the incompatibility measure.

    % Reference: V. A. Morozov, "Methods for Solving Incorrectly Posed
    % Problems", Springer, 1984; Chapter 26.

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

    % Set defaults.
    thr = 100*sqrt(eps);  % Relative stopping criterion.
    it_max = 50;      % Max number of iterations.
    % lambda_0 = eps;

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

    % Iterate; avoid negative values of lambda.
    lambda = lambda_0^2; step = 1; it = 0;
    while (abs(step) > thr*lambda & it < it_max), it = it+1;
      f = s2./(s2 + lambda);
      if (ps==1)
        r  = (1-f).*(beta - s.*omega);
        dr = f.*(beta - omega)./(s2 + lambda);
      else
        r  = (1-f).*(beta - sigma.*omega);
        dr = f.*(beta - sigma.*omega)./(s2 + lambda);
      end
      res = sqrt(r'*r + delta_0^2);
      step = -(res - delta)*res/(dr'*r);
      lambda = lambda + step;
      if (lambda <= 0), lambda = eps*max(s2); end
    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