svds: Find a few singular values and vectors on large, sparse matrix

Description

Compute a few singular triplets from a specified region (the largest, the smallest, the closest to a point) on a matrix using PRIMME [1]. Only the matrix-vector product of the matrix is required. The used method is usually faster than a direct method (such as svd) if seeking few singular values and the matrix-vector product is cheap. For accelerating the convergence consider to use preconditioning and/or educated initial guesses.

Usage

svds(
  A,
  NSvals,
  which = "L",
  tol = 1e-06,
  u0 = NULL,
  v0 = NULL,
  orthou = NULL,
  orthov = NULL,
  prec = NULL,
  isreal = NULL,
  ...
)

Value

list with the next elements

d: the singular values $\sigma_i$
u: the left singular vectors $u_i$
v: the right singular vectors $v_i$
rnorms: the residual vector norms $\sqrt{\|Av - \sigma u\|^2+\|A^*u - \sigma v\|^2}$
stats$numMatvecs: matrix-vector products performed
stats$numPreconds: number of preconditioner applications performed
stats$elapsedTime: time expended by the eigensolver
stats$timeMatvec: time expended in the matrix-vector products
stats$timePrecond: time expended in applying the preconditioner
stats$timeOrtho: time expended in orthogonalizing
stats$estimateANorm: estimation of the norm of A

Arguments

A

matrix or a function with signature f(x, trans) that returns A %*% x when trans == "n" and t(Conj(A)) %*% x when trans == "c".

NSvals

number of singular triplets to seek.

which

which singular values to find:

"L": the largest values;

"S"

the smallest values;

vector of numbers

the closest values to these points.

tol

a triplet $(\sigma,u,v)$ is marked as converged when $\sqrt{\|Av - \sigma u\|^2+\|A^*u - \sigma v\|^2} \le tol\|A\|$ is smaller than $tol*||A||$, or close to the minimum tolerance that the selected method can achieve.

matrix whose columns are educated guesses of the left singular vectors to find.

matrix whose columns are educated guesses of the right singular vectors to find.

orthou

find left singular vectors orthogonal to the space spanned by the columns of this matrix; useful to avoid finding some triplets or to find new solutions.

orthov

find right singular vectors orthogonal to the space spanned by the columns of this matrix.

prec

preconditioner used to accelerated the convergence; it is a named list of matrices or functions such as solve(prec[[mode]],x) or prec[[mode]](x) return an approximation of $OP^{-1} x$, where

`mode`	$OP$
`"AHA"`	$A^*A$
`"AAH"`	$A A^*$
`"aug"`	$[0 A; A^* 0]$

The three values haven't to be set. It is recommended to set "AHA" for matrices with nrow > ncol; "AAH" for matrices with nrow < ncol; and additionally "aug" for tol < 1e-8.

isreal

whether A %*% x always returns real number and not complex.

...

other PRIMME options (see details).

Details

Optional arguments to pass to PRIMME eigensolver (see further details at [2]):

aNorm

estimation of norm-2 of A, used in convergence test (if not provided, it is estimated as the largest eigenvalue in magnitude seen)

maxBlockSize

maximum block size (like in subspace iteration or LOBPCG)

printLevel

message level reporting, from 0 (no output) to 5 (show all)

locking

1, hard locking; 0, soft locking

maxBasisSize

maximum size of the search subspace

minRestartSize

minimum Ritz vectors to keep in restarting

maxMatvecs

maximum number of matrix vector multiplications

iseed

an array of four numbers used as a random seed

method

which equivalent eigenproblem to solve

"primme_svds_normalequation": $A^*A$ or $AA^*$

"primme_svds_augmented"

$[0 A^*;A 0]$

"primme_svds_hybrid"

first normal equations and then augmented (default)

locking

1, hard locking; 0, soft locking

primmeStage1, primmeStage2

list with options for the first and the second stage solver; see eigs_sym

If method is "primme_svds_normalequation", the minimum tolerance that can be achieved is $\|A\|\epsilon/\sigma$, where $\epsilon$ is the machine precision. If method is "primme_svds_augmented" or "primme_svds_hybrid", the minimum tolerance is $\|A\|\epsilon$. However it may not return triplets with singular values smaller than $\|A\|\epsilon$.

References

[1] L. Wu, E. Romero and A. Stathopoulos, PRIMME_SVDS: A High-Performance Preconditioned SVD Solver for Accurate Large-Scale Computations, J. Sci. Comput., Vol. 39, No. 5, (2017), S248--S271.

[2] https://www.cs.wm.edu/~andreas/software/doc/svdsc.html#parameters-guide

Examples

Run this code

A <- diag(1:5,10,5)  # the singular values of this matrix are 1:10 and the
                        # left and right singular vectors are the columns of
                        # diag(1,100,10) and diag(10), respectively
r <- svds(A, 3);
r$d # the three largest singular values on A
r$u # the corresponding approximate left singular vectors
r$v # the corresponding approximate right singular vectors
r$rnorms # the corresponding residual norms
r$stats$numMatvecs # total matrix-vector products spend

r <- svds(A, 3, "S") # compute the three smallest values

r <- svds(A, 3, 2.5) # compute the three closest values to 2.5

A <- diag(1:500,500,100)   # we use a larger matrix to amplify the difference
r <- svds(A, 3, 2.5, tol=1e-3); # compute the values with 
r$rnorms                               # residual norm <= 1e-3*||A||

# Build the diagonal squared preconditioner
# and see how reduce the number matrix-vector products
P <- diag(colSums(A^2))
svds(A, 3, "S", tol=1e-3)$stats$numMatvecs
svds(A, 3, "S", tol=1e-3, prec=list(AHA=P))$stats$numMatvecs

# Passing A and the preconditioner as functions
Af <- function(x,mode) if (mode == "n") A%*%x else crossprod(A,x);
P = colSums(A^2);
PAHAf <- function(x) x / P;
r <- svds(Af, 3, "S", tol=1e-3, prec=list(AHA=PAHAf), m=500, n=100)

# Passing initial guesses
v0 <- diag(1,100,4) + matrix(rnorm(400), 100, 4)/100;
svds(A, 4, "S", tol=1e-3)$stats$numMatvecs
svds(A, 4, "S", tol=1e-3, v0=v0)$stats$numMatvecs

# Passing orthogonal constrain, in this case, already compute singular vectors
r <- svds(A, 4, "S", tol=1e-3); r$d
svds(A, 4, "S", tol=1e-3, orthov=r$v)$d

Run the code above in your browser using DataLab

Data Engineering and BI courses are free this week!