function [S,A,U,ll,info]=icaML(X,K,par,debug_draw)
% icaML : ICA by ML (Infomax) with square mixing matrix and no noise.
%
% function [S,A,U,ll,info]=icaML(X,[K],[par]) Independent component analysis (ICA) using
% maximum likelihood, square mixing matrix and
% no noise [1] (Infomax). Source prior is assumed
% to be p(s)=1/pi*exp(-ln(cosh(s))). For optimization
% the BFGS algorithm is used [2]. See code for
% references.
%
% X : Mixed signals
% K : Number of source components.
% For K=0 (default) number of sources are equal to number of
% observations.
% For K < number of observations, SVD is used to reduce the
% dimension.
% par: Vector with 4 elements:
% (1) : Expected length of initial step
% Stopping criteria:
% (2) : Gradient ||g||_inf <= par(2)
% (3) : Parameter changes ||dW||_2 <= par(3)*(par(3) + ||W||_2)
% (4) : Maximum number of iterations
% Any illegal element in opts is replaced by its
% default value, [1 1e-4*||g(x0)||_inf 1e-8 100]
% debug_draw : Draw debug information
%
% S : Estimated source signals with variance
% scaled to one.
% A : Estimated mixing matrix
% U : Principal directions of preprocessing PCA.
% If K (the number of sources) is equal to the number
% of observations then no PCA is performed and U=eye(K).
% ll : Log likelihood for estimated sources
% info : Performance information, vector with 6 elements:
% (1:3) : final values of [ll ||g||_inf ||dx||_2]
% (4:5) : no. of iteration steps and evaluations of (ll,g)
% (6) : 1 means stopped by small gradient
% 2 means stopped by small x-step
% 3 means stopped by max number of iterations.
% 4 means stopped by zero step.
%
%
% Eks. Separate with number of sources equal number of
% observations.
%
% [S,A] = icaML(X);
%
% Eks. Separate with number of sources equal to k, using SVD
% as pre-processing.
%
% [S,A,U] = icaML(X,k);
%
% - version 1.5 (Revised 9/9-2003)
% - IMM, Technical University of Denmark
% - version 1.4
% - by Thomas Kolenda 2002 - IMM, Technical University of Denmark
% Revised: 9/9-2003, Mads, mad@imm.dtu.dk
% * Fixed help message to inform the user about U.
% * Removed the automatic use of PCA in the quadratic case.
% Bibtex references:
% [1]
% @article{Bell95,
% author = "A. Bell and T.J. Sejnowski",
% title = "An Information-Maximization Approach to Blind Separation and Blind Deconvolution",
% journal = "Neural Computation",
% year = "1995",
% volume = "7",
% pages = "1129-1159",
% }
%
% [2]
% @techreport{Nielsen01:unopt,
% author = "H.B. Nielsen",
% title = "UCMINF - an Algorithm for Unconstrained, Nonlinear Optimization ",
% institution = "IMM, Technical University of Denmark",
% number = "IMM-TEC-0019",
% year = "2001",
% url = "http://www.imm.dtu.dk/pubdb/views/edoc_download.php/642/ps/imm642.ps",
% }
% algorithm parameter settings
MaxNrIte = 1000;
try
debug = debug_draw;
catch
debug = 0;
end
if debug==1,fprintf('\n** Start icaML ***************************************\n');tic;end
% Scale X to avoid numerical problems
Xorg = X;
scaleX = max([abs(max(X(:))),abs(min(X(:)))]);
X = X./scaleX;
% set number of source parameters
if nargin<2, K=size(X,1); end
if ((K>0) & (K<size(X,1))),
[U,X] = callSVD(X,K,debug);
else
U=eye(size(X,1));
if debug==1,disp('Don''t use SVD');end
end
% initialize optimize parameters
try
ucminf_opt(1) = par(1);
ucminf_opt(2) = par(2);
ucminf_opt(3) = par(3);
ucminf_opt(4) = par(4);
catch
if debug==1,disp('Use default parameters to optimize');end
ucminf_opt = [1 1e-4 1e-8 MaxNrIte];
end
% initialize variables
[M,N] = size(X);
W = eye(M);
if debug==1,fprintf('Number of samples %d - Number of sources %d\n',N,M);end
% optimize
if debug==1,fprintf('Optimize ICA ... ');end
par.X=X; par.M=M; par.N=N;
[W,info] = ucminf( 'ica_MLf' , par , W(:) , ucminf_opt );
W = reshape(W,M,M);
% estimates
A=pinv(W);
S=W*X;
if debug==1,disp('done optimize ICA!');end
% sort components according to energy
Avar=diag(A'*A)/M;
Svar=diag(S*S')/N;
vS=var(S');
sig=Avar.*Svar;
[a,indx]=sort(sig);
S=S(indx(M:-1:1),:);
A=A(:,indx(M:-1:1));
% scale back
A=A.*scaleX;
% log likelihood
if nargout>3
ll= N*log(abs(det(inv(A)))) - sum(sum(log( cosh(S) ))) - N*M*log(pi);
end
if debug==1,fprintf('** End of icaML ****** time %0.2f sec******************\n\n',toc/60);end
function [f,dW]=ica_MLf(W,par)
% returns the negative log likelihood and its gradient w.r.t. W
X=par.X; M=par.M; N=par.N;
W=reshape(W,M,M);
S=W*X;
% negative log likelihood function
f=-( N*log(abs(det(W))) - sum(sum(log( cosh(S) ))) - N*M*log(pi) );
if nargout>1
% gradient w.r.t. W
dW=-(N*inv(W') - tanh(S)*X');
dW=dW(:);
end
function [U,DV] = callSVD(X,K,draw)
% Reduce dimension with SVD
[M,N] = size(X);
if N>M % Transpose if matrix is flat, to speed up the svd and later transpose back again
if draw==1,disp('Do Transpose SVD');end
[V,D,U] = svd(X',0);
else;
if draw==1,disp('Do SVD');end
[U,D,V] = svd(X,0);
end;
DV = D(1:K,1:K)*V(:,1:K)';
function [X, info, perf, D] = ucminf(fun,par, x0, opts, D0)
%UCMINF BFGS method for unconstrained nonlinear optimization:
% Find xm = argmin{f(x)} , where x is an n-vector and the scalar
% function F with gradient g (with elements g(i) = DF/Dx_i )
% must be given by a MATLAB function with with declaration
% function [F, g] = fun(x, par)
% par holds parameters of the function. It may be dummy.
%
% Call: [X, info {, perf {, D}}] = ucminf(fun,par, x0, opts {,D0})
%
% Input parameters
% fun : String with the name of the function.
% par : Parameters of the func
- 1
- 2
- 3
- 4
前往页