CALCULATING AN INVARIANT SUBSPACE OF

DIAGONALLY DOMINANT MATRICES – PART II

Noah H. Rhee

University of Missouri – Kansas City

Abstract. In part I of this paper, a modification of the algorithm of Blevins and Stewart for calculating an invariant subspace of diagonally dominant matrices was given. Here, we will discuss the numerical aspects of this algorithm. In particular, numerical results show that the modified algorithm seems to converge faster than the original algorithm.

The notation and terminology we will use will be the same as in Part I [4] of the article.

4. Numerical Results. In practice ϵ = ||E||_{2} is difficult to calculate. But we note that ϵ̃ = ||E_{11}|| + ||E_{22}|| is easy to compute and ϵ ≤ ϵ̃. This can be seen by using [1] as follows. Using the definition of spectral norm and relation (9) we obtain

ϵ ≤ ||E||_{2} = sup_{||P ||=1} ||EP ||

= sup_{||P ||=1} ||P E_{11} − E_{22}P ||

≤ sup_{||P ||=1}(||P ||||E_{11}||_{2} + ||E_{22}||_{2}||P ||)

= ||E_{11}||_{2} + ||E_{22}||_{2} ≤ ||E_{11}|| + ||E_{22}|| = ϵ̃.

Let

(27)  Γ̃ = \{P ∈ ℝ^{(n−l)×l} : ||P || ≤ \frac{2γ}{δ − ϵ̃}\}.

Note that the following variant of Theorem 9 is also valid.

\underline{Theorem 13}. Let η = ||E_{12}||, γ = ||E_{21}||, ϵ̃ = ||E_{11}|| + ||E_{22}||, and δ = ||D^{−1}||_{2}^{−1}. Then if

(28)  δ − ϵ̃ > 2\sqrt{ηγ},

then, for any P_{0} ∈ Γ, the sequence P_{k}, defined by

P_{k+1} = Φ(P_{k}), k = 0, 1, …,

converges to the unique solution P^{∗} of equation (7) in Γ̃. Moreover,

(29)  ||P^{∗} − P_{k}|| ≤ \frac{ρ̃}{1 − ρ̃}||P_{k} − P_{k−1}||, k = 1, 2, …,

and

(30)  ||P_{k+1} − P_{k}|| ≤ ρ̃||P_{k} − P_{k−1}||, k = 1, 2, …,

where

(31)  ρ̃ = \frac{ϵ̃}{δ} + \frac{4ηγ}{δ(δ − ϵ̃)} < 1.

\underline{Remark 14}. Condition (28) is stronger than condition (18), since ϵ̃ ≥ ϵ. However, the advantage of Theorem 13 over Theorem 9 is that ϵ̃ is easy to compute. Note that η and γ are easy to compute. Hence, if we know how to compute δ, we can check convergence criterion (28), and compute ρ̃ in relation (31). But, from relation (11), it is easy to see [1] that

(32)  δ = ||D^{−1}||^{−1} = min {|d_{j} − d_{i}| : l + 1 ≤ i ≤ n, 1 ≤ j ≤ l}.

Under condition (28), Algorithm 11 becomes, with a slight change in step (1)

\underline{Algorithm 15}.

(1) \overline{P}_{0} = 0 and \overline{P}_{1} = Φ(\overline{P}_{0})

(2) For k = 1, 2, …, \overline{P}_{k+1} = Ψ(\overline{P}_{k}). If ||\overline{P}_{k+1}|| > 2γ/(δ − ϵ̃) or ||\overline{P}_{k+1} − \overline{P}_{k}|| > ρ̃||\overline{P}_{k} − \overline{P}_{k−1}|| (k ≥ 2), go to (3); otherwise go to (2).

(3) For j = k, k + 1, …, \overline{P}_{j+1} = Φ(\overline{P}_{j}). \underline{Remark 16}. The change in step (1) is due to the fact that Φ(0) is cheaper to compute

than Ψ(0).

Now for the given diagonally dominant matrix A with which the first l diagonal elements form a cluster that is well separated from the other diagonal elements, it is easy to check whether the condition (28) is satisfied or not. Once condition (28) is satisfied, we can locate an invariant subspace of A either by using Algorithm 5 or Algorithm 15. The iteration can be terminated with P_{k+1} (\overline{P}_{k+1}) if

||P_{k+1} − P_{k}||(||\overline{P}_{k+1} − \overline{P}_{k}||) ≤ eps

in both algorithms, where eps is a prescribed error tolerance. Note that from relation (29) this implies

||P^{∗} − P_{k+1}|| ≤ \frac{ρ̃}{1 − ρ̃}eps.

How to use Algorithm 5 is straightforward. Now we briefly discuss the practical details in using Algorithm 15. There are two things that need to be explained; (i) how to form \overline{P}_{k+1} from \overline{P}_{k} by using the function Ψ in step (2) and (ii) how to check two conditions in step (2).

(i) To compute \overline{P}_{k+1} from \overline{P}_{k} we have to solve (see relation (23))

(D + L)\overline{P}_{k+1} = E_{21} − \overline{P}_{k}E_{12}\overline{P}_{k} − U\overline{P}_{k}.

Let

\overline{Q}_{k} = E_{21} − \overline{P}_{k}E_{12}\overline{P}_{k} − U\overline{P}_{k}.

Then the above linear system becomes

(D + L)\overline{P}_{k+1} = \overline{Q}_{k}.

Let \overline{P}_{k+1} = (p_{ij}), \overline{Q}_{k} = (q_{ij}), l + 1 ≤ i ≤ n, 1 ≤ j ≤ l. Then from equation (25) we can express p_{ij} by the following equations

p_{ij} = \frac{q_{ij} + ∑_{k=l+1}^{i−1} l_{ik}p_{kj} − ∑_{k=j+1}^{l} p_{ik}l_{kj}}{d_{j} − d_{i}} , l + 1 ≤ i ≤ n, 1 ≤ j ≤ l.

Here the empty sums are regarded as zero. By using the ordering defined in relation (26), the above equations can be solved by forward substitution, since the matrix representation of L + D is a lower triangular matrix.

(ii) To check the terminating condition we need to form

σ_{k} = ||\overline{P}_{k+1} − \overline{P}_{k}||, k = 0, 1, ….

To check the second condition in step (2), all we need to check is whether σ_{k}/σ_{k−1} is bigger than ρ̃ or not for k = 2, 3, . . .. The first condition in step (2) can be checked by directly computing ||\overline{P}_{k+1}|| and compare it with 2γ/(δ − ϵ̃). However, a more economic way of checking this condition is to form σ_{0} + σ_{1} + ⋯ + σ_{k}, and compare this with 2γ/(δ − ϵ̃). Note that ||\overline{P}_{k+1}|| ≤ σ_{0} + σ_{1} + ⋯ + σ_{k}. Hence, if

(33)  σ_{0} + σ_{1} + ⋯ + σ_{k} > \frac{2γ}{δ − ϵ̃}

and

(34)  σ_{0} + σ_{1} + ⋯ + σ_{k−1} ≤ \frac{2γ}{δ − ϵ̃},

then ||\overline{P}_{k}|| ≤ 2γ/(δ − ϵ̃) and ||\overline{P}_{k+1}|| may be greater than 2γ/(δ − ϵ̃). So if conditions (33) and (34) are satisfied for k = 1, 2, …, we go to step (3). We also mention that in using Algorithm 15 we stored A^{T} instead of A to make the program column oriented. Now we conclude this paper with several numerical examples, which demonstrate faster convergence of Algorithm 15.

5. Numerical Examples. We will denote Algorithm 5 by algorithm (BS) and Algorithm 15 by algorithm (MBS). In both algorithms eps was set to 10^{−15}. All computations were done on a VAX 6000 Model 500 using double precision arithmetic, that is, the machine precision of 1.11 × 10^{−16}.

\underline{Example 16}. We generated a matrix A of order n = 40, whose off diagonal entries were determined by (2 random − 1) × 0.025, where random denotes a random number between 0 and 1. The diagonal entries were set d_{1} = d_{2} = d_{3} = 3 and d_{i} = i, 4 ≤ i ≤ n. l was chosen to be 3. Both algorithms were used to locate the invariant subspace spanned by eigenvectors whose eigenvalues are close to 3. The convergence condition (28) was satisfied. In fact, we have

\begin{matrix}δ − ϵ̃ = 0.4272, 2\sqrt{ηγ} = 0.3101 \\ \frac{2γ}{δ − ϵ̃} = 0.7579, ρ̃ = 0.7979 \end{matrix}

Then we obtained the following results:

(BS)  (MBS)
k  ||P_{k+1} − P_{k}||  ||\overline{P}_{k+1} − \overline{P}_{k}||
0  3.1075 · 10^{−2}  3.1075 · 10^{−2}
1  7.8481 · 10^{−4}  7.8347 · 10^{−4}
2  1.8771 · 10^{−5}  9.5336 · 10^{−6}
3  5.3214 · 10^{−7}  2.2330 · 10^{−7}
4  1.7232 · 10^{−8}  4.5688 · 10^{−9}
5  3.8518 · 10^{−10}  1.0147 · 10^{−10}
6  8.5870 · 10^{−12}  8.5379 · 10^{−13}
7  1.7331 · 10^{−13}  4.7367 · 10^{−15}
8  3.2774 · 10^{−15}  1.4187 · 10^{−17}
9  9.2689 · 10^{−17}

We note that 1 iteration was saved by using algorithm (MBS).

\underline{Example 17}. Let A be of order n = 40, whose upper (lower) off diagonal entries are uniformly 0.01 (0.02). The diagonal entries were set d_{1} = d_{2} = d_{3} = 3 and d_{i} = i, 4 ≤ i ≤ n. l was chosen to be 3. Both algorithms were used to locate the invariant subspace spanned by eigenvectors whose eigenvalues are close to 3. The convergence condition (28) was satisfied. In fact, we have

\begin{matrix}δ − ϵ̃ = 0.3842, 2\sqrt{ηγ} = 0.2980 \\ \frac{2γ}{δ − ϵ̃} = 1.0969, ρ̃ = 0.8469 \end{matrix}

Then we obtained the following results:

(BS)  (MBS)
k  ||P_{k+1} − P_{k}||  ||\overline{P}_{k+1} − \overline{P}_{k}||
0  4.4067 · 10^{−2}  4.4067 · 10^{−2}
1  1.3768 · 10^{−3}  1.3574 · 10^{−3}
2  8.4031 · 10^{−5}  8.0349 · 10^{−5}
3  3.4503 · 10^{−6}  9.4730 · 10^{−7}
4  1.8811 · 10^{−7}  1.1611 · 10^{−8}
5  9.0326 · 10^{−9}  9.6739 · 10^{−11}
6  4.9280 · 10^{−10}  4.9305 · 10^{−13}
7  2.6768 · 10^{−11}  5.1917 · 10^{−15}
8  1.5416 · 10^{−12}  3.1615 · 10^{−17}
9  9.0982 · 10^{−14}
10  5.4875 · 10^{−15}
11  3.3363 · 10^{−16}

We note that 3 iterations were saved by using algorithm (MBS). So in this example algorithm (MBS) was substantially faster than algorithm (BS).

\underline{Example 18}. Let A be of order n = 400, whose upper (lower) diagonal entries are uniformly 0.001 (0.002). The diagonal entries were set d_{1} = d_{2} = d_{3} = 3 and d_{i} = 4, 4 ≤ i ≤ n. l was chosen to be 3. Both algorithms were used to locate the invariant subspace spanned by eigenvectors whose eigenvalues are close to 3. The convergence condition (28) was satisfied. In fact, we have

\begin{matrix}δ − ϵ̃ = 0.3692, 2\sqrt{ηγ} = 0.0971 \\ \frac{2γ}{δ − ϵ̃} = 0.3739, ρ̃ = 0.6566 \end{matrix}

It was observed that ||P_{58} − P_{57}|| ≤ eps, while ||\overline{P}_{15} − \overline{P}_{14}|| ≤ eps. Hence, 43 iterations were saved by using algorithm (MBS). Accordingly the execution time of algorithm (BS) was 33.86 seconds, while the execution time of algorithm (MBS) was 9.24 seconds. This is an extreme example. However, this example demonstrates the point that to compute Ψ(P) and Φ(P) takes essentially the same time and also the point that in some cases, algorithm (MBS) is much faster than algorithm (BS).

\underline{Remark 19}. In the above three examples, step (3) of algorithm (MBS) was never used.

\underline{References}

1. M. M. Blevins and G. W. Stewart, “Calculating the Eigenvectors of Diagonally Dominant Matrices,” Journal of the Association for Computing Machinery, 21 (1974), 261– 271.

2. W. W. Hager, Applied Numerical Linear Algebra, Prentice Hall, Englewood Cliffs, N.J., 1988.

3. J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.

4. N. H. Rhee, “Calculating an Invariant Subspace of Diagonally Dominant Matrices – Part I,” Missouri Journal of Mathematical Sciences, 5 (1993), 71–82.

5. G. W. Stewart, “Error Bounds for Invariant Subspaces of Closed Operators,” SIAM J. Numer. Anal., 8 (1971), 196–208.

6. J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford Univ. Press, London/New York, 1965.