NAG Library Routine Document

1Purpose

f02wdf returns the Householder $QU$ factorization of a real rectangular $m$ by $n$ $\left(m\ge n\right)$ matrix $A$. Further, on request or if $A$ is not of full rank, part or all of the singular value decomposition of $A$ is returned.

2Specification

Fortran Interface
 Subroutine f02wdf ( m, n, a, lda, b, tol, svd, z, sv, r, ldr, pt, ldpt, work,
 Integer, Intent (In) :: m, n, lda, ldr, ldpt, lwork Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: irank Real (Kind=nag_wp), Intent (In) :: tol Real (Kind=nag_wp), Intent (Inout) :: a(lda,n), b(m), r(ldr,n), pt(ldpt,n) Real (Kind=nag_wp), Intent (Out) :: z(n), sv(n), work(lwork) Logical, Intent (In) :: wantb, wantr, wantpt Logical, Intent (Inout) :: svd
#include nagmk26.h
 void f02wdf_ (const Integer *m, const Integer *n, double a[], const Integer *lda, const logical *wantb, double b[], const double *tol, logical *svd, Integer *irank, double z[], double sv[], const logical *wantr, double r[], const Integer *ldr, const logical *wantpt, double pt[], const Integer *ldpt, double work[], const Integer *lwork, Integer *ifail)

3Description

The real $m$ by $n$ $\left(m\ge n\right)$ matrix $A$ is first factorized as
 $A=Q U 0 ,$
where $Q$ is an $m$ by $m$ orthogonal matrix and $U$ is an $n$ by $n$ upper triangular matrix.
If either $U$ is singular or svd is supplied as .TRUE., then the singular value decomposition (SVD) of $U$ is obtained so that $U$ is factorized as
 $U=RDPT,$
where $R$ and $P$ are $n$ by $n$ orthogonal matrices and $D$ is the $n$ by $n$ diagonal matrix
 $D=diagsv1,sv2,…,svn,$
with $s{v}_{1}\ge s{v}_{2}\ge \cdots \ge s{v}_{n}\ge 0\text{.}$
Note that the SVD of $A$ is then given by
 $A=Q1 D 0 PT where Q1=Q R 0 0 I ,$
the diagonal elements of $D$ being the singular values of $A$.
The option to form a vector ${Q}^{\mathrm{T}}b$, or if appropriate ${Q}_{1}^{\mathrm{T}}b$, is also provided.
The rank of the matrix $A$, based upon a user-supplied argument tol, is also returned.
The $QU$ factorization of $A$ is obtained by Householder transformations. To obtain the SVD of $U$ the matrix is first reduced to bidiagonal form by means of plane rotations and then the $QR$ algorithm is used to obtain the SVD of the bidiagonal form.

4References

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

5Arguments

1:     $\mathbf{m}$ – IntegerInput
On entry: $m$, the number of rows of the matrix $A$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
2:     $\mathbf{n}$ – IntegerInput
On entry: $n$, the number of columns of the matrix $A$.
Constraint: $1\le {\mathbf{n}}\le {\mathbf{m}}$.
3:     $\mathbf{a}\left({\mathbf{lda}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: the leading $m$ by $n$ part of a must contain the matrix to be factorized.
On exit: the leading $m$ by $n$ part of a, together with the $n$-element vector z, contains details of the Householder $QU$ factorization.
Details of the storage of the $QU$ factorization are given in Section 9.4.
4:     $\mathbf{lda}$ – IntegerInput
On entry: the first dimension of the array a as declared in the (sub)program from which f02wdf is called.
Constraint: ${\mathbf{lda}}\ge {\mathbf{m}}$.
5:     $\mathbf{wantb}$ – LogicalInput
On entry: must be .TRUE. if ${Q}^{\mathrm{T}}b$ or ${Q}_{1}^{\mathrm{T}}b$ is required.
If on entry ${\mathbf{wantb}}=\mathrm{.FALSE.}$, b is not referenced.
6:     $\mathbf{b}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: if wantb is supplied as .TRUE., b must contain the $m$ element vector $b$. Otherwise, b is not referenced.
On exit: contains ${Q}_{1}^{\mathrm{T}}b$ if svd is returned as .TRUE. and ${Q}^{\mathrm{T}}b$ if svd is returned as .FALSE..
7:     $\mathbf{tol}$ – Real (Kind=nag_wp)Input
On entry: must specify a relative tolerance to be used to determine the rank of $A$. tol should be chosen as approximately the largest relative error in the elements of $A$. For example, if the elements of $A$ are correct to about $4$ significant figures, tol should be set to about $5×{10}^{-4}$. See Section 9.3 for a description of how tol is used to determine rank.
If tol is outside the range $\left(\epsilon ,1.0\right)$, where $\epsilon$ is the machine precision, the value $\epsilon$ is used in place of tol. For most problems this is unreasonably small.
8:     $\mathbf{svd}$ – LogicalInput/Output
On entry: must be .TRUE. if the singular values are to be found even if $A$ is of full rank.
If before entry, ${\mathbf{svd}}=\mathrm{.FALSE.}$ and $A$ is determined to be of full rank, only the $QU$ factorization of $A$ is computed.
On exit: is returned as .FALSE. if only the $QU$ factorization of $A$ has been obtained and is returned as .TRUE. if the singular values of $A$ have been obtained.
9:     $\mathbf{irank}$ – IntegerOutput
On exit: returns the rank of the matrix $A$. (It should be noted that it is possible for irank to be returned as $n$ and svd to be returned as .TRUE., even if svd was supplied as .FALSE.. This means that the matrix $U$ only just failed the test for nonsingularity.)
10:   $\mathbf{z}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: the $n$-element vector z contains some details of the Householder transformations. See Section 9.4 for further information.
11:   $\mathbf{sv}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: if svd is returned as .TRUE., sv contains the $n$ singular values of $A$ arranged in descending order.
12:   $\mathbf{wantr}$ – LogicalInput
On entry: must be .TRUE. if the orthogonal matrix $R$ is required when the singular values are computed.
If on entry ${\mathbf{wantr}}=\mathrm{.FALSE.}$, r is not referenced.
13:   $\mathbf{r}\left({\mathbf{ldr}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the second dimension of the array r must be at least ${\mathbf{n}}$ if ${\mathbf{wantr}}=\mathrm{.TRUE.}$, and at least $1$ otherwise.
On exit: if svd is returned as .TRUE. and wantr was supplied as .TRUE., the leading $n$ by $n$ part of r will contain the left-hand orthogonal matrix of the svd of $U$.
14:   $\mathbf{ldr}$ – IntegerInput
On entry: the first dimension of the array r as declared in the (sub)program from which f02wdf is called.
Constraints:
• if ${\mathbf{wantr}}=\mathrm{.TRUE.}$, ${\mathbf{ldr}}\ge {\mathbf{n}}$;
• otherwise ${\mathbf{ldr}}\ge 1$.
15:   $\mathbf{wantpt}$ – LogicalInput
On entry: must be .TRUE. if the orthogonal matrix ${P}^{\mathrm{T}}$ 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 argument pt.
16:   $\mathbf{pt}\left({\mathbf{ldpt}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: if svd is returned as .TRUE. and wantpt was supplied as .TRUE., the leading $n$ by $n$ part of pt contains the orthogonal matrix ${P}^{\mathrm{T}}$.
If svd is returned as .TRUE., but wantpt was supplied as .FALSE., the leading $n$ by $n$ part of pt is used for internal workspace.
17:   $\mathbf{ldpt}$ – IntegerInput
On entry: the first dimension of the array pt as declared in the (sub)program from which f02wdf is called.
Constraint: ${\mathbf{ldpt}}\ge {\mathbf{n}}$.
18:   $\mathbf{work}\left({\mathbf{lwork}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: if svd is returned as .FALSE., ${\mathbf{work}}\left(1\right)$ contains the condition number ${‖U‖}_{E}{‖{U}^{-1}‖}_{E}$ of the upper triangular matrix $U$.
If svd is returned as .TRUE., ${\mathbf{work}}\left(1\right)$ will contain the total number of iterations taken by the $QR$ algorithm.
The rest of the array is used as workspace and so contains no meaningful information.
19:   $\mathbf{lwork}$ – IntegerInput
On entry: the dimension of the array work as declared in the (sub)program from which f02wdf is called.
Constraint: ${\mathbf{lwork}}\ge 3×{\mathbf{n}}$.
20:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{n}}<1$, or ${\mathbf{m}}<{\mathbf{n}}$, or ${\mathbf{lda}}<{\mathbf{m}}$, or ${\mathbf{ldr}}<{\mathbf{n}}$ when ${\mathbf{wantr}}=\mathrm{.TRUE.}$, or ${\mathbf{ldpt}}<{\mathbf{n}}$ or ${\mathbf{lwork}}<3×{\mathbf{n}}$.
(The routine only checks ldr if wantr is supplied as .TRUE..)
${\mathbf{ifail}}>1$
The $QR$ algorithm has failed to converge to the singular values in $50×{\mathbf{n}}$ iterations. In this case ${\mathbf{sv}}\left(1\right),{\mathbf{sv}}\left(2\right),\dots ,{\mathbf{sv}}\left({\mathbf{ifail}}-1\right)$ may not have been correctly found and the remaining singular values may not be the smallest singular values. The matrix $A$ has nevertheless been factorized as $A={Q}_{1}C{P}^{\mathrm{T}}$, where $C$ is an upper bidiagonal matrix with ${\mathbf{sv}}\left(1\right),{\mathbf{sv}}\left(2\right),\dots ,{\mathbf{sv}}\left(n\right)$ as its diagonal elements and ${\mathbf{work}}\left(2\right),{\mathbf{work}}\left(3\right),\dots ,{\mathbf{work}}\left(n\right)$ as its superdiagonal elements.
This failure cannot occur if svd is returned as .FALSE. and in any case is extremely rare.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7Accuracy

The computed factors $Q$, $U$, $R$, $D$ and ${P}^{\mathrm{T}}$ satisfy the relations
 $Q U 0 =A+E,$
 $Q R 0 0 I D 0 PT=A+F$
where ${‖E‖}_{2}\le {c}_{1}\epsilon {‖A‖}_{2}$, ${‖F‖}_{2}\le {c}_{2}\epsilon {‖A‖}_{2}$,
$\epsilon$ being the machine precision and ${c}_{1}$ and ${c}_{2}$ are modest functions of $m$ and $n$. Note that ${‖A‖}_{2}=s{v}_{1}$.

8Parallelism and Performance

f02wdf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f02wdf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9.1Timing

The time taken by f02wdf to obtain the Householder $QU$ factorization is approximately proportional to ${n}^{2}\left(3m-n\right)$.
The additional time taken to obtain the singular value decomposition is approximately proportional to ${n}^{3}$, where the constant of proportionality depends upon whether or not the orthogonal matrices $R$ and ${P}^{\mathrm{T}}$ are required.

9.2General 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.
Unless otherwise stated in the Users' Note for your implementation, the routine may be called with the same array for arguments z and sv, in which case, if svd is returned as .TRUE., the singular values will overwrite the original contents of z; also, if ${\mathbf{wantpt}}=\mathrm{.FALSE.}$, it may be called with the same array for arguments r and pt. However this is not standard Fortran, and may not work on all systems.
This routine is called by the least squares routine f04jgf.

9.3Determining the Rank of $A$

Following the $QU$ factorization of $A$, if svd is supplied as .FALSE., then the condition number of $U$ given by
 $CU=UF U-1F$
is found, where ${‖.‖}_{F}$ denotes the Frobenius norm, and if $C\left(U\right)$ is such that
 $CU×tol>1.0$
then $U$ is regarded as singular and the singular values of $A$ are computed. If this test is not satisfied, then the rank of $A$ is set to $n$. Note that if svd is supplied as .TRUE. then this test is omitted.
When the singular values are computed, then the rank of $A$, $r$, is returned as the largest integer such that
 $svr>tol×sv1,$
unless $s{v}_{1}=0$ in which case $r$ is returned as zero. That is, singular values which satisfy $s{v}_{i}\le {\mathbf{tol}}×s{v}_{1}$ are regarded as negligible because relative perturbations of order tol can make such singular values zero.

9.4Storage Details of the $QU$ Factorization

The $k$th Householder transformation matrix, ${T}_{k}$, used in the $QU$ factorization is chosen to introduce the zeros into the $k$th column and has the form
 $Tk=I-2 0 u 0 uT , uTu=1,$
where $u$ is an $\left(m-k+1\right)$ element vector.
In place of $u$ the routine actually computes the vector $z$ given by
 $z=2u1u.$
The first element of $z$ is stored in ${\mathbf{z}}\left(k\right)$ and the remaining elements of $z$ are overwritten on the subdiagonal elements of the $k$th column of a. The upper triangular matrix $U$ is overwritten on the $n$ by $n$ upper triangular part of a.

10Example

This example obtains the rank and the singular value decomposition of the $6$ by $4$ matrix $A$ given by
 $A= 22.25 31.75 -38.25 65.50 20.00 26.75 28.50 -26.50 -15.25 24.25 27.75 18.50 27.25 10.00 3.00 2.00 -17.25 -30.75 11.25 7.50 17.25 30.75 -11.25 -7.50$
the value tol to be taken as $5×{10}^{-4}$.

10.1Program Text

Program Text (f02wdfe.f90)

10.2Program Data

Program Data (f02wdfe.d)

10.3Program Results

Program Results (f02wdfe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017