An algorithm for the inversion of classes of sparse matrices

An algorithm for the inversion of classes of sparse matrices D J U R O M. M I ~ L J E N O V I ( ~

Mathematical Institute, Beograd, Yugoslavia

INTRODUCTION

REVIEW OF THE CHOLESKY METHOD

Many problems in science are solved using the system of N linear equations with N unknown variables

Let the matrix A be symmetric and non-singular. Then there exists a lower triangular matrix B and a diagonal matrix D which satisfies matrix equation

Ax = b

(1)

It can be done by using two groups of methods. The first group is based on direct elimination procedures, like Gauss method. The second requires the inverse of the matrix A. Thus, if one has to solve the system in equation (1) only once or a few times, then some of the direct methods could be applied. Very often we are concerned with numbers of right-handed side vectors b and it could then be better to use an inverse matrix approach, which is concurrent with the direct methods if the problem in equation (1) has to be solved for approximately N right-hand sides. Usually we have not such a large number of different cases but for special structures of the matrix A, methods for inversion becomes applicable. Sparseness of special structures (which is described later on) and/or bandedness reduce the number of multiplications to be performed to the level adoptable for the wide class of problems. STATEMENT OF THE PROBLEM

Let the physical system consist of N objects with mutual interactions described by symmetric and diagonal-dominant matrix A. The first property means that the mutual influence of any two objects is the same in both directions. The second property finds that the influence of one object on itself is greater than or equal to its influence on all other objects together. The objects in BE analysis are elements and the matrix A is slightly more general, i.e. this matrix is nearly symmetric as the consequence of the collocation method used. The term nearly symmetric is used for the matrix A with elements ai/close to the elements a/i, for all i and]; so, we do not require equality. The standard approach is to solve the system of linear equations obtained with a suitable numerical method, like the Gauss elimination. In applications one needs to solve a number of systems because of different load cases, insufficient accuracy of results obtained etc. So, a more flexible (and time consuming) matrix inversion method becomes a good alternative. Here, we develop one efficient matrix inversion algorithm for the class of the problems previously described. The Cholesky decomposition method is found to be the best one because the theorems develop naturally to cover assumptions imposed on the problem. 0141-1195/83/010039-04 $2.00 © 1983 CML Publications

A = BDB T

(2)

If the matrix A is positive definite, one can choose D = I and decomposition becomes simpler A = BB T

(3)

These procedures are usually applied for solving the system of linear equations Ax = b

(4)

The numerical technique for solving equations (4) based on decomposition equation (2) is called the Cholesky method; similarly, using the decomposition equation (3) gives us the well-known square-root method. Here, we demonstrate the applicability of these methods for inversion of the matrix A. It can be used in an efficient way for a class of sparse and/or banded matrices, as the positions of zero-elements in the decomposition matrices can be preserved. In the following sections various cases of the Cholesky method are presented. We assume the matrix A to be positive definite, for simplicity. The more general case of a non-singular matrix can be treated in a similar way, using decomposition equation (2); note that elements of matrix D can be chosen as plus or minus 1 by using the principle of inertia. The Einstein summation convention is used together with details about running indices, if necessary.

Theorem 1 (standard case): Let the matrix A be symmetric A =[ai/],

ai/= a/t

(5)

Also, suppose the lower triangular matrix B is given B = [bt/],

bi] = 0

] > i

(6)

where all diagonal terms are positive. Then if the matrix A is positive definite there exist the unique decomposition A = BB T

(7)

Proof." Here we demonstrate the way of calculating the elements of matrix B. So, it will be easy to validate the uniqueness of the solution obtained. The condition in equation (7) in the scalar form is ai/= bikb/k

k = 1. . . . . N

Adv. Eng. Software, 1983, Vol. 5, No. 1

(8)

39

Using the property of indices of the matrix B the range for index k is

or in the matrix element form a i / = bigbik

k ~< min(i,/)

(231)

k = I..... N

(9)

Obviously, it is sufficient to deal with the elements from the lower triangular portion of the matrix A, i.e. only for j <~ i. Applying the previous restriction, the range of the index k is further reduced to

It is sufficient to work with the lower triangular portion of the matrix A, because of symmetry. The starting expression can be transformed by using the property equation (21) of matrix B and the summation is reduced to the range

k
k = 1. . . . . min(i,/)

(10)

So, the following basic formula is obtained aij = bikb/k

(11)

k <~/

Hence, for i = 1 we have f = 1 and ala = b]l

(12)

For positive definite matrices all diagonal minors are positive; thus, a ~ is positive number and

Now, let the row index i be fixed and let the column index j take values from the set 1. . . . . i - m , i.e. we are only interested in zero elements of the matrix A at the row i. In this case the summation can be further reduced and we obtain the formula aij = bikbik = 0

Now, let us perform calculations for the fixed row i, from left to right, i.e. for j = 1, 2 . . . . . i--1. For j = 1 we have

(25)

k = 1. . . . . i - - m

For/" = 1 we have

(13)

bH = SQR T(alx)

(24)

bllbil = 0

(26)

i> m

and the Main Property of the Cholesky decomposition gives us bil = 0

(27)

i> m

(14)

all = b i l b l l

The expression (25) has one summand less. Hence and. consequently, we find

For 1 < / <

i we have aij = b i k b j k

bikb/k = 0

(15)

bil = a i l / b l l

(16)

k <~j

bij = 0 m
(17)

The new value bi/ is expressed in terms of all known quantities bij = (aij -- b i m b j m ) / b j j

m
(18)

F o r / = i one obtains a i i = bikbik = bi~

k <~ i

(19)

or, equivalently, in developed form bii = S Q R T ( a i i -- b2m)

m < i

(20)

The condition of positive definitiveness allows us to perform calculations in the previous expression, which completes the proof. T h e o r e m 2 (banded case): If the symmetric, banded and positive definite matrix A, with semi-band-width equal to m is transformed using the Cholesky method then the corresponding triangular Cholesky matrix B is with the same semi-band-width. R e m a r k : The so-called Main Property of the Cholesky decomposition states that a i i > O. Proof." Let the matrix B = [bi/] be lower triangular, i.e. bo=O

/>i

(21)

40

Adv. Eng. Software, 1983, Vol. 5, No. 1

(22)

/ = 1. . . . . i - - m

Q.E.D. (29)

Theorem 3 (variable band-width): Let the symmetric and positive definite matrix be A with m i zeros at the beginning position in the row corresponding to i. Let the lower triangular matrix B correspond to the Cholesky decomposition of the matrix A. Then matrix B has m i zeros at the same positions as the zero elements of the matrix A. Proof." It is easy to see that the proof of Theorem 2 gives all necessary information for this, more general, case. In fact, one needs to perform the same calculations for each row of the matrix A, in this case having a variable bandwidth, and the conclusion becomes apparent. Now, we can state the main conclusion of interest as the Theorem. Theorem 4 (general case): Let A be the symmetric and positive definite matrix of a special structure (Fig. 1) with

The assumptions about matrix A are sufficient for the existence of the Cholesky decomposition matrix B. So, let the matrix A be represented as the product of the Cholesky matrix B and its transpose B T, i.e. A = BB T

(28)

Successive usage f o r / = 2, ] = 3 . . . . . / = i - - m gives us the corresponding elements of the matrix B, which are equal to zero. So, we obtain the desired result for the f'Lxed row i

or in the developed form aij = bimbjm + bijbjj

k = 2 ..... i--m

Figure 1.

a semi-band-width equal to m and the lower triangular matrix of dimensions m by m. Using the Cholesky method it is possible to split the matrix A into the pair B and B T of triangular matrices of the form presented in Fig. 2, with positive diagonal elements, where the matrix B has the same semi-band-width as the matrix A. Remark: Test examples assures us that the last m rows of the matrix B are really fully populated. Thus, the previous Theorem is, in terms of finding locations of elements equal to zero, the best possible.

which is the second order process (with quadratic convergence). In matrix form, this formula gives us E~ as the first approximation to A-~ for Ek+l = Ek(2/-- AEk)

C = LL ~

(30)

Fourth step: invert the C matrix by using L in order to obtain the approximation for B -~, or, consequently for A -1. Fifth step: use some iterative procedure to improve the approximation to A -1 obtained. The previous algorithm needs some additional explanations. (a) Symmetricization of the matrix A is not necessary, i.e. we assume that the matrix A is symmetric. In that case B = A. (b) If we symmetricize A we obtain matrix B as B = (A + AT)/2

(31)

(c) An element of the matrix B is small if the following sparseness condition is satisfied

Ibi/l
]
S E ( 0 , 1)

D1 = I - - AE1

Figure 2.

(35)

Later on, an error is Dk+l = D~c+I = D ~ = D~k

(f)

kEN

(36)

sO, the convergence is very fast if D1 is close to the zero matrix. Usually, it is enough for one or two iterations to obtain the inverse matrix with a high accuracy. Our first trial was to perform the simplified version of the algorithm by using the Gauss method for the inverse of the sparse matrix C. The direct elimination phase was performed successfully because of the sparseness of the matrix C but the substitution phase became cumbersome as the sparseness was destroyed and we were forced to work with fully populated matrices. So, we can not suggest this approach.

NUMERICAL CALCULATIONS For our applications the matrix C is the structure given by Fig. 1. This structure is a consequence of taking into account only the influences between elements close to one another; the rest of the influences are equal to zero. In order to preserve simplicity, let us assume that the diagonal part of matrix C has semid-band-width equal to rn and that the lower triangular matrix of order m x m is located at the left-hand corner of the matrix. Later on, we shall relax this assumption by applying Theorem 3. Here we want to take advantage of the sparseness of matrix C for numerical calculations. The easiest way is to split this matrix into four blocks. The Cholesky matrix can be split in the same way and one obtains:

(32)

Usually, S E ( 0 , 0.2) and this corresponds to the semiband-width of the C matrix. (d) We demonstrate below that the matrix L is also sparse, with a special structure. (e) In order to improve the approximation of the matrix A -1 one may use the matrix equivalent of the Newton iterative formula for the reciprocal of the number a, i.e. Xk+l = xk(2 -- aXk)

(34)

The error for the first approximation is given by

MATRIX INVERSION ALGORITHM In order to perform the inversion of the matrix A as fast as possible, we propose the following procedure: First step: symmetricize the matrix A, which gives us matrix B. Second step: transform the matrix B into a sparse matrix C by making all small elements in B equal to zero. Third step: find the Cholesky decomposition of the matrix C, i.e.

k EN

(33)

tC21 C22J where Cn and C22 are symmetric and banded matrices, U, V are lower triangular matrices and C12 = C21. Using the symmetry of matrix C we have: Ctl = UU T Cl2 = UB T

(38)

C22 = BBT + VV T The first of these equations allows us to calculate the matrix U. It can be done efficiently by using Theorem 2 because of the banded nature of C1~. The second equation is necessary for determining matrix B. Note that matrix U is banded and triangular (Theorem 2) and with a fully populated lower triangular inverse, but the matrix C12 is sparse; so, a lot of elements in B are equal to zero. With our assumption about bandness it is possible to demonstrate that the only non-zero elements of B are located in the lower m rows at the righthand comer in the upper triangular form (Fig. 3). The remaining equation is suitable for the evaluation of matrix V. The matrix C22 is symmetric and so is the product BBT, with the special structure depicted in Fig. 4. The off-diagonal blocks of this matrix product had

Adv. Eng. Software, 1983, Vol. 5, No. 1

41

APPLICATION POSSIBILITIES Here, we present the results for decomposition of the symmetric and banded matrices, which are suitable for processes of sub-decomposition to be applied. It can be used independently as in finite elements analysis bearing in mind the sparseness inside the band, which may affect applicability. Also, in BE applications, the important class of problems with planes of symmetry can be treated. The sparse matrix C, an approximation of the corresponding BE matrix A, is really banded and theorem 2 is a powerful tool for handling this type of problem successfully. The main difference, in comparison with equations (37)-(41) is that the matrix A22--BB T is sparse, with the same band-width as matrix A22 , SO, the structure of the matrix is preserved. This feature is valuable for calculations on smaller computer systems because of smaller storage requirements. Note, we used the fact that B obtained here has one upl~r triangular block in the upper right-hand corner and BB" consists of a fully populated matrix in the upper lefthand comer.

@

Figure 3.

@

SPARSENESS OF THE MATRIX

I/

111

The Cholesky matrix is very sparse. For the symmetric N x N matrix of the type in Fig. 1, with semi-band-width equal to m = sN, s E (0, 0.5), we have to calculate the following number of elements for the lower triangular matrix

Figure 4.

destroyed the banded nature of the matrix C22 and the resulting Cholesky matrix V is not banded but with the structure as presented in Fig. 2. BBT can also be calculated in the following way. The second equation of equation (38) gives us: B T = U-1C12

M(s) = 4s(1 - s)

M(O.05) = 0.19

(39)

M(0.10) = 0.36

T -1 C12CllC12

(41)

Note that C11is banded but C~I is fully populated; also, the matrix C12 is sparse. So, it is possible to fred which elements of BBT are non-zero and calculate only them from the right-hand-side expression, taking into account the sparseness of C12 also. It is important to note that the matrix C22--BB T has not the same form as the starting one (from Fig. 1); the difference between these two matrices is in the form of non-diagonal blocks, which are now fully populated matrices. Applying the same procedure for the later matrices we Fred that the corresponding B matrix in the second step of calculations has the same form as the starting matrix, i.e. the matrix from Fig. 1. Hence, this slightly more general case can be treated in exactly the same way. The importance of the splitting procedure is in the possibility of performing calculations in a step-by-step fashion; so, we have an opportunity to solve large problems using small computers.

42

Adv. Eng. Software, 1983, Vol. 5, No. 1

(43)

M(0.15) = 0.51 (40)

Thus, the product BB T is given by: BBT=

(42)

So, for standard cases we have s = 0.05, s = 0.1, s = 0.15 and the corresponding values of this quadratic function are:

and transposing this equation we obtain: B = c T ( u T ) -a

s E (0, 0.5)

These numbers, multiplied by 100 gives us the percentage of element numbers from the symmetric portion of the starting matrix which have to be calculated. ACKNOWLEDGEMENTS

The author is grateful to Dr C. A. Brebbia, University of Southampton, UK for his useful discussions and comments about these results. He also thanks the Mathematical Institute, Beograd, Yugoslavia, for all their support. BIBLIOGRAPHY

1 Bertolino,M. Numerical Analysis, Naucna knjiga, Beograd, 1977. 2 Brebbia, C. A. The Boundary Element Method for Engineers, Pentech Press, London, 1978. 3 Brebbia, C. A. and Walker, S. Boundary Element Technique in Engineering, Newnes-Butterworth, 1980. 4 Brebbia, C. A. (ed.) Progress in Boundary Elements, Vol. 1, Pentech Press, London, 1981. 5 Brebbia, C. A. (ed.) Boundary Element Methods, Springer Verlag, Heidelberg, 1981. 6 Brebbia,C. A. Private communication. 7 Gantmacher, F. R. The Theory of Matrices, Chelsea Publishing Co., New York, 1959.