f01 Chapter Contents
f01 Chapter Introduction (PDF version)
NAG C Library Manual

NAG Library Chapter Introduction

f01 – Matrix Operations, Including Inversion

+− Contents

    1  Scope of the Chapter
+− 2  Background to the Problems
+− 3  Recommendations on Choice and Use of Available Functions
    4  Decision Trees
    5  Functionality Index
    6  Functions Withdrawn or Scheduled for Withdrawal
    7  References

1  Scope of the Chapter

This chapter provides facilities for three types of problem:
(i) Matrix Inversion
(ii) Matrix Factorizations
(iii) Matrix Functions
These problems are discussed separately in Section 2.1, Section 2.2 and Section 2.3.

2  Background to the Problems

2.1  Matrix Inversion

(i) Non-singular square matrices of order n.
If A, a square matrix of order n, is nonsingular (has rank n), then its inverse X exists and satisfies the equations AX=XA=I (the identity or unit matrix).
It is worth noting that if AX-I=R, so that R is the ‘residual’ matrix, then a bound on the relative error is given by R, i.e.,
X-A-1 A-1 ≤R.
(ii) General real rectangular matrices.
A real matrix A has no inverse if it is square (n by n) and singular (has rank <n), or if it is of shape (m by n) with m≠n, but there is a Generalized or Pseudo Inverse Z which satisfies the equations
AZA=A,  ZAZ=Z,  AZT=AZ,  ZAT=ZA
(which of course are also satisfied by the inverse X of A if A is square and nonsingular).
(a) if m≥n and rankA=n then A can be factorized using a QR factorization, given by
A=Q R 0 ,
where Q is an m by m orthogonal matrix and R is an n by n, nonsingular, upper triangular matrix. The pseudo-inverse of A is then given by
Z=R-1Q~T,
where Q~ consists of the first n columns of Q.
(b) if m≤n and rankA=m then A can be factorized using an RQ factorization, given by
A=R 0PT
where P is an n by n orthogonal matrix and R is an m by m, nonsingular, upper triangular matrix. The pseudo-inverse of A is then given by
Z=P~R-1,
where P~ consists of the first m columns of P.
(c) if m≥n and rankA=r≤n then A can be factorized using a QR factorization, with column interchanges, as
A=Q R 0 PT,
where Q is an m by m orthogonal matrix, R is an r by n upper trapezoidal matrix and P is an n by n permutation matrix. The pseudo-inverse of A is then given by
Z=PRTRRT-1Q~T,
where Q~ consists of the first r columns of Q.
(d) if rankA=r≤k=minm,n, then A can be factorized as the singular value decomposition
A=QDPT,
where Q is an m by m orthogonal matrix, P is an n by n orthogonal matrix and D is an m by n diagonal matrix with non-negative diagonal elements. The first k columns of Q and P are the left- and right-hand singular vectors of A respectively and the k diagonal elements of D are the singular values of A. D may be chosen so that
d1≥d2≥⋯≥dk≥0
and in this case if rankA=r then
d1≥d2≥⋯≥dr>0,  dr+1=⋯=dk=0.
If Q~ and P~ consist of the first r columns of Q and P respectively and Σ is an r by r diagonal matrix with diagonal elements d1,d2,…,dr then A is given by
A=Q~ΣP~T
and the pseudo-inverse of A is given by
Z=P~Σ-1Q~T.
Notice that
ATA=PDTDPT
which is the classical eigenvalue (spectral) factorization of ATA.
(e) if A is complex then the above relationships are still true if we use ‘unitary’ in place of ‘orthogonal’ and conjugate transpose in place of transpose. For example, the singular value decomposition of A is
A=QDPH,
where Q and P are unitary, PH the conjugate transpose of P and D is as in (d) above.

2.2  Matrix Factorizations

The functions in this section perform matrix factorizations which are required for the solution of systems of linear equations with various special structures. A few functions which perform associated computations are also included.
Other functions for matrix factorizations are to be found in Chapters f07, f08 and f11.
This section also contains a few functions associated with eigenvalue problems (see Chapter f02). (Historical note: this section used to contain many more such functions, but they have now been superseded by functions in Chapter f08.)

2.3  Matrix Functions

Given a square matrix A, the matrix function fA is a matrix with the same dimensions as A which provides a generalization of the scalar function f.
If A has a full set of eigenvectors V then A can be factorized as
A = V D VT ,
where D is the diagonal matrix whose diagonal elements, di, are the eigenvalues of A. fA is given by
fA = V fD VT ,
where fD is the diagonal matrix whose ith diagonal element is fdi.
In general, A may not have a full set of eigenvectors. The matrix function can then be defined via a Cauchy integral. For A∈ℂn×n,
fA = 1 2π i ∫ Γ fz zI-A-1 dz ,
where Γ is a closed contour surrounding the eigenvalues of A, and f is analytic within Γ.
Algorithms for computing matrix functions are usually tailored to a specific function. Currently Chapter f01 contains routines for calculating the exponential, logarithm, sine, cosine, sinh and cosh of both real and complex matrices. In addition there are routines to compute a general function of real symmetric and complex Hermitian matrices and a general function of general real and complex matrices.

3  Recommendations on Choice and Use of Available Functions

3.1  Matrix Inversion

Note:  before using any function for matrix inversion, consider carefully whether it is really needed.
Although the solution of a set of linear equations Ax=b can be written as x=A-1b, the solution should never be computed by first inverting A and then computing A-1b; the functions in Chapters f04 or f07 should always be used to solve such sets of equations directly; they are faster in execution, and numerically more stable and accurate. Similar remarks apply to the solution of least squares problems which again should be solved by using the functions in Chapters f04 and f08 rather than by computing a pseudo-inverse.
(a) Non-singular square matrices of order n 
This chapter describes techniques for inverting a general real matrix A and matrices which are positive definite (have all eigenvalues positive) and are either real and symmetric or complex and Hermitian. It is wasteful and uneconomical not to use the appropriate function when a matrix is known to have one of these special forms. A general function must be used when the matrix is not known to be positive definite. In most functions the inverse is computed by solving the linear equations Axi=ei, for i=1,2,…,n, where ei is the ith column of the identity matrix.
The residual matrix R=AX-I, where X is a computed inverse of A, conveys useful information. Firstly R is a bound on the relative error in X and secondly R<12  guarantees the convergence of the iterative process in the ‘corrected’ inverse functions.
The decision trees for inversion show which functions in Chapter f04 and Chapter f07 should be used for the inversion of other special types of matrices not treated in the chapter.
(b) General real rectangular matrices
For real matrices nag_dgeqrf (f08aec) returns the QR factorization of the matrix and nag_dgeqp3 (f08bfc) returns the QR factorization with column interchanges. The corresponding complex functions are nag_zgeqrf (f08asc) and nag_zgeqp3 (f08btc) respectively. Functions are also provided to form the orthogonal matrices and transform by the orthogonal matrices following the use of the above functions.
nag_dgesvd (f08kbc) and nag_zgesvd (f08kpc) compute the singular value decomposition as described in Section 2 for real and complex matrices respectively. If A has rank r≤k=minm,n then the k-r smallest singular values will be negligible and the pseudo-inverse of A can be obtained as Z=PΣ-1QT as described in Section 2. If the rank of A is not known in advance it can be estimated from the singular values (see Section 2.4 in the f04 Chapter Introduction). For large sparse matrices, leading terms in the singular value decomposition can be computed using functions from Chapter f12.

3.2  Matrix Factorizations

Each of these functions serves a special purpose required for the solution of sets of simultaneous linear equations or the eigenvalue problem. For further details you should consult Sections 3 or 4 in the f02 Chapter Introduction or Sections 3 or 4 in the f04 Chapter Introduction.
nag_sparse_nsym_fac (f11dac) is provided for factorizing general real sparse matrices. A more recent algorithm for the same problem is available through nag_superlu_lu_factorize (f11mec). For factorizing real symmetric positive definite sparse matrices, see nag_sparse_sym_chol_fac (f11jac). These functions should be used only when A is not banded and when the total number of nonzero elements is less than 10% of the total number of elements. In all other cases either the band functions or the general functions should be used.

3.3  Matrix Functions

nag_real_gen_matrix_exp (f01ecc) and nag_matop_complex_gen_matrix_exp (f01fcc) compute the matrix exponential, eA, of a real and complex square matrix A respectively. nag_real_symm_matrix_exp (f01edc) and nag_matop_complex_herm_matrix_exp (f01fdc) compute the matrix exponential of a real symmetric and complex Hermitian matrix respectively. If the matrix is real symmetric, or complex Hermitian then it is recommended that nag_real_symm_matrix_exp (f01edc), or nag_matop_complex_herm_matrix_exp (f01fdc) be used as they are more efficient and, in general, more accurate than nag_real_gen_matrix_exp (f01ecc) and nag_matop_complex_gen_matrix_exp (f01fcc).
nag_matop_real_gen_matrix_log (f01ejc) and nag_matop_complex_gen_matrix_log (f01fjc) compute the principal matrix logarithm, logA, of a real and complex square matrix A respectively.
nag_matop_real_gen_matrix_fun_std (f01ekc) and nag_matop_complex_gen_matrix_fun_std (f01fkc) compute the matrix exponential, sine, cosine, sinh or cosh of a real and complex square matrix A respectively. If the matrix exponential is required then it is recommended that nag_real_gen_matrix_exp (f01ecc) or nag_matop_complex_gen_matrix_exp (f01fcc) be used as they are, in general, more accurate than nag_matop_real_gen_matrix_fun_std (f01ekc) and nag_matop_complex_gen_matrix_fun_std (f01fkc).
nag_matop_real_gen_matrix_fun_usd (f01emc) computes the matrix function, fA, of a real square matrix. nag_matop_complex_gen_matrix_fun_usd (f01fmc) computes the matrix function of a complex square matrix. The derivatives of f are required for these computations. nag_matop_real_gen_matrix_fun_usd (f01emc) and nag_matop_complex_gen_matrix_fun_usd (f01fmc) use derivatives you have supplied. nag_matop_real_symm_matrix_fun (f01efc) and nag_matop_complex_herm_matrix_fun (f01ffc) compute the matrix function, fA, of a real symmetric and complex Hermitian matrix A respectively. If the matrix is real symmetric or complex Hermitian then it is recommended that nag_matop_real_symm_matrix_fun (f01efc) or nag_matop_complex_herm_matrix_fun (f01ffc) be used as they are more efficient and, in general, more accurate than nag_matop_real_gen_matrix_fun_usd (f01emc) and nag_matop_complex_gen_matrix_fun_usd (f01fmc).

4  Decision Trees

The decision trees show the functions in this chapter and in Chapter f04 that should be used for inverting matrices of various types.
(i) Matrix Inversion:

Tree 1

Is A an n by n matrix of rank n? _
yes		 Is A a real matrix? _
yes		 see Tree 2
| no
|
| see Tree 3
no
|
see Tree 4

Tree 2: Inverse of a real n by n matrix of full rank

Is A a band matrix? _
yes		 See Note 1.
no
|
Is A symmetric? _
yes		 Is A positive definite? _
yes		 Is one triangle of A stored as a linear array? _
yes		 f07gdc and f07gjc
| | no
|
| | f07fdc and f07fjc
| no
|
| Is one triangle of A stored as a linear array? _
yes		 f07pdc and f07pjc
| no
|
| f07mdc and f07mjc
no
|
Is A triangular? _
yes		 Is A stored as a linear array? _
yes		 f07ujc
| no
|
| f07tjc
no
|
f07adc and f07ajc

Tree 3: Inverse of a complex n by n matrix of full rank

Is A a band matrix? _
yes		 See Note 1.
no
|
Is A Hermitian? _
yes		 Is A positive definite? _
yes		 Is one triangle of A stored as a linear array? _
yes		 f07grc and f07gwc
| | no
|
| | f07frc and f07fwc
| no
|
| Is one triangle A stored as a linear array? _
yes		 f07prc and f07pwc
| no
|
| f07mrc and f07mwc
no
|
Is A symmetric? _
yes		 Is one triangle of A stored as a linear array? _
yes		 f07qrc and f07qwc
| no
|
| f07nrc and f07nwc
no
|
Is A triangular? _
yes		 Is A stored as a linear array? _
yes		 f07uwc
| no
|
| f07twc
no
|
f07anc or f07arc and f07awc

Tree 4: Pseudo-inverses

Is A a complex matrix? _
yes		 Is A of full rank? _
yes		 Is A an m by n matrix with m<n? _
yes		 f08avc and f08awc or f08axc
| | no
|
| | f08asc and f08auc or f08atc
| no
|
| f08kpc
no
|
Is A of full rank? _
yes		 Is A an m by n matrix with m<n? _
yes		 f08ahc and f08ajc or f08akc
| no
|
| f08aec and f08agc or f08afc
no
|
Is A an m by n matrix with m<n? _
yes		 f08kbc
no
|
f02wec
(ii) Matrix Factorizations: see the decision trees in Section 4 in the f02 and f04 Chapter Introductions.
(iii) Matrix Functions: not required.
Note 1: the inverse of a band matrix A does not in general have the same shape as A, and no functions are provided specifically for finding such an inverse. The matrix must either be treated as a full matrix, or the equations AX=B must be solved, where B has been initialized to the identity matrix I. In the latter case, see the decision trees in Section 4 in the f04 Chapter Introduction.
Note 2: by ‘guaranteed accuracy’ we mean that the accuracy of the inverse is improved by use of the iterative refinement technique using additional precision.

5  Functionality Index

Matrix Function 
    complex Hermitian n by n matrix 
        matrix exponential nag_matop_complex_herm_matrix_exp (f01fdc)
        matrix function nag_matop_complex_herm_matrix_fun (f01ffc)
    complex n by n matrix 
        matrix exponential nag_matop_complex_gen_matrix_exp (f01fcc)
    real n by n matrix 
        matrix exponential nag_real_gen_matrix_exp (f01ecc)
        matrix exponential, sine, cosine, sinh or cosh nag_matop_real_gen_matrix_fun_std (f01ekc)
        matrix function, using user-supplied derivatives nag_matop_real_gen_matrix_fun_usd (f01emc)
        matrix logarithm nag_matop_real_gen_matrix_log (f01ejc)
    real symmetric n by n matrix 
        matrix exponential nag_real_symm_matrix_exp (f01edc)
        matrix exponential, sine, cosine, sinh or cosh nag_matop_complex_gen_matrix_fun_std (f01fkc)
        matrix function nag_matop_real_symm_matrix_fun (f01efc)
        matrix function, using user-supplied derivatives nag_matop_complex_gen_matrix_fun_usd (f01fmc)
        matrix logarithm nag_matop_complex_gen_matrix_log (f01fjc)
Matrix Transformations, 
    real band symmetric positive definite matrix, 
        variable bandwidth, LDLT factorization nag_real_cholesky_skyline (f01mcc)

6  Functions Withdrawn or Scheduled for Withdrawal

Withdrawn
Function		Mark of
Withdrawal
Replacement Function(s)
nag_complex_cholesky (f01bnc)25nag_zpotrf (f07frc)
nag_real_qr (f01qcc)25nag_dgeqrf (f08aec)
nag_real_apply_q (f01qdc)25nag_dormqr (f08agc)
nag_real_form_q (f01qec)25nag_dorgqr (f08afc)
nag_complex_qr (f01rcc)25nag_zgeqrf (f08asc)
nag_complex_apply_q (f01rdc)25nag_zunmqr (f08auc)
nag_complex_form_q (f01rec)25nag_zungqr (f08atc)

7  References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Higham N J (2008) Functions of Matrices: Theory and Computation SIAM, Philadelphia, PA, USA
Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford
Wilkinson J H (1977) Some recent advances in numerical linear algebra The State of the Art in Numerical Analysis (ed D A H Jacobs) Academic Press
Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag
f01 Chapter Contents
f01 Chapter Introduction (PDF version)
NAG C Library Manual
© The Numerical Algorithms Group Ltd, Oxford, UK. 2012