naginterfaces.library.rnla.svd_​rowext_​real

naginterfaces.library.rnla.svd_rowext_real(jobu, jobvt, a, k, statecomm, rtol_abs=- 1.0, rtol_rel=- 1.0)

svd_rowext_real computes the singular value decomposition (SVD) of a real $m\times n$ matrix $A$, optionally computing the left and/or right singular vectors by using a randomised numerical linear algebra (RNLA) method.

For full information please refer to the NAG Library document for f10ca

https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/f10/f10caf.html

Parameters

jobustr, length 1

Specifies options for computing part of or none of the matrix $U$.

$jobu=\text{'S'}$

The first $k$ columns of $U$ (the left singular vectors) are returned in the array $u$.

$jobu=\text{'N'}$

No columns of $U$ (no left singular vectors) are computed.

jobvtstr, length 1

Specifies options for computing part of or none of the matrix $V^T$.

$jobvt=\text{'S'}$

The first $k$ rows of $V^T$ (the right singular vectors) are returned in the array $vt$.

$jobvt=\text{'N'}$

No rows of $V^T$ (no right singular vectors) are computed.

afloat, array-like, shape $(m,n)$

The $m\times n$ matrix $A$.

kint

$k$, number of columns in random projection, $Y=A\Omega$.

statecommdict, RNG communication object, modified in place

RNG communication structure.

This argument must have been initialized by a prior call to rand.init_repeat or rand.init_nonrepeat.

rtol_absfloat, optional

The absolute tolerance, ${ϵ}_a$ used in defining the threshold on estimating the rank of $A$. If $rtol\_abs<0.0$ then $0.0$ is used unless $rtol\_rel\le 0.0$ in which case $\sqrt{\text{machine precision}}$ is used.

rtol_relfloat, optional

The relative tolerance, ${ϵ}_r$ used in defining the threshold on estimating the rank of $A$. If $rtol\_rel<0.0$ then $0.0$ is used unless $rtol\_abs\le 0.0$ in which case $\sqrt{\text{machine precision}}$ is used.

Returns

sfloat, ndarray, shape $\left(r\right)$

The $r$ largest singular values of $A$ in descending order.

uNone or float, ndarray, shape $(:,:)$

If $jobu=\text{'S'}$, $u$ contains the first $r$ columns of $U$ (the left singular vectors, stored column-wise).

vtNone or float, ndarray, shape $(:,:)$

If $jobvt=\text{'S'}$, $vt$ contains the first $r$ rows of $V^T$ (the right singular vectors).

rint

$r$, contains estimated rank of array $a$.

Raises

NagValueError

(errno $1$)

On entry, $jobu=⟨value⟩$.

Constraint: $jobu=\text{'S'}$ or $\text{'N'}$.

(errno $2$)

On entry, $jobvt=⟨value⟩$.

Constraint: $jobvt=\text{'S'}$ or $\text{'N'}$.

(errno $3$)

On entry, $m=⟨value⟩$.

Constraint: $m>0$.

(errno $4$)

On entry, $n=⟨value⟩$.

Constraint: $n>0$.

(errno $7$)

On entry, $k=⟨value⟩$.

Constraint: $0<k\le min(m,n)$.

(errno $8$)

On entry, $statecomm$[‘state’] vector has been corrupted or not initialized.

(errno $9$)

On entry, $jobu=⟨value⟩$ and $\text{ldu}=⟨value⟩$.

Constraint: if $jobu=\text{'S'}$, $\text{ldu}\ge m$.

(errno $10$)

On entry, $jobvt=⟨value⟩$ and $\text{ldvt}=⟨value⟩$.

Constraint: if $jobvt=\text{'S'}$, $\text{ldvt}\ge k$.

(errno $21$)

On exit, $r=k$, the rank of $A$ may be larger than $r$.

Increase $k$ to obtain a more accurate rank estimate.

Smallest diagonal element of $R$, from $QR$ of $Y$, $=$$⟨value⟩$.

Tolerance used to determine rank $=$$⟨value⟩$.

(errno $22$)

$A$ has effective rank of zero.

First diagonal element of $R$, from $QR$ of $Y$, $=$$⟨value⟩$.

Tolerance used to determine rank $=$$⟨value⟩$.

Notes

The SVD is written as

$$A=U\Sigma V^T\text{,}$$

where $\Sigma$ is an $m\times n$ matrix which is zero except for its $min(m,n)$ diagonal elements, $U$ is an $m\times m$ orthogonal matrix, and $V$ is an $n\times n$ orthogonal matrix. The diagonal elements of $\Sigma$ are the singular values of $A$; they are real and non-negative, and are returned in descending order. The first $min(m,n)$ columns of $U$ and $V$ are the left and right singular vectors of $A$.

Note that the function returns $V^T$, not $V$.

If the rank of $A$ is $r<min(m,n)$, then $\sigma$ has $r$ nonzero elements, and only $r$ columns of $U$ and $V$ are well-defined. In this case we can reduce $\sigma$ to an $r\times r$ matrix, $U$ to an $m\times r$ matrix and $V^T$ to an $r\times n$ matrix.

svd_rowext_real is designed for efficiently computing the SVD in the case $r<min(m,n)$. The input argument $k$ should be greater than $r$ by a small oversampling parameter, $\delta$, such that $k=r+\delta$. A reasonable value for $\delta$, to compute the SVD to within machine precision, is $\delta =5$. The value of $\delta$ should not vary based on $m$ or $n$. If $r$ is not known then the function can be used iteratively to refine the estimate and accuracy of the computed SVD; using a larger value of $\delta$ than necessary increases the computational cost of the function.

As a by-product of computing the SVD, the function estimates $r$.

If the input argument $k$ is less than $r$ the accuracy depends on the $(k+1)$th singular value, ${\sigma }_{k+1}$. See Accuracy for more details.

A call to svd_rowext_real consists of the following:

1. A random projection is applied, $Y=A\Omega$, where $\Omega$ is an $n\times k$ matrix. (Note that the product $A\Omega$ is computed using a Fast Fourier Transform, so can be computed in $O\left(mn\log \left(k\right)\right)$ time.) See randproj_dct_real() for more details on the random projection.

2. A pivoted $QR$ decomposition of $Y$ is calculated (see lapackeig.dgeqp3 for more details). The rank estimate is then $r$ such that, on the diagonal of $R$,

   $$\left|R_{rr}\right|>{ϵ}_a+{ϵ}_r\times R_{11}\text{,}$$

   where ${ϵ}_a$ and ${ϵ}_r$ are the absolute and relative error tolerances, respectively, and $r$ is the largest diagonal index for which the above relation holds.

3. Obtain the SVD from the $QR$ decomposition of $Y$ (or, depending on the rank, an approximation to the SVD) of $A$. This is referred to as row extraction.

Further details of the randomized SVD procedure can be found in Sections 4 and 5 of Halko et al. (2011).

References

Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, LAPACK Users’ Guide, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug

Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore

Halko, N, 2012, Randomized methods for computing low-rank approximations of matrices, PhD thesis

Halko, N, Martinsson, P G and Tropp, J A, 2011, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Rev. (53(2)), 217–288, https://epubs.siam.org/doi/abs/10.1137/090771806