Direct Methods

In this section, we briefly discuss direct methods for the computation of eigenvalues of symmetric matrices that can be stored in the computer as full matrices. These direct methods are sometimes called transformation methods and are built up around similarity transformations.

They are implemented in LAPACK [12], ScaLAPACK [52],
and the `eig` command in
MATLAB.^{}

All the subspace-based algorithms we are going to discuss later in this chapter need to use dense, tridiagonal, or banded matrix routines as inner iterations to get Ritz approximations to subspace eigenvalues.

The most common direct methods have two phases:

- Find an orthogonal matrix such that is a tridiagonal matrix.
- Compute the eigendecomposition of the tridiagonal matrix .

The initial reduction to tridiagonal form is made by a sequence of
orthogonal Householder reflections and takes
floating point operations,
or
if eigenvectors are also desired.
The symmetric or Hermitian driver routines of LAPACK start with this
reduction: it is computed by the computational routines `xSYTRD`
for the real and `xHETRD` for the complex case (`PxSYTRD`
and `PxHETRD` in ScaLAPACK).
There are also implementations for packed and banded storage.

There is a choice of methods for finding the eigendecomposition of a tridiagonal matrix:

*QR algorithm:*- This algorithm finds all the eigenvalues
and optionally all the eigenvectors.
It takes floating point operations for finding all the eigenvalues
of a tridiagonal matrix. Since
reducing a dense matrix to tridiagonal form costs
floating point operations, is negligible for large enough .
For finding all the eigenvectors as well, QR iteration takes a little
over floating point operations on average.
It is implemented as
`xSTEQR`and`xSTERF`in LAPACK.This is the algorithm underlying the MATLAB command

`eig`.^{} *Divide-and-conquer method:*-
It divides the tridiagonal matrix into two halves, solves the eigenproblems
of each of the halves, and glues the solutions together
by solving a special rational equation. It is implemented in LAPACK
as
`xSTEVD`.`xSTEVD`can be many times faster than`xSTEQR`for large matrices but needs more workspace ( or ). *Bisection method and inverse iteration:*-
Bisection may be used to find just a subset of the
eigenvalues, say those in an interval .
It needs only floating point operations, where is the
number of eigenvalues desired.
Thus the bisection method could be much faster than
the QR method when .
It can be highly accurate, but may be adjusted to run faster
if lower accuracy is acceptable.
Inverse iteration can then be used to find the corresponding eigenvectors. In the best case, when the eigenvalues are well separated, inverse iteration also costs only floating point operations. This is much less than either QR or divide-and-conquer, even when all eigenvalues and eigenvectors are desired (). On the

other hand, when many eigenvalues are clustered close together, Gram-Schmidt orthogonalization will be needed to make sure that we do not get several identical eigenvectors. This will add floating point operations to the operation count in the worst case.Bisection method and inverse iteration are implemented in the LAPACK routines

`xSTEBZ`and`xSTEIN`, and they are available as options in`xSYEVX`. In ScaLAPACK, they are subroutines`PxSTEBZ`and`PxSTEIN`. *Relatively robust representation algorithm:*-
This algorithm uses an LDL factorization of a number of translates
of , for one shift near each cluster of eigenvalues.
For each translate the algorithm computes very accurate eigenpairs
for the tiny eigenvalues.
It is implemented as
`xSTEGR`in LAPACK.`xSTEGR`is faster than all the routines except in a few cases and uses the least workspace. See [358,128,360].

Finally, a classical method for solving Hermitian eigenvalue problems
is the *Jacobi method*.
This method constructs an orthogonal transformation to diagonal form,

by applying a sequence of elementary orthogonal rotations, each time reducing the sum of squares of the nondiagonal elements of the matrix, until it is of diagonal form to working accuracy.

The Jacobi algorithm has been very popular, since it is implemented by a very simple program and gives eigenvectors that are orthogonal to working accuracy. However, it cannot compete with the QR method in terms of operation counts: Jacobi needs multiplications for sweeps ( to 5 usually), which is more than the needed for tridiagonal reduction.

There is one important advantage to the Jacobi algorithm. Properly implemented it can deliver eigenvalue approximations with a small error in the relative sense, in contrast to algorithms based on tridiagonalization, which only guarantee that the error is bounded relative to the norm of the matrix (4.3). See [124].