F08 Chapter Introduction : NAG Library, Mark 24

This chapter provides routines for the solution of linear least squares problems, eigenvalue problems and singular value problems, as well as associated computations. It provides routines for:

– solution of linear least squares problems
– solution of symmetric eigenvalue problems
– solution of nonsymmetric eigenvalue problems
– solution of singular value problems
– solution of generalized symmetric-definite eigenvalue problems
– solution of generalized nonsymmetric eigenvalue problems
– solution of generalized singular value problems
– solution of generalized linear least squares problems
– matrix factorizations associated with the above problems
– estimating condition numbers of eigenvalue and eigenvector problems
– estimating the numerical rank of a matrix
– solution of the Sylvester matrix equation

Routines are provided for both real and complex data.

For a general introduction to the solution of linear least squares problems, you should turn first to Chapter F04. The decision trees, at the end of Chapter F04, direct you to the most appropriate routines in Chapters F04 or F08. Chapters F04 and F08 contain Black Box (or driver) routines which enable standard linear least squares problems to be solved by a call to a single routine.

For a general introduction to eigenvalue and singular value problems, you should turn first to Chapter F02. The decision trees, at the end of Chapter F02, direct you to the most appropriate routines in Chapters F02 or F08. Chapters F02 and F08 contain Black Box (or driver) routines which enable standard types of problem to be solved by a call to a single routine. Often routines in Chapter F02 call Chapter F08 routines to perform the necessary computational tasks.

The routines in this chapter (Chapter F08) handle only dense, band, tridiagonal and Hessenberg matrices (not matrices with more specialised structures, or general sparse matrices). The tables in Section 3 and the decision trees in Section 4 direct you to the most appropriate routines in Chapter F08.

The routines in this chapter have all been derived from the LAPACK project (see Anderson et al. (1999)). They have been designed to be efficient on a wide range of high-performance computers, without compromising efficiency on conventional serial machines.

It is not expected that you will need to read all of the following sections, but rather you will pick out those sections relevant to your particular problem.

This section is only a brief introduction to the numerical solution of linear least squares problems, eigenvalue and singular value problems. Consult a standard textbook for a more thorough discussion, for example Golub and Van Loan (1996).

The linear least squares problem is

\underset{x}{minimize} {‖b - A x‖}_{2},

(1)

where

A

is an

m

by

n

matrix,

b

is a given

m

element vector and

x

is an

n

-element solution vector.

In the most usual case

m \geq n

and

rank (A) = n

, so that

A

has full rank and in this case the solution to problem (1) is unique; the problem is also referred to as finding a least squares solution to an overdetermined system of linear equations.

When

m < n

and

rank (A) = m

, there are an infinite number of solutions

x

which exactly satisfy

b - A x = 0

. In this case it is often useful to find the unique solution

x

which minimizes

{‖x‖}_{2}

, and the problem is referred to as finding a minimum norm solution to an underdetermined system of linear equations.

In the general case when we may have

rank (A) < \min (m, n)

– in other words,

A

may be rank-deficient – we seek the minimum norm least squares solution

x

which minimizes both

{‖x‖}_{2}

and

{‖b - A x‖}_{2}

.

This chapter (Chapter F08) contains driver routines to solve these problems with a single call, as well as computational routines that can be combined with routines in Chapter F07 to solve these linear least squares problems. Chapter F04 also contains Black Box routines to solve these linear least squares problems in standard cases. The next two sections discuss the factorizations that can be used in the solution of linear least squares problems.

A number of routines are provided for factorizing a general rectangular

m

by

n

matrix

A

, as the product of an orthogonal matrix (unitary if complex) and a triangular (or possibly trapezoidal) matrix.

A real matrix

Q

is orthogonal if

Q^{T} Q = I

; a complex matrix

Q

is unitary if

Q^{H} Q = I

. Orthogonal or unitary matrices have the important property that they leave the

2

-norm of a vector invariant, so that

{‖x‖}_{2} = {‖Q x‖}_{2},

if

Q

is orthogonal or unitary. They usually help to maintain numerical stability because they do not amplify rounding errors.

Orthogonal factorizations are used in the solution of linear least squares problems. They may also be used to perform preliminary steps in the solution of eigenvalue or singular value problems, and are useful tools in the solution of a number of other problems.

The most common, and best known, of the factorizations is the

Q R

factorization given by

A = Q (\begin{matrix} R \\ 0 \end{matrix}),   if ​ m \geq n,

where

R

is an

n

by

n

upper triangular matrix and

Q

is an

m

by

m

orthogonal (or unitary) matrix. If

A

is of full rank

n

, then

R

is nonsingular. It is sometimes convenient to write the factorization as

A = (Q_{1} Q_{2}) (\begin{matrix} R \\ 0 \end{matrix})

which reduces to

A = Q_{1} R,

where

Q_{1}

consists of the first

n

columns of

Q

, and

Q_{2}

the remaining

m - n

columns.

If

m < n

,

R

is trapezoidal, and the factorization can be written

A = Q (R_{1} R_{2}),   if ​ m < n,

where

R_{1}

is upper triangular and

R_{2}

is rectangular.

The

Q R

factorization can be used to solve the linear least squares problem (1) when

m \geq n

and

A

is of full rank, since

{‖b - A x‖}_{2} = {‖Q^{T} b - Q^{T} A x‖}_{2} = {‖(\begin{matrix} c_{1} - R x \\ c_{2} \end{matrix})‖}_{2},

where

c \equiv (\begin{matrix} c_{1} \\ c_{2} \end{matrix}) = (\begin{matrix} Q_{1}^{T} b \\ Q_{2}^{T} b \end{matrix}) = Q^{T} b;

and

c_{1}

is an

n

-element vector. Then

x

is the solution of the upper triangular system

R x = c_{1} .

The residual vector

r

is given by

r = b - A x = Q (\begin{matrix} 0 \\ c_{2} \end{matrix}) .

The residual sum of squares

{‖r‖}_{2}^{2}

may be computed without forming

r

explicitly, since

{‖r‖}_{2} = {‖b - A x‖}_{2} = {‖c_{2}‖}_{2} .

The

L Q

factorization is given by

A = (L 0) Q = (L 0) (\begin{matrix} Q_{1} \\ Q_{2} \end{matrix}) = L Q_{1},   if ​ m \leq n,

where

L

is

m

by

m

lower triangular,

Q

is

n

by

n

orthogonal (or unitary),

Q_{1}

consists of the first

m

rows of

Q

, and

Q_{2}

the remaining

n - m

rows.

The

L Q

factorization of

A

is essentially the same as the

Q R

factorization of

A^{T}

(

A^{H}

if

A

is complex), since

A = (L 0) Q \Leftrightarrow A^{T} = Q^{T} (\begin{matrix} L^{T} \\ 0 \end{matrix}) .

The

L Q

factorization may be used to find a minimum norm solution of an underdetermined system of linear equations

A x = b

where

A

is

m

by

n

with

m < n

and has rank

m

. The solution is given by

x = Q^{T} (\begin{matrix} L^{- 1} b \\ 0 \end{matrix}) .

To solve a linear least squares problem (1) when

A

is not of full rank, or the rank of

A

is in doubt, we can perform either a

Q R

factorization with column pivoting or a singular value decomposition.

The

Q R

factorization with column pivoting is given by

A = Q (\begin{matrix} R \\ 0 \end{matrix}) P^{T}, m \geq n,

where

Q

and

R

are as before and

P

is a (real) permutation matrix, chosen (in general) so that

|r_{11}| \geq |r_{22}| \geq \dots \geq |r_{n n}|

and moreover, for each

k

,

|r_{k k}| \geq {‖R_{k : j, j}‖}_{2}, j = k + 1, \dots, n .

If we put

R = (\begin{matrix} R_{11} & R_{12} \\ 0 & R_{22} \end{matrix})

where

R_{11}

is the leading

k

by

k

upper triangular sub-matrix of

R

then, in exact arithmetic, if

rank (A) = k

, the whole of the sub-matrix

R_{22}

in rows and columns

k + 1

to

n

would be zero. In numerical computation, the aim must be to determine an index

k

, such that the leading sub-matrix

R_{11}

is well-conditioned, and

R_{22}

is negligible, so that

R = (\begin{matrix} R_{11} & R_{12} \\ 0 & R_{22} \end{matrix}) ≃ (\begin{matrix} R_{11} & R_{12} \\ 0 & 0 \end{matrix}) .

Then

k

is the effective rank of

A

. See Golub and Van Loan (1996) for a further discussion of numerical rank determination.

The so-called basic solution to the linear least squares problem (1) can be obtained from this factorization as

x = P (\begin{matrix} R_{11}^{- 1} {\hat{c}}_{1} \\ 0 \end{matrix}),

where

{\hat{c}}_{1}

consists of just the first

k

elements of

c = Q^{T} b

.

The

Q R

factorization with column pivoting does not enable us to compute a minimum norm solution to a rank-deficient linear least squares problem, unless

R_{12} = 0

. However, by applying for further orthogonal (or unitary) transformations from the right to the upper trapezoidal matrix

(\begin{matrix} R_{11} & R_{12} \end{matrix})

,

R_{12}

can be eliminated:

(\begin{matrix} R_{11} & R_{12} \end{matrix}) Z = (\begin{matrix} T_{11} & 0 \end{matrix}) .

This gives the complete orthogonal factorization

AP = Q (\begin{matrix} T_{11} & 0 \\ 0 & 0 \end{matrix}) Z^{T}

from which the minimum norm solution can be obtained as

x = P Z (\begin{matrix} T_{11}^{- 1} & {\hat{c}}_{1} \\ 0 \end{matrix}) .

The

Q L

and

R Q

factorizations are given by

A = Q (\begin{matrix} 0 \\ L \end{matrix}),  if ​ m \geq n,

and

A = (\begin{matrix} 0 & R \end{matrix}) Q,  if ​ m \leq n .

The factorizations are less commonly used than either the

Q R

or

L Q

factorizations described above, but have applications in, for example, the computation of generalized

Q R

factorizations.

The singular value decomposition (SVD) of an

m

by

n

matrix

A

is given by

A = U Σ V^{T}, (A = U Σ V^{H} in the complex case)

where

U

and

V

are orthogonal (unitary) and

Σ

is an

m

by

n

diagonal matrix with real diagonal elements,

σ_{i}

, such that

σ_{1} \geq σ_{2} \geq \dots \geq σ_{\min (m, n)} \geq 0 .

The

σ_{i}

are the singular values of

A

and the first

\min (m, n)

columns of

U

and

V

are the left and right singular vectors of

A

. The singular values and singular vectors satisfy

A v_{i} = σ_{i} u_{i} and A^{T} u_{i} = σ_{i} v_{i} (or ​ A^{H} u_{i} = σ_{i} v_{i})

where

u_{i}

and

v_{i}

are the

i

th columns of

U

and

V

respectively.

The computation proceeds in the following stages.

The matrix $A$ is reduced to bidiagonal form $A = U_{1} B V_{1}^{T}$ if $A$ is real ( $A = U_{1} B V_{1}^{H}$ if $A$ is complex), where $U_{1}$ and $V_{1}$ are orthogonal (unitary if $A$ is complex), and $B$ is real and upper bidiagonal when $m \geq n$ and lower bidiagonal when $m < n$ , so that $B$ is nonzero only on the main diagonal and either on the first superdiagonal (if $m \geq n$ ) or the first subdiagonal (if $m < n$ ).
The SVD of the bidiagonal matrix $B$ is computed as $B = U_{2} Σ V_{2}^{T}$ , where $U_{2}$ and $V_{2}$ are orthogonal and $Σ$ is diagonal as described above. The singular vectors of $A$ are then $U = U_{1} U_{2}$ and $V = V_{1} V_{2}$ .

If

m ≫ n

, it may be more efficient to first perform a

Q R

factorization of

A

, and then compute the SVD of the

n

by

n

matrix

R

, since if

A = Q R

and

R = U Σ V^{T}

, then the SVD of

A

is given by

A = (Q U) Σ V^{T}

.

Similarly, if

m ≪ n

, it may be more efficient to first perform an

L Q

factorization of

A

.

This chapter supports two primary algorithms for computing the SVD of a bidiagonal matrix. They are:

(i)	the divide and conquer algorithm;
(ii)	the $Q R$ algorithm.

The divide and conquer algorithm is much faster than the

Q R

algorithm if singular vectors of large matrices are required.

The SVD may be used to find a minimum norm solution to a (possibly) rank-deficient linear least squares problem (1). The effective rank,

k

, of

A

can be determined as the number of singular values which exceed a suitable threshold. Let

\hat{Σ}

be the leading

k

by

k

sub-matrix of

Σ

, and

\hat{V}

be the matrix consisting of the first

k

columns of

V

. Then the solution is given by

x = \hat{V} {\hat{Σ}}^{- 1} {\hat{c}}_{1},

where

{\hat{c}}_{1}

consists of the first

k

elements of

c = U^{T} b = U_{2}^{T} U_{1}^{T} b

.

The simple type of linear least squares problem described in Section 2.1 can be generalized in various ways.

Linear least squares problems with equality constraints:
$find x to minimize S = {‖c - A x‖}_{2}^{2} subject to B x = d,$
where $A$ is $m$ by $n$ and $B$ is $p$ by $n$ , with $p \leq n \leq m + p$ . The equations $B x = d$ may be regarded as a set of equality constraints on the problem of minimizing $S$ . Alternatively the problem may be regarded as solving an overdetermined system of equations
$(\begin{matrix} A \\ B \end{matrix}) x = (\begin{matrix} c \\ d \end{matrix}),$
where some of the equations (those involving $B$ ) are to be solved exactly, and the others (those involving $A$ ) are to be solved in a least squares sense. The problem has a unique solution on the assumptions that $B$ has full row rank $p$ and the matrix $(\begin{matrix} A \\ B \end{matrix})$ has full column rank $n$ . (For linear least squares problems with inequality constraints, refer to Chapter E04.)
General Gauss–Markov linear model problems:
$minimize {‖y‖}_{2} subject to d = A x + B y,$
where $A$ is $m$ by $n$ and $B$ is $m$ by $p$ , with $n \leq m \leq n + p$ . When $B = I$ , the problem reduces to an ordinary linear least squares problem. When $B$ is square and nonsingular, it is equivalent to a weighted linear least squares problem:
$find x to minimize {‖B^{- 1} (d - A x)‖}_{2} .$
The problem has a unique solution on the assumptions that $A$ has full column rank $n$ , and the matrix $(A, B)$ has full row rank $m$ . Unless $B$ is diagonal, for numerical stability it is generally preferable to solve a weighted linear least squares problem as a general Gauss–Markov linear model problem.

The generalized $Q R$ (GQR) factorization of an

n

by

m

matrix

A

and an

n

by

p

matrix

B

is given by the pair of factorizations

A = Q R and B = Q T Z,

where

Q

and

Z

are respectively

n

by

n

and

p

by

p

orthogonal matrices (or unitary matrices if

A

and

B

are complex).

R

has the form

R = \begin{array}{rc} m \\ m & ( & R_{11} & ) \\ n - m & 0 \end{array}, if ​ n \geq m,

or

R = \begin{array}{rc} n & m - n \\ n & ( & R_{11} & R_{12} & ) \end{array}, if ​ n < m,

where

R_{11}

is upper triangular.

T

has the form

T = \begin{array}{rc} p - n & n \\ n & ( & 0 & T_{12} & ) \end{array}, if ​ n \leq p,

or

T = \begin{array}{rc} p \\ n - p & ( & T_{11} & ) \\ p & T_{21} \end{array}, if ​ n > p,

where

T_{12}

or

T_{21}

is upper triangular.

Note that if

B

is square and nonsingular, the GQR factorization of

A

and

B

implicitly gives the

Q R

factorization of the matrix

B^{- 1} A

:

B^{- 1} = Z^{T} (T^{- 1} R)

without explicitly computing the matrix inverse

B^{- 1}

or the product

B^{- 1} A

.

The GQR factorization can be used to solve the general (Gauss–Markov) linear model problem (GLM) (see Section 2.5, but note that

A

and

B

are dimensioned differently there as

m

by

n

and

p

by

n

respectively). Using the GQR factorization of

A

and

B

, we rewrite the equation

d = A x + B y

as

\begin{matrix} Q^{T} d & = Q^{T} A x + Q^{T} B y \\ = R x + T Z y . \end{matrix}

We partition this as

\begin{matrix} (\begin{matrix} d_{1} \\ d_{2} \end{matrix}) \end{matrix} = \begin{array}{rc} m \\ m & ( & R_{11} & ) \\ n - m & 0 \end{array} x + \begin{array}{rc} p - n + m & n - m \\ m & ( & T_{11} & T_{12} & ) \\ n - m & 0 & T_{22} \end{array} (\begin{matrix} y_{1} \\ y_{2} \end{matrix})

where

(\begin{matrix} d_{1} \\ d_{2} \end{matrix}) \equiv Q^{T} d,  and (\begin{matrix} y_{1} \\ y_{2} \end{matrix}) \equiv Z y .

The GLM problem is solved by setting

y_{1} = 0 and y_{2} = T_{22}^{- 1} d_{2}

from which we obtain the desired solutions

x = R_{11}^{- 1} (d_{1} - T_{12} y_{2}) and y = Z^{T} (\begin{matrix} 0 \\ y_{2} \end{matrix}) .

The generalized $R Q$ (GRQ) factorization of an

m

by

n

matrix

A

and a

p

by

n

matrix

B

is given by the pair of factorizations

A = R Q, B = Z T Q

where

Q

and

Z

are respectively

n

by

n

and

p

by

p

orthogonal matrices (or unitary matrices if

A

and

B

are complex).

R

has the form

R = \begin{array}{rc} n - m & m \\ m & ( & 0 & R_{12} & ) \end{array},  if ​ m \leq n,

or

R = \begin{array}{rc} n \\ m - n & ( & R_{11} & ) \\ n & R_{21} \end{array},  if ​ m > n,

where

R_{12}

or

R_{21}

is upper triangular.

T

has the form

T = \begin{array}{rc} n \\ n & ( & T_{11} & ) \\ p - n & 0 \end{array},   if ​ p \geq n,

or

T = \begin{array}{rc} p & n - p \\ p & ( & T_{11} & T_{12} & ) \end{array},   if ​ p < n,

where

T_{11}

is upper triangular.

Note that if

B

is square and nonsingular, the GRQ factorization of

A

and

B

implicitly gives the

R Q

factorization of the matrix

A B^{- 1}

:

A B^{- 1} = (R T^{- 1}) Z^{T}

without explicitly computing the matrix

B^{- 1}

or the product

A B^{- 1}

.

The GRQ factorization can be used to solve the linear equality-constrained least squares problem (LSE) (see Section 2.5). We use the GRQ factorization of

B

and

A

(note that

B

and

A

have swapped roles), written as

B = T Q and A = Z R Q .

We write the linear equality constraints

B x = d

as

T Q x = d,

which we partition as:

\begin{array}{rc} n - p & p \\ p & ( & 0 & T_{12} & ) \end{array} (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) = d where (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) \equiv Q_{x} .

Therefore

x_{2}

is the solution of the upper triangular system

T_{12} x_{2} = d .

Furthermore,

\begin{matrix} {‖A x - c‖}_{2} & = & {‖Z^{T} A x - Z^{T} c‖}_{2} \\ = & {‖R Q x - Z^{T} c‖}_{2} \end{matrix} .

We partition this expression as:

\begin{array}{rc} n - p & p \\ n - p & ( & R_{11} & R_{12} & ) \\ p + m - n & 0 & R_{22} \end{array} (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) - (\begin{matrix} c_{1} \\ c_{2} \end{matrix}),

where

(\begin{matrix} c_{1} \\ c_{2} \end{matrix}) \equiv Z^{T} c

.

To solve the LSE problem, we set

R_{11} x_{1} + R_{12} x_{2} - c_{1} = 0

which gives

x_{1}

as the solution of the upper triangular system

R_{11} x_{1} = c_{1} - R_{12} x_{2} .

Finally, the desired solution is given by

x = Q^{T} (\begin{matrix} x_{1} \\ x_{2} \end{matrix}) .

The generalized (or quotient) singular value decomposition of an

m

by

n

matrix

A

and a

p

by

n

matrix

B

is given by the pair of factorizations

A = U Σ_{1} [0, R] Q^{T} and B = V Σ_{2} [0, R] Q^{T} .

The matrices in these factorizations have the following properties:

–	$U$ is $m$ by $m$ , $V$ is $p$ by $p$ , $Q$ is $n$ by $n$ , and all three matrices are orthogonal. If $A$ and $B$ are complex, these matrices are unitary instead of orthogonal, and $Q^{T}$ should be replaced by $Q^{H}$ in the pair of factorizations.
–	$R$ is $r$ by $r$ , upper triangular and nonsingular. $[0, R]$ is $r$ by $n$ (in other words, the $0$ is an $r$ by $n - r$ zero matrix). The integer $r$ is the rank of $(\begin{matrix} A \\ B \end{matrix})$ , and satisfies $r \leq n$ .
–	$Σ_{1}$ is $m$ by $r$ , $Σ_{2}$ is $p$ by $r$ , both are real, non-negative and diagonal, and $Σ_{1}^{T} Σ_{1} + Σ_{2}^{T} Σ_{2} = I$ . Write $Σ_{1}^{T} Σ_{1} = diag (α_{1}^{2}, \dots, α_{r}^{2})$ and $Σ_{2}^{T} Σ_{2} = diag (β_{1}^{2}, \dots, β_{r}^{2})$ , where $α_{i}$ and $β_{i}$ lie in the interval from $0$ to $1$ . The ratios $α_{1} / β_{1}, \dots, α_{r} / β_{r}$ are called the generalized singular values of the pair $A$ , $B$ . If $β_{i} = 0$ , then the generalized singular value $α_{i} / β_{i}$ is infinite.

Σ_{1}

and

Σ_{2}

have the following detailed structures, depending on whether

m - r \geq 0

or

m - r < 0

. In the first case,

m - r \geq 0

, then

Σ_{1} = \begin{array}{rc} k & l \\ k & ( & I & 0 & ) \\ l & 0 & C \\ m - k - l & 0 & 0 \end{array} and Σ_{2} = \begin{array}{rc} k & l \\ l & ( & 0 & S & ) \\ p - l & 0 & 0 \end{array} .

Here

l

is the rank of

B

,

k = r - l

,

C

and

S

are diagonal matrices satisfying

C^{2} + S^{2} = I

, and

S

is nonsingular. We may also identify

α_{1} = \dots = α_{k} = 1

,

α_{k + i} = c_{i i}

, for

i = 1, 2, \dots, l

,

β_{1} = \dots = β_{k} = 0

, and

β_{k + i} = s_{i i}

, for

i = 1, 2, \dots, l

. Thus, the first

k

generalized singular values

α_{1} / β_{1}, \dots, α_{k} / β_{k}

are infinite, and the remaining

l

generalized singular values are finite.

In the second case, when

m - r < 0

,

Σ_{1} = \begin{array}{rc} k & m - k & k + l - m \\ k & ( & I & 0 & 0 & ) \\ m - k & 0 & C & 0 \end{array}

and

Σ_{2} = \begin{array}{rc} k & m - k & k + l - m \\ m - k & ( & 0 & S & 0 & ) \\ k + l - m & 0 & 0 & I \\ p - l & 0 & 0 & 0 \end{array} .

Again,

l

is the rank of

B

,

k = r - l

,

C

and

S

are diagonal matrices satisfying

C^{2} + S^{2} = I

, and

S

is nonsingular, and we may identify

α_{1} = \dots = α_{k} = 1

,

α_{k + i} = c_{i i}

, for

i = 1, 2, \dots, m - k

,

α_{m + 1} = \dots = α_{r} = 0

,

β_{1} = \dots = β_{k} = 0

,

β_{k + i} = s_{i i}

, for

i = 1, 2, \dots, m - k

and

β_{m + 1} = \dots = β_{r} = 1

. Thus, the first

k

generalized singular values

α_{1} / β_{1}, \dots, α_{k} / β_{k}

are infinite, and the remaining

l

generalized singular values are finite.

Here are some important special case of the generalized singular value decomposition. First, if

B

is square and nonsingular, then

r = n

and the generalized singular value decomposition of

A

and

B

is equivalent to the singular value decomposition of

A B^{- 1}

, where the singular values of

A B^{- 1}

are equal to the generalized singular values of the pair

A

,

B

:

A B^{- 1} = (U Σ_{1} R Q^{T}) {(V Σ_{2} R Q^{T})}^{- 1} = U (Σ_{1} Σ_{2}^{- 1}) V^{T} .

Second, if the columns of

{(A^{T} B^{T})}^{T}

are orthonormal, then

r = n

,

R = I

and the generalized singular value decomposition of

A

and

B

is equivalent to the CS (Cosine–Sine) decomposition of

{(A^{T} B^{T})}^{T}

:

(\begin{matrix} A \\ B \end{matrix}) = (\begin{matrix} U & 0 \\ 0 & V \end{matrix}) (\begin{matrix} Σ_{1} \\ Σ_{2} \end{matrix}) Q^{T} .

Third, the generalized eigenvalues and eigenvectors of

A^{T} A - λ B^{T} B

can be expressed in terms of the generalized singular value decomposition: Let

X = Q (\begin{matrix} I & 0 \\ 0 & R^{- 1} \end{matrix}) .

Then

X^{T} A^{T} A X = (\begin{matrix} 0 & 0 \\ 0 & Σ_{1}^{T} Σ_{1} \end{matrix}) and X^{T} B^{T} B X = (\begin{matrix} 0 & 0 \\ 0 & Σ_{2}^{T} Σ_{2} \end{matrix}) .

Therefore, the columns of

X

are the eigenvectors of

A^{T} A - λ B^{T} B

, and ‘nontrivial’ eigenvalues are the squares of the generalized singular values (see also Section 2.8). ‘Trivial’ eigenvalues are those corresponding to the leading

n - r

columns of

X

, which span the common null space of

A^{T} A

and

B^{T} B

. The ‘trivial eigenvalues’ are not well defined.

The symmetric eigenvalue problem is to find the eigenvalues,

λ

, and corresponding eigenvectors,

z \neq 0

, such that

A z = λ z, A = A^{T},   where ​ A ​ is real.

For the Hermitian eigenvalue problem we have

A z = λ z, A = A^{H},   where ​ A ​ is complex.

For both problems the eigenvalues

λ

are real.

When all eigenvalues and eigenvectors have been computed, we write

A = Z Λ Z^{T} (or ​ A = Z Λ Z^{H} ​ if complex),

where

Λ

is a diagonal matrix whose diagonal elements are the eigenvalues, and

Z

is an orthogonal (or unitary) matrix whose columns are the eigenvectors. This is the classical spectral factorization of

A

.

The basic task of the symmetric eigenproblem routines is to compute values of

λ

and, optionally, corresponding vectors

z

for a given matrix

A

. This computation proceeds in the following stages.

The real symmetric or complex Hermitian matrix $A$ is reduced to real tridiagonal form $T$ . If $A$ is real symmetric this decomposition is $A = Q T Q^{T}$ with $Q$ orthogonal and $T$ symmetric tridiagonal. If $A$ is complex Hermitian, the decomposition is $A = Q T Q^{H}$ with $Q$ unitary and $T$ , as before, real symmetric tridiagonal.
Eigenvalues and eigenvectors of the real symmetric tridiagonal matrix $T$ are computed. If all eigenvalues and eigenvectors are computed, this is equivalent to factorizing $T$ as $T = S Λ S^{T}$ , where $S$ is orthogonal and $Λ$ is diagonal. The diagonal entries of $Λ$ are the eigenvalues of $T$ , which are also the eigenvalues of $A$ , and the columns of $S$ are the eigenvectors of $T$ ; the eigenvectors of $A$ are the columns of $Z = Q S$ , so that $A = Z Λ Z^{T}$ ( $Z Λ Z^{H}$ when $A$ is complex Hermitian).

This chapter supports four primary algorithms for computing eigenvalues and eigenvectors of real symmetric matrices and complex Hermitian matrices. They are:

(i)	the divide-and-conquer algorithm;
(ii)	the $Q R$ algorithm;
(iii)	bisection followed by inverse iteration;
(iv)	the Relatively Robust Representation (RRR).

The divide-and-conquer algorithm is generally more efficient than the traditional

Q R

algorithm for computing all eigenvalues and eigenvectors, but the RRR algorithm tends to be fastest of all. For further information and references see Anderson et al. (1999).

This section is concerned with the solution of the generalized eigenvalue problems

A z = λ B z

,

A B z = λ z

, and

B A z = λ z

, where

A

and

B

are real symmetric or complex Hermitian and

B

is positive definite. Each of these problems can be reduced to a standard symmetric eigenvalue problem, using a Cholesky factorization of

B

as either

B = L L^{T}

or

B = U^{T} U

(

L L^{H}

or

U^{H} U

in the Hermitian case).

With

B = L L^{T}

, we have

A z = λ B z \Rightarrow (L^{- 1} A L^{- T}) (L^{T} z) = λ (L^{T} z) .

Hence the eigenvalues of

A z = λ B z

are those of

C y = λ y

, where

C

is the symmetric matrix

C = L^{- 1} A L^{- T}

and

y = L^{T} z

. In the complex case

C

is Hermitian with

C = L^{- 1} A L^{- H}

and

y = L^{H} z

.

Table 1 summarises how each of the three types of problem may be reduced to standard form

C y = λ y

, and how the eigenvectors

z

of the original problem may be recovered from the eigenvectors

y

of the reduced problem. The table applies to real problems; for complex problems, transposed matrices must be replaced by conjugate-transposes.

	Type of problem	Factorization of $B$	Reduction	Recovery of eigenvectors
1.	$A z = λ B z$	$B = L L^{T}$ , $B = U^{T} U$	$C = L^{- 1} A L^{- T}$ , $C = U^{- T} A U^{- 1}$	$z = L^{- T} y$ , $z = U^{- 1} y$
2.	$A B z = λ z$	$B = L L^{T}$ , $B = U^{T} U$	$C = L^{T} A L$ , $C = U A U^{T}$	$z = L^{- T} y$ , $z = U^{- 1} y$
3.	$B A z = λ z$	$B = L L^{T}$ , $B = U^{T} U$	$C = L^{T} A L$ , $C = U A U^{T}$	$z = L y$ , $z = U^{T} y$

Table 1
Reduction of generalized symmetric-definite eigenproblems to standard problems

When the generalized symmetric-definite problem has been reduced to the corresponding standard problem

C y = λ y

, this may then be solved using the routines described in the previous section. No special routines are needed to recover the eigenvectors

z

of the generalized problem from the eigenvectors

y

of the standard problem, because these computations are simple applications of Level 2 or Level 3 BLAS (see Chapter F06).

Routines which handle symmetric matrices are usually designed so that they use either the upper or lower triangle of the matrix; it is not necessary to store the whole matrix. If either the upper or lower triangle is stored conventionally in the upper or lower triangle of a two-dimensional array, the remaining elements of the array can be used to store other useful data. However, that is not always convenient, and if it is important to economize on storage, the upper or lower triangle can be stored in a one-dimensional array of length

n (n + 1) / 2

; that is, the storage is almost halved.

This storage format is referred to as packed storage; it is described in Section 3.3.2 in the F07 Chapter Introduction.

Routines designed for packed storage are usually less efficient, especially on high-performance computers, so there is a trade-off between storage and efficiency.

A band matrix is one whose elements are confined to a relatively small number of subdiagonals or superdiagonals on either side of the main diagonal. Algorithms can take advantage of bandedness to reduce the amount of work and storage required. The storage scheme for band matrices is described in Section 3.3.4 in the F07 Chapter Introduction.

If the problem is the generalized symmetric definite eigenvalue problem

A z = λ B z

and the matrices

A

and

B

are additionally banded, the matrix

C

as defined in Section 2.8 is, in general, full. We can reduce the problem to a banded standard problem by modifying the definition of

C

thus:

C = X^{T} A X,   where X = U^{- 1} Q or ​ L^{- T} Q,

where

Q

is an orthogonal matrix chosen to ensure that

C

has bandwidth no greater than that of

A

.

A further refinement is possible when

A

and

B

are banded, which halves the amount of work required to form

C

. Instead of the standard Cholesky factorization of

B

as

U^{T} U

or

L L^{T}

, we use a split Cholesky factorization

B = S^{T} S

, where

S = (\begin{array}{l} U_{11} \\ M_{21} & L_{22} \end{array})

with

U_{11}

upper triangular and

L_{22}

lower triangular of order approximately

n / 2

;

S

has the same bandwidth as

B

.

The nonsymmetric eigenvalue problem is to find the eigenvalues,

λ

, and corresponding eigenvectors,

v \neq 0

, such that

A v = λ v .

More precisely, a vector

v

as just defined is called a right eigenvector of

A

, and a vector

u \neq 0

satisfying

u^{T} A = λ u^{T} (u^{H} A = λ u^{H} when ​ u ​ is complex)

is called a left eigenvector of

A

.

A real matrix

A

may have complex eigenvalues, occurring as complex conjugate pairs.

This problem can be solved via the Schur factorization of

A

, defined in the real case as

A = Z T Z^{T},

where

Z

is an orthogonal matrix and

T

is an upper quasi-triangular matrix with

1

by

1

and

2

by

2

diagonal blocks, the

2

by

2

blocks corresponding to complex conjugate pairs of eigenvalues of

A

. In the complex case, the Schur factorization is

A = Z T Z^{H},

where

Z

is unitary and

T

is a complex upper triangular matrix.

The columns of

Z

are called the Schur vectors. For each

k

(

1 \leq k \leq n

), the first

k

columns of

Z

form an orthonormal basis for the invariant subspace corresponding to the first

k

eigenvalues on the diagonal of

T

. Because this basis is orthonormal, it is preferable in many applications to compute Schur vectors rather than eigenvectors. It is possible to order the Schur factorization so that any desired set of

k

eigenvalues occupy the

k

leading positions on the diagonal of

T

.

The two basic tasks of the nonsymmetric eigenvalue routines are to compute, for a given matrix

A

, all

n

values of

λ

and, if desired, their associated right eigenvectors

v

and/or left eigenvectors

u

, and the Schur factorization.

These two basic tasks can be performed in the following stages.

A general matrix $A$ is reduced to upper Hessenberg form $H$ which is zero below the first subdiagonal. The reduction may be written $A = Q H Q^{T}$ with $Q$ orthogonal if $A$ is real, or $A = Q H Q^{H}$ with $Q$ unitary if $A$ is complex.
The upper Hessenberg matrix $H$ is reduced to Schur form $T$ , giving the Schur factorization $H = S T S^{T}$ (for $H$ real) or $H = S T S^{H}$ (for $H$ complex). The matrix $S$ (the Schur vectors of $H$ ) may optionally be computed as well. Alternatively $S$ may be postmultiplied into the matrix $Q$ determined in stage $1$ , to give the matrix $Z = Q S$ , the Schur vectors of $A$ . The eigenvalues are obtained from the diagonal elements or diagonal blocks of $T$ .
Given the eigenvalues, the eigenvectors may be computed in two different ways. Inverse iteration can be performed on $H$ to compute the eigenvectors of $H$ , and then the eigenvectors can be multiplied by the matrix $Q$ in order to transform them to eigenvectors of $A$ . Alternatively the eigenvectors of $T$ can be computed, and optionally transformed to those of $H$ or $A$ if the matrix $S$ or $Z$ is supplied.

The accuracy with which eigenvalues can be obtained can often be improved by balancing a matrix. This is discussed further in Section 2.14.6 below.

The generalized nonsymmetric eigenvalue problem is to find the eigenvalues,

λ

, and corresponding eigenvectors,

v \neq 0

, such that

A v = λ B v .

More precisely, a vector

v

as just defined is called a right eigenvector of the matrix pair

(A, B)

, and a vector

u \neq 0

satisfying

u^{T} A = λ u^{T} B (u^{H} A = λ u^{H} B ​ when ​ u ​ is complex)

is called a left eigenvector of the matrix pair

(A, B)

.

If

B

is singular then the problem has one or more infinite eigenvalues

λ = \infty

, corresponding to

B v = 0

. Note that if

A

is nonsingular, then the equivalent problem

μ A v = B v

is perfectly well defined and an infinite eigenvalue corresponds to

μ = 0

. To deal with both finite (including zero) and infinite eigenvalues, the routines in this chapter do not compute

λ

explicitly, but rather return a pair of numbers

(α, β)

such that if

β \neq 0

λ = α / β

and if

α \neq 0

and

β = 0

then

λ = \infty

.

β

is always returned as real and non-negative. Of course, computationally an infinite eigenvalue may correspond to a small

β

rather than an exact zero.

For a given pair

(A, B)

the set of all the matrices of the form

(A - λ B)

is called a matrix pencil and

λ

and

v

are said to be an eigenvalue and eigenvector of the pencil

(A - λ B)

. If

A

and

B

are both singular and share a common null space then

\det (A - λ B) \equiv 0

so that the pencil

(A - λ B)

is singular for all

λ

. In other words any

λ

can be regarded as an eigenvalue. In exact arithmetic a singular pencil will have

α = β = 0

for some

(α, β)

. Computationally if some pair

(α, β)

is small then the pencil is singular, or nearly singular, and no reliance can be placed on any of the computed eigenvalues. Singular pencils can also manifest themselves in other ways; see, in particular, Sections 2.3.5.2 and 4.11.1.4 of Anderson et al. (1999) for further details.

The generalized eigenvalue problem can be solved via the generalized Schur factorization of the pair

(A, B)

defined in the real case as

A = Q S Z^{T}, B = Q T Z^{T},

where

Q

and

Z

are orthogonal,

T

is upper triangular with non-negative diagonal elements and

S

is upper quasi-triangular with

1

by

1

and

2

by

2

diagonal blocks, the

2

by

2

blocks corresponding to complex conjugate pairs of eigenvalues. In the complex case, the generalized Schur factorization is

A = Q S Z^{H}, B = Q T Z^{H},

where

Q

and

Z

are unitary and

S

and

T

are upper triangular, with

T

having real non-negative diagonal elements. The columns of

Q

and

Z

are called respectively the left and right generalized Schur vectors and span pairs of deflating subspaces of

A

and

B

, which are a generalization of invariant subspaces.

It is possible to order the generalized Schur factorization so that any desired set of

k

eigenvalues correspond to the

k

leading positions on the diagonals of the pair

(S, T)

.

The two basic tasks of the generalized nonsymmetric eigenvalue routines are to compute, for a given pair

(A, B)

, all

n

values of

λ

and, if desired, their associated right eigenvectors

v

and/or left eigenvectors

u

, and the generalized Schur factorization.

These two basic tasks can be performed in the following stages.

The matrix pair $(A, B)$ is reduced to generalized upper Hessenberg form $(H, R)$ , where $H$ is upper Hessenberg (zero below the first subdiagonal) and $R$ is upper triangular. The reduction may be written as $A = Q_{1} H Z_{1}^{T}, B = Q_{1} R Z_{1}^{T}$ in the real case with $Q_{1}$ and $Z_{1}$ orthogonal, and $A = Q_{1} H Z_{1}^{H}, B = Q_{1} R Z_{1}^{H}$ in the complex case with $Q_{1}$ and $Z_{1}$ unitary.
The generalized upper Hessenberg form $(H, R)$ is reduced to the generalized Schur form $(S, T)$ using the generalized Schur factorization $H = Q_{2} S Z_{2}^{T}$ , $R = Q_{2} T Z_{2}^{T}$ in the real case with $Q_{2}$ and $Z_{2}$ orthogonal, and $H = Q_{2} S Z_{2}^{H}, R = Q_{2} T Z_{2}^{H}$ in the complex case. The generalized Schur vectors of $(A, B)$ are given by $Q = Q_{1} Q_{2}$ , $Z = Z_{1} Z_{2}$ . The eigenvalues are obtained from the diagonal elements (or blocks) of the pair $(S, T)$ .
Given the eigenvalues, the eigenvectors of the pair $(S, T)$ can be computed, and optionally transformed to those of $(H, R)$ or $(A, B)$ .

The accuracy with which eigenvalues can be obtained can often be improved by balancing a matrix pair. This is discussed further in Section 2.14.8 below.

The Sylvester equation is a matrix equation of the form

A X + X B = C,

where

A

,

B

, and

C

are given matrices with

A

being

m

by

m

,

B

an

n

by

n

matrix and

C

, and the solution matrix

X

,

m

by

n

matrices. The solution of a special case of this equation occurs in the computation of the condition number for an invariant subspace, but a combination of routines in this chapter allows the solution of the general Sylvester equation.

Routines are also provided for solving a special case of the generalized Sylvester equations

A R - L B = C, D R - L E = F,

where

(A, D)

,

(B, E)

and

(C, F)

are given matrix pairs, and

R

and

L

are the solution matrices.

In this section we discuss the effects of rounding errors in the solution process and the effects of uncertainties in the data, on the solution to the problem. A number of the routines in this chapter return information, such as condition numbers, that allow these effects to be assessed. First we discuss some notation used in the error bounds of later sections.

The bounds usually contain the factor

p (n)

(or

p (m, n)

), which grows as a function of the matrix dimension

n

(or matrix dimensions

m

and

n

). It measures how errors can grow as a function of the matrix dimension, and represents a potentially different function for each problem. In practice, it usually grows just linearly;

p (n) \leq 10 n

is often true, although generally only much weaker bounds can be actually proved. We normally describe

p (n)

as a ‘modestly growing’ function of

n

. For detailed derivations of various

p (n)

, see Golub and Van Loan (1996) and Wilkinson (1965).

For linear equation (see Chapter F07) and least squares solvers, we consider bounds on the relative error

‖x - \hat{x}‖ / ‖x‖

in the computed solution

\hat{x}

, where

x

is the true solution. For eigenvalue problems we consider bounds on the error

|λ_{i} - {\hat{λ}}_{i}|

in the

i

th computed eigenvalue

{\hat{λ}}_{i}

, where

λ_{i}

is the true

i

th eigenvalue. For singular value problems we similarly consider bounds

|σ_{i} - {\hat{σ}}_{i}|

.

Bounding the error in computed eigenvectors and singular vectors

{\hat{v}}_{i}

is more subtle because these vectors are not unique: even though we restrict

{‖{\hat{v}}_{i}‖}_{2} = 1

and

{‖v_{i}‖}_{2} = 1

, we may still multiply them by arbitrary constants of absolute value

1

. So to avoid ambiguity we bound the angular difference between

{\hat{v}}_{i}

and the true vector

v_{i}

, so that

\begin{array}{lcl} θ (v_{i}, {\hat{v}}_{i}) & = & acute angle between ​ v_{i} ​ and ​ {\hat{v}}_{i} \\ = & \arccos |v_{i}^{H} {\hat{v}}_{i}| . \end{array}

(2)

Here

\arccos (θ)

is in the standard range:

0 \leq \arccos (θ) < π

. When

θ (v_{i}, {\hat{v}}_{i})

is small, we can choose a constant

α

with absolute value

1

so that

{‖α v_{i} - {\hat{v}}_{i}‖}_{2} \approx θ (v_{i}, {\hat{v}}_{i})

.

In addition to bounds for individual eigenvectors, bounds can be obtained for the spaces spanned by collections of eigenvectors. These may be much more accurately determined than the individual eigenvectors which span them. These spaces are called invariant subspaces in the case of eigenvectors, because if

v

is any vector in the space,

A v

is also in the space, where

A

is the matrix. Again, we will use angle to measure the difference between a computed space

\hat{S}

and the true space

S

:

\begin{array}{lcl} θ (S, \hat{S}) & = & acute angle between ​ S ​ and ​ \hat{S} \\ = & \max_{\begin{matrix} s \in S \\ s \neq 0 \end{matrix}} \min_{\begin{matrix} \hat{s} \in \hat{S} \\ \hat{s} \neq 0 \end{matrix}} θ (s, \hat{s}) or \max_{\begin{matrix} \hat{s} \in \hat{S} \\ \hat{s} \neq 0 \end{matrix}} \min_{\begin{matrix} s \in S \\ s \neq 0 \end{matrix}} θ (s, \hat{s}) \end{array}

(3)

θ (S, \hat{S})

may be computed as follows. Let

S

be a matrix whose columns are orthonormal and

span S

. Similarly let

\hat{S}

be an orthonormal matrix with columns spanning

\hat{S}

. Then

θ (S, \hat{S}) = \arccos σ_{\min} (S^{H} \hat{S}) .

Finally, we remark on the accuracy of the bounds when they are large. Relative errors like

‖\hat{x} - x‖ / ‖x‖

and angular errors like

θ ({\hat{v}}_{i}, v_{i})

are only of interest when they are much less than

1

. Some stated bounds are not strictly true when they are close to

1

, but rigorous bounds are much more complicated and supply little extra information in the interesting case of small errors. These bounds are indicated by using the symbol

≲

, or ‘approximately less than’, instead of the usual

\leq

. Thus, when these bounds are close to 1 or greater, they indicate that the computed answer may have no significant digits at all, but do not otherwise bound the error.

A number of routines in this chapter return error estimates and/or condition number estimates directly. In other cases Anderson et al. (1999) gives code fragments to illustrate the computation of these estimates, and a number of the Chapter F08 example programs, for the driver routines, implement these code fragments.

The conventional error analysis of linear least squares problems goes as follows. The problem is to find the

x

minimizing

{‖A x - b‖}_{2}

. Let

\hat{x}

be the solution computed using one of the methods described above. We discuss the most common case, where

A

is overdetermined (i.e., has more rows than columns) and has full rank.

Then the computed solution

\hat{x}

has a small normwise backward error. In other words

\hat{x}

minimizes

{‖(A + E) \hat{x} - (b + f)‖}_{2}

, where

\max (\frac{{‖E‖}_{2}}{{‖A‖}_{2}}, \frac{{‖f‖}_{2}}{{‖b‖}_{2}}) \leq p (n) ε

and

p (n)

is a modestly growing function of

n

and

ε

is the machine precision. Let

κ_{2} (A) = σ_{\max} (A) / σ_{\min} (A)

,

ρ = {‖A x - b‖}_{2}

, and

\sin (θ) = ρ / {‖b‖}_{2}

. Then if

p (n) ε

is small enough, the error

\hat{x} - x

is bounded by

\frac{{‖x - \hat{x}‖}_{2}}{{‖x‖}_{2}} ≲ p (n) ε \{\frac{2 κ_{2} (A)}{\cos (θ)} + \tan (θ) κ_{2}^{2} (A)\} .

If

A

is rank-deficient, the problem can be regularized by treating all singular values less than a user-specified threshold as exactly zero. See Golub and Van Loan (1996) for error bounds in this case, as well as for the underdetermined case.

The solution of the overdetermined, full-rank problem may also be characterised as the solution of the linear system of equations

(\begin{matrix} I & A \\ A^{T} & 0 \end{matrix}) (\begin{matrix} r \\ x \end{matrix}) = (\begin{matrix} b \\ 0 \end{matrix}) .

By solving this linear system (see Chapter F07) component-wise error bounds can also be obtained Arioli et al. (1989).

The usual error analysis of the SVD algorithm is as follows (see Golub and Van Loan (1996)).

The computed SVD,

\hat{U} \hat{Σ} {\hat{V}}^{T}

, is nearly the exact SVD of

A + E

, i.e.,

A + E = (\hat{U} + δ \hat{U}) \hat{Σ} (\hat{V} + δ \hat{V})

is the true SVD, so that

\hat{U} + δ \hat{U}

and

\hat{V} + δ \hat{V}

are both orthogonal, where

{‖E‖}_{2} / {‖A‖}_{2} \leq p (m, n) ε

,

‖δ \hat{U}‖ \leq p (m, n) ε

, and

‖δ \hat{V}‖ \leq p (m, n) ε

. Here

p (m, n)

is a modestly growing function of

m

and

n

and

ε

is the machine precision. Each computed singular value

{\hat{σ}}_{i}

differs from the true

σ_{i}

by an amount satisfying the bound

|{\hat{σ}}_{i} - σ_{i}| \leq p (m, n) ε σ_{1} .

Thus large singular values (those near

σ_{1}

) are computed to high relative accuracy and small ones may not be.

The angular difference between the computed left singular vector

{\hat{u}}_{i}

and the true

u_{i}

satisfies the approximate bound

θ ({\hat{u}}_{i}, u_{i}) ≲ \frac{p (m, n) ε {‖A‖}_{2}}{{gap}_{i}}

where

{gap}_{i} = \min_{j \neq i} |σ_{i} - σ_{j}|

is the absolute gap between

σ_{i}

and the nearest other singular value. Thus, if

σ_{i}

is close to other singular values, its corresponding singular vector

u_{i}

may be inaccurate. The same bound applies to the computed right singular vector

{\hat{v}}_{i}

and the true vector

v_{i}

. The gaps may be easily obtained from the computed singular values.

Let

\hat{S}

be the space spanned by a collection of computed left singular vectors

\{{\hat{u}}_{i}, i \in I\}

, where

I

is a subset of the integers from

1

to

n

. Let

S

be the corresponding true space. Then

θ (\hat{S}, S) ≲ \frac{p (m, n) ε {‖A‖}_{2}}{{gap}_{I}} .

where

{gap}_{I} = \min \{|σ_{i} - σ_{j}| for ​ i \in I, j \notin I\}

is the absolute gap between the singular values in

I

and the nearest other singular value. Thus, a cluster of close singular values which is far away from any other singular value may have a well determined space

\hat{S}

even if its individual singular vectors are ill-conditioned. The same bound applies to a set of right singular vectors

\{{\hat{v}}_{i}, i \in I\}

.

In the special case of bidiagonal matrices, the singular values and singular vectors may be computed much more accurately (see Demmel and Kahan (1990)). A bidiagonal matrix

B

has nonzero entries only on the main diagonal and the diagonal immediately above it (or immediately below it). Reduction of a dense matrix to bidiagonal form

B

can introduce additional errors, so the following bounds for the bidiagonal case do not apply to the dense case.

Using the routines in this chapter, each computed singular value of a bidiagonal matrix is accurate to nearly full relative accuracy, no matter how tiny it is, so that

|{\hat{σ}}_{i} - σ_{i}| \leq p (m, n) ε σ_{i} .

The computed left singular vector

{\hat{u}}_{i}

has an angular error at most about

θ ({\hat{u}}_{i}, u_{i}) ≲ \frac{p (m, n) ε}{{relgap}_{i}}

where

{relgap}_{i} = \min_{j \neq i} |σ_{i} - σ_{j}| / (σ_{i} + σ_{j})

is the relative gap between

σ_{i}

and the nearest other singular value. The same bound applies to the right singular vector

{\hat{v}}_{i}

and

v_{i}

. Since the relative gap may be much larger than the absolute gap, this error bound may be much smaller than the previous one. The relative gaps may be easily obtained from the computed singular values.

The usual error analysis of the symmetric eigenproblem is as follows (see Parlett (1998)).

The computed eigendecomposition

\hat{Z} \hat{Λ} {\hat{Z}}^{T}

is nearly the exact eigendecomposition of

A + E

, i.e.,

A + E = (\hat{Z} + δ \hat{Z}) \hat{Λ} {(\hat{Z} + δ \hat{Z})}^{T}

is the true eigendecomposition so that

\hat{Z} + δ \hat{Z}

is orthogonal, where

{‖E‖}_{2} / {‖A‖}_{2} \leq p (n) ε

and

{‖δ \hat{Z}‖}_{2} \leq p (n) ε

and

p (n)

is a modestly growing function of

n

and

ε

is the machine precision. Each computed eigenvalue

{\hat{λ}}_{i}

differs from the true

λ_{i}

by an amount satisfying the bound

|{\hat{λ}}_{i} - λ_{i}| \leq p (n) ε {‖A‖}_{2} .

Thus large eigenvalues (those near

\max_{i} |λ_{i}| = {‖A‖}_{2}

) are computed to high relative accuracy and small ones may not be.

The angular difference between the computed unit eigenvector

{\hat{z}}_{i}

and the true

z_{i}

satisfies the approximate bound

θ ({\hat{z}}_{i}, z_{i}) ≲ \frac{p (n) ε {‖A‖}_{2}}{{gap}_{i}}

if

p (n) ε

is small enough, where

{gap}_{i} = \min_{j \neq i} |λ_{i} - λ_{j}|

is the absolute gap between

λ_{i}

and the nearest other eigenvalue. Thus, if

λ_{i}

is close to other eigenvalues, its corresponding eigenvector

z_{i}

may be inaccurate. The gaps may be easily obtained from the computed eigenvalues.

Let

\hat{S}

be the invariant subspace spanned by a collection of eigenvectors

\{{\hat{z}}_{i}, i \in I\}

, where

I

is a subset of the integers from

1

to

n

. Let

S

be the corresponding true subspace. Then

θ (\hat{S}, S) ≲ \frac{p (n) ε {‖A‖}_{2}}{{gap}_{I}}

where

{gap}_{I} = \min \{|λ_{i} - λ_{j}| for ​ i \in I, j \notin I\}

is the absolute gap between the eigenvalues in

I

and the nearest other eigenvalue. Thus, a cluster of close eigenvalues which is far away from any other eigenvalue may have a well determined invariant subspace

\hat{S}

even if its individual eigenvectors are ill-conditioned.

In the special case of a real symmetric tridiagonal matrix

T

, routines in this chapter can compute the eigenvalues and eigenvectors much more accurately. See Anderson et al. (1999) for further details.

The three types of problem to be considered are

A - λ B

,

A B - λ I

and

B A - λ I

. In each case

A

and

B

are real symmetric (or complex Hermitian) and

B

is positive definite. We consider each case in turn, assuming that routines in this chapter are used to transform the generalized problem to the standard symmetric problem, followed by the solution of the symmetric problem. In all cases

{gap}_{i} = \min_{j \neq i} |λ_{i} - λ_{j}|

is the absolute gap between

λ_{i}

and the nearest other eigenvalue.

$A - λ B$ . The computed eigenvalues ${\hat{λ}}_{i}$ can differ from the true eigenvalues $λ_{i}$ by an amount
$|{\hat{λ}}_{i} - λ_{i}| ≲ p (n) ε {‖B^{- 1}‖}_{2} {‖A‖}_{2} .$
The angular difference between the computed eigenvector ${\hat{z}}_{i}$ and the true eigenvector $z_{i}$ is
$θ ({\hat{z}}_{i}, z_{i}) ≲ \frac{p (n) ε {‖B^{- 1}‖}_{2} {‖A‖}_{2} {(κ_{2} (B))}^{1 / 2}}{{gap}_{i}} .$
$A B - λ I$ or $B A - λ I$ . The computed eigenvalues ${\hat{λ}}_{i}$ can differ from the true eigenvalues $λ_{i}$ by an amount
$|{\hat{λ}}_{i} - λ_{i}| ≲ p (n) ε {‖B‖}_{2} {‖A‖}_{2} .$
The angular difference between the computed eigenvector ${\hat{z}}_{i}$ and the true eigenvector $z_{i}$ is
$θ ({\hat{z}}_{i}, z_{i}) ≲ \frac{p (n) ε {‖B‖}_{2} {‖A‖}_{2} {(κ_{2} (B))}^{1 / 2}}{{gap}_{i}} .$

These error bounds are large when

B

is ill-conditioned with respect to inversion (

κ_{2} (B)

is large). It is often the case that the eigenvalues and eigenvectors are much better conditioned than indicated here. One way to get tighter bounds is effective when the diagonal entries of

B

differ widely in magnitude, as for example with a graded matrix.

$A - λ B$ . Let $D = diag (b_{11}^{- 1 / 2}, \dots, b_{n n}^{- 1 / 2})$ be a diagonal matrix. Then replace $B$ by $D B D$ and $A$ by $D A D$ in the above bounds.
$A B - λ I$ or $B A - λ I$ . Let $D = diag (b_{11}^{- 1 / 2}, \dots, b_{n n}^{- 1 / 2})$ be a diagonal matrix. Then replace $B$ by $D B D$ and $A$ by $D^{- 1} A D^{- 1}$ in the above bounds.

Further details can be found in Anderson et al. (1999).

The nonsymmetric eigenvalue problem is more complicated than the symmetric eigenvalue problem. In this section, we just summarise the bounds. Further details can be found in Anderson et al. (1999).

We let

{\hat{λ}}_{i}

be the

i

th computed eigenvalue and

λ_{i}

the

i

th true eigenvalue. Let

{\hat{v}}_{i}

be the corresponding computed right eigenvector, and

v_{i}

the true right eigenvector (so

A v_{i} = λ_{i} v_{i}

). If

I

is a subset of the integers from

1

to

n

, we let

λ_{I}

denote the average of the selected eigenvalues:

λ_{I} = (\sum_{i \in I} λ_{i}) / (\sum_{i \in I} 1)

, and similarly for

{\hat{λ}}_{I}

. We also let

S_{I}

denote the subspace spanned by

\{v_{i}, i \in I\}

; it is called a right invariant subspace because if

v

is any vector in

S_{I}

then

A v

is also in

S_{I}

.

{\hat{S}}_{I}

is the corresponding computed subspace.

The algorithms for the nonsymmetric eigenproblem are normwise backward stable: they compute the exact eigenvalues, eigenvectors and invariant subspaces of slightly perturbed matrices

(A + E) E

, where

‖E‖ \leq p (n) ε ‖A‖

. Some of the bounds are stated in terms of

{‖E‖}_{2}

and others in terms of

{‖E‖}_{F}

; one may use

p (n) ε

for either quantity.

Routines are provided so that, for each (

{\hat{λ}}_{i}, {\hat{v}}_{i}

) pair the two values

s_{i}

and

{sep}_{i}

, or for a selected subset

I

of eigenvalues the values

s_{I}

and

{sep}_{I}

can be obtained, for which the error bounds in Table 2 are true for sufficiently small

‖E‖

, (which is why they are called asymptotic):

Simple eigenvalue	$\|{\hat{λ}}_{i} - λ_{i}\| ≲ {‖E‖}_{2} / s_{i}$
Eigenvalue cluster	$\|{\hat{λ}}_{I} - λ_{I}\| ≲ {‖E‖}_{2} / s_{I}$
Eigenvector	$θ ({\hat{v}}_{i}, v_{i}) ≲ {‖E‖}_{F} / {sep}_{i}$
Invariant subspace	$θ ({\hat{S}}_{I}, S_{I}) ≲ {‖E‖}_{F} / {sep}_{I}$

Table 2
Asymptotic error bounds for the nonsymmetric eigenproblem

If the problem is ill-conditioned, the asymptotic bounds may only hold for extremely small

‖E‖

. The global error bounds of Table 3 are guaranteed to hold for all

{‖E‖}_{F} < s \times sep / 4

:

Simple eigenvalue	$\|{\hat{λ}}_{i} - λ_{i}\| \leq n {‖E‖}_{2} / s_{i}$	Holds for all $E$
Eigenvalue cluster	$\|{\hat{λ}}_{I} - λ_{I}\| \leq 2 {‖E‖}_{2} / s_{I}$	Requires ${‖E‖}_{F} < s_{I} \times {sep}_{I} / 4$
Eigenvector	$θ ({\hat{v}}_{i}, v_{i}) \leq \arctan (2 {‖E‖}_{F} / ({sep}_{i} - 4 {‖E‖}_{F} / s_{i}))$	Requires ${‖E‖}_{F} < s_{i} \times {sep}_{i} / 4$
Invariant subspace	$θ ({\hat{S}}_{I}, S_{I}) \leq \arctan (2 {‖E‖}_{F} / ({sep}_{I} - 4 {‖E‖}_{F} / s_{I}))$	Requires ${‖E‖}_{F} < s_{I} \times {sep}_{I} / 4$

Table 3
Global error bounds for the nonsymmetric eigenproblem

There are two preprocessing steps one may perform on a matrix

A

in order to make its eigenproblem easier. The first is permutation, or reordering the rows and columns to make

A

more nearly upper triangular (closer to Schur form):

A^{'} = P A P^{T}

, where

P

is a permutation matrix. If

A^{'}

is permutable to upper triangular form (or close to it), then no floating point operations (or very few) are needed to reduce it to Schur form. The second is scaling by a diagonal matrix

D

to make the rows and columns of

A^{'}

more nearly equal in norm:

A^{''} = D A^{'} D^{- 1}

. Scaling can make the matrix norm smaller with respect to the eigenvalues, and so possibly reduce the inaccuracy contributed by roundoff (see Chapter 11 of Wilkinson and Reinsch (1971)). We refer to these two operations as balancing.

Permuting has no effect on the condition numbers or their interpretation as described previously. Scaling, however, does change their interpretation and further details can be found in Anderson et al. (1999).

The algorithms for the generalized nonsymmetric eigenvalue problem are normwise backward stable: they compute the exact eigenvalues (as the pairs

(α, β)

), eigenvectors and deflating subspaces of slightly perturbed pairs

(A + E, B + F)

, where

{‖(E, F)‖}_{F} \leq p (n) ε {‖(A, B)‖}_{F} .

Asymptotic and global error bounds can be obtained, which are generalizations of those given in Tables 2 and 3. See Section 4.11 of Anderson et al. (1999) for details. Routines are provided to compute estimates of reciprocal conditions numbers for eigenvalues and eigenspaces.

As with the standard nonsymmetric eigenvalue problem, there are two preprocessing steps one may perform on a matrix pair

(A, B)

in order to make its eigenproblem easier; permutation and scaling, which together are referred to as balancing, as indicated in the following two steps.

The balancing routine first attempts to permute

A

and

B

to block upper triangular form by a similarity transformation:

\begin{matrix} P A P^{T} = F = (\begin{matrix} F_{11} & F_{12} & F_{13} \\ F_{22} & F_{23} \\ F_{33} \end{matrix}), \\ P B P^{T} = G = (\begin{matrix} G_{11} & G_{12} & G_{13} \\ G_{22} & G_{23} \\ G_{33} \end{matrix}), \end{matrix}

where

P

is a permutation matrix,

F_{11}

,

F_{33}

,

G_{11}

and

G_{33}

are upper triangular. Then the diagonal elements of the matrix

(F_{11}, G_{11})

and

(G_{33}, H_{33})

are generalized eigenvalues of

(A, B)

. The rest of the generalized eigenvalues are given by the matrix pair

(F_{22}, G_{22})

. Subsequent operations to compute the eigenvalues of

(A, B)

need only be applied to the matrix

(F_{22}, G_{22})

; this can save a significant amount of work if

(F_{22}, G_{22})

is smaller than the original matrix pair

(A, B)

. If no suitable permutation exists (as is often the case), then there is no gain in efficiency or accuracy.

The balancing routine applies a diagonal similarity transformation to

(F, G)

, to make the rows and columns of

(F_{22}, G_{22})

as close as possible in the norm:

\begin{matrix} D F D^{- 1} = (\begin{matrix} I \\ D_{22} \\ I \end{matrix}) (\begin{matrix} F_{11} & F_{12} & F_{13} \\ F_{22} & F_{23} \\ F_{33} \end{matrix}) (\begin{matrix} I \\ D_{22}^{- 1} \\ I \end{matrix}), \\ D G D^{- 1} = (\begin{matrix} I \\ D_{22} \\ I \end{matrix}) (\begin{matrix} G_{11} & G_{12} & G_{13} \\ G_{22} & G_{23} \\ G_{33} \end{matrix}) (\begin{matrix} I \\ D_{22}^{- 1} \\ I \end{matrix}) . \end{matrix}

This transformation usually improves the accuracy of computed generalized eigenvalues and eigenvectors. However, there are exceptional occasions when this transformation increases the norm of the pencil; in this case accuracy could be lower with diagonal balancing.

See Anderson et al. (1999) for further details.

Error bounds for other problems such as the generalized linear least squares problem and generalized singular value decomposition can be found in Anderson et al. (1999).

A number of the routines in this chapter use what is termed a block partitioned algorithm. This means that at each major step of the algorithm a block of rows or columns is updated, and much of the computation is performed by matrix-matrix operations on these blocks. The matrix-matrix operations are performed by calls to the Level 3 BLAS (see Chapter F06), which are the key to achieving high performance on many modern computers. In the case of the

Q R

algorithm for reducing an upper Hessenberg matrix to Schur form, a multishift strategy is used in order to improve performance. See Golub and Van Loan (1996) or Anderson et al. (1999) for more about block partitioned algorithms and the multishift strategy.

The performance of a block partitioned algorithm varies to some extent with the block size – that is, the number of rows or columns per block. This is a machine-dependent parameter, which is set to a suitable value when the library is implemented on each range of machines. You do not normally need to be aware of what value is being used. Different block sizes may be used for different routines. Values in the range

16

to

64

are typical.

On more conventional machines there is often no advantage from using a block partitioned algorithm, and then the routines use an unblocked algorithm (effectively a block size of

1

), relying solely on calls to the Level 2 BLAS (see Chapter F06 again).

The only situation in which you need some awareness of the block size is when it affects the amount of workspace to be supplied to a particular routine. This is discussed in Section 3.4.3.

The tables in the following sub-sections show the routines which are provided for performing different computations on different types of matrices. Each entry in the table gives the NAG routine name and the LAPACK double precision name (see Section 3.2).

Black box (or driver) routines are provided for the solution of most problems. In a number of cases there are simple drivers, which just return the solution to the problem, as well as expert drivers, which return additional information, such as condition number estimates, and may offer additional facilities such as balancing. The following sub-sections give tables for the driver routines.

It is possible to solve problems by calling two or more routines in sequence. Some common sequences of routines are indicated in the tables in the following sub-sections; an asterisk (

*

) against a routine name means that the sequence of calls is illustrated in the example program for that routine.

Routines are provided for

Q R

factorization (with and without column pivoting), and for

L Q

,

Q L

and

R Q

factorizations (without pivoting only), of a general real or complex rectangular matrix. A routine is also provided for the

R Q

factorization of a real or complex upper trapezoidal matrix. (LAPACK refers to this as the

R Z

factorization.)

The factorization routines do not form the matrix

Q

explicitly, but represent it as a product of elementary reflectors (see Section 3.3.6). Additional routines are provided to generate all or part of

Q

explicitly if it is required, or to apply

Q

in its factored form to another matrix (specifically to compute one of the matrix products

Q C

,

Q^{T} C

,

C Q

or

C Q^{T}

with

Q^{T}

replaced by

Q^{H}

if

C

and

Q

are complex).

	Factorize without pivoting	Factorize with pivoting	Generate matrix $Q$	Apply matrix $Q$
$Q R$ factorization, real matrices	F08AEF (DGEQRF)	F08BFF (DGEQP3)	F08AFF (DORGQR)	F08AGF (DORMQR)
$L Q$ factorization, real matrices	F08AHF (DGELQF)		F08AJF (DORGLQ)	F08AKF (DORMLQ)
$Q L$ factorization, real matrices	F08CEF (DGEQLF)		F08CFF (DORGQL)	F08CGF (DORMQL)
$R Q$ factorization, real matrices	F08CHF (DGERQF)		F08CJF (DORGRQ)	F08CKF (DORMRQ)
$R Q$ factorization, real upper trapezoidal matrices	F08BHF (DTZRZF)			F08BKF (DORMRZ)
$Q R$ factorization, complex matrices	F08ASF (ZGEQRF)	F08BTF (ZGEQP3)	F08ATF (ZUNGQR)	F08AUF (ZUNMQR)
$L Q$ factorization, complex matrices	F08AVF (ZGELQF)		F08AWF (ZUNGLQ)	F08AXF (ZUNMLQ)
$Q L$ factorization, complex matrices	F08CSF (ZGEQLF)		F08CTF (ZUNGQL)	F08CUF (ZUNMQL)
$R Q$ factorization, complex matrices	F08CVF (ZGERQF)		F08CWF (ZUNGRQ)	F08CXF (ZUNMRQ)
$R Q$ factorization, complex upper trapezoidal matrices	F08BVF (ZTZRZF)			F08BXF (ZUNMRZ)

To solve linear least squares problems, as described in Sections 2.2.1 or 2.2.3, routines based on the

Q R

factorization can be used:

real data, full-rank problem	F08AEF*, F06YJF, F08AGF
complex data, full-rank problem	F08ASF*, F06ZJF, F08AUF
real data, rank-deficient problem	F08BFF*, F06YJF, F08AGF
complex data, rank-deficient problem	F08BTF*, F06ZJF, F08AUF

To find the minimum norm solution of under-determined systems of linear equations, as described in Section 2.2.2, routines based on the

L Q

factorization can be used:

real data, full-rank problem	F08AHF*, F06YJF, F08AKF
complex data, full-rank problem	F08AVF*, F06ZJF, F08AXF

Routines are provided for the generalized

Q R

and

R Q

factorizations of real and complex matrix pairs.

	Factorize
Generalized $Q R$ factorization, real matrices	F08ZEF (DGGQRF)
Generalized $R Q$ factorization, real matrices	F08ZFF (DGGRQF)
Generalized $Q R$ factorization, complex matrices	F08ZSF (ZGGQRF)
Generalized $R Q$ factorization, complex matrices	F08ZTF (ZGGRQF)

Routines are provided to reduce a general real or complex rectangular matrix

A

to real bidiagonal form

B

by an orthogonal transformation

A = Q B P^{T}

(or by a unitary transformation

A = Q B P^{H}

if

A

is complex). Different routines allow a full matrix

A

to be stored conventionally (see Section 3.3.1), or a band matrix to use band storage (see Section 3.3.4 in the F07 Chapter Introduction).

The routines for reducing full matrices do not form the matrix

Q

or

P

explicitly; additional routines are provided to generate all or part of them, or to apply them to another matrix, as with the routines for orthogonal factorizations. Explicit generation of

Q

or

P

is required before using the bidiagonal

Q R

algorithm to compute left or right singular vectors of

A

.

The routines for reducing band matrices have options to generate

Q

or

P

if required.

Further routines are provided to compute all or part of the singular value decomposition of a real bidiagonal matrix; the same routines can be used to compute the singular value decomposition of a real or complex matrix that has been reduced to bidiagonal form.

	Reduce to bidiagonal form	Generate matrix $Q$ or $P^{T}$	Apply matrix $Q$ or $P$	Reduce band matrix to bidiagonal form	SVD of bidiagonal form ( $QR$ algorithm)	SVD of bidiagonal form (divide and conquer)
real matrices	F08KEF (DGEBRD)	F08KFF (DORGBR)	F08KGF (DORMBR)	F08LEF (DGBBRD)	F08MEF (DBDSQR)	F08MDF (DBDSDC)
complex matrices	F08KSF (ZGEBRD)	F08KTF (ZUNGBR)	F08KUF (ZUNMBR)	F08LSF (ZGBBRD)	F08MSF (ZBDSQR)

Given the singular values, F08FLF (DDISNA) is provided to compute the reciprocal condition numbers for the left or right singular vectors of a real or complex matrix.

To compute the singular values and vectors of a rectangular matrix, as described in Section 2.3, use the following sequence of calls:

Rectangular matrix (standard storage)

real matrix, singular values and vectors	F08KEF, F08KFF*, F08MEF
complex matrix, singular values and vectors	F08KSF, F08KTF*, F08MSF

Rectangular matrix (banded)

real matrix, singular values and vectors	F08LEF, F08KFF, F08MEF
complex matrix, singular values and vectors	F08LSF, F08KTF, F08MSF

To use the singular value decomposition to solve a linear least squares problem, as described in Section 2.4, the following routines are required:

real data	F06YAF, F08KEF, F08KFF, F08KGF, F08MEF
complex data	F06ZAF, F08KSF, F08KTF, F08KUF, F08MSF

Routines are provided to compute the generalized SVD of a real or complex matrix pair

(A, B)

in upper trapezoidal form. Routines are also provided to reduce a general real or complex matrix pair to the required upper trapezoidal form.

	Reduce to trapezoidal form	Generalized SVD of trapezoidal form
real matrices	F08VEF (DGGSVP)	F08YEF (DTGSJA)
complex matrices	F08VSF (ZGGSVP)	F08YSF (ZTGSJA)

Routines are provided to reduce a real symmetric or complex Hermitian matrix

A

to real tridiagonal form

T

by an orthogonal similarity transformation

A = Q T Q^{T}

(or by a unitary transformation

A = Q T Q^{H}

if

A

is complex). Different routines allow a full matrix

A

to be stored conventionally (see Section 3.3.1 in the F07 Chapter Introduction) or in packed storage (see Section 3.3.2 in the F07 Chapter Introduction); or a band matrix to use band storage (see Section 3.3.4 in the F07 Chapter Introduction).

The routines for reducing full matrices do not form the matrix

Q

explicitly; additional routines are provided to generate

Q

, or to apply it to another matrix, as with the routines for orthogonal factorizations. Explicit generation of

Q

is required before using the

Q R

algorithm to find all the eigenvectors of

A

; application of

Q

to another matrix is required after eigenvectors of

T

have been found by inverse iteration, in order to transform them to eigenvectors of

A

.

The routines for reducing band matrices have an option to generate

Q

if required.

	Reduce to tridiagonal form	Generate matrix $Q$	Apply matrix $Q$
real symmetric matrices	F08FEF (DSYTRD)	F08FFF (DORGTR)	F08FGF (DORMTR)
real symmetric matrices (packed storage)	F08GEF (DSPTRD)	F08GFF (DOPGTR)	F08GGF (DOPMTR)
real symmetric band matrices	F08HEF (DSBTRD)
complex Hermitian matrices	F08FSF (ZHETRD)	F08FTF (ZUNGTR)	F08FUF (ZUNMTR)
complex Hermitian matrices (packed storage)	F08GSF (ZHPTRD)	F08GTF (ZUPGTR)	F08GUF (ZUPMTR)
complex Hermitian band matrices	F08HSF (ZHBTRD)

Given the eigenvalues, F08FLF (DDISNA) is provided to compute the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix.

A variety of routines are provided to compute eigenvalues and eigenvectors of the real symmetric tridiagonal matrix

T

, some computing all eigenvalues and eigenvectors, some computing selected eigenvalues and eigenvectors. The same routines can be used to compute eigenvalues and eigenvectors of a real symmetric or complex Hermitian matrix which has been reduced to tridiagonal form.

Eigenvalues and eigenvectors of real symmetric tridiagonal matrices:

The original (non-reduced) matrix is Real or Complex Hermitian

all eigenvalues (root-free $Q R$ algorithm)	F08JFF
all eigenvalues (root-free $Q R$ algorithm called by divide-and-conquer)	F08JCF or F08JHF
all eigenvalues (RRR)	F08JLF
selected eigenvalues (bisection)	F08JJF

The original (non-reduced) matrix is Real

all eigenvalues and eigenvectors ( $Q R$ algorithm)	F08JEF
all eigenvalues and eigenvectors (divide-and-conquer)	F08JCF or F08JHF
all eigenvalues and eigenvectors (RRR)	F08JLF
all eigenvalues and eigenvectors (positive definite case)	F08JGF
selected eigenvectors (inverse iteration)	F08JKF

The original (non-reduced) matrix is Complex Hermitian

all eigenvalues and eigenvectors ( $Q R$ algorithm)	F08JSF
all eigenvalues and eigenvectors (divide and conquer)	F08JVF
all eigenvalues and eigenvectors (RRR)	F08JYF
all eigenvalues and eigenvectors (positive definite case)	F08JUF
selected eigenvectors (inverse iteration)	F08JXF

The following sequences of calls may be used to compute various combinations of eigenvalues and eigenvectors, as described in Section 2.7.

Sequences for computing eigenvalues and eigenvectors

Real Symmetric matrix (standard storage)

all eigenvalues and eigenvectors (using divide-and-conquer)	F08FCF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08FEF, F08FFF*, F08JEF
all eigenvalues and eigenvectors (RRR)	F08FEF, F08FGF, F08JLF
selected eigenvalues and eigenvectors (bisection and inverse iteration)	F08FEF, F08FGF, F08JJF, F08JKF*

Real Symmetric matrix (packed storage)

all eigenvalues and eigenvectors (using divide-and-conquer)	F08GCF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08GEF, F08GFF*, F08JEF
all eigenvalues and eigenvectors (RRR)	F08GEF, F08GGF, F08JLF
selected eigenvalues and eigenvectors (bisection and inverse iteration)	F08GEF, F08GGF, F08JJF, F08JKF*

Real Symmetric banded matrix

all eigenvalues and eigenvectors (using divide-and-conquer)	F08HCF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08HEF*, F08JEF

Complex Hermitian matrix (standard storage)

all eigenvalues and eigenvectors (using divide-and-conquer)	F08FQF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08FSF, F08FTF*, F08JSF
all eigenvalues and eigenvectors (RRR)	F08FSF, F08FUF, F08JYF
selected eigenvalues and eigenvectors (bisection and inverse iteration)	F08FSF, F08FUF, F08JJF, F08JXF*

Complex Hermitian matrix (packed storage)

all eigenvalues and eigenvectors (using divide-and-conquer)	F08GQF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08GSF, F08GTF*, F08JSF
all eigenvalues and eigenvectors (RRR)	F08GSF, F08GUF and F08JYF
selected eigenvalues and eigenvectors (bisection and inverse iteration)	F08GSF, F08GUF, F08JJF, F08JXF*

Complex Hermitian banded matrix

all eigenvalues and eigenvectors (using divide-and-conquer)	F08HQF
all eigenvalues and eigenvectors (using $Q R$ algorithm)	F08HSF*, F08JSF

Routines are provided for reducing each of the problems

A x = λ B x

,

A B x = λ x

or

B A x = λ x

to an equivalent standard eigenvalue problem

C y = λ y

. Different routines allow the matrices to be stored either conventionally or in packed storage. The positive definite matrix

B

must first be factorized using a routine from Chapter F07. There is also a routine which reduces the problem

A x = λ B x

where

A

and

B

are banded, to an equivalent banded standard eigenvalue problem; this uses a split Cholesky factorization for which a routine in Chapter F08 is provided.

	Reduce to standard problem	Reduce to standard problem (packed storage)	Reduce to standard problem (band matrices)
real symmetric matrices	F08SEF (DSYGST)	F08TEF (DSPGST)	F08UEF (DSBGST)
complex Hermitian matrices	F08SSF (ZHEGST)	F08TSF (ZHPGST)	F08USF (ZHBGST)

The equivalent standard problem can then be solved using the routines discussed in Section 3.1.2.5. For example, to compute all the eigenvalues, the following routines must be called:

real symmetric-definite problem	F07FDF, F08SEF*, F08FEF, F08JFF
real symmetric-definite problem, packed storage	F07GDF, F08TEF*, F08GEF, F08JFF
real symmetric-definite banded problem	F08UFF, F08UEF, F08HEF, F08JFF
complex Hermitian-definite problem	F07FRF, F08SSF*, F08FSF, F08JFF
complex Hermitian-definite problem, packed storage	F07GRF, F08TSF*, F08GSF, F08JFF
complex Hermitian-definite banded problem	F08UTF, F08USF, F08HSF, F08JFF

If eigenvectors are computed, the eigenvectors of the equivalent standard problem must be transformed back to those of the original generalized problem, as indicated in Section 2.8; routines from Chapter F06 may be used for this.

Routines are provided to reduce a general real or complex matrix

A

to upper Hessenberg form

H

by an orthogonal similarity transformation

A = Q H Q^{T}

(or by a unitary transformation

A = Q H Q^{H}

if

A

is complex).

These routines do not form the matrix

Q

explicitly; additional routines are provided to generate

Q

, or to apply it to another matrix, as with the routines for orthogonal factorizations. Explicit generation of

Q

is required before using the

Q R

algorithm on

H

to compute the Schur vectors; application of

Q

to another matrix is needed after eigenvectors of

H

have been computed by inverse iteration, in order to transform them to eigenvectors of

A

.

Routines are also provided to balance the matrix before reducing it to Hessenberg form, as described in Section 2.14.6. Companion routines are required to transform Schur vectors or eigenvectors of the balanced matrix to those of the original matrix.

	Reduce to Hessenberg form	Generate matrix $Q$	Apply matrix $Q$	Balance	Backtransform vectors after balancing
real matrices	F08NEF (DGEHRD)	F08NFF (DORGHR)	F08NGF (DORMHR)	F08NHF (DGEBAL)	F08NJF (DGEBAK)
complex matrices	F08NSF (ZGEHRD)	F08NTF (ZUNGHR)	F08NUF (ZUNMHR)	F08NVF (ZGEBAL)	F08NWF (ZGEBAK)

Routines are provided to compute the eigenvalues and all or part of the Schur factorization of an upper Hessenberg matrix. Eigenvectors may be computed either from the upper Hessenberg form by inverse iteration, or from the Schur form by back-substitution; these approaches are equally satisfactory for computing individual eigenvectors, but the latter may provide a more accurate basis for a subspace spanned by several eigenvectors.

Additional routines estimate the sensitivities of computed eigenvalues and eigenvectors, as discussed in Section 2.14.5.

	Eigenvalues and Schur factorization ( $QR$ algorithm)	Eigenvectors from Hessenberg form (inverse iteration)	Eigenvectors from Schur factorization	Sensitivities of eigenvalues and eigenvectors
real matrices	F08PEF (DHSEQR)	F08PKF (DHSEIN)	F08QKF (DTREVC)	F08QLF (DTRSNA)
complex matrices	F08PSF (ZHSEQR)	F08PXF (ZHSEIN)	F08QXF (ZTREVC)	F08QYF (ZTRSNA)

Finally routines are provided for reordering the Schur factorization, so that eigenvalues appear in any desired order on the diagonal of the Schur form. The routines F08QFF (DTREXC) and F08QTF (ZTREXC) simply swap two diagonal elements or blocks, and may need to be called repeatedly to achieve a desired order. The routines F08QGF (DTRSEN) and F08QUF (ZTRSEN) perform the whole reordering process for the important special case where a specified cluster of eigenvalues is to appear at the top of the Schur form; if the Schur vectors are reordered at the same time, they yield an orthonormal basis for the invariant subspace corresponding to the specified cluster of eigenvalues. These routines can also compute the sensitivities of the cluster of eigenvalues and the invariant subspace.

	Reorder Schur factorization	Reorder Schur factorization, find basis for invariant subspace and estimate sensitivities
real matrices	F08QFF (DTREXC)	F08QGF (DTRSEN)
complex matrices	F08QTF (ZTREXC)	F08QUF (ZTRSEN)

The following sequences of calls may be used to compute various combinations of eigenvalues, Schur vectors and eigenvectors, as described in Section 2.11:

real matrix, all eigenvalues and Schur factorization	F08NEF, F08NFF*, F08PEF
real matrix, all eigenvalues and selected eigenvectors	F08NEF, F08NGF, F08PEF, F08PKF
real matrix, all eigenvalues and eigenvectors (with balancing)	F08NHF*, F08NEF, F08NFF, F08NJF, F08PEF, F08PKF
complex matrix, all eigenvalues and Schur factorization	F08NSF, F08NTF*, F08PSF
complex matrix, all eigenvalues and selected eigenvectors	F08NSF, F08NUF, F08PSF, F08PXF*
complex matrix, all eigenvalues and eigenvectors (with balancing)	F08NVF*, F08NSF, F08NTF, F08NWF, F08PSF, F08PXF

Routines are provided to reduce a real or complex matrix pair

(A_{1}, R_{1})

, where

A_{1}

is general and

R_{1}

is upper triangular, to generalized upper Hessenberg form by orthogonal transformations

A_{1} = Q_{1} H Z_{1}^{T}

,

R_{1} = Q_{1} R Z_{1}^{T}

, (or by unitary transformations

A_{1} = Q_{1} H Z_{1}^{H}

,

R = Q_{1} R_{1} Z_{1}^{H}

, in the complex case). These routines can optionally return

Q_{1}

and/or

Z_{1}

. Note that to transform a general matrix pair

(A, B)

to the form

(A_{1}, R_{1})

a

Q R

factorization of

B

(

B = \tilde{Q} R_{1}

) should first be performed and the matrix

A_{1}

obtained as

A_{1} = {\tilde{Q}}^{T} A

(see Section 3.1.2.1 above).

Routines are also provided to balance a general matrix pair before reducing it to generalized Hessenberg form, as described in Section 2.14.8. Companion routines are provided to transform vectors of the balanced pair to those of the original matrix pair.

	Reduce to generalized Hessenberg form	Balance	Backtransform vectors after balancing
real matrices	F08WEF (DGGHRD)	F08WHF (DGGBAL)	F08WJF (DGGBAK)
complex matrices	F08WSF (ZGGHRD)	F08WVF (ZGGBAL)	F08WWF (ZGGBAK)

Routines are provided to compute the eigenvalues (as the pairs

(α, β)

) and all or part of the generalized Schur factorization of a generalized upper Hessenberg matrix pair. Eigenvectors may be computed from the generalized Schur form by back-substitution.

Additional routines estimate the sensitivities of computed eigenvalues and eigenvectors.

	Eigenvalues and generalized Schur factorization ( $QZ$ algorithm)	Eigenvectors from generalized Schur factorization	Sensitivities of eigenvalues and eigenvectors
real matrices	F08XEF (DHGEQZ)	F08YKF (DTGEVC)	F08YLF (DTGSNA)
complex matrices	F08XSF (ZHGEQZ)	F08YXF (ZTGEVC)	F08YYF (ZTGSNA)

Finally, routines are provided for reordering the generalized Schur factorization so that eigenvalues appear in any desired order on the diagonal of the generalized Schur form. F08YFF (DTGEXC) and F08YTF (ZTGEXC) simply swap two diagonal elements or blocks, and may need to be called repeatedly to achieve a desired order. F08YGF (DTGSEN) and F08YUF (ZTGSEN) perform the whole reordering process for the important special case where a specified cluster of eigenvalues is to appear at the top of the generalized Schur form; if the Schur vectors are reordered at the same time, they yield an orthonormal basis for the deflating subspace corresponding to the specified cluster of eigenvalues. These routines can also compute the sensitivities of the cluster of eigenvalues and the deflating subspace.

	Reorder generalized Schur factorization	Reorder generalized Schur factorization, find basis for deflating subspace and estimate sensitivites
real matrices	F08YFF (DTGEXC)	F08YGF (DTGSEN)
complex matrices	F08YTF (ZTGEXC)	F08YUF (ZTGSEN)

The following sequences of calls may be used to compute various combinations of eigenvalues, generalized Schur vectors and eigenvectors

real matrix pair, all eigenvalues (with balancing)	F08AEF, F08AGF, F08WEF, F08WHF, F08XEF*
real matrix pair, all eigenvalues and generalized Schur factorization	F08AEF, F08AFF, F08AGF, F08WEF, F08XEF
real matrix pair, all eigenvalues and eigenvectors (with balancing)	F06QFF, F06QHF, F08AEF, F08AFF, F08AGF, F08WEF, F08WHF, F08XEF, F08YKF*, F08WJF
complex matrix pair, all eigenvalues (with balancing)	F08ASF, F08AUF, F08WSF, F08WVF, F08XSF*
complex matrix pair, all eigenvalues and generalized Schur factorization	F08ASF, F08ATF, F08AUF, F08WSF, F08XSF
complex matrix pair, all eigenvalues and eigenvectors (with balancing)	F06TFF, F06THF, F08ASF, F08ATF, F08AUF, F08WSF, F08WVF, F08XSF, F08YXF*, F08WWF

Routines are provided to solve the real or complex Sylvester equation

A X \pm X B = C

, where

A

and

B

are upper quasi-triangular if real, or upper triangular if complex. To solve the general form of the Sylvester equation in which

A

and

B

are general square matrices,

A

and

B

must be reduced to upper (quasi-) triangular form by the Schur factorization, using routines described in Section 3.1.2.7. For more details, see the documents for the routines listed below.

	Solve the Sylvester equation
real matrices	F08QHF (DTRSYL)
complex matrices	F08QVF (ZTRSYL)

Routines are also provided to solve the real or complex generalized Sylvester equations

A R - L B = C,   ​ D R - L E = F,

where the pairs

(A, D)

and

(B, E)

are in generalized Schur form. To solve the general form of the generalized Sylvester equation in which

(A, D)

and

(B, E)

are general matrix pairs,

(A, D)

and

(B, E)

must first be reduced to generalized Schur form.

	Solve the generalized Sylvester equation
real matrices	F08YHF (DTGSYL)
complex matrices	F08YVF (ZTGSYL)

As well as the NAG routine name (beginning F08), the tables in Section 3.1 show the LAPACK routine names in double precision.

The routines may be called either by their NAG or LAPACK names. When using the NAG Library, the double precision form of the LAPACK name must be used (beginning with D- or Z-).

References to Chapter F08 routines in the manual normally include the LAPACK double precision names, for example F08AEF (DGEQRF). The LAPACK routine names follow a simple scheme (which is similar to that used for the BLAS in Chapter F06). Each name has the structure XYYZZZ, where the components have the following meanings:

–

the initial letter X indicates the data type (real or complex) and precision:

S	– real, single precision (in Fortran 77, REAL)
D	– real, double precision (in Fortran 77, DOUBLE PRECISION)
C	– complex, single precision (in Fortran 77, COMPLEX)
Z	– complex, double precision (in Fortran 77, COMPLEX*16 or DOUBLE COMPLEX)

–

the second and third letters YY indicate the type of the matrix

A

or matrix pair

(A, B)

(and in some cases the storage scheme):

BD	– bidiagonal
DI	– diagonal
GB	– general band
GE	– general
GG	– general pair ( $B$ may be triangular)
HG	– generalized upper Hessenberg
HS	– upper Hessenberg
OP	– (real) orthogonal (packed storage)
UP	– (complex) unitary (packed storage)
OR	– (real) orthogonal
UN	– (complex) unitary
PT	– symmetric or Hermitian positive definite tridiagonal
SB	– (real) symmetric band
HB	– (complex) Hermitian band
SP	– symmetric (packed storage)
HP	– Hermitian (packed storage)
ST	– (real) symmetric tridiagonal
SY	– symmetric
HE	– Hermitian
TG	– triangular pair (one may be quasi-triangular)
TR	– triangular (or quasi-triangular)

–

the last three letters ZZZ indicate the computation performed. For example, QRF is a

Q R

factorization.

Thus the routine DGEQRF performs a

Q R

factorization of a real general matrix; the corresponding routine for a complex general matrix is ZGEQRF.

In this chapter the following storage schemes are used for matrices:

– conventional storage in a two-dimensional array;
– packed storage for symmetric or Hermitian matrices;
– packed storage for orthogonal or unitary matrices;
– band storage for general, symmetric or Hermitian band matrices;
– storage of bidiagonal, symmetric or Hermitian tridiagonal matrices in two one-dimensional arrays.

These storage schemes are compatible with those used in Chapters F06 and F07, but different schemes for packed, band and tridiagonal storage are used in a few older routines in Chapters F01, F02, F03 and F04.

Please see Section 3.3.1 in the F07 Chapter Introduction for full details.

Please see Section 3.3.2 in the F07 Chapter Introduction for full details.

Please see Section 3.3.4 in the F07 Chapter Introduction for full details.

A symmetric tridiagonal or bidiagonal matrix is stored in two one-dimensional arrays, one of length

n

containing the diagonal elements, and one of length

n - 1

containing the off-diagonal elements. (Older routines in Chapter F02 store the off-diagonal elements in elements

2 : n

of a vector of length

n

.)

Please see Section 3.3.6 in the F07 Chapter Introduction for full details.

A real orthogonal or complex unitary matrix (usually denoted

Q

) is often represented in the NAG Library as a product of elementary reflectors – also referred to as elementary Householder matrices (usually denoted

H_{i}

). For example,

Q = H_{1} H_{2} \dots H_{k} .

You need not be aware of the details, because routines are provided to work with this representation, either to generate all or part of

Q

explicitly, or to multiply a given matrix by

Q

or

Q^{T}

(

Q^{H}

in the complex case) without forming

Q

explicitly.

Nevertheless, the following further details may occasionally be useful.

An elementary reflector (or elementary Householder matrix)

H

of order

n

is a unitary matrix of the form

H = I - τ v v^{H}

(4)

where

τ

is a scalar, and

v

is an

n

-element vector, with

{|τ|}^{2} {‖v‖}_{2}^{2} = 2 \times Re (τ)

;

v

is often referred to as the Householder vector. Often

v

has several leading or trailing zero elements, but for the purpose of this discussion assume that

H

has no such special structure.

There is some redundancy in the representation

(4)

, which can be removed in various ways. The representation used in Chapter F08 and in LAPACK (which differs from those used in some of the routines in Chapters F01, F02, F04 and F06) sets

v_{1} = 1

; hence

v_{1}

need not be stored. In real arithmetic,

1 \leq τ \leq 2

, except that

τ = 0

implies

H = I

.

In complex arithmetic,

τ

may be complex, and satisfies

1 \leq Re (τ) \leq 2

and

|τ - 1| \leq 1

. Thus a complex

H

is not Hermitian (as it is in other representations), but it is unitary, which is the important property. The advantage of allowing

τ

to be complex is that, given an arbitrary complex vector

x, H

can be computed so that

H^{H} x = β {(1, 0, \dots, 0)}^{T}

with real

β

. This is useful, for example, when reducing a complex Hermitian matrix to real symmetric tridiagonal form, or a complex rectangular matrix to real bidiagonal form.

Most routines in this chapter have one or more option parameters, of type CHARACTER. The descriptions in Section 5 of the routine documents refer only to upper case values (for example

UPLO ='U'

or

UPLO ='L'

); however in every case, the corresponding lower case characters may be supplied (with the same meaning). Any other value is illegal.

A longer character string can be passed as the actual parameter, making the calling program more readable, but only the first character is significant. (This is a feature of Fortran 77.) For example:

CALL SSYTRD ('Upper',...)

It is permissible for the problem dimensions (for example, M or N) to be passed as zero, in which case the computation (or part of it) is skipped. Negative dimensions are regarded as an error.

A number of routines implementing block algorithms require workspace sufficient to hold one block of rows or columns of the matrix if they are to achieve optimum levels of performance – for example, workspace of size

n \times n b

, where

n b

is the optimal block size. In such cases, the actual declared length of the work array must be passed as a separate argument LWORK, which immediately follows WORK in the argument-list.

The routine will still perform correctly when less workspace is provided: it simply uses the largest block size allowed by the amount of workspace supplied, as long as this is likely to give better performance than the unblocked algorithm. On exit,

WORK (1)

contains the minimum value of LWORK which would allow the routine to use the optimal block size; this value of LWORK can be used for subsequent runs.

If LWORK indicates that there is insufficient workspace to perform the unblocked algorithm, this is regarded as an illegal value of LWORK, and is treated like any other illegal parameter value (see Section 3.4.4).

If you are in doubt how much workspace to supply and are concerned to achieve optimal performance, supply a generous amount (assume a block size of

64

, say), and then examine the value of

WORK (1)

on exit.

Routines in this chapter do not use the usual NAG Library error-handling mechanism, involving the parameter IFAIL. Instead they have a diagnostic parameter INFO. (Thus they preserve complete compatibility with the LAPACK specification.)

Whereas IFAIL is an Input/Output parameter and must be set before calling a routine, INFO is purely an Output parameter and need not be set before entry.

INFO indicates the success or failure of the computation, as follows:

$INFO = 0$ : successful termination;
$INFO > 0$ : failure in the course of computation, control returned to the calling program.

If the routine document specifies that the routine may terminate with

INFO > 0

, then it is essential to test INFO on exit from the routine. (This corresponds to a soft failure in terms of the usual NAG error-handling terminology.) No error message is output.

All routines check that input parameters such as N or LDA or option parameters of type CHARACTER have permitted values. If an illegal value of the

i

th parameter is detected, INFO is set to

- i

, a message is output, and execution of the program is terminated. (This corresponds to a hard failure in the usual NAG terminology.) In some implementations, especially when linking to vendor versions of LAPACK, execution of the program may continue, in which case, it is essential to test INFO on exit from the routine.

The following decision trees are principally for the computation (general purpose) routines. See Section 3.1.1.1 for tables of the driver (black box) routines.

Note: the routines for band matrices only handle the problem

A x = λ B x

; the other routines handle all three types of problems (

A x = λ B x

,

A B x = λ x

or

B A x = λ x

) except that, if the problem is

B A x = λ x

and eigenvectors are required, F06PHF must be used instead of F06PLF and F06YFF instead of F06YJF.

F08PAZ	nagf_lapack_dgees_dummy_select See the description of the argument SELECT in F08PAF (DGEES) and F08PBF (DGEESX).
F08PNZ	nagf_lapack_zgees_dummy_select See the description of the argument SELECT in F08PNF (ZGEES) and F08PPF (ZGEESX).
F08XAZ	nagf_lapack_dgges_dummy_selctg See the description of the argument SELCTG in F08XAF (DGGES) and F08XBF (DGGESX).
F08XNZ	nagf_lapack_zgges_dummy_selctg See the description of the argument SELCTG in F08XNF (ZGGES) and F08XPF (ZGGESX).

None.

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia

Arioli M, Duff I S and de Rijk P P M (1989) On the augmented system approach to sparse least-squares problems Numer. Math. 55 667–684

Demmel J W and Kahan W (1990) Accurate singular values of bidiagonal matrices SIAM J. Sci. Statist. Comput. 11 873–912

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Moler C B and Stewart G W (1973) An algorithm for generalized matrix eigenproblems SIAM J. Numer. Anal. 10 241–256

Parlett B N (1998) The Symmetric Eigenvalue Problem SIAM, Philadelphia

Stewart G W and Sun J-G (1990) Matrix Perturbation Theory Academic Press, London

Ward R C (1981) Balancing the generalized eigenvalue problem SIAM J. Sci. Stat. Comp. 2 141–152

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

Operation	real	complex
solve LLS using $Q R$ or $L Q$ factorization solve LLS using complete orthogonal factorization solve LLS using SVD solve LLS using divide-and-conquer SVD	F08AAF (DGELS) F08BAF (DGELSY) F08KAF (DGELSS) F08KCF (DGELSD)	F08ANF (ZGELS) F08BNF (ZGELSY) F08KNF (ZGELSS) F08KQF (ZGELSD)

Operation	real	complex
solve LSE problem using GRQ solve GLM problem using GQR	F08ZAF (DGGLSE) F08ZBF (DGGGLM)	F08ZNF (ZGGLSE) F08ZPF (ZGGGLM)

Function and storage scheme	real	complex
simple driver divide-and-conquer driver expert driver RRR driver	F08FAF (DSYEV) F08FCF (DSYEVD) F08FBF (DSYEVX) F08FDF (DSYEVR)	F08FNF (ZHEEV) F08FQF (ZHEEVD) F08FPF (ZHEEVX) F08FRF (ZHEEVR)
packed storage simple driver divide-and-conquer driver expert driver	F08GAF (DSPEV) F08GCF (DSPEVD) F08GBF (DSPEVX)	F08GNF (ZHPEV) F08GQF (ZHPEVD) F08GPF (ZHPEVX)
band matrix simple driver divide-and-conquer driver expert driver	F08HAF (DSBEV) F08HCF (DSBEVD) F08HBF (DSBEVX)	F08HNF (ZHBEV) F08HQF (ZHBEVD) F08HPF (ZHBEVX)
tridiagonal matrix simple driver divide-and-conquer driver expert driver RRR driver	F08JAF (DSTEV) F08JCF (DSTEVD) F08JBF (DSTEVX) F08JDF (DSTEVR)

Function and storage scheme	real	complex
simple driver for Schur factorization expert driver for Schur factorization simple driver for eigenvalues/vectors expert driver for eigenvalues/vectors	F08PAF (DGEES) F08PBF (DGEESX) F08NAF (DGEEV) F08NBF (DGEEVX)	F08PNF (ZGEES) F08PPF (ZGEESX) F08NNF (ZGEEV) F08NPF (ZGEEVX)

Function and storage scheme	real	complex
simple driver divide-and-conquer driver simple driver for one-sided Jacobi SVD expert driver for one-sided Jacobi SVD	F08KBF (DGESVD) F08KDF (DGESDD) F08KJF (DGESVJ) F08KHF (DGEJSV)	F08KPF (ZGESVD) F08KRF (ZGESDD)

Are eigenvalues only required?	_ yes	Are all the eigenvalues required?	_ yes	Is $A$ tridiagonal?	_ yes	F08JCF or F08JFF
\|		\|		no \|
\|		\|		Is $A$ band matrix?	_ yes	(F08HEF and F08JFF) or F08HCF
\|		\|		no \|
\|		\|		Is one triangle of $A$ stored as a linear array?	_ yes	(F08GEF and F08JFF) or F08GCF
\|		\|		no \|
\|		\|		(F08FEF and F08JFF) or F08FAF or F08FCF
\|		no \|
\|		Is $A$ tridiagonal?	_ yes	F08JJF
\|		no \|
\|		Is $A$ a band matrix?	_ yes	F08HEF and F08JJF
\|		no \|
\|		Is one triangle of $A$ stored as a linear array?	_ yes	F08GEF and F08JJF
\|		no \|
\|		(F08FEF and F08JJF) or F08FBF
no \|
Are all eigenvalues and eigenvectors required?	_ yes	Is $A$ tridiagonal?	_ yes	F08JEF, F08JCF, F08JHF or F08JLF
\|		no \|
\|		Is $A$ a band matrix?	_ yes	(F08HEF and F08JEF) or F08HCF
\|		no \|
\|		Is one triangle of $A$ stored as a linear array?	_ yes	(F08GEF, F08GFF and F08JEF) or F08GCF
\|		no \|
\|		(F08FEF, F08FFF and F08JEF) or F08FAF or F08FCF
no \|
Is $A$ tridiagonal?	_ yes	F08JJF and F08JKF
no \|
Is one triangle of $A$ stored as a linear array?	_ yes	F08GEF, F08JJF, F08JKF and F08GGF
no \|
(F08FEF, F08JJF, F08JKF and F08FGF) or F08FBF

Are eigenvalues required?	_ yes	Is $A$ an upper Hessenberg matrix?	_ yes	F08PEF
\|		no \|
\|		F08NAF or F08NBF or (F08NHF, F08NEF and F08PEF)
no \|
Is the Schur factorization of $A$ required?	_ yes	Is $A$ an upper Hessenberg matrix?	_ yes	F08PEF
\|		no \|
\|		F08NBF or (F08NEF, F08NFF, F08PEF or F08NJF)
no \|
Are all eigenvectors required?	_ yes	Is $A$ an upper Hessenberg matrix?	_ yes	F08PEF or F08QKF
\|		no \|
\|		F08NAF or F08NBF or (F08NHF, F08NEF, F08NFF, F08PEF, F08QKF or F08NJF)
no \|
Is $A$ an upper Hessenberg matrix?	_ yes	F08PEF or F08PKF
no \|
F08NHF, F08NEF, F08PEF, F08PKF, F08NGF or F08NJF

Are eigenvalues only required?	_ yes	Are $A$ and $B$ in generalized upper Hessenberg form?	_ yes	F08XEF
\|		no \|
\|		F08WHF, F08AEF, F08AGF, F08WEF and F08XEF
no \|
Is the generalized Schur factorization of $A$ and $B$ required?	_ yes	Are $A$ and $B$ in generalized upper Hessenberg form?	_ yes	F08XEF
\|		no \|
\|		F08AEF, F08AGF, F06QHF, F06QFF, F08AFF, F08WEF, F08XEF and F08YKF
no \|
Are $A$ and $B$ in generalized upper Hessenberg form?	_ yes	F08XEF and F08YKF
no \|
F08WHF, F08AEF, F08AGF, F06QHF, F06QFF, F08AFF, F08WEF, F08XEF, F08YKF and F08WJF

Are eigenvalues only required?	_ yes	Are all eigenvalues required?	_ yes	Are $A$ and $B$ stored with one triangle as a linear array?	_ yes	F07GRF, F08TSF, F08GSF and F08JFF
\|		\|		no \|
\|		\|		F07FRF, F08SSF, F08FSF and F08JFF
\|		no \|
\|		Are $A$ and $B$ stored with one triangle as a linear array?	_ yes	F07GRF, F08TSF, F08GSF and F08JJF
\|		no \|
\|		F07FRF, F08SSF, F08GSF and F08JJF
no \|
Are all eigenvalues and eigenvectors required?	_ yes	Are $A$ and $B$ stored with one triangle as a linear array?	_ yes	F07GRF, F08TSF, F08GSF, F08GTF and F06PSF
\|		no \|
\|		F07FRF, F08SSF, F08FSF, F08FTF, F08JSF and F06ZJF
no \|
Are $A$ and $B$ stored with one triangle as a linear array?	_ yes	F07GRF, F08TSF, F08GSF, F08JJF, F08JXF, F08GUF and F06SLF
no \|
F07FRF, F08SSF, F08FSF, F08JJF, F08JXF, F08FUF and F06ZJF

Is $A$ a complex matrix?	_ yes	Is $A$ banded?	_ yes	F08LSF and F08MSF
\|		no \|
\|		Are singular values only required?	_ yes	F08KSF and F08MSF
\|		no \|
\|		F08KSF, F08KTF and F08MSF
no \|
Is $A$ bidiagonal?	_ yes	F08MEF
no \|
Is $A$ banded?	_ yes	F08LEF and F08MEF
no \|
Are singular values only required?	_ yes	F08KEF and F08MEF
no \|
F08KEF, F08KFF and F08MEF

NAG Library Chapter IntroductionF08 – Least Squares and Eigenvalue Problems (LAPACK)

+− Contents

1 Scope of the Chapter

2 Background to the Problems

2.1 Linear Least Squares Problems

2.2 Orthogonal Factorizations and Least Squares Problems

2.2.1 QR factorization

2.2.2 LQ factorization

2.2.3 QR factorization with column pivoting

2.2.4 Complete orthogonal factorization

2.2.5 Other factorizations

2.3 The Singular Value Decomposition

2.4 The Singular Value Decomposition and Least Squares Problems

2.5 Generalized Linear Least Squares Problems

2.6 Generalized Orthogonal Factorization and Generalized Linear Least Squares Problems

2.6.1 Generalized QR Factorization

2.6.2 Generalized RQ Factorization

2.6.3 Generalized Singular Value Decomposition (GSVD)

2.7 Symmetric Eigenvalue Problems

2.8 Generalized Symmetric-definite Eigenvalue Problems

2.9 Packed Storage for Symmetric Matrices

2.10 Band Matrices

2.11 Nonsymmetric Eigenvalue Problems

2.12 Generalized Nonsymmetric Eigenvalue Problem

2.13 The Sylvester Equation and the Generalized Sylvester Equation

2.14 Error and Perturbation Bounds and Condition Numbers

2.14.1 Least squares problems

2.14.2 The singular value decomposition

2.14.3 The symmetric eigenproblem

2.14.4 The generalized symmetric-definite eigenproblem

2.14.5 The nonsymmetric eigenproblem

2.14.6 Balancing and condition for the nonsymmetric eigenproblem

2.14.7 The generalized nonsymmetric eigenvalue problem

2.14.8 Balancing the generalized eigenvalue problem

2.14.9 Other problems

2.15 Block Partitioned Algorithms

3 Recommendations on Choice and Use of Available Routines

3.1 Available Routines

3.1.1 Driver routines

3.1.1.1 Linear least squares problems (LLS)

3.1.1.2 Generalized linear least squares problems (LSE and GLM)

3.1.1.3 Symmetric eigenvalue problems (SEP)

3.1.1.4 Nonsymmetric eigenvalue problem (NEP)

3.1.1.5 Singular value decomposition (SVD)

3.1.1.6 Generalized symmetric definite eigenvalue problems (GSEP)

3.1.1.7 Generalized nonsymmetric eigenvalue problem (GNEP)

3.1.1.8 Generalized singular value decomposition (GSVD)

3.1.2 Computational routines

3.1.2.1 Orthogonal factorizations

3.1.2.2 Generalized orthogonal factorizations

3.1.2.3 Singular value problems

3.1.2.4 Generalized singular value decomposition

3.1.2.5 Symmetric eigenvalue problems

3.1.2.6 Generalized symmetric-definite eigenvalue problems

3.1.2.7 Nonsymmetric eigenvalue problems

3.1.2.8 Generalized nonsymmetric eigenvalue problems

3.1.2.9 The Sylvester equation and the generalized Sylvester equation

3.2 NAG Names and LAPACK Names

3.3 Matrix Storage Schemes

3.3.1 Conventional storage

3.3.2 Packed storage

3.3.3 Band storage

3.3.4 Tridiagonal and bidiagonal matrices

3.3.5 Real diagonal elements of complex matrices

3.3.6 Representation of orthogonal or unitary matrices

3.4 Parameter Conventions

3.4.1 Option parameters

3.4.2 Problem dimensions

3.4.3 Length of work arrays

3.4.4 Error-handling and the diagnostic parameter INFO

4 Decision Trees

4.1 General Purpose Routines (eigenvalues and eigenvectors)

Tree 1: Real Symmetric Eigenvalue Problems

Tree 2: Real Generalized Symmetric-definite Eigenvalue Problems

Tree 3: Real Nonsymmetric Eigenvalue Problems

Tree 4: Real Generalized Nonsymmetric Eigenvalue Problems

Tree 5: Complex Hermitian Eigenvalue Problems

Tree 6: Complex Generalized Hermitian-definite Eigenvalue Problems

Tree 7: Complex non-Hermitian Eigenvalue Problems

Tree 8: Complex Generalized non-Hermitian Eigenvalue Problems

NAG Library Chapter Introduction

F08 – Least Squares and Eigenvalue Problems (LAPACK)

2.2.1 $Q R$ factorization

2.2.2 $L Q$ factorization

2.2.3 $Q R$ factorization with column pivoting

2.6.1 Generalized $Q R$ Factorization

2.6.2 Generalized $R Q$ Factorization