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_qu_svd (f02wd)

Purpose

nag_eigen_real_gen_qu_svd (f02wd) returns the Householder QUQU factorization of a real rectangular mm by nn (mn)(mn) matrix AA. Further, on request or if AA is not of full rank, part or all of the singular value decomposition of AA is returned.

Syntax

[a, b, svd, irank, z, sv, r, pt, work, ifail] = f02wd(a, wantb, b, tol, svd, wantr, wantpt, lwork, 'm', m, 'n', n)
[a, b, svd, irank, z, sv, r, pt, work, ifail] = nag_eigen_real_gen_qu_svd(a, wantb, b, tol, svd, wantr, wantpt, lwork, 'm', m, 'n', n)

Description

The real mm by nn (mn)(mn) matrix AA is first factorized as
A = Q
(U)
0
,
A=Q U 0 ,
where QQ is an mm by mm orthogonal matrix and UU is an nn by nn upper triangular matrix.
If either UU is singular or svd is supplied as true, then the singular value decomposition (SVD) of UU is obtained so that UU is factorized as
U = RDPT,
U=RDPT,
where RR and PP are nn by nn orthogonal matrices and DD is the nn by nn diagonal matrix
D = diag(sv1,sv2,,svn),
D=diag(sv1,sv2,,svn),
with sv1sv2svn0.sv1sv2svn0.
Note that the SVD of AA is then given by
A = Q1
(D)
0
PT  where  Q1 = Q
(R0)
0 I
,
A=Q1 D 0 PT  where  Q1=Q R 0 0 I ,
the diagonal elements of DD being the singular values of AA.
The option to form a vector QTbQTb, or if appropriate Q1T b Q1T b, is also provided.
The rank of the matrix AA, based upon a user-supplied parameter tol, is also returned.
The QUQU factorization of AA is obtained by Householder transformations. To obtain the SVD of UU the matrix is first reduced to bidiagonal form by means of plane rotations and then the QRQR algorithm is used to obtain the SVD of the bidiagonal form.

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:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldamldam.
The leading mm by nn part of a must contain the matrix to be factorized.
2:     wantb – logical scalar
Must be true if QTbQTb or Q1Tb Q1Tb is required.
If on entry wantb = falsewantb=false, b is not referenced.
3:     b(m) – double array
m, the dimension of the array, must satisfy the constraint mnmn.
If wantb is supplied as true, b must contain the mm element vector bb. Otherwise, b is not referenced.
4:     tol – double scalar
Must specify a relative tolerance to be used to determine the rank of AA. tol should be chosen as approximately the largest relative error in the elements of AA. For example, if the elements of AA are correct to about 44 significant figures, tol should be set to about 5 × 1045×10-4. See Section [Determining the Rank of ] for a description of how tol is used to determine rank.
If tol is outside the range (ε,1.0)(ε,1.0), where εε is the machine precision, the value εε is used in place of tol. For most problems this is unreasonably small.
5:     svd – logical scalar
Must be true if the singular values are to be found even if AA is of full rank.
If before entry, svd = falsesvd=false and AA is determined to be of full rank, only the QUQU factorization of AA is computed.
6:     wantr – logical scalar
Must be true if the orthogonal matrix RR is required when the singular values are computed.
If on entry wantr = falsewantr=false, r is not referenced.
7:     wantpt – logical scalar
Must be true if the orthogonal matrix PTPT is required when the singular values are computed.
Note that if svd is returned as true, pt is referenced even if wantpt is supplied as false, but see parameter pt.
8:     lwork – int64int32nag_int scalar
The dimension of the array work as declared in the (sub)program from which nag_eigen_real_gen_qu_svd (f02wd) is called.
Constraint: lwork3 × nlwork3×n.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The dimension of the array b and the first dimension of the array a. (An error is raised if these dimensions are not equal.)
mm, the number of rows of the matrix AA.
Constraint: mnmn.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
nn, the number of columns of the matrix AA.
Constraint: 1nm1nm.

Input Parameters Omitted from the MATLAB Interface

lda ldr ldpt

Output Parameters

1:     a(lda,n) – double array
ldamldam.
The leading mm by nn part of a, together with the nn-element vector z, contains details of the Householder QUQU factorization.
Details of the storage of the QUQU factorization are given in Section [Storage Details of the Factorization].
2:     b(m) – double array
Contains Q1TbQ1Tb if svd is returned as true and QTbQTb if svd is returned as false.
3:     svd – logical scalar
Is returned as false if only the QUQU factorization of AA has been obtained and is returned as true if the singular values of AA have been obtained.
4:     irank – int64int32nag_int scalar
Returns the rank of the matrix AA. (It should be noted that it is possible for irank to be returned as nn and svd to be returned as true, even if svd was supplied as false. This means that the matrix UU only just failed the test for nonsingularity.)
5:     z(n) – double array
The nn-element vector z contains some details of the Householder transformations. See Section [Storage Details of the Factorization] for further information.
6:     sv(n) – double array
If svd is returned as true, sv contains the nn singular values of AA arranged in descending order.
7:     r(ldr,n) – double array
The first dimension, ldr, of the array r will be
  • if wantr = truewantr=true, ldrnldrn;
  • otherwise 11.
The second dimension of the array will be nn if wantr = truewantr=true, and at least 11 otherwise
If svd is returned as true and wantr was supplied as true, the leading nn by nn part of r will contain the left-hand orthogonal matrix of the svd of UU.
8:     pt(ldpt,n) – double array
ldptnldptn.
If svd is returned as true and wantpt was supplied as true, the leading nn by nn part of pt contains the orthogonal matrix PTPT.
If svd is returned as true, but wantpt was supplied as false, the leading nn by nn part of pt is used for internal workspace.
9:     work(lwork) – double array
If svd is returned as false, work(1)work1 contains the condition number UEU1EUEU-1E of the upper triangular matrix UU.
If svd is returned as true, work(1)work1 will contain the total number of iterations taken by the QRQR algorithm.
The rest of the array is used as workspace and so contains no meaningful information.
10:   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:
  ifail = 1ifail=1
On entry,n < 1n<1,
orm < nm<n,
orlda < mlda<m,
orldr < nldr<n when wantr = truewantr=true,
orldpt < nldpt<n 
orlwork < 3 × nlwork<3×n.
(The function only checks ldr if wantr is supplied as true.)
  ifail > 1ifail>1
The QRQR algorithm has failed to converge to the singular values in 50 × n50×n iterations. In this case sv(1),sv(2),,sv(ifail1)sv1,sv2,,svifail-1 may not have been correctly found and the remaining singular values may not be the smallest singular values. The matrix AA has nevertheless been factorized as A = Q1CPTA=Q1CPT, where CC is an upper bidiagonal matrix with sv(1),sv(2),,sv(n)sv1,sv2,,svn as its diagonal elements and work(2),work(3),,work(n)work2,work3,,workn as its superdiagonal elements.
This failure cannot occur if svd is returned as false and in any case is extremely rare.

Accuracy

The computed factors QQ, UU, RR, DD and PTPT satisfy the relations
Q
(U)
0
= A + E,
Q U 0 =A+E,
Q
(R0)
0 I
(D)
0
PT = A + F
Q R 0 0 I D 0 PT=A+F
where E2c1ε A2E2c1ε A2, F2c2ε A2F2c2ε A2,
εε being the machine precision and c1c1 and c2c2 are modest functions of mm and nn. Note that A2 = sv1A2=sv1.

Further Comments

Timing

The time taken by nag_eigen_real_gen_qu_svd (f02wd) to obtain the Householder QUQU factorization is approximately proportional to n2(3mn)n2(3m-n).
The additional time taken to obtain the singular value decomposition is approximately proportional to n3n3, where the constant of proportionality depends upon whether or not the orthogonal matrices RR and PTPT are required.

General Remarks

Singular vectors associated with a zero or multiple singular value, are not uniquely determined, even in exact arithmetic, and very different results may be obtained if they are computed on different machines.
This function is called by the least squares function nag_linsys_real_gen_solve (f04jg).

Determining the Rank of A

Following the QUQU factorization of AA, if svd is supplied as false, then the condition number of UU given by
C(U) = UF U1F
C(U)=UF U-1F
is found, where .F.F denotes the Frobenius norm, and if C(U)C(U) is such that
C(U) × tol > 1.0
C(U)×tol>1.0
then UU is regarded as singular and the singular values of AA are computed. If this test is not satisfied, then the rank of AA is set to nn. Note that if svd is supplied as true then this test is omitted.
When the singular values are computed, then the rank of AA, rr, is returned as the largest integer such that
svr > tol × sv1,
svr>tol×sv1,
unless sv1 = 0sv1=0 in which case rr is returned as zero. That is, singular values which satisfy svitol × sv1svitol×sv1 are regarded as negligible because relative perturbations of order tol can make such singular values zero.

Storage Details of the QU Factorization

The kkth Householder transformation matrix, TkTk, used in the QUQU factorization is chosen to introduce the zeros into the kkth column and has the form
Tk = I2
(0)
u
(0uT)
,  uTu = 1,
Tk=I-2 0 u 0 uT ,  uTu=1,
where uu is an (mk + 1)(m-k+1) element vector.
In place of uu the function actually computes the vector zz given by
z = 2u1u.
z=2u1u.
The first element of zz is stored in z(k)zk and the remaining elements of zz are overwritten on the subdiagonal elements of the kkth column of a. The upper triangular matrix UU is overwritten on the nn by nn upper triangular part of a.

Example

function nag_eigen_real_gen_qu_svd_example
a = [22.25, 31.75, -38.25, 65.5;
     20, 26.75, 28.5, -26.5;
     -15.25, 24.25, 27.75, 18.5;
     27.25, 10, 3, 2;
     -17.25, -30.75, 11.25, 7.5;
     17.25, 30.75, -11.25, -7.5];
wantb = false;
b = [0;
     0;
     0;
     0;
     0;
     0];
tol = 0.0005;
svd = true;
wantr = true;
wantpt = true;
lwork = int64(24);
[aOut, bOut, svdOut, irank, z, sv, r, pt, work, ifail] = ...
    nag_eigen_real_gen_qu_svd(a, wantb, b, tol, svd, wantr, wantpt, lwork)
 

aOut =

  -49.6519  -44.4092   20.3542   -8.8818
    0.4028  -48.2767   -9.5887  -20.3761
   -0.3071    0.8369   52.9270  -48.8806
    0.5488   -0.3907   -0.8364  -50.6742
   -0.3474   -0.2585   -0.1851    0.6321
    0.3474    0.2585    0.1851   -0.6321


bOut =

     0
     0
     0
     0
     0
     0


svdOut =

     1


irank =

                    4


z =

    1.4481
    1.1153
    1.4817
    1.4482


sv =

   91.0000
   68.2500
   45.5000
   22.7500


r =

    0.5639   -0.6344    0.4229    0.3172
    0.3512   -0.3951   -0.6791   -0.5093
    0.6403    0.5692    0.3095   -0.4126
    0.3855    0.3427   -0.5140    0.6854


pt =

   -0.3077   -0.4615    0.4615   -0.6923
    0.4615    0.6923    0.3077   -0.4615
   -0.4615    0.3077    0.6923    0.4615
   -0.6923    0.4615   -0.4615   -0.3077


work =

    1.0000
  -49.6523
   44.2299
   28.9880
         0
         0
    0.8037
    0.9584
         0
         0
   -0.5951
   -0.2854
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0


function f02wd_example
a = [22.25, 31.75, -38.25, 65.5;
     20, 26.75, 28.5, -26.5;
     -15.25, 24.25, 27.75, 18.5;
     27.25, 10, 3, 2;
     -17.25, -30.75, 11.25, 7.5;
     17.25, 30.75, -11.25, -7.5];
wantb = false;
b = [0;
     0;
     0;
     0;
     0;
     0];
tol = 0.0005;
svd = true;
wantr = true;
wantpt = true;
lwork = int64(24);
[aOut, bOut, svdOut, irank, z, sv, r, pt, work, ifail] = ...
    f02wd(a, wantb, b, tol, svd, wantr, wantpt, lwork)
 

aOut =

  -49.6519  -44.4092   20.3542   -8.8818
    0.4028  -48.2767   -9.5887  -20.3761
   -0.3071    0.8369   52.9270  -48.8806
    0.5488   -0.3907   -0.8364  -50.6742
   -0.3474   -0.2585   -0.1851    0.6321
    0.3474    0.2585    0.1851   -0.6321


bOut =

     0
     0
     0
     0
     0
     0


svdOut =

     1


irank =

                    4


z =

    1.4481
    1.1153
    1.4817
    1.4482


sv =

   91.0000
   68.2500
   45.5000
   22.7500


r =

    0.5639   -0.6344    0.4229    0.3172
    0.3512   -0.3951   -0.6791   -0.5093
    0.6403    0.5692    0.3095   -0.4126
    0.3855    0.3427   -0.5140    0.6854


pt =

   -0.3077   -0.4615    0.4615   -0.6923
    0.4615    0.6923    0.3077   -0.4615
   -0.4615    0.3077    0.6923    0.4615
   -0.6923    0.4615   -0.4615   -0.3077


work =

    1.0000
  -49.6523
   44.2299
   28.9880
         0
         0
    0.8037
    0.9584
         0
         0
   -0.5951
   -0.2854
         0
         0
         0
         0
         0
         0
         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