svds {PRIMME}

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,
  ...
)

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.

u0

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

v0

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.

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

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

See Also

svd for computing all singular triplets; eigs_sym for computing a few eigenvalues and vectors from a symmetric/Hermitian matrix.

Examples

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