hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
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, kk, of the largest singular values and corresponding vectors of an mm by nn matrix AA. The value of kk should be small relative to mm and nn, for example kO(min (m,n))kO(min(m,n)). The full singular value decomposition (SVD) of an mm by nn matrix AA is given by
A = UΣVT ,
A=UΣVT ,
where UU and VV are orthogonal and ΣΣ is an mm by nn diagonal matrix with real diagonal elements, σiσi, such that
σ1 σ2 σmin (m,n) 0 .
σ1 σ2 σ min(m,n) 0 .
The σiσi are the singular values of AA and the first min (m,n)min(m,n) columns of UU and VV are the left and right singular vectors of AA.
If UkUk, VkVk denote the leading kk columns of UU and VV respectively, and if ΣkΣk denotes the leading principal submatrix of Σ Σ, then
Ak Uk Σk VTk
Ak Uk Σk VTk
is the best rank-kk approximation to AA in both the 22-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 uiui and vivi are the iith columns of UkUk and VkVk respectively.
Thus, for mnmn, the largest singular values and corresponding right singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix ATAATA. For m < nm<n, the largest singular values and corresponding left singular vectors are computed by finding eigenvalues and eigenvectors for the symmetric matrix AATAAT. 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 AA 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 AxAx or ATxATx for a given real vector xx (of length nn or mm 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
mm, the number of rows of the matrix AA.
Constraint: m0m0.
If m = 0m=0, an immediate return is effected.
2:     n – int64int32nag_int scalar
nn, the number of columns of the matrix AA.
Constraint: n0n0.
If n = 0n=0, an immediate return is effected.
3:     k – int64int32nag_int scalar
kk, the number of singular values to be computed.
Constraint: 0 < k < min (m,n)10<k<min(m,n)-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 ATAATA (mnmn) or AATAAT (m < nm<n).
At present there is no a priori analysis to guide the selection of ncv relative to k. However, it is recommended that ncv2 × k + 1ncv2×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)k<ncvmin(m,n).
5:     av – function handle or string containing name of m-file
av must return the vector result of the matrix-vector product AxAx or ATxATx, as indicated by the input value of iflag, for the given vector xx.
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 = 1iflag=1, ax must return the mm-vector result of the matrix-vector product AxAx.
If iflag = 2iflag=2, ax must return the nn-vector result of the matrix-vector product ATxATx.
2:     m – int64int32nag_int scalar
The number of rows of the matrix AA.
3:     n – int64int32nag_int scalar
The number of columns of the matrix AA.
4:     x( : :) – double array
The vector to be pre-multiplied by the matrix AA or ATAT.
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 AxAx or ATxATx. 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 = 1iflag=1, contains the mm-vector result of the matrix-vector product AxAx.
If iflag = 2iflag=2, contains the nn-vector result of the matrix-vector product ATxATx.
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)ldumax(1,m).
The left singular vectors corresponding to the singular values stored in sigma.
The iith element of the jjth left singular vector ujuj is stored in u(i,j)uij, for i = 1,2,,mi=1,2,,m and j = 1,2,,nconvj=1,2,,nconv.
4:     v(ldv,ncv) – double array
ldvmax (1,n)ldvmax(1,n).
The right singular vectors corresponding to the singular values stored in sigma.
The iith element of the jjth right singular vector vjvj is stored in v(i,j)vij, for i = 1,2,,ni=1,2,,n and j = 1,2,,nconvj=1,2,,nconv.
5:     resid(ncv) – double array
The residual Avjσjuj A vj - σj uj , for mnmn, or ATujσjvj AT uj - σj vj , for m < nm<n, for each of the converged singular values σjσj and corresponding left and right singular vectors ujuj and vjvj.
6:     user – Any MATLAB object
7:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).
nag_eigen_real_gen_partialsvd (f02wg) returns with ifail = 0ifail=0 if at least kk 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 < 0ifail<0
iflag was set to a negative value following a call to av.
  ifail = 1ifail=1
On entry,m < 0m<0.
  ifail = 2ifail=2
On entry,n < 0n<0.
  ifail = 3ifail=3
On entry,k0k0.
  ifail = 4ifail=4
On entry,ncvk + 1ncvk+1,
orncv > min (m,n)ncv>min(m,n).
  ifail = 5ifail=5
On entry,ldu < mldu<m.
  ifail = 6ifail=6
On entry,ldv < nldv<n.
  ifail = 7ifail=7
No converged singular values were found to sufficient accuracy.
  ifail = 8ifail=8
The internal maximum number of iterations has been reached. Some singular values may have converged.
  ifail = 9ifail=9
  ifail = 10ifail=10
  ifail = 20ifail=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.

Further Comments

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



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013