NAG Library Routine Document
f04jgf (real_gen_solve)
1
Purpose
f04jgf finds the solution of a linear least squares problem, $Ax=b$, where $A$ is a real $m$ by $n\left(m\ge n\right)$ matrix and $b$ is an $m$ element vector. If the matrix of observations is not of full rank, then the minimal least squares solution is returned.
2
Specification
Fortran Interface
Subroutine f04jgf ( 
m, n, a, lda, b, tol, svd, sigma, irank, work, lwork, ifail) 
Integer, Intent (In)  ::  m, n, lda, 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)  Real (Kind=nag_wp), Intent (Out)  ::  sigma, work(lwork)  Logical, Intent (Out)  ::  svd 

C Header Interface
#include <nagmk26.h>
void 
f04jgf_ (const Integer *m, const Integer *n, double a[], const Integer *lda, double b[], const double *tol, logical *svd, double *sigma, Integer *irank, double work[], const Integer *lwork, Integer *ifail) 

3
Description
The minimal least squares solution of the problem $Ax=b$ is the vector $x$ of minimum (Euclidean) length which minimizes the length of the residual vector $r=bAx$.
The real
$m$ by
$n\left(m\ge n\right)$ matrix
$A$ is factorized as
where
$Q$ is an
$m$ by
$m$ orthogonal matrix and
$R$ is an
$n$ by
$n$ upper triangular matrix. If
$R$ is of full rank, then the least squares solution is given by
If
$R$ is not of full rank, then the singular value decomposition of
$R$ is obtained so that
$R$ is factorized as
where
$U$ and
$V$ are
$n$ by
$n$ orthogonal matrices and
$D$ is the
$n$ by
$n$ diagonal matrix
with
${\sigma}_{1}\ge {\sigma}_{2}\ge \dots \ge {\sigma}_{n}\ge 0$, these being the singular values of
$A$. If the singular values
${\sigma}_{k+1},\dots ,{\sigma}_{n}$ are negligible, but
${\sigma}_{k}$ is not negligible, relative to the data errors in
$A$, then the rank of
$A$ is taken to be
$k$ and the minimal least squares solution is given by
where
$S=\mathrm{diag}\left({\sigma}_{1},{\sigma}_{2},\dots ,{\sigma}_{k}\right)$.
The routine also returns the value of the standard error
${r}^{\mathrm{T}}r$ being the residual sum of squares and
$k$ the rank of
$A$.
4
References
Lawson C L and Hanson R J (1974) Solving Least Squares Problems Prentice–Hall
5
Arguments
 1: $\mathbf{m}$ – IntegerInput

On entry:
$m$, the number of rows of
a.
Constraint:
${\mathbf{m}}\ge {\mathbf{n}}$.
 2: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of columns of ${\mathbf{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 $m$ by $n$ matrix $A$.
On exit: if
svd is returned as .FALSE.,
a is overwritten by details of the
$QR$ factorization of
$A$.
If
svd is returned as .TRUE., the first
$n$ rows of
a are overwritten by the righthand singular vectors, stored by rows; and the remaining rows of the array are used as workspace.
 4: $\mathbf{lda}$ – IntegerInput

On entry: the first dimension of the array
a as declared in the (sub)program from which
f04jgf is called.
Constraint:
${\mathbf{lda}}\ge {\mathbf{m}}$.
 5: $\mathbf{b}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the righthand side vector $b$.
On exit: the first
$n$ elements of
b contain the minimal least squares solution vector
$x$. The remaining
$mn$ elements are used for workspace.
 6: $\mathbf{tol}$ – Real (Kind=nag_wp)Input

On entry: 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 then
tol should be set to about
$5\times {10}^{4}$. See
Section 9 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.
 7: $\mathbf{svd}$ – LogicalOutput

On exit: is returned as .FALSE. if the least squares solution has been obtained from the
$QR$ factorization of
$A$. In this case
$A$ is of full rank.
svd is returned as .TRUE. if the least squares solution has been obtained from the singular value decomposition of
$A$.
 8: $\mathbf{sigma}$ – Real (Kind=nag_wp)Output

On exit: the standard error, i.e., the value $\sqrt{{r}^{\mathrm{T}}r/\left(mk\right)}$ when $m>k$, and the value zero when $m=k$. Here $r$ is the residual vector $bAx$ and $k$ is the rank of $A$.
 9: $\mathbf{irank}$ – IntegerOutput

On exit:
$k$, 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.. This means that the matrix
$R$ only just failed the test for nonsingularity.
 10: $\mathbf{work}\left({\mathbf{lwork}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: if
svd is returned as .FALSE., then the first
$n$ elements of
work contain information on the
$QR$ factorization of
$A$ (see argument
a above), and
${\mathbf{work}}\left(n+1\right)$ contains the condition number
${\Vert R\Vert}_{E}{\Vert {R}^{1}\Vert}_{E}$ of the upper triangular matrix
$R$.
If
svd is returned as .TRUE., then the first
$n$ elements of
work contain the singular values of
$A$ arranged in descending order and
${\mathbf{work}}\left(n+1\right)$ contains the total number of iterations taken by the
$QR$ algorithm. The rest of
work is used as workspace.
 11: $\mathbf{lwork}$ – IntegerInput

On entry: the dimension of the array
work as declared in the (sub)program from which
f04jgf is called.
Constraint:
${\mathbf{lwork}}\ge 4\times {\mathbf{n}}$.
 12: $\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).
6
Error 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{lda}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lda}}\ge {\mathbf{m}}$.
On entry,
lwork is too small. Minimum size required:
$\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{n}}\ge 1$.
 ${\mathbf{ifail}}=2$

The $QR$ algorithm has failed to converge to the singular values in $50\times {\mathbf{n}}$ iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
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.
7
Accuracy
The computed factors
$Q$,
$R$,
$U$,
$D$ and
${V}^{\mathrm{T}}$ satisfy the relations
where
$\epsilon $ being the
machine precision, and
${c}_{1}$ and
${c}_{2}$ being modest functions of
$m$ and
$n$. Note that
${\Vert A\Vert}_{2}={\sigma}_{1}$.
For a fuller discussion, covering the accuracy of the solution
$x$ see
Lawson and Hanson (1974), especially pages 50 and 95.
8
Parallelism and Performance
f04jgf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f04jgf 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 implementationspecific information.
If the least squares solution is obtained from the $QR$ factorization then the time taken by the routine is approximately proportional to ${n}^{2}\left(3mn\right)$. If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to ${n}^{2}\left(3m+19n\right)$. The approximate proportionality factor is the same in each case.
This routine is column biased and so is suitable for use in paged environments.
Following the
$QR$ factorization of
$A$ the condition number
is determined and if
$c\left(R\right)$ is such that
then
$R$ is regarded as singular and the singular values of
$A$ are computed. If this test is not satisfied,
$R$ is regarded as nonsingular and the rank of
$A$ is set to
$n$. When the singular values are computed the rank of
$A$, say
$k$, is returned as the largest integer such that
unless
${\sigma}_{1}=0$ in which case
$k$ is returned as zero. That is, singular values which satisfy
${\sigma}_{i}\le {\mathbf{tol}}\times {\sigma}_{1}$ are regarded as negligible because relative perturbations of order
tol can make such singular values zero.
10
Example
This example obtains a least squares solution for
$r=bAx$, where
and the value
tol is to be taken as
$5\times {10}^{4}$.
10.1
Program Text
Program Text (f04jgfe.f90)
10.2
Program Data
Program Data (f04jgfe.d)
10.3
Program Results
Program Results (f04jgfe.r)