[Master Index] [Index for PublicToolbox/regutools]

# maxent

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

### Function Synopsis

    [x_lambda,rho,eta,data,X] = maxent(A,b,lambda,w,x0)

### Help Text

    MAXENT Maximum entropy regularization.

     [x_lambda,rho,eta] = maxent(A,b,lambda)

     Maximum entropy regularization:
        min { || A x - b ||^2 + lambda^2*x'*log(diag(w)*x) } ,
     where -x'*log(diag(w)*x) is the entropy of the solution x.
     If no weights w are specified, unit weights are used.

     If lambda is a vector, then x_lambda is a matrix such that
        x_lambda = [x_lambda(1), x_lambda(2), ... ] .

     This routine uses a nonlinear conjugate gradient algorithm with "soft"
     line search and a step-length control that ensures a positive solution.
     If the starting vector x0 is not specified, then the default is
        x0 = norm(b)/norm(A,1)*ones(n,1) .

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

    function [x_lambda,rho,eta,data,X] = maxent(A,b,lambda,w,x0)
    %MAXENT Maximum entropy regularization.
    %
    % [x_lambda,rho,eta] = maxent(A,b,lambda)
    %
    % Maximum entropy regularization:
    %    min { || A x - b ||^2 + lambda^2*x'*log(diag(w)*x) } ,
    % where -x'*log(diag(w)*x) is the entropy of the solution x.
    % If no weights w are specified, unit weights are used.
    %
    % If lambda is a vector, then x_lambda is a matrix such that
    %    x_lambda = [x_lambda(1), x_lambda(2), ... ] .
    %
    % This routine uses a nonlinear conjugate gradient algorithm with "soft"
    % line search and a step-length control that ensures a positive solution.
    % If the starting vector x0 is not specified, then the default is
    %    x0 = norm(b)/norm(A,1)*ones(n,1) .

    % Per Christian Hansen, UNI-C and Tommy Elfving, Dept. of Mathematics,
    % Linkoping University, 06/10/92.

    % Reference: R. Fletcher, "Practical Methods for Optimization",
    % Second Edition, Wiley, Chichester, 1987.

    % Set defaults.
    flat = 1e-3;     % Measures a flat minimum.
    flatrange = 10;  % How many iterations before a minimum is considered flat.
    maxit = 150;     % Maximum number of CG iterations;
    minstep = 1e-12; % Determines the accuracy of x_lambda.
    sigma = 0.5;     % Threshold used in descent test.
    tau0 = 1e-3;     % Initial threshold used in secant root finder.

    % Initialization.
    [m,n] = size(A); x_lambda = zeros(n,length(lambda)); F = zeros(maxit,1);
    if (min(lambda) <= 0)
      error('Regularization parameter lambda must be positive')
    end
    if (nargin ==3), w  = ones(n,1); end
    if (nargin < 5), x0 = ones(n,1); end

    % Treat each lambda separately.
    for j=1:length(lambda);

      % Prepare for nonlinear CG iteration.
      l2 = lambda(j)^2;
      x  = x0; Ax = A*x;
      g  = 2*A'*(Ax - b) + l2*(1 + log(w.*x));
      p  = -g;
      r  = Ax - b;

      % Start the nonlinear CG iteration here.
      delta_x = x; dF = 1; it = 0; phi0 = p'*g;
      while (norm(delta_x) > minstep*norm(x) & dF > flat & it < maxit & phi0 < 0)
        it = it + 1;

        % Compute some CG quantities.
        Ap = A*p; gamma = Ap'*Ap; v = A'*Ap;

        % Determine the steplength alpha by "soft" line search in which
        % the minimum of phi(alpha) = p'*g(x + alpha*p) is determined to
        % a certain "soft" tolerance.
        % First compute initial parameters for the root finder.
        alpha_left = 0; phi_left = phi0;
        if (min(p) >= 0)
          alpha_right = -phi0/(2*gamma);
          h = 1 + alpha_right*p./x;
        else
          % Step-length control to ensure a positive x + alpha*p.
          I = find(p < 0);
          alpha_right = min(-x(I)./p(I));
          h = 1 + alpha_right*p./x; delta = eps;
          while (min(h) <= 0)
            alpha_right = alpha_right*(1 - delta);
            h = 1 + alpha_right*p./x;
            delta = delta*2;
          end
        end
        z = log(h);
        phi_right = phi0 + 2*alpha_right*gamma + l2*p'*z;
        alpha = alpha_right; phi = phi_right;

        if (phi_right <= 0)

          % Special treatment of the case when phi(alpha_right) = 0.
          z = log(1 + alpha*p./x);
          g_new = g + l2*z + 2*alpha*v; t = g_new'*g_new;
          beta = (t - g'*g_new)/(phi - phi0);

        else

          % The regular case: improve the steplength alpha iteratively
          % until the new step is a descent step.
          t = 1; u = 1; tau = tau0;
          while (u > -sigma*t)

            % Use the secant method to improve the root of phi(alpha) = 0
            % to within an accuracy determined by tau.
            while (abs(phi/phi0) > tau)
              alpha = (alpha_left*phi_right - alpha_right*phi_left)/...
                      (phi_right - phi_left);
              z = log(1 + alpha*p./x);
              phi = phi0 + 2*alpha*gamma + l2*p'*z;
              if (phi > 0)
                alpha_right = alpha; phi_right = phi;
              else
                alpha_left  = alpha; phi_left  = phi;
              end
            end

            % To check the descent step, compute u = p'*g_new and
            % t = norm(g_new)^2, where g_new is the gradient at x + alpha*p.
            g_new = g + l2*z + 2*alpha*v; t = g_new'*g_new;
            beta = (t - g'*g_new)/(phi - phi0);
            u = -t + beta*phi;
            tau = tau/10;

          end  % End of improvement iteration.

        end  % End of regular case.

        % Update the iteration vectors.
        g = g_new; delta_x = alpha*p;
        x = x + delta_x;
        p = -g + beta*p;
        r = r + alpha*Ap;
        phi0 = p'*g;

        % Compute some norms and check for flat minimum.
        rho(j,1) = norm(r); eta(j,1) = x'*log(w.*x);
        F(it) = rho(j,1)^2 + l2*eta(j,1);
        if (it <= flatrange)
          dF = 1;
        else
          dF = abs(F(it) - F(it-flatrange))/abs(F(it));
        end

        data(it,:) = [F(it),norm(delta_x),norm(g)];
        X(:,it) = x;

      end  % End of iteration for x_lambda(j).

      x_lambda(:,j) = x;

    end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON