HOMOTOPY METHOD FOR THE SINGULAR SYMMETRIC TRIDIAGONAL EIGENPROBLEM

Kuiyuan Li

University of West Florida

Tien-Yien Li

Michigan State University

Abstract. In this paper, a homotopy algorithm for finding some or all finite eigenvalues and corresponding eigenvectors of a real symmetric matrix pencil (A, B) is presented, where A is a symmetric tridiagonal matrix and B is a diagonal matrix with b_{i} ≥ 0, i = 1, 2, …, n. It is shown that there are exactly m (m is the number of finite eigenvalues of (A, B)) disjoint, smooth homotopy paths connecting the trivial eigenpairs to the desired eigenpairs. And the eigenvalue curves are monotonic and easy to follow. The performance of the parallel version of our algorithm is presented.

1. Introduction. Consider the generalized real symmetric eigenvalue problem:

(1)  Ax = λBx

where A, B are real symmetric and B is positive semidefinite. By MDR reduction [1], A can be reduced to a symmetric tridiagonal matrix and B to a positive semidefinite diagonal matrix simultaneously. Hence, from now on, we will assume A is a real symmetric tridiagonal matrix and B is a diagonal matrix with b_{i} ≥ 0. Then, (1) is called singular symmetric tridiagonal eigenproblem.

If all b_{i}’s are positive and well conditioned, then (1) can be reduced to a standard symmetric tridiagonal eigenproblem. There are many efficient algorithms for that problem. However, if B is singular, this technique is not available. To the best of our knowledge, there are few efficient algorithms for solving (1) but Fix Heiberger [3], Bunse Gerstner [1] and the QZ method [4]. However, the complexity of these methods is O(n^{3}).

In this paper, we will present a new homotopy method for finding some or all finite eigenpairs of (1).

Assume in (1)

(2)  A = ( \begin{matrix}α_{1} & β_{2} \\ β_{2} & α_{2} & β_{3} \\ & ... & ... & ... \\ & & β_{n−1} & α_{n−1} & β_{n} \\ & & & β_{n} & α_{n} \end{matrix} )

and B = diag(b_{1}, b_{2}, …, b_{n}) with b_{i} ≥ 0.

If β_{i} = 0 for some i, 2 ≤ i ≤ n, then ℝ^{n} can clearly be decomposed into two complementary subspaces invariant under A. Thus the generalized eigenproblem Ax = λBx is decomposed in an obvious way into two smaller subproblems. Hence, we will assume that each β_{i} ≠ 0. That is, A is unreduced.

Let D be an n × n symmetric tridiagonal matrix and consider the homotopy H: ℝ^{n} × ℝ × [0, 1] → ℝ^{n} × ℝ, defined by

(3)  \begin{aligned}H(x, λ, t) & = (1 − t) ( λBx − Dx_{\frac{x^{T} Bx−1}{2}} ) + t ( λBx − Ax_{\frac{x^{T} Bx−1}{2}} ) \\ & = ( λBx − [(1 − t)D + tA]x_{\frac{x^{T} Bx−1}{2}} ) \\ & = ( λBx − A(t)x_{\frac{x^{T} Bx−1}{2}} ) \end{aligned}

where A(t) = (1 − t)D + tA. The pencil (D, B) is called an initial pencil.

In section 2, we will show that the solution set of H(x, λ, t) = 0 in (3) consists of exactly m (m is the number of finite eigenvalues of (A, B)) disjoint smooth curves (x(t), λ(t)), each joins an eigenpair of (D, B) to one of (A, B). We call each of these curves a homotopy curve or an eigenpath and its component λ(t) an eigenvalue curve. We will also show that each eigenvalue curve is monotonic in t. In section 3, we will give the details of our algorithm. In section 4, some numerical results will be presented.

The homotopy method may become an efficient method for this problem since it can be used to find some or all finite eigenpairs without any waste on computing the infinite eigenvalues and its complexity is O(n^{2}). It is also a parallel scheme since the homotopy curves can be followed independently.

2. Preliminary Analysis. If

(4) B = ( \begin{matrix}B_{1} \\ & O_{1} \\ & & B_{2} \\ & & & . & . & . \\ & & & & & & B_{r} \end{matrix} ) , we let A = ( \begin{matrix}A_{1} & ∗ \\ ∗ & A_{1}^{0} & ∗ \\ & ... & ... & ... \\ & & & A_{r−1}^{0} & ∗ \\ & & & ∗ & A_{r} \end{matrix} ) ,

where O_{i}’s are zero matrices with dim(O_{i}) = dim(A_{i}^{0}), and B_{i}’s are positive definite diagonal matrices with dim(B_{i}) = dim(A_{i}), i = 1, 2, …, r.

\underline{Theorem 2.1}. Let n(A, B) denote the number of the eigenvalues of the pencil (A, B), then

n(A, B) = rank(B) − ∑_{i=1}^{r}(1 − sign(| det(A_{i}^{0})|)).

\underline{Proof}. See [6].

Choose k such that β_{k+1} is an off–diagonal element of A_{i} for some i or the off–diagonal element of A which joins block matrices A_{i−1}^{0} and A_{i}, then let

(5)  D = ( \begin{matrix}D_{1} & 0 \\ 0 & D_{2} \end{matrix} ) ,

where

D_{1} = ( \begin{matrix}α_{1} & β_{2} \\ β_{2} & α_{2} & β_{3} \\ & ... & ... & ... \\ & & β_{k−1} & α_{k−1} & β_{k} \\ & & & β_{k} & α_{k} \end{matrix} ) , D_{2} = ( \begin{matrix}α_{k+1} & β_{k+2} \\ β_{k+2} & α_{k+2} & β_{k+3} \\ & ... & ... & ... \\ & & β_{n−1} & α_{n−1} & β_{n} \\ & & & β_{n} & α_{n} \end{matrix} ) ,

and then let

B = ( \begin{matrix}B_{1} & 0 \\ 0 & B_{2} \end{matrix} )

with dim(B_{i}) = dim(D_{i}), i = 1, 2.

Let A(t) and B be the same as in (3), then we get the following corollary immediately. \underline{Corollary 2.1}.

n(A(t), B) = rank(B) − ∑_{i=1}^{r}(1 − sign(| det(A_{i}^{0})|)) = n(A, B) for t ∈ (0, 1].

n(D_{1}, B_{1}) + n(D_{2}, B_{2}) = n(A, B).

Therefore, there are exactly m = n(A, B) finite homotopy curves.

\underline{Theorem 2.2}. Let H and D be given in (3) and (5), then the homotopy paths are distinct and smooth.

\underline{Proof}. Differentiating H(x, λ, t) = 0 with respect to t, we get

(6)  H_{x}ẋ + H_{λ}λ̇ + H_{t} = 0.

By the simple computation, we have

(7)  ( \begin{matrix}λB − A(t) & Bx \\ x^{T} B & 0 \end{matrix} ) ( ẋ ; λ̇ ) = ( \begin{matrix}(A − D)x \\ 0 \end{matrix} ) .

Now, we show

(8)  H_{(x,λ)}(x, λ, t) = ( \begin{matrix}λB − A(t) & Bx \\ x^{T} B & 0 \end{matrix} )

is nonsingular if (λ(t), x(t)) is an eigenpair of (A(t), B).

Since (λ(t)B−A(t)) is an unreduced tridiagonal matrix for any t in (0, 1] and (λ(t), x(t)) is one of its eigenpairs, dim(ker(λ(t)B − A(t))) = 1. If H_{(x,λ)}(x, λ, t)y = 0 for some y, then

(9)  \{ \begin{matrix}(λ(t)B − A(t))y_{1} + Bxy_{2} = 0 \\ x^{T} By_{1} = 0 \end{matrix}

where y = (y_{1}, y_{2})^{T}.

Clearly, x(t)^{T} Bx(t) ≠ 0. Otherwise, x_{1} = 0 since b_{1} > 0 and b_{i} ≥ 0. Hence, (λ(t)B − A(t))x = 0 implies that the last n − 1 columns of λ(t)B − A(t) are linearly dependent. This contradicts to the fact that λ(t)B − A(t) is an unreduced tridiagonal matrix.

Since x(t)^{T} Bx(t) ≠ 0 and x^{T} (t)(λ(t)B − A(t)) = 0, y_{2} = 0 from (9). Hence, y_{1} ∈ ker(λ(t)B − A(t)). Since dim(ker(λ(t)B − A(t))) = 1 and x^{T} By_{1} = 0, y_{1} = 0. Hence, y = 0, i.e., H_{(x,λ)}(x, λ, t) is nonsingular.

By the nonsingularity of H_{(x,λ)}(x, λ, t) and the implicit function theorem, the homotopy curves are distinct and smooth.

It follows from Corollary 2.1 and Theorem 2.2 that none of the homotopy paths come from infinity at t = 0 and diverges to infinity at t = 1.

Let (A)^{1} denote the lower (n − 1) × (n − 1) submatrix of A, (A)_{1} denote the upper (n − 1) × (n − 1) submatrix of A.

\underline{Lemma 2.1}. The eigenvalues of (A, B) strictly separate the eigenvalues of ((A)_{1}, (B)_{1}), where A is an unreduced symmetric tridiagonal matrix and B is a diagonal matrix with b_{1} > 0 and b_{i} ≥ 0.

\underline{Proof}. See [15].

\underline{Theorem 2.3}. The eigenvalue curve λ(t) is either constant or strictly monotonic. Furthermore, if all eigenvalues of the initial pencil (D, B) are distinct, then

(i) λ̇(t)λ̈(t) > 0 for t small, if λ̇(t) ≢ 0.

(ii) λ(t) is bounded by two consecutive eigenvalues of (D, B).

\underline{Proof}. Let f_{1} = det(D_{1} − λB_{1}), f_{2} = det(D_{2} − λB_{2}), f_{3} = det((D_{1})_{1} − λ(B_{1})_{1}) and f_{4} = det((D_{2})^{1} − λ(B_{2})^{1}). Since

det(A(t) − λ(t)B) = det ( \begin{matrix}D_{1} − λ(t)B_{1} & tβ_{k+1} \\ tβ_{k+1} & D_{2} − λ(t)B_{2} \end{matrix} ) = 0,

we have

(10)  f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t)) = 0.

If there exists a t_{0} in [0, 1] for which f_{3}(λ(t_{0}))f_{4}(λ(t_{0})) = 0 then either f_{3}(λ(t_{0})) = 0 or f_{4}(λ(t_{0})) = 0; say f_{3}(λ(t_{0})) = 0. It follows from Lemma 2.1 that f_{1}(λ(t_{0})) ≠ 0. Hence, f_{2}(λ(t_{0})) = 0 in (10); accordingly,

f_{1}(λ(t_{0}))f_{2}(λ(t_{0})) − t^{2}β_{k+1}^{2}f_{3}(λ(t_{0}))f_{4}(λ(t_{0})) ≡ 0.

This implies det(A(t) − λ(t_{0})B) = 0. Thus, λ(t) = λ(t_{0}) for all t in [0, 1].

Assume f_{3}(λ(t))f_{4}(λ(t)) ≠ 0 for any t in [0, 1]. Write λ̇(t) = \frac{d}{dt}λ(t). Differentiating (10) with respect to t yields,

(11)  \frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))]λ̇(t) = 2tβ_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))

so,

\frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))] ≠ 0 for any t ∈ (0, 1].

Hence,

λ̇(t) = \frac{2tβ_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))}{\frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))]} ≠ 0 for any t ∈ (0, 1].

Therefore, the eigenvalue curve is strictly monotonic.

Furthermore, if all the eigenvalues of (D, B) are distinct, then we claim that

(12)  \frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))]|_{t=0} ≠ 0.

Otherwise,

\begin{matrix}\frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))]|_{t=0} = \frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t))]|_{t=0} \\ = \frac{d}{dλ}f_{1}(λ(t))|_{t=0}f_{2}(λ(0)) + f_{1}(λ(0))\frac{d}{dλ}f_{2}(λ(t))|_{t=0} = 0. \end{matrix}

Since λ(0) is an eigenvalue of (D, B), we have f_{1}(λ(0))f_{2}(λ(0)) = 0. If f_{1}(λ(0)) = 0 then f_{2}(λ(0)) ≠ 0 since all eigenvalues of (D, B) are distinct. Hence, \frac{d}{dλ}f_{1}(λ(t))|_{t=0} = 0. Consequently, λ(0) is a multiple eigenvalue of (D_{1}, B_{1}) which contradicts to the fact that D_{1} is unreduced. The proof of (12) for the case f_{2}(λ(0)) = 0 follows by the same argument.

Hence,

λ̇(t) = \frac{2tβ_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))}{\frac{d}{dλ}[f_{1}(λ(t))f_{2}(λ(t)) − t^{2}β_{k+1}^{2}f_{3}(λ(t))f_{4}(λ(t))]}

is well defined for t in [0, 1] and λ̇(0) = 0.

Now, let g = f_{1}f_{2}, f = β_{k+1}^{2}f_{3}f_{4} and h = (g − t^{2}f)_{λ} = h(t, λ(t)), then λ̇ = 2tf/h

(13)  λ̈(t) = −4t^{2} \frac{f^{2}(g − t^{2}f)λλ}{h^{3}} + 8t^{2} \frac{f_{λ}f}{h^{2}} + \frac{2f}{h} .

Since f^{2}(g −t^{2}f)λλ/h^{3} and f_{λ}f/h^{2} are continuous in [0, 1], they are bounded. By (13),

λ̈(0) = lim_{t→0} λ̈(t) = lim_{t→0} 2f/h.

Hence, λ̇(t)λ̈(t) > 0 for small t, since λ̇(t) = 2tf/h.

Now we show that each eigenvalue curve λ_{i}(t) is bounded by two consecutive eigenvalues of (D, B). If there exists a t_{0} ≠ 0 such that λ_{i}(t_{0}) = λ_{j}(0) for some j, then

f_{1}(λ_{i}(t_{0}))f_{2}(λ_{i}(t_{0})) − t_{0}^{2}β_{k+1}^{2}f_{3}(λ_{i}(t_{0}))f_{4}(λ_{i}(t_{0})) = 0

from (10). Since λ_{i}(t_{0}) = λ_{j}(0) for some j, f_{1}(λ_{i}(t_{0}))f_{2}(λ_{i}(t_{0})) = 0. Therefore, f_{3}(λ_{i}(t_{0}))f_{4}(λ_{i}(t_{0})) = 0, since t_{0}^{2}β_{k+1}^{2} ≠ 0. Hence, λ_{i}(t) is constant.

From Theorem 2.2 and Theorem 2.3, every eigenvalue curve must be one of those in Figure 1. These results provide very important information in designing our code.

Figure 1.

3. Algorithm. The algorithm to follow the eigenpath (x(t), λ(t)) has the basic features given below.

(i) Initiating at t = 0 (ii) Prediction (iii) Correction (iv) Step–size selection

(v) Terminating at t = 1. In this section, we give a detailed description of these features.

(i) Initiating at t = 0. Choose the initial pencil (D, B) as mentioned in the last section and also try to make

the sizes of the blocks D_{1} and D_{2} roughly the same.

When the initial matrix D is decided, we calculate the eigenvalues of (D_{1}, B_{1}) and (D_{2}, B_{2}) by the multisection method.

(ii) Prediction.

Assume that after i steps the approximate value (x̃(t_{i}), λ̃(t_{i})) on the eigenpath (x(t), λ(t)) at t_{i} is known and the next step–size h is determined; that is, t_{i+1} = t_{i} + h. We want to find an approximate value (x̃(t_{i+1}), λ̃(t_{i+1})) of (x(t_{i+1}), λ(t_{i+1})) on the eigenpath at t_{i+1}. Notice that (x̃(t_{i+1}), λ̃(t_{i+1})) is an approximate eigenpair of (A(t_{i+1}), B).

We use the Euler predictor to predict the eigenvalue at t_{i+1}, namely,

λ^{0}(t_{i+1}) = λ(t_{i}) + λ̇(t_{i})h.

When i > 1, we know precisely which two consecutive eigenvalues of the initial pencil (D, B) bound λ(t). Since λ(t) is strictly monotonic, λ(t_{i}) must be either a lower bound or an upper bound of λ(t). If λ^{0}(t_{i+1}) is not in that interval, we let λ^{0}(t_{i+1}) be the middle point of the interval. To predict the eigenvector, we use the inverse power method of A(t_{i+1}) on x(t_{i}) with shift λ^{0}(t_{i+1}). At t_{i} = 0, since we skip the calculations of eigenvectors of (D, B), x(0) is not available. We choose a random vector to substitute for x(0).

(iii) Correction.

Simple computation shows that Newton’s method for the nonlinear problem of n + 1 equations

(14)  F (λ, x) = \{ \begin{matrix}λBx − Ax = 0 \\ \frac{x^{T} Bx−1}{2} = 0 \end{matrix}

in the n + 1 variables λ, x_{1}, x_{2}, …, x_{n} at (λ^{(k)}, x^{(k)}) is the inverse iteration,

(15)  (A − λ^{(k)}B)y = x^{(k)}

and

(16)  \{ \begin{matrix}λ^{(k+1)} = λ^{(k)} + \frac{(x^{(k)})^{T} x^{(k)}+1}{2(x^{(k)})^{T} y} \\ x^{(k+1)} = (λ^{(k+1)} − λ^{(k)})y. \end{matrix}

Newton’s method is used as a corrector. The stop point (x^{j}(t_{i+1}), λ^{j}(t_{i+1})) of Newton’s iteration will be taken as an approximate eigenpair (x̃(t_{i+1}), λ(t_{i+1})) of A(t_{i+1}). As mentioned in (ii), we know the upper and lower bounds of λ(t). When (x̃(t_{i+1}), λ̃(t_{i+1})) is taken as an approximate eigenpair of (A(t_{i+1}), B), we check if λ̃(t_{i+1}) is still in that interval. If the checking fails, we reduce the step size to h/2 and repeat the whole process once again beginning with the eigenvalue prediction in (ii).

(iv) Step–size selection.

In the first attempt, we always choose step–size h = 1 − t_{i} at t_{i} < 1. If after the prediction and correction steps the checking step fails, we reduce the step size to h/2 as mentioned in (iii). This extremely liberal choice of step–size can be justified because of the closeness of the matrix D to A as well as the effective checking algorithm. Indeed, since the initial matrix D is chosen to be so close to A, from our experiences, the majority of the eigenpairs of A can be reached in one step, i.e., h = 1.

Very small step–size can also cause the inefficiency of the algorithm. Therefore, we impose a minimum γ on step–size h. If h < γ, we simply give up following the eigenpath and the corresponding eigenpair of A will be calculated at the end of the algorithm by the method of bisection with inverse iterations (see (iv)). We usually choose γ ≈ 0.25.

(v) Terminating at t = 1.

At t = 1, when an approximate eigenvalue λ̃(1) is reached, we compute the generalized Sturm sequence at λ̃(1) ± ϵλ̃(1) with ϵ = machine precision to ensure the correct order. If the checking fails, we have jumped into a wrong eigenpath. More precisely, suppose we are following the ith eigenpair, the checking algorithm detects that we have reached the jth eigenpair instead. In this situation, we will save the jth eigenpair before the step–size is cut. By saving the jth eigenpair, the computation of following the jth eigenpair is no longer needed.

As mentioned in (iv), we may give up following some eigenpaths to avoid adopting a step–size that is too small. Without extra computation, we know exactly which eigenpairs are lost at t = 1. In order to find these eigenpairs, we first use the bisection to find the missing eigenvalues and then use inverse iteration to find the eigenvectors.

4. Numerical Results. In this section, we present our numerical results. Our homotopy continuation algorithm is in its preliminary stage, and much development and testing are necessary. But the numerical results on the examples we have looked at seem remarkable. Our testing examples are:

Type 1. A is an unreduced symmetric tridiagonal matrix with both diagonal and offdiagonal elements being uniformly distributed random numbers between 0 and 1. B is a diagonal matrix with the first n/2 diagonal elements being uniformly distributed random numbers between 0 and 1, and the last n/2 being zeros.

Type 2. A is an unreduced symmetric tridiagonal matrix with both diagonal and offdiagonal elements being uniformly distributed random numbers between 0 and 1. B is a diagonal matrix with the first 3n/10 and the last 3n/10 diagonal elements being uniformly distributed random numbers between 0 and 1, and the rest being zeros.

Type 3. A is the Toeplitz matrix [1,2,1]. B is a diagonal matrix with the first n/2 diagonal elements being 1, and the rest being zeros.

Type 4. A is the Toeplitz matrix [1,2,1]. B is a diagonal matrix with the first 3n/10 and the last 3n/10 diagonal elements being 1, and the rest being zeros.

Type 5. A is an unreduced symmetric tridiagonal matrix with both diagonal and offdiagonal elements being uniformly distributed random numbers between 0 and 1. B is a diagonal matrix with all diagonal elements being random numbers between 0 and 1.

Table 1: Execution time (second), speed-up and efficiency of STH.

The homotopy algorithm is to a large degree parallel since each eigenpath can be followed independently. This inherent nature of the homotopy method makes the parallel implementation much simpler than other methods.

In our parallel algorithm, after all the eigenvalues of D are computed and put in increasing order, we assign each processor to trace roughly m/p eigencurves, where m = n(A, B) = n(D_{1}, B_{1}) + n(D_{2}, B_{2}) and p is the number of processors being used. Let the first processor trace the first m/p smallest eigencurves from the smallest to the largest and let the second processor trace the second m/p smallest eigencurves, and so on.

We present the numerical results of the parallel implementation of our algorithm. All examples were executed on BUTTERFLY GP 1000, a shared memory multiprocessor machine.

The speed–up is defined as

S_{p} = \frac{execution time using one processor}{execution time using p processors}

and the efficiency is the ratio of the speed–up over p.

Table 1 shows the execution time and the speed–up S_{p} as well as the efficiency S_{p}/p of our algorithm STH on all five type matrices. The numerical result shows the homotopy method may become an excellent candidate for a variety of advanced architectures.

\underline{References}

1. A. Bunse–Gerstner, “An Algorithm for the Symmetric Generalized Eigenvalue Problem,” Linear Alg. Appl., 58 (1984), 43–68.

2. M. T. Chu, T. Y. Li, and T. Sauer, “Homotopy Method for General λ–Matrix Problem,” SIAM J. Matrix Anal. and Appl., 9 (1988), 528–536.

3. G. Fix and R. Heiberger, “An Algorithm for the Ill–Conditioned Generalized Eigenvalue Problem,” SIAM J. Numer. Analy., 9 (1972), 78–88.

4. G. Golub and C. V. Loan, Matrix Computations, The Johns Hopkins Press, Baltimore, MD, 1983.

5. B. S. Garbow, et.al., Matrix Eigensystem Routines - EISPACK Guide Extension, New York, Springer–Verlag, 1977.

6. K. Li, Homotopy Methods and Algorithms for Real Symmetric Eigenproblems, Ph.D. thesis, Michigan State University, 1991.

7. K. Li and T. Y. Li, “An Algorithm for Symmetric Tridiagonal Eigen–problem – Divide and Conquer with Homotopy Continuation,” SIAM J. Sci. Comput., 14 (1993), 735– 751.

8. T. Y. Li and N. H. Rhee, “Homotopy Algorithm for Symmetric Eigenvalue Problems,” Numer. Math., 55 (1989), 256–280.

9. T. Y. Li and T. Sauer, “Homotopy Method for Generalized Eigenvalue Problems,” Lin. Alg. Appl., 91 (1987), 65–74.

10. T. Y. Li, T. Sauer and J. Yorke, “Numerical Solution of a Class of Deficient Polynomial Systems,” SIAM J. Numer. Anal., 24 (1987), 435–451.

11. T. Y. Li, Z. Zeng, and L. Cong, “Solving Eigenvalue Problems of Real Nonsymmetric Matrices with Real Homotopy,” SIAM Numer. Anal., to appear.

12. B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice–Hall, Englewood Cliffs, N.J., 1980.

13. G. Peters and J. H. Wilkinson, “Ax = λBx and the Generalized Eigenproblem,” SIAM J. Numer. Anal., 7 (1970), 479–492.

14. J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, New York, 1965.

15. Z. Zhang, “Generalized Sturm Sequences,” preprint.