NAG Library Chapter Introduction
f08 – Least Squares and Eigenvalue Problems (LAPACK)
1 Scope of the Chapter
This chapter provides functions for the solution of linear least squares problems, eigenvalue problems and singular value problems, as well as associated computations. It provides functions 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 linear least squares problems
- – solution of generalized symmetric-definite eigenvalue problems
- – solution of generalized nonsymmetric eigenvalue problems
- – solution of generalized singular value 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
Functions 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 functions in
Chapters f04 or
f08.
Chapters f04 and
f08 contain
Black Box (or
driver) functions which enable standard linear least squares problems to be solved by a call to a single function.
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 functions in
Chapters f02 or
f08.
Chapters f02 and
f08 contain
Black Box (or
driver) functions which enable standard types of problem to be solved by a call to a single function. Often functions in
Chapter f02 call
Chapter f08 functions to perform the necessary computational tasks.
The functions 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 functions in
Chapter f08.
The functions 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.
2 Background to the Problems
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 (2012).
2.1 Linear Least Squares Problems
The
linear least squares problem is
where
is an
by
matrix,
is a given
element vector and
is an
-element solution vector.
In the most usual case
and
, so that
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 and , there are an infinite number of solutions which exactly satisfy . In this case it is often useful to find the unique solution which minimizes , 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 – in other words, may be rank-deficient – we seek the minimum norm least squares solution which minimizes both and .
This chapter (
Chapter f08) contains driver functions to solve these problems with a single call, as well as computational functions that can be combined with functions in
Chapter f07 to solve these linear least squares problems.
The next two sections discuss the factorizations that can be used in the solution of linear least squares problems.
2.2 Orthogonal Factorizations and Least Squares Problems
A number of functions are provided for factorizing a general rectangular by matrix , as the product of an orthogonal matrix (unitary if complex) and a triangular (or possibly trapezoidal) matrix.
A real matrix
is
orthogonal if
; a complex matrix
is
unitary if
. Orthogonal or unitary matrices have the important property that they leave the
-norm of a vector invariant, so that
if
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.
2.2.1 factorization
The most common, and best known, of the factorizations is the
factorization given by
where
is an
by
upper triangular matrix and
is an
by
orthogonal (or unitary) matrix. If
is of full rank
, then
is nonsingular. It is sometimes convenient to write the factorization as
which reduces to
where
consists of the first
columns of
, and
the remaining
columns.
If
,
is trapezoidal, and the factorization can be written
where
is upper triangular and
is rectangular.
The
factorization can be used to solve the linear least squares problem
(1) when
and
is of full rank, since
where
and
is an
-element vector. Then
is the solution of the upper triangular system
The residual vector
is given by
The residual sum of squares
may be computed without forming
explicitly, since
2.2.2 factorization
The
factorization is given by
where
is
by
lower triangular,
is
by
orthogonal (or unitary),
consists of the first
rows of
, and
the remaining
rows.
The
factorization of
is essentially the same as the
factorization of
(
if
is complex), since
The
factorization may be used to find a minimum norm solution of an underdetermined system of linear equations
where
is
by
with
and has rank
. The solution is given by
2.2.3 factorization with column pivoting
To solve a linear least squares problem
(1) when
is not of full rank, or the rank of
is in doubt, we can perform either a
factorization with column pivoting or a singular value decomposition.
The
factorization with column pivoting is given by
where
and
are as before and
is a (real) permutation matrix, chosen (in general) so that
and moreover, for each
,
If we put
where
is the leading
by
upper triangular sub-matrix of
then, in exact arithmetic, if
, the whole of the sub-matrix
in rows and columns
to
would be zero. In numerical computation, the aim must be to determine an index
, such that the leading sub-matrix
is well-conditioned, and
is negligible, so that
Then
is the effective rank of
. See
Golub and Van Loan (2012) 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
where
consists of just the first
elements of
.
2.2.4 Complete orthogonal factorization
The
factorization with column pivoting does not enable us to compute a
minimum norm solution to a rank-deficient linear least squares problem, unless
. However, by applying for further orthogonal (or unitary) transformations from the right to the upper trapezoidal matrix
,
can be eliminated:
This gives the
complete orthogonal factorization
from which the minimum norm solution can be obtained as
2.2.5 Updating a factorization
Section 2.2.1 gave the forms of the
factorization of an
by
matrix
for the two cases
and
. Taking first the case
, the least squares solution of
is the solution of
If the original system is now augmented by the addition of
rows so that we require the solution of
where
is
by
, then this is equivalent to finding the least squares solution of
This now requires the factorization of the by triangular-rectangular matrix .
For the case
, the least squares solution of the augmented system reduces to
where
is pentagonal.
In both cases can be written as a special case of a triangular-pentagonal matrix consisting of an upper triangular part on top of a rectangular part which is itself on top of a trapezoidal part. In the first case there is no trapezoidal part, in the second case a zero upper triangular part can be added, and more generally the two cases can be combined.
2.2.6 Other factorizations
The
and
factorizations are given by
and
The factorizations are less commonly used than either the or factorizations described above, but have applications in, for example, the computation of generalized factorizations.
2.3 The Singular Value Decomposition
The
singular value decomposition (SVD) of an
by
matrix
is given by
where
and
are orthogonal (unitary) and
is an
by
diagonal matrix with real diagonal elements,
, such that
The
are the
singular values of
and the first
columns of
and
are the
left and
right singular vectors of
. The singular values and singular vectors satisfy
where
and
are the
th columns of
and
respectively.
The computation proceeds in the following stages.
| 1. |
The matrix is reduced to bidiagonal form if is real ( if is complex), where and are orthogonal (unitary if is complex), and is real and upper bidiagonal when and lower bidiagonal when , so that is nonzero only on the main diagonal and either on the first superdiagonal (if ) or the first subdiagonal (if ). |
| 2. |
The SVD of the bidiagonal matrix is computed as , where and are orthogonal and is diagonal as described above. The singular vectors of are then and . |
If , it may be more efficient to first perform a factorization of , and then compute the SVD of the by matrix , since if and , then the SVD of is given by .
Similarly, if , it may be more efficient to first perform an factorization of .
This chapter supports two primary algorithms for computing the SVD of a bidiagonal matrix. They are:
| (i) |
the divide and conquer algorithm; |
| (ii) |
the algorithm. |
The divide and conquer algorithm is much faster than the algorithm if singular vectors of large matrices are required.
2.4 The Singular Value Decomposition and Least Squares Problems
The SVD may be used to find a minimum norm solution to a (possibly) rank-deficient linear least squares problem
(1). The effective rank,
, of
can be determined as the number of singular values which exceed a suitable threshold. Let
be the leading
by
sub-matrix of
, and
be the matrix consisting of the first
columns of
. Then the solution is given by
where
consists of the first
elements of
.
2.5 Generalized Linear Least Squares Problems
The simple type of linear least squares problem described in
Section 2.1 can be generalized in various ways.
| 1. |
Linear least squares problems with equality constraints:
where is by and is by , with . The equations may be regarded as a set of equality constraints on the problem of minimizing . Alternatively the problem may be regarded as solving an overdetermined system of equations
where some of the equations (those involving ) are to be solved exactly, and the others (those involving ) are to be solved in a least squares sense. The problem has a unique solution on the assumptions that has full row rank and the matrix has full column rank . (For linear least squares problems with inequality constraints, refer to Chapter e04.) |
| 2. |
General Gauss–Markov linear model problems:
where is by and is by , with . When , the problem reduces to an ordinary linear least squares problem. When is square and nonsingular, it is equivalent to a weighted linear least squares problem:
The problem has a unique solution on the assumptions that has full column rank , and the matrix has full row rank . Unless 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. |
2.6 Generalized Orthogonal Factorization and Generalized Linear Least Squares Problems
2.6.1 Generalized Factorization
The
generalized (GQR) factorization of an
by
matrix
and an
by
matrix
is given by the pair of factorizations
where
and
are respectively
by
and
by
orthogonal matrices (or unitary matrices if
and
are complex).
has the form
or
where
is upper triangular.
has the form
or
where
or
is upper triangular.
Note that if
is square and nonsingular, the GQR factorization of
and
implicitly gives the
factorization of the matrix
:
without explicitly computing the matrix inverse
or the product
(remembering that the inverse of an invertible upper triangular matrix and the product of two upper triangular matrices is an upper triangular matrix).
The GQR factorization can be used to solve the general (Gauss–Markov) linear model problem (GLM) (see
Section 2.5, but note that
and
are dimensioned differently there as
by
and
by
respectively). Using the GQR factorization of
and
, we rewrite the equation
as
We partition this as
where
The GLM problem is solved by setting
from which we obtain the desired solutions
2.6.2 Generalized Factorization
The
generalized (GRQ) factorization of an
by
matrix
and a
by
matrix
is given by the pair of factorizations
where
and
are respectively
by
and
by
orthogonal matrices (or unitary matrices if
and
are complex).
has the form
or
where
or
is upper triangular.
has the form
or
where
is upper triangular.
Note that if
is square and nonsingular, the GRQ factorization of
and
implicitly gives the
factorization of the matrix
:
without explicitly computing the matrix
or the product
(remembering that the inverse of an invertible upper triangular matrix and the product of two upper triangular matrices is an upper triangular matrix).
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
and
(note that
and
have swapped roles), written as
We write the linear equality constraints
as
which we partition as:
Therefore
is the solution of the upper triangular system
We partition this expression as:
where
.
To solve the LSE problem, we set
which gives
as the solution of the upper triangular system
Finally, the desired solution is given by
2.6.3 Generalized Singular Value Decomposition (GSVD)
The
generalized (or quotient) singular value decomposition of an
by
matrix
and a
by
matrix
is given by the pair of factorizations
The matrices in these factorizations have the following properties:
| – |
is by , is by , is by , and all three matrices are orthogonal. If and are complex, these matrices are unitary instead of orthogonal, and should be replaced by in the pair of factorizations. |
| – |
is by , upper triangular and nonsingular. is by (in other words, the is an by zero matrix). The integer is the rank of , and satisfies . |
| – |
is by , is by , both are real, non-negative and diagonal, and . Write and , where and lie in the interval from to . The ratios are called the generalized singular values of the pair , . If , then the generalized singular value is infinite. |
and
have the following detailed structures, depending on whether
or
. In the first case,
, then
Here is the rank of , , and are diagonal matrices satisfying , and is nonsingular. We may also identify , , for ,
, and , for . Thus, the first generalized singular values are infinite, and the remaining generalized singular values are finite.
In the second case, when
,
and
Again, is the rank of , , and are diagonal matrices satisfying , and is nonsingular, and we may identify , , for , , , , for and . Thus, the first generalized singular values are infinite, and the remaining generalized singular values are finite.
Here are some important special cases of the generalized singular value decomposition. First, if
is square and nonsingular, then
and the generalized singular value decomposition of
and
is equivalent to the singular value decomposition of
, where the singular values of
are equal to the generalized singular values of the pair
,
:
Second, for the matrix
, where
if the columns of
are orthonormal, then
,
and the generalized singular value decomposition of
and
is equivalent to the CS (Cosine–Sine) decomposition of
:
Third, the generalized eigenvalues and eigenvectors of
can be expressed in terms of the generalized singular value decomposition: Let
Therefore, the columns of
are the eigenvectors of
, and ‘nontrivial’ eigenvalues are the squares of the generalized singular values (see also
Section 2.8). ‘Trivial’ eigenvalues are those corresponding to the leading
columns of
, which span the common null space of
and
. The ‘trivial eigenvalues’ are not well defined.
2.6.4 The Full CS Decomposition of Orthogonal Matrices
In
Section 2.6.3 the CS (Cosine-Sine) decomposition of an orthogonal matrix partitioned into two submatrices
and
was given by
The full CS decomposition of an
by
orthogonal matrix
partitions
into four submatrices and factorizes as
where,
is a
by
submatrix (which implies the dimensions of
,
and
);
,
,
and
are orthogonal matrices of dimensions
,
,
and
respectively;
is the
by
single-diagonal matrix
is the
by
single-diagonal matrix
is the
by
single-diagonal matrix
and,
is the
by
single-diagonal matrix
where
and the missing zeros remind us that either the column or the row is missing. The
by
diagonal matrices
and
are such that
.
This is equivalent to the simultaneous singular value decomposition of the four submatrices , , and .
2.7 Symmetric Eigenvalue Problems
The
symmetric eigenvalue problem is to find the
eigenvalues,
, and corresponding
eigenvectors,
, such that
For the
Hermitian eigenvalue problem we have
For both problems the eigenvalues
are real.
When all eigenvalues and eigenvectors have been computed, we write
where
is a diagonal matrix whose diagonal elements are the eigenvalues, and
is an orthogonal (or unitary) matrix whose columns are the eigenvectors. This is the classical
spectral factorization of
.
The basic task of the symmetric eigenproblem functions is to compute values of
and, optionally, corresponding vectors
for a given matrix
. This computation proceeds in the following stages.
| 1. |
The real symmetric or complex Hermitian matrix is reduced to real tridiagonal form
. If is real symmetric this decomposition is with orthogonal and symmetric tridiagonal. If is complex Hermitian, the decomposition is with unitary and , as before, real symmetric tridiagonal. |
| 2. |
Eigenvalues and eigenvectors of the real symmetric tridiagonal matrix are computed. If all eigenvalues and eigenvectors are computed, this is equivalent to factorizing as , where is orthogonal and is diagonal. The diagonal entries of are the eigenvalues of , which are also the eigenvalues of , and the columns of are the eigenvectors of ; the eigenvectors of are the columns of , so that ( when 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 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
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).
2.8 Generalized Symmetric-definite Eigenvalue Problems
This section is concerned with the solution of the generalized eigenvalue problems , , and , where and are real symmetric or complex Hermitian and is positive definite. Each of these problems can be reduced to a standard symmetric eigenvalue problem, using a Cholesky factorization of as either or ( or in the Hermitian case).
With
, we have
Hence the eigenvalues of
are those of
, where
is the symmetric matrix
and
. In the complex case
is Hermitian with
and
.
Table 1 summarises how each of the three types of problem may be reduced to standard form
, and how the eigenvectors
of the original problem may be recovered from the eigenvectors
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 |
Reduction |
Recovery of eigenvectors |
| 1. |
|
,
|
,
|
,
|
| 2. |
|
,
|
,
|
,
|
| 3. |
|
,
|
,
|
,
|
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
, this may then be solved using the functions described in the previous section. No special functions are needed to recover the eigenvectors
of the generalized problem from the eigenvectors
of the standard problem, because these computations are simple applications of Level 2 or Level 3 BLAS
(see
Chapter f16).
2.9 Packed Storage for Symmetric Matrices
Functions 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 ; 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.
Functions designed for packed storage are usually less efficient, especially on high-performance computers, so there is a trade-off between storage and efficiency.
2.10 Band Matrices
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
and the matrices
and
are additionally banded, the matrix
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
thus:
where
is an orthogonal matrix chosen to ensure that
has bandwidth no greater than that of
.
A further refinement is possible when
and
are banded, which halves the amount of work required to form
. Instead of the standard Cholesky factorization of
as
or
, we use a
split Cholesky factorization
, where
with
upper triangular and
lower triangular of order approximately
;
has the same bandwidth as
.
2.11 Nonsymmetric Eigenvalue Problems
The
nonsymmetric eigenvalue problem is to find the
eigenvalues,
, and corresponding
eigenvectors,
, such that
More precisely, a vector
as just defined is called a
right eigenvector of
, and a vector
satisfying
is called a
left eigenvector of
.
A real matrix may have complex eigenvalues, occurring as complex conjugate pairs.
This problem can be solved via the
Schur factorization of
, defined in the real case as
where
is an orthogonal matrix and
is an upper quasi-triangular matrix with
by
and
by
diagonal blocks, the
by
blocks corresponding to complex conjugate pairs of eigenvalues of
. In the complex case, the Schur factorization is
where
is unitary and
is a complex upper triangular matrix.
The columns of are called the Schur vectors. For each (), the first columns of form an orthonormal basis for the invariant subspace corresponding to the first eigenvalues on the diagonal of . 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 eigenvalues occupy the leading positions on the diagonal of .
The two basic tasks of the nonsymmetric eigenvalue functions are to compute, for a given matrix , all values of and, if desired, their associated right eigenvectors and/or left eigenvectors , and the Schur factorization.
These two basic tasks can be performed in the following stages.
| 1. |
A general matrix is reduced to upper Hessenberg form
which is zero below the first subdiagonal. The reduction may be written with orthogonal if is real, or with unitary if is complex. |
| 2. |
The upper Hessenberg matrix is reduced to Schur form , giving the Schur factorization (for real) or (for complex). The matrix (the Schur vectors of ) may optionally be computed as well. Alternatively may be postmultiplied into the matrix determined in stage , to give the matrix , the Schur vectors of . The eigenvalues are obtained from the diagonal elements or diagonal blocks of . |
| 3. |
Given the eigenvalues, the eigenvectors may be computed in two different ways. Inverse iteration can be performed on to compute the eigenvectors of , and then the eigenvectors can be multiplied by the matrix in order to transform them to eigenvectors of . Alternatively the eigenvectors of can be computed, and optionally transformed to those of or if the matrix or 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.
2.12 Generalized Nonsymmetric Eigenvalue Problem
The
generalized nonsymmetric eigenvalue problem is to find the eigenvalues,
, and corresponding
eigenvectors,
, such that
More precisely, a vector
as just defined is called a
right eigenvector of the matrix pair
, and a vector
satisfying
is called a
left eigenvector of the matrix pair
.
If
is singular then the problem has one or more
infinite eigenvalues
, corresponding to
. Note that if
is nonsingular, then the equivalent problem
is perfectly well defined and an infinite eigenvalue corresponds to
. To deal with both finite (including zero) and infinite eigenvalues, the functions in this chapter do not compute
explicitly, but rather return a pair of numbers
such that if
and if
and
then
.
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
the set of all the matrices of the form
is called a matrix
pencil and
and
are said to be an eigenvalue and eigenvector of the pencil
. If
and
are both singular and share a common null space then
so that the pencil
is
singular for all
. In other words any
can be regarded as an eigenvalue. In exact arithmetic a singular pencil will have
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
defined in the real case as
where
and
are orthogonal,
is upper triangular with non-negative diagonal elements and
is upper quasi-triangular with
by
and
by
diagonal blocks, the
by
blocks corresponding to complex conjugate pairs of eigenvalues. In the complex case, the generalized Schur factorization is
where
and
are unitary and
and
are upper triangular, with
having real non-negative diagonal elements. The columns of
and
are called respectively the
left and
right generalized Schur vectors and span pairs of
deflating subspaces of
and
, which are a generalization of invariant subspaces.
It is possible to order the generalized Schur factorization so that any desired set of eigenvalues correspond to the leading positions on the diagonals of the pair .
The two basic tasks of the generalized nonsymmetric eigenvalue functions are to compute, for a given pair , all values of and, if desired, their associated right eigenvectors and/or left eigenvectors , and the generalized Schur factorization.
These two basic tasks can be performed in the following stages.
| 1. |
The matrix pair is reduced to generalized upper Hessenberg form ,
where is upper Hessenberg (zero below the first subdiagonal) and is upper triangular. The reduction may be written as in the real case with and orthogonal, and in the complex case with and unitary. |
| 2. |
The generalized upper Hessenberg form is reduced to the generalized
Schur form using the generalized Schur factorization , in the real case with and orthogonal, and in the complex case. The generalized Schur vectors of are given by , . The eigenvalues are obtained from the diagonal elements (or blocks) of the pair . |
| 3. |
Given the eigenvalues, the eigenvectors of the pair can be computed, and optionally transformed to those of or . |
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.
2.13 The Sylvester Equation and the Generalized Sylvester Equation
The Sylvester equation is a matrix equation of the form
where
,
, and
are given matrices with
being
by
,
an
by
matrix and
, and the solution matrix
,
by
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 functions in this chapter allows the solution of the general Sylvester equation.
Functions are also provided for solving a special case of the generalized Sylvester equations
where
,
and
are given matrix pairs, and
and
are the solution matrices.
2.14 Error and Perturbation Bounds and Condition Numbers
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 functions 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
(or
), which grows as a function of the matrix dimension
(or matrix dimensions
and
). 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;
is often true, although generally only much weaker bounds can be actually proved. We normally describe
as a ‘modestly growing’ function of
. For detailed derivations of various
, see
Golub and Van Loan (2012) and
Wilkinson (1965).
For linear equation (see
Chapter f07) and least squares solvers, we consider bounds on the relative error
in the computed solution
, where
is the true solution. For eigenvalue problems we consider bounds on the error
in the
th computed eigenvalue
, where
is the true
th eigenvalue. For singular value problems we similarly consider bounds
.
Bounding the error in computed eigenvectors and singular vectors
is more subtle because these vectors are not unique: even though we restrict
and
, we may still multiply them by arbitrary constants of absolute value
. So to avoid ambiguity we bound the
angular difference between
and the true vector
, so that
Here
is in the standard range:
. When
is small, we can choose a constant
with absolute value
so that
.
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
is any vector in the space,
is also in the space, where
is the matrix. Again, we will use angle to measure the difference between a computed space
and the true space
:
may be computed as follows. Let
be a matrix whose columns are orthonormal and
. Similarly let
be an orthonormal matrix with columns spanning
. Then
Finally, we remark on the accuracy of the bounds when they are large. Relative errors like
and angular errors like
are only of interest when they are much less than
. Some stated bounds are not strictly true when they are close to
, 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
. 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 functions 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 functions, implement these code fragments.
2.14.1 Least squares problems
The conventional error analysis of linear least squares problems goes as follows. The problem is to find the minimizing . Let be the solution computed using one of the methods described above. We discuss the most common case, where is overdetermined (i.e., has more rows than columns) and has full rank.
Then the computed solution
has a small normwise backward error. In other words
minimizes
, where
and
is a modestly growing function of
and
is the
machine precision. Let
,
, and
. Then if
is small enough, the error
is bounded by
If
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 (2012) 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
By solving this linear system (see
Chapter f07) component-wise error bounds can also be obtained (see
Arioli et al. (1989)).
2.14.2 The singular value decomposition
The usual error analysis of the SVD algorithm is as follows (see
Golub and Van Loan (2012)).
The computed SVD,
, is nearly the exact SVD of
, i.e.,
is the true SVD, so that
and
are both orthogonal, where
,
, and
. Here
is a modestly growing function of
and
and
is the
machine precision. Each computed singular value
differs from the true
by an amount satisfying the bound
Thus large singular values (those near
) are computed to high relative accuracy and small ones may not be.
The angular difference between the computed left singular vector
and the true
satisfies the approximate bound
where
is the
absolute gap between
and the nearest other singular value. Thus, if
is close to other singular values, its corresponding singular vector
may be inaccurate. The same bound applies to the computed right singular vector
and the true vector
. The gaps may be easily obtained from the computed singular values.
Let
be the space spanned by a collection of computed left singular vectors
, where
is a subset of the integers from
to
. Let
be the corresponding true space. Then
where
is the absolute gap between the singular values in
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
even if its individual singular vectors are ill-conditioned. The same bound applies to a set of right singular vectors
.
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
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
can introduce additional errors, so the following bounds for the bidiagonal case do not apply to the dense case.
Using the functions 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
The computed left singular vector
has an angular error at most about
where
is the
relative gap between
and the nearest other singular value. The same bound applies to the right singular vector
and
. 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.
2.14.3 The symmetric eigenproblem
The usual error analysis of the symmetric eigenproblem is as follows (see
Parlett (1998)).
The computed eigendecomposition
is nearly the exact eigendecomposition of
, i.e.,
is the true eigendecomposition so that
is orthogonal, where
and
and
is a modestly growing function of
and
is the
machine precision. Each computed eigenvalue
differs from the true
by an amount satisfying the bound
Thus large eigenvalues (those near
) are computed to high relative accuracy and small ones may not be.
The angular difference between the computed unit eigenvector
and the true
satisfies the approximate bound
if
is small enough, where
is the
absolute gap between
and the nearest other eigenvalue. Thus, if
is close to other eigenvalues, its corresponding eigenvector
may be inaccurate. The gaps may be easily obtained from the computed eigenvalues.
Let
be the invariant subspace spanned by a collection of eigenvectors
, where
is a subset of the integers from
to
. Let
be the corresponding true subspace. Then
where
is the absolute gap between the eigenvalues in
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
even if its individual eigenvectors are ill-conditioned.
In the special case of a real symmetric tridiagonal matrix
, functions in this chapter can compute the eigenvalues and eigenvectors much more accurately. See
Anderson et al. (1999) for further details.
2.14.4 The generalized symmetric-definite eigenproblem
The three types of problem to be considered are
,
and
. In each case
and
are real symmetric (or complex Hermitian) and
is positive definite. We consider each case in turn, assuming that functions 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
is the
absolute gap between
and the nearest other eigenvalue.
| 1. |
. The computed eigenvalues can differ from the true eigenvalues by an amount
The angular difference between the computed eigenvector and the true eigenvector is
|
| 2. |
or . The computed eigenvalues can differ from the true eigenvalues by an amount
The angular difference between the computed eigenvector and the true eigenvector is
|
These error bounds are large when
is ill-conditioned with respect to inversion (
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
differ widely in magnitude, as for example with a
graded matrix.
| 1. |
. Let
be a diagonal matrix. Then replace by and by in the above bounds. |
| 2. |
or . Let be a diagonal matrix. Then replace by and by in the above bounds. |
Further details can be found in
Anderson et al. (1999).
2.14.5 The nonsymmetric eigenproblem
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 be the th computed eigenvalue and the th true eigenvalue. Let be the corresponding computed right eigenvector, and the true right eigenvector (so ). If is a subset of the integers from to , we let denote the average of the selected eigenvalues:
, and similarly for . We also let denote the subspace spanned by ; it is called a right invariant subspace because if is any vector in then is also in . 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 , where . Some of the bounds are stated in terms of and others in terms of ; one may use for either quantity.
Functions are provided so that, for each (
) pair the two values
and
, or for a selected subset
of eigenvalues the values
and
can be obtained, for which the error bounds in
Table 2 are true for sufficiently small
, (which is why they are called asymptotic):
| Simple eigenvalue |
|
| Eigenvalue cluster |
|
| Eigenvector |
|
| Invariant subspace |
|
Table 2
Asymptotic error bounds for the nonsymmetric eigenproblem
If the problem is ill-conditioned, the asymptotic bounds may only hold for extremely small
. The global error bounds of
Table 3 are guaranteed to hold for all
:
Simple
eigenvalue |
|
Holds for all |
Eigenvalue
cluster |
|
Requires |
| Eigenvector |
|
Requires |
Invariant
subspace |
|
Requires |
Table 3
Global error bounds for the nonsymmetric eigenproblem
2.14.6 Balancing and condition for the nonsymmetric eigenproblem
There are two preprocessing steps one may perform on a matrix
in order to make its eigenproblem easier. The first is
permutation, or reordering the rows and columns to make
more nearly upper triangular (closer to Schur form):
, where
is a permutation matrix. If
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
to make the rows and columns of
more nearly equal in norm:
. 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).
2.14.7 The generalized nonsymmetric eigenvalue problem
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
, where
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. Functions are provided to compute estimates of reciprocal conditions numbers for eigenvalues and eigenspaces.
2.14.8 Balancing the generalized eigenvalue problem
As with the standard nonsymmetric eigenvalue problem, there are two preprocessing steps one may perform on a matrix pair
in order to make its eigenproblem easier; permutation and scaling, which together are referred to as balancing, as indicated in the following two steps.
| 1. |
The balancing function first attempts to permute and to block upper triangular form by a similarity transformation:
where is a permutation matrix, , , and are upper triangular. Then the diagonal elements of the matrix and are generalized eigenvalues of . The rest of the generalized eigenvalues are given by the matrix pair . Subsequent operations to compute the eigenvalues of need only be applied to the matrix ; this can save a significant amount of work if is smaller than the original matrix pair . If no suitable permutation exists (as is often the case), then there is no gain in efficiency or accuracy. |
| 2. |
The balancing function applies a diagonal similarity transformation to , to make the rows and columns of as close as possible in the norm:
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. |
2.14.9 Other problems
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).
2.15 Block Partitioned Algorithms
A number of the functions 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 f16),
which are the key to achieving high performance on many modern computers. In the case of the
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 (2012) 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 argument, 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 functions. Values in the range to are typical.
On more conventional machines there is often no advantage from using a block partitioned algorithm, and then the functions use an
unblocked algorithm (effectively a block size of
), relying solely on calls to the Level 2 BLAS
(see
Chapter f16 again).
3 Recommendations on Choice and Use of Available Functions
3.1 Available Functions
The tables in the following sub-sections show the functions which are provided for performing different computations on different types of matrices. Each entry in the table gives the NAG function short name.
Black box (or driver) functions 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 functions.
3.1.1 Driver functions
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)
| 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
|
f08pac
f08pbc
f08nac
f08nbc
|
f08pnc
f08ppc
f08nnc
f08npc
|
3.1.1.5 Singular value decomposition (SVD)
| 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
|
f08kbc
f08kdc
f08kjc
f08khc |
f08kpc
f08krc |
3.1.1.6 Generalized symmetric definite eigenvalue problems (GSEP)
3.1.1.7 Generalized nonsymmetric eigenvalue problem (GNEP)
| 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
|
f08xac
f08xbc
f08wac
f08wbc
|
f08xnc
f08xpc
f08wnc
f08wpc
|
3.1.1.8 Generalized singular value decomposition (GSVD)
| Function and storage scheme |
real |
complex |
| singular values/vectors |
f08vac |
f08vnc |
3.1.2 Computational functions
It is possible to solve problems by calling two or more functions in sequence. Some common sequences of functions are indicated in the tables in the following sub-sections; an asterisk () against a function name means that the sequence of calls is illustrated in the example program for that function.
3.1.2.1 Orthogonal factorizations
Functions are provided for factorization (with and without column pivoting), and for , and factorizations (without pivoting only), of a general real or complex rectangular matrix. A function is also provided for the factorization of a real or complex upper trapezoidal matrix. (LAPACK refers to this as the factorization.)
The factorization functions do not form the matrix
explicitly, but represent it as a product of elementary reflectors (see
Section 3.3.6). Additional functions are provided to generate all or part of
explicitly if it is required, or to apply
in its factored form to another matrix (specifically to compute one of the matrix products
,
,
or
with
replaced by
if
and
are complex).
| |
Factorize
without
pivoting |
Factorize
with
pivoting |
Factorize
(blocked) |
Generate
matrix |
Apply
matrix |
Apply
(blocked) |
factorization,
real matrices |
f08aec |
f08bfc |
f08abc |
f08afc |
f08agc |
f08acc |
factorization,
real triangular-pentagonal |
|
|
f08bbc |
|
|
f08bcc |
factorization,
real matrices |
f08ahc |
|
|
f08ajc |
f08akc |
|
factorization,
real matrices |
f08cec |
|
|
f08cfc |
f08cgc |
|
factorization,
real matrices |
f08chc |
|
|
f08cjc |
f08ckc |
|
factorization,
real upper trapezoidal matrices |
f08bhc |
|
|
|
f08bkc |
|
factorization,
complex matrices |
f08asc |
f08btc |
f08apc |
f08atc |
f08auc |
f08aqc |
factorization,
complex triangular-pentagonal |
|
|
f08bpc |
|
|
f08bqc |
factorization,
complex matrices |
f08avc |
|
|
f08awc |
f08axc |
|
factorization,
complex matrices |
f08csc |
|
|
f08ctc |
f08cuc |
|
factorization,
complex matrices |
f08cvc |
|
|
f08cwc |
f08cxc |
|
factorization,
complex upper trapezoidal matrices |
f08bvc |
|
|
|
f08bxc |
|
To solve linear least squares problems, as described in
Sections 2.2.1 or
2.2.3, functions based on the
factorization can be used:
| real data, full-rank problem |
f08aac, f08aec and f08agc, f08abc and f08acc, f16yjc |
| complex data, full-rank problem |
f08anc, f08asc and f08auc, f08apc and f08aqc, f16zjc |
| real data, rank-deficient problem |
f08bfc*, f16yjc, f08agc |
| complex data, rank-deficient problem |
f08btc*, f16zjc, f08auc |
To find the minimum norm solution of under-determined systems of linear equations, as described in
Section 2.2.2, functions based on the
factorization can be used:
3.1.2.2 Generalized orthogonal factorizations
Functions are provided for the generalized
and
factorizations of real and complex matrix pairs.
| |
Factorize |
| Generalized factorization, real matrices |
f08zec |
| Generalized factorization, real matrices |
f08zfc |
| Generalized factorization, complex matrices |
f08zsc |
| Generalized factorization, complex matrices |
f08ztc |
3.1.2.3 Singular value problems
Functions are provided to reduce a general real or complex rectangular matrix
to real bidiagonal form
by an orthogonal transformation
(or by a unitary transformation
if
is complex). Different functions allow a full matrix
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 functions for reducing full matrices do not form the matrix or explicitly; additional functions are provided to generate all or part of them, or to apply them to another matrix, as with the functions for orthogonal factorizations. Explicit generation of or is required before using the bidiagonal algorithm to compute left or right singular vectors of .
The functions for reducing band matrices have options to generate or if required.
Further functions are provided to compute all or part of the singular value decomposition of a real bidiagonal matrix; the same functions can be used to compute the singular value decomposition of a real or complex matrix that has been reduced to bidiagonal form.
Given the singular values,
f08flc 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)
Rectangular matrix (banded)
To use the singular value decomposition to solve a linear least squares problem, as described in
Section 2.4, the following functions are required:
3.1.2.4 Generalized singular value decomposition
Functions are provided to compute the generalized SVD of a real or complex matrix pair
in upper trapezoidal form. Functions are also provided to reduce a general real or complex matrix pair to the required upper trapezoidal form.
Functions are provided for the full CS decomposition of orthogonal and unitary matrices expressed as
by
partitions of submatrices. For real orthogonal matrices the CS decomposition is performed by
f08rac, while for unitary matrices the equivalent function is
f08rnc.
3.1.2.5 Symmetric eigenvalue problems
Functions are provided to reduce a real symmetric or complex Hermitian matrix
to real tridiagonal form
by an orthogonal similarity transformation
(or by a unitary transformation
if
is complex). Different functions allow a full matrix
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 functions for reducing full matrices do not form the matrix explicitly; additional functions are provided to generate , or to apply it to another matrix, as with the functions for orthogonal factorizations. Explicit generation of is required before using the algorithm to find all the eigenvectors of ; application of to another matrix is required after eigenvectors of have been found by inverse iteration, in order to transform them to eigenvectors of .
The functions for reducing band matrices have an option to generate
if required.
Given the eigenvalues,
f08flc is provided to compute the reciprocal condition numbers for the eigenvectors of a real symmetric or complex Hermitian matrix.
A variety of functions are provided to compute eigenvalues and eigenvectors of the real symmetric tridiagonal matrix , some computing all eigenvalues and eigenvectors, some computing selected eigenvalues and eigenvectors. The same functions 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 Symmetric or Complex Hermitian
| all eigenvalues (root-free algorithm) |
f08jfc |
| all eigenvalues (root-free algorithm called by divide-and-conquer) |
f08jcc or f08jhc |
| all eigenvalues (RRR) |
f08jlc |
| selected eigenvalues (bisection) |
f08jjc |
The original (non-reduced) matrix is Real Symmetric
| all eigenvalues and eigenvectors ( algorithm) |
f08jec |
| all eigenvalues and eigenvectors (divide-and-conquer) |
f08jcc or f08jhc |
| all eigenvalues and eigenvectors (RRR) |
f08jlc |
| all eigenvalues and eigenvectors (positive definite case) |
f08jgc |
| selected eigenvectors (inverse iteration) |
f08jkc |
The original (non-reduced) matrix is Complex Hermitian
| all eigenvalues and eigenvectors ( algorithm) |
f08jsc |
| all eigenvalues and eigenvectors (divide and conquer) |
f08jvc |
| all eigenvalues and eigenvectors (RRR) |
f08jyc |
| all eigenvalues and eigenvectors (positive definite case) |
f08juc |
| selected eigenvectors (inverse iteration) |
f08jxc |
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)
Real Symmetric matrix (packed storage)
Real Symmetric banded matrix
| all eigenvalues and eigenvectors (using divide-and-conquer) |
f08hcc |
| all eigenvalues and eigenvectors (using algorithm) |
f08hec*, f08jec |
Complex Hermitian matrix (standard storage)
Complex Hermitian matrix (packed storage)
Complex Hermitian banded matrix
| all eigenvalues and eigenvectors (using divide-and-conquer) |
f08hqc |
| all eigenvalues and eigenvectors (using algorithm) |
f08hsc*, f08jsc |
3.1.2.6 Generalized symmetric-definite eigenvalue problems
Functions are provided for reducing each of the problems
,
or
to an equivalent standard eigenvalue problem
. Different functions allow the matrices to be stored either conventionally or in packed storage. The positive definite matrix
must first be factorized using a function from
Chapter f07. There is also a function which reduces the problem
where
and
are banded, to an equivalent banded standard eigenvalue problem; this uses a split Cholesky factorization for which a function in
Chapter f08 is provided.
| |
Reduce to
standard problem |
Reduce to
standard problem
(packed storage) |
Reduce to
standard problem
(band matrices) |
| real symmetric matrices |
f08sec |
f08tec |
f08uec |
| complex Hermitian matrices |
f08ssc |
f08tsc |
f08usc |
The equivalent standard problem can then be solved using the functions discussed in
Section 3.1.2.5. For example, to compute all the eigenvalues, the following functions must be called:
| real symmetric-definite problem |
f07fdc, f08sec*, f08fec, f08jfc |
| real symmetric-definite problem, packed storage |
f07gdc, f08tec*, f08gec, f08jfc |
| real symmetric-definite banded problem |
f08ufc*, f08uec*, f08hec, f08jfc |
| complex Hermitian-definite problem |
f07frc, f08ssc*, f08fsc, f08jfc |
| complex Hermitian-definite problem, packed storage |
f07grc, f08tsc*, f08gsc, f08jfc |
| complex Hermitian-definite banded problem |
f08utc*, f08usc*, f08hsc, f08jfc |
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; functions from
Chapter f16
may be used for this.
3.1.2.7 Nonsymmetric eigenvalue problems
Functions are provided to reduce a general real or complex matrix to upper Hessenberg form by an orthogonal similarity transformation (or by a unitary transformation if is complex).
These functions do not form the matrix explicitly; additional functions are provided to generate , or to apply it to another matrix, as with the functions for orthogonal factorizations. Explicit generation of is required before using the algorithm on to compute the Schur vectors; application of to another matrix is needed after eigenvectors of have been computed by inverse iteration, in order to transform them to eigenvectors of .
Functions are also provided to balance the matrix before reducing it to Hessenberg form, as described in
Section 2.14.6. Companion functions are required to transform Schur vectors or eigenvectors of the balanced matrix to those of the original matrix.
Functions 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 functions estimate the sensitivities of computed eigenvalues and eigenvectors, as discussed in
Section 2.14.5.
| |
Eigenvalues and
Schur factorization
( algorithm) |
Eigenvectors from
Hessenberg form
(inverse iteration) |
Eigenvectors from
Schur factorization |
Sensitivities of
eigenvalues and
eigenvectors |
| real matrices |
f08pec |
f08pkc |
f08qkc |
f08qlc |
| complex matrices |
f08psc |
f08pxc |
f08qxc |
f08qyc |
Finally functions are provided for reordering the Schur factorization, so that eigenvalues appear in any desired order on the diagonal of the Schur form. The functions
f08qfc and
f08qtc simply swap two diagonal elements or blocks, and may need to be called repeatedly to achieve a desired order. The functions
f08qgc and
f08quc 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 functions 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 |
f08qfc |
f08qgc |
| complex matrices |
f08qtc |
f08quc |
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 |
f08nec, f08nfc*, f08pec |
| real matrix, all eigenvalues and selected eigenvectors |
f08nec, f08ngc, f08pec, f08pkc
|
| real matrix, all eigenvalues and eigenvectors (with balancing) |
f08nhc*, f08nec, f08nfc, f08njc, f08pec, f08pkc |
| complex matrix, all eigenvalues and Schur factorization |
f08nsc, f08ntc*, f08psc |
| complex matrix, all eigenvalues and selected eigenvectors |
f08nsc, f08nuc, f08psc, f08pxc* |
| complex matrix, all eigenvalues and eigenvectors (with balancing) |
f08nvc*, f08nsc, f08ntc, f08nwc, f08psc, f08pxc |
3.1.2.8 Generalized nonsymmetric eigenvalue problems
Functions are provided to reduce a real or complex matrix pair
,
where
is general and
is upper triangular, to generalized upper
Hessenberg form by orthogonal transformations
,
,
(or by unitary transformations
,
,
in the complex case). These functions can optionally return
and/or
. Note that to transform a general matrix pair
to the form
a
factorization of
(
) should first be performed and the matrix
obtained as
(see
Section 3.1.2.1 above).
Functions are also provided to balance a general matrix pair before reducing it to generalized Hessenberg form, as described in
Section 2.14.8. Companion functions are provided to transform vectors of the balanced pair to those of the original matrix pair.
Functions 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 functions estimate the sensitivities of computed eigenvalues and eigenvectors.
| |
Eigenvalues and
generalized Schur
factorization
( algorithm) |
Eigenvectors from
generalized Schur
factorization |
Sensitivities of
eigenvalues and
eigenvectors |
| real matrices |
f08xec |
f08ykc |
f08ylc |
| complex matrices |
f08xsc |
f08yxc |
f08yyc |
Finally, functions are provided for reordering the generalized Schur factorization so that eigenvalues appear in any desired order on the diagonal of the generalized Schur form.
f08yfc and
f08ytc simply swap two diagonal elements or blocks, and may need to be called repeatedly to achieve a desired order.
f08ygc and
f08yuc 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 functions 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 |
f08yfc |
f08ygc |
| complex matrices |
f08ytc |
f08yuc |
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) |
f08aec, f08agc (or f08abc, f08acc), f08wec, f08whc, f08xec* |
| real matrix pair, all eigenvalues and generalized Schur factorization |
f08aec, f08afc, f08agc (or f08abc, f08acc), f08wec, f08xec |
| real matrix pair, all eigenvalues and eigenvectors (with balancing) |
f16qfc, f16qhc, f08aec, f08afc, f08agc (or f08abc, f08acc), f08wec, f08whc, f08xec, f08ykc*, f08wjc |
| complex matrix pair, all eigenvalues (with balancing) |
f08asc, f08auc (or f08apc, f08aqc), f08wsc, f08wvc, f08xsc* |
| complex matrix pair, all eigenvalues and generalized Schur factorization |
f08asc, f08atc, f08auc (or f08apc, f08aqc), f08wsc, f08xsc |
| complex matrix pair, all eigenvalues and eigenvectors (with balancing) |
f16tfc, f16thc, f08asc, f08atc, f08auc (or f08apc, f08aqc), f08wsc, f08wvc, f08xsc, f08yxc*, f08wwc |
3.1.2.9 The Sylvester equation and the generalized Sylvester equation
Functions are provided to solve the real or complex Sylvester equation
, where
and
are upper quasi-triangular if real, or upper triangular if complex. To solve the general form of the Sylvester equation in which
and
are general square matrices,
and
must be reduced to upper (quasi-) triangular form by the Schur factorization, using functions described in
Section 3.1.2.7. For more details, see the documents for the functions listed below.
| |
Solve the Sylvester equation |
| real matrices |
f08qhc |
| complex matrices |
f08qvc |
Functions are also provided to solve the real or complex generalized Sylvester equations
where the pairs
and
are in generalized Schur form. To solve the general form of the generalized Sylvester equation in which
and
are general matrix pairs,
and
must first be reduced to generalized Schur form.
| |
Solve the generalized Sylvester equation |
| real matrices |
f08yhc |
| complex matrices |
f08yvc |
3.2 NAG Names and LAPACK Names
The functions may be called either by their NAG short names or by their NAG long names which contain their double precision LAPACK names.
References to
Chapter f08 functions in the manual normally include the LAPACK double precision names, for example
nag_dgeqrf (f08aec). The LAPACK routine names follow a simple scheme. 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 |
| d |
– real, double precision |
| c |
– complex, single precision |
| z |
– complex, double precision |
|
| – |
the second and third letters yy indicate the type of the matrix or matrix pair (and in some cases the storage scheme):
| bd |
– bidiagonal |
| di |
– diagonal |
| gb |
– general band |
| ge |
– general |
| gg |
– general pair ( may be triangular) |
| hb |
– (complex) Hermitian band |
| he |
– Hermitian |
| hg |
– generalized upper Hessenberg |
| hp |
– Hermitian (packed storage) |
| hs |
– upper Hessenberg |
| op |
– (real) orthogonal (packed storage) |
| or |
– (real) orthogonal |
| pt |
– symmetric or Hermitian positive definite tridiagonal |
| sb |
– (real) symmetric band |
| sp |
– symmetric (packed storage) |
| st |
– (real) symmetric tridiagonal |
| sy |
– symmetric |
| tg |
– triangular pair (one may be quasi-triangular) |
| tp |
– triangular-pentagonal |
| tr |
– triangular (or quasi-triangular) |
| un |
– (complex) unitary |
| up |
– (complex) unitary (packed storage) |
|
| – |
the last three letters zzz indicate the computation performed. For example, qrf is a factorization. |
Thus the function
nag_dgeqrf performs a
factorization of a real general matrix; the corresponding function for a complex general matrix is
nag_zgeqrf.
3.3 Matrix Storage Schemes
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 f07 and
f16,
but different schemes for packed, band and tridiagonal storage are used in a few older functions in
Chapters f01,
f02,
f03 and
f04.
3.3.1 Conventional storage
Please see
Section 3.3.1 in the f07 Chapter Introduction for full details.
3.3.2 Packed storage
Please see
Section 3.3.2 in the f07 Chapter Introduction for full details.
3.3.3 Band storage
Please see
Section 3.3.4 in the f07 Chapter Introduction for full details.
3.3.4 Tridiagonal and bidiagonal matrices
A symmetric tridiagonal or bidiagonal matrix is stored in two one-dimensional arrays, one of length
containing the diagonal elements, and one of length
containing the off-diagonal elements. (Older functions in
Chapter f02 store the off-diagonal elements in elements
of a vector of length
.)
3.3.5 Real diagonal elements of complex matrices
Please see
Section 3.3.6 in the f07 Chapter Introduction for full details.
3.3.6 Representation of orthogonal or unitary matrices
A real orthogonal or complex unitary matrix (usually denoted
) is often represented in the NAG C Library as a product of
elementary reflectors – also referred to as
elementary Householder matrices (usually denoted
). For example,
You need not be aware of the details, because functions are provided to work with this representation, either to generate all or part of
explicitly, or to multiply a given matrix by
or
(
in the complex case) without forming
explicitly.
Nevertheless, the following further details may occasionally be useful.
An elementary reflector (or elementary Householder matrix)
of order
is a unitary matrix of the form
where
is a scalar, and
is an
-element vector, with
;
is often referred to as the
Householder vector. Often
has several leading or trailing zero elements, but for the purpose of this discussion assume that
has no such special structure.
There is some redundancy in the representation
, 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 functions in
Chapters f01,
f02 and
f04)
sets
; hence
need not be stored. In real arithmetic,
, except that
implies
.
In complex arithmetic,
may be complex, and satisfies
and
. Thus a complex
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
can be computed so that
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.
3.4 Argument Conventions
3.4.1 Option arguments
In addition to the order argument of type Nag_OrderType, most functions in this Chapter have one or more option arguments of various types; only options of the correct type may be supplied.
For example,
nag_dsytrd(Nag_RowMajor,Nag_Upper,...)
3.4.2 Problem dimensions
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.
4 Decision Trees
The following decision trees are principally for the computation (general purpose) functions.
4.1 General Purpose Functions (eigenvalues and eigenvectors)
Tree 1: Real Symmetric Eigenvalue Problems
Tree 2: Real Generalized Symmetric-definite Eigenvalue Problems
| Are eigenvalues only required? |
|
Are all the eigenvalues required? |
|
Are and band matrices? |
|
f08ufc, f08uec, f08hec and f08jfc |
| yes | yes | yes |
| | no | | | no | | | no | |
| |
|
Are and stored with one triangle as a linear array? |
|
f07gdc, f08tec, f08gec and f08jfc |
| | | yes |
| | | | no | |
| |
|
f07fdc, f08sec, f08fec and f08jfc |
| | |
| | |
| |
Are and band matrices? |
|
f08ufc, f08uec, f08hec and f08jjc |
| | yes |
| | | no | |
| |
Are and stored with one triangle as a linear array? |
|
f07gdc, f08tec, f08gec and f08jjc |
| | yes |
| | | no | |
| |
f07fdc, f08sec, f08gec and f08jjc |
| |
| |
| Are all eigenvalues and eigenvectors required? |
|
Are and stored with one triangle as a linear array? |
|
f07gdc, f08tec, f08gec, f08gfc, f08jec and f16plc |
| yes | yes |
| | no | | | no | |
| |
f07fdc, f08sec, f08fec, f08ffc, f08jec and f16yjc |
| |
| |
| Are and band matrices? |
|
f08ufc, f08uec, f08hec, f08jkc and f16yjc |
| yes |
| | no | |
| Are and stored with one triangle as a linear array? |
|
f07gdc, f08tec, f08gec, f08jjc, f08jkc, f08ggc and f16plc |
| yes |
| | no | |
| f07fdc, f08sec, f08fec, f08jjc, f08jkc, f08fgc and f16yjc |
|
Note: the functions for band matrices only handle the problem
; the other functions handle all three types of problems (
,
or
) except that, if the problem is
and eigenvectors are required,
f16phc must be used instead of
f16plc and
f16yfc instead of
f16yjc.
Tree 3: Real Nonsymmetric Eigenvalue Problems
| Are eigenvalues required? |
|
Is an upper Hessenberg matrix? |
|
f08pec |
| yes | yes |
| | no | | | no | |
| |
f08nac or f08nbc or (f08nhc, f08nec and f08pec) |
| |
| |
| Is the Schur factorization of required? |
|
Is an upper Hessenberg matrix? |
|
f08pec |
| yes | yes |
| | no | | | no | |
| |
f08nbc or (f08nec, f08nfc, f08pec or f08njc) |
| |
| |
| Are all eigenvectors required? |
|
Is an upper Hessenberg matrix? |
|
f08pec or f08qkc |
| yes | yes |
| | no | | | no | |
| |
f08nac or f08nbc or (f08nhc, f08nec, f08nfc, f08pec, f08qkc or f08njc) |
| |
| |
| Is an upper Hessenberg matrix? |
|
f08pec or f08pkc |
| yes |
| | no | |
| f08nhc, f08nec, f08pec, f08pkc, f08ngc or f08njc |
|
Tree 4: Real Generalized Nonsymmetric Eigenvalue Problems
| Are eigenvalues only required? |
|
Are and in generalized upper Hessenberg form? |
|
f08xec |
| yes | yes |
| | no | | | no | |
| |
f08whc, (f08aec and f08agc or f08abc and f08acc), f08wec and f08xec |
| |
| |
| Is the generalized Schur factorization of and required? |
|
Are and in generalized upper Hessenberg form? |
|
f08xec |
| yes | yes |
| | no | | | no | |
| |
(f08aec and f08agc or f08abc and f08acc), f16qhc, f16qfc, f08afc, f08wec, f08xec and f08ykc |
| |
| |
| Are and in generalized upper Hessenberg form? |
|
f08xec and f08ykc |
| yes |
| | no | |
| f08whc, (f08aec and f08agc or f08abc and f08acc), f16qhc, f16qfc, f08afc, f08wec, f08xec, f08ykc and f08wjc |
|
Tree 5: Complex Hermitian Eigenvalue Problems
Tree 6: Complex Generalized Hermitian-definite Eigenvalue Problems
| Are eigenvalues only required? |
|
Are all eigenvalues required? |
|
Are and stored with one triangle as a linear array? |
|
f07grc, f08tsc, f08gsc and f08jfc |
| yes | yes | yes |
| | no | | | no | | | no | |
| |
|
f07frc, f08ssc, f08fsc and f08jfc |
| | |
| | |
| |
Are and stored with one triangle as a linear array? |
|
f07grc, f08tsc, f08gsc and f08jjc |
| | yes |
| | | no | |
| |
f07frc, f08ssc, f08gsc and f08jjc |
| |
| |
| Are all eigenvalues and eigenvectors required? |
|
Are and stored with one triangle as a linear array? |
|
f07grc, f08tsc, f08gsc, f08gtc and f16psc |
| yes | yes |
| | no | | | no | |
| |
f07frc, f08ssc, f08fsc, f08ftc, f08jsc and f16zjc |
| |
| |
| Are and stored with one triangle as a linear array? |
|
f07grc, f08tsc, f08gsc, f08jjc, f08jxc, f08guc and f16slc |
| yes |
| | no | |
| f07frc, f08ssc, f08fsc, f08jjc, f08jxc, f08fuc and f16zjc |
|
Tree 7: Complex non-Hermitian Eigenvalue Problems
| Are eigenvalues only required? |
|
Is an upper Hessenberg matrix? |
|
f08psc |
| yes | yes |
| | no | | | no | |
| |
f08nvc, f08nsc and f08psc |
| |
| |
| Is the Schur factorization of required? |
|
Is an upper Hessenberg matrix? |
|
f08psc |
| yes | yes |
| | no | | | no | |
| |
f08nsc, f08ntc, f08psc and f08nwc |
| |
| |
| Are all eigenvectors required? |
|
Is an upper Hessenberg matrix? |
|
f08psc and f08qxc |
| yes | yes |
| | no | | | no | |
| |
f08nvc, f08nsc, f08ntc, f08psc, f08qxc and f08nwc |
| |
| |
| Is an upper Hessenberg matrix? |
|
f08psc and f08pxc |
| yes |
| | no | |
| f08nvc, f08nsc, f08psc, f08pxc, f08nuc and f08nwc |
|
Tree 8: Complex Generalized non-Hermitian Eigenvalue Problems
| Are eigenvalues only required? |
|
Are and in generalized upper Hessenberg form? |
|
f08xsc |
| yes | yes |
| | no | | | no | |
| |
f08wvc, (f08asc and f08auc or f08apc and f08aqc), f08wsc and f08xsc |
| |
| |
| Is the generalized Schur factorization of and required? |
|
Are and in generalized upper Hessenberg form? |
|
f08xsc |
| yes | yes |
| | no | | | no | |
| |
(f08asc and f08auc or f08apc and f08aqc), f16thc, f16tfc, f08atc, f08wsc, f08xsc and f08yxc |
| |
| |
| Are and in generalized upper Hessenberg form? |
|
f08xsc and f08yxc |
| yes |
| | no | |
| f08wvc, (f08asc and f08auc or f08apc and f08aqc), f16thc, f16tfc, f08atc, f08wsc, f08xsc, f08yxc and f08wwc |
|
4.2 General Purpose Functions (singular value decomposition)
Tree 9
5 Functionality Index
| Backtransformation of eigenvectors from those of balanced forms, | | |
| Backtransformation of generalized eigenvectors from those of balanced forms, | | |
| Eigenvalue problems for condensed forms of matrices, | | |
| complex Hermitian matrix, | | |
| eigenvalues and eigenvectors, | | |
| all eigenvalues and eigenvectors by a divide-and-conquer algorithm, using packed storage | | nag_zhbevd (f08hqc) |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_zhbevx (f08hpc) |
| all eigenvalues and eigenvectors by a divide-and-conquer algorithm, using packed storage | | nag_zhpevd (f08gqc) |
| all eigenvalues and eigenvectors by root-free QR algorithm, using packed storage | | nag_zhpev (f08gnc) |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_zheevx (f08fpc) |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration, using packed storage | | nag_zhpevx (f08gpc) |
| all eigenvalues and eigenvectors using Relatively Robust Representations or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_zheevr (f08frc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm | | nag_zhbev (f08hnc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection | | nag_zhbevx (f08hpc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, using packed storage | | nag_zhbevd (f08hqc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm | | nag_zheev (f08fnc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection | | nag_zheevx (f08fpc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection, using packed storage | | nag_zhpevx (f08gpc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, using packed storage | | nag_zhpev (f08gnc) |
| complex upper Hessenberg matrix, reduced from complex general matrix, | | |
| singular value decomposition, | | |
| eigenvalues and eigenvectors, | | |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_dsbevx (f08hbc) |
| all eigenvalues and eigenvectors by a divide-and-conquer algorithm, using packed storage | | nag_dspevd (f08gcc) |
| all eigenvalues and eigenvectors by root-free QR algorithm, using packed storage | | nag_dspev (f08gac) |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_dsyevx (f08fbc) |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration, using packed storage | | nag_dspevx (f08gbc) |
| all eigenvalues and eigenvectors using Relatively Robust Representations or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_dsyevr (f08fdc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm | | nag_dsbev (f08hac) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection | | nag_dsbevx (f08hbc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm | | nag_dsyev (f08fac) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection | | nag_dsyevx (f08fbc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection, using packed storage | | nag_dspevx (f08gbc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, using packed storage | | nag_dspev (f08gac) |
| real symmetric tridiagonal matrix, | | |
| eigenvalues and eigenvectors, | | |
| after reduction from complex Hermitian matrix, | | |
| all eigenvalues and eigenvectors by root-free QR algorithm or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_dstevx (f08jbc) |
| all eigenvalues and eigenvectors using Relatively Robust Representations or selected eigenvalues and eigenvectors by bisection and inverse iteration | | nag_dstevr (f08jdc) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm | | nag_dstev (f08jac) |
| all eigenvalues by the Pal–Walker–Kahan variant of the QL or QR algorithm, or selected eigenvalues by bisection | | nag_dstevx (f08jbc) |
| real upper Hessenberg matrix, reduced from real general matrix, | | |
| Eigenvalue problems for nonsymmetric matrices, | | |
| all eigenvalues, Schur form, Schur vectors and reciprocal condition numbers | | nag_zgeesx (f08ppc) |
| all eigenvalues and left/right eigenvectors, plus balancing transformation and reciprocal condition numbers | | nag_zgeevx (f08npc) |
| all eigenvalues, real Schur form, Schur vectors and reciprocal condition numbers | | nag_dgeesx (f08pbc) |
| all eigenvalues and left/right eigenvectors, plus balancing transformation and reciprocal condition numbers | | nag_dgeevx (f08nbc) |
| Eigenvalues and generalized Schur factorization, | | |
| General Gauss–Markov linear model, | | |
| Generalized eigenvalue problems for condensed forms of matrices, | | |
| complex Hermitian-definite eigenproblems, | | |
| all eigenvalues and eigenvectors by a divide-and-conquer algorithm, packed storage format | | nag_zhpgvd (f08tqc) |
| all eigenvalues and eigenvectors by reduction to tridiagonal form, packed storage format | | nag_zhpgv (f08tnc) |
| selected eigenvalues and eigenvectors by reduction to tridiagonal form, packed storage format | | nag_zhpgvx (f08tpc) |
| real symmetric-definite eigenproblems, | | |
| all eigenvalues and eigenvectors by a divide-and-conquer algorithm, packed storage format | | nag_dspgvd (f08tcc) |
| all eigenvalues and eigenvectors by reduction to tridiagonal form, packed storage format | | nag_dspgv (f08tac) |
| selected eigenvalues and eigenvectors by reduction to tridiagonal form, packed storage format | | nag_dspgvx (f08tbc) |
| Generalized eigenvalue problems for nonsymmetric matrix pairs, | | |
| complex nonsymmetric matrix pairs, | | |
| all eigenvalues, generalized Schur form, Schur vectors and reciprocal condition numbers | | nag_zggesx (f08xpc) |
| all eigenvalues and left/right eigenvectors, plus the balancing transformation and reciprocal condition numbers | | nag_zggevx (f08wpc) |
| real nonsymmetric matrix pairs, | | |
| all eigenvalues, generalized real Schur form and left/right Schur vectors | | nag_dgges (f08xac) |
| all eigenvalues, generalized real Schur form and left/right Schur vectors, plus reciprocal condition numbers | | nag_dggesx (f08xbc) |
| all eigenvalues and left/right eigenvectors, plus the balancing transformation and reciprocal condition numbers | | nag_dggevx (f08wbc) |
| Generalized QR factorization, | | |
| Generalized RQ factorization, | | |
| Generalized singular value decomposition, | | |
| after reduction from complex general matrix, | | |
| after reduction from real general matrix, | | |
| reduction of a pair of general matrices to triangular or trapezoidal form, | | |
| minimum norm solution using the singular value decomposition (divide-and-conquer) | | nag_zgelsd (f08kqc) |
| minimum norm solution using the singular value decomposition (divide-and-conquer) | | nag_dgelsd (f08kcc) |
| least squares problems with linear equality constraints, | | |
| minimum norm solution subject to linear equality constraints using a generalized RQ factorization | | nag_zgglse (f08znc) |
| minimum norm solution subject to linear equality constraints using a generalized RQ factorization | | nag_dgglse (f08zac) |
| Left and right eigenvectors of a pair of matrices, | | |
| LQ factorization and related operations, | | |
| Operations on eigenvectors of a real symmetric or complex Hermitian matrix, or singular vectors of a general matrix, | | |
| Operations on generalized Schur factorization of a general matrix pair, | | |
| re-order Schur factorization, compute generalized eigenvalues and condition numbers | | nag_ztgsen (f08yuc) |
| re-order Schur factorization, compute generalized eigenvalues and condition numbers | | nag_dtgsen (f08ygc) |
| Operations on Schur factorization of a general matrix, | | |
| re-order Schur factorization, compute basis of invariant subspace, and estimate sensitivities | | nag_ztrsen (f08quc) |
| re-order Schur factorization, compute basis of invariant subspace, and estimate sensitivities | | nag_dtrsen (f08qgc) |
| Overdetermined and underdetermined linear systems, | | |
| Performs a reduction of eigenvalue problems to condensed forms, and related operations, | | |
| QL factorization and related operations, | | |
| QR factorization and related operations, | | |
| triangular-pentagonal matrices, | | |
| triangular-pentagonal matrices, | | |
| Reduction of a pair of general matrices to generalized upper Hessenberg form, | | |
| Reduction of eigenvalue problems to condensed forms, and related operations, | | |
| complex general matrix to upper Hessenberg form, | | |
| complex Hermitian matrix to real symmetric tridiagonal form, | | |
| complex rectangular matrix to real bidiagonal form, | | |
| real general matrix to upper Hessenberg form, | | |
| real rectangular matrix to bidiagonal form, | | |
| real symmetric matrix to symmetric tridiagonal form, | | |
| Reduction of generalized eigenproblems to standard eigenproblems, | | |
| complex Hermitian-definite generalized eigenproblem Ax = λBx, ABx = λx or BAx = λx | | nag_zhegst (f08ssc) |
| complex Hermitian-definite generalized eigenproblem Ax = λBx, ABx = λx or BAx = λx, packed storage | | nag_zhpgst (f08tsc) |
| real symmetric-definite generalized eigenproblem Ax = λBx, ABx = λx or BAx = λx | | nag_dsygst (f08sec) |
| real symmetric-definite generalized eigenproblem Ax = λBx, ABx = λx or BAx = λx, packed storage | | nag_dspgst (f08tec) |
| RQ factorization and related operations, | | |
| Singular value decomposition, | | |
| preconditioned Jacobi SVD using fast scaled rotations and de Rijks pivoting | | nag_dgejsv (f08khc) |
| Solve generalized Sylvester equation, | | |
| Solve reduced form of Sylvester matrix equation, | | |
| Split Cholesky factorization, | | |
6 Auxiliary Functions Associated with Library Function Arguments
None.
7 Functions Withdrawn or Scheduled for Withdrawal
None.
8 References
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 (2012) Matrix Computations (4th 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