ICA工具箱（Infomax实现）

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