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_complex_triang_svd (f02xu)

Purpose

nag_eigen_complex_triang_svd (f02xu) returns all, or part, of the singular value decomposition of a complex upper triangular matrix.

Syntax

[a, b, q, sv, rwork, ifail] = f02xu(a, b, wantq, wantp, 'n', n, 'ncolb', ncolb)
[a, b, q, sv, rwork, ifail] = nag_eigen_complex_triang_svd(a, b, wantq, wantp, 'n', n, 'ncolb', ncolb)

Description

The nn by nn upper triangular matrix RR is factorized as
R = QSPH,
R=QSPH,
where QQ and PP are nn by nn unitary matrices and SS is an nn by nn diagonal matrix with real non-negative diagonal elements, sv1,sv2,,svnsv1,sv2,,svn, ordered such that
sv1sv2svn0.
sv1sv2svn0.
The columns of QQ are the left-hand singular vectors of RR, the diagonal elements of SS are the singular values of RR and the columns of PP are the right-hand singular vectors of RR.
Either or both of QQ and PHPH may be requested and the matrix CC given by
C = QHB,
C=QHB,
where BB is an nn by ncolbncolb given matrix, may also be requested.
nag_eigen_complex_triang_svd (f02xu) obtains the singular value decomposition by first reducing RR to bidiagonal form by means of Givens plane rotations and then using the QRQR algorithm to obtain the singular value decomposition of the bidiagonal form.
Good background descriptions to the singular value decomposition are given in Dongarra et al. (1979), Hammarling (1985) and Wilkinson (1978).
Note that if KK is any unitary diagonal matrix so that
KKH = I,
KKH=I,
then
A = (QK)S(PK)H
A=(QK)S(PK)H
is also a singular value decomposition of AA.

References

Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979) LINPACK Users' Guide SIAM, Philadelphia
Hammarling S (1985) The singular value decomposition in multivariate statistics SIGNUM Newsl. 20(3) 2–25
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:     a(lda, : :) – complex array
The first dimension of the array a must be at least max (1,n)max(1,n)
The second dimension of the array must be at least max (1,n)max(1,n)
The leading nn by nn upper triangular part of the array a must contain the upper triangular matrix RR.
2:     b(ldb, : :) – complex array
The first dimension, ldb, of the array b must satisfy
  • if ncolb > 0ncolb>0, ldbmax (1,n)ldbmax(1,n);
  • otherwise ldb1ldb1.
The second dimension of the array must be at least max (1,ncolb)max(1,ncolb)
If ncolb > 0ncolb>0, the leading nn by ncolbncolb part of the array b must contain the matrix to be transformed.
3:     wantq – logical scalar
Must be true if the matrix QQ is required.
If wantq = falsewantq=false then the array q is not referenced.
4:     wantp – logical scalar
Must be true if the matrix PHPH is required, in which case PHPH is returned in the array a, otherwise wantp must be false.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array a The second dimension of the array a.
nn, the order of the matrix RR.
If n = 0n=0, an immediate return is effected.
Constraint: n0n0.
2:     ncolb – int64int32nag_int scalar
Default: The second dimension of the array b.
ncolbncolb, the number of columns of the matrix BB.
If ncolb = 0ncolb=0, the array b is not referenced.
Constraint: ncolb0ncolb0.

Input Parameters Omitted from the MATLAB Interface

lda ldb ldq cwork

Output Parameters

1:     a(lda, : :) – complex array
The first dimension of the array a will be max (1,n)max(1,n)
The second dimension of the array will be max (1,n)max(1,n)
ldamax (1,n)ldamax(1,n).
If wantp = truewantp=true, the nn by nn part of a will contain the nn by nn unitary matrix PHPH, otherwise the nn by nn upper triangular part of a is used as internal workspace, but the strictly lower triangular part of a is not referenced.
2:     b(ldb, : :) – complex array
The first dimension, ldb, of the array b will be
  • if ncolb > 0ncolb>0, ldbmax (1,n)ldbmax(1,n);
  • otherwise ldb1ldb1.
The second dimension of the array will be max (1,ncolb)max(1,ncolb)
Stores the nn by ncolbncolb matrix QHBQHB.
3:     q(ldq, : :) – complex array
The first dimension, ldq, of the array q will be
  • if wantq = truewantq=true, ldqmax (1,n)ldqmax(1,n);
  • otherwise ldq1ldq1.
The second dimension of the array will be max (1,n)max(1,n) if wantq = truewantq=true, and at least 11 otherwise
If wantq = truewantq=true, the leading nn by nn part of the array q will contain the unitary matrix QQ. Otherwise the array q is not referenced.
4:     sv(n) – double array
The nn diagonal elements of the matrix SS.
5:     rwork( : :) – double array
Note: the dimension of the array rwork must be at least max (1,2 × (n1))max(1,2×(n-1)) if ncolb = 0ncolb=0 and wantq = falsewantq=false and wantp = falsewantp=false, max (1,3 × (n1))max(1,3×(n-1)) if ncolb = 0ncolb=0 and wantq = falsewantq=false and wantp = truewantp=true or ncolb > 0ncolb>0 and wantp = falsewantp=false or wantq = truewantq=true and wantp = falsewantp=false, and at least max (1,5 × (n1))max(1,5×(n-1)) otherwise.
rwork(n) contains the total number of iterations taken by the QRQR algorithm.
The rest of the array is used as workspace.
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  ifail = 1ifail=-1
On entry,n < 0n<0,
orlda < nlda<n,
orncolb < 0ncolb<0,
orldb < nldb<n and ncolb > 0ncolb>0,
orldq < nldq<n and wantq = truewantq=true
W ifail > 0ifail>0
The QRQR algorithm has failed to converge in 50 × n50×n iterations. In this case sv(1),sv(2),,sv(ifail)sv1,sv2,,svifail may not have been found correctly and the remaining singular values may not be the smallest. The matrix RR will nevertheless have been factorized as R = QEPHR=QEPH, where EE is a bidiagonal matrix with sv(1),sv(2),,sv(n)sv1,sv2,,svn as the diagonal elements and rwork(1),rwork(2),,rwork(n1)rwork1,rwork2,,rworkn-1 as the superdiagonal elements.
This failure is not likely to occur.

Accuracy

The computed factors QQ, SS and PP satisfy the relation
QSPH = A + E,
QSPH=A+E,
where
Ecε A,
Ecε A,
εε is the machine precision, cc is a modest function of nn and . . denotes the spectral (two) norm. Note that A = sv1A=sv1.

Further Comments

For given values of ncolb, wantq and wantp, the number of floating point operations required is approximately proportional to n3n3.
>Following the use of nag_eigen_complex_triang_svd (f02xu) the rank of RR may be estimated as follows:
tol = eps;
irank = 1;
while (irank <= numel(sv) && sv(irank) >= tol*sv(1) )
  irank = irank + 1;
end
returns the value kk in irank, where kk is the smallest integer for which sv(k) < tol × sv(1)svk<tol×sv1, where toltol is typically the machine precision, so that irank is an estimate of the rank of SS and thus also of RR.

Example

function nag_eigen_complex_triang_svd_example
a = [1,  1 + 1i,  1 + 1i;
      0 + 0i,  -2 + 0i,  -1 - 1i;
      0 + 0i,  0 + 0i,  -3 + 0i];
b = [ 1 + 1i;  -1 + 0i;  -1 + 1i];
wantq = true;
wantp = true;
[aOut, bOut, q, sv, rwork, ifail] = nag_eigen_complex_triang_svd(a, b, wantq, wantp)
 

aOut =

   0.1275 + 0.0000i   0.3899 + 0.2046i   0.5289 + 0.7142i
   0.2265 + 0.0000i   0.3397 + 0.7926i   0.0000 - 0.4529i
   0.9656 + 0.0000i  -0.1311 - 0.2129i  -0.0698 + 0.0119i


bOut =

   1.9656 + 0.7935i
  -0.1132 + 0.3397i
   0.0915 + 0.6086i


q =

   0.5005 + 0.0000i   0.4529 + 0.0000i   0.7378 + 0.0000i
  -0.5152 + 0.1514i  -0.1132 + 0.5661i   0.4190 - 0.4502i
  -0.4041 + 0.5457i  -0.0000 - 0.6794i   0.2741 + 0.0468i


sv =

    3.9263
    2.0000
    0.7641


rwork =

         0
         0
    1.0000
    0.7071
         0
         0
         0
         0
         0
         0


ifail =

                    0


function f02xu_example
a = [1,  1 + 1i,  1 + 1i;
      0 + 0i,  -2 + 0i,  -1 - 1i;
      0 + 0i,  0 + 0i,  -3 + 0i];
b = [ 1 + 1i;  -1 + 0i;  -1 + 1i];
wantq = true;
wantp = true;
[aOut, bOut, q, sv, rwork, ifail] = f02xu(a, b, wantq, wantp)
 

aOut =

   0.1275 + 0.0000i   0.3899 + 0.2046i   0.5289 + 0.7142i
   0.2265 + 0.0000i   0.3397 + 0.7926i   0.0000 - 0.4529i
   0.9656 + 0.0000i  -0.1311 - 0.2129i  -0.0698 + 0.0119i


bOut =

   1.9656 + 0.7935i
  -0.1132 + 0.3397i
   0.0915 + 0.6086i


q =

   0.5005 + 0.0000i   0.4529 + 0.0000i   0.7378 + 0.0000i
  -0.5152 + 0.1514i  -0.1132 + 0.5661i   0.4190 - 0.4502i
  -0.4041 + 0.5457i  -0.0000 - 0.6794i   0.2741 + 0.0468i


sv =

    3.9263
    2.0000
    0.7641


rwork =

         0
         0
    1.0000
    0.7071
         0
         0
         0
         0
         0
         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