naginterfaces.library.linsys.real_​gen_​sparse_​lsqsol

naginterfaces.library.linsys.real_gen_sparse_lsqsol(n, b, aprod, damp=0.0, atol=- 1.0, btol=- 1.0, conlim=0.0, itnlim=0, msglvl=0, data=None, io_manager=None)

real_gen_sparse_lsqsol solves sparse nonsymmetric equations, sparse linear least squares problems and sparse damped linear least squares problems, using a Lanczos algorithm.

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

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

Parameters

nint

$n$, the number of columns of the matrix $A$.

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

The right-hand side vector $b$.

aprodcallable (mode, x, y) = aprod(mode, x, y, data=None)

$aprod$ must perform the operations $y:=y+Ax$ and $x:=x+A^{T}y$ for given vectors $x$ and $y$.

Parameters

modeint

Specifies which operation is to be performed.

$mode=1$

$aprod$ must compute $y+Ax$.

$mode=2$

$aprod$ must compute $x+A^{T}y$.

xfloat, ndarray, shape $\left(n\right)$

The vector $x$.

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

The vector $y$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

modeint

May be used as a flag to indicate a failure in the computation of $y+Ax$ or $x+A^{T}y$. If $mode$ is negative on exit from $aprod$, real_gen_sparse_lsqsol will exit immediately with $\text{errno}$ set to $mode$.

xfloat, array-like, shape $\left(n\right)$

If $mode=1$, $x$ must be unchanged.

If $mode=2$, $x$ must contain $x+A^{T}y$.

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

If $mode=1$, $y$ must contain $y+Ax$.

If $mode=2$, $y$ must be unchanged.

dampfloat, optional

The value $\lambda$. If either problem (1) or problem (2) is to be solved, $damp$ must be supplied as zero.

atolfloat, optional

The tolerance, ${\text{tol}}_1$, of the convergence criteria (6) and (7); it should be an estimate of the largest relative error in the elements of $A$. For example, if the elements of $A$ are correct to about $4$ significant figures, then $atol$ should be set to about $5\times {10}^{-4}$. If $atol$ is supplied as less than $ϵ$, where $ϵ$ is the machine precision, the value $ϵ$ is used instead of $atol$.

btolfloat, optional

The tolerance, ${\text{tol}}_2$, of the convergence criterion (6); it should be an estimate of the largest relative error in the elements of $B$. For example, if the elements of $B$ are correct to about $4$ significant figures, then $btol$ should be set to about $5\times {10}^{-4}$. If $btol$ is supplied as less than $ϵ$, the value $ϵ$ is used instead of $btol$.

conlimfloat, optional

The value $c_{lim}$ of equation (8); it should be an upper limit on the condition number of $\overline{A}$. $conlim$ should not normally be chosen much larger than $1.0/atol$. If $conlim$ is supplied as zero, the value $1.0/ϵ$ is used instead of $conlim$.

itnlimint, optional

An upper limit on the number of iterations. If $itnlim\le 0$, the value $n$ is used in place of $itnlim$, but for ill-conditioned problems a higher value of $itnlim$ is likely to be necessary.

msglvlint, optional

The level of printing from real_gen_sparse_lsqsol. If $msglvl\le 0$, then no printing occurs, but otherwise messages will be output on the file object associated with the advisory I/O unit (see FileObjManager). A description of the printed output is given in Further Comments. The level of printing is determined as follows:

$msglvl\le 0$

No printing.

$msglvl=1$

A brief summary is printed just prior to return from real_gen_sparse_lsqsol.

$msglvl\ge 2$

A summary line is printed periodically to monitor the progress of real_gen_sparse_lsqsol, together with a brief summary just prior to return from real_gen_sparse_lsqsol.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

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

$b$ is overwritten.

xfloat, ndarray, shape $\left(n\right)$

The solution vector $x$.

sefloat, ndarray, shape $\left(n\right)$

The estimates of the standard errors of the components of $x$. Thus $se[i-1]$ contains an estimate of $\sqrt{{\nu }_{ii}}$, where ${\nu }_{ii}$ is the $i$th diagonal element of the estimated variance-covariance matrix $V$. The estimates returned in $se$ will be the lower bounds on the actual estimated standard errors, but will usually be correct to at least one significant figure.

itnlimint

Unchanged unless $itnlim\le 0$ on entry, in which case it is set to $n$.

itnint

The number of iterations performed.

anormfloat

An estimate of $\Vert \overline{A}\Vert$ for the matrix $\overline{A}$ of (4).

acondfloat

An estimate of $cond\left(\overline{A}\right)$ which is a lower bound.

rnormfloat

An estimate of $\parallel \overline{r}\parallel$ for the residual, $\overline{r}$, of (5) corresponding to the solution $x$ returned in $x$. Note that $\parallel \overline{r}\parallel$ is the function being minimized.

arnormfloat

An estimate of the $\Vert {\overline{A}}^T\overline{r}\Vert$ corresponding to the solution $x$ returned in $x$.

xnormfloat

An estimate of $\parallel x\parallel$ for the solution $x$ returned in $x$.

informint

The reason for termination of real_gen_sparse_lsqsol.

$inform=0$

The exact solution is $x=0$. No iterations are performed in this case.

$inform=1$

The termination criterion of (6) has been satisfied with ${\text{tol}}_1$ and ${\text{tol}}_2$ as the values supplied in $atol$ and $btol$ respectively.

$inform=2$

The termination criterion of (7) has been satisfied with ${\text{tol}}_1$ as the value supplied in $atol$.

$inform=3$

The termination criterion of (6) has been satisfied with ${\text{tol}}_1$ and/or ${\text{tol}}_2$ as the value $ϵ$, where $ϵ$ is the machine precision. One or both of the values supplied in $atol$ and $btol$ must have been less than $ϵ$ and was too small for this machine.

$inform=4$

The termination criterion of (7) has been satisfied with ${\text{tol}}_1$ as the value $ϵ$. The value supplied in $atol$ must have been less than $ϵ$ and was too small for this machine.

The values $inform=5$, $6$ or $7$ correspond to failure with $errno$ = 2, 3 and 4 respectively (see Exceptions) and when $\text{errno}$ is negative $inform$ will be set to the same negative value.

Raises

NagValueError

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 1$.

(errno $1$)

On entry, $m=⟨value⟩$.

Constraint: $m\ge 1$.

(errno $2$)

Termination criteria not satisfied, but the condition number bound supplied in $conlim$ has been met.

(errno $3$)

Termination criteria not satisfied, but the condition number is bounded by $1$/machine.precision.

(errno $4$)

Maximum number of iterations reached.

Warns

NagAlgorithmicWarning

(errno $i<0$)

$mode$ in $aprod$ has been set to: $⟨value⟩$.

Notes

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

real_gen_sparse_lsqsol can be used to solve a system of linear equations

$$Ax=b$$

where $A$ is an $n\times n$ sparse nonsymmetric matrix, or can be used to solve linear least squares problems, so that real_gen_sparse_lsqsol minimizes the value $\rho$ given by

$$\rho =\parallel r\parallel \text{,}r=b-Ax$$

where $A$ is an $m\times n$ sparse matrix and $\parallel r\parallel$ denotes the Euclidean length of $r$ so that ${\parallel r\parallel }^2=r^{T}r$. A damping argument, $\lambda$, may be included in the least squares problem in which case real_gen_sparse_lsqsol minimizes the value $\rho$ given by

$${\rho }^2={\parallel r\parallel }^2+{\lambda }^2{\parallel x\parallel }^2\text{.}$$

$\lambda$ is supplied as the argument $damp$ and should of course be zero if the solution to problems (1) and (2) is required. Minimizing $\rho$ in (3) is often called ridge regression.

real_gen_sparse_lsqsol is based upon algorithm LSQR (see Paige and Saunders (1982a) and Paige and Saunders (1982b)) and solves the problems by an algorithm based upon the Lanczos process. The function does not require $A$ explicitly, but $A$ is specified via $aprod$ which must perform the operations $(y+Ax)$ and $(x+A^{T}y)$ for a given $n$-element vector $x$ and $m$ element vector $y$. An argument to $aprod$ specifies which of the two operations is required on a given entry.

The function also returns estimates of the standard errors of the sample regression coefficients ( $x_{\text{i}}$, for $\text{i}=1,2,\dots ,n$) given by the diagonal elements of the estimated variance-covariance matrix $V$. When problem (2) is being solved and $A$ is of full rank, then $V$ is given by

$$V=s^2{\left(A^{T}A\right)}^{-1}\text{,}{s}^2={\rho }^2/(m-n)\text{,}m>n$$

and when problem (3) is being solved then $V$ is given by

$$V=s^2{(A^{T}A+{\lambda }^{2}I)}^{-1}\text{,}{s}^2={\rho }^2/m\text{,}\lambda \ne 0\text{.}$$

Let $\overline{A}$ denote the matrix

$$\begin{array}{c}\overline{A}=A\text{,}\lambda =0\text{;}\overline{A}=\left(\begin{array}{c}A\\ \lambda I\end{array}\right)\text{,}\lambda \ne 0\text{,}\end{array}$$

let $\overline{r}$ denote the residual vector

$$\begin{array}{c}\overline{r}=r\text{,}\lambda =0\text{;}\overline{r}=\left(\begin{array}{c}b\\ 0\end{array}\right)-\overline{A}x\text{,}\lambda \ne 0\end{array}$$

corresponding to an iterate $x$, so that $\rho =\parallel \overline{r}\parallel$ is the function being minimized, and let $\parallel A\parallel$ denote the Frobenius (Euclidean) norm of $A$. Then the function accepts $x$ as a solution if it is estimated that one of the following two conditions is satisfied:

$$\rho \le {\text{tol}}_1\Vert \overline{A}\Vert .\parallel x\parallel +{\text{tol}}_2\parallel b\parallel$$

$$\Vert {\overline{A}}^T\overline{r}\Vert \le {\text{tol}}_1\Vert \overline{A}\Vert \rho$$

where ${\text{tol}}_1$ and ${\text{tol}}_2$ are user-supplied tolerances which estimate the relative errors in $A$ and $b$ respectively. Condition (6) is appropriate for compatible problems where, in theory, we expect the residual to be zero and will be satisfied by an acceptable solution $x$ to a compatible problem. Condition (7) is appropriate for incompatible systems where we do not expect the residual to be zero and is based on the observation that, in theory,

$${\overline{A}}^T\overline{r}=0$$

when $x$ is a solution to the least squares problem, and so (7) will be satisfied by an acceptable solution $x$ to a linear least squares problem.

The function also includes a test to prevent convergence to solutions, $x$, with unacceptably large elements. This can happen if $A$ is nearly singular or is nearly rank deficient. If we let the singular values of $\overline{A}$ be

$${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_n\ge 0$$

then the condition number of $\overline{A}$ is defined as

$$cond\left(\overline{A}\right)={\sigma }_1/{\sigma }_k$$

where ${\sigma }_k$ is the smallest nonzero singular value of $\overline{A}$ and hence $k$ is the rank of $\overline{A}$. When $k<n$, then $\overline{A}$ is rank deficient, the least squares solution is not unique and real_gen_sparse_lsqsol will normally converge to the minimal length solution. In practice $\overline{A}$ will not have exactly zero singular values, but may instead have small singular values that we wish to regard as zero.

The function provides for this possibility by terminating if

$$cond\left(\overline{A}\right)\ge c_{lim}$$

where $c_{lim}$ is a user-supplied limit on the condition number of $\overline{A}$. For problem (1) termination with this condition indicates that $A$ is nearly singular and for problem (2) indicates that $A$ is nearly rank deficient and so has near linear dependencies in its columns. In this case inspection of $\parallel r\parallel$, $\Vert A^{T}r\Vert$ and $\parallel x\parallel$, which are all returned by the function, will indicate whether or not an acceptable solution has been found. Condition (8), perhaps in conjunction with $\lambda \ne 0$, can be used to try and ‘regularize’ least squares solutions. A full discussion of the stopping criteria is given in Section 6 of Paige and Saunders (1982a).

Introduction of a nonzero damping argument $\lambda$ tends to reduce the size of the computed solution and to make its components less sensitive to changes in the data, and real_gen_sparse_lsqsol is applicable when a value of $\lambda$ is known a priori. To have an effect, $\lambda$ should normally be at least $\sqrt{ϵ}\parallel A\parallel$ where $ϵ$ is the machine precision. For further discussion see Paige and Saunders (1982b) and the references given there.

Whenever possible the matrix $A$ should be scaled so that the relative errors in the elements of $A$ are all of comparable size. Such a scaling helps to prevent the least squares problem from being unnecessarily sensitive to data errors and will normally reduce the number of iterations required. At the very least, in the absence of better information, the columns of $A$ should be scaled to have roughly equal column length.

References

Paige, C C and Saunders, M A, 1982, LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Software (8), 43–71

Paige, C C and Saunders, M A, 1982, Algorithm 583 LSQR: Sparse linear equations and least squares problems, ACM Trans. Math. Software (8), 195–209