icamf.m

Independent Component Analysis源代码,最大似然方法

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); 
    end
end
fprintf(' %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 parameters
switch 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) % em

EM_step=0;

if par.draw
    h = waitbar(0,'Mean Field ICA - running EM optimization...') ;
end

while 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
    end

end % EM_step

if par.draw close(h), end


function [A,Sigma] = icaMFaem(A,Sigma,X,par) % aem

EM_step=0; 
fold = inf ; 
alpha = 1.05 ; % increase factor for overrelaxed
eta = 1 / alpha ;

if par.draw
    h = waitbar(0,'Mean Field ICA - running AEM optimization...') ;
end

while 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
    end

end % EM_step

if par.draw close(h), end