NAG Library Chapter Introduction
f01 – Matrix Operations, Including Inversion
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) |
Nonsingular square matrices of order .
If , a square matrix of order , is nonsingular (has rank ), then its inverse exists and satisfies the equations (the identity or unit matrix).
It is worth noting that if , so that is the ‘residual’ matrix, then a bound on the relative error is given by , i.e.,
|
(ii) |
General real rectangular matrices.
A real matrix has no inverse if it is square ( by ) and singular (has rank ), or if it is of shape ( by ) with , but there is a Generalized or Pseudo-inverse
which satisfies the equations
(which of course are also satisfied by the inverse of if is square and nonsingular).
(a) |
if and then can be factorized using a factorization, given by
where is an by orthogonal matrix and is an by , nonsingular, upper triangular matrix. The pseudo-inverse of is then given by
where consists of the first columns of . |
(b) |
if and then can be factorized using an RQ factorization, given by
where is an by orthogonal matrix and is an by , nonsingular, upper triangular matrix. The pseudo-inverse of is then given by
where consists of the first columns of . |
(c) |
if and then can be factorized using a factorization, with column interchanges, as
where is an by orthogonal matrix, is an by upper trapezoidal matrix and is an by permutation matrix. The pseudo-inverse of is then given by
where consists of the first columns of . |
(d) |
if , then can be factorized as the singular value decomposition
where is an by orthogonal matrix, is an by orthogonal matrix and is an by diagonal matrix with non-negative diagonal elements . The first columns of and are the left- and right-hand singular vectors of respectively and the diagonal elements of are the singular values of . may be chosen so that
and in this case if then
If and consist of the first columns of and respectively and is an by diagonal matrix with diagonal elements then is given by
and the pseudo-inverse of is given by
Notice that
which is the classical eigenvalue (spectral) factorization of . |
(e) |
if 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 is
where and are unitary, the conjugate transpose of and 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 , the matrix function is a matrix with the same dimensions as which provides a generalization of the scalar function .
If
has a full set of eigenvectors
then
can be factorized as
where
is the diagonal matrix whose diagonal elements,
, are the eigenvalues of
.
is given by
where
is the diagonal matrix whose
th diagonal element is
.
In general,
may not have a full set of eigenvectors. The matrix function can then be defined via a Cauchy integral. For
,
where
is a closed contour surrounding the eigenvalues of
, and
is analytic within
.
Some matrix functions are defined implicitly. A matrix logarithm is a solution
to the equation
In general
is not unique, but if
has no eigenvalues on the closed negative real line then a unique
principal logarithm exists whose eigenvalues have imaginary part between
and
. Similarly, a matrix square root is a solution
to the equation
If has no eigenvalues on the closed negative real line then a unique principal square root exists with eigenvalues in the right half-plane. If has a vanishing eigenvalue then cannot be computed. If the vanishing eigenvalue is defective (its algebraic multiplicity exceeds its geometric multiplicity, or equivalently it occurs in a Jordan block of size greater than ) then the square root cannot be computed. If the vanishing eigenvalue is semisimple (its algebraic and geometric multiplicities are equal, or equivalently it occurs only in Jordan blocks of size ) then a square root can be computed.
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, cosh, square root and general real power 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.
The Fréchet derivative of a matrix function
in the direction of the matrix
is the linear function mapping
to
such that
The Fréchet derivative measures the first-order effect on
of perturbations in
.
Chapter f01 contains functions for calculating the Fréchet derivative of the exponential, logarithm and real powers of both real and complex matrices.
The condition number of a matrix function is a measure of its sensitivity to perturbations in the data. The absolute condition number measures these perturbations in an absolute sense, and is defined by
The relative condition number, which is usually of more interest, measures these perturbations in a relative sense, and is defined by
The absolute and relative condition numbers can be expressed in terms of the norm of the Fréchet derivative by
Chapter f01 contains routines for calculating the condition number of the matrix exponential, logarithm, sine, cosine, sinh, cosh, square root and general real power of both real and complex matrices. It also contains routines for estimating the condition number of a general function of a real or complex matrix.
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
can be written as
, the solution should
never be computed by first inverting
and then computing
; 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) |
Nonsingular square matrices of order This chapter describes techniques for inverting a general real matrix 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 , for , where is the th column of the identity matrix.
The residual matrix , where is a computed inverse of , conveys useful information. Firstly
is a bound on the relative error in and secondly guarantees the convergence of the iterative process in the ‘corrected’ inverse functions.
The decision trees for inversion show which functions in
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 factorization of the matrix and nag_dgeqp3 (f08bfc) returns the 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 has rank then the smallest singular values will be negligible and the pseudo-inverse of can be obtained as as described in Section 2. If the rank of 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
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,
, of a real and complex square matrix
respectively. If estimates of the condition number of the matrix exponential are required then
nag_matop_real_gen_matrix_cond_exp (f01jgc) and
nag_matop_complex_gen_matrix_cond_exp (f01kgc) should be used. If Fréchet derivatives are required then
nag_matop_real_gen_matrix_frcht_exp (f01jhc) and
nag_matop_complex_gen_matrix_frcht_exp (f01khc) should be used.
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,
, of a real and complex square matrix
respectively. If estimates of the condition number of the matrix logarithm are required then
nag_matop_real_gen_matrix_cond_log (f01jjc) and
nag_matop_complex_gen_matrix_cond_log (f01kjc) should be used. If Fréchet derivatives are required then
nag_matop_real_gen_matrix_frcht_log (f01jkc) and
nag_matop_complex_gen_matrix_frcht_log (f01kkc) should be used.
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
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). If estimates of the condition number of the matrix function are required then
nag_matop_real_gen_matrix_cond_std (f01jac) and
nag_matop_complex_gen_matrix_cond_std (f01kac) should be used.
nag_matop_real_gen_matrix_fun_num (f01elc) and
nag_matop_real_gen_matrix_fun_usd (f01emc) compute the matrix function,
, of a real square matrix.
nag_matop_complex_gen_matrix_fun_num (f01flc) and
nag_matop_complex_gen_matrix_fun_usd (f01fmc) compute the matrix function of a complex square matrix. The derivatives of
are required for these computations.
nag_matop_real_gen_matrix_fun_num (f01elc) and
nag_matop_complex_gen_matrix_fun_num (f01flc) use numerical differentiation to obtain the derivatives of
.
nag_matop_real_gen_matrix_fun_usd (f01emc) and
nag_matop_complex_gen_matrix_fun_usd (f01fmc) use derivatives you have supplied. If estimates of the condition number are required but you are not supplying derivatives then
nag_matop_real_gen_matrix_cond_num (f01jbc) and
nag_matop_complex_gen_matrix_cond_num (f01kbc) should be used.
If estimates of the condition number of the matrix function are required and you are supplying derivatives of
, then
nag_matop_real_gen_matrix_cond_usd (f01jcc) and
nag_matop_complex_gen_matrix_cond_usd (f01kcc) should be used.
nag_matop_real_gen_matrix_actexp (f01gac) and
nag_matop_complex_gen_matrix_actexp (f01hac) compute the matrix function
for explicitly stored dense real and complex matrices
and
respectively while
nag_matop_real_gen_matrix_actexp_rcomm (f01gbc) and
nag_matop_complex_gen_matrix_actexp_rcomm (f01hbc) compute the same using reverse communication. In the latter case, control is returned to you. You should calculate any required matrix-matrix products and then call the function again.
nag_matop_real_gen_matrix_sqrt (f01enc) and
nag_matop_complex_gen_matrix_sqrt (f01fnc) compute the principal square root
of a real and complex square matrix
respectively. If
is complex and upper triangular then
nag_matop_complex_tri_matrix_sqrt (f01fpc) should be used. If
is real and upper quasi-triangular then
nag_matop_real_tri_matrix_sqrt (f01epc) should be used. If estimates of the condition number of the matrix square root are required then
nag_matop_real_gen_matrix_cond_sqrt (f01jdc) and
nag_matop_complex_gen_matrix_cond_sqrt (f01kdc) should be used.
4 Decision Trees
The decision trees show the functions in this chapter and in
Chapter f07 and
Chapter f08 that should be used for inverting matrices of various types. They also show which function should be used to calculate various matrix functions.
(i) Matrix Inversion:
Tree 1
Is an by matrix of rank ? |
_ yes |
Is 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 band matrix? |
_ yes |
See Note 1. |
no | |
|
Is symmetric? |
_ yes |
Is positive definite? |
_ yes |
Is one triangle of stored as a linear array? |
_ yes |
f07gdc and f07gjc |
| |
|
| |
|
no | |
|
| |
| | |
|
f07fdc and f07fjc |
| |
|
no | |
|
| |
|
Is one triangle of stored as a linear array? |
_ yes |
f07pdc and f07pjc |
| |
|
no | |
|
| |
|
f07mdc and f07mjc |
no | |
|
Is triangular? |
_ yes |
Is 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
Tree 4: Pseudo-inverses
Note 1: the inverse of a band matrix
does not in general have the same shape as
, and no functions are provided specifically for finding such an inverse. The matrix must either be treated as a full matrix, or the equations
must be solved, where
has been initialized to the identity matrix
. 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.
(ii)
Matrix Factorizations: see the decision trees in Section 4 in the
f02 and
f04 Chapter Introductions.
(iii)
Matrix Functions:
Tree 5: Matrix functions of an n by n real matrix
Is required? |
_ yes |
Is stored in dense format? |
_ yes |
f01gac |
| |
|
no | |
|
| |
|
f01gbc |
no | |
|
Is real symmetric? |
_ yes |
Is required? |
_ yes |
f01edc |
| |
|
no | |
|
| |
|
f01efc |
no | |
|
Is or
or
or
required? |
_ yes |
Is the condition number of the matrix function required? |
_ yes |
f01jac |
| |
|
no | |
|
| |
|
f01ekc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix logarithm required? |
_ yes |
f01jjc |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix logarithm required? |
_ yes |
f01jkc |
| |
|
no | |
|
| |
|
f01ejc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix exponential required? |
_ yes |
f01jgc |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix exponential required? |
_ yes |
f01jhc |
| |
|
no | |
|
| |
|
f01ecc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix square root required? |
_ yes |
f01jdc |
| |
|
no | |
|
| |
|
Is the matrix upper quasi-triangular? |
_ yes |
f01epc |
| |
|
no | |
|
| |
|
f01enc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix power required? |
_ yes |
f01jec |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix power required? |
_ yes |
f01jfc |
| |
|
no | |
|
| |
|
f01eqc |
no | |
|
will be computed. Will derivatives of be supplied by the user? |
_ yes |
Is the condition number of the matrix function required? |
_ yes |
f01jcc |
| |
|
no | |
|
| |
|
f01emc |
no | |
|
Is the condition number of the matrix function required? |
_ yes |
f01jbc |
no | |
|
f01elc |
Tree 6: Matrix functions of an n by n complex matrix
Is required? |
_ yes |
Is stored in dense format? |
_ yes |
f01hac |
| |
|
no | |
|
| |
|
f01hbc |
no | |
|
Is complex Hermitian? |
_ yes |
Is required? |
_ yes |
f01fdc |
| |
|
no | |
|
| |
|
f01ffc |
no | |
|
Is or
or
or
required? |
_ yes |
Is the condition number of the matrix function required? |
_ yes |
f01kac |
| |
|
no | |
|
| |
|
f01fkc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix logarithm required? |
_ yes |
f01kjc |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix logarithm required? |
_ yes |
f01kkc |
| |
|
no | |
|
| |
|
f01fjc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix exponential required? |
_ yes |
f01kgc |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix exponential required? |
_ yes |
f01khc |
| |
|
no | |
|
| |
|
f01fcc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix square root required? |
_ yes |
f01kdc |
| |
|
no | |
|
| |
|
Is the matrix upper triangular? |
_ yes |
f01fpc |
| |
|
no | |
|
| |
|
f01fnc |
no | |
|
Is required? |
_ yes |
Is the condition number of the matrix power required? |
_ yes |
f01kec |
| |
|
no | |
|
| |
|
Is the Fréchet derivative of the matrix power required? |
_ yes |
f01kfc |
| |
|
no | |
|
| |
|
f01fqc |
no | |
|
will be computed. Will derivatives of be supplied by the user? |
_ yes |
Is the condition number of the matrix function required? |
_ yes |
f01kcc |
| |
|
no | |
|
| |
|
f01fmc |
no | |
|
Is the condition number of the matrix function required? |
_ yes |
f01kbc |
no | |
|
f01flc |
5 Functionality Index
complex Hermitian n by n matrix, | | |
real symmetric n by n matrix, | | |
real band symmetric positive definite matrix, | | |
6 Auxiliary Functions Associated with Library Function Arguments
None.
7 Functions Withdrawn or Scheduled for Withdrawal
The following lists all those functions that have been withdrawn since Mark 23 of the Library or are scheduled for withdrawal at one of the next two marks.
8 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