PDF version (NAG web site
, 64bit version, 64bit version)
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=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 function also returns the value of the standard error
${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 righthand 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\times {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\times {\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 righthand 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
$mn$ 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(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$.
 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
${\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.
 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\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.
 ${\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
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.
Further Comments
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(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 function 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.
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}$.
Open in the MATLAB editor:
f04jg_example
function f04jg_example
fprintf('f04jg example results\n\n');
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
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015