# Why such a poor result from sparse PCA R package?

I'm preparing to use R to perform sparse analysis on my data. I tried to get started with an ad hoc example, but the reconstruction result turned out really poor. I'm wondering if I was making any mistake in the code, or was I encountering any intrinsic drawbacks of the method?

> prod
[,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    2
[2,]    1    2    1    2    3
[3,]    0    1    0    1    1
[4,]    0    0    1    1    1
[5,]    0    1    0    1    1


and I use this

rst = nsprcomp(prod, retx=T, nneg=T)


to perform non-negtive sparse PCA.

I reconstruct the data like this:

recon = rst$x %*% t(rst$rotation) + matrix(1,5,1) %*% rst$center  The final result turned out to be really poor: > abs(prod - recon) [,1] [,2] [,3] [,4] [,5] [1,] 0.386523767 0.2432819 0.1186152 0.3535836 0.3770144 [2,] 0.184944985 1.1214659 0.6230812 1.2173770 1.8819417 [3,] 0.096185913 0.4301642 0.2249522 0.4825192 0.6942987 [4,] 0.009206957 0.5044193 0.2917920 0.6059221 0.8703587 [5,] 0.096185913 0.4301642 0.2249522 0.4825192 0.6942987  Any idea or comment is appreciated. ## Edit The package is called nsprcomp • can you explain your example? PCA etc dimension reduction require a statistical relationship (correlation/linear relation between coordinates) to perform well... is this the case with your data? – seanv507 Aug 2 '13 at 8:49 • What is the package that you are using that includes the nsprcomp function? – Marc in the box Aug 2 '13 at 9:12 • @Marcinthebox Thx. See the update please. – SolessChong Aug 2 '13 at 10:12 ## 2 Answers nsprcomp computes the scores matrix Z (rst$x in your example) as $Z=XW$, where $X$ is the data matrix (prod in your example) and $W$ is the matrix of principal axes (rst$rotation in your example). This is in accordance with standard PCA and the predict.prcomp interface. However, non-negative sparse PCA usually results in principal axes which are not pairwise orthogonal, and therefore a reconstruction$\hat{X} = ZW^t = XWW^t$doesn't recover$X$even if$W$has full rank, because$W$is not an orthogonal matrix. If you reconstruct using the pseudo-inverse$W^\dagger=(W^tW)^{-1}W^t$instead,$\hat{X}_2 = ZW^\dagger = XW(W^tW)^{-1}W^t$corresponds to an orthogonal projection of$X$onto the principal subspace spanned by$W$, and recovers$X$if$W$has full rank: library(MASS) recon2 = predict(rst) %*% ginv(rst$rotation) + matrix(1,5,1) %*% rst$center abs(prod - recon2) [,1] [,2] [,3] [,4] [,5] [1,] 4.440892e-16 4.440892e-16 2.220446e-16 8.881784e-16 2.220446e-16 [2,] 4.440892e-16 1.332268e-15 0.000000e+00 1.998401e-15 4.440892e-16 [3,] 1.110223e-16 4.440892e-16 2.220446e-16 4.440892e-16 2.220446e-16 [4,] 3.330669e-16 4.440892e-16 2.220446e-16 2.220446e-16 0.000000e+00 [5,] 1.110223e-16 4.440892e-16 2.220446e-16 4.440892e-16 2.220446e-16  I will amend the documentation accordingly. @Marc: nneg=TRUE enforces non-negative loadings, i.e. the principal axes are constrained to the non-negative orthant. Your code looks OK for the reconstruction, but this does not seem to be appropriate when the argument nneg=TRUE ("a logical value indicating whether the principal axes should be constrained to the non-negative orthant"). When this argument is set to FALSE, then the reconstruction works in the typical way: rst.f = nsprcomp(prod, retx=T, nneg=FALSE) recon = rst.f$x %*% t(rst.f$rotation) + matrix(1,5,1) %*% rst.f$center
abs(prod - recon)
[,1]         [,2]         [,3]         [,4] [,5]
[1,] 1.110223e-16 1.110223e-16 0.000000e+00 2.220446e-16    0
[2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00    0
[3,] 5.551115e-17 0.000000e+00 2.220446e-16 2.220446e-16    0
[4,] 0.000000e+00 4.440892e-16 0.000000e+00 2.220446e-16    0
[5,] 5.551115e-17 0.000000e+00 2.220446e-16 2.220446e-16    0


Sorry that I can't comment more on the use of the "non-negative orthant" option.

I might also point out that your example is not truely using a "sparse" matrix. But this method does appear do deal nicely with such objects:

require(Matrix)
prod <- replace(prod, prod==0, NaN)
tmp <- cbind(expand.grid(i=seq(nrow(prod)), j=seq(ncol(prod))), p=c(prod))[-which(is.na(c(prod))),]
prod2 <- sparseMatrix(i=tmp$i, j=tmp$j, x=tmp$p) rst.f = nsprcomp(prod2, retx=T, nneg=FALSE) recon = rst.f$x %*% t(rst.f$rotation) + matrix(1,5,1) %*% rst.f$center
abs(prod2 - recon)
5 x 5 sparse Matrix of class "dgCMatrix"

[1,] 2.220446e-16 1.110223e-16 .            .            .
[2,] 2.220446e-16 4.440892e-16 .            .            .
[3,] 1.110223e-16 2.220446e-16 2.220446e-16 2.220446e-16 .
[4,] .            .            .            .            .
[5,] 1.110223e-16 2.220446e-16 2.220446e-16 2.220446e-16 .