dspca.m

稀疏PCA的优化解算法

function [U,X,x,F,k,dualitygap,cputime,perceigs]=DSPCA(A,rho,gapchange,maxiter,info,algo,U0,F0,k0)
% Wrapper function for sparse_rank_one_mex
% DSPCA finds a sparse rank-one approximation to 
% a given symmetric matrix A, by solving the SDP
% min_U lambda_max(A+X) : X = X', abs(X(i,j)) <= rho, 1<=i,j<= n
% and its dual:
% max_X Tr(UA) - rho sum_ij |U_ij| : U=U', U \succeq 0, Tr(U)=1
% ***	inputs: ***
% A                 nxn symmetric matrix (left unchanged)
% n                 problem size
% rho              non-negative scalar 
% gapchange required change in gap from first gap (default: 1e-4) 
% MaxIter		maximum number of iterations
% info             controls verbosity: 0 silent, n>0 frequency of progress report
% k0               0 if cold start, k0 if WarmStart (total number of iterations in previous run)
% F0               Average gradient (for warm start, Fmat is updated)
% X0               symmetric matrix that solves the above SDP 
 
% ***	outputs: ***
% U                 symmetric matrix that solves the above SDP 
% X                 dual variable, solves the dual SDP 
% x                 largest eigenvector of U 
% F                 Average gradient
% k                 number of iterations run
% dualitygap  vector of duality gaps at designated iterations
% cputime      vector of cumulative cpu times at designated iterations
% perceigs     vector of percentage of eigenvalues used (in partial
%                       eigenvalue decomposition) at designated iterations

% Refer to d扐spremont, A., El Ghaoui, L., Jordan, M. & Lanckriet, G. R. G. (2005), 
% 慉 direct formulation for sparse PCA using semidefinite programming? 
% Advances in Neural Information Processing Systems 17, 41?8.

% Compute sparsest pseudo eigenvector of A
if nargin < 5
    disp('Error: DSPCA must take at least 5 arguments');
    U=[];X=[];x=[];F=[];k=[];
    return;
elseif nargin == 5
    algo=1;
elseif nargin == 6
    % argument number is fine - nothing to do
elseif nargin == 9
    if isempty(algo)
        algo=1;
    end
    if isempty(U0) || isempty(F0) || isempty(k0)
        disp('Error: DSPCA needs values for U,F,and k0 for warm Startup');
        U=[];X=[];x=[];F=[];k=[];
        return;
    end
else
    disp('Error: DSPCA can take either 5,6, or 9 arguments');
    U=[];X=[];x=[];F=[];k=[];
    return;
end

% Test if A is 2 dimensional
if (length(size(A))~= 2)
    disp('Error: A is not 2 dimensional');
    U=[];X=[];x=[];F=[];k=[];
    return;
end
% Test if A is square
n=size(A,1);
m=size(A,2);
if (m~=n) 
    disp('Error: A is not square');
    U=[];X=[];x=[];F=[];k=[];
    return;
end

% Test if A is symmetric
checksymm = sum(sum(A-A'));
if (checksymm > 1e-10)
    disp('Error: A is not symmetric');
    U=[];X=[];x=[];F=[];k=[];
    return;
end

% check that variables make sense
if rho <= 0
    disp('Error: rho must be > 0');
    U=[];X=[];x=[];F=[];k=[];
    return;
end
if gapchange <= 0
    disp('Error: gapchange must be > 0');
    U=[];X=[];x=[];F=[];k=[];
    return;
end
if maxiter < 1
    disp('Error: maxiter must be >= 1');
    U=[];X=[];x=[];F=[];k=[];
    return;
end 
if info < 0
    disp('Error: info must be >= 0');
    U=[];X=[];x=[];F=[];k=[];
    return;
end
if (nargin == 6 || nargin == 9) && (algo~=1 && algo ~=2 && algo ~=3)
    disp('Error: the algorithm must be 1,2, or 3');
    U=[];X=[];x=[];F=[];k=[];
    return;
end
    
% Parameters that can be modified here and affect the matrix exponential computation using partial eigenvalue decomposition

% Run the mex file
if nargin < 14
    [U,X,x,F,k,dualitygap,cputime,perceigs] = sparse_rank_one_mex(A,rho,gapchange,maxiter,info,algo,numeigs,addeigs,maxit_eigs,checkgap,propeigs,check_for_more_eigs,check_for_unnec_eigs);
else
    [U,X,x,F,k,dualitygap,cputime,perceigs] = sparse_rank_one_mex(A,rho,gapchange,maxiter,info,algo,numeigs,addeigs,maxit_eigs,checkgap,propeigs,check_for_more_eigs,check_for_unnec_eigs,U0,F0,k0);
end
clear sparse_rank_one_mex;