MATLAB Interface¶
-
function [varargout] = primme_svds(varargin)
primme_svds()
finds a few singular values and vectors of a matrixA
by calling PRIMME.A
is typically large and sparse.S = primme_svds(A)
returns a vector with the 6 largest singular values ofA
.S = primme_svds(AFUN,M,N)
accepts the function handleAFUN
to perform the matrix vector products with an M-by-N matrixA
.AFUN(X,'notransp')
returnsA*X
whileAFUN(X,'transp')
returnsA’*X
. In all the following,A
can be replaced byAFUN,M,N
.S = primme_svds(A,k)
computes thek
largest singular values ofA
.S = primme_svds(A,k,sigma)
computes thek
singular values closest to the scalar shiftsigma
.- If
sigma
is a vector, find the singular valueS(i)
closest to eachsigma(i)
, fori<=k
. - If
sigma
is'L'
, it computes the largest singular values. - if
sigma
is'S'
, it computes the smallest singular values.
S = primme_svds(A,k,sigma,OPTIONS)
specifies extra solver parameters. Some default values are indicated in brackets {}:aNorm
: estimation of the 2-norm ofA
{0.0 (estimate the norm internally)}tol
: convergence toleranceNORM([A*V-U*S;A'*U-V*S]) <= tol * NORM(A)
(seeeps
) {1e-10
for double precision and1e-3
for single precision}maxit
: maximum number of matvecs withA
andA'
(seemaxMatvecs
) {inf}p
: maximum basis size (seemaxBasisSize
)disp
: level of reporting 0-3 (see HIST) {0: no output}isreal
: if 0, the matrix is complex; else it's real {0: complex}isdouble
: if 0, the matrix is single; else it's double {1: double}method
: which equivalent eigenproblem to solve- '
primme_svds_normalequations
':A'*A
orA*A'
- '
primme_svds_augmented
':[0 A';A 0]
- '
primme_svds_hybrid
': first normal equations and then augmented (default)
- '
u0
: initial guesses to the left singular vectors (seeinitSize
) {[]}v0
: initial guesses to the right singular vectors {[]}orthoConst
: external orthogonalization constraints (seenumOrthoConst
) {[]}locking
: 1, hard locking; 0, soft lockingmaxBlockSize
: maximum block sizeiseed
: random seedprimme
: options for first stage solverprimmeStage2
: options for second stage solverconvTestFun
: function handler with an alternative convergence criterion. IfFUN(SVAL,LSVEC,RSVEC,RNORM)
returns a nonzero value, the triplet(SVAL,LSVEC,RSVEC)
with residual normRNORM
is considered converged.
The available options for
OPTIONS.primme
andprimmeStage2
are the same asprimme_eigs()
, plus the option'method'
.S = primme_svds(A,k,sigma,OPTIONS,P)
applies a preconditionerP
as follows:- If
P
is a matrix it appliesP\X
andP'\X
to approximateA\X
andA'\X
. - If
P
is a function handle,PFUN
,PFUN(X,'notransp')
returnsP\X
andPFUN(X,'transp')
returnsP’\X
, approximatingA\X
andA'\X
respectively. - If
P
is astruct
, it can have one or more of the following fields: P.AHA\X
orP.AHA(X)
returns an approximation of(A'*A)\X
,P.AAH\X
orP.AAH(X)
returns an approximation of(A*A')\X
,P.aug\X
orP.aug(X)
returns an approximation of[zeros(N,N) A';A zeros(M,M)]\X
.
- If
- If
P
is[]
then no preconditioner is applied.
S = primme_svds(A,k,sigma,OPTIONS,P1,P2
) applies a factorized preconditioner:- If both
P1
andP2
are nonempty, apply(P1*P2)\X
to approximateA\X
. - If
P1
is[]
andP2
is nonempty, then(P2'*P2)\X
approximatesA'*A
.P2
can be the R factor of an (incomplete) QR factorization ofA
or the L factor of an (incomplete) LL' factorization ofA'*A
(RIF). - If both
P1
andP2
are[]
then no preconditioner is applied.
[U,S,V] = primme_svds(...)
returns also the corresponding singular vectors. IfA
is M-by-N andk
singular triplets are computed, thenU
is M-by-k with orthonormal columns,S
is k-by-k diagonal, andV
is N-by-k with orthonormal columns.[S,R] = primme_svds(...)
[U,S,V,R] = primme_svds(...)
returns the residual norm of eachk
triplet,NORM([A*V(:,i)-S(i,i)*U(:,i); A'*U(:,i)-S(i,i)*V(:,i)])
.[U,S,V,R,STATS] = primme_svds(...)
returns how many timesA
andP
were used and elapsed time. The application ofA
is counted independently from the application ofA'
.[U,S,V,R,STATS,HIST] = primme_svds(...)
returns the convergence history, instead of printing it. Every row is a record, and the columns report:HIST(:,1)
: number of matvecsHIST(:,2)
: timeHIST(:,3)
: number of converged/locked tripletsHIST(:,4)
: stageHIST(:,5)
: block indexHIST(:,6)
: approximate singular valueHIST(:,7)
: residual normHIST(:,8)
: QMR residual norm
OPTS.disp
controls the granularity of the record. IfOPTS.disp == 1
,HIST
has one row per converged triplet and only the first four columns are reported; ifOPTS.disp == 2
,HIST
has one row per outer iteration and only the first seven columns are reported; and otherwiseHIST
has one row per QMR iteration and all columns are reported.Examples:
A = diag(1:50); A(200,1) = 0; % rectangular matrix of size 200x50 s = primme_svds(A,10) % the 10 largest singular values s = primme_svds(A,10,'S') % the 10 smallest singular values s = primme_svds(A,10,25) % the 10 closest singular values to 25 opts = struct(); opts.tol = 1e-4; % set tolerance opts.method = 'primme_svds_normalequations' % set svd solver method opts.primme.method = 'DEFAULT_MIN_TIME' % set first stage eigensolver method opts.primme.maxBlockSize = 2; % set block size for first stage [u,s,v] = primme_svds(A,10,'S',opts); % find 10 smallest svd triplets opts.orthoConst = {u,v}; [s,rnorms] = primme_svds(A,10,'S',opts) % find another 10 % Compute the 5 smallest singular values of a rectangular matrix using % Jacobi preconditioner on (A'*A) A = sparse(diag(1:50) + diag(ones(49,1), 1)); A(200,50) = 1; % size(A)=[200 50] P = diag(sum(abs(A).^2)); precond.AHA = @(x)P\x; s = primme_svds(A,5,'S',[],precond) % find the 5 smallest values % Estimation of the smallest singular value A = diag([1 repmat(2,1,1000) 3:100]); [~,sval,~,rnorm] = primme_svds(A,1,'S',struct('convTestFun',@(s,u,v,r)r<s*.1)); sval - rnorm % approximate smallest singular value
See also: MATLAB svds,
primme_eigs()
- If