DURAND–KERNER ROOT–FINDING METHOD FOR THE GENERALIZED TRIDIAGONAL EIGENPROBLEM

Kuiyuan Li

Abstract. In this paper we present a method, parallel in nature, for finding all eigenvalues of a symmetric generalized tridiagonal matrix. Our method employs the determinant evaluation and the Durand–Kerner root–finding scheme. It will be shown that the method converges quadratically and is reliable, efficient, and easy to implement in practice.

1. Introduction. In this paper we consider the generalized eigenvalue problem

Ax = λMx,  (1)

where A and M are both real, symmetric, and tridiagonal and one of them, say M, is positive definite. These assumptions imply real eigenvalues and the existence of an M — orthogonal basis of eigenvectors. Eigenvalue problems with such a special structure arise in many applications, such as finding numerical solutions of the Sturm–Liouville and the radial Schrödinger equations [1, 3, 8, 18] and the finite element approximations for free longitudinal vibrations problems of non–uniform rod [16, 19].

There are some traditional approaches to solving (1). The first approach is to reduce (1) to a standard eigenvalue problem [7, 14, 20]:

L^{−1}AL^{−T} (L^{T} x) = λ(L^{T} x),  (2)

where M = LL^{T} is the Cholesky factorization of M. Then the eigenvalue problem (2) can be solved by many efficient algorithms, such as the QR algorithm [14], the bisection/multisection algorithms [13], the divide–conquer algorithm [4, 5], and homotopy algorithm [11]. However, this approach is less attractive because it can not take advantage of the tridiagonal form of A and M and a full matrix L^{−1}AL^{−T} is generated in the process. Furthermore, the accuracy of this method also depends on the conditioning of M, since the inverse of L is explicitly required. The complexity of this method is O(n^{3}).

The second approach, the QZ method [20], disregards the symmetry and the tridiagonal structure of the problem and after a direct phase requiring O(n^{3}) operations, enters an iterative phase requiring O(n^{3}) operations to simultaneously reduce A and M to triangular–quasi–triangular form. This algorithm is not sensitive to the conditioning of M, though it is even more expensive than the first approach.

Another approach is the bisection/multisection method [15]. This parallel method is very easy to implement in practice and is reliable and accurate. However, the speed of the convergence is slow.

In this paper, we propose a method that computes all eigenvalues of (1) through finding zeroes of the polynomial equation

p(λ) = \frac{1}{α} det(a − λM),  (3)

where α is (−1)^{n} times det(M). Since M is positive definite, α ≠ 0.

Our method can be divided into two stages: root isolation and root extraction. The root isolation can be done by bisections or multisections of order n_{p}, the number of processors. A multisection splits a given interval [a, b] into n_{p} subintervals [μ_{i}, μ_{i+1}], where μ_{i} = a + i(b − a)/n_{p}, i = 0, 1, 2, …, n_{p}, by computing Sturm sequences at μ_{i}. As a result, the number of roots of p(λ) which are smaller than each μ_{i} are obtained and so are the number of roots in each subinterval. The process is recursively applied to every subinterval containing more than one root. This step results in a list of intervals consisting of single roots. When multiple roots are present, the isolation process continues until the required precision is reached. The root extraction shall be done by the Durand–Kerner root–finding method which shall be discussed in detail in section 3. We shall also show how we can virtually eliminate the overflow and underflow problems in section 2 and 3. In section 4, we shall show that our method is reliable and backward stable.

There are other approaches to this problem. In [2, 17], divide–and–conquer methods are proposed; the method based on the determinant evaluation and the Laguerre iterations is suggested in [10, 12]; and the parallel implementations of the multisection method and Lanczos method are discussed in [15].

2. The Sturm Sequences. Without loss of generality, we assume the matrix pencil (A, M) is unreducible, that is, β_{i}^{2} + δ_{i}^{2} ≠ 0, i = 2, …, n, where

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

(4)

It is well known that polynomials defined by

\{ \begin{matrix}ρ_{0}(λ) = 1, ρ_{1}(λ) = α_{1} − λγ_{1}, \\ ρ_{i}(λ) = (α_{i} − λγ_{i})ρ_{i−1}(λ) − (β_{i} − λδ_{i})^{2}ρ_{i−2}(λ), & i = 2, 3, …, n, \end{matrix}  (5)

form a Sturm sequence and

p(λ) = \frac{1}{α}ρ_{n}(λ).

The number of sign changes between consecutive terms of {ρ_{i}}_{i=1}^{n}, denoted by κ(λ) is equal to the number of roots of p(λ) that are less than λ.

The Sturm sequence (5) may suffer from severe overflow and underflow problems and require frequent scaling [10]. To remedy this, we propose following alternative recurrences. Let

ξ_{i} = \frac{ρ_{i}}{ρ_{i−1}} , i = 1, 2, …, n.  (6)

Dividing both sides of (5) by ρ_{i−1} yields,

\{ \begin{matrix}ξ_{1} = α_{1} − λγ_{1}, \\ ξ_{i} = α_{i} − λγ_{i} − \frac{(β_{i}−λδ_{i})^{2}}{ξ_{i−1}} , & i = 2, 3, …, n. \end{matrix}  (7)

To avoid possible breakdown, namely, when ξ_{i} = 0 for certain i, the following adjustment is asserted:

• If ξ_{1} = 0, set ξ_{1} = γ_{1}ϵ^{2}.
• If ξ_{i} = 0, i > 1, set ξ_{i} = \frac{(|β_{i}|+|λδ_{i}|)^{2}ϵ^{2}}{ξ_{i−1}} .

That is, if ξ_{i} = 0, we perturb the corresponding entries beyond the last significant digit stored in the machine.

Obviously,

ρ_{i} = ∏_{k=1}^{i} ξ_{k}  (8)

and κ(λ), the number of roots of p(λ) that are less than λ, is now equal to the number of negative terms in {ξ_{i}}_{i=1}^{n}.

(7) virtually eliminates hazards of overflow and underflow problems to compute κ(λ) because of its self–scaling.

3. The Durand–Kerner Method. Since

p(λ) = \frac{1}{α} det(A − λM),

it can be rewritten as

p(λ) = λ^{n} + a_{1}λ^{n−1} + ⋯ + a_{n−1}λ + a_{n}.  (9)

Assume λ_{1}, λ_{2}, …, λ_{n} are zeroes of p(λ).

Let

h(λ) = (λ − μ_{1})(λ − μ_{2}) ⋯ (λ − μ_{n}).

The following iterative formula

μ_{i}^{(k+1)} = μ_{i}^{(k)} − \frac{p(μ_{i}^{(k)})}{h′(μ_{i}^{(k)})}, i = 1, 2, …, n  (10)

is called the Durand–Kerner root–finding iteration [6, 9]. Numbers μ_{i}^{(k)}, i = 1, 2, …, n, are a set of approximations to n zeroes of p(λ). Numbers μ_{i}^{(k+)}, i = 1, 2, …, n, are a set of, normally, better approximations.

\underline{Theorem 3.1}. Let μ_{i}^{(0)}, i = 1, 2, …, n be arbitrary real numbers and μ_{i}^{(1)}, i = 1, 2, …, n be computed from (10). Then

∑_{i=1}^{n} μ_{i}^{(1)} = ∑_{i=1}^{n} λ_{i},  (11)

where λ_{i}’s are eigenvalues.

\underline{Proof}. Since λ_{1}, λ_{2}, …, λ_{n} are roots of p(λ) = 0,
p(λ) = (λ − λ_{1})(λ − λ_{2}) ⋯ (λ − λ_{n}).  (12)

Comparing (9) with (12) yields

a_{1} = − ∑_{i=1}^{n} λ_{i}.

On the other hand, by the results from [6, 9],

a_{1} = − ∑_{i=1}^{n} μ_{i}^{(1)}.

Therefore,

∑_{i=1}^{n} μ_{i}^{(1)} = ∑_{i=1}^{n} λ_{i}.

Theorem 3.1 implies that after one iteration, the sum of n approximations is equal to the sum of all roots of p(λ), even if initial approximations are arbitrary. The sum of n improved approximations is always constant and equal to the sum of all roots of p(λ).

It can be shown [6] that the iteration (10) is the Newton method applied to the function f(x) = (f_{1}(x), f_{2}(x), …, f_{n}(x))^{T}, where x = (x_{1}, x_{2}, …, x_{n})^{T} and

f_{1}(x) = x_{1} + x_{2} + ⋯ + x_{n} + a_{1},

⋮

f_{k}(x) = ∑_{i_{1}<i_{2}<···<i_{k}} x_{i_{1}}x_{i_{2}} ⋯ x_{i_{k}} + (−1)^{k−1}a_{k},

⋮

f_{n}(x) = ∏_{i=1}^{n} x_{i} + (−1)^{n−1}a_{n}.

The Newton iteration implies that the convergence is quadratic whenever the Jacobian of f(x) is nonsingular at (λ_{1}, λ_{2}, …, λ_{n})^{T}, i.e., whenever those roots λ_{i} are simple and there are r_{i} components of x which converge to λ_{i}, if the multiplicity of λ_{i} is r_{i}. In the case of multiple zeroes, the convergence is only linear. However, the case of multiple zeroes can be eliminated since our method consists of two parts: root isolation and root extraction. Without loss of generality, assume λ_{1} is a multiple root with multiplicity m and λ_{m+1}, …, λ_{n} are distinct. After the root isolation, λ_{1} is obtained. Therefore, we use the modified Durand–Kerner method:

\begin{matrix}μ_{i}^{(k+1)} = λ_{1}, i = 1, 2, …, m. \\ μ_{i}^{(k+1)} = μ_{i}^{(k)} − \frac{p(μ_{i}^{(k)})}{h′(μ_{i}^{(k)})}, i = m + 1, m + 2, …, n. \end{matrix}

The convergence is still quadratic.

Although (10) always converges in practice [6, 9] and no counterexample has been given so far, no proof of this property has been given. Many authors conjecture that this scheme converges for almost any starting points. However, to our special problem, we have the following result.

\underline{Theorem 3.2}. Assume

μ_{1} < λ_{1} < μ_{2} < λ_{2} < ⋯ < μ_{[ \frac{n}{2} ]} < λ_{[ \frac{n}{2} ]} < λ_{[ \frac{n}{2} ]+1} < μ_{[ \frac{n}{2} ]+1} < ⋯ < λ_{n} < μ_{n},

where λ_{1}, λ_{2}, …, λ_{n} are roots of p(λ) = 0 and μ_{1}, μ_{2}, …, μ_{n} are starting points for (10). Then,

μ_{i} < μ_{i}^{(1)}, if i ≤ \frac{n}{2}

and

μ_{i} > μ_{i}^{(1)}, if i > \frac{n}{2},

where μ_{1}^{(1)}, μ_{2}^{(1)}, …, μ_{n}^{(1)} are computed from (10).

\underline{Proof}.

μ_{i}^{(1)} = μ_{i} − \frac{p(μ_{i})}{h′(μ_{i})}

= μ_{i} − \frac{∏_{j=1}^{n}(μ_{i} − λ_{j})}{∏_{{j=1 \atop j≠i}} (μ_{i} − μ_{j})}

= μ_{i} − (μ_{i} − λ_{i})\frac{∏_{{j=1 \atop j≠i}} (μ_{i} − λ_{j})}{∏_{{j=1 \atop j≠i}} (μ_{i} − μ_{j})}.

Clearly,

\frac{∏_{{j=1 \atop j≠i}} (μ_{i} − λ_{j})}{∏_{{j=1 \atop j≠i}} (μ_{i} − μ_{j})} > 0, for i = 1, 2, …, n,

therefore,

μ_{i} < μ_{i}^{(1)}, if i ≤ \frac{n}{2}

and

μ_{i} > μ_{i}^{(1)}, if i > \frac{n}{2}.

Theorem 3.1 and 3.2 suggest that we should choose starting points like those we mentioned in Theorem 3.2. After each iteration, we should have better approximations. In case μ_{i}^{(k+1)} > μ_{i+1}^{(k)} happens for i ≤ \frac{n}{2} or μ_{i}^{(k+1)} < μ_{i−1}^{(k)} for i ≥ \frac{n}{2}, we give up μ_{i}^{(k+1)} and do bisection at this point to insure the convergence.

From (8), (10) can be rewritten as

μ_{i}^{(k+1)} = μ_{i}^{(k)} − \frac{∏_{j=1}^{n} ξ_{j}}{α ∏n ; {j=1 \atop j≠i} (μ_{j}^{(k)} − μ_{i}^{(k)})}  (13)

= μ_{i}^{(k)} − \frac{ξ_{i} ∏n ; {j=1 \atop j≠i} ξ_{j}}{α ∏n ; {j=1 \atop j≠i} η_{ji}}

= μ_{i}^{(k)} − \frac{ξ_{i}}{α} ∏_{{j=1 \atop j≠i}}^{n} \frac{ξ_{j}}{η_{ji}}

i = 1, 2, …, n,

where η_{ji} = μ_{j}^{(k)} − μ_{i}^{(k)}.

Directly computing those products in (13) may lead to severe overflow and underflow problems. To avoid the problems, we do the following:

(1). Compute v_{j} = \frac{ξ_{j}}{η_{ji}}, j = 1, 2, …, n and j ≠ i. (2). If |v_{j}| ≤ 1, let v_{j} in S. If |v_{j}| > 1, let v_{j} in B. Then S = {s_{1}, s_{2}, … } and

B = {b_{1}, b_{2}, … }.

(3). Let w = s_{1}b_{1}s_{2}b_{2} · · ·. Hence,

μ_{i}^{(k+1)} = μ_{i}^{(k)} − wξ_{i}/α, i = 1, 2, …, n.  (14)

In this way, overflow and underflow problems can be eliminated.

Since μ_{i}^{(k+1)} and μ_{j}^{(k+1)}, i ≠ j, can be computed independently, the method is parallel in nature.

4. Error Analysis. Let fl(•) denote the floating point computation of •. In practice, when we intend to evaluate zeros of p(λ), actually zeros of fl[p(λ)] are evaluated. We now do some error analysis and thereby establish the backward stability of our algorithm.

The most frequently used model for floating point arithmetic is

fl(x ∘ y) = (x ∘ y) · (1 + e), |e| ≤ ϵ,

where ∘ is either +, −, ×, or / and ϵ above is the machine precision, which is approximately 6 × 10^{−8} for single precision and 2.2 × 10^{−16} for double precision under IEEE standards.

\underline{Proposition 4.1}. [12]
fl[f(λ)] = (1 + ϵ_{1}) det[(A + δA) − λ(M + δM)],  (15)

where

|ϵ_{1}| ≤ nϵ,

and both δA and δM are symmetric tridiagonal matrices satisfying entrywise inequalities

|δA| ≤ 2.51ϵ|A|, |δM| ≤ 3.51ϵ|M|.  (16)
Assume λ_{min}(M) > 3.51ϵ∥M∥_{∞}. Then by Proposition 4.1,

\begin{aligned}λ_{min}(M + δM) & > λ_{min} − ∥δM∥_{∞} \\ & ≥ λ_{min} − 3.51ϵ∥M∥_{∞} > 0. \end{aligned}

Therefore, (A+δA, M +δM) is a symmetric, positive definite, tridiagonal pencil, if λ_{min}(M) > 3.51ϵ∥M∥_{∞} and this implies that, unless the pencil (A, M) is extremely ill–conditioned, for which λ_{min} < 3.51ϵ∥M∥_{∞}, the floating point computation of our algorithm is still performed on a symmetric definite pencil.

This establishes the backward stability of our algorithm.

5. Conclusion. In this paper, we have shown that the Durand–Kerner root– finding method converges for the generalized tridiagonal eigenvalue problem. The asymptotic convergence rate is quadratic. After each iteration, the sum of the improved approximations is the sum of all the eigenvalues. This indicates why our method converges so fast in practice. We have also shown how overflow and underflow problems can be eliminated. These good properties and schemes make our method an excellent candidate for the generalized tridiagonal eigenvalue problem.

\underline{References}

1. A. E. Berger, J. M. Solomon, and M. Ciment, “Higher Order Accurate Tridiagonal Difference Methods for Diffusion Convection Equations,” Advances in Computer Methods for Partial Differential Equations III, R. Vichnevetsky and R. S. Stepleman, editors, IMACS, 1978, 322–330.

2. C. F. Borges and W. B. Gragg, A Parallel Divide and Conquer Algorithm for the Generalized Real Symmetric Definite Tridiagonal Eigenproblem, Naval Postgraduate School, preprint, 1992.

3. R. F. Boisvert, “Families of Higher Order Accurate Discretizations of Some Elliptic Problems,” SIAM J. Sci. Stat. Comput., 2 (1981), 268–284.

4. J. J. M. Cuppen, “A Divide and Conquer Method for the Symmetric Tridiagonal Eigenproblem,” Numer. Math., 36 (1981), 177–195.

5. J. J. Dongarra and D. C. Sorensen, “A Fully Parallel Algorithm for the Symmetric Eigenvalue Problem,” SIAM J. Sci. Stat. Comput., 8 (1987), 139–154.

6. P. Fraigniaud, “The Durand–Kerner Polynomials Roots–Finding Method in Case of Multiple Roots,” BIT, 31 (1991), 112–123.

7. G. H. Golub and Van Loan, Matrix Computations, 2nd ed., The Johns Hopkins University Press, Baltimore, Maryland, 1989.

8. B. Grieves and D. Dunn, “Symmetric Matrix Methods for Schrödinger Eigenvectors,” J. Phys. A: Math Gen., Vol. 23, 5479–5491.

9. G. Kjellberg, “Two Observations on Durand–Kerner’s Root–Finding Method,” BIT, 24 (1984), 556–559.

10. K. Li, “A Fully Parallel Method for Tridiagonal Eigenvalue Problem,” Internat. J. of Math. & Math. Sci., 17 (1994), 741–752.

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

12. K. Li, T. Y. Li, and Z. Zeng, “An Algorithm for Generalized Symmetric Tridiagonal Eigenvalue Problem,” Numerical Algorithms, 8 (1994), 269–291.

13. S–S. Lo, B. Phillipe, and A. Sameh, “A Multiprocessor Algorithm for the Symmetric Tridiagonal Eigenvalue Problems,” SIAM J. Sci. Stat. Comput., 8 (1987), 155–165.

14. B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice–Hall, Inc., Englewood Cliffs, New Jersey, 1980.

15. B. Philippe and B. Vital, “Parallel Implementations for Solving Generalized Eigenvalue Problems with Symmetric Sparse Matrices,” Applied Numerical Mathematics, 12 (1993), 391–402.

16. Y. M. Ram and G. M. L. Gladwell, “Constructing a Finite–Element Model of a Vibratory Rod From Eigendata,” J. of Sound and Vibration, 165 (1993), 36–65.

17. C. J. Ribbens and C. Beattie, “Parallel Solutions of Generalized Symmetric Tridiagonal Eigenvalue Problem on Shared Memory Multiprocessors,” Technical Report TR 92-47, Department of Computer Science, Virginia Polytechnic Institute and State University, 1992.

18. G. Strang and G. Fix, An Analysis of the Finite Element Method, Prentice– Hall, Inc., Englewood Cliffs, New Jersey, 1973.

19. W. Weaver, Jr. and P. R. Johnson, Finite Elements for Structural Analysis, Prentice–Hall, New Jersey, 1984.

20. J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford, 1965.

Kuiyuan Li Department of Mathematics and Statistics University of West Florida Pensacola, FL 32514 email: kli@math-stat.math.uwf.edu