cfa.m

一个很好的Matlab编制的数据降维处理软件

function mappedX = cfa(X, no_dims, no_analyzers, no_iterations)
%CFA Performs manifold charting on dataset X 
%
%   mappedX = cfa(X, no_dims, no_analyzers, no_iterations)
%
% Performs manifold charting on dataset X to reduce its dimensionality to
% no_dims dimensions. The variable no_analyzers determines the number of
% local factor analyzers that is used in the mixture of factor analyzers
% (default = 40). The variable no_iterations sets the number of
% iterations that is employed in the EM algorithm (default = 200).
%

% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.4b.
% The toolbox can be obtained from http://www.cs.unimaas.nl/l.vandermaaten
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you 
% maintain the name of the original author.
%
% (C) Laurens van der Maaten
% Maastricht University, 2007

    if ~exist('no_dims', 'var')
        no_dims = 2;
    end
    if ~exist('no_analyzers', 'var')
        no_analyzers = 40;
    end
    if ~exist('no_iterations', 'var')
        no_iterations = 200;
    end
    
    % Make sure data is zero-mean, unit-variance
    X = X - repmat(mean(X, 1), [size(X, 1) 1]);
    X = X ./ repmat(var(X, 1), [size(X, 1) 1]);
    
    % Initialize some parameters
    min_var = 0;                     % minimum STD of Gaussians
    X = X';
    [D n] = size(X);
    c = no_analyzers;
    d = no_dims;
   
    % Randomly initialize the parameters of the CFA model
    SigmaN = rand(d, d, n);
    Z = rand(d, n);
    Pi = repmat(1 / c, [1 c]);
    Kappa = rand(d, c);
    Mu = rand(D, c);
    SigmaC = rand(d, d, c);
    Lambda = rand(D, d, c);
    Psi = zeros(D, D, c);
    for j=1:c
        tmp = zeros(D, D);
        tmp(1:size(tmp, 2) + 1:end) = rand(D, 1);
        Psi(:,:,j) = tmp;
    end
    
    % Perform the EM algorithm for the optimization
    iter = 0;
    while iter < no_iterations

        % E-step
        % ====================================================
        
        % Do some precomputations for speed
        invPsi = zeros(size(Psi));
        detSigmaC = zeros(1, c);
        detPsi = zeros(1, c);
        logPi = zeros(1, c);
        for j=1:c
            invPsi(:,:,j) = inv(Psi(:,:,j));
            detSigmaC(j) = det(SigmaC(:,:,j));
            detPsi(j) = det(Psi(:,:,j));
            logPi(j) = log(Pi(j));
        end
        
        % Compute matrices Epsilon, Vc, and m
        Eps = zeros(n, c);
        Vc = zeros(d, d, c);
        m = zeros(d, c);
        const = ((D + d) / 2) * log(2 * pi);
        for j=1:c
            
            % Precomputations
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            tmpProd = Lambda(:,:,j)' * invPsi(:,:,j) * Lambda(:,:,j);
            
            % Compute Epsilon
            for i=1:n            
                normX = (Xnc(:,i) - Lambda(:,:,j) * Znc(:,i));
                Eps(i, j) = -logPi(j) + const ...
                            + (.5 * log(detSigmaC(j))) + (.5 * detPsi(j)) ...
                            + (.5 * trace(SigmaC(:,:,j) * (SigmaN(:,:,i) + Znc(:,i) * Znc(:,i)'))) ...
                            + (.5 * trace(SigmaN(:,:,i) * tmpProd)) ...
                            + (.5 * normX' * invPsi(:,:,j) * normX);
            end
            
            % Compute Vc and m
            Vc(:,:,j) = inv(SigmaC(:,:,j)) + tmpProd;
            m(:,j) = Kappa(:,j) + inv(Vc(:,:,j)) * Lambda(:,:,j)' * invPsi(:,:,j) * mean(Xnc, 2);
        end
        
        % Update estimate of Q
        Q = (repmat(sum(exp(-Eps), 2), [1 c]) .^ -1) .* exp(-Eps);
        
        % Update estimate of SigmaN
        for i=1:n
            tmp = zeros(d, d);
            for j=1:c
                tmp = tmp + Q(i, j) * Vc(:,:,j);
            end
            
            % Enforce minimum variance
            tmp(1:size(tmp, 2) + 1:end) = max(min_var, tmp(1:size(tmp, 2) + 1:end));
            SigmaN(:,:,i) = inv(tmp);           % code above gave us inv(SigmaN)
        end
        
        % Update estimate of Z (= mappedX)
        for i=1:n
            tmp = zeros(d, 1);
            for j=1:c
                tmp = tmp + (Q(i, j) * Vc(:,:,j) * m(:,j));
            end
            Z(:,i) = SigmaN(:,:,i) * tmp;
        end
        
        % M-step
        % ====================================================
        
        % Update estimate of Pi
        Pi = sum(Q, 1) ./ n;
        
        % Update estimate of Mu and Kappa
        for j=1:c
            tmpQ = Q(:,c) ./ sum(Q(:,c));
            Mu(:,c) = sum(repmat(tmpQ', [D 1]) .* X, 2);
            Kappa(:,c) = sum(repmat(tmpQ', [d 1]) .* Z, 2);
        end
        
        % Update estimate of SigmaC
        for j=1:c
            tmp = 0;
            tmpQ = Q(:,c) ./ sum(Q(:,c));
            for i=1:n
                Znc = Z(:,i) - Kappa(:,j);
                tmp = tmp + (tmpQ(i) * (SigmaN(:,:,i) + Znc * Znc'));
            end
            
            % Enforce some variance
            tmp(1:size(tmp, 2) + 1:end) = max(min_var, tmp(1:size(tmp, 2) + 1:end));
            SigmaC(:,:,j) = tmp;
        end
        
        % Update estimate of Lambda
        for j=1:c
            Sc = zeros(D, d);
            tmpQ = Q(:,j) ./ sum(Q(:,j));
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            for i=1:n
                Sc = Sc + tmpQ(i) * (Xnc(:,i) * Znc(:,i)');
            end
            Lambda(:,:,j) = Sc * inv(SigmaC(:,:,j));
        end
        
        % Update estimate of Psi
        for j=1:c
            tmpPsi = zeros(D, D);
            tmpQ = Q(:,j) ./ sum(Q(:,j));
            Xnc = X - repmat(Mu(:,j), [1 n]);
            Znc = Z - repmat(Kappa(:,j), [1 n]);
            tmpProd = Lambda(:,:,j) * SigmaN(:,:,i) * Lambda(:,:,j)';
            for i=1:n
                tmpPsi(1:size(tmpPsi, 2) + 1:end) = tmpPsi(1:size(tmpPsi, 2) + 1:end) + ...
                    tmpQ(i) * (((Xnc(:,i) - Lambda(:,:,j) * Znc(:,i)) .^ 2)' + tmpProd(1:size(tmpProd, 2) + 1:end));
            end
            Psi(:,:,c) = tmpPsi;
        end
        
        % Update number of iterations
        iter = iter + 1;
        if rem(iter, 5) == 0
            fprintf('.');
        end
    end
    
    % Transpose to get lowdimensional data representation
    mappedX = Z';