NAG FL Interface
f04jgf (real_​gen_​solve)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

f04jgf finds the solution of a linear least squares problem, Ax=b, where A is a real m×n(mn) 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.

2 Specification

Fortran Interface
Subroutine f04jgf ( m, n, a, lda, b, tol, svd, sigma, irank, work, lwork, ifail)
Integer, Intent (In) :: m, n, lda, lwork
Integer, Intent (Inout) :: ifail
Integer, Intent (Out) :: irank
Real (Kind=nag_wp), Intent (In) :: tol
Real (Kind=nag_wp), Intent (Inout) :: a(lda,n), b(m)
Real (Kind=nag_wp), Intent (Out) :: sigma, work(lwork)
Logical, Intent (Out) :: svd
C Header Interface
#include <nag.h>
void  f04jgf_ (const Integer *m, const Integer *n, double a[], const Integer *lda, double b[], const double *tol, logical *svd, double *sigma, Integer *irank, double work[], const Integer *lwork, Integer *ifail)
The routine may be called by the names f04jgf or nagf_linsys_real_gen_solve.

3 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=b-Ax.
The real m×n(mn) matrix A is factorized as
A=Q ( R 0 )  
where Q is an m×m orthogonal matrix and R is an n×n upper triangular matrix. If R is of full rank, then the least squares solution is given by
x= ( R-1 0 ) QTb.  
If R is not of full rank, then the singular value decomposition of R is obtained so that R is factorized as
R=UDVT,  
where U and V are n×n orthogonal matrices and D is the n×n diagonal matrix
D=diag(σ1,σ2,,σn),  
with σ1σ2σn0, these being the singular values of A. If the singular values σk+1,,σn are negligible, but σ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
x=V ( S-1 0 0 0 ) ( UT 0 0 I ) QTb,  
where S=diag(σ1,σ2,,σk).
The routine also returns the value of the standard error
σ = rTr m-k , if ​m>k, = 0, if ​m=k,  
rTr being the residual sum of squares and k the rank of A.

4 References

Lawson C L and Hanson R J (1974) Solving Least Squares Problems Prentice–Hall

5 Arguments

1: m Integer Input
On entry: m, the number of rows of a.
Constraint: mn.
2: n Integer Input
On entry: n, the number of columns of a.
Constraint: 1nm.
3: a(lda,n) Real (Kind=nag_wp) array Input/Output
On entry: the m×n matrix A.
On exit: 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.
4: lda Integer Input
On entry: the first dimension of the array a as declared in the (sub)program from which f04jgf is called.
Constraint: ldam.
5: b(m) Real (Kind=nag_wp) array Input/Output
On entry: the right-hand side vector b.
On exit: the first n elements of b contain the minimal least squares solution vector x. The remaining m-n elements are used for workspace.
6: tol Real (Kind=nag_wp) Input
On entry: 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×10−4. See Section 9 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.
7: svd Logical Output
On exit: 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.
8: sigma Real (Kind=nag_wp) Output
On exit: the standard error, i.e., the value rTr/(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.
9: irank Integer Output
On exit: 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.
10: work(lwork) Real (Kind=nag_wp) array Output
On exit: 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(n+1) contains the condition number RER-1E 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(n+1) contains the total number of iterations taken by the QR algorithm. The rest of work is used as workspace.
11: lwork Integer Input
On entry: the dimension of the array work as declared in the (sub)program from which f04jgf is called.
Constraint: lwork4×n.
12: 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, lda=value and m=value.
Constraint: ldam.
On entry, lwork is too small. Minimum size required: value.
On entry, m=value and n=value.
Constraint: mn.
On entry, n=value.
Constraint: n1.
ifail=2
The QR algorithm has failed to converge to the singular values in 50×n iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.
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

The computed factors Q, R, U, D and VT satisfy the relations
Q ( R 0 )=A+E,  
Q ( U 0 0 I ) ( D 0 ) VT=A+F,  
where
E2c1ε A2,  
F2c2εA2,  
ε being the machine precision, and c1 and c2 being modest functions of m and n. Note that A2=σ1.
For a fuller discussion, covering the accuracy of the solution x see Lawson and Hanson (1974), especially pages 50 and 95.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
f04jgf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f04jgf 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

If the least squares solution is obtained from the QR factorization then the time taken by the routine is approximately proportional to n2(3m-n). If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to n2(3m+19n). The approximate proportionality factor is the same in each case.
This routine is column biased and so is suitable for use in paged environments.
Following the QR factorization of A the condition number
c(R)=RE R-1E  
is determined and if c(R) is such that
c(R)×tol>1.0  
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
σk>tol×σ1,  
unless σ1=0 in which case k is returned as zero. That is, singular values which satisfy σitol×σ1 are regarded as negligible because relative perturbations of order tol can make such singular values zero.

10 Example

This example obtains a least squares solution for r=b-Ax, where
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.30 -0.30 0.30 0.30 0.40 -0.40 0.40 0.40 ) ,  b=( 1 2 3 4 5 6 )  
and the value tol is to be taken as 5×10−4.

10.1 Program Text

Program Text (f04jgfe.f90)

10.2 Program Data

Program Data (f04jgfe.d)

10.3 Program Results

Program Results (f04jgfe.r)