cca.m

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

                                        tempFidx(:,1); ... % symmetric
                                        col+D^2+1; ... % for P being p.s.d
                                        row+D^2+1; ... % being symmetric
                                    ],...
                                    [tempFidx(:,3); ... % for cost function
                                        tempFidx(:,3); ... % symmetric
                                        1;                  % for P being p.s.d
                                        1;                  % symmetric
                                    ], dimF, dimF);

                end
            end
        end
        idx=idx+1;
        % for the F matrix corresponding to t
        FI{idx} = sparse(D^2+1, D^2+1, 1, dimF, dimF);

        % now for F0
        if TRACE==1
            F0 = sparse( [[1:D^2] dimF-1 dimF], [[1:D^2] dimF-1 dimF], [ones(1, D^2) -1 1], dimF, dimF);
        else
            F0 = sparse( [[1:D^2]], [[1:D^2]], [ones(1, D^2)], dimF, dimF);
        end

        % now for c
        b = reshape(-b, D, D);
        b = b*2 - diag(diag(b)); 
        c = zeros(idx-1,1);
        kdx=0;
        %keyboard;
        for col=1:D
            for row=col:D
              kdx = kdx+1;
              c(kdx) = b(row, col);
            end
        end
        %keyboard;
        c = [c; 1]; % remember: we use only half of P
    return;


    function [A, b, c] = sdpToSeDuMi(F0, FI, cc)
    % convert the canonical SDP dual formulation:
    % (see  Vandenberche and Boyd 1996, SIAM Review)
    %  max -Tr(F0 Z)
    % s.t. Tr(Fi Z) = cci and Z is positive definite
    %
    % in which cc = (cc1, cc2, cc3,..) and FI = {F1, F2, F3,...}
    % 
    % to SeDuMi format (formulated as vector decision variables ):
    % min c'x
    % s.t. Ax = b and x is positive definite (x is a vector, so SeDuMi
    % really means that vec2mat(x) is positive definite)
    %
    % by [e-mail], June, 10, 2004

        if nargin < 3
            error('Cannot convert SDP formulation to SeDuMi formulation in sdpToSeDumi!');
        end

        [m, n] = size(F0);
        if m ~= n
            error('F0 matrix must be squared matrix in sdpToSeDumi(F0, FI, b)');
        end

        p = length(cc);
        if p ~= length(FI)
            error('FI matrix cellarray must have the same length as b in sdpToSeDumi(F0,FI,b)');
        end

        % should check every element in the cell array FI...later..

        % x = reshape(Z, n*n, 1);  % optimization variables from matrix to vector

        % converting objective function of the canonical SDP
        c = reshape(F0', n*n,1);

        % converting equality constraints of the canonical SDP
        zz= 0;
        for idx=1:length(FI)
            zz= zz + nnz(FI{idx});
        end
        A = spalloc( n*n, p, zz);
        for idx = 1:p
            temp = reshape(FI{idx}, n*n,1);
            lst = find(temp~=0);
            A(lst, idx) = temp(lst);
        end
        % The SeDuMi solver actually expects the transpose of A as in following
        % dual problem
        % max b'y
        % s.t. c - A'y is positive definite
        % Therefore, we transpose A
        % A = A';

        % b doesn't need to be changed
        b = cc;
    return;


    % Check OPTS that is passed into
    function OPTS = check_opt(OPTS)
        if isfield(OPTS,'method') == 0  
            OPTS.method = 'cca';
            disp('Options does''t have method field, so running CCA');
        end

        if strncmpi(OPTS.method, 'MVU',3)==1
            OPTS.CCA = 0; OPTS.MVU = 1;
        else
            OPTS.CCA = 1; OPTS.MVU = 0;
        end

        if isfield(OPTS, 'relative')==0
            OPTS.relative = 0;
        end

        if OPTS.CCA==1 && OPTS.relative ==1
            disp('Running CCA, so the .relative flag set to 0');
            OPTS.relative = 0;
        end

        if isfield(OPTS, 'regularizer')==0
            OPTS.regularizer = 0;
        end
    return



    function [tnn vidx]= triangNN(snn, TRI)
    % function [TNN VIDX]= triangNN(SNN) triangulates a sparse graph coded by spare matrix
    % SNN. TNN records the original edges in SNN as well as those that are
    % triangulated. Each edge is associated with a scaling factor that is specific
    % to a vertex. And VIDX records the id of the vertex.
    %
    % by [e-mail]  Aug. 15, 2006.

        N = size(snn,1);
        %fprintf('The graph has %d vertices\n', N);
        % figure out maximum degree a vertex has
        connectivs = sum(snn,1);
        maxDegree =  max(connectivs);
        tnn = spalloc(N, N, round(maxDegree*N)); % prealloc estimated storage for speedup 

        % triangulation
        for idx=1:N
            lst = find(snn(:, idx)>0);
            for jdx=1:length(lst)
                col = min (idx, lst(jdx));
                row = max(idx, lst(jdx));
                tnn(row, col) = tnn(row, col)+1;
                if TRI == 1
                    for kdx = jdx+1:length(lst)
                        col = min(lst(jdx), lst(kdx));
                        row = max(lst(jdx), lst(kdx));
                        tnn(row, col) = tnn(row, col)+1;
                    end
                end
            end
        end

        numEdges = nnz(tnn);
        %fprintf('Total %d edges\n', numEdges);
        numVertexIdx = full(sum(tnn(:)));
        %fprintf('%d vertex entries are needed\n', numVertexIdx);

        rowIdx = zeros(numVertexIdx,1);
        colIdx = zeros(numVertexIdx,1);
        vidx = zeros(numVertexIdx,1);

        whichEdge = 0;
        for idx=1:N
            lst = find(snn(:, idx)>0);
            for jdx=1:length(lst)
                col = min(lst(jdx), idx);
                row  = max(lst(jdx), idx);
                whichEdge = whichEdge+1;
                rowIdx(whichEdge) = row;
                colIdx(whichEdge) = col;
                vidx(whichEdge)  = idx;
                if TRI==1
                    for kdx = jdx+1:length(lst)
                        col = min(lst(jdx), lst(kdx));
                        row = max(lst(jdx), lst(kdx));
                        whichEdge = whichEdge+1;
                        rowIdx(whichEdge) = row;
                        colIdx(whichEdge) = col;
                        vidx(whichEdge)  = idx;
                    end
                end
            end
        end
        linearIdx  = sub2ind([N N],rowIdx, colIdx);
        [sa, sIdx] = sort(linearIdx);
        vidx = vidx(sIdx);
    return


    % turn sparse graph snn into row and col indices
    function [edgesrow, edgescol, value] = sparse_nn(snn)
        N = size(snn,1);
        edgescol = zeros(N+1,1);
        nnzer = nnz(snn);
        edgesrow = zeros(nnzer,1);
        value = zeros(nnzer,1);

        edgescol(1) = 0;
        for jdx=1:N
            lst = find(snn(:, jdx)>0);
            %lst = lst(find(lst>jdx));
            edgescol(jdx+1) = edgescol(jdx)+length(lst);
            edgesrow(edgescol(jdx)+1:edgescol(jdx+1)) = lst-1;
            value(edgescol(jdx)+1:edgescol(jdx+1)) = snn(lst, jdx);
        end
    return