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,w,x0) 
% 
% 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 insures 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, IMM 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 insure 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