...matlab,

MATLAB is a registered trademark of The MathWorks, Inc.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... Software)

http://gams.nist.gov, developed by NIST (National Institute of Standards and Technology).

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... NETLIB.

http://www.netlib.org

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... constant.

The word pencil arises historically, because the set of all matrices of the form $A - \lambda B$ for constant $A$, constant $B$, and varying $\lambda$ forms a ``line'' of matrices in matrix space, and this resembles a ``pencil'' or beam of light.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... form,

It is usually called staircase form in the literature, but we find Jordan-Schur more appropriate.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... random,''

For example, choose each entry independently and randomly from $(-1,1)$ or any other open interval of real numbers.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...

In the case of Hermitian eigenproblems, we mentioned the ability to count the number of eigenvalues in an interval $[\alpha, \beta]$. This can sometimes be done for sparse Hermitian matrices much more cheaply than computing the eigenvalues in $[\alpha, \beta]$. Although such counting algorithms exist for the NHEP [32], their costs are similar to transformation methods that actually compute the eigenvalues, so we do not pursue them.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

....

For example, $p(\lambda)$ is identically zero for the 1-by-1 pencil $A - \lambda B = 0 - \lambda \cdot 0$.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... form,

It is usually called staircase form in the literature, but we find Weierstrass-Schur more appropriate.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

....

To see this in the generic case when all the eigenvalues are distinct, multiply $AX = BX \Lambda$ by $Y^*$ and $Y^*A = \Lambda Y^* B$ by $X$ to get $\Lambda Y^*BX = Y^*AX = Y^*BX \Lambda$. $\Lambda Y^*BX = Y^*BX \Lambda$ implies that $Y^*BX = \Lambda_B$ must be diagonal, whence $Y^*AX = \Lambda_A = \Lambda_B \Lambda$ is also diagonal.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... random,''

For example, choose each entry independently and randomly from $(-1,1)$ or any other open interval of real numbers.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... eigenvector.

These blocks look like $I - \lambda J(0)$, where $J(0)$ is a Jordan block with 0 eigenvalue.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... users.

In the case of GHEPs, we mentioned the ability to count the number of eigenvalues in an interval $[\alpha, \beta]$. This can sometimes be done for sparse Hermitian pencils much more cheaply than computing the eigenvalues in $[\alpha, \beta]$. Although such counting algorithms exist for the regular GNHEP [32], their costs are similar to transformation methods that actually compute the eigenvalues, so we do not pursue them.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... out.

Note that ``direct'' methods must still iterate, since finding eigenvalues is mathematically equivalent to finding zeros of polynomials, for which no noniterative methods can exist. We can call a method direct if experience shows that it (nearly) never fails to converge in a fixed number of iterations.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... MATLAB.

It has been announced that an LAPACK-based numerics library will be part of the next major release of MATLAB; see The Newsletter for MATLAB, Simulink and Toolbox Users, Winter 2000, available at http://www.mathworks.com/company/newsletter/clevescorner/winter2000.cleve.shtml

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...eig.

MATLAB checks to see whether the argument of eig is symmetric or not and uses the symmetric QR algorithm when appropriate.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

....

There is now a much better algorithm for constructing and applying these transformations; see D. Sorensen, Deflation for Implicitly Restarted Arnoldi Methods, CAAM Technical Report, TR 98-12, Rice University, 1998 (revised August 2000).

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... matrix

Deflation with approximate eigenvectors may introduce an error of order $\epsilon^2$ on the eigenvalues, provided that the computed eigenvalues are well separated from the remaining ones [353, Section 5.1].

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...

When $B$ is ill-conditioned, incorporating pivoting into the Cholesky decomposition can improve the numerical stability of the process [472].

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... large.

Exactly how large is too large is an ambiguous notion and should generally be dealt with on a case-by-case basis. Usually one may regard any number over 1000 as too large for condition numbers.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...kapj82,sun98,hihi98.

Kahan, Parlett, and Jiang [256, Theorem 1$^{\prime}$] gave only optimal norms for $\wtd\beta E-\wtd\alpha F\equiv Z$. Take $E=\wtd\beta Z$ and $F=-\wtd\alpha Z$ to see that $(E,F)$ has the same 2-norm and Frobenius norm as $Z$.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... MATLAB.

It has been announced that an LAPACK-based numerics library will be part of the next major release of MATLAB; see The Newsletter for MATLAB, Simulink and Toolbox Users, Winter 2000, available at http://www.mathworks.com/company/newsletter/clevescorner/winter2000.cleve.shtml

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... called.

This code makes very mild assumptions about floating point arithmetic. It will work on machines with a guard digit in add/subtract or on those binary machines without guard digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably fail on hexadecimal or decimal machines without guard digits, but we know of none.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... method

Note that a direct method must still iterate, since finding eigenvalues is mathematically equivalent to finding zeros of polynomials, for which no noniterative methods can exist. We call a method direct if experience shows that it (nearly) never fails to converge in a fixed number of iterations.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... MATLAB.

It has been announced that an LAPACK-based numerics library will be part of the next major relase of MATLAB; see The Newsletter for MATLAB, Simulink and Toolbox Users, Winter $2000$, available at http://www.mathworks.com/company/newsletter/clevescorner/winter2000.cleve.shtml

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

....

There is now a much better algorithm for constructing and applying these transformations; see D. Sorensen, Deflation for Implictly Restarted Arnoldi Methods, CAAM Technical Report, TR 98-12, Rice University, 1998 (revised August 2000).

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...lanc50.

Historical note: It appears that Lanczos did not know Krylov's work. Ostrowski informed Lanczos about the early parallel work of Krylov.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...

A version of the algorithm that enhances numerical stability at the cost of some extra inner products is described in R. W. Freund, Using Krylov subspace methods with inexact deflation for reduced-order model. Bell Labs Numerical Analysis Manuscript, August 2000.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... MATLAB.

It has been announced that an LAPACK-based numerics library will be part of the next major release of MATLAB; see The Newsletter for MATLAB, Simulink and Toolbox Users, Winter $2000$ available at http://www.mathworks.com/company/newsletter/clevescorner/winter2000.cleve.shtml

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... product.

The $B$ inner product between two vectors $x$ and $y$ is defined as $(x,y)_B = y^{\ast} B x$. If $B$ is symmetric and indefinite, then $(x,y)_B$ is a pseudo-inner product. It violates the condition $(x,x)_B \geq 0$ for all $x$ in the definition of a standard inner product.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... problem:

The term ``linear'' is in quotation marks because $\lambda$ appears in the eigenvectors as well.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

...equivalent

Two matrix polynomials $M_1(\lambda)$ and $M_2(\lambda)$ of size $n \times n$ are called equivalent if $M_1(\lambda) = E(\lambda)M_2(\lambda)F(\lambda)$ for some $n \times n$ matrix polynomials $E(\lambda)$ and $F(\lambda)$ with constant nonzero determinants (unimodular).

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... generalization,

See the book's homepage, ETHOME, for access to the software described in this section.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

... vectors,

Over complex numbers, we always take the real part of $\tr(A^H B)$.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.