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

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_linsys_real_gen_solve (f04jg)

## Purpose

nag_linsys_real_gen_solve (f04jg) finds the solution of a linear least squares problem, Ax = b$Ax=b$, where A$A$ is a real m$m$ by n(mn)$n\left(m\ge n\right)$ matrix and b$b$ is an m$m$ element vector. If the matrix of observations is not of full rank, then the minimal least squares solution is returned.

## Syntax

[a, b, svd, sigma, irank, work, ifail] = f04jg(a, b, tol, lwork, 'm', m, 'n', n)
[a, b, svd, sigma, irank, work, ifail] = nag_linsys_real_gen_solve(a, b, tol, lwork, 'm', m, 'n', n)

## Description

The minimal least squares solution of the problem Ax = b$Ax=b$ is the vector x$x$ of minimum (Euclidean) length which minimizes the length of the residual vector r = bAx$r=b-Ax$.
The real m$m$ by n(mn)$n\left(m\ge n\right)$ matrix A$A$ is factorized as
A = Q
 ( U ) 0
$A=Q U 0$
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. If U$U$ is of full rank, then the least squares solution is given by
 x = (U − 10)QTb. $x=(U-10)QTb.$
If U$U$ is not of full rank, then the singular value decomposition of U$U$ is obtained so that U$U$ is factorized as
 U = RDPT, $U=RDPT,$
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
 D = diag(σ1,σ2, … ,σn), $D=diag(σ1,σ2,…,σn),$
with σ1σ2σn0${\sigma }_{1}\ge {\sigma }_{2}\ge \dots \ge {\sigma }_{n}\ge 0$, these being the singular values of A$A$. If the singular values σk + 1,,σn${\sigma }_{k+1},\dots ,{\sigma }_{n}$ are negligible, but σk${\sigma }_{k}$ is not negligible, relative to the data errors in A$A$, then the rank of A$A$ is taken to be k$k$ and the minimal least squares solution is given by
x = P
 ( S − 1 0 ) 0 0
 ( RT 0 ) 0 I
QTb,
$x=P S-1 0 0 0 RT 0 0 I QTb,$
where S = diag(σ1,σ2,,σk)$S=\mathrm{diag}\left({\sigma }_{1},{\sigma }_{2},\dots ,{\sigma }_{k}\right)$.
This function obtains the factorizations by a call to nag_eigen_real_gen_qu_svd (f02wd).
The function also returns the value of the standard error
 σ = sqrt((rTr)/(m − k)), if ​m > k, = 0, if ​m = k,
$σ = rTr m-k , if ​m>k, = 0, if ​m=k,$
rTr${r}^{\mathrm{T}}r$ being the residual sum of squares and k$k$ the rank of A$A$.

## References

Lawson C L and Hanson R J (1974) Solving Least Squares Problems Prentice–Hall

## Parameters

### Compulsory Input Parameters

1:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldam$\mathit{lda}\ge {\mathbf{m}}$.
The m$m$ by n$n$ matrix A$A$.
2:     b(m) – double array
m, the dimension of the array, must satisfy the constraint mn${\mathbf{m}}\ge {\mathbf{n}}$.
The right-hand side vector b$b$.
3:     tol – double scalar
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 then tol should be set to about 5 × 104$5×{10}^{-4}$. See Section [Further Comments] for a description of how tol is used to determine rank. If tol is outside the range (ε,1.0)$\left(\epsilon ,1.0\right)$, where ε$\epsilon$ is the machine precision, then the value ε$\epsilon$ is used in place of tol. For most problems this is unreasonably small.
4:     lwork – int64int32nag_int scalar
The dimension of the array work as declared in the (sub)program from which nag_linsys_real_gen_solve (f04jg) is called.
Constraint: lwork4 × n${\mathbf{lwork}}\ge 4×{\mathbf{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.)
m$m$, the number of rows of a.
Constraint: mn${\mathbf{m}}\ge {\mathbf{n}}$.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
n$n$, the number of columns of a${\mathbf{a}}$.
Constraint: 1nm$1\le {\mathbf{n}}\le {\mathbf{m}}$.

lda

### Output Parameters

1:     a(lda,n) – double array
ldam$\mathit{lda}\ge {\mathbf{m}}$.
If svd is returned as false, a stores details of the QU$QU$ factorization of A$A$ (see nag_eigen_real_gen_qu_svd (f02wd) for further details).
If svd is returned as true, the first n$n$ rows of a store the right-hand singular vectors, stored by rows; and the remaining rows of the array are used as workspace.
2:     b(m) – double array
The first n$n$ elements of b contain the minimal least squares solution vector x$x$. The remaining mn$m-n$ elements are used for workspace.
3:     svd – logical scalar
Is returned as false if the least squares solution has been obtained from the QU$QU$ factorization of A$A$. In this case A$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$A$.
4:     sigma – double scalar
The standard error, i.e., the value sqrt(rTr / (mk))$\sqrt{{r}^{\mathrm{T}}r/\left(m-k\right)}$ when m > k$m>k$, and the value zero when m = k$m=k$. Here r$r$ is the residual vector bAx$b-Ax$ and k$k$ is the rank of A$A$.
5:     irank – int64int32nag_int scalar
k$k$, the rank of the matrix A$A$. It should be noted that it is possible for irank to be returned as n$n$ and svd to be returned as true. This means that the matrix U$U$ only just failed the test for nonsingularity.
6:     work(lwork) – double array
If svd is returned as false, then the first n$n$ elements of work contain information on the QU$QU$ factorization of A$A$ (see parameter a above and nag_eigen_real_gen_qu_svd (f02wd)), and work(n + 1)${\mathbf{work}}\left(n+1\right)$ contains the condition number UEU1E${‖U‖}_{E}{‖{U}^{-1}‖}_{E}$ of the upper triangular matrix U$U$.
If svd is returned as true, then the first n$n$ elements of work contain the singular values of A$A$ arranged in descending order and work(n + 1)${\mathbf{work}}\left(n+1\right)$ contains the total number of iterations taken by the QR$QR$ algorithm. The rest of work is used as workspace.
7:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n < 1${\mathbf{n}}<1$, or m < n${\mathbf{m}}<{\mathbf{n}}$, or lda < m$\mathit{lda}<{\mathbf{m}}$, or lwork < 4 × n${\mathbf{lwork}}<4×{\mathbf{n}}$.
ifail = 2${\mathbf{ifail}}=2$
The QR$QR$ algorithm has failed to converge to the singular values in 50 × n$50×{\mathbf{n}}$ iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.

## Accuracy

The computed factors Q$Q$, U$U$, R$R$, D$D$ and PT${P}^{\mathrm{T}}$ satisfy the relations
Q
 ( U ) 0
= A + E,Q
 ( R 0 ) 0 I
 ( D ) 0
PT = A + F,
$Q U 0 =A+E,Q R 0 0 I D 0 PT=A+F,$
where
 ‖E‖2 ≤ c1ε ‖A‖2, $‖E‖2≤c1ε ‖A‖2,$
 ‖F‖2 ≤ c2ε‖A‖2, $‖F‖2≤c2ε‖A‖2,$
ε$\epsilon$ being the machine precision, and c1${c}_{1}$ and c2${c}_{2}$ being modest functions of m$m$ and n$n$. Note that A2 = σ1${‖A‖}_{2}={\sigma }_{1}$.
For a fuller discussion, covering the accuracy of the solution x$x$ see Lawson and Hanson (1974), especially pages 50 and 95.

If the least squares solution is obtained from the QU$QU$ factorization then the time taken by the function is approximately proportional to n2(3mn)${n}^{2}\left(3m-n\right)$. If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to n2(3m + 19n)${n}^{2}\left(3m+19n\right)$. The approximate proportionality factor is the same in each case.
This function is column biased and so is suitable for use in paged environments.
Following the QU$QU$ factorization of A$A$ the condition number
 c(U) = ‖U‖E ‖U − 1‖E $c(U)=‖U‖E ‖U-1‖E$
is determined and if c(U)$c\left(U\right)$ is such that
 c(U) × tol > 1.0 $c(U)×tol>1.0$
then U$U$ is regarded as singular and the singular values of A$A$ are computed. If this test is not satisfied, U$U$ is regarded as nonsingular and the rank of A$A$ is set to n$n$. When the singular values are computed the rank of A$A$, say k$k$, is returned as the largest integer such that
 σk > tol × σ1, $σk>tol×σ1,$
unless σ1 = 0${\sigma }_{1}=0$ in which case k$k$ is returned as zero. That is, singular values which satisfy σitol × σ1${\sigma }_{i}\le {\mathbf{tol}}×{\sigma }_{1}$ are regarded as negligible because relative perturbations of order tol can make such singular values zero.

## Example

```function nag_linsys_real_gen_solve_example
a = [0.05, 0.05, 0.25, -0.25;
0.25, 0.25, 0.05, -0.05;
0.35, 0.35, 1.75, -1.75;
1.75, 1.75, 0.35, -0.35;
0.3, -0.3, 0.3, 0.3;
0.4, -0.4, 0.4, 0.4];
b = [1;
2;
3;
4;
5;
6];
tol = 0.0005;
lwork = int64(32);
[aOut, bOut, svd, sigma, irank, work, ifail] = nag_linsys_real_gen_solve(a, b, tol, lwork)
```
```

aOut =

-0.5000   -0.5000   -0.5000    0.5000
0.5000    0.5000   -0.5000    0.5000
-0.5000    0.5000   -0.5000   -0.5000
0.5000   -0.5000   -0.5000   -0.5000
0.1604   -0.5793   -0.0473    0.1714
0.2138   -0.7724   -0.0630    0.2228

bOut =

4.9667
-2.8333
4.5667
3.2333
1.1258
0.8292

svd =

1

sigma =

0.9092

irank =

3

work =

3.0000
2.0000
1.0000
0.0000
1.0000
-1.7333
0.4000
-7.8000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ifail =

0

```
```function f04jg_example
a = [0.05, 0.05, 0.25, -0.25;
0.25, 0.25, 0.05, -0.05;
0.35, 0.35, 1.75, -1.75;
1.75, 1.75, 0.35, -0.35;
0.3, -0.3, 0.3, 0.3;
0.4, -0.4, 0.4, 0.4];
b = [1;
2;
3;
4;
5;
6];
tol = 0.0005;
lwork = int64(32);
[aOut, bOut, svd, sigma, irank, work, ifail] = f04jg(a, b, tol, lwork)
```
```

aOut =

-0.5000   -0.5000   -0.5000    0.5000
0.5000    0.5000   -0.5000    0.5000
-0.5000    0.5000   -0.5000   -0.5000
0.5000   -0.5000   -0.5000   -0.5000
0.1604   -0.5793   -0.0473    0.1714
0.2138   -0.7724   -0.0630    0.2228

bOut =

4.9667
-2.8333
4.5667
3.2333
1.1258
0.8292

svd =

1

sigma =

0.9092

irank =

3

work =

3.0000
2.0000
1.0000
0.0000
1.0000
-1.7333
0.4000
-7.8000
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ifail =

0

```