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)

▸▿ Contents

1 Purpose

2 Syntax

3 Description

4 References

▸▿ 5 Parameters

5.1 Compulsory Input Parameters

5.2 Optional Input Parameters

5.3 Output Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

Purpose

nag_matop_real_gen_pseudinv (f01bl) calculates the rank and pseudo-inverse of an

m

n

real matrix,

m \geq n

, using a

Q R

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:

At Mark 22:

m was made optional

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

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 Further Comments).
2: $a (lda, n)$ – double array: lda, the first dimension of the array, must satisfy the constraint $lda \geq m$ .
The $m$ by $n$ rectangular matrix $A$ .

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.
$m$ and $n$ , the number of rows and columns in the matrix $A$ .

Constraint: $m \geq n$ .

Output Parameters

1: $a (lda, n)$ – double array: The transpose of the pseudo-inverse of $A$ .
2: $aijmax (n)$ – double array: $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 Further Comments).
3: $irank$ – int64int32nag_int scalar: $r$ , the rank of $A$ 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 = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$: Inverse not found, due to an incorrect determination of irank (see Further Comments).

$ifail = 2$

Invalid tolerance, due to

(i)	t is negative, $irank = - 1$ ;
(ii)	t too large, $irank = 0$ ;
(iii)	t too small, $irank > 0$ .

$ifail = 3$

On entry,

m < n

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

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.

Further Comments

The time taken by nag_matop_real_gen_pseudinv (f01bl) 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 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

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)

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.

Open in the MATLAB editor: f01bl_example

function f01bl_example


fprintf('f01bl example results\n\n');

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];

anorm = norm(a);
t = anorm*eps('double');

[ainv, aijmax, irank, inc, ifail] = f01bl(t, a);

fprintf('Maximum element in reduced matrix at Kth stages\n');
fprintf('   K      Modulus\n');
ind = [1:double(irank)+1];
fprintf('%4d   %12.4e\n',[ind ; aijmax']);
fprintf('\nRank = %5d\n\n', irank);
fprintf('Pseudo-inverse\n');
disp(ainv');

f01bl example results

Maximum element in reduced matrix at Kth stages
   K      Modulus
   1     9.0000e+00
   2     9.3101e+00
   3     8.7461e+00
   4     5.6832e+00
   5     2.8449e-16

Rank =     4

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

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox