pca_full.m

Matlab package for PCA for datasets with missing values

                        M, sXv, ndata );
        %cost = cf_full( X, A, S, Mu, V, Va, Vmu, M, sXv, ndata );
        %cost = cf_ppca( X, M, A.me, S.me, Mu, V, S.va, S.Iv );
        lc.cost = [ lc.cost cost ];
    end
    
    DisplayProgress( dsph, lc )
    angleA = subspace(A,Aold);
    PrintStep( opts.verbose, lc, angleA )

    convmsg = converg_check( opts, lc, angleA );
    if ~isempty(convmsg)
        if opts.verbose, fprintf( '%s', convmsg ), end
        break
    end
    Aold = A;
    
    time = clock;
    if etime(time,time_autosave) > opts.autosave
        time_autosave = time;
        save( opts.filename,...
              'A', 'S', 'Mu', 'V', 'Av', 'Muv', 'Sv', 'Isv',...
              'Va', 'Vmu',...
              'lc', 'Ir', 'Ic', 'n1x', 'n2x', 'n1', 'n2' )
    end
end

%[ A, S, Sv ] = RotateToPCA( A, S, Sv );

if n1 < n1x
    [ A, Av ] = addmrows( A, Av, Ir, n1x, Va );
    [ Mu, Muv ] = addmrows( Mu, Muv, Ir, n1x, Vmu );
    %A = addmrows( A, Ir, n1x, NaN );
    %Mu = addmrows( Mu, Ir, n1x, NaN );
end
if n2 < n2x
    [ S, Sv, Isv ] = addmcols( S, Sv, Ic, n2x, Isv );
end

if nargout == 1
    A = struct( 'A', A );
    A.S = S;
    A.Mu = Mu;
    A.V = V;
    A.Va = Va;
    A.Vmu = Vmu;
    A.Av = Av;
    A.Sv = Sv;
    A.Isv = Isv;
    A.Muv = Muv;
    A.lc = lc;

else
    cv.A = Av;
    cv.S = Sv;
    cv.Isv = Isv;
    cv.Mu = Muv;

    hp.Va = Va;
    hp.Vmu = Vmu;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X,Xprobe] = SubtractMu( Mu, X, M, Xprobe, Mprobe, update_bias )

n2 = size(X,2);

if ~update_bias
    return
end
   
if issparse(X)
    X = subtract_mu( X, Mu );
    if ~isempty(Xprobe)
        Xprobe = subtract_mu( Xprobe, Mu );
    end
else
    X = X - repmat(Mu,1,n2).*M;
    if ~isempty(Xprobe)
        Xprobe = Xprobe - repmat( Mu, 1, n2 ).*Mprobe;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the PCA rotation: This has to be checked
function [ dMu, A, Av, S, Sv ] = ...
    RotateToPCA( A, Av, S, Sv, Isv, obscombj, update_bias );

n1 = size(A,1);
n2 = size(S,2);

% TODO: Take into account the prior for Mu, A???
if update_bias
    mS = mean(S,2);
    dMu = A*mS;
    S = S - repmat(mS,1,n2);
else
    dMu = 0;
end

covS = S*S';
if isempty(Isv)
    for j = 1:n2
        covS = covS + Sv{j};
    end
else
    nobscomb = length(obscombj);
    for j = 1:nobscomb
        covS = covS + ( length(obscombj{j})*Sv{j} );
    end
end
    
covS = covS / n2;
[VS,D] = eig(covS);
RA = VS*sqrt(D);
A = A*RA;
covA = A'*A;
if ~isempty(Av)
    for i = 1:n1
        Av{i} = RA'*Av{i}*RA;
        covA = covA + Av{i};
    end
end
covA = covA / n1;
[VA,DA] = eig(covA);
[DA,I] = sort( -diag(DA) );
DA = -DA;
VA = VA(:,I);
A = A*VA;

if ~isempty(Av)
    for i = 1:n1
        Av{i} = VA'*Av{i}*VA;
    end
end
R = VA'*diag(1./sqrt(diag(D)))*VS';

S = R*S;
for j = 1:length(Sv)
    Sv{j} = R*Sv{j}*R';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ A, S, Mu, V, Av, Sv, Muv ] = ...
    InitParms( init, n1, n2, ncomp, nobscomb, Isv )

if ischar(init)
    if strcmpi( init, 'random' )
        init = struct([]);
    else
        % Load from a file
        init = load( init );
    end
end

if isstruct(init)
    if isfield( init, 'A' )
        A = init.A;
    else
        A = orth(randn(n1,ncomp));
    end
    if isfield( init, 'Av' ) && ~isempty(init.Av)
        if iscell(init.Av)
            Av = init.Av;
        else
            for i = 1:n1, Av{i} = diag(init.Av(i,:)); end
        end
    else
        Av = cell(1,n1);
        for i = 1:n1
            Av{i} = eye(ncomp);
        end
    end

    if isfield( init, 'Mu' )
        Mu = init.Mu;
    else
        Mu = [];
    end
    if isfield( init, 'Muv' ) && ~isempty(init.Muv)
        Muv = init.Muv;
    else
        Muv = ones(n1,1);
    end

    if isfield( init, 'V' )
        V = init.V;
    else
        V = 1;
    end
    
    if isfield( init, 'S' )
        S = init.S;
    else
        S = randn(ncomp,n2);
    end

    % TODO: Check this
    %Sv = {};
    if isfield( init, 'Sv' ) && ~isempty(init.Sv)
       if nobscomb < n2
            [B,I] = unique(Isv,'first');
            if ~iscell(init.Sv)
                Sv = cell(1,nobscomb);
                for j = 1:nobscomb
                    Sv{j} = diag(init.Sv(:,Isv(I(j))));
                end
            
            elseif isfield( init, 'Isv' ) && ~isempty(init.Isv)
                Sv = { init.Sv{ init.Isv(I) } };
            else
                for j = 1:nobscomb
                    Sv{j} = init.Sv{Isv(I(j))};
                end
            end
        else
            if ~iscell(init.Sv)
                Sv = cell(1,n2);
                for j = 1:n2, Sv{j} = diag(init.Sv(:,j)); end
            
            elseif isfield( init, 'Isv' ) && ~isempty(init.Isv)
                Sv = { init.Sv{ init.Isv } };
            elseif length(init.Sv) == n2
                Sv = init.Sv;
            end
       end
    else
        Sv = cell(1,nobscomb);
        for j = 1:nobscomb, Sv{j} = eye(ncomp); end
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PrintFirstStep( verbose, rms, prms )

if ~verbose
    return
end
fprintf( 'Step 0: rms = %.6f', rms )
if ~isnan(prms)
    fprintf( ' (%.6f)', prms )
end
fprintf( '\n' )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PrintStep( verbose, lc, Aangle )

if ~verbose
    return
end

iter = length(lc.rms)-1;
steptime = lc.time(end)-lc.time(end-1);

fprintf( 'Step %i: ', iter )
if ~isnan(lc.cost(end))
    fprintf( 'cost = %.6f, ', lc.cost(end) );
end
fprintf( 'rms = %.6f', lc.rms(end) );
if ~isnan(lc.prms(end))
    fprintf( ' (%.6f)', lc.prms(end) );
end
fprintf( ', angle = %.2e', Aangle )
if steptime > 1
    fprintf( ' (%i sec)\n', round(steptime) )
else
    fprintf( ' (%.0e sec)\n', steptime )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PrintProgressBar( verbose, str )

if verbose == 2
    fprintf( [ str '\n' ] )
    %fprintf( '\n|                                                  |\r|' )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function PrintProgress( verbose, i, n, str )

if verbose == 2
    fprintf( '\r%s %i/%i', str, i, n )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dsph = DisplayInit( display, lc )

dsph.display = display;
if ~dsph.display
    return
end
dsph.fig = figure;
subplot(2,1,1)
dsph.rms = plot( 0:length(lc.rms)-1, lc.rms );
title( 'RMS training error' )
subplot(2,1,2)
title( 'RMS test error' )
dsph.prms = plot( 0:length(lc.prms)-1, lc.prms );

drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DisplayProgress( dsph, lc )

if dsph.display
    set( dsph.rms, 'xdata', 0:length(lc.rms)-1,...
                   'ydata', lc.rms )
    set( dsph.prms, 'xdata', 0:length(lc.prms)-1,...
                   'ydata', lc.prms )
    drawnow
end