mcholmz3.m

该代码为盲源分离中的相对牛顿法

%
%  [L,D,E,pneg]=mchol2(G)
%
%  Given a symmetric matrix G, find a matrix E of "small" norm and c
%  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])
%  Modification:  Michael Zibulevsky   ([e-mail])
%
function [L,D,E,pneg]=mcholmz3(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]);


C=diag(diag(G));

L=zeros(n);
D=zeros(n,1);  %  use vestors as diagonal matrices
E=zeros(n,1);  %%****************** **************


%
%  Loop through, calculating column j of L for j=1:n
%


for j=1:n,

    bb=[1:j-1];
    ee=[j+1:n];


    if (j > 1),
        L(j,bb)=C(j,bb)./D(bb)';     %  Calculate the jth row of L.  
    end;


    
    if (j >= 2)
        if (j < n), 
            C(ee,j)=G(ee,j)- C(ee,bb)*L(j,bb)';    %  Update the jth column of C.
        end;
    else
        C(ee,j)=G(ee,j);
    end;
    
    
    %%%%     Update theta. 
    
    if (j == n)
        theta(j)=0;
    else
        theta(j)=max(abs(C(ee,j)));
    end;
    
    %%%%%  Update D 
    
    D(j)=max([eps,abs(C(j,j)),theta(j)^2/beta2]');
    
    
    %%%%%%  Update E. 
    
    E(j)=D(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;
    
    ind=[j*(n+1)+1 : n+1 : n*n]';
    C(ind)=C(ind)-(1/D(j))*C(ee,j).^2;


end;

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

ind=[1 : n+1 : n*n]';
L(ind)=1;


%
%  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;


return


%%%%%%%%%%%%%%%%%%%%%% Test %%%%%%%%%%%%%

n=3; G=rand(n);G=G+G';eigG=eig(G), 
[L,D,E,pneg]=mchol(G)
[L1,D1,E1,pneg1]=mcholmz1(G);
LL=L-L1, DD=D-D1, EE=E-E1, pp=pneg-pneg1

n=3; G=rand(n);G=G+G';eigG=eig(G),
[L,D,E,pneg]=mchol(G)
[L2,D2,E2,pneg2]=mcholmz2(G);
LL=L-L2, DD=D-diag(D2), EE=E-diag(E2), pp=pneg-pneg2