Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_eigen_real_gen_partialsvd (f02wg)

## Purpose

nag_eigen_real_gen_partialsvd (f02wg) returns leading terms in the singular value decomposition (SVD) of a real general matrix and computes the corresponding left and right singular vectors.

## Syntax

[nconv, sigma, u, v, resid, user, ifail] = f02wg(m, n, k, ncv, av, 'user', user)
[nconv, sigma, u, v, resid, user, ifail] = nag_eigen_real_gen_partialsvd(m, n, k, ncv, av, 'user', user)

## Description

nag_eigen_real_gen_partialsvd (f02wg) computes a few, $k$, of the largest singular values and corresponding vectors of an $m$ by $n$ matrix $A$. The value of $k$ should be small relative to $m$ and $n$, for example $k\sim O\left(\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)\right)$. The full singular value decomposition (SVD) of an $m$ by $n$ matrix $A$ is given by
 $A=UΣVT ,$
where $U$ and $V$ are orthogonal and $\Sigma$ is an $m$ by $n$ diagonal matrix with real diagonal elements, ${\sigma }_{i}$, such that
 $σ1 ≥ σ2 ≥⋯≥ σ minm,n ≥ 0 .$
The ${\sigma }_{i}$ are the singular values of $A$ and 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$.
If ${U}_{k}$, ${V}_{k}$ denote the leading $k$ columns of $U$ and $V$ respectively, and if ${\Sigma }_{k}$ denotes the leading principal submatrix of $\Sigma$, then
 $Ak ≡ Uk Σk VTk$
is the best rank-$k$ approximation to $A$ in both the $2$-norm and the Frobenius norm.
The singular values and singular vectors satisfy
 $Avi = σi ui and ATui = σi vi so that ATA νi = σi2 νi ​ and ​ A AT ui = σ i 2 u i ,$
where ${u}_{i}$ and ${v}_{i}$ are the $i$th columns of ${U}_{k}$ and ${V}_{k}$ respectively.
Thus, for $m\ge n$, the largest singular values and corresponding right singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix ${A}^{\mathrm{T}}A$. For $m, the largest singular values and corresponding left singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix $A{A}^{\mathrm{T}}$. These eigenvalues and eigenvectors are found using functions from Chapter F12. You should read the F12 Chapter Introduction for full details of the method used here.
The real matrix $A$ is not explicitly supplied to nag_eigen_real_gen_partialsvd (f02wg). Instead, you are required to supply a function, av, that must calculate one of the requested matrix-vector products $Ax$ or ${A}^{\mathrm{T}}x$ for a given real vector $x$ (of length $n$ or $m$ respectively).

## References

Wilkinson J H (1978) Singular Value Decomposition – Basic Aspects Numerical Software – Needs and Availability (ed D A H Jacobs) Academic Press

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{m}$int64int32nag_int scalar
$m$, the number of rows of the matrix $A$.
Constraint: ${\mathbf{m}}\ge 0$.
If ${\mathbf{m}}=0$, an immediate return is effected.
2:     $\mathrm{n}$int64int32nag_int scalar
$n$, the number of columns of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
If ${\mathbf{n}}=0$, an immediate return is effected.
3:     $\mathrm{k}$int64int32nag_int scalar
$k$, the number of singular values to be computed.
Constraint: $0<{\mathbf{k}}<\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)-1$.
4:     $\mathrm{ncv}$int64int32nag_int scalar
The dimension of the arrays sigma and resid and the second dimension of the arrays u and v. this is the number of Lanczos basis vectors to use during the computation of the largest eigenvalues of ${A}^{\mathrm{T}}A$ ($m\ge n$) or $A{A}^{\mathrm{T}}$ ($m).
At present there is no a priori analysis to guide the selection of ncv relative to k. However, it is recommended that ${\mathbf{ncv}}\ge 2×{\mathbf{k}}+1$. If many problems of the same type are to be solved, you should experiment with varying ncv while keeping k fixed for a given test problem. This will usually decrease the required number of matrix-vector operations but it also increases the internal storage required to maintain the orthogonal basis vectors. The optimal ‘cross-over’ with respect to CPU time is problem dependent and must be determined empirically.
Constraint: ${\mathbf{k}}<{\mathbf{ncv}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
5:     $\mathrm{av}$ – function handle or string containing name of m-file
av must return the vector result of the matrix-vector product $Ax$ or ${A}^{\mathrm{T}}x$, as indicated by the input value of iflag, for the given vector $x$.
av is called from nag_eigen_real_gen_partialsvd (f02wg) with the argument user as supplied to nag_eigen_real_gen_partialsvd (f02wg). You are free to use these arrays to supply information to av.
[iflag, ax, user] = av(iflag, m, n, x, user)

Input Parameters

1:     $\mathrm{iflag}$int64int32nag_int scalar
If ${\mathbf{iflag}}=1$, ax must return the $m$-vector result of the matrix-vector product $Ax$.
If ${\mathbf{iflag}}=2$, ax must return the $n$-vector result of the matrix-vector product ${A}^{\mathrm{T}}x$.
2:     $\mathrm{m}$int64int32nag_int scalar
The number of rows of the matrix $A$.
3:     $\mathrm{n}$int64int32nag_int scalar
The number of columns of the matrix $A$.
4:     $\mathrm{x}\left(:\right)$ – double array
The vector to be pre-multiplied by the matrix $A$ or ${A}^{\mathrm{T}}$.
5:     $\mathrm{user}$ – Any MATLAB object
av is called from nag_eigen_real_gen_partialsvd (f02wg) with the object supplied to nag_eigen_real_gen_partialsvd (f02wg).

Output Parameters

1:     $\mathrm{iflag}$int64int32nag_int scalar
May be used as a flag to indicate a failure in the computation of $Ax$ or ${A}^{\mathrm{T}}x$. If iflag is negative on exit from av, nag_eigen_real_gen_partialsvd (f02wg) will exit immediately with ifail set to iflag.
2:     $\mathrm{ax}\left(:\right)$ – double array
If ${\mathbf{iflag}}=1$, contains the $m$-vector result of the matrix-vector product $Ax$.
If ${\mathbf{iflag}}=2$, contains the $n$-vector result of the matrix-vector product ${A}^{\mathrm{T}}x$.
3:     $\mathrm{user}$ – Any MATLAB object

### Optional Input Parameters

1:     $\mathrm{user}$ – Any MATLAB object
user is not used by nag_eigen_real_gen_partialsvd (f02wg), but is passed to av. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

### Output Parameters

1:     $\mathrm{nconv}$int64int32nag_int scalar
The number of converged singular values found.
2:     $\mathrm{sigma}\left({\mathbf{ncv}}\right)$ – double array
The nconv converged singular values are stored in the first nconv elements of sigma.
3:     $\mathrm{u}\left(\mathit{ldu},{\mathbf{ncv}}\right)$ – double array
The left singular vectors corresponding to the singular values stored in sigma.
The $\mathit{i}$th element of the $\mathit{j}$th left singular vector ${u}_{\mathit{j}}$ is stored in ${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$, for $\mathit{i}=1,2,\dots ,m$ and $\mathit{j}=1,2,\dots ,{\mathbf{nconv}}$.
4:     $\mathrm{v}\left(\mathit{ldv},{\mathbf{ncv}}\right)$ – double array
The right singular vectors corresponding to the singular values stored in sigma.
The $\mathit{i}$th element of the $\mathit{j}$th right singular vector ${v}_{\mathit{j}}$ is stored in ${\mathbf{v}}\left(\mathit{i},\mathit{j}\right)$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,{\mathbf{nconv}}$.
5:     $\mathrm{resid}\left({\mathbf{ncv}}\right)$ – double array
The residual $‖A{v}_{j}-{\sigma }_{j}{u}_{j}‖$, for $m\ge n$, or $‖{A}^{\mathrm{T}}{u}_{j}-{\sigma }_{j}{v}_{j}‖$, for $m, for each of the converged singular values ${\sigma }_{j}$ and corresponding left and right singular vectors ${u}_{j}$ and ${v}_{j}$.
6:     $\mathrm{user}$ – Any MATLAB object
7:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).
nag_eigen_real_gen_partialsvd (f02wg) returns with ${\mathbf{ifail}}={\mathbf{0}}$ if at least $k$ singular values have converged and the corresponding left and right singular vectors have been computed.

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
Constraint: ${\mathbf{m}}\ge 0$.
${\mathbf{ifail}}=2$
Constraint: ${\mathbf{n}}\ge 0$.
${\mathbf{ifail}}=3$
Constraint: ${\mathbf{k}}>0$.
${\mathbf{ifail}}=4$
Constraint: ${\mathbf{k}}<{\mathbf{ncv}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
${\mathbf{ifail}}=5$
Constraint: $\mathit{ldu}\ge {\mathbf{m}}$.
${\mathbf{ifail}}=6$
Constraint: $\mathit{ldv}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=8$
The maximum number of iterations has been reached.
${\mathbf{ifail}}=9$
No shifts could be applied during a cycle of the implicitly restarted Lanczos iteration.
${\mathbf{ifail}}=10$
Could not build a full Lanczos factorization.
${\mathbf{ifail}}=11$
The number of eigenvalues found to sufficient accuracy is zero.
${\mathbf{ifail}}=20$
An error occurred during an internal call. Consider increasing the size of ncv relative to k.
${\mathbf{ifail}}<0$
On output from user-defined function av, iflag was set to a negative value, .
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

See The singular value decomposition in the F08 Chapter Introduction.

None.

## Example

This example finds the four largest singular values ($\sigma$) and corresponding right and left singular vectors for the matrix $A$, where $A$ is the $m$ by $n$ real matrix derived from the simplest finite difference discretization of the two-dimensional kernal $k\left(s,t\right)dt$ where
 $ks,t = st-1 if ​0≤s≤t≤ 1 ts-1 if ​0≤t
```function f02wg_example

fprintf('f02wg example results\n\n');

m = int64(100);
n = int64(500);
k = int64(4);
ncv = int64(10);

[nconv, sigma, u, v, resid, user, ifail] = ...
f02wg(m, n, k, ncv, @av);

res = [sigma(1:nconv) resid(1:nconv)];
fprintf('  Singular Value   Residual\n');
fprintf('   %10.5f  %12.2e\n',res');

function [iflag, ax, user] = av(iflag, m, n, x, user)
if (iflag == 1)
ax = zeros(m, 1);
else
ax = zeros(n, 1);
end

% Computes  w <- A*x or w <- Trans(A)*x.
h = 1/(double(m)+1);
k = 1/(double(n)+1);
if (iflag == 1)
t = 0;
for j = 1:n
t = t + k;
s = 0;
ktx = k*(t-1)*x(j);
for i = 1:min(j,m)
s = s + h;
ax(i) = ax(i) + ktx*s;
end
ktx = k*t*x(j);
for i = j+1:m
s = s + h;
ax(i) = ax(i) + ktx*(s-1);
end
end
else
t = 0;
for j = 1:n
t = t + k;
s = 0;
for i = 1:min(j,m)
s = s + h;
ax(j) = ax(j) + k*s*(t-1)*x(i);
end
for i = j+1:m
s = s + h;
ax(j) = ax(j) + k*t*(s-1)*x(i);
end
end
end
```
```f02wg example results

Singular Value   Residual
0.00830      1.82e-18
0.01223      3.95e-18
0.02381      1.95e-17
0.11274      2.53e-17
```