NAG Library Routine Document

F04JGF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

F04JGF finds the solution of a linear least squares problem,

A x = b

, where

A

is a real

m

n (m \geq 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.

2 Specification

SUBROUTINE F04JGF (

M, N, A, LDA, B, TOL, SVD, SIGMA, IRANK, WORK, LWORK, IFAIL)

INTEGER	M, N, LDA, IRANK, LWORK, IFAIL
REAL (KIND=nag_wp)	A(LDA,N), B(M), TOL, SIGMA, WORK(LWORK)
LOGICAL	SVD

3 Description

The minimal least squares solution of the problem

A x = b

is the vector

x

of minimum (Euclidean) length which minimizes the length of the residual vector

r = b - A x

The real

m

n (m \geq n)

matrix

A

is factorized as

A = Q (\begin{array}{l} U \\ 0 \end{array})

where

Q

is an

m

m

orthogonal matrix and

U

is an

n

n

upper triangular matrix. If

U

is of full rank, then the least squares solution is given by

x = (U^{- 1} 0) Q^{T} b .

U

is not of full rank, then the singular value decomposition of

U

is obtained so that

U

is factorized as

U = R D P^{T},

where

R

and

P

are

n

n

orthogonal matrices and

D

is the

n

n

diagonal matrix

D = diag (σ_{1}, σ_{2}, \dots, σ_{n}),

with

σ_{1} \geq σ_{2} \geq \dots \geq σ_{n} \geq 0

, these being the singular values of

A

. If the singular values

σ_{k + 1}, \dots, σ_{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 = P (\begin{array}{l} S^{- 1} & 0 \\ 0 & 0 \end{array}) (\begin{array}{l} R^{T} & 0 \\ 0 & I \end{array}) Q^{T} b,

where

S = diag (σ_{1}, σ_{2}, \dots, σ_{k})

This routine obtains the factorizations by a call to F02WDF.

The routine also returns the value of the standard error

\begin{array}{ccl} σ & = & \sqrt{\frac{r^{T} r}{m - k}}, & if ​ m > k, \\ = & 0, & if ​ m = k, \end{array}

r^{T} r

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 Parameters

1: M – INTEGERInput: On entry: $m$ , the number of rows of A.
Constraint: $M \geq N$ .
2: N – INTEGERInput: On entry: $n$ , the number of columns of $A$ .
Constraint: $1 \leq N \leq M$ .
3: A(LDA,N) – REAL (KIND=nag_wp) arrayInput/Output: On entry: the $m$ by $n$ matrix $A$ .

On exit: if SVD is returned as .FALSE., A is overwritten by details of the $Q U$ factorization of $A$ (see F02WDF for further details).
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 – INTEGERInput: On entry: the first dimension of the array A as declared in the (sub)program from which F04JGF is called.
Constraint: $LDA \geq M$ .
5: B(M) – REAL (KIND=nag_wp) arrayInput/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 \times 10^{- 4}$ . See Section 8 for a description of how TOL is used to determine rank. If TOL is outside the range $(ε, 1.0)$ , where $ε$ is the machine precision, then the value $ε$ is used in place of TOL. For most problems this is unreasonably small.
7: SVD – LOGICALOutput: On exit: is returned as .FALSE. if the least squares solution has been obtained from the $Q U$ 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 $\sqrt{r^{T} r / (m - k)}$ when $m > k$ , and the value zero when $m = k$ . Here $r$ is the residual vector $b - A x$ and $k$ is the rank of $A$ .
9: IRANK – INTEGEROutput: 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 $U$ only just failed the test for nonsingularity.
10: WORK(LWORK) – REAL (KIND=nag_wp) arrayOutput: On exit: if SVD is returned as .FALSE., then the first $n$ elements of WORK contain information on the $Q U$ factorization of $A$ (see parameter A above and F02WDF), and $WORK (n + 1)$ contains the condition number ${‖U‖}_{E} {‖U^{- 1}‖}_{E}$ of the upper triangular matrix $U$ .
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 $Q R$ algorithm. The rest of WORK is used as workspace.
11: LWORK – INTEGERInput: On entry: the dimension of the array WORK as declared in the (sub)program from which F04JGF is called.
Constraint: $LWORK \geq 4 \times N$ .
12: IFAIL – INTEGERInput/Output: On entry: IFAIL must be set to $0$ , $- 1 or 1$ . If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $- 1 or 1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is $0$ . 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

- 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,	$N < 1$ ,
or	$M < N$ ,
or	$LDA < M$ ,
or	$LWORK < 4 \times N$ .

$IFAIL = 2$: The $Q R$ 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.

7 Accuracy

The computed factors

Q

U

R

D

and

P^{T}

satisfy the relations

Q (\begin{array}{l} U \\ 0 \end{array}) = A + E, Q (\begin{array}{l} R & 0 \\ 0 & I \end{array}) (\begin{array}{l} D \\ 0 \end{array}) P^{T} = A + F,

where

{‖E‖}_{2} \leq c_{1} ε {‖A‖}_{2},

{‖F‖}_{2} \leq c_{2} ε {‖A‖}_{2},

ε

being the machine precision, and

c_{1}

and

c_{2}

being modest functions of

m

and

n

. Note that

{‖A‖}_{2} = σ_{1}

For a fuller discussion, covering the accuracy of the solution

x

see Lawson and Hanson (1974), especially pages 50 and 95.

8 Further Comments

If the least squares solution is obtained from the

Q U

factorization then the time taken by the routine is approximately proportional to

n^{2} (3 m - n)

. If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to

n^{2} (3 m + 19 n)

. 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

Q U

factorization of

A

the condition number

c (U) = {‖U‖}_{E} {‖U^{- 1}‖}_{E}

is determined and if

c (U)

is such that

c (U) \times TOL > 1.0

then

U

is regarded as singular and the singular values of

A

are computed. If this test is not satisfied,

U

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 \times σ_{1},

unless

σ_{1} = 0

in which case

k

is returned as zero. That is, singular values which satisfy

σ_{i} \leq TOL \times σ_{1}

are regarded as negligible because relative perturbations of order TOL can make such singular values zero.

9 Example

This example obtains a least squares solution for

r = b - A x

, where

A = (\begin{array}{r} 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 \end{array}), b = (\begin{array}{l} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{array})

and the value TOL is to be taken as

5 \times 10^{- 4}

NAG Library Routine DocumentF04JGF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

F04JGF