A ROOT–FINDING METHOD FOR THE GENERALIZED TRIDIAGONAL EIGENPROBLEM

Kuiyuan Li and Matoteng Ncube

Abstract. In this paper a globally convergent, parallel method for the generalized real symmetric eigenproblem is presented. The method can be easily implemented on computers.

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

Ax = λMx,  (1)

where A and M are both real symmetric tridiagonal and one of them, say M, is positive definite. 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 numerical solutions of Sturm–Liouville and radial Schrödinger equations [1,9] and finite element approximation for free longitudinal vibrations problems of non–uniform rod [10].

There are some approaches to solving (1). The first approach is to reduce (1) to a standard eigenvalue problem [11]:

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 [11], the bisection/multisection algorithms [7], the divide–conquer algorithm [2], and the homotopy algorithm [5]. However, this approach is less attractive because it cannot 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 [11], 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 very sensitive to the conditioning of M though it is even more expensive than the first approach.

The third approach is the bisection/multisection method [8]. This parallel method is very easy to implement in practice and is reliable and accurate. However, the rate of the convergence is linear.

Another approach is the homotopy–like method [6]. The method employs the determinant evaluation, split–and–merge strategy and the Laguerre’s iteration. Since Laguerre’s iteration requires us to compute the first and second derivatives of f(λ) = det(A − λM), more computations are needed.

Recently, the Durand–Kerner root finding method [4] was proposed for this problem. The method does not require us to compute f′(λ) or f′′(λ). Although it always converges in practice [4], and no counterexample has been given so far, no proof of this property has been given.

In this paper, we present a parallel method which converges monotonically and quadratically without computing f′(λ) or f′′(λ).

2. Proposed Root Finding Method. Let p(x) = \frac{1}{α} det(A − λM), where α = (−1)^{n} det(M). Since M is positive definite, α ≠ 0. Clearly, to find the eigenvalues, one only needs to find the zeros of p(x).

For this purpose, we propose the following root–finding method. Let

x_{i}^{(k+1)} = x_{i}^{(k)} − \frac{p(x_{i}^{(k)})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})}, (3)

for i = 1, 2, …, n and

y_{i}^{(k+1)} = y_{i}^{(k)} − \frac{p(y_{i}^{(k)})}{(y_{i}^{(k)} − x_{1}^{(k)}) ⋯ (y_{i}^{(k)} − x_{i−1}^{(k)})(y_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (y_{i}^{(k)} − y_{n}^{(k)})}, (4)

for i = 1, 2, …, n, where the numbers x_{i}^{(k)} and y_{i}^{(k)}, i = 1, 2, …, n, are sets of approximations to n zeroes of p(x). The numbers x_{i}^{(k+1)} and y_{i}^{(k+1)}, i = 1, 2, …, n, are sets of better approximations than x_{i}^{(k)} and y_{i}^{(k)}, i = 1, 2, …, n, respectively.

Clearly, (3) and (4) give us a fully parallel method for finding all zeros of p(x), that is, x_{i}^{(k+1)} and y_{i}^{(k+1)}, i = 1, 2, …, n can be computed independently.

\underline{Theorem 2.1}. Assume x_{i}^{(k)} < λ_{i} < y_{i}^{(k)} < x_{i+1}^{(k)} < λ_{i+1} < y_{i+1}^{(k)}, i = 1, 2, …, n− 1 and for some k, then for any i,

x_{i}^{(k)} < x_{i}^{(k+1)} < λ_{i}, for k > 0,

and

y_{i}^{(k)} > y_{i}^{(k+1)} > λ_{i}, for k > 0.

\underline{Proof}.

x_{i}^{(k+1)} = x_{i}^{(k)} − \frac{p(x_{i}^{(k)})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})}

= x_{i}^{(k)} − \frac{(x_{i}^{(k)} − λ_{1})(x_{i}^{(k)} − λ_{2}) ⋯ (x_{i}^{(k)} − λ_{n})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})}

= x_{i}^{(k)} − \frac{(x_{i}^{(k)} − λ_{1}) ⋯ (x_{i}^{(k)} − λ_{i−1})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})} \frac{(x_{i}^{(k)} − λ_{i+1}) ⋯ (x_{i}^{(k)} − λ_{n})}{(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})} · (x_{i}^{(k)} − λ_{i}).

Since

\frac{(x_{i}^{(k)} − λ_{1}) ⋯ (x_{i}^{(k)} − λ_{i−1})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})} > 0,

\frac{(x_{i}^{(k)} − λ_{i+1}) ⋯ (x_{i}^{(k)} − λ_{n})}{(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})} > 0,

and

x_{i}^{(k)} − λ_{i} < 0,
x_{i}^{(k+1)} > x_{i}^{(k)}.

On the other hand,

\frac{(x_{i}^{(k)} − λ_{1}) ⋯ (x_{i}^{(k)} − λ_{i−1})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})} < 1,

\frac{(x_{i}^{(k)} − λ_{i+1}) ⋯ (x_{i}^{(k)} − λ_{n})}{(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})} < 1,

and

x_{i}^{(k)} − λ_{i} < 0,

therefore,

x_{i}^{(k+1)} < x_{i}^{(k)} − (x_{i}^{(k)} − λ_{i}) = λ_{i}.

Similarly, we can prove that y_{i}^{(k)} > y_{i}^{(k+1)} > λ_{i}, for k > 0. Theorem 2.1 implies that x_{i}^{(k)} and y_{i}^{(k)} converge to λ_{i} strictly monotonically. \underline{Theorem 2.2}. Assume x_{i}^{(k)} < λ_{i} < y_{i}^{(k)} < x_{i+1}^{(k)} < λ_{i+1} < y_{i+1}^{(k)}, i = 1, 2, …, n−

1, then λ_{i} − x_{i}^{(k+1)} < y_{i}^{(k)} − λ_{i} and y_{i}^{(k+1)} − λ_{i} < λ_{i} − x_{i}^{(k)} for sufficiently large k.

\underline{Proof}. Let

ϵ = max_{1≤j≤n}{λ_{j} − x_{j}^{(k)}, y_{j}^{(k)} − λ_{j}}.

For any i,

|x_{i}^{(k+1)} − λ_{i}|

= |(x_{i}^{(k)} − λ_{i}) − \frac{(x_{i}^{(k)} − λ_{1})(x_{i}^{(k)} − λ_{2}) ⋯ (x_{i}^{(k)} − λ_{n})}{(x_{i}^{(k)} − x_{1}^{(k)}) ⋯ (x_{i}^{(k)} − x_{i−1}^{(k)})(x_{i}^{(k)} − y_{i+1}^{(k)}) ⋯ (x_{i}^{(k)} − y_{n}^{(k)})}|

= |(x_{i}^{(k)} − λ_{i}) − ∏_{j=1}^{i−1} \frac{x_{i}^{(k)} − λ_{j}}{x_{i}^{(k)} − x_{j}^{(k)}} ∏_{j=i+1}^{n} \frac{x_{i}^{(k)} − λ_{j}}{x_{i}^{(k)} − y_{j}^{(k)}} · (x_{i}^{(k)} − λ_{i})|

= |(x_{i}^{(k)} − λ_{i}) − ∏_{j=1}^{i−1}(1 + \frac{x_{j}^{(k)} − λ_{j}}{x_{i}^{(k)} − x_{j}^{(k)}} ) ∏_{j=i+1}^{n}(1 + \frac{x_{j}^{(k)} − λ_{j}}{x_{i}^{(k)} − y_{j}^{(k)}} ) · (x_{i}^{(k)} − λ_{i})|

= |(x_{i}^{(k)} − λ_{i})| | ∑_{j=1}^{i−1} \frac{x_{j}^{(k)} − λ_{j}}{x_{i}^{(k)} − x_{j}^{(k)}} + ∑_{j=i+1}^{n} \frac{x_{j}^{(k)} − λ_{j}}{x_{i}^{(k)} − y_{j}^{(k)}} + O(ϵ^{2})|

≤ |x_{i}^{(k)} − λ_{i}| |ϵ · (∑_{j=1}^{i−1} \frac{1}{x_{i}^{(k)} − x_{j}^{(k)}} + ∑_{j=i+1}^{n} \frac{1}{y_{j}^{(k)} − x_{i}^{(k)}} )O(ϵ^{2})|

≤ ϵ(ϵ(∑_{j=1}^{i−2} \frac{1}{λ_{i−1} − λ_{j}} + \frac{1}{x_{i} − λ_{i−1}} + ∑_{j=i+1}^{n} \frac{1}{λ_{j} − λ_{i}} ) + O(ϵ^{2}))

≤ ϵ^{2}(∑_{j=1}^{i−2} \frac{1}{λ_{i−1} − λ_{j}} + \frac{1}{λ_{i} − λ_{i−1} − ϵ} + ∑_{j=i+1}^{n} \frac{1}{λ_{j} − λ_{i}} ) + O(ϵ^{3})

≤ cϵ^{2} ≤ c(y_{i}^{(k)} − λ_{i})^{2}.

Since y_{i}^{(k)} converges to λ_{i}, for sufficiently large k, c(y_{i}^{(k)} − λ_{i}) < 1. Therefore, λ_{i} − x_{i}^{(k)} < y_{i}^{(k)} − λ_{i}.

Similarly, we can show that y_{i}^{(k+1)} − λ_{i} < λ_{i} − x_{i}^{(k)}.

It seems that to compute the approximation of λ_{i}, two sequences are needed, and therefore, more computations are needed. However, Theorem 2.2 shows that more computations may not be really needed. We do two iterations at each step.

Theorem 2.2 also implies that the root finding method converges quadratically.

\underline{Remark}. The method we discussed is a parallel method, but at the same time it is possible to use a Gauss–Seidel–like procedure, that is, we use the new approximations, x_{j}^{(m+1)}, where j = 1, 2, …, i−1, to replace x_{j}^{(m)} and y_{j}^{(m+1)}, where j = n, n − 1, …, i + 1, to replace y_{j}^{(m)} in (3) and (4). We compute x_{j}^{(m+1)} and y_{j}^{(m+1)} in the following order: x_{1}^{(m+1)}, y_{n}^{(m+1)}, x_{2}^{(m+1)}, y_{n−1}^{(m+1)}, …, x_{n}^{(m+1)}, and y_{1}^{(m+1)}. In this case, the sequences converge to λ_{i} super quadratically.

3. Numerical Results. In this section, we present our numerical results. Our method for finding all eigenvalues can be divided into two stages: root

isolation and root extraction. The root isolation can be done by using the Sturm sequence. This technique is discussed in detail in [4]. This step results in a list of intervals containing single roots. When multiple roots are present, the isolation process continues until the required precision is reached. The root extraction step uses the proposed root–finding method.

\underline{Experiment 1}. We implemented our algorithm LGT and the routine RSG in EISPACK [3] on equation (1). The matrices A and M are obtained from piecewise linear finite element discretization [9] of the Sturm–Liouville problem

−\frac{d}{dx}(p(x)\frac{du}{dx}) + q(x)u = λu,

where u = u(x), 0 < x < π, u(0) = u′(π) = 0, and p(x) > 0. When [0, π] is divided into n + 1 subintervals of equal length, equation (1) is obtained. Here, both A and M are symmetric tridiagonal, and positive definite. To make the problem simpler, we use p(x) = 1 and q(x) = 6.

The experiment was conducted on a Sun–1000 workstation with IEEE standard and machine precision ϵ ≈ 2.2 × 10^{−16}.

For the generalized symmetric tridiagonal eigenproblem, the subroutine RSG in EISPACK is somewhat less attractive since it transforms the generalized eigenproblem Ax = λMx to a standard eigenproblem T̃y = λy of a dense matrix T̃. It cannot take advantage of the tridiagonal structure of the pencil (A, M) and clearly requires O(n^{3}) flops. The RSG is also sensitive to the conditioning of M. We conduct the comparison here mainly because it is the only algorithm available in EISPACK for eigenproblems of symmetric definite pencils.

Table 1 shows the results of our algorithm LGT and the algorithm RSG. The execution time apparently shows that the complexities of LGT and RSG are O(n^{2}) and O(n^{3}), respectively for the generalized tridiagonal problem.

Order Execution Time of All Eigenvalues

n  RSG  LGT
100  2.25  0.31
200  19.02  1.21
300  70.29  2.79
400  194.71  4.87

Table 1: Execution time (seconds) of computed eigenvalues.

\underline{Experiment 2}. Our algorithm LGT and the routine BIS, the bisection method are applied to equation (1), where A and M are symmetric random matrices with m_{i,i} = 2(m_{i,i−1} + m_{i,i+1}) so that M is positive definite. For each n, we ran 50 problems to find the average time. Both the LGT and BIS can fully take advantage of the symmetric tridiagonal forms of A and M. Since the bisection converges only linearly, Table 2 shows that our algorithm leads in speed by a considerable margin.

Order Execution Time of All Eigenvalues

n  BIS  LGT
100  0.76  0.34
200  3.09  1.26
300  6.47  2.99
400  11.46  5.04
1000  77.19  27.24

Table 2: Execution time (seconds) of computed eigenvalues.

\underline{References}

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

2. 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.

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

4. K. Li, “Durand–Kerner Root–Finding Method for the Generalized Tridiagonal Eigenproblem,” Missouri Journal of Mathematical Sciences, (to appear).

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

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

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

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

9. G. Strang and G. Fix, An Analysis of the Finite Element Method, Prentice– Hall, Inc., Englewood Cliffs, N.J., 1973.

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

11. 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

Matoteng Ncube Department of Mathematics and Statistics University of West Florida Pensacola, FL 32514