function [U, S, V] = mySVD(X,ReducedDim) %mySVD Accelerated singular value decomposition. % [U,S,V] = mySVD(X) produces a diagonal matrix S, of the % dimension as the rank of X and with nonnegative diagonal elements in % decreasing order, and unitary matrices U and V so that % X = U*S*V'. % % [U,S,V] = mySVD(X,ReducedDim) produces a diagonal matrix S, of the % dimension as ReducedDim and with nonnegative diagonal elements in % decreasing order, and unitary matrices U and V so that % Xhat = U*S*V' is the best approximation (with respect to F norm) of X % among all the matrices with rank no larger than ReducedDim. % % Based on the size of X, mySVD computes the eigvectors of X*X^T or X^T*X % first, and then convert them to the eigenvectors of the other. % % See also SVD. % % version 2.0 --Feb/2009 % version 1.0 --April/2004 % % Written by Deng Cai (dengcai AT gmail.com) % MAX_MATRIX_SIZE = 1600; % You can change this number according your machine computational power EIGVECTOR_RATIO = 0.1; % You can change this number according your machine computational power if ~exist('ReducedDim','var') ReducedDim = 0; end [nSmp, mFea] = size(X); if mFea/nSmp > 1.0713 ddata = X*X'; ddata = max(ddata,ddata'); dimMatrix = size(ddata,1); if (ReducedDim > 0) && (dimMatrix > MAX_MATRIX_SIZE) && (ReducedDim < dimMatrix*EIGVECTOR_RATIO) option = struct('disp',0); [U, eigvalue] = eigs(ddata,ReducedDim,'la',option); eigvalue = diag(eigvalue); else if issparse(ddata) ddata = full(ddata); end [U, eigvalue] = eig(ddata); eigvalue = diag(eigvalue); [dump, index] = sort(-eigvalue); eigvalue = eigvalue(index); U = U(:, index); end clear ddata; maxEigValue = max(abs(eigvalue)); eigIdx = find(abs(eigvalue)/maxEigValue < 1e-10); eigvalue(eigIdx) = []; U(:,eigIdx) = []; if (ReducedDim > 0) && (ReducedDim < length(eigvalue)) eigvalue = eigvalue(1:ReducedDim); U = U(:,1:ReducedDim); end eigvalue_Half = eigvalue.^.5; S = spdiags(eigvalue_Half,0,length(eigvalue_Half),length(eigvalue_Half)); if nargout >= 3 eigvalue_MinusHalf = eigvalue_Half.^-1; V = X'*(U.*repmat(eigvalue_MinusHalf',size(U,1),1)); end else ddata = X'*X; ddata = max(ddata,ddata'); dimMatrix = size(ddata,1); if (ReducedDim > 0) && (dimMatrix > MAX_MATRIX_SIZE) && (ReducedDim < dimMatrix*EIGVECTOR_RATIO) option = struct('disp',0); [V, eigvalue] = eigs(ddata,ReducedDim,'la',option); eigvalue = diag(eigvalue); else if issparse(ddata) ddata = full(ddata); end [V, eigvalue] = eig(ddata); eigvalue = diag(eigvalue); [dump, index] = sort(-eigvalue); eigvalue = eigvalue(index); V = V(:, index); end clear ddata; maxEigValue = max(abs(eigvalue)); eigIdx = find(abs(eigvalue)/maxEigValue < 1e-10); eigvalue(eigIdx) = []; V(:,eigIdx) = []; if (ReducedDim > 0) && (ReducedDim < length(eigvalue)) eigvalue = eigvalue(1:ReducedDim); V = V(:,1:ReducedDim); end eigvalue_Half = eigvalue.^.5; S = spdiags(eigvalue_Half,0,length(eigvalue_Half),length(eigvalue_Half)); eigvalue_MinusHalf = eigvalue_Half.^-1; U = X*(V.*repmat(eigvalue_MinusHalf',size(V,1),1)); end