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$, 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.

## 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$ is the vector $x$ of minimum (Euclidean) length which minimizes the length of the residual vector $r=b-Ax$.
The real $m$ by $n\left(m\ge n\right)$ matrix $A$ is factorized as
 $A=Q R 0$
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
 $x= R-1 0 QTb.$
If $R$ is not of full rank, then the singular value decomposition of $R$ is obtained so that $R$ is factorized as
 $R=UDVT,$
where $U$ and $V$ are $n$ by $n$ orthogonal matrices and $D$ is the $n$ by $n$ diagonal matrix
 $D=diagσ1,σ2,…,σn,$
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
 $x=V S-1 0 0 0 UT 0 0 I QTb,$
where $S=\mathrm{diag}\left({\sigma }_{1},{\sigma }_{2},\dots ,{\sigma }_{k}\right)$.
The function also returns the value of the standard error
 $σ = rTr m-k , if ​m>k, = 0, if ​m=k,$
${r}^{\mathrm{T}}r$ being the residual sum of squares and $k$ the rank of $A$.

## References

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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge {\mathbf{m}}$.
The $m$ by $n$ matrix $A$.
2:     $\mathrm{b}\left({\mathbf{m}}\right)$ – double array
The right-hand side vector $b$.
3:     $\mathrm{tol}$ – double scalar
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×{10}^{-4}$. See Further Comments 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, then the value $\epsilon$ is used in place of tol. For most problems this is unreasonably small.
4:     $\mathrm{lwork}$int64int32nag_int scalar
The dimension of the array work.
Constraint: ${\mathbf{lwork}}\ge 4×{\mathbf{n}}$.

### Optional Input Parameters

1:     $\mathrm{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$, the number of rows of a.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
2:     $\mathrm{n}$int64int32nag_int scalar
Default: the second dimension of the array a.
$n$, the number of columns of ${\mathbf{a}}$.
Constraint: $1\le {\mathbf{n}}\le {\mathbf{m}}$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
If svd is returned as false, a stores details of the $QR$ factorization of $A$.
If svd is returned as true, the first $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:     $\mathrm{b}\left({\mathbf{m}}\right)$ – double array
The first $n$ elements of b contain the minimal least squares solution vector $x$. The remaining $m-n$ elements are used for workspace.
3:     $\mathrm{svd}$ – logical scalar
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$.
4:     $\mathrm{sigma}$ – double scalar
The standard error, i.e., the value $\sqrt{{r}^{\mathrm{T}}r/\left(m-k\right)}$ when $m>k$, and the value zero when $m=k$. Here $r$ is the residual vector $b-Ax$ and $k$ is the rank of $A$.
5:     $\mathrm{irank}$int64int32nag_int scalar
$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.
6:     $\mathrm{work}\left({\mathbf{lwork}}\right)$ – double array
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 ${‖R‖}_{E}{‖{R}^{-1}‖}_{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.
7:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
Constraint: $\mathit{lda}\ge {\mathbf{m}}$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
Constraint: ${\mathbf{n}}\ge 1$.
On entry, lwork is too small.
${\mathbf{ifail}}=2$
The $QR$ algorithm has failed to converge to the singular values in $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.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

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

If the least squares solution is obtained from the $QR$ factorization then the time taken by the function is approximately proportional to ${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 ${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 $QR$ factorization of $A$ the condition number
 $cR=RE R-1E$
is determined and if $c\left(R\right)$ is such that
 $cR×tol>1.0$
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
 $σk>tol×σ1,$
unless ${\sigma }_{1}=0$ in which case $k$ is returned as zero. That is, singular values which satisfy ${\sigma }_{i}\le {\mathbf{tol}}×{\sigma }_{1}$ are regarded as negligible because relative perturbations of order tol can make such singular values zero.

## Example

This example obtains a least squares solution for $r=b-Ax$, where
 $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.30 -0.30 0.30 0.30 0.40 -0.40 0.40 0.40 , b= 1 2 3 4 5 6$
and the value tol is to be taken as $5×{10}^{-4}$.
function f04jg_example

fprintf('f04jg example results\n\n');

% Solve linear least squares problem Ax = b for general A
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);
[a, x, svd, sigma, irank, work, ifail] = ...
f04jg(a, b, tol, lwork);

disp('Solution');
disp(x(1:4)');
fprintf('Standard error = %6.3f, rank = %2d\n\n', sigma, irank);
if svd==1
fprintf('solution obtained from SVD of A\n');
else
fprintf('solution obtained from QU factorization of A\n');
end

f04jg example results

Solution
4.9667   -2.8333    4.5667    3.2333

Standard error =  0.909, rank =  3

solution obtained from SVD of A