mchol.m

这是在网上下的一个东东

%
%  [L,D,E,pneg]=mchol(G)
%
%  Given a symmetric matrix G, find a matrix E of "small" norm and 
%  L, and D such that  G+E is Positive Definite, and 
%
%      G+E = L*D*L'
%
%  Also, calculate a direction pneg, such that if G is not PSD, then
%
%      pneg'*G*pneg < 0
%
%  Reference: Gill, Murray, and Wright, "Practical Optimization", p111.
%  Author: Brian Borchers ([e-mail])
%
function [L,D,E,pneg]=mchol(G)
%
%  n gives the size of the matrix.
%
n=size(G,1);
%
%  gamma, zi, nu, and beta2 are quantities used by the algorithm.  
%
gamma=max(diag(G));
zi=max(max(G-diag(diag(G))));
nu=max([1,sqrt(n^2-1)]);
beta2=max([gamma, zi/nu, 1.0E-15]);
%
%  Initialize diag(C) to diag(G).
%
C=diag(diag(G));
%
%  Loop through, calculating column j of L for j=1:n
%
for j=1:n,
%
%  Calculate the jth row of L.  
%
  if (j > 1),
    L(j,1:j-1)=C(j,1:j-1)./diag(D(1:j-1,1:j-1))';
  end;
%
%  Update the jth column of C.
%
 if (j >= 2),
   if (j < n), 
     C(j+1:n,j)=G(j+1:n,j)-(L(j,1:j-1)*C(j+1:n,1:j-1)')';
   end;
 else
   C(j+1:n,j)=G(j+1:n,j);
 end;
%
% Update theta. 
%
 if (j == n)
   theta(j)=0;
  else
   theta(j)=max(abs(C(j+1:n,j)));
  end;
%
%  Update D
%
 D(j,j)=max([eps,abs(C(j,j)),theta(j)^2/beta2]');
%
% Update E.
%
 E(j,j)=D(j,j)-C(j,j);
%
%  Update C again...
% 
 for i=j+1:n,
   C(i,i)=C(i,i)-C(i,j)^2/D(j,j);
 end;

end;
%
% Put 1's on the diagonal of L
%
for j=1:n,
  L(j,j)=1;
end;

%
%  if needed, finded a descent direction.  
%
if (nargout == 4)
  [m,col]=min(diag(C));
  rhs=zeros(n,1);
  rhs(col)=1;
  pneg=L'\rhs;
end;