naginterfaces.library.linsys.real_​gen_​solve

naginterfaces.library.linsys.real_gen_solve(a, b, tol, lwork)

real_gen_solve finds the solution of a linear least squares problem, $Ax=b$, where $A$ is a real $m\times n(m\ge n)$ 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.

For full information please refer to the NAG Library document for f04jg

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/f04/f04jgf.html

Parameters

afloat, array-like, shape $(m,n)$

The $m\times n$ matrix $A$.

bfloat, array-like, shape $\left(m\right)$

The right-hand side vector $b$.

tolfloat

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 $(ϵ,1.0)$, where $ϵ$ is the machine precision, the value $ϵ$ is used in place of $tol$. For most problems this is unreasonably small.

lworkint

The dimension of the array $work$.

Returns

afloat, ndarray, shape $(m,n)$

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 right-hand singular vectors, stored by rows; and the remaining rows of the array are used as workspace.

bfloat, ndarray, shape $\left(m\right)$

The first $n$ elements of $b$ contain the minimal least squares solution vector $x$. The remaining $m-n$ elements are used for workspace.

svdbool

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

sigmafloat

The standard error, i.e., the value $\sqrt{r^{T}r/(m-k)}$ 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$.

irankint

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

workfloat, ndarray, shape $\left(lwork\right)$

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 $work\left[n\right]$ contains the condition number ${\parallel R\parallel }_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 $work\left[n\right]$ contains the total number of iterations taken by the $QR$ algorithm.

The rest of $work$ is used as workspace.

Raises

NagValueError

(errno $1$)

On entry, $lwork$ is too small. Minimum size required: $⟨value⟩$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 1$.

(errno $1$)

On entry, $m=⟨value⟩$ and $n=⟨value⟩$.

Constraint: $m\ge n$.

(errno $2$)

The $QR$ algorithm has failed to converge to the singular values in $50\times n$ iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

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\times n(m\ge n)$ matrix $A$ is factorized as

$$\begin{array}{c}A=Q\left(\begin{array}{c}R\\ 0\end{array}\right)\end{array}$$

where $Q$ is an $m\times m$ orthogonal matrix and $R$ is an $n\times n$ upper triangular matrix. If $R$ is of full rank, then the least squares solution is given by

$$x=\left(\begin{array}{cc}{R}^{-1}& 0\end{array}\right)Q^{T}b\text{.}$$

If $R$ is not of full rank, then the singular value decomposition of $R$ is obtained so that $R$ is factorized as

$$R=UD{V}^T\text{,}$$

where $U$ and $V$ are $n\times n$ orthogonal matrices and $D$ is the $n\times n$ diagonal matrix

$$D=diag({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n)\text{,}$$

with ${\sigma }_1\ge {\sigma }_2\ge \cdots \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

$$\begin{array}{c}x=V\left(\begin{array}{cc}{S}^{-1}& 0\\ 0& 0\end{array}\right)\left(\begin{array}{cc}{U}^T& 0\\ 0& I\end{array}\right)Q^{T}b\text{,}\end{array}$$

where $S=diag({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_k)$.

The function also returns the value of the standard error

$$\begin{array}{c}\begin{array}{cccc}\sigma & =& \sqrt{\frac{r^{T}r}{m-k}}\text{,}& \text{if}m>k\text{,}\\ & & & \\ & =& 0\text{,}& \text{if}m=k\text{,}\end{array}\end{array}$$

$r^{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