function [S,A,loglikelihood,Sigma,chi]=icaMF(X,par);% ICAMF  Mean field independent component analysis (ICA)%    [S,A,LOGLIKELIHOOD,SIGMA,CHI]=ICAMF(X) performs linear %    instanteneous mixing ICA with additive noise X = A S + eta.  %    Empirical Bayes (maximum Likelihood) is used to estimate the mixing %    matrix A and the noise covariance SIGMA [1]. The sufficient statistics %    (mean and covarinces of the sources) are estimated using either %    variational with linear response correction for covariances [1] or%    with the expectation consistent framework (formerly known as adaptive%    TAP) [1-4].%   %    Input and output arguments: %      %    X        : Mixed signal%    A        : Estimated mixing matrix%    S        : Estimated source mean values                                   %    SIGMA    : Estimated noise covariance%    LL       : Normalized marginal log likelihood, log P(X|A,SIGMA) / N %    CHI      : Estimated source covariances%%    Outputs are sorted in descending order according to 'energy'%    sum(A.^2,1))'.*sum(S.^2,2).%%    [S,A,LL,SIGMA,CHI]=ICA_ADATAP(X,PAR) allows for specification of %    parameters to specify priors on A, S and SIGMA, the method to be used%    and additional parameters. PAR.method overrides other choices of %    priors. PAR.Aprior, PAR.Sprior, PAR.Sigmaprior and PAR.method are %    character strings that can take the following values%    %    PAR.Aprior :         constant      constant mixing matrix %                         free          unconstrained optimization%                         positive      positively constrained optimization %                                       (uses QUADPROG if negative data or%                                       sources occur)%%    PAR.Sprior :         bigauss       sum of two Gaussians, default%                                       variance 1 centered at +/-1.  %                         binary        binary +/-1 %                         binary_01     binary 0/1             %                         combi         combinations of priors%                                       E.g. 2 binary_01's and 1%                                       exponential are specified by%                                       PAR.S(1).prior='binary_01';%                                       PAR.S(2).prior='binary_01';%                                       PAR.S(3).prior='exponential'; %                                       Parameters to the priors can also%                                       be specified, see mog prior below.%                         exponential   exponential (positive)%                         Gauss         Gaussian for factor analysis and %                                       probabilistic PCA, see below.%                         heavy_tail    heavy_tailed (non-analytic %  				                        power law tail)%                         Laplace       Laplace (double exponential)%                         mog           mixture of Gaussians. Parameters %                                       for mixing proportions, mean and %                                       variances are specified with %                                       PAR.S.pi, PAR.S.mu and PAR.S.Sigma. %                                       Parameters can also be specific to%                                       the source using PAR.S(source).pi%                                       etc. Default is shared and heavy-%                                       tailed: PAR.S.pi = [ 0.5 0.5 ]', %                                       PAR.S.mu = [ 0 0 ]' and%                                       PAR.S.Sigma = [ 1 0.01 ]'.%%    PAR.Sigmaprior :     constant      constant noise covariance%                         free          full noise covariance matrix  %                         isotropic     white noise, with same (scalar) %                                       noise variance on all sensors%                         diagonal      white noise, different noise %                                       variance for each sensor.%%%    PRIOR.method :       constant      constant A and Sigma %                                       for test sets. A and Sigma%                                       should initialized, see %                                       below. Sets PAR.Aprior='constant'%                                       and PAR.Sigmaprior='constant'.%                         fa            factor analysis (FA) sets%                                       PAR.Aprior='free', %                                       PAR.Sprior='Gauss' and%                                       PAR.Sigmaprior='diagonal'.%                         neg_kurtosis  negative kurtosis set%                                       PAR.Sprior='bigauss'.%                         pos_kurtosis  positive kurtosis sets%                                       PAR.Sprior='mog'. %                         positive      positive ICA sets%                                       PAR.Sprior='exponential' and%                                       PAR.Aprior='positive'.%                         ppca          probabilistic PCA (pPCA)%                                       sets PAR.Aprior='free', %                                       PAR.Sprior='Gauss' and %                                       PAR.Sigmaprior='isotropic'.%                         default       PAR.Aprior='free', %                                       PAR.Sprior='mog', and %                                       PAR.Sigmaprior='isotropic'%%    Additional fields that can be set in PAR (e.g. PAR.sources) are%   %       sources           Number of sources (default is quadratic %                         size(X,1)). This parameter is overridden %                         by size(PAR.A_init,2) if PAR.A_init is defined.%       optimizer         optimization method for A and SIGMA%%                         aem           adaptive overrelaxed em by%                                       Salakhutdinov and Roweis, default.%                         em            expectation maximization%                         bfgs          bfgs quasi-Newton with line search.%                         conjgrad      conjugated gradient.%%       solver            method for estimating the statistics of the %                         sources%%                         ec            expectation consistent framework,%                                       default.%                         variational   variational (Bayes) with linear %                                       response correction for source%                                       covariances and likelihood %                                       ( = - free energy).%                         %       A_init            Initial mixing matrix, default is small random %                         values.        %       S_init            Initial source values, default is zero only used%                         for estimating initial noise covariance.%       Sigma_init        Initial noise covariance (scalar for isotropic). %                         Default is empirical variance.%       max_ite           Number of A and SIGMA updates, default is 50.%       S_max_ite         Maximum number of source mean iteration steps, %                         default 100.%       S_tol             Required squared deviation of solution of mean%                         field equations, default 1e-10.%       draw              toggle whether to plot run-time information,%                         defualt is 1.%%%    The memory requirements are O(M^2*N+D*N), where D is the %    number of sensors, N the number of samples (X is D*N) and M is %    the number of sources. The computational complexity is typically%    dominated by estimation of sources covariances O(M^3*N). Additional %    complexity can come when D is large for positive mixing matrix and/or%    full noise covariances.  Use of QUADPROG (see above) also tends to %    slow down the program.%%    It is recommended as a preprocessing step to normalize the %    data, e.g. by dividing by the largest absolute element of the data: %    X=X/max(max(abs(X))).%%    See also icaMF_bic%% References [1-4]:%%    Mean Field Approaches to Independent Component Analysis%    Pedro A.d.F.R. H鴍en-S鴕ensen, Ole Winther and Lars Kai Hansen%    Neural Computation 14, 889-918 (2002).%%    Expectation Consistent Approximate Inference%    Manfred Opper and Ole Winther%    Journal of Machine Learning Research (2005).%  %    Tractable Approximations for Probabilistic Models: The Adaptive %    Thouless-Anderson-Palmer Mean Field Approach%    M. Opper and O. Winther%    Phys. Rev. Lett. 86, 3695-3699 (2001).%%    Adaptive and Self-averaging Thouless-Anderson-Palmer Mean Field %    Theory for Probabilistic Modeling%    M. Opper and O. Winther%    Physical Review E 64, 056131 (2001). % % Program uses notation of reference [2].% - by Ole Winther 2002, 2005 - IMM, Technical University of Denmark% - http://isp.imm.dtu.dk/staff/winther/% - version 3.0[A,Sigma,par] = initializeMF(X,par) ;fprintf(' Mean Field ICA empirical Bayes algorithm\n')fprintf(' %s mean field solver, %s source prior\n',par.solver,par.Sprior);  if strcmp(par.Sprior,'combi')    for i=1:par.sources        fprintf('   %s source prior\n',par.S(i).prior);     endendfprintf(' %s parameter optimizer for %s A and %s noise\n',...    par.optimizer,par.Aprior,par.Sigmaprior);fprintf(' sensors: D = %i\n sources: M = %i\n samples: N = %i\n\n',par.D,par.M,par.N); if par.draw   if size(Sigma,1) * size(Sigma,2) == par.D    fprintf(' initial noise variance = %g \n\n',sum(Sigma)/par.D);  else           fprintf(' initial noise variance = %g \n\n',trace(Sigma)/size(Sigma,1));  end  end% run mean field ICA empirical Bayes algorithm% choose optimization method for parametersswitch lower(par.optimizer)    case 'constant'        case 'conjgrad'                 if par.draw            h = waitbar(1,'Mean Field ICA - conjugated gradient optimization...') ;        end                popt = ASigma2popt(A,Sigma,par) ;                popt = minimize(popt,'icaMFconjgrad',par.max_ite,X,par) ; % conjugated gradients            [A,Sigma] = popt2ASigma(popt,par) ;                if par.draw close(h), end    case 'bfgs'         if par.draw            h = waitbar(1,'Mean Field ICA - BFGS optimization...') ;        end                popt = ASigma2popt(A,Sigma,par) ;              par.X = X ;        ucminf_opt = [10  1e-8  1e-12  par.max_ite];        [popt,info] = ucminf( 'icaMFbfgs' , par , popt , ucminf_opt ); % BFGS           [A,Sigma] = popt2ASigma(popt,par) ;        if par.draw close(h), end            case 'em'            [A,Sigma] = icaMFem(A,Sigma,X,par) ;        case 'aem' % adaptive overrelaxed em by R Salakhutdinov + S Roweis        [A,Sigma] = icaMFaem(A,Sigma,X,par) ;        end% solve mean field equations and exit![theta,J,par] = XASigma2thetaJ(X,A,Sigma,par) ;[S,chi,f,ite,dSdS] = par.mf_solver(theta,J,par) ; loglikelihood = - f  ;% return variables sorted according to 'energy'  [Y,I]=sort((sum(A.^2,1))'.*sum(S.^2,2)) ;I=flipud(I);S=S(I,:); A=A(:,I); for i=1:par.N  chi(:,:,i)=squeeze(chi(I,I,i));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                        end of main program                                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                          help functions                                           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              free energy (minus log likelihood) and derivative wrt parameters                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [f,df] = icaMFconjgrad(popt,X,par) % for minimize[A,Sigma] = popt2ASigma(popt,par);[theta,J,par] = XASigma2thetaJ(X,A,Sigma,par) ;[S,chi,f] = par.mf_solver(theta,J,par) ;df = dpoptf(X,A,Sigma,S,chi,par) ;function [f,df] = icaMFbfgs(popt,par) % for ucminf[A,Sigma] = popt2ASigma(popt,par);[theta,J,par] = XASigma2thetaJ(par.X,A,Sigma,par) ;[S,chi,f] = par.mf_solver(theta,J,par) ;df = dpoptf(par.X,A,Sigma,S,chi,par) ;function [A,Sigma] = icaMFem(A,Sigma,X,par) % emEM_step=0;if par.draw    h = waitbar(0,'Mean Field ICA - running EM optimization...') ;endwhile EM_step<par.max_ite %& (dA_rel>tol | dSigma_rel>tol)     EM_step=EM_step+1;        if par.draw waitbar(EM_step/par.max_ite,h), end        [theta,J,par] = XASigma2thetaJ(X,A,Sigma,par) ;    if par.draw        [S,chi,f,S_ite,dSdS]=par.mf_solver(theta,J,par) ;    else        [S,chi] = par.mf_solver(theta,J,par) ;    end        [Anew,Sigmanew] = dASigmaf(X,A,Sigma,S,chi,par) ;    dpoptrel = sum(abs([ (Anew(:) - A(:))'  (Sigmanew(:) - Sigma(:))' ])) / ...               ( sum(abs([A(:)' Sigma(:)'])) + eps ) ;       % em update    A = Anew ;    Sigma = Sigmanew ;                   if par.draw        fprintf(' Iteration %d\n S_ite = %d   dSdS = %g dpar/|par| = %g\n',...            EM_step,S_ite,dSdS,dpoptrel);        fprintf(' loglikelihood = %6.3f  ',-f);        if size(Sigma,1)*size(Sigma,2)==par.D             fprintf(' noise variance = %g \n\n',sum(Sigma)/D);        else                   fprintf(' noise variance = %g \n\n',trace(Sigma)/size(Sigma,1));        end    endend % EM_stepif par.draw close(h), endfunction [A,Sigma] = icaMFaem(A,Sigma,X,par) % aemEM_step=0; fold = inf ; alpha = 1.05 ; % increase factor for overrelaxedeta = 1 / alpha ;if par.draw    h = waitbar(0,'Mean Field ICA - running AEM optimization...') ;endwhile EM_step<par.max_ite %& (dA_rel>tol | dSigma_rel>tol)     EM_step=EM_step+1;        if par.draw waitbar(EM_step/par.max_ite,h), end        [theta,J,par] = XASigma2thetaJ(X,A,Sigma,par) ;        [S,chi,f,S_ite,dSdS] = par.mf_solver(theta,J,par) ;        [Aem,Sigmaem] = dASigmaf(X,A,Sigma,S,chi,par) ;       if f > fold & eta > 1        % we took a step that decreased the likelihood        % backtrack and take regular em step instead        eta = 1;        A = Aold ; Sigma = Sigmaold ; f = fold ;        Aem = Aemold ; Sigmaem = Sigmaemold ;    else        % increase step size        eta = alpha * eta ;        Aemold = Aem ; Sigmaemold = Sigmaem ;    end       % aem update    Aold = A ; Sigmaold = Sigma ;  fold = f ;     switch par.Aprior        case 'free'            A = A + eta * ( Aem - A ) ;        case 'positive'            A = A .* ( Aem ./ A ).^eta ;            maxratio = 2; minratio = 0.5 ; % maximim increase/decrease factors            A = min(A,maxratio*Aold); A = max(A,minratio*Aold);    end           switch par.Sigmaprior         case 'free'            logSigma = logm(Sigma) ;             Sigma = expm( logSigma + eta * ( logm(Sigmaem) - logSigma) ) ;        case {'isotropic','diagonal'}            Sigma = Sigma .* ( Sigmaem ./ Sigma ).^eta ;    end      dpoptrel = sum(abs([ (A(:) - Aold(:))'  (Sigma(:) - Sigmaold(:))' ])) / ...                ( sum(abs([Aold(:)' Sigmaold(:)'])) + eps ) ;        if par.draw        fprintf(' Iteration %d\n S_ite = %d   dSdS = %g  eta  = %g  dpar/|par| = %g\n',...            EM_step,S_ite,dSdS,eta,dpoptrel);        fprintf(' loglikelihood = %6.3f  ',-f);        if size(Sigma,1)*size(Sigma,2)==par.D             fprintf(' noise variance = %g \n\n',sum(Sigma)/par.D);        else                   fprintf(' noise variance = %g \n\n',trace(Sigma)/size(Sigma,1));        end    endend % EM_stepif par.draw close(h), end