cca.m

一个基于matlab的数据降维工具箱,包括MDS,LEE等方法

function [Z, ccaEigen, ccaDetails] = cca(X, Y, EDGES, OPTS)
%
% Function [Z, CCAEIGEN, CCADETAILS] = CCA(X, Y, EDGES, OPTS) computes a low
% dimensional embedding Z in R^d that maximally preserves angles among  input 
% data X that lives in R^D, with the algorithm Conformal Component Analysis.
%
% The embedding Z is constrained to be Z = L*Y where Y is a partial basis that 
% spans the space of R^d. Such Y can be computed from graph Laplacian (such as 
% the outputs of Laplacian eigenmap and Locally Linear Embedding, ie, LLE). 
% The parameterization matrix L is found by this function as to maximally
% prserve angles between edges coded in the sparse matrix EDGES.
%
% A basic usage of this function is given below:
%
% Inputs:
%   X: input data stored in matrix  (D x N) where D is the dimensionality
%
%   Y: partial basis stored in matrix (d x N)
%
%   EDGES: a sparse matrix of (N x N). In each column i, the row indices j to
%   nonzero entrices define data points that are in the nearest neighbors of
%   data point i.
%
%   OPTS:
%     OPTS.method: 'CCA' 
%
% Outputs:
%   Z: low dimensional embedding (d X N)
%   CCAEIGN: eigenspectra of the matrix P = L'*L. If P is low-rank (say d' < d),
%   then Z can be cutoff at d' dimension as dimensionality reduced further.
%
% The CCA() function is fairly versatile. For more details, consult the file
% README.
%
%  by [e-mail] Aug 18, 2006
%  Feel free to use it for educational and research purpose.

% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.3b.
% 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

    % sanity check
    if nargin ~= 4
        error('Incorrect number of inputs supplied to cca().');
    end
    N = size(X,2);
    if (N~=size(Y,2)) | (N ~= size(EDGES,1)) | (N~=size(EDGES,2))
        disp('Unmatched matrix dimensions in cca().');
        fprintf('# of data points: %d\n', N);
        fprintf('# of data points in Y: %d\n', size(Y,2));
        fprintf('Size of the sparse matrix for edges: %d x %d\n', size(EDGES,1), size(EDGES,2));
        error('All above 4 numbers should be the same.');
    end
    % check necessary programs
    if exist('mexCCACollectData') ~= 3
        error('Missing mexCCACollectData mex file on the path');
    end
    if exist('csdp') ~= 2
        error('You will need CSDP solver to run cca(). Please make sure csdp.m is in your path');
    end
    % check options
    OPTS = check_opt(OPTS);

    D = size(X,1); d = size(Y,1);

    % 
    if OPTS.CCA == 1
        %disp('-------- Running CCA ---------------------------');
    else
        %disp('-------- Running MVU ---------------------------');
    end
    %disp('Step I.  collect data needed for SDP formulation');
    t0 = clock;
    [tnn, vidx] = triangNN(EDGES, OPTS.CCA);
    [erow, ecol, evalue] = sparse_nn(tnn);

    irow = int32(erow); icol = int32(ecol);
    ividx = int32(vidx); ivalue = int32(evalue);
    [A,B, g] = mexCCACollectData(X,Y, irow, icol, int32(OPTS.relative), ivalue, ividx );
    clear erow ecol irow icol tnn ividx ivalue evalue vidx;
    lst = find(g~=0);
    g = g(lst); B = B(:, lst);
    if OPTS.CCA == 1
        BG = B*spdiags(1./sqrt(g),0, length(g),length(g));
        Q = A - BG*BG';
        BIAS = OPTS.regularizer*reshape(eye(d), d^2,1);
    else
        Q = A; BIAS = 2*sum(B,2)+OPTS.regularizer*reshape(eye(d), d^2,1);
    end
    [V, E] = eig(Q+eye(size(Q))); % adding an identity matrix to Q for numerical
    E = E-eye(size(Q));           % stability
    E(E<0) = 0;
    if ~isreal(diag(E))
        warning('\tThe quadratic matrix is not positive definite..forced to be positive definite...\n');
        E=real(E);
        V = real(V);
        S = sqrt(E)*V';
        %keyboard;
    else
        S = sqrt(E)*V';
    end
    %fprintf('\tElapsed time: %ds\n', etime(clock,t0));

    % formulate the SDP problem
    %disp('Step II.  formulate the SDP problem');
    t0 = clock;
    [AA, bb, cc] = formulateSDP(S, d, BIAS, (OPTS.CCA==1));
    sizeSDP = d^2+1 + d + 2*(OPTS.CCA==1);
    csdppars.s = sizeSDP;
    csdpopts.printlevel = 0;
    %fprintf('\tElapsed time: %ds\n', etime(clock,t0));

    % solve it via csdp
    %disp('Step III. solve the SDP problem');
    t0 = clock;
    [xx, yy, zz, info] = csdp(AA, bb, cc, csdppars,csdpopts);
    %fprintf('\tSDP solver exit flag is %d\n', info);
    %fprintf('\tElapsed time: %ds\n', etime(clock,t0));
    ccaDetails.sdpflag = info;
    % the negate of yy is our solution
    yy = -yy;
    idx = 0;
    P = zeros(d);
    for col=1:d
        for row = col:d
            idx=idx+1;
            P(row, col) = yy(idx);
        end
    end
    % convert P to a positive definite matrix
    P = P+P' - diag(diag(P));

    % transform the original projection to the new
    [V, E] = eig(P);
    E(E<0) = 0; % make sure there is no very small negative eigenvalue
    L = diag(sqrt(diag(E))) * V';
    newY = L*Y;

    % eigenvalue of the new projection, doing PCA using covariance matrix
    % because the dimension of newY or Y is definitely less than the number of
    % points
    [newV, newE] = eig(newY *newY');
    newE = diag(newE);
    [dummy, idx] = sort(newE);
    newE = newE(idx(end:-1:1));
    newY = newV'*newY;
    Z = newY(idx(end:-1:1),:);

    ccaEigen = newE;
    ccaDetails.cost = P(:)'*Q*P(:) - BIAS'*P(:) + sum(g(:))*(OPTS.MVU==1);
    if OPTS.CCA == 1
        ccaDetails.c = spdiags(1./sqrt(g),0, length(g),length(g))*B'*P(:);
    else
        ccaDetails.c = [];
    end
    ccaDetails.P = P;
    ccaDetails.opts = OPTS;
    return;


    %%%%%%%%%%%%%%%%%%%% FOLLOWING IS SUPPORTING MATLAB FUNCTIONS
    function [A, b, c]=formulateSDP(S, D, bb, TRACE)
    [F0, FI, c] = localformulateSDP(S, D, bb, TRACE);
    [A, b, c] = sdpToSeDuMi(F0, FI, c);
    return

    function [F0, FI, c] = localformulateSDP(S, D, b, TRACE)
    % formulate SDP problem
    % each FI that corresponds to the LMI for the quadratic cost function has
    % precisely 2*D^2 nonzero elements. But we need only D^2 storage for
    % indexing these elements since the FI are symmetric
        tempFidx = zeros(D^2, 3);
        dimF = (D^2+1) + D + 2*TRACE;
        idx= 0;
        tracearray = ones(TRACE,1);
        for col=1:D
            for row=col:D
                idx = idx+1;
                lindx1 = sub2ind([D D], row, col);
                lindx2 = sub2ind([D D], col, row);
                tempFidx(:,1) = [1:D^2]';
                tempFidx(:,2) = D^2+1;
                if col==row
                    tempFidx(:,3) = S(:, lindx1) ;
                    FI{idx} = sparse([tempFidx(:,1); ...  % for cost function
                                        tempFidx(:,2); ... % symmetric
                                        row+D^2+1; ... % for P being p.s.d
                                        tracearray*(D^2+1+D+1); % for trace
                                        tracearray*(D^2+1+D+2); % for negate trace
                                    ], ...
                                    [tempFidx(:,2); ...  % for cost function
                                        tempFidx(:,1); ... % symmetric
                                        row+D^2+1; ... % for P being p.s.d
                                        tracearray*(D^2+1+D+1); % for trace
                                        tracearray*(D^2+1+D+2); % for negate trace
                                    ],...
                                    [tempFidx(:,3); ... % for cost function
                                        tempFidx(:,3); ... % symmetric
                                        1;                  % for P being p.s.d
                                        tracearray*1; % for trace
                                        tracearray*(-1); % for negate trace                               
                                    ], dimF, dimF);
                else

                    tempFidx(:,3) = S(:, lindx1) + S(:, lindx2);
                    FI{idx} = sparse([tempFidx(:,1); ...  % for cost function
                                        tempFidx(:,2); ... % symmetric
                                        row+D^2+1; ... % for P being p.s.d
                                        col+D^2+1; ... % symmetric
                                    ], ...
                                    [tempFidx(:,2); ...  % for cost function