NAG FL Interface
f04qaf (real_​gen_​sparse_​lsqsol)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

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

2 Specification

Fortran Interface
Subroutine f04qaf ( m, n, b, x, se, aprod, damp, atol, btol, conlim, itnlim, msglvl, itn, anorm, acond, rnorm, arnorm, xnorm, work, ruser, lruser, iuser, liuser, inform, ifail)
Integer, Intent (In) :: m, n, msglvl, lruser, liuser
Integer, Intent (Inout) :: itnlim, iuser(liuser), ifail
Integer, Intent (Out) :: itn, inform
Real (Kind=nag_wp), Intent (In) :: damp, atol, btol, conlim
Real (Kind=nag_wp), Intent (Inout) :: b(m), ruser(lruser)
Real (Kind=nag_wp), Intent (Out) :: x(n), se(n), anorm, acond, rnorm, arnorm, xnorm, work(n,2)
External :: aprod
C Header Interface
#include <nag.h>
void  f04qaf_ (const Integer *m, const Integer *n, double b[], double x[], double se[],
void (NAG_CALL *aprod)(Integer *mode, const Integer *m, const Integer *n, double x[], double y[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser),
const double *damp, const double *atol, const double *btol, const double *conlim, Integer *itnlim, const Integer *msglvl, Integer *itn, double *anorm, double *acond, double *rnorm, double *arnorm, double *xnorm, double work[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser, Integer *inform, Integer *ifail)
The routine may be called by the names f04qaf or nagf_linsys_real_gen_sparse_lsqsol.

3 Description

f04qaf can be used to solve a system of linear equations
Ax=b (1)
where A is an n×n sparse nonsymmetric matrix, or can be used to solve linear least squares problems, so that f04qaf minimizes the value ρ given by
ρ=r,  r=b-Ax (2)
where A is an m×n sparse matrix and r denotes the Euclidean length of r so that r2=rTr. A damping argument, λ, may be included in the least squares problem in which case f04qaf minimizes the value ρ given by
ρ2=r2+λ2x2. (3)
λ is supplied as the argument damp and should of course be zero if the solution to problems (1) or (2) is required. Minimizing ρ in (3) is often called ridge regression.
f04qaf 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 routine does not require A explicitly, but A is specified via aprod which must perform the operations (y+Ax) and (x+ATy) for a given n-element vector x and m element vector y. A argument to aprod specifies which of the two operations is required on a given entry.
The routine also returns estimates of the standard errors of the sample regression coefficients (xi, for i=1,2,,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=s2(ATA)−1,  s2=ρ2/(m-n),  m>n  
and when problem (3) is being solved then V is given by
V=s2(ATA+λ2I) −1,  s2=ρ2/m,  λ0.  
Let A¯ denote the matrix
A¯=A,  λ=0;  A¯=( A λI ) ,  λ0, (4)
let r¯ denote the residual vector
r¯=r,  λ=0;  r¯=( b 0 )-A¯x,  λ0 (5)
corresponding to an iterate x, so that ρ=r¯ is the function being minimized, and let A denote the Frobenius (Euclidean) norm of A. Then the routine accepts x as a solution if it is estimated that one of the following two conditions is satisfied:
ρtol1A¯.x+tol2b (6)
A¯Tr¯ tol1 A¯ ρ (7)
where tol1 and tol2 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,
A¯T 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 routine 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 A¯ be
σ1σ2σn0  
then the condition number of A¯ is defined as
cond(A¯)=σ1/σk  
where σk is the smallest nonzero singular value of A¯ and hence k is the rank of A¯. When k<n, then A¯ is rank deficient, the least squares solution is not unique and f04qaf will normally converge to the minimal length solution. In practice A¯ will not have exactly zero singular values, but may instead have small singular values that we wish to regard as zero.
The routine provides for this possibility by terminating if
cond(A¯)clim (8)
where clim is a user-supplied limit on the condition number of 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 r, ATr and x, which are all returned by the routine, will indicate whether or not an acceptable solution has been found. Condition (8), perhaps in conjunction with λ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 λ tends to reduce the size of the computed solution and to make its components less sensitive to changes in the data, and f04qaf is applicable when a value of λ is known a priori. To have an effect, λ should normally be at least εA 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.

4 References

Paige C C and Saunders M A (1982a) 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 (1982b) Algorithm 583 LSQR: Sparse linear equations and least squares problems ACM Trans. Math. Software 8 195–209

5 Arguments

1: m Integer Input
On entry: m, the number of rows of the matrix A.
Constraint: m1.
2: n Integer Input
On entry: n, the number of columns of the matrix A.
Constraint: n1.
3: b(m) Real (Kind=nag_wp) array Input/Output
On entry: the right-hand side vector b.
On exit: b is overwritten.
4: x(n) Real (Kind=nag_wp) array Output
On exit: the solution vector x.
5: se(n) Real (Kind=nag_wp) array Output
On exit: the estimates of the standard errors of the components of x. Thus se(i) contains an estimate of νii, where νii is the ith 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.
6: aprod Subroutine, supplied by the user. External Procedure
aprod must perform the operations yy+Ax and xx+ATy for given vectors x and y.
The specification of aprod is:
Fortran Interface
Subroutine aprod ( mode, m, n, x, y, ruser, lruser, iuser, liuser)
Integer, Intent (In) :: m, n, lruser, liuser
Integer, Intent (Inout) :: mode, iuser(liuser)
Real (Kind=nag_wp), Intent (Inout) :: x(n), y(m), ruser(lruser)
C Header Interface
void  aprod (Integer *mode, const Integer *m, const Integer *n, double x[], double y[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser)
1: mode Integer Input/Output
On entry: specifies which operation is to be performed.
mode=1
aprod must compute y+Ax.
mode=2
aprod must compute x+ATy.
On exit: may be used as a flag to indicate a failure in the computation of y+Ax or x+ATy. If mode is negative on exit from aprod, f04qaf will exit immediately with ifail set to mode.
2: m Integer Input
On entry: m, the number of rows of A.
3: n Integer Input
On entry: n, the number of columns of A.
4: x(n) Real (Kind=nag_wp) array Input/Output
On entry: the vector x.
On exit: if mode=1, x must be unchanged.
If mode=2, x must contain x+ATy.
5: y(m) Real (Kind=nag_wp) array Input/Output
On entry: the vector y.
On exit: if mode=1, y must contain y+Ax.
If mode=2, y must be unchanged.
6: ruser(lruser) Real (Kind=nag_wp) array User Workspace
7: lruser Integer Input
8: iuser(liuser) Integer array User Workspace
9: liuser Integer Input
aprod is called with the arguments ruser, lruser, iuser and liuser as supplied to f04qaf. You should use the arrays ruser, lruser, iuser and liuser to supply information to aprod.
aprod must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which f04qaf is called. Arguments denoted as Input must not be changed by this procedure.
Note: aprod should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by f04qaf. If your code inadvertently does return any NaNs or infinities, f04qaf is likely to produce unexpected results.
7: damp Real (Kind=nag_wp) Input
On entry: the value λ. If either problem (1) or problem (2) is to be solved, damp must be supplied as zero.
8: atol Real (Kind=nag_wp) Input
On entry: the tolerance, tol1, 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×10−4. If atol is supplied as less than ε, where ε is the machine precision, the value ε is used instead of atol.
9: btol Real (Kind=nag_wp) Input
On entry: the tolerance, tol2, 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×10−4. If btol is supplied as less than ε, the value ε is used instead of btol.
10: conlim Real (Kind=nag_wp) Input
On entry: the value clim of equation (8); it should be an upper limit on the condition number of 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.
11: itnlim Integer Input/Output
On entry: an upper limit on the number of iterations. If itnlim0, the value n is used in place of itnlim, but for ill-conditioned problems a higher value of itnlim is likely to be necessary.
On exit: unchanged unless itnlim0 on entry, in which case it is set to n.
12: msglvl Integer Input
On entry: the level of printing from f04qaf. If msglvl0, then no printing occurs, but otherwise messages will be output on the current advisory message unit (see x04abf). A description of the printed output is given in Section 9.1. The level of printing is determined as follows:
msglvl0
No printing.
msglvl=1
A brief summary is printed just prior to return from f04qaf.
msglvl2
A summary line is printed periodically to monitor the progress of f04qaf, together with a brief summary just prior to return from f04qaf.
13: itn Integer Output
On exit: the number of iterations performed.
14: anorm Real (Kind=nag_wp) Output
On exit: an estimate of A¯ for the matrix A¯ of (4).
15: acond Real (Kind=nag_wp) Output
On exit: an estimate of cond(A¯) which is a lower bound.
16: rnorm Real (Kind=nag_wp) Output
On exit: an estimate of r¯ for the residual, r¯, of (5) corresponding to the solution x returned in x. Note that r¯ is the function being minimized.
17: arnorm Real (Kind=nag_wp) Output
On exit: an estimate of the A¯Tr¯ corresponding to the solution x returned in x.
18: xnorm Real (Kind=nag_wp) Output
On exit: an estimate of x for the solution x returned in x.
19: work(n,2) Real (Kind=nag_wp) array Workspace
20: ruser(lruser) Real (Kind=nag_wp) array User Workspace
ruser is not used by f04qaf, but is passed directly to aprod and may be used to pass information to this routine.
21: lruser Integer Input
On entry: the dimension of the array ruser as declared in the (sub)program from which f04qaf is called.
22: iuser(liuser) Integer array User Workspace
iuser is not used by f04qaf, but is passed directly to aprod and may be used to pass information to this routine.
23: liuser Integer Input
On entry: the dimension of the array iuser as declared in the (sub)program from which f04qaf is called.
24: inform Integer Output
On exit: the reason for termination of f04qaf.
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 tol1 and tol2 as the values supplied in atol and btol respectively.
inform=2
The termination criterion of (7) has been satisfied with tol1 as the value supplied in atol.
inform=3
The termination criterion of (6) has been satisfied with tol1 and/or tol2 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 tol1 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 and 7 correspond to failure with ifail=2, 3 and 4 respectively (see Section 6) and when ifail is negative inform will be set to the same negative value.
25: ifail Integer Input/Output
On entry: ifail must be set to 0, −1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of −1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value −1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, m=value.
Constraint: m1.
On entry, n=value.
Constraint: n1.
ifail=2
Termination criteria not satisfied, but the condition number bound supplied in conlim has been met.
The condition of (8) has been satisfied for the value of clim supplied in conlim. If this failure is unexpected you should check that aprod is working correctly. Although conditions (6) or (7) have not been satisfied, the values returned in rnorm, arnorm and xnorm may nevertheless indicate that an acceptable solution has been reached.
ifail=3
Termination criteria not satisfied, but the condition number is bounded by 1/x02ajf.
The condition of (8) has been satisfied for the value clim=1.0/ε, where ε is the machine precision. The matrix A¯ is nearly singular or rank deficient and the problem is too ill-conditioned for this machine. If this failure is unexpected, you should check that aprod is working correctly.
ifail=4
Maximum number of iterations reached.
The limit on the number of iterations has been reached. The number of iterations required by f04qaf and the condition of the matrix A¯ can depend strongly on the scaling of the problem. Poor scaling of the rows and columns of A should be avoided whenever possible.
ifail<0
mode in aprod has been set to: value.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

When the problem is compatible, the computed solution x will satisfy the equation
r=b-Ax,  
where an estimate of r is returned in the argument rnorm. When the problem is incompatible, the computed solution x will satisfy the equation
A¯T r¯ = e ,  
where an estimate of e is returned in the argument arnorm. See also Section 6.2 of Paige and Saunders (1982b).

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
f04qaf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The time taken by f04qaf is likely to be principally determined by the time taken in aprod, which is called twice on each iteration, once with mode=1 and once with mode=2. The time taken per iteration by the remaining operations in f04qaf is approximately proportional to max(m,n).
The Lanczos process will usually converge more quickly if A is pre-conditioned by a nonsingular matrix M that approximates A in some sense and is also chosen so that equations of the form My=c can efficiently be solved for y. For a discussion of preconditioning, see the F11 Chapter Introduction. In the context of f04qaf, problem (1) is equivalent to
(AM-1)y=b,  Mx=y  
and problem (2) is equivalent to minimizing
ρ=r,   r=b-(AM-1)y,   Mx=y.  
Note that the normal matrix (AM-1)T(AM-1)=M-T(ATA)M-1 so that the preconditioning AM-1 is equivalent to the preconditioning M-T(ATA)M-1 of the normal matrix ATA.
Pre-conditioning can be incorporated into f04qaf simply by coding aprod to compute y+AM-1x and x+M-TATy in place of y+Ax and x+ATy respectively, and then solving the equations Mx=y for x on return from f04qaf. The quantity y+AM-1x should be computed by solving Mz=x for z and then computing y+Az, and x+M-TATy should be computed by solving MTz=ATy for z and then forming x+z.

9.1 Description of the Printed Output

When msglvl>0, then f04qaf will produce output (except in the case where the routine fails with ifail=1) on the advisory message channel (see x04abf).
When msglvl2 then a summary line is printed periodically giving the following information:
Output Meaning
ITN Iteration number, k.
X(1) The first element of the current iterate xk.
FUNCTION The current value of the function, ρ, being minimized.
COMPAT An estimate of r¯k/b, where r¯k is the residual corresponding to xk. This value should converge to zero (in theory) if and only if the problem is compatible. COMPAT decreases monotonically.
INCOMPAT An estimate of A¯Tr¯k / (A¯r¯k) which should converge to zero if and only if at the solution ρ is nonzero. INCOMPAT is not usually monotonic.
NRM(ABAR) A monotonically increasing estimate of A¯.
COND(ABAR) A monotonically increasing estimate of the condition number cond(A¯).

10 Example

This example solves the linear least squares problem
minρ =r,  r=b-Ax  
where A is the 13×12 matrix and b is the 13 element vector given by
A=( 1 0 0 −1 0 0 0 0 0 0 0 0 0 1 0 0 −1 0 0 0 0 0 0 0 0 0 1 −1 0 0 0 0 0 0 0 0 −1 0 −1 4 −1 0 0 −1 0 0 0 0 0 −1 0 −1 4 −1 0 0 −1 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 −1 0 0 0 0 0 0 0 −1 0 0 −1 4 −1 0 −1 0 0 0 0 0 −1 0 0 −1 4 −1 0 −1 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 −1 0 0 1 0 0 0 0 0 0 0 0 0 −1 0 0 1 1 1 1 0 0 1 1 0 0 1 1 1 ) ,  b=-h2 ( 0 0 0 1 1 0 0 1 1 0 0 0 -h−3 )  
with h=0.1.
Such a problem can arise by considering the Neumann problem on a rectangle
δu δn =0 δu δn =0 2 u=g(x,y) δu δn =0 c u=1 δu δn =0  
where C is the boundary of the rectangle, and discretizing as illustrated below with the square mesh
Figure 1
Figure 1
The 12×12 symmetric part of A represents the difference equations and the final row comes from the normalizing condition. The example program has g(x,y)=1 at all the internal mesh points, but apart from this is written in a general manner so that the number of rows (NROWS) and columns (NCOLS) in the grid can readily be altered.

10.1 Program Text

Program Text (f04qafe.f90)

10.2 Program Data

Program Data (f04qafe.d)

10.3 Program Results

Program Results (f04qafe.r)