Integer type:** int32**** int64**** nag_int** show int32 show int32 show int64 show int64 show nag_int show nag_int

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

The real m$m$ by n$n$ (m ≥ n)$(m\ge n)$ matrix A$A$ is first factorized as

where Q$Q$ is an m$m$ by m$m$ orthogonal matrix and U$U$ is an n$n$ by n$n$ upper triangular matrix.

A = Q
$$A=Q\left(\begin{array}{l}U\\ 0\end{array}\right)\text{,}$$ |

If either U$U$ is singular or svd is supplied as **true**, then the singular value decomposition (SVD) of U$U$ is obtained so that U$U$ is factorized as

where R$R$ and P$P$ are n$n$ by n$n$ orthogonal matrices and D$D$ is the n$n$ by n$n$ diagonal matrix

with sv_{1} ≥ sv_{2} ≥ ⋯ ≥ sv_{n} ≥ 0.$s{v}_{1}\ge s{v}_{2}\ge \cdots \ge s{v}_{n}\ge 0\text{.}$

U = RDP ^{T},
$$U=RD{P}^{\mathrm{T}}\text{,}$$ |

D = diag(sv _{1},sv_{2}, … ,sv_{n}),
$$D=\mathrm{diag}(s{v}_{1},s{v}_{2},\dots ,s{v}_{n})\text{,}$$ |

Note that the SVD of A$A$ is then given by

the diagonal elements of D$D$ being the singular values of A$A$.

A = Q _{1}
^{T} where Q_{1} = Q
$$A={Q}_{1}\left(\begin{array}{l}D\\ 0\end{array}\right){P}^{\mathrm{T}}\text{\hspace{1em} where \hspace{1em}}{Q}_{1}=Q\left(\begin{array}{ll}R& 0\\ 0& I\end{array}\right)\text{,}$$ |

The option to form a vector Q^{T}b${Q}^{\mathrm{T}}b$, or if appropriate
Q_{1}^{T}
b${Q}_{1}^{\mathrm{T}}b$, is also provided.

The QU$QU$ factorization of A$A$ is obtained by Householder transformations. To obtain the SVD of U$U$ the matrix is first reduced to bidiagonal form by means of plane rotations and then the QR$QR$ algorithm is used to obtain the SVD of the bidiagonal form.

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

- 1: a(lda,n) – double array
- lda, the first dimension of the array, must satisfy the constraint lda ≥ m$\mathit{lda}\ge {\mathbf{m}}$.
- 2: wantb – logical scalar
- 3: b(m) – double array
- 4: tol – double scalar
- Must specify a relative tolerance to be used to determine the rank of A$A$. tol should be chosen as approximately the largest relative error in the elements of A$A$. For example, if the elements of A$A$ are correct to about 4$4$ significant figures, tol should be set to about 5 × 10
^{ − 4}$5\times {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)$(\epsilon ,1.0)$, where ε$\epsilon $ is the machine precision, the value ε$\epsilon $ 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 A$A$ is of full rank.If before entry, svd = false${\mathbf{svd}}=\mathbf{false}$**and**A$A$ is determined to be of full rank, only the QU$QU$ factorization of A$A$ is computed. - 6: wantr – logical scalar
- 7: wantpt – logical scalar
- 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.

- 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.)m$m$, the number of rows of the matrix A$A$.- 2: n – int64int32nag_int scalar
*Default*: The second dimension of the array a.n$n$, the number of columns of the matrix A$A$.

- lda ldr ldpt

- 1: a(lda,n) – double array
- lda ≥ m$\mathit{lda}\ge {\mathbf{m}}$.The leading m$m$ by n$n$ part of a, together with the n$n$-element vector z, contains details of the Householder QU$QU$ factorization.Details of the storage of the QU$QU$ factorization are given in Section [Storage Details of the Factorization].
- 2: b(m) – double array
- 3: svd – logical scalar
- Is returned as
**false**if only the QU$QU$ factorization of A$A$ has been obtained and is returned as**true**if the singular values of A$A$ have been obtained. - 4: irank – int64int32nag_int scalar
- 5: z(n) – double array
- The n$n$-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
- 7: r(ldr,n) – double array
- 8: pt(ldpt,n) – double array
- ldpt ≥ n$\mathit{ldpt}\ge {\mathbf{n}}$.
- 9: work(lwork) – double array
- If svd is returned as
**false**, work(1)${\mathbf{work}}\left(1\right)$ contains the condition number ‖U‖_{E}‖U^{ − 1}‖_{E}${\Vert U\Vert}_{E}{\Vert {U}^{-1}\Vert}_{E}$ of the upper triangular matrix U$U$.If svd is returned as**true**, work(1)${\mathbf{work}}\left(1\right)$ will contain the total number of iterations taken by the QR$QR$ algorithm.The rest of the array is used as workspace and so contains no meaningful information. - 10: ifail – int64int32nag_int scalar
- ifail = 0${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Errors or warnings detected by the function:

On entry, n < 1${\mathbf{n}}<1$, or m < n${\mathbf{m}}<{\mathbf{n}}$, or lda < m$\mathit{lda}<{\mathbf{m}}$, or ldr < n$\mathit{ldr}<{\mathbf{n}}$ when wantr = true${\mathbf{wantr}}=\mathbf{true}$, or ldpt < n$\mathit{ldpt}<{\mathbf{n}}$ or lwork < 3 × n${\mathbf{lwork}}<3\times {\mathbf{n}}$.

- The QR$QR$ algorithm has failed to converge to the singular values in 50 × n$50\times {\mathbf{n}}$ iterations. In this case sv(1),sv(2), … ,sv(ifail − 1)${\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$A$ has nevertheless been factorized as A = Q
_{1}CP^{T}$A={Q}_{1}C{P}^{\mathrm{T}}$, where C$C$ is an upper bidiagonal matrix with sv(1),sv(2), … ,sv(n)${\mathbf{sv}}\left(1\right),{\mathbf{sv}}\left(2\right),\dots ,{\mathbf{sv}}\left(n\right)$ as its diagonal elements and work(2),work(3), … ,work(n)${\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.

The computed factors Q$Q$, U$U$, R$R$, D$D$ and P^{T}${P}^{\mathrm{T}}$ satisfy the relations

where ‖E‖_{2} ≤ c_{1}ε ‖A‖_{2}${\Vert E\Vert}_{2}\le {c}_{1}\epsilon {\Vert A\Vert}_{2}$, ‖F‖_{2} ≤ c_{2}ε ‖A‖_{2}${\Vert F\Vert}_{2}\le {c}_{2}\epsilon {\Vert A\Vert}_{2}$,

Q
$$Q\left(\begin{array}{l}U\\ 0\end{array}\right)=A+E\text{,}$$ |

Q
^{T} = A + F
$$Q\left(\begin{array}{ll}R& 0\\ 0& I\end{array}\right)\left(\begin{array}{l}D\\ 0\end{array}\right){P}^{\mathrm{T}}=A+F$$ |

ε$\epsilon $ being the machine precision and c_{1}${c}_{1}$ and c_{2}${c}_{2}$ are modest functions of m$m$ and n$n$. Note that ‖A‖_{2} = sv_{1}${\Vert A\Vert}_{2}=s{v}_{1}$.

The time taken by nag_eigen_real_gen_qu_svd (f02wd) to obtain the Householder QU$QU$ factorization is approximately proportional to n^{2}(3m − n)${n}^{2}(3m-n)$.

The **additional** time taken to obtain the singular value decomposition is approximately proportional to n^{3}${n}^{3}$, where the constant of proportionality depends upon whether or not the orthogonal matrices R$R$ and P^{T}${P}^{\mathrm{T}}$ are required.

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).

Following the QU$QU$ factorization of A$A$, if svd is supplied as
**false**, then the condition number of U$U$ given by

is found, where ‖.‖_{F}${\Vert .\Vert}_{F}$ denotes the Frobenius norm, and if C(U)$C\left(U\right)$ is such that

then U$U$ is regarded as singular and the singular values of A$A$ are computed. If this test is not satisfied, then the rank of A$A$ is set to n$n$. Note that if svd is supplied as **true** then this test is omitted.

C(U) = ‖U‖ _{F} ‖U^{ − 1}‖_{F}
$$C\left(U\right)={\Vert U\Vert}_{F}{\Vert {U}^{-1}\Vert}_{F}$$ |

$$C\left(U\right)\times {\mathbf{tol}}>1.0$$ |

When the singular values are computed, then the rank of A$A$, r$r$, is returned as the largest integer such that

unless sv_{1} = 0$s{v}_{1}=0$ in which case r$r$ is returned as zero. That is, singular values which satisfy sv_{i} ≤ tol × sv_{1}$s{v}_{i}\le {\mathbf{tol}}\times s{v}_{1}$ are regarded as negligible because relative perturbations of order tol can make such singular values zero.

$$s{v}_{r}>{\mathbf{tol}}\times s{v}_{1}\text{,}$$ |

The k$k$th Householder transformation matrix, T_{k}${T}_{k}$, used in the QU$QU$ factorization is chosen to introduce the zeros into the k$k$th column and has the form

where u$u$ is an (m − k + 1)$(m-k+1)$ element vector.

T _{k} = I − 2
^{T}u = 1,
$${T}_{k}=I-2\left(\begin{array}{c}0\\ u\end{array}\right)\left(\begin{array}{cc}0& {u}^{\mathrm{T}}\end{array}\right)\text{, \hspace{1em}}{u}^{\mathrm{T}}u=1\text{,}$$ |

In place of u$u$ the function actually computes the vector z$z$ given by

The first element of z$z$ is stored in z(k)${\mathbf{z}}\left(k\right)$ and the remaining elements of z$z$ are overwritten on the subdiagonal elements of the k$k$th column of a. The upper triangular matrix U$U$ is overwritten on the n$n$ by n$n$ upper triangular part of a.

z = 2u _{1}u.
$$z=2{u}_{1}u\text{.}$$ |

Open in the MATLAB editor: nag_eigen_real_gen_qu_svd_example

function nag_eigen_real_gen_qu_svd_examplea = [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

Open in the MATLAB editor: f02wd_example

function f02wd_examplea = [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

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