rsvd (version 2.0.0)

rrpca: Randomized robust principal component analysis

Description

Estimate a low-rank component and an entrywise sparse component whose sum reconstructs a real matrix. Principal component pursuit uses penalties to encourage this separation, with randomized or deterministic singular value decompositions for the low-rank updates.

Usage

rrpca(
  A,
  lambda = NULL,
  maxiter = 50,
  tol = 1e-05,
  p = 10,
  q = 2,
  trace = FALSE,
  rand = TRUE
)

Value

A list of class "rrpca" with the following components.

L

A dense \(m\) by \(n\) matrix containing the estimated low-rank component.

S

A dense \(m\) by \(n\) matrix containing the estimated sparse component. It is stored densely even when most entries are zero.

err

A numeric vector of relative Frobenius reconstruction residuals, one per completed iteration. For a zero input, the function returns zero components and err = 0 immediately.

Rows of L and S correspond to retained input rows after missing rows have been omitted. No fitted loading basis or new-data prediction method is returned.

Arguments

A

A real numeric matrix or data frame to decompose. A sparse matrix from Matrix is accepted and converted to a dense matrix. Rows containing missing values are omitted with a warning; all remaining entries must be finite. Complex input is not supported. In the dimensions below, \(m\) and \(n\) are the numbers of rows and columns after omission; both must be positive.

lambda

A positive numeric penalty on the entrywise sparse component. The default, NULL, uses 1 / sqrt(max(m, n)). Larger values penalize the sparse component more strongly.

maxiter

A positive integer maximum number of iterations (default 50). Reaching this limit does not imply convergence.

tol

A positive numeric tolerance for both the relative reconstruction residual and stationarity diagnostic (default 1e-5).

p

A nonnegative integer oversampling parameter for rsvd (default 10); used only in randomized SVD steps.

q

A nonnegative integer initial number of power iterations for rsvd (default 2); increased adaptively if randomized accuracy stalls. Ignored when rand = FALSE.

trace

Logical; print iteration progress, including both convergence diagnostics (default FALSE).

rand

Logical; start with randomized SVD steps when TRUE (default), with deterministic SVD used for large target ranks or stalled accuracy. With FALSE, use the deterministic SVD throughout.

Author

N. Benjamin Erichson

Details

A matrix can contain patterns shared across many entries together with corruption confined to a few entries. For example, frames of a nearly static video can share a low-rank background while a small moving object affects only a subset of the pixels in each frame. Principal component pursuit seeks a low-rank term \(L\) and a sparse term \(S\) satisfying \(A = L + S\).

Among such splits, the objective minimizes the sum of singular values of \(L\) plus lambda times the sum of absolute entries of \(S\). These penalties are called the nuclear norm and the entrywise one-norm, respectively. They encourage low rank in \(L\) and few nonzero entries in \(S\). Successful denoising depends on the data and the assumptions that make low-rank structure distinguishable from sparse corruption.

Computation

The inexact augmented Lagrange multiplier method (IALM) first updates the sparse component by shrinking entries toward zero and setting small entries to zero. It then applies the same idea to singular values to update the low-rank component. These steps are entrywise soft thresholding and singular-value thresholding. The surviving singular values determine the rank, so the function selects it during iteration.

Randomized singular-value thresholding increases its target rank if all returned singular values exceed the shrinkage threshold, allowing it to examine further directions that might survive thresholding. Large target ranks use deterministic SVD. If accuracy stalls, the algorithm makes bounded increases in the number of power iterations before falling back to deterministic SVD. These refinements can increase runtime.

Input and working arrays are dense, including for sparse input. Apply any desired centering or scaling before calling rrpca; the function uses the supplied values without this preprocessing. Use set.seed before fitting for reproducible random draws when rand = TRUE.

Stopping diagnostics

Many different splits can reconstruct the input, so reconstruction accuracy alone is insufficient to establish convergence. The algorithm also monitors the change in the low-rank component.

The relative reconstruction residual combines errors across all entries and compares their Frobenius norm with that of the input, $$\frac{\|A-L-S\|_F}{\|A\|_F}.$$ A small value means that \(L+S\) closely reproduces \(A\) relative to its magnitude. With the sparse component updated first, the stationarity diagnostic measures the change in \(L\) scaled by the current penalty, $$\frac{\mu\|L_{new}-L_{old}\|_F}{\max(1,\|Z\|_F)}.$$ Here \(Z\) is the updated Lagrange multiplier and \(\mu\) is the penalty used for that iteration. Both diagnostics must be at most tol for the stopping test to be satisfied. The penalty increases only when the stationarity diagnostic is at most tol. These normalized tests are invariant to rescaling the input.

The iteration limit maxiter may be reached before both diagnostics satisfy tol. The stopping diagnostics and randomized approximations provide no general certificate of optimality. Use trace = TRUE to display both diagnostics during fitting. The returned err records reconstruction residuals only; a small final value cannot establish convergence on its own.

Using the decomposition

rrpca fits all entries of the supplied matrix jointly. Its L and S describe those observations, and a fitted PCA rotation or a predict.rrpca method is unavailable. Projecting new rows onto a basis extracted from L gives an ordinary low-rank projection with a generally dense residual. Estimating a robust sparse component for those new rows with a fixed learned basis requires a separate optimization problem. For ordinary PCA with new-data scores, use rpca and its predict method.

References

• N. B. Erichson, S. Voronin, S. L. Brunton and J. N. Kutz. 2019. Randomized Matrix Decompositions Using R. Journal of Statistical Software, 89(11), 1--48. tools:::Rd_expr_doi("10.18637/jss.v089.i11").

• Z. Lin, M. Chen and Y. Ma. 2010. The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices. https://arxiv.org/abs/1009.5055.

• S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein. 2011. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. Foundations and Trends in Machine Learning, 3(1), 1--122. tools:::Rd_expr_doi("10.1561/2200000016").

See Also

rpca, rsvd.

Examples

# A small video: each column is one frame, with pixels stored in rows.
background <- outer(1:32, 1:32, function(x, y) 0.1 + 0.001 * (x + y))
video <- matrix(rep(as.vector(background), 30), nrow = 32 * 32)
for (i in 1:30) {
  foreground <- matrix(0, 32, 32)
  foreground[i:(i + 2), 15:17] <- 0.8
  video[, i] <- video[, i] + as.vector(foreground)
}

set.seed(42)
out <- rrpca(video, maxiter = 100, tol = 1e-6)
tail(out$err, 1) # Reconstruction residual; see Details for convergence.

# Compare one frame and its estimated background and foreground.
# Each panel uses its own intensity range.
old_par <- par(mfrow = c(1, 4), mar = c(0, 0, 2, 0))
image(background, col = gray.colors(64), axes = FALSE,
      asp = 1, useRaster = TRUE, main = "Background")
image(matrix(video[, 10], 32, 32), col = gray.colors(64),
      axes = FALSE, asp = 1, useRaster = TRUE, main = "Frame")
image(matrix(out$L[, 10], 32, 32), col = gray.colors(64),
      axes = FALSE, asp = 1, useRaster = TRUE, main = "Low-rank part")
image(matrix(out$S[, 10], 32, 32), col = gray.colors(64),
      axes = FALSE, asp = 1, useRaster = TRUE, main = "Sparse part")
par(old_par)