4
$\begingroup$

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?

The data is ad hoc:

> 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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – seanv507 Aug 2 '13 at 8:49
  • $\begingroup$ What is the package that you are using that includes the nsprcomp function? $\endgroup$ – Marc in the box Aug 2 '13 at 9:12
  • $\begingroup$ @Marcinthebox Thx. See the update please. $\endgroup$ – SolessChong Aug 2 '13 at 10:12
10
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
3
$\begingroup$

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 .
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.