CALCULATING AN INVARIANT SUBSPACE OF

DIAGONALLY DOMINANT MATRICES – PART I

Noah H. Rhee

University of Missouri – Kansas City

Abstract. In this paper a modification of the algorithm of Blevins and Stewart for calculating an invariant subspace of diagonally dominant matrices is given. The relation between the algorithm of Blevins and Stewart and our modified algorithm is essentially the same as the relation between the Jacobi method and the Gauss–Seidel method in solving linear systems iteratively.

1. Introduction. Let A be a real matrix of order n, and write A in the form

A = D + E,

where

D = diag(d_{1}, d_{2}, …, d_{n})

is the diagonal part of A. Following the definition in [1] we shall say that A is diagonally dominant if

||E|| = σ||D||,

where σ < 1. Here || · || denotes the Frobenius matrix norm defined by

||C||^{2} = ∑_{i=1}^{m} ∑_{j=1}^{n} c_{ij}^{2},

where C = (c_{ij}) is any real m × n matrix. In this paper we are concerned with algorithms for computing an invariant subspace of A when σ << 1.

The calculation of an invariant subspace of diagonally dominant matrices arises in the problem of refining systems of approximate eigenvectors. [1]

Throughout this paper ℝ^{n} denotes real Euclidean space of dimension n and ℝ^{m×n} denotes the set of m × n real matrices. Linear operators on a vector space, as opposed to their matrix representations, are denoted by bold face letters. The superscript T will denote the transpose of a given vector or a matrix.

For the linear operator T from ℝ^{m×n} to ℝ^{k×l}, we define the spectral norm of T by

||T||_{2} = sup_{||P ||=1} ||TP ||.

We note that for any matrices B and C regarded as operators, ||BC|| ≤ ||B||_{2}||C|| and ||BC|| ≤ ||B||||C||_{2}, whenever the product BC is defined. We also note that ||B||_{2} ≤ ||B||.

In section 2, we briefly summarize the main results of Blevins and Stewart [1]. In section 3, we discuss some further results and give a modified algorithm of Blevins and Stewart. In part 2 of this article, we will give some numerical results which demonstrate the faster convergence of the modified algorithm.

2. The Main Results of Blevins and Stewart. Let A ∈ ℝ^{n×n} be a diagonally dominant matrix. There is a natural correspondence between the eigenvectors of A and those of D, since A is a diagonally dominant matrix. If d_{i} is distinct from the other elements of D and if σ is sufficiently small, then the Gerschgorin theorem shows that there is a unique eigenvalue of A which is near d_{i}. Since the eigenvector corresponding to the distinct eigenvalue is a continuous function of the entries of the matrix (see [5], p. 67), if the ith components of the corresponding eigenvectors of D and A are normalized to unity the other components of the eigenvectors of A must be small. For definiteness we shall compare the eigenvectors corresponding to the diagonal element d_{1}. If we write these eigenvectors of D and A in the forms (1, 0)^{T} and (1, p^{T})^{T}, respectively, then p ∈ ℝ^{n−1} is small.

The first step is to find the equation satisfied by p. Because (1, p^{T})^{T} is an eigenvector of A, the vector A(1, p^{T})^{T} lies in the same direction as (1, p^{T})^{T}. But the matrix

( \begin{matrix}−p^{T} \\ I_{n−1} \end{matrix} )

has n − 1 linearly independent columns, all of which are orthogonal to (1, p^{T})^{T}. Here I_{k} denotes the identity matrix in ℝ^{k×k}. It follows that

(1)  ( \begin{matrix}−p^{T} \\ I_{n−1} \end{matrix} )^{T} A ( \begin{matrix}1 \\ p \end{matrix} ) = 0.

If we partition A in the form

A = ( \begin{matrix}d_{1} & ϵ_{12}^{T} \\ ϵ_{21} & D_{2} + E_{22} \end{matrix} ) ,

then (1) becomes

(2)  T p = ϵ_{21} − pϵ_{12}^{T}p,

where T is the matrix T = d_{1}I_{n−1} − (D_{2} + E_{22}).

Any method for solving the nonlinear equation (2) for p is effectively a method for computing the eigenvector corresponding to the eigenvalue near to d_{1}. In theory, this is not difficult to do. In [4] it is shown that under the following condition,

(3)  ||T^{−1}||_{2}^{2}||ϵ_{21}||||ϵ_{12}|| < \frac{1}{4},

the sequence defined by

(4)  T p_{k+1} = ϵ_{21} − p_{k}ϵ_{12}^{T}p_{k}, p_{0} = 0

converges to a solution of (2). As a practical method, however, this iteration has two drawbacks. First, condition (3) may not be satisfied, that is, d_{1} may not be sufficiently well separated from the other diagonal elements of A. Second, the solution of equation (4) for p_{k+1} is expensive.

Blevins and Stewart disposed of the first problem by showing how to find an invariant subspace corresponding to a cluster of close eigenvalues. Then they disposed of the second problem by using an approximated inverse of T.

To explain the idea of Blevins and Stewart we define the following. \underline{Definition 1}. Let ℜ(A) denote the column space of the matrix A. \underline{Definition 2}. For a subspace Ω ⊆ ℝ^{n}, let dim(Ω) denote the dimension of the subspace

Ω.

\underline{Definition 3}. A subspace Ω ⊆ ℝ^{n} is an invariant subspace of A if AΩ ⊆ Ω.

We note that an eigenvector of A spans an invariant subspace of dimension unity. In order to compute an invariant subspace Ω of A, Blevins and Stewart set up an equation analogous to (2) for a basis for Ω. The following lemma indicates how this was done.

\underline{Lemma 4}. Let A ∈ ℝ^{n×n} and let X = [X_{1} X_{2}] be such that X_{1} ∈ ℝ^{n×l}, X_{2} ∈ R^{n×(n−l)}, X_{1}^{T} X_{2} = 0 and the columns of X_{1}(X_{2}) are linearly independent. Then a necessary and sufficient condition that ℜ(X_{1}) be an invariant subspace of A is

(5)  X_{2}^{T} AX_{1} = 0.

\underline{Proof}. This lemma is a modification of the Lemma 3.1 in [1]. Since

dim(ℜ(X_{1})) + dim(ℜ(X_{2})) = n

and X_{1}^{T} X_{2} = 0, ℜ(X_{2}) is the orthogonal complement of ℜ(X_{1}) in ℝ^{n}. The remaining proof is identical with the proof given in [1].

For definiteness suppose that the first l diagonal elements of A form a cluster that is well separated from the other diagonal elements, and partition A in the form

A = ( \begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} ) = ( \begin{matrix}D_{1} + E_{11} & E_{12} \\ E_{21} & D_{2} + E_{22} \end{matrix} ) ,

where D_{1} ∈ ℝ^{l×l}. And they attempted to find a basis for an invariant subspace of A that in some sense corresponds to the matrix D_{1}.

Because E_{21} is presumed small, equation (5) is very nearly satisfied by the matrix

X = ( \begin{matrix}I_{l} & 0 \\ 0 & I_{n−l} \end{matrix} ) .

This suggests that we seek X in the form

X = [X_{1} X_{2}] = ( \begin{matrix}I_{l} & −P^{T} \\ P & I_{n−l} \end{matrix} ) ,

where P ∈ ℝ^{(n−l)×l} and ||P || is small. Note that the columns of X_{1}(X_{2}) are linearly independent and X_{1}^{T} X_{2} = 0. Hence, by Lemma 4 the necessary and sufficient condition that ℜ(X_{1}) be an invariant subspace of A is that

X_{2}^{T} AX_{1} = 0.

With the partitions of A and X, the above equation becomes

(6)  P A_{11} − A_{22}P = E_{21} − P E_{12}P.
Then they defined the linear operator T: ℝ^{(n−l)×l} → ℝ^{(n−l)×l} by
TP = P A_{11} − A_{22}P.

So equation (6) becomes

(7)  TP = E_{21} − P E_{12}P.

Equation (7) is perfectly analogous to equation (2), and it has a small solution under analogous conditions. Since l is unrestricted, one is free to augment D_{1} until the conditions for the existence of a solution are satisfied, which disposes of the problem of very close diagonal elements.

Equation (7) can be solved by an iterative process analogous to (4); however, each iteration requires the solution of an equation of the form TP = Q, which is prohibitively expensive, even when l = 1. However, T can be written T = D + E, where

(8)  DP = P D_{1} − D_{2}P.

and

(9)  EP = P E_{11} − E_{22}P.

Then they proposed the following iterative process to solve equation (7).

\underline{Algorithm 5}.

(1) P_{0} = 0

(2) For k = 0, 1, …

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

where

(10)  Φ(P) = D^{−1}(E_{21} − P E_{12}P − EP).

Since E_{11} and E_{22} are small, the operator D is near the operator T. Since D_{1} and D_{2} are diagonal matrices, the equation of the form DP = Q can be easily solved. In fact, the matrix representation of the linear operator D is a diagonal matrix. This can be seen as follows.

Suppose DP = Q. Let P = (p_{ij}) and Q = (q_{ij}), l + 1 ≤ i ≤ n, 1 ≤ j ≤ l. Here we shifted the index i by l for notational convenience. From equation (8) we have

P D_{1} − D_{2}P = Q.

This equation is equivalent to

(11)  (d_{j} − d_{i})p_{ij} = q_{ij}, l + 1 ≤ i ≤ n, 1 ≤ j ≤ l.

Note that q_{ij} depends only on p_{ij}, l + 1 ≤ i ≤ n, 1 ≤ j ≤ l. This means that the matrix representation of the linear operator D is a diagonal matrix.

The iteration defined by Algorithm 5 converges to a solution of equation (7) under rather general conditions. They proved the following Theorem.

\underline{Theorem 6}. Let η = ||E_{12}||, γ = ||E_{21}||, ϵ = ||E||_{2}, and δ = ||D^{−1}||_{2}^{−1}. Then if

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

the sequence P_{k} defined by Algorithm 5 converges to a solution P^{∗} of equation (7) satisfying

(13)  ||P^{∗}|| ≤ \frac{2γ}{δ − ϵ}.

Moreover

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

and

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

where

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

To find out eigenvalues and eigenvectors from the given invariant subspace, they also proved the following theorem.

\underline{Theorem 7}. Let P^{∗} be determined as in Theorem 6. Then the eigenvalues of A corresponding to the invariant subspace determined by P^{∗} are the eigenvalues of the matrix A_{11} + E_{12}P^{∗}. Moreover, if the columns of Z form a complete set of eigenvectors for A_{11} + E_{12}P^{∗}, the corresponding eigenvectors of A are the columns of the matrix

( \begin{matrix}I \\ P^{∗} \end{matrix} ) Z = ( \begin{matrix}Z \\ P^{∗}Z \end{matrix} ) .

3. A Modified Algorithm. The idea of proving Theorem 6 in [1] is as follows. First, it was shown that all the iterates P_{k} generated by Algorithm 5 remain in the region defined by (13). Second, it was shown that the function Φ is a contraction, with constant ρ, in that region. Then the result with error bounds (14) and (15) follow from a variant of the contraction mapping theorem.

But we can prove a stronger result by using the standard contraction mapping theorem. \underline{Theorem 8}. Contraction Mapping Theorem.

Suppose that G: D ⊆ ℝ^{n} → ℝ^{n} maps a closed set D_{0} ⊆ D into itself and that

||Gx − Gy|| ≤ ρ||x − y||, x, y ∈ D_{0},

for some ρ < 1. Then, for any x_{0} ∈ D_{0}, the sequence generated by

x_{k+1} = Gx_{k}, k = 0, 1, …,

converges to the unique fixed point x^{∗} of G (that is, G(x^{∗}) = x^{∗}), in D_{0} and

||x^{∗} − x_{k}|| ≤ \frac{ρ}{1 − ρ}||x_{k} − x_{k−1}||, k = 1, 2, ….

\underline{Proof}. (See [3], p. 385) . Let

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

By using the above standard Contraction Mapping Theorem, we can improve Theorem 6 as follows.

\underline{Theorem 9}. Let η = ||E_{12}||, γ = ||E_{21}||, ϵ = ||E||_{2}, and δ = ||D^{−1}||_{2}^{−1}. If

(18)  δ − ϵ > 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

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

and

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

where

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

\underline{Proof}. The fact that Φ is a contraction in Γ has been established in [1], with the contraction constant ρ described in the relation (21). Clearly Γ is a closed set in ℝ^{(n−l)×l}. Now we establish the fact that Φ(Γ) ⊆ Γ. Suppose that ||P || ≤ 2γ/(δ − ϵ) and

||Φ(P)|| > 2γ/(δ − ϵ).

Then using relation (10) we obtain

\begin{aligned}\frac{2γ}{δ − ϵ} & < ||Φ(P)|| ≤ \frac{1}{δ} (γ + η||P ||^{2} + ϵ||P ||) \\ & ≤ \frac{1}{δ} [γ + η (\frac{2γ}{δ − ϵ})^{2} + ϵ (\frac{2γ}{δ − ϵ})] . \end{aligned}

Equivalently,

2δ(δ − ϵ) < (δ − ϵ)^{2} + 4ηγ + 2ϵ(δ − ϵ).

However, it can be seen easily by using the assumption δ−ϵ > 2\sqrt{ηγ} that the right hand side of the above inequality is less than 2δ(δ − ϵ). So we obtain a contradiction. It follows that Φ(Γ) ⊆ Γ. Hence, by the Contraction Mapping Theorem the sequence {P_{k}} converges to the unique fixed point P^{∗} of Φ in Γ. This means that P^{∗} is the unique solution of P^{∗} = Φ(P^{∗}) in Γ. By using the definition of Φ in (10) we see that P^{∗} is the unique solution of equation (7) in Γ.

\underline{Remark 10}. The differences between Theorem 6 and Theorem 9 are the following. First, Theorem 9 gives the information that P^{∗} is the unique solution of (7) in Γ, whereas Theorem 6 does not. Second, in Theorem 6 P_{0} has to be 0, whereas in Theorem 9, P_{0} can be any element in Γ. This fact will be used crucially in our modified algorithm.

Now we are going to suggest a modified algorithm. Recall that ((8), (9))

T = D + E,

where

DP = P D_{1} − D_{2}P,

and

EP = P E_{11} − E_{22}P.

Now we let

T = L + D + U,

where

LP = P L_{1} − L_{2}P,

and

UP = P U_{1} − U_{2}P.

Here L_{1}(U_{1}) and L_{2}(U_{2}) are strictly lower (upper) triangular parts of A_{11} and A_{22}, respectively. We note that L_{1}(U_{1}) and L_{2}(U_{2}) can also be viewed as the strictly lower (upper) triangular parts of E_{11} and E_{22}, respectively.

Then the corresponding algorithm to Algorithm 5, based on the above LDU splitting of T, is

\underline{Algorithm 10}.

(1) \overline{P}_{0} = 0
(2) For k = 0, 1, …,
\overline{P}_{k+1} = Ψ(\overline{P}_{k}),

where

(22)  Ψ(\overline{P}) = (D + L)^{−1} (E_{21} − \overline{P} E_{12}\overline{P} − U\overline{P} ) .
Of course, we determine \overline{P}_{k+1} by solving
(23)  (D + L)\overline{P}_{k+1} = E_{21} − \overline{P}_{k}E_{12}\overline{P}_{k} − U\overline{P}_{k},

which is easy to solve, since (D + L)P = P (D_{1} + L_{1}) − (D_{2} + L_{2})P and D_{i} + L_{i} is a lower triangular matrix for i = 1, 2. In fact, there is a lower triangular matrix representation of the operator D + L. This can be seen as follows.

Suppose (D + L)P = Q. Then we have

(24)  P (D_{1} + L_{1}) − (D_{2} + L_{2})P = Q.

As before we let P = (p_{ij}) and Q = (q_{ij}), l + 1 ≤ i ≤ n, 1 ≤ j ≤ l. We also let L_{1} = (l_{ij}), l ≤ i, j ≤ l, and L_{2} = (l_{ij}), l + 1 ≤ i, j ≤ n. If we use the fact that L_{1} and L_{2} are strictly lower triangular matrices, equation (24) is equivalent to

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

Here the empty sums are zero.

Now we identify any matrix V = (v_{ij}) ∈ ℝ^{(n−l)×l}, l + 1 ≤ i ≤ n, 1 ≤ j ≤ l, with the long vector

V = (v_{l+1,l}, v_{l+2,l}, …, v_{n,l}; v_{l+1,l−1}, v_{l+2,l−1}, …, v_{n,l−1}; …; v_{l+1,1}, v_{l+2,1}, …, v_{n,1})^{T}

∈ ℝ^{l(n−l)}. (26)

That is, we order v_{ij} by the following principles. For any two elements of v_{i_{1}j_{1}} and v_{i_{2}j_{2}} of V, v_{i_{1}j_{1}} precedes v_{i_{2}j_{2}} if and only if one of the following conditions holds:

(1) j_{1} > j_{2}

(2) j_{1} = j_{2} and i_{1} < i_{2}.

From equation (25) we note that q_{ij} depends only on p_{ik} for which j ≤ k ≤ l and p_{kj} for which l + 1 ≤ k ≤ i. Hence, according to our ordering q_{ij} depends only on p_{ij} and all the elements that precede it. This means that if we identify ℝ^{(n−l)×l} to ℝ^{l(n−l)} with the above ordering, the matrix representation of D + L is a lower triangular matrix in ℝ^{l(n−l)×l(n−l)}. So equation (23) is easy to solve. Henceforth we shall use the ordering defined by relation (26).

We note that the cost of calculating \overline{P}_{k+1} from \overline{P}_{k} by using Algorithm 10 is exactly the same as calculating P_{k+1} from P_{k} using Algorithm 5 except calculating \overline{P}_{1} from \overline{P}_{0} = 0, which is more expensive than calculating P_{1} from P_{0} = 0.

The relation between Algorithm 5 and Algorithm 10 is essentially the same as the relation between the Jacobi method and the Gauss–Seidel method in solving linear systems iteratively. It is well known that in many applications the Gauss–Seidel method is faster than the Jacobi method (see [2], p. 324). In part 2 we will see numerically that Algorithm 10 is faster than Algorithm 5.

However, to establish a rigorous convergence criterion of Algorithm 10 we need to compute \overline{δ} = ||(D + L)^{−1}||_{2}^{−1}, which is prohibitively expensive to calculate. So we suggest the following algorithm.

\underline{Algorithm 11}.

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

(2) For k = 0, 1, …, \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 ≥ 1), go to (3); otherwise go to (2).

(3) For j = k, k + 1, …, \overline{P}_{j+1} = Φ(\overline{P}_{j}). \underline{Remark 11}. In step (2) the condition ||\overline{P}_{k+1}|| > 2γ/(δ − ϵ) means that \{\overline{P}_{k}\}_{k=0}^{∞} does

not converge to P^{∗}, and the condition ||\overline{P}_{k+1} − \overline{P}_{k}|| > ρ||\overline{P}_{k} − \overline{P}_{k−1}|| (k ≥ 1) means that \{P_{k}\} generated by Algorithm 5 will converge to P^{∗} faster than \{\overline{P}_{k}\} generated by Algorithm 10.

\underline{Remark 12}. The convergence of the sequence of \{\overline{P}_{k}\}_{k=0}^{∞}, generated by Algorithm 11, is guaranteed by step (3) together with Theorem 9. So Algorithm 11 enjoys possible faster convergence of Algorithm 10 as well as the guaranteed convergence of Algorithm 5.

\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. G. W. Stewart, “Error Bounds for Invariant Subspaces of Closed Operators,” SIAM J. Numer. Anal., 8 (1971), 196–208.

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