# NAG CL Interfacef10cac (svd_​rowext_​real)

## 1Purpose

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

## 2Specification

 #include
 void f10cac (Nag_ComputeUType jobu, Nag_ComputeVTType jobvt, Integer m, Integer n, const double a[], Integer pda, Integer k, double rtol_abs, double rtol_rel, Integer state[], double s[], double u[], Integer pdu, double vt[], Integer pdvt, Integer *r, NagError *fail)
The function may be called by the names: f10cac or nag_rnla_svd_rowext_real.

## 3Description

The SVD is written as
 $A = UΣVT ,$
where $\Sigma$ is an $m$ by $n$ matrix which is zero except for its $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$ diagonal elements, $U$ is an $m$ by $m$ orthogonal matrix, and $V$ is an $n$ by $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 $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$ columns of $U$ and $V$ are the left and right singular vectors of $A$.
Note that the function returns ${V}^{\mathrm{T}}$, not $V$.
If the rank of $A$ is $r<\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$, 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$ by $r$ matrix, $U$ to an $m$ by $r$ matrix and ${V}^{\mathrm{T}}$ to an $r$ by $n$ matrix.
f10cac is designed for efficiently computing the SVD in the case $r<\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$. 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 $\left(k+1\right)$th singular value, ${\sigma }_{k+1}$. See Section 7 for more details.
A call to f10cac consists of the following:
1. 1.A random projection is applied, $Y=A\Omega$, where $\Omega$ is an $n$ by $k$ matrix. (Note that the product $A\Omega$ is computed using a Fast Fourier Transform, so can be computed in $\mathit{O}\left(mn\mathrm{log}\left(k\right)\right)$ time.) See f10dac for more details on the random projection.
2. 2.A pivoted $QR$ decomposition of $Y$ is calculated (see f08bec for more details). The rank estimate is then $r$ such that, on the diagonal of $R$,
 $R rr > εa + εr × R 11 ,$
where ${\epsilon }_{a}$ and ${\epsilon }_{r}$ are the absolute and relative error tolerances, respectively, and $r$ is the largest diagonal index for which the above relation holds.
3. 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).

## 4References

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

## 5Arguments

1: $\mathbf{jobu}$Nag_ComputeUType Input
On entry: specifies options for computing part of or none of the matrix $U$.
${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$
The first $k$ columns of $U$ (the left singular vectors) are returned in the array u.
${\mathbf{jobu}}=\mathrm{Nag_NotU}$
No columns of $U$ (no left singular vectors) are computed.
Constraint: ${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$ or $\mathrm{Nag_NotU}$.
2: $\mathbf{jobvt}$Nag_ComputeVTType Input
On entry: specifies options for computing part of or none of the matrix ${V}^{\mathrm{T}}$.
${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$
The first $k$ rows of ${V}^{\mathrm{T}}$ (the right singular vectors) are returned in the array vt.
${\mathbf{jobvt}}=\mathrm{Nag_NotVT}$
No rows of ${V}^{\mathrm{T}}$ (no right singular vectors) are computed.
Constraint: ${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$ or $\mathrm{Nag_NotVT}$.
3: $\mathbf{m}$Integer Input
On entry: $m$, the number of rows of the matrix $A$.
Constraint: ${\mathbf{m}}>0$.
4: $\mathbf{n}$Integer Input
On entry: $n$, the number of columns of the matrix $A$.
Constraint: ${\mathbf{n}}>0$.
5: $\mathbf{a}\left[\mathit{dim}\right]$const double Input
Note: the dimension, dim, of the array a must be at least ${\mathbf{pda}}×{\mathbf{n}}$.
The $\left(i,j\right)$th element of the matrix $A$ is stored in ${\mathbf{a}}\left[\left(j-1\right)×{\mathbf{pda}}+i-1\right]$.
On entry: the $m$ by $n$ matrix $A$.
6: $\mathbf{pda}$Integer Input
On entry: the stride separating matrix row elements in the array a.
Constraint: ${\mathbf{pda}}\ge {\mathbf{m}}$.
7: $\mathbf{k}$Integer Input
On entry: $k$, number of columns in random projection, $Y=A\Omega$.
Constraint: $0<{\mathbf{k}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
8: $\mathbf{rtol_abs}$double Input
On entry: the absolute tolerance, ${\epsilon }_{a}$ used in defining the threshold on estimating the rank of $A$. If ${\mathbf{rtol_abs}}<0.0$ then $0.0$ is used unless ${\mathbf{rtol_rel}}\le 0.0$ in which case is used.
9: $\mathbf{rtol_rel}$double Input
On entry: the relative tolerance, ${\epsilon }_{r}$ used in defining the threshold on estimating the rank of $A$. If ${\mathbf{rtol_rel}}<0.0$ then $0.0$ is used unless ${\mathbf{rtol_abs}}\le 0.0$ in which case is used.
10: $\mathbf{state}\left[\mathit{dim}\right]$Integer Communication Array
Note: the dimension, $\mathit{dim}$, of this array is dictated by the requirements of associated functions that must have been previously called. This array MUST be the same array passed as argument state in the previous call to nag_rand_init_repeatable (g05kfc) or nag_rand_init_nonrepeatable (g05kgc).
On entry: contains information on the selected base generator and its current state.
On exit: contains updated information on the state of the generator.
11: $\mathbf{s}\left[{\mathbf{k}}\right]$double Output
On exit: the first r elements of s contain the r largest singular values of $A$ in descending order. The remaining values are set to zero.
12: $\mathbf{u}\left[\mathit{dim}\right]$double Output
Note: the dimension, dim, of the array u must be at least ${\mathbf{k}}$ when ${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$.
The $\left(i,j\right)$th element of the matrix $U$ is stored in ${\mathbf{u}}\left[\left(j-1\right)×{\mathbf{pdu}}+i-1\right]$.
On exit: if ${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$, u contains the first r columns of $U$ (the left singular vectors, stored column-wise); the remaining elements of u are set to zero.
If ${\mathbf{jobu}}=\mathrm{Nag_NotU}$, u is not referenced.
13: $\mathbf{pdu}$Integer Input
On entry: the stride separating matrix row elements in the array u.
Constraint: if ${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$, ${\mathbf{pdu}}\ge {\mathbf{m}}$
14: $\mathbf{vt}\left[\mathit{dim}\right]$double Output
Note: the dimension, dim, of the array vt must be at least
• ${\mathbf{pdvt}}×{\mathbf{n}}$ when ${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$;
• $1$ otherwise.
The $\left(i,j\right)$th element of the matrix is stored in ${\mathbf{vt}}\left[\left(j-1\right)×{\mathbf{pdvt}}+i-1\right]$.
On exit: if ${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$, vt contains the first r rows of ${V}^{\mathrm{T}}$ (the right singular vectors); the remaining elements of vt are set to zero.
If ${\mathbf{jobvt}}=\mathrm{Nag_NotVT}$, vt is not referenced.
15: $\mathbf{pdvt}$Integer Input
On entry: the stride separating matrix row elements in the array vt.
Constraint: if ${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$, ${\mathbf{pdvt}}\ge {\mathbf{k}}$
16: $\mathbf{r}$Integer * Output
On exit: $r$, contains estimated rank of array a.
17: $\mathbf{fail}$NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

## 6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
On entry, argument $〈\mathit{\text{value}}〉$ had an illegal value.
NE_INT
On entry, ${\mathbf{k}}=〈\mathit{\text{value}}〉$.
Constraint: $0<{\mathbf{k}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
On entry, ${\mathbf{m}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{m}}>0$.
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}>0$.
On entry, ${\mathbf{pda}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{pda}}\ge {\mathbf{m}}$.
On entry, state vector has been corrupted or not initialized.
NE_INT_2
On entry, ${\mathbf{jobu}}=〈\mathit{\text{value}}〉$ and ${\mathbf{pdu}}=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{jobu}}=\mathrm{Nag_SingularVecsU}$, ${\mathbf{pdu}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{jobvt}}=〈\mathit{\text{value}}〉$ and ${\mathbf{pdvt}}=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{jobvt}}=\mathrm{Nag_SingularVecsVT}$, ${\mathbf{pdvt}}\ge {\mathbf{k}}$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_REAL
$A$ has effective rank of zero.
First diagonal element of $R$, from $QR$ of $Y$, $=$$〈\mathit{\text{value}}〉$.
Tolerance used to determine rank $=$$〈\mathit{\text{value}}〉$.
On exit, ${\mathbf{r}}={\mathbf{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$, $=$$〈\mathit{\text{value}}〉$.
Tolerance used to determine rank $=$$〈\mathit{\text{value}}〉$.

## 7Accuracy

The error is approximately,
 $A-UσV* < 1 + 1 + 4k n-k ε ,$
where,
 $ε= 1+7n/k · σk+1 .$
The norm on the left-hand side of the first equation is the spectral norm, and ${\sigma }_{k+1}$ is the $\left(k+1\right)$th singular value of $A$. More details on the error bound can be found in Sections 5 and 11 of Halko et al. (2011).

## 8Parallelism and Performance

f10cac is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f10cac makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

The total number of floating-point operations is $\mathit{O}\left(\mathrm{log}\left(k\right)mn+{k}^{2}\left(m+n\right)\right)$. The first term corresponds to applying the random projection, i.e., computing $Y=A\Omega$. The second term corresponds to the $QR$ decomposition of $Y$ and the steps required to obtain the SVD of the original matrix $A$.
Deterministic SVD solvers, such as f08kbc, require $\mathit{O}\left(m{n}^{2}\right)$ operations when $m\ge n$ and $\mathit{O}\left(n{m}^{2}\right)$ operations when $n\ge m$.
The default values for rtol_abs and rtol_rel assume that you need an accurate approximation to $A$. If you only need to use a small number of singular values or singular vectors, larger values for these tolerances are appropriate. Increasing tolerances sufficiently will decrease r, the estimated rank. Decreasing r means that k can then be decreased to reduce the run-time of the routine.

## 10Example

This example finds the singular values, the left and right singular vectors, and the rank of the $6$ by $6$ matrix
 $A = 2.27 0.28 -0.48 1.07 -2.35 0.62 -1.54 -1.67 -3.09 1.22 2.93 -7.39 1.15 0.94 0.99 0.79 -1.45 1.03 -1.94 -0.78 -0.21 0.63 2.30 -2.57 -1.94 -0.78 -0.21 0.63 2.30 -2.57 -1.94 -0.78 -0.21 0.63 2.30 -2.57 ,$
using the randomised solver, f10cac, and a deterministic solver, f08kbc for comparison.

### 10.1Program Text

Program Text (f10cace.c)

### 10.2Program Data

Program Data (f10cace.d)

### 10.3Program Results

Program Results (f10cace.r)