Introduction. Let us consider a system of linear algebraic equations

K. M O S Z Y ´ N S K I (Warszawa)

CONCERNING DECOMPOSITION OF A SYSTEM OF LINEAR ALGEBRAIC EQUATIONS

Introduction. Let us consider a system of linear algebraic equations

(1) Ax = b,

where A is an N ×N real, invertible matrix. In [1] a method of decomposition of (1) was proposed. The purpose of such a decomposition is to enable parallelization of the algorithm, and if possible to make the problem better conditioned.

Let R = U A − AU . The general idea of the method mentioned above is based on the following observation: if an N × N matrix U of rank r < N commutes sufficiently well with A, i.e. R is sufficiently small , then U defines an approximate decomposition of (1).

Let U = QF , where Q is an N × r matrix and F is an r × N matrix, both of rank r. In [1] it is proposed to replace (1) by one of following systems, which can be solved by iteration:

(2) Q ^T AQy n+1 + Q ^T RQy n + Q ^T RSz n = Q ^T U b, S ^T ASz n+1 − S ^T RQy n − S ^T RSz n = G(I − U )b, or

(3) AQy n+1 + F RQy n + F RSz n = F U b, GASz n+1 − GRQy n − GRSz n = G(I − U )b,

where I − U = SG with an N × s matrix S and an s × N matrix G, and, in general, N − r ≤ s ≤ N . Moreover, x n = Qy n + Sz n converges to the solution x = A ⁻¹ b of the system (1).

We may easily transform (2) and (3) to a more convenient form not containing R (see [1]):

1991 Mathematics Subject Classification: 65F10, 65F35.

Key words and phrases: decomposition of a matrix, nearly commuting matrices, paral- lelization.

[191]

(4) Q ^T AQv n+1 = Q ^T U r n , S ^T ASw n+1 = S ^T (I − U )r n , or

(5) F AQv n+1 = F U r n , GASw n+1 = G(I − U )r n , where

v n+1 = y n+1 − F x _n , x n = Qy n + Sz n , w n+1 = z n+1 − Gx _n , r n = b − Ax n .

If U is a projector, i.e. if U ² = U , each of the systems (4) and (5) contains exactly N equations (s = N − r), hence such a choice is preferable.

This paper concerns the following problem:

Given U = QF , where Q and F are N ×r and r×N matrices respectively, both of rank r ≤ N , we have to construct an N × N matrix V , satisfying the following conditions:

1. rank(V ) = rank(U ) = r;

2. V ² = V ;

3. Im(V ) = Im(U );

4. If at least one of the processes (4) and (5) converges and R is suffi- ciently small , then after replacing U by V , at least one of (4) and (5) will converge as well.

The matrix V

Lemma 1. Let U be an N × N matrix of rank r ≤ N . Assume that there are r linearly independent columns u p

, . . . , u p

of U and r linearly independent columns w _p

1

, . . . , w p

of U ^T such that (w _p

i

, u p

) 6= 0 for i = 1, . . . , r. Then there exist four matrices Q, Q ⁰ , F , F ⁰ of dimensions N × r, N × r, r × N , r × N respectively, such that

(6) U = QF, U ^T = Q ⁰ F ⁰

and

(7) Q ^T Q ⁰ = Q ^0T Q = I r .

P r o o f. Observe that in this case a kind of Gram–Schmidt process of biorthogonalization can be applied to the double system of vectors u p

, . . . . . . , u p

, w p

, . . . , w p

.

We start with

u p

= γ 1,1 q 1 , w p

= γ _1,1 ⁰ q ⁰ ₁ , γ 1,1 γ ⁰ _1,1 = (w p

, u p

) 6= 0,

and then we proceed with the formulas u p

=

k

X

j=1

γ k,j q j , γ k,j = (q _j ⁰ , u p

), j = 1, . . . , k − 1,

w _p

k

=

k

X

j=1

γ _k,j ⁰ q _j ⁰ , γ _k,j ⁰ = (q j , w _p

k

), j = 1, . . . , k − 1, γ k,k γ _k,k ⁰ = (w p

, u p

) −

k−1

X

j=1

γ k,j γ ⁰ _k,j

for k = 1, . . . , r. In this way we get u s =

r

X

j=1

γ s,j q j , γ s,j = (u s , q ⁰ _j ),

w s =

r

X

j=1

γ _s,j ⁰ q ⁰ _j , γ _s,j ⁰ = (w s , q j ),

for all s = 1, . . . , N , j = 1, . . . , N . The last formulas can be written in the form

U = QF and U ^T = Q ⁰ F ⁰ ,

where Q and Q ⁰ are the N × r matrices with columns q j and q ⁰ _j respec- tively, and F and F ⁰ are the r × N matrices of the coefficients γ i,j and γ ⁰ _i,j , respectively.

Assume now that the decompositions from Lemma 1: U = QF and U ^T = Q ⁰ F ⁰ are possible, and are given. Define

V = QQ ^0T and R ⁰ = V A − AV.

Proposition 1. V is a projector.

P r o o f. V V = QQ ^0T QQ ^0T = QI r Q ^0T = V .

Since V is a projector of rank r, I − V is a projector of rank N − r.

Hence I − V = S ⁰ G ⁰ , where S ⁰ and G ⁰ are N × (N − r) and (N − r) × N matrices respectively. This decomposition may be obtained for example by usual Gram–Schmidt orthogonalization, applied to the columns of I − V .

Proposition 2. U = U V = V U .

P r o o f. We have V U = QQ ^0T QF = QI r F = QF = U . Moreover,

U V = (Q ⁰ F ⁰ ) ^T QQ ^0T = F ^0T Q ^0T QQ ^0T = F ^0T I r Q ^0T = (Q ⁰ F ⁰ ) ^T = U .

Proposition 3.

R(I − V ) = U R ⁰ , (I − V )R = R ⁰ U.

P r o o f. By Proposition 2 it follows that R = U A − AU = U V A − AU V ; since AU = U A − R, we have

R = U V A − U AV + RV = U R ⁰ + RV

and so U R ⁰ = R(I − V ). Similarly, R = U A − AU = V U A − AV U , and U A = AU + R, hence

R = V AU + V R − AV U = R ⁰ U + V R, whence R ⁰ U = (I − V )R.

Proposition 4. Q ^0T R ⁰ Q = 0.

P r o o f. Observe that R ⁰ = V A − AV = A(I − V ) − (I − V )A. This yields

V R ⁰ V = V A(I − V )V − V (I − V )AV = 0,

because V is a projector and V (I − V ) = (I − V )V = 0. On the other hand, 0 = V R ⁰ V = QQ ^0T R ⁰ QQ ^0T and 0 = Q ^T V R ⁰ V Q ⁰ = Q ^T QQ ^0T R ⁰ QQ ^0T Q ⁰ . Now, Q ^T Q and Q ^0T Q ⁰ are the Gram matrices of the bases q 1 , . . . , q r and q ⁰ ₁ , . . . , q ⁰ _r , and hence are invertible. Finally, we deduce that Q ^0T R ⁰ Q = 0.

Proposition 5. G ⁰ R ⁰ S ⁰ = 0.

P r o o f. Since V is a projector, we have

(I − V )R ⁰ (I − V ) = (I − V )(V A − AV )(I − V ) = 0, because (I − V )V = V (I − V ) = 0. Therefore

S ⁰ G ⁰ R ⁰ S ⁰ G ⁰ = 0 and

S ^0T S ⁰ G ⁰ R ⁰ S ⁰ G ⁰ G ^0T = 0.

We conclude that G ⁰ R ⁰ S ⁰ = 0, the Gram matrices S ^0T S ⁰ and G ⁰ G ^0T being invertible.

Proposition 6.

R ⁰ Q = (I − V )RF ^T (F F ^T ) ⁻¹ = O(R).

P r o o f. From Proposition 3, (I − V )R = R ⁰ QF , and hence

R ⁰ Q = (I − V )RF ^T (F F ^T ) ⁻¹ .

Proposition 7. If F Q is invertible (that is, U is in some sense close to a projector ), then

Q ^0T R ⁰ = (F Q) ⁻¹ (Q ^T Q) ⁻¹ Q ^T R(I − V ) = O(R).

P r o o f. By Proposition 3, R(I − V ) = U R ⁰ , and by Proposition 2, U = U V = QF QQ ^0T . Hence R(I − V ) = QF QQ ^0T R ⁰ , which implies

Q ^T R(I − V ) = Q ^T QF QQ ^0T R ⁰ . The assertion follows by invertibility of F Q and Q ^T Q.

Theorem 1. Assume that the hypotheses of Lemma 1 are satisfied , the matrix U depends continuously on R, where R = U A − AU , and F Q is invertible for R small. Then the process (5), with U replaced by V , converges for R small enough. This process can now be written as follows:

Q ^0T AQv n+1 = Q ^0T r n , G ⁰ AS ⁰ w n+1 = G ⁰ (I − V )r n , (8) v n+1 = y n+1 − Q ^0T x n , x n = Qy n + S ⁰ z n ,

w n+1 = z n+1 − G ⁰ x n , r n = b − Ax n .

P r o o f. Let us return to the equation (3), equivalent to (5). Now, if U is replaced by V , in view of Propositions 1–7, the equation (3) admits the following form:

Q ^0T AQy n+1 + Q ^0T R ⁰ S ⁰ z n = Q ^0T b, G ⁰ AS ⁰ z n+1 − G ⁰ R ⁰ Qy n = G ⁰ (I − V )b.

By Propositions 1–7, the coefficients of all terms containing y n and z n are of order O(R); hence the convergence follows by standard arguments.

Case of A symmetric. Put now V = QQ ^T , R = U A − AU , R ⁰ = V A − AV , and U = QF with Q ^T Q = I r . A decomposition of this kind may be obtained for example by application of the Gram–Schmidt process to the columns of U .

Proposition 8. V is an orthogonal projector.

P r o o f. V V = QQ ^T QQ ^T = QI r Q ^T = V . Moreover, V ^T = (QQ ^T ) ^T = QQ ^T = V .

Since I − V is of rank N − r, we may decompose (by the Gram–Schmidt process)

I − V = S ⁰ G ⁰ , where S ^0T S ⁰ = I N −r . Proposition 9. If A = A ^T , then R ^0T = −R ⁰ .

P r o o f. We have R ^0T = (V A − AV ) ^T = A ^T V ^T − V ^T A ^T = AV − V A

= −R ⁰ .

Proposition 10. R ⁰ Q = (I − V )RF ^T (F F ^T ) ⁻¹ = O(R).

P r o o f. We have Q ^T U = Q ^T QF = F , hence U = QQ ^T U = V U , R = U A − AU = V U A − AV U = V (AU + R) − AV U = R ⁰ U + V R, and so (I − V )R = R ⁰ U = R ⁰ QF . Since F F ^T is invertible, we get R ⁰ Q = (I − V )RF ^T (F F ^T ) ⁻¹ .

Proposition 11. If A = A ^T , then Q ^T R ⁰ = −(F F ^T ) ⁻¹ F R ^T (I − V ) = O(R).

P r o o f. We have

R ⁰ Q = (I − V )RF ^T (F F ^T ) ⁻¹ , whence by Proposition 9,

−Q ^T R ^0T = Q ^T R ⁰ = −(F F ^T ) ⁻¹ F R ^T (I − V ).

Proposition 12. If A = A ^T , then Q ^T R ⁰ Q = 0.

P r o o f. Observe that (I − V )Q = Q − QQ ^T Q = Q − Q = 0 and Q ^T R ⁰ Q = −(F F ^T ) ⁻¹ F R ^T (I − V )Q = 0.

Proposition 13. S ^0T R ⁰ S ⁰ = 0.

P r o o f. We have

(I − V )R ⁰ (I − V ) = (I − V )(V A − AV )(I − V )

= (I − V )V A(I − V ) − (I − V )AV (I − V ) = 0 because V (I − V ) = (I − V )V = 0, where V is an orthogonal projector.

Since I − V is symmetric, it follows that I − V = (I − V ) ^T = G ^0T S ^0T and (I −V )R ⁰ (I −V ) = G ^0T S ^0T R ⁰ S ⁰ G ⁰ . Observe that G ⁰ G ^0T is invertible, whence G ⁰ (I − V )R ⁰ (I − V )G ^0T = 0, which completes the proof.

Theorem 2. Assume that A = A ^T , and that U = QF , where Q ^T Q

= I r , depends continuously on R = U A − AU . Then the process (4), with U replaced by V = QQ ^T , which is now of the following form:

Q ^T AQv n+1 = Q ^T r n , S ^0T AS ⁰ w n+1 = S ^0T S ⁰ G ⁰ r n , (9) v n+1 = y n+1 − Q ^T x n , x n = Qy n + S ⁰ z n ,

w n+1 = z n+1 − G ⁰ x n , r n = b − Ax x , converges for R small enough.

P r o o f. We recall the equation (2), equivalent to (4), which now takes the form

Q ^T AQy n+1 + Q ^T R ⁰ S ⁰ z n = Q ^T (I − V )b,

S ^0T AS ⁰ z n+1 − S ^0T R ⁰ Qy n = G ⁰ (I − V )b.

From Propositions 8–12 it follows that the terms containing y n and z n are of order O(R); hence, for R small the convergence follows by standard ar- guments.

Example. Assume that an N × N matrix A and an M × M matrix C, with M < N , are two finite-dimensional approximations of a certain linear operator. For simplicity, assume both matrices A and C to be symmetric and invertible.

Let

p : R ^M → R ^N and r : R ^N → R ^M

be linear extension and restriction operators, respectively (see [2]). Put U = pCr : R ^N → R ^N .

If p and r are properly chosen (see [2]), then we may expect that R = U A − AU will be small for sufficiently large N and M , M < N . We may also expect (at least in certain situations—see the Laplace operator for ex- ample), that in general the matrix C will correspond to a lower part of the spectrum of the original operator than the matrix A. This phenomenon may be explained as follows: approximation on a rough grid in general does not allow passing higher frequency oscillations.

We may apply our algorithm (9) to the matrix A and U . Application of the Gram–Schmidt process to the columns of the matrix pC will give pC = QΓ with Q ^T Q = I M . Hence we get

pCr = QF

with F = Γ r. We can construct in an arbitrary way an N × (N − M ) matrix Q in order to get an N × N orthogonal matrix e

[Q| e Q].

We have V = QQ ^T and

I − V = [Q| e Q][Q| e Q] ^T − QQ ^T = QQ ^T + e Q e Q ^T − QQ ^T = e Q e Q ^T . In other words, S ⁰ = e Q and G ⁰ = e Q ^T .

Now the system (9) can be written in the following form:

Q ^T AQv n+1 = Q ^T r n , Q e ^T A e Qz n+1 = e Q ^T r n ,

v n+1 = y n+1 − Q ^T x n , x n = Qy n + e Qz n ,

w n+1 = z n+1 − e Q ^T x n , r n = b − Ax n .

References

[1] K. M o s z y ´ n s k i, On solving linear algebraic equations with an ill-conditioned matrix , Appl. Math. (Warsaw) 22 (1995), 499–513.

[2] R. T e m a m, Numerical Analysis, Reidel, 1973.

KRZYSZTOF MOSZY ´ NSKI

DEPARTMENT OF MATHEMATICS, COMPUTER SCIENCE AND MECHANICS UNIVERSITY OF WARSAW

BANACHA 2

02-097 WARSZAWA, POLAND E-mail: KMOSZYNS@MIMUW.EDU.PL

Received on 28.9.1994