Linear Least-Squares Systems

Solution of Linear Systems, Linear Least-Squares Systems
Iterative Methods for Least-Squares: Conjugate Gradient

Get these three pieces of code!
LSQRF90_INT.F90
LSQRF90.F90
LSQRF90_EX.F90

Compile these in the order shown. This order is important since Fortran modules are used. Missing routines are found in the IMSL Fortran 90 MP Library product.

Solving Linear Equations and Least-Squares Systems
The conjugate gradient routine LSQRF90 is used to solve large sparse least-squares systems. These are represented by Axb where A is an m x n coefficient data matrix, b is a given right-hand side m-vector, and x is the solution n-vector being computed. There are no restrictions on the relative sizes of m and n. When A is sparse, this algorithm can substantially reduce computer time and storage requirements, compared with using a routine for solving a general system.

The user provides code to compute Av and A^Tu for given n-vectors v and m-vectors u. This routine uses forward or reverse communication, so A can be in any convenient storage scheme. In the document that factors RELPR = epsilon(0D0), a Fortran 90 elemental function giving the largest relative floating point spacing.

Solve a real least-squares linear system using a conjugate gradient method with forward or reverse communication. A preconditioner can also be used.

CALL LSQRF90 (M,N,DAMP,U,V,X,ISTOP, &
   IMODE,APROD,ITNLIM,ATOL,BTOL,CONLIM,NOUT,SE,DOPE)

Required Arguments

M	(Input) Number of linear equations
N	(Input) Number of unknowns
DAMP	(Input) The damping (regularizing, ridge-regression or Levenberg-Marquardt) parameter. If the system Axb is incompatible, a value of DAMP in the interval from 0 to sqrt( RELPR) x ||A|| is reasonable. Increasing the value of DAMP decreases the solution vector length, ||x|| , increases the residual vector length, ||b - Ax||,and tends to reduce the number of iterations required by LSQRF90. The work per iteration and the storage required are the same for all values of DAMP.
U(1:M)	(Input/Output) Assumed size array of length M initially containing the right-hand side, b. This array is used for computing matrix-vector products, U <--- U + A*V either internal to the APROD subroutine or during reverse communication. It is over-written during this process.
V (1:N)	(Input/Output) Assumed size array of length N. This array is used for computing matrix products, V <--- V + transpose(A )*U either internal to the APROD subroutine or during reverse communication.
X(1:N)	(Output) Assumed size array of length N containing the solution, x . The solution is defined only after all iterations are completed. This is an issue with reverse communication, so X should not be used in that case until IMODE=-1.
ISTOP	(Output) The reason for termination: 0 The exact solution is x=0. No iterations were performed. 1 The equations Axb are probably compatible. The value ||b-Ax|| is sufficiently small, given the values of ATOL and BTOL. 2 The system Axb is probably not compatible. A least-squares solution has been obtained which is sufficiently accurate, given the value of ATOL. 3 The matrix C= has a condition number that exceeds CONLIM. The system Axb appears to be ill-conditioned. This might indicate an error in the evaluation of the matrix-vector products in the user's code. 4 The equations Ax = b are probably compatible. The value ||b - Ax|| is as small as seems reasonable on this machine. 5 The system Ax = b is probably not compatible. A least-squares solution has been obtained which is as accurate as seems reasonable on this machine. 6 The condition number of Cseems to be so large that there is not much point in doing further iterations, given the precision of the arithmetic. This might indicate an error in the evaluation of the matrix-vector products in the user's code. 7 The iteration limit ITNLIM was reached.

Optional Arguments

Note: Exactly one of IMODE and APROD arguments must be used.

IMODE	(Input/Output) This is the communication direction flag for Reverse Communication, where the matrix-vector products are provided in the calling program unit instead of the subroutine APROD. When IMODE is present, always initialize IMODE = -2. The alternate use of APROD is known as Forward Communication. Exactly one of the arguments IMODE or APROD must be present when solving a system. Construct a loop around the call to LSQRF90 that responds to requests for matrix-vector products as indicated. IMODE=-2 SOLVE: DO call LSQRF90 (M,N,...,IMODE=IMODE,...) IF(IMODE < 0) EXIT SOLVE IF IMODE = 1, COMPUTE U = U + A*V and CYCLE SOLVE IF IMODE = 2, COMPUTE V = V + TRANSPOSE(A )*U and CYCLE SOLVE. END DO SOLVE Note: The loop can be EXITED whenever desired. To properly shut down the process and set the output parameters, reset ITNLIM = 1 and re-enter LSQRF90. An exit with IMODE = -1 will follow shortly, and then the loop will exit as programmed.
APROD	(External Subroutine Name) If Reverse Communication is used it is not necessary to write the subroutine APROD. When APROD is present, the matrix-vector products are obtained within LSQRF90 by subroutine calls of the form: CALL APROD( MODE,M,N,V,U) which must perform the following functions: IF MODE = 1, COMPUTE U = U + A*V IF MODE = 2, COMPUTE V = V + TRANSPOSE(A )*U To initiate shut-down inside of LSQRF90 set MODE = 0 after performing the required function. This results in the same effect as setting ITNLIM=1 when using reverse communication. An orderly termination of LSQRF90 will shortly follow. There is a subroutine with these arguments provided. It is a "poison-pill" routine, APROD. If a user does not write a routine of their own and uses the optional argument 'APROD = APROD', a fatal error message will be printed when this routine is called inside of LSQRF90.
ITNLIM	(Input) Upper limit for the number of iterations permitted Suggested value: = N/2 for well-conditioned systems, and
= 4*N	otherwise

Algorithm

Subroutine LSQRF90 solves the least-squares linear system Ax b using the conjugate gradient method as developed for least-squares by Paige and Saunders. This is summarized by Åke Björck, Numerical Methods for Least-Squares Problems, SIAM Publications, Philadelphia, PA, (1996), pages 306-309. The routine LSQRF90 was originally published as Algorithm 583, "LSQR: Sparse linear equations and least-squares," ACM Trans. Math. Software, (1982), pages 195-209. We provide a more flexible user interface with fewer arguments in the calling sequence.

The preconditioner matrix, R , approximates some factor of A = QR, where Q is orthogonal. This matrix is generally chosen so that the linear systems Rz = xare easy to solve. These properties are in conflict. Whatever preconditioner matrix is chosen, the operations for the matrix-vector products are not much different from the case where a preconditioner is not used. The basic steps are illustrated by considering the problem with the conditioned matrix AR^-1. Instead of directly computing this matrix, the same effect is achieved in the user's interface code by the steps:

To compute u <-- u AR^-1v (MODE=1),
Solve Rz = v for z. Then u <-- u + Az.

To compute v <-- v + (AR^-1)Tu, (MODE=2),
set z = u^TA and solve z^T=w^TR for w. Then v <-- v + w.

Following the last matrix-vector product or iteration, compute x <--R^-1x. For systems with large condition numbers, the use of a preconditioner may be essential for computing a solution.

Example

The program LSQRF90_EX illustrates several uses of routine LSQRF90 including an example using a diagonal matrix preconditioner. Values of M and N are input by the user. Then a random rectangular matrix and a right-hand side are generated. Matrix-vector products are provided to the routine LSQRF90 using Reverse Communication and then Forward Communication. In Forward Communication, with the example subroutine MYPROD, the matrix data is communicated using a Fortran 90 module which gives an associated pointer to the matrix data array of the main program unit. This obviates the need for COMMON to communicate data to subroutines.

A final example illustrates reverse communication and sparse matrices represented in "coordinate form." Thus A is equivalent to the list {i, j, a_ij}, including nz non-zero a_ij and the corresponding row and column indices. We construct this matrix by selecting entries of the random rectangular matrix with values between 1/4 and 3/4.