NAG FL Interface
f01blf (real_gen_pseudinv)
1
Purpose
f01blf calculates the rank and pseudoinverse of an $m$ by $n$ real matrix, $m\ge n$, using a $QR$ factorization with column interchanges.
2
Specification
Fortran Interface
Subroutine f01blf ( 
m, n, t, a, lda, aijmax, irank, inc, d, u, ldu, du, ifail) 
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) 

C++ Header Interface
#include <nag.h> extern "C" {
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=QU$, where
$F$ is
$A$ with its columns permuted,
$Q$ is the first
$r$ columns of an
$m$ by
$m$ orthogonal matrix and
$U$ is an
$r$ by
$n$ uppertrapezoidal matrix of rank
$r$. The pseudoinverse of
$F$ is given by
$X$ where
If the matrix is found to be of maximum rank,
$r=n$,
$U$ is a nonsingular
$n$ by
$n$ uppertriangular matrix and the pseudoinverse of
$F$ simplifies to
$X={U}^{1}{Q}^{\mathrm{T}}$. The transpose of the pseudoinverse of
$A$ is overwritten on
$A$.
4
References
Peters G and Wilkinson J H (1970) The least squares problem and pseudoinverses Comput. J. 13 309–316
Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag
5
Arguments

1:
$\mathbf{m}$ – Integer
Input

2:
$\mathbf{n}$ – Integer
Input

On entry: $m$ and $n$, the number of rows and columns in the matrix $A$.
Constraint:
${\mathbf{m}}\ge {\mathbf{n}}$.

3:
$\mathbf{t}$ – Real (Kind=nag_wp)
Input

On entry: the tolerance used to decide when elements can be regarded as zero (see
Section 9).

4:
$\mathbf{a}\left({\mathbf{lda}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) array
Input/Output

On entry: the $m$ by $n$ rectangular matrix $A$.
On exit: the transpose of the pseudoinverse of $A$.

5:
$\mathbf{lda}$ – Integer
Input

On entry: the first dimension of the array
a as declared in the (sub)program from which
f01blf is called.
Constraint:
${\mathbf{lda}}\ge {\mathbf{m}}$.

6:
$\mathbf{aijmax}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) array
Output

On exit:
${\mathbf{aijmax}}\left(i\right)$ 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
${\mathbf{aijmax}}\left(1\right)/{\mathbf{aijmax}}\left(r\right)$ usually gives an indication of the condition number of the original matrix (see
Section 9).

7:
$\mathbf{irank}$ – Integer
Output

On exit:
$r$, the rank of
$A$ as determined using the tolerance
t.

8:
$\mathbf{inc}\left({\mathbf{n}}\right)$ – Integer array
Output

On exit: the record of the column interchanges in the Householder factorization.

9:
$\mathbf{d}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) array
Workspace

10:
$\mathbf{u}\left({\mathbf{ldu}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) array
Workspace


11:
$\mathbf{ldu}$ – Integer
Input

On entry: the first dimension of the array
u as declared in the (sub)program from which
f01blf is called.
Constraint:
${\mathbf{ldu}}\ge {\mathbf{n}}$.

12:
$\mathbf{du}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) array
Workspace


13:
$\mathbf{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 $\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{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:
 ${\mathbf{ifail}}=1$

Inverse not found. Incorrect
irank.
 ${\mathbf{ifail}}=2$

Invalid tolerance, ${\mathbf{t}}<0.0$: ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
Invalid tolerance,
t too large:
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=3$

On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}\ge {\mathbf{n}}$.
 ${\mathbf{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.
 ${\mathbf{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.
 ${\mathbf{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 pseudoinverse 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
$QU=F+E$ where
in which
$c$ is a modest function of
$m$ and
$n$,
$\eta $ is the value of
t, and
$\epsilon $ is the
machine precision.
8
Parallelism and Performance
f01blf is not threaded in any implementation.
The time taken by f01blf is approximately proportional to $mnr$.
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$ is ${10}^{p}$ we expect to get $p$ figures wrong in the pseudoinverse. An estimate of the condition number is usually given by ${\mathbf{aijmax}}\left(1\right)/{\mathbf{aijmax}}\left(r\right)$.
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 pseudoinverse. Data and results are given for an example which is a
$6$ by
$5$ matrix of deficient rank in which the last column is a linear combination of the other four. Setting
t to
$\epsilon $ times the norm of the matrix, the rank is correctly determined as
$4$ and the pseudoinverse is computed to full implementation accuracy.
10.1
Program Text
10.2
Program Data
10.3
Program Results