f01 Chapter Introduction : NAG Library, Mark 25

This chapter provides facilities for four types of problem:

(i)	Matrix Inversion
(ii)	Matrix Factorizations
(iii)	Matrix Arithmetic and Manipulation
(iv)	Matrix Functions

These problems are discussed separately in Section 2.1, Section 2.2, Section 2.3 and Section 2.4.

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.)

The intention of functions in this section (sub-chapter f01v) is to cater for some of the commonly occurring operations in matrix manipulation, i.e., conversion between different storage formats,such as conversion between rectangular band matrix storage and packed band matrix storage. For vector or matrix-vector or matrix-matrix operations refer to Chapters f06 and f16.

Given a square matrix

A

, the matrix function

f (A)

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 V^{- 1},

where

D

is the diagonal matrix whose diagonal elements,

d_{i}

, are the eigenvalues of

A

.

f (A)

is given by

f (A) = V f (D) V^{- 1},

where

f (D)

is the diagonal matrix whose

i

th diagonal element is

f (d_{i})

.

In general,

A

may not have a full set of eigenvectors. The matrix function can then be defined via a Cauchy integral. For

A \in ℂ^{n \times n}

,

f (A) = \frac{1}{2 π i} \int_{Γ} f (z) {(z I - A)}^{- 1} d z,

where

Γ

is a closed contour surrounding the eigenvalues of

A

, and

f

is analytic within

Γ

.

Some matrix functions are defined implicitly. A matrix logarithm is a solution

X

to the equation

e^{X} = A .

In general

X

is not unique, but if

A

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

X

to the equation

X^{2} = A .

If

A

has no eigenvalues on the closed negative real line then a unique principal square root exists with eigenvalues in the right half-plane. If

A

has a vanishing eigenvalue then

\log (A)

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

1

) 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

1

) 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

f (A)

in the direction of the matrix

E

is the linear function mapping

E

to

L_{f} (A, E)

such that

f (A + E) - f (A) - L_{f} (A, E) = O (‖E‖) .

The Fréchet derivative measures the first-order effect on

f (A)

of perturbations in

A

. 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

{cond}_{abs} (f, A) : = \lim_{ε \to 0} \sup_{\{‖E‖ \to 0\}} \frac{‖f (A + E) - f (A)‖}{ε} .

The relative condition number, which is usually of more interest, measures these perturbations in a relative sense, and is defined by

{cond}_{rel} (f, A) = {cond}_{abs} (f, A) \frac{‖A‖}{‖f (A)‖} .

The absolute and relative condition numbers can be expressed in terms of the norm of the Fréchet derivative by

{cond}_{abs} (f, A) = \max_{E \neq 0} \frac{‖L (A, E)‖}{‖E‖},

{cond}_{rel} (f, A) = \frac{‖A‖}{‖f (A)‖} \max_{E \neq 0} \frac{‖L (A, E)‖}{‖E‖} .

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.

Note: before using any function for matrix inversion, consider carefully whether it is really needed.

Although the solution of a set of linear equations

A x = b

can be written as

x = A^{- 1} b

, the solution should never be computed by first inverting

A

and then computing

A^{- 1} b

; 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

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

A x_{i} = e_{i}

, for

i = 1, 2, \dots, n

, where

e_{i}

is the

i

th column of the identity matrix.

The residual matrix

R = A X - 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‖ < \frac{1}{2}

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

Q R

factorization of the matrix and nag_dgeqp3 (f08bfc) returns the

Q R

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 \leq k = \min (m, n)

then the

k - r

smallest singular values will be negligible and the pseudo-inverse of

A

can be obtained as

A^{+} = V Σ^{- 1} U^{T}

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.

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.

The functions in the f01v (LAPACK) section are designed to allow conversion between full storage format and one of the packed storage schemes required by some of the functions in Chapters f02, f04, f06, f07 and f08.

Functions with NAG short name beginning f01v may be called either by their NAG short names or by their NAG long names. The NAG long names for a function is simply the LAPACK name (in lower case) prepended by nag_, for example, nag_dtrttf is the long name for f01vec.

When using the NAG Library, the double precision form of the LAPACK name must be used (beginning with d- or z-).

References to Chapter f01 functions in the manual normally include the LAPACK double precision names, for example, nag_dtrttf (f01vec).

The LAPACK function names follow a simple scheme (which is similar to that used for the BLAS in Chapter f16). Most names have the structure nag_xyytzz, where the components have the following meanings:

– the initial letter, x, indicates the data type (real or complex) and precision:

d – real, double precision
z – complex, double precision

– the fourth letter, t, indicates that the function is performing a storage scheme transformation (conversion)

– the letters yy indicate the original storage scheme used to store a triangular part of the matrix

A

, while the letters zz indicate the target storage scheme of the conversion (yy cannot equal zz since this would do nothing):

tf – Rectangular Full Packed Format (RFP)
tp – Packed Format
tr – Full Format

nag_real_gen_matrix_exp (f01ecc) and nag_matop_complex_gen_matrix_exp (f01fcc) compute the matrix exponential,

e^{A}

, of a real and complex square matrix

A

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,

e^{A}

, 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,

\log (A)

, of a real and complex square matrix

A

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

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). 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,

f (A)

, 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

f

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

f

. 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

f

, then nag_matop_real_gen_matrix_cond_usd (f01jcc) and nag_matop_complex_gen_matrix_cond_usd (f01kcc) should be used.

If the matrix

A

is real symmetric or complex Hermitian then it is recommended that to compute the matrix function,

f (A)

, nag_matop_real_symm_matrix_fun (f01efc) and nag_matop_complex_herm_matrix_fun (f01ffc) are used respectively as they are more efficient and, in general, more accurate than nag_matop_real_gen_matrix_fun_num (f01elc), nag_matop_real_gen_matrix_fun_usd (f01emc), nag_matop_complex_gen_matrix_fun_num (f01flc) and nag_matop_complex_gen_matrix_fun_usd (f01fmc).

nag_matop_real_gen_matrix_actexp (f01gac) and nag_matop_complex_gen_matrix_actexp (f01hac) compute the matrix function

e^{t A} B

for explicitly stored dense real and complex matrices

A

and

B

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. See Section 3.2.2 in the Essential Introduction for further information.

nag_matop_real_gen_matrix_sqrt (f01enc) and nag_matop_complex_gen_matrix_sqrt (f01fnc) compute the principal square root

A^{1 / 2}

of a real and complex square matrix

A

respectively. If

A

is complex and upper triangular then nag_matop_complex_tri_matrix_sqrt (f01fpc) should be used. If

A

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.

nag_matop_real_gen_matrix_pow (f01eqc) and nag_matop_complex_gen_matrix_pow (f01fqc) compute the matrix power

A^{p}

, where

p \in ℝ

, of real and complex matrices respectively. If estimates of the condition number of the matrix power are required then nag_matop_real_gen_matrix_cond_pow (f01jec) and nag_matop_complex_gen_matrix_cond_pow (f01kec) should be used. If Fréchet derivatives are required then nag_matop_real_gen_matrix_frcht_pow (f01jfc) and nag_matop_complex_gen_matrix_frcht_pow (f01kfc) should be used.

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:

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

A X = 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.

(ii) Matrix Factorizations: see the decision trees in Section 4 in the f02 and f04 Chapter Introductions.

(iii) Matrix Arithmetic and Manipulation: not appropriate.

(iv) Matrix Functions:

None.

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.

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

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


see Tree 4

Is $A$ a band matrix?			See Note 1.
		yes
	no
Is $A$ symmetric?			Is $A$ positive definite?			Is one triangle of $A$ stored as a linear array?			f07gdc and f07gjc
		yes			yes			yes
	no			no			no
						f07fdc and f07fjc


			Is one triangle of $A$ stored as a linear array?			f07pdc and f07pjc
					yes
				no
			f07mdc and f07mjc


Is $A$ triangular?			Is $A$ stored as a linear array?			f07ujc
		yes			yes
	no			no
			f07tjc


f07adc and f07ajc

Is $A$ a band matrix?			See Note 1.
		yes
	no
Is $A$ Hermitian?			Is $A$ positive definite?			Is one triangle of $A$ stored as a linear array?			f07grc and f07gwc
		yes			yes			yes
	no			no			no
						f07frc and f07fwc


			Is one triangle $A$ stored as a linear array?			f07prc and f07pwc
					yes
				no
			f07mrc and f07mwc


Is $A$ symmetric?			Is one triangle of $A$ stored as a linear array?			f07qrc and f07qwc
		yes			yes
	no			no
			f07nrc and f07nwc


Is $A$ triangular?			Is $A$ stored as a linear array?			f07uwc
		yes			yes
	no			no
			f07twc


f07anc or f07arc and f07awc

Is $A$ a complex matrix?			Is $A$ of full rank?			Is $A$ an $m$ by $n$ matrix with $m < n$ ?			f08avc and f08awc or f08axc
		yes			yes			yes
	no			no			no
						f08asc and f08auc or f08atc


			f08kpc


Is $A$ of full rank?			Is $A$ an $m$ by $n$ matrix with $m < n$ ?			f08ahc and f08ajc or f08akc
		yes			yes
	no			no
			f08aec and f08agc or f08afc


f08kbc

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

NAG Library Chapter Introduction

f01 – Matrix Operations, Including Inversion

▸▿ Contents

1 Scope of the Chapter

2 Background to the Problems

2.1 Matrix Inversion

2.2 Matrix Factorizations

2.3 Matrix Arithmetic and Manipulation

2.4 Matrix Functions

3 Recommendations on Choice and Use of Available Functions

3.1 Matrix Inversion

3.2 Matrix Factorizations

3.3 Matrix Arithmetic and Manipulation

3.3.1 NAG Names and LAPACK Names

3.4 Matrix Functions

4 Decision Trees

Tree 1

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

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

Tree 4: Pseudo-inverses

Tree 5: Matrix functions $f (A)$ of an n by n real matrix $A$

Tree 6: Matrix functions $f (A)$ of an n by n complex matrix $A$

5 Functionality Index

6 Auxiliary Functions Associated with Library Function Arguments

7 Functions Withdrawn or Scheduled for Withdrawal

8 References

Is $e^{t A} B$ required?			Is $A$ stored in dense format?			f01gac
		yes			yes
	no			no
			f01gbc


Is $A$ real symmetric?			Is $e^{A}$ required?			f01edc
		yes			yes
	no			no
			f01efc


Is $\cos (A)$ or $\cosh (A)$ or $\sin (A)$ or $\sinh (A)$ required?			Is the condition number of the matrix function required?			f01jac
		yes			yes
	no			no
			f01ekc


Is $\log (A)$ required?			Is the condition number of the matrix logarithm required?			f01jjc
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix logarithm required?			f01jkc
					yes
				no
			f01ejc


Is $\exp (A)$ required?			Is the condition number of the matrix exponential required?			f01jgc
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix exponential required?			f01jhc
					yes
				no
			f01ecc


Is $A^{1 / 2}$ required?			Is the condition number of the matrix square root required?			f01jdc
		yes			yes
	no			no
			Is the matrix upper quasi-triangular?			f01epc
					yes
				no
			f01enc


Is $A^{p}$ required?			Is the condition number of the matrix power required?			f01jec
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix power required?			f01jfc
					yes
				no
			f01eqc


$f (A)$ will be computed. Will derivatives of $f$ be supplied by the user?			Is the condition number of the matrix function required?			f01jcc
		yes			yes
	no			no
			f01emc


Is the condition number of the matrix function required?			f01jbc
		yes
	no
f01elc

Is $e^{t A} B$ required?			Is $A$ stored in dense format?			f01hac
		yes			yes
	no			no
			f01hbc


Is $A$ complex Hermitian?			Is $e^{A}$ required?			f01fdc
		yes			yes
	no			no
			f01ffc


Is $\cos (A)$ or $\cosh (A)$ or $\sin (A)$ or $\sinh (A)$ required?			Is the condition number of the matrix function required?			f01kac
		yes			yes
	no			no
			f01fkc


Is $\log (A)$ required?			Is the condition number of the matrix logarithm required?			f01kjc
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix logarithm required?			f01kkc
					yes
				no
			f01fjc


Is $\exp (A)$ required?			Is the condition number of the matrix exponential required?			f01kgc
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix exponential required?			f01khc
					yes
				no
			f01fcc


Is $A^{1 / 2}$ required?			Is the condition number of the matrix square root required?			f01kdc
		yes			yes
	no			no
			Is the matrix upper triangular?			f01fpc
					yes
				no
			f01fnc


Is $A^{p}$ required?			Is the condition number of the matrix power required?			f01kec
		yes			yes
	no			no
			Is the Fréchet derivative of the matrix power required?			f01kfc
					yes
				no
			f01fqc


$f (A)$ will be computed. Will derivatives of $f$ be supplied by the user?			Is the condition number of the matrix function required?			f01kcc
		yes			yes
	no			no
			f01fmc


Is the condition number of the matrix function required?			f01kbc
		yes
	no
f01flc

NAG Library Chapter Introductionf01 – Matrix Operations, Including Inversion

▸▿ Contents

1 Scope of the Chapter

2 Background to the Problems

2.1 Matrix Inversion

2.2 Matrix Factorizations

2.3 Matrix Arithmetic and Manipulation

2.4 Matrix Functions

3 Recommendations on Choice and Use of Available Functions

3.1 Matrix Inversion

3.2 Matrix Factorizations

3.3 Matrix Arithmetic and Manipulation

3.3.1 NAG Names and LAPACK Names

3.4 Matrix Functions

4 Decision Trees

Tree 1

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

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

Tree 4: Pseudo-inverses

Tree 5: Matrix functions fA of an n by n real matrix A

Tree 6: Matrix functions fA of an n by n complex matrix A

5 Functionality Index

6 Auxiliary Functions Associated with Library Function Arguments

7 Functions Withdrawn or Scheduled for Withdrawal

8 References

NAG Library Chapter Introduction

f01 – Matrix Operations, Including Inversion

Tree 5: Matrix functions $f (A)$ of an n by n real matrix $A$

Tree 6: Matrix functions $f (A)$ of an n by n complex matrix $A$