RobustPCA: Decompose a Matrix into Low-Rank and Sparse Components

Description

Suppose we have a data matrix, which is the superposition of a low-rank component and a sparse component. Candes, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11. prove that we can recover each component individually under some suitable assumptions. It is possible to recover both the low-rank and the sparse components exactly by solving a very convenient convex program called Principal Component Pursuit; among all feasible decompositions, simply minimize a weighted combination of the nuclear norm and of the L1 norm. This package implements this decomposition algorithm resulting with Robust PCA approach.

Details

Index of help topics:

F2norm                  Frobenius norm of a matrix
rpca                    Decompose a matrix into a low-rank component
                        and a sparse component by solving Principal
                        Components Pursuit
rpca-package            RobustPCA: Decompose a Matrix into Low-Rank and
                        Sparse Components
thresh.l1               Shrinkage operator
thresh.nuclear          Thresholding operator

This package contains rpca function, which decomposes a rectangular matrix M into a low-rank component, and a sparse component, by solving a convex program called Principal Component Pursuit:

\textrm{minimize}\quad \|L\|_{*} + \lambda \|S\|_{1}

\textrm{subject to}\quad L+S = M

where \|L\|_{*} is the nuclear norm of L (sum of singular values).

Note

Use citation("rpca") to cite this R package.

Author(s)

Maciek Sykulski [aut, cre]

Maintainer: Maciek Sykulski <macieksk@gmail.com>

References

Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11.

Yuan, X., & Yang, J. (2009). Sparse and low-rank matrix decomposition via alternating direction methods. preprint, 12.

See Also

rpca

Frobenius norm of a matrix

Description

Frobenius norm of a matrix.

Usage

F2norm(M)

Arguments

Value

Frobenius norm of M.

Examples

## The function is currently defined as
function (M) 
sqrt(sum(M^2))

F2norm(matrix(runif(100),nrow=5))

Decompose a matrix into a low-rank component and a sparse component by solving Principal Components Pursuit

Description

This function decomposes a rectangular matrix M into a low-rank component, and a sparse component, by solving a convex program called Principal Component Pursuit.

Usage

rpca(M, 
     lambda = 1/sqrt(max(dim(M))), mu = prod(dim(M))/(4 * sum(abs(M))), 
     term.delta = 10^(-7), max.iter = 5000, trace = FALSE,
     thresh.nuclear.fun = thresh.nuclear, thresh.l1.fun = thresh.l1, 
     F2norm.fun = F2norm)

Arguments

Details

These functions decompose a rectangular matrix M into a low-rank component, and a sparse component, by solving a convex program called Principal Component Pursuit:

\textrm{minimize}\quad \|L\|_{*} + \lambda \|S\|_{1}

\textrm{subject to}\quad L+S = M

where \|L\|_{*} is the nuclear norm of L (sum of singular values).

Value

The function returns two matrices S and L, which have the property that L+S \simeq M, where the quality of the approximation depends on the argument term.delta, and the convergence of the algorithm.

Author(s)

Maciek Sykulski [aut, cre]

References

Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11.

Yuan, X., & Yang, J. (2009). Sparse and low-rank matrix decomposition via alternating direction methods. preprint, 12.

Examples

data(iris)
M <- as.matrix(iris[,1:4])
Mcent <- sweep(M,2,colMeans(M))

res <- rpca(Mcent)

## Check convergence and number of iterations
with(res$convergence,list(converged,iterations))
## Final delta F2 norm divided by F2norm(Mcent)
with(res$convergence,final.delta)

## Check properites of the decomposition
with(res,c(
all(abs( L+S - Mcent ) < 10^-5),
all( L == L.svd$u%*%(L.svd$d*L.svd$vt) )
))
# [1] TRUE TRUE

## The low rank component has rank 2
length(res$L.svd$d)
## However, the sparse component is not sparse 
## - thus this data set is not the best example here.
mean(res$S==0)

## Plot the first (the only) two principal components
## of the low-rank component L
rpc<-res$L.svd$u%*%diag(res$L.svd$d)
plot(jitter(rpc[,1:2],amount=.001),col=iris[,5])

## Compare with classical principal components
pc <- prcomp(M,center=TRUE)
plot(pc$x[,1:2],col=iris[,5])
points(rpc[,1:2],col=iris[,5],pch="+")

## "Sparse" elements distribution
plot(density(abs(res$S),from=0))
curve(dexp(x,rate=1/mean(abs(res$S))),add=TRUE,lty=2)

## Plot measurements against measurements corrected by sparse components
par(mfcol=c(2,2))
for(i in 1:4) {
plot(M[,i],M[,i]-res$S[,i],col=iris[,5],xlab=colnames(M)[i])
}

Shrinkage operator

Description

Shrinkage operator: S[x] = sgn(x) max(|x| - thr, 0). For description see section 5 of Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?.

Usage

thresh.l1(x, thr)

Arguments

Value

S[x] = sgn(x) max(|x| - thr, 0)

References

Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11

Yuan, X., & Yang, J. (2009). Sparse and low-rank matrix decomposition via alternating direction methods. preprint, 12.

See Also

thresh.nuclear

Examples

## The function is currently defined as
function(x,thr){sign(x)*pmax(abs(x)-thr,0)}

summary(thresh.l1(runif(100),0.3))

Thresholding operator

Description

Thresholding operator, an application of the shrinkage operator on a singular value decomposition: D[X] = U S[Sigma] V . For description see section 5 of Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?.

Usage

thresh.nuclear(M, thr)

Arguments

Value

Returned is a thresholded Singular Value Decomposition with thr subtracted from singular values, and values smaller than 0 dropped together with their singular vectors.

References

Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM), 58(3), 11

Yuan, X., & Yang, J. (2009). Sparse and low-rank matrix decomposition via alternating direction methods. preprint, 12.

See Also

thresh.l1

Examples

## The function is currently defined as
function (M, thr) {
    s <- La.svd.cmp(M)
    dd <- thresh.l1(s$d, thr)
    id <- which(dd != 0)
    s$d <- dd[id]
    s$u <- s$u[, id, drop = FALSE]
    s$vt <- s$vt[id, , drop = FALSE]
    s$L <- s$u %*% (s$d * s$vt)
    s
  }

l<-thresh.nuclear(matrix(runif(600),nrow=20),2)
l$d