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$k$, of the largest singular values and corresponding vectors of an m$m$ by n$n$ matrix A$A$. The value of k$k$ should be small relative to m$m$ and n$n$, for example kO(min (m,n))$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$m$ by n$n$ matrix A$A$ is given by
 A = UΣVT , $A=UΣVT ,$
where U$U$ and V$V$ are orthogonal and Σ$\Sigma$ is an m$m$ by n$n$ diagonal matrix with real diagonal elements, σi${\sigma }_{i}$, such that
 σ1 ≥ σ2 ≥ ⋯ ≥ σmin (m,n) ≥ 0 . $σ1 ≥ σ2 ≥⋯≥ σ min(m,n) ≥ 0 .$
The σi${\sigma }_{i}$ are the singular values of A$A$ and the first min (m,n)$\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$ columns of U$U$ and V$V$ are the left and right singular vectors of A$A$.
If Uk${U}_{k}$, Vk${V}_{k}$ denote the leading k$k$ columns of U$U$ and V$V$ respectively, and if Σk${\Sigma }_{k}$ denotes the leading principal submatrix of Σ $\Sigma$, then
 Ak ≡ Uk Σk VTk $Ak ≡ Uk Σk VTk$
is the best rank-k$k$ approximation to A$A$ in both the 2$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 = σi2 ui , $Avi = σi ui and ATui = σi vi so that ATA νi = σi2 νi ​ and ​ A AT ui = σ i 2 u i ,$
where ui${u}_{i}$ and vi${v}_{i}$ are the i$i$th columns of Uk${U}_{k}$ and Vk${V}_{k}$ respectively.
Thus, for mn$m\ge n$, the largest singular values and corresponding right singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix ATA${A}^{\mathrm{T}}A$. For m < n$m, the largest singular values and corresponding left singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix AAT$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$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$Ax$ or ATx${A}^{\mathrm{T}}x$ for a given real vector x$x$ (of length n$n$ or m$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:     m – int64int32nag_int scalar
m$m$, the number of rows of the matrix A$A$.
Constraint: m0${\mathbf{m}}\ge 0$.
If m = 0${\mathbf{m}}=0$, an immediate return is effected.
2:     n – int64int32nag_int scalar
n$n$, the number of columns of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
If n = 0${\mathbf{n}}=0$, an immediate return is effected.
3:     k – int64int32nag_int scalar
k$k$, the number of singular values to be computed.
Constraint: 0 < k < min (m,n)1$0<{\mathbf{k}}<\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)-1$.
4:     ncv – int64int32nag_int scalar
The dimension of the arrays sigma and resid and the second dimension of the arrays u and v as declared in the (sub)program from which nag_eigen_real_gen_partialsvd (f02wg) is called. This is the number of Lanczos basis vectors to use during the computation of the largest eigenvalues of ATA${A}^{\mathrm{T}}A$ (mn$m\ge n$) or AAT$A{A}^{\mathrm{T}}$ (m < n$m).
At present there is no a priori analysis to guide the selection of ncv relative to k. However, it is recommended that ncv2 × k + 1${\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: k < ncvmin (m,n)${\mathbf{k}}<{\mathbf{ncv}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
5:     av – function handle or string containing name of m-file
av must return the vector result of the matrix-vector product Ax$Ax$ or ATx${A}^{\mathrm{T}}x$, as indicated by the input value of iflag, for the given vector x$x$.
av is called from nag_eigen_real_gen_partialsvd (f02wg) with the parameter 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:     iflag – int64int32nag_int scalar
If iflag = 1${\mathbf{iflag}}=1$, ax must return the m$m$-vector result of the matrix-vector product Ax$Ax$.
If iflag = 2${\mathbf{iflag}}=2$, ax must return the n$n$-vector result of the matrix-vector product ATx${A}^{\mathrm{T}}x$.
2:     m – int64int32nag_int scalar
The number of rows of the matrix A$A$.
3:     n – int64int32nag_int scalar
The number of columns of the matrix A$A$.
4:     x( : $:$) – double array
The vector to be pre-multiplied by the matrix A$A$ or AT${A}^{\mathrm{T}}$.
5:     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:     iflag – int64int32nag_int scalar
May be used as a flag to indicate a failure in the computation of Ax$Ax$ or ATx${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:     ax( : $:$) – double array
If iflag = 1${\mathbf{iflag}}=1$, contains the m$m$-vector result of the matrix-vector product Ax$Ax$.
If iflag = 2${\mathbf{iflag}}=2$, contains the n$n$-vector result of the matrix-vector product ATx${A}^{\mathrm{T}}x$.
3:     user – Any MATLAB object

Optional Input Parameters

1:     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.

Input Parameters Omitted from the MATLAB Interface

ldu ldv iuser ruser

Output Parameters

1:     nconv – int64int32nag_int scalar
The number of converged singular values found.
2:     sigma(ncv) – double array
The nconv converged singular values are stored in the first nconv elements of sigma.
3:     u(ldu,ncv) – double array
ldumax (1,m)$\mathit{ldu}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The left singular vectors corresponding to the singular values stored in sigma.
The i$\mathit{i}$th element of the j$\mathit{j}$th left singular vector uj${u}_{\mathit{j}}$ is stored in u(i,j)${\mathbf{u}}\left(\mathit{i},\mathit{j}\right)$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$ and j = 1,2,,nconv$\mathit{j}=1,2,\dots ,{\mathbf{nconv}}$.
4:     v(ldv,ncv) – double array
ldvmax (1,n)$\mathit{ldv}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The right singular vectors corresponding to the singular values stored in sigma.
The i$\mathit{i}$th element of the j$\mathit{j}$th right singular vector vj${v}_{\mathit{j}}$ is stored in v(i,j)${\mathbf{v}}\left(\mathit{i},\mathit{j}\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$ and j = 1,2,,nconv$\mathit{j}=1,2,\dots ,{\mathbf{nconv}}$.
5:     resid(ncv) – double array
The residual Avjσjuj $‖A{v}_{j}-{\sigma }_{j}{u}_{j}‖$, for mn$m\ge n$, or ATujσjvj $‖{A}^{\mathrm{T}}{u}_{j}-{\sigma }_{j}{v}_{j}‖$, for m < n$m, for each of the converged singular values σj${\sigma }_{j}$ and corresponding left and right singular vectors uj${u}_{j}$ and vj${v}_{j}$.
6:     user – Any MATLAB object
7:     ifail – int64int32nag_int scalar
${\mathrm{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$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:
ifail < 0${\mathbf{ifail}}<0$
iflag was set to a negative value following a call to av.
ifail = 1${\mathbf{ifail}}=1$
 On entry, m < 0${\mathbf{m}}<0$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, n < 0${\mathbf{n}}<0$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, k ≤ 0${\mathbf{k}}\le 0$.
ifail = 4${\mathbf{ifail}}=4$
 On entry, ncv ≤ k + 1${\mathbf{ncv}}\le {\mathbf{k}}+1$, or ncv > min (m,n)${\mathbf{ncv}}>\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$.
ifail = 5${\mathbf{ifail}}=5$
 On entry, ldu < m$\mathit{ldu}<{\mathbf{m}}$.
ifail = 6${\mathbf{ifail}}=6$
 On entry, ldv < n$\mathit{ldv}<{\mathbf{n}}$.
ifail = 7${\mathbf{ifail}}=7$
No converged singular values were found to sufficient accuracy.
ifail = 8${\mathbf{ifail}}=8$
The internal maximum number of iterations has been reached. Some singular values may have converged.
ifail = 9${\mathbf{ifail}}=9$
ifail = 10${\mathbf{ifail}}=10$
ifail = 20${\mathbf{ifail}}=20$
An error occurred during an internal call. Consider increasing the size of ncv relative to k.

Accuracy

See Section [The singular value decomposition] in the F08 Chapter Introduction.

None.

Example

```function nag_eigen_real_gen_partialsvd_example
m = int64(100);
n = int64(500);
k = int64(4);
ncv = int64(10);
[nconv, sigma, u, v, resid, user, ifail] = nag_eigen_real_gen_partialsvd(m, n, k, ncv, @av);
nconv, sigma, resid, user, ifail

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:double(n)
t = t + k;
s = 0;
for i = 1:double(min(j,m))
s = s + h;
ax(i) = ax(i) + k*s*(t-1)*x(j);
end
for i = double(j+1):double(m)
s = s + h;
ax(i) = ax(i) + k*t*(s-1)*x(j);
end
end
else
t = 0;
for j = 1:double(n)
t = t + k;
s = 0;
for i = 1:double(min(j,m))
s = s + h;
ax(j) = ax(j) + k*s*(t-1)*x(i);
end
for i = double(j + 1):double(m)
s = s + h;
ax(j) = ax(j) + k*t*(s-1)*x(i);
end
end
end
```
```

nconv =

4

sigma =

0.0083
0.0122
0.0238
0.1127
0
0
0
0
0
0

resid =

1.0e-16 *

0.0154
0.0422
0.0943
0.3472
0
0
0
0
0
0

user =

0

ifail =

0

```
```function f02wg_example
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);
nconv, sigma, resid, user, ifail

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:double(n)
t = t + k;
s = 0;
for i = 1:double(min(j,m))
s = s + h;
ax(i) = ax(i) + k*s*(t-1)*x(j);
end
for i = double(j+1):double(m)
s = s + h;
ax(i) = ax(i) + k*t*(s-1)*x(j);
end
end
else
t = 0;
for j = 1:double(n)
t = t + k;
s = 0;
for i = 1:double(min(j,m))
s = s + h;
ax(j) = ax(j) + k*s*(t-1)*x(i);
end
for i = double(j + 1):double(m)
s = s + h;
ax(j) = ax(j) + k*t*(s-1)*x(i);
end
end
end
```
```

nconv =

4

sigma =

0.0083
0.0122
0.0238
0.1127
0
0
0
0
0
0

resid =

1.0e-16 *

0.0154
0.0422
0.0943
0.3472
0
0
0
0
0
0

user =

0

ifail =

0

```