Integer, Intent (In)	::	m, n, lda, ldu
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	irank, inc(n)
Real (Kind=nag_wp), Intent (In)	::	t
Real (Kind=nag_wp), Intent (Inout)	::	a(lda,n), u(ldu,n)
Real (Kind=nag_wp), Intent (Out)	::	aijmax(n), d(m), du(n)

C Header Interface

#include <nag.h>

void	f01blf_ (const Integer m, const Integer n, const double t, double a[], const Integer lda, double aijmax[], Integer irank, Integer inc[], double d[], double u[], const Integer ldu, double du[], Integer *ifail)

The routine may be called by the names f01blf or nagf_matop_real_gen_pseudinv.

3 Description

Householder's factorization with column interchanges is used in the decomposition

F = Q U

, where

F

A

with its columns permuted,

Q

is the first

r

columns of an

m

m

orthogonal matrix and

U

is an

r

n

upper-trapezoidal matrix of rank

r

. The pseudo-inverse of

F

is given by

X

where

X = U^{T} {(U U^{T})}^{- 1} Q^{T} .

If the matrix is found to be of maximum rank,

r = n

U

is a nonsingular

n

n

upper-triangular matrix and the pseudo-inverse of

F

simplifies to

X = U^{- 1} Q^{T}

. The transpose of the pseudo-inverse of

A

is overwritten on

A

4 References

Peters G and Wilkinson J H (1970) The least squares problem and pseudo-inverses Comput. J. 13 309–316

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

5 Arguments

1: $m$ – Integer Input
2: $n$ – Integer Input: On entry: $m$ and $n$ , the number of rows and columns in the matrix $A$ .

Constraint: $m \geq n$ .
3: $t$ – Real (Kind=nag_wp) Input: On entry: the tolerance used to decide when elements can be regarded as zero (see Section 9).
4: $a (lda, n)$ – Real (Kind=nag_wp) array Input/Output: On entry: the $m$ by $n$ rectangular matrix $A$ .

On exit: the transpose of the pseudo-inverse of $A$ .
5: $lda$ – Integer Input: On entry: the first dimension of the array a as declared in the (sub)program from which f01blf is called.

Constraint: $lda \geq m$ .
6: $aijmax (n)$ – Real (Kind=nag_wp) array Output: On exit: $aijmax (i)$ contains the element of largest modulus in the reduced matrix at the $i$ th stage. If $r < n$ , then only the first $r + 1$ elements of aijmax have values assigned to them; the remaining elements are unused. The ratio $aijmax (1) / aijmax (r)$ usually gives an indication of the condition number of the original matrix (see Section 9).
7: $irank$ – Integer Output: On exit: $r$ , the rank of $A$ as determined using the tolerance t.
8: $inc (n)$ – Integer array Output: On exit: the record of the column interchanges in the Householder factorization.
9: $d (m)$ – Real (Kind=nag_wp) array Workspace
10: $u (ldu, n)$ – Real (Kind=nag_wp) array Workspace
11: $ldu$ – Integer Input: On entry: the first dimension of the array u as declared in the (sub)program from which f01blf is called.

Constraint: $ldu \geq n$ .
12: $du (n)$ – Real (Kind=nag_wp) array Workspace
13: $ifail$ – Integer Input/Output: On entry: ifail must be set to $0$ , $- 1 or 1$ . If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface 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 argument, 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$: Inverse not found. Incorrect irank.

$ifail = 2$: Invalid tolerance, $t < 0.0$ : $t = 〈value〉$ .

Invalid tolerance, t too large: $t = 〈value〉$ .

$ifail = 3$: On entry, $m = 〈value〉$ and $n = 〈value〉$ .
Constraint: $m \geq n$ .

$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

For most matrices the pseudo-inverse is the best possible having regard to the condition of

A

and the choice of t. Note that only the singular value decomposition method can be relied upon to give maximum accuracy for the precision of computation used and correct determination of the condition of a matrix (see Wilkinson and Reinsch (1971)).

The computed factors

Q

and

U

satisfy the relation

Q U = F + E

where

{‖E‖}_{2} < c ε {‖A‖}_{2} + η \sqrt{(m - r) (n - r)}

in which

c

is a modest function of

m

and

n

η

is the value of t, and

ε

is the machine precision.

8 Parallelism and Performance

f01blf is not threaded in any implementation.

9 Further Comments

The time taken by f01blf is approximately proportional to

m n r

The most difficult practical problem is the determination of the rank of the matrix (see pages 314–315 of Peters and Wilkinson (1970)); only the singular value decomposition method gives a reliable indication of rank deficiency (see pages 134–151 of Wilkinson and Reinsch (1971) and f08kbf). In f01blf a tolerance, t, is used to recognize ‘zero’ elements in the remaining matrix at each step in the factorization. The value of t should be set at

n

times the bound on possible errors in individual elements of the original matrix. If the elements of

A

vary widely in their orders of magnitude, of course this presents severe difficulties. Sound decisions can only be made by somebody who appreciates the underlying physical problem.

If the condition number of

A

10^{p}

we expect to get

p

figures wrong in the pseudo-inverse. An estimate of the condition number is usually given by

aijmax (1) / aijmax (r)

10 Example

A complete program follows which outputs the maximum of the moduli of the ‘remaining’ elements at each step in the factorization, the rank, as determined by the given value of t, and the transposed pseudo-inverse. Data and results are given for an example which is a

6

5

matrix of deficient rank in which the last column is a linear combination of the other four. Setting t to

ε

times the norm of the matrix, the rank is correctly determined as

4

and the pseudo-inverse is computed to full implementation accuracy.

Interfaces: FL CL AD

NAG FL Interface Introduction

F01 (Matop) Chapter Contents

F01 (Matop) Chapter Introduction

f01bl: FL CL AD

NAG FL Interfacef01blf (real_​gen_​pseudinv)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
f01blf (real_gen_pseudinv)