fast_newt_dir.m

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

function X=fast_newt_dir(D,B),

%%%%%%%%%% Fast solution of Newton matrix equation  D .* X + X' = B %%%%%%%%%%
%          using eigenvalue decomposition of 2x2 matrices                         
delta=1e-8;

normB=norm(B(:)); %% Scale to avoid too small/large numbers
B=B/normB;

b1=B;
b2=B';

a11=D;
a22=D';

Tr=a11+a22;    % Trace
dt=a11.*a22-1; % Determinant

T2= Tr/2;
sq=sqrt(T2.^2-dt);
lam1=(T2-sq);
lam2=(T2+sq);

ss=lam1-a11;

norms=sqrt(1+ss.^2);

s11=1./norms;
s21=ss./norms;
s12=s21;
s22=-s11;

c1=(s11.*b1+s21.*b2)./max(abs(lam1),delta);
c2=(s12.*b1+s22.*b2)./max(abs(lam2),delta);
%c1=(s11.*b1+s21.*b2)./abs(lam1);
%c2=(s12.*b1+s22.*b2)./abs(lam2);
%c1=(s11.*b1+s21.*b2)./lam1;
%c2=(s12.*b1+s22.*b2)./lam2;

x1=s11.*c1+s12.*c2;
%x2=s21*c1+s22*c2

X=x1;

% Diagonal should be treated separately
[n,n]=size(D); idiag=[1:(n+1):n^2]';
X(idiag)=B(idiag)./(D(idiag)+1);

lam1(idiag)= D(idiag)+1;
lam2(idiag)=lam1(idiag);


lammin=min([lam1(:);lam2(:)]);
lammax=max([lam1(:);lam2(:)]);
fprintf(' eigval=%.2g _ %.2g     ',lammin,lammax)
%if lammin < delta, fprintf('\nNegative or small eig. value of the Hessian: lammin=%.2e\n',lammin);end
%keyboard

X=X*normB;  %% Scale back to avoid too small/large numbers