hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_matop_real_gen_pseudinv (f01bl)

Purpose

nag_matop_real_gen_pseudinv (f01bl) calculates the rank and pseudo-inverse of an mm by nn real matrix, mnmn, using a QRQR factorization with column interchanges.

Syntax

[a, aijmax, irank, inc, ifail] = f01bl(t, a, 'm', m, 'n', n)
[a, aijmax, irank, inc, ifail] = nag_matop_real_gen_pseudinv(t, a, 'm', m, 'n', n)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: m has been made optional
.

Description

Householder's factorization with column interchanges is used in the decomposition F = QUF=QU, where FF is AA with its columns permuted, QQ is the first rr columns of an mm by mm orthogonal matrix and UU is an rr by nn upper-trapezoidal matrix of rank rr. The pseudo-inverse of FF is given by XX where
X = UT(UUT)1QT.
X=UT(UUT)-1QT.
If the matrix is found to be of maximum rank, r = nr=n, UU is a nonsingular nn by nn upper-triangular matrix and the pseudo-inverse of FF simplifies to X = U1QTX=U-1QT. The transpose of the pseudo-inverse of AA is overwritten on AA.

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

Parameters

Compulsory Input Parameters

1:     t – double scalar
The tolerance used to decide when elements can be regarded as zero (see Section [Further Comments]).
2:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldamldam.
The mm by nn rectangular matrix AA.

Optional Input Parameters

1:     m – int64int32nag_int scalar
2:     n – int64int32nag_int scalar
Default: For m, the first dimension of the array a. The second dimension of the array a.
mm and nn, the number of rows and columns in the matrix AA.
Constraint: mnmn.

Input Parameters Omitted from the MATLAB Interface

lda d u ldu du

Output Parameters

1:     a(lda,n) – double array
ldamldam.
The transpose of the pseudo-inverse of AA.
2:     aijmax(n) – double array
aijmax(i)aijmaxi contains the element of largest modulus in the reduced matrix at the iith stage. If r < nr<n, then only the first r + 1r+1 elements of aijmax have values assigned to them; the remaining elements are unused. The ratio aijmax(1) / aijmax(r)aijmax1/aijmaxr usually gives an indication of the condition number of the original matrix (see Section [Further Comments]).
3:     irank – int64int32nag_int scalar
rr, the rank of AA as determined using the tolerance t.
4:     inc(n) – int64int32nag_int array
The record of the column interchanges in the Householder factorization.
5:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
Inverse not found, due to an incorrect determination of irank (see Section [Further Comments]).
  ifail = 2ifail=2
Invalid tolerance, due to
(i) t is negative, irank = 1irank=-1;
(ii) t too large, irank = 0irank=0;
(iii) t too small, irank > 0irank>0.
  ifail = 3ifail=3
On entry,m < nm<n.

Accuracy

For most matrices the pseudo-inverse is the best possible having regard to the condition of AA 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 QQ and UU satisfy the relation QU = F + EQU=F+E where
E2 < cεA2 + η×sqrt((mr)(nr))
E2<cεA2+η(m-r)(n-r)
in which cc is a modest function of mm and nn, ηη is the value of t, and εε is the machine precision.

Further Comments

The time taken by nag_matop_real_gen_pseudinv (f01bl) is approximately proportional to mnrmnr.
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 nag_lapack_dgesvd (f08kb)). In nag_matop_real_gen_pseudinv (f01bl) 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 nn times the bound on possible errors in individual elements of the original matrix. If the elements of AA 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 AA is 10p10p we expect to get pp figures wrong in the pseudo-inverse. An estimate of the condition number is usually given by aijmax(1) / aijmax(r)aijmax1/aijmaxr.

Example

function nag_matop_real_gen_pseudinv_example
t = 2.99772935794301e-15;
a = [7, -2, 4, 9, 1.8;
     3, 8, -4, 6, 1.3;
     9, 6, 1, 5, 2.1;
     -8, 7, 5, 2, 0.6;
     4, -1, 2, 8, 1.3;
     1, 6, 3, -5, 0.5];
[aOut, aijmax, irank, inc, ifail] = nag_matop_real_gen_pseudinv(t, a)
 

aOut =

    0.0178   -0.0216    0.0520    0.0237    0.0072
   -0.0118    0.0434   -0.0813    0.0357   -0.0014
    0.0472    0.0294    0.0139   -0.0138    0.0077
   -0.0566    0.0291    0.0474    0.0305    0.0050
   -0.0037   -0.0138    0.0166    0.0357    0.0035
    0.0384    0.0343    0.0576   -0.0571    0.0073


aijmax =

    9.0000
    9.3101
    8.7461
    5.6832
    0.0000


irank =

                    4


inc =

                    4
                    2
                    4
                    4
                    0


ifail =

                    0


function f01bl_example
t = 2.99772935794301e-15;
a = [7, -2, 4, 9, 1.8;
     3, 8, -4, 6, 1.3;
     9, 6, 1, 5, 2.1;
     -8, 7, 5, 2, 0.6;
     4, -1, 2, 8, 1.3;
     1, 6, 3, -5, 0.5];
[aOut, aijmax, irank, inc, ifail] = f01bl(t, a)
 

aOut =

    0.0178   -0.0216    0.0520    0.0237    0.0072
   -0.0118    0.0434   -0.0813    0.0357   -0.0014
    0.0472    0.0294    0.0139   -0.0138    0.0077
   -0.0566    0.0291    0.0474    0.0305    0.0050
   -0.0037   -0.0138    0.0166    0.0357    0.0035
    0.0384    0.0343    0.0576   -0.0571    0.0073


aijmax =

    9.0000
    9.3101
    8.7461
    5.6832
    0.0000


irank =

                    4


inc =

                    4
                    2
                    4
                    4
                    0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013