- Email: [email protected]
Iterative inversion of single exit compartmental matrices
OOIO-4825/79/lOO-0317 $02.00/O.
Compur. Biol. Med. Vol. 9, pp. 317-330. @ Pergmon Press Ltd. 1979. Printed in Great Britain.
ITERATIVE INVERSION COMPARTMENTAL DAVID
OF SINGLE EXIT MATRICES
H. ANDERSON
Department of Mathematics, Southern Methodist University, Dallas, Texas 75275 and Department of Medical Computer Science, The University of Texas Health Science Center, Dallas, Texas 75235, U.S.A. (Received
25
July 1978;
in
revised form 5 February 1979)
Abstract - This article develops an iterative method for solving a system of linear algebraic equations in which the coefficient matrix is a compartmental matrix of the single-exit type. Matrix properties, invertibility, ill-conditioning, algorithm acceleration and error analysis are discussed. Matrix Single-exit compartmental systems Simultaneous linear algebraic equations Numerical iterative methods Error analysis and invertibility Acceleration of convergence rate of convergence
1. INTRODUCTION In the mathematical modeling and simulation of phenomena in such areas as biology, medicine and chemistry, there naturally arises an important class of square matrices in which each matrix has the following properties: (a) each off-diagonal element of the matrix is nonnegative; (b) each diagonal element is negative ; (c) all columns except one add to zero ; (d) that exceptional column sums to a negative number. A simple example is -0.3 0.1
K=
0.1 -0.1
0.2
0.0
0.0
0.0
0.0
0.2
0.2
0.2
-0.4 0.0
0.1 -0.5
J
One way such matrices come about in practice is through linear compartmental models associated with radioactive tracer studies dealing with human physiology [ 11. For example, suppose there are 3 identifiable compartments in a certain clinical study and of these there is only a single compartment, the first, through which exit to the outside -environment is possible (this “single exit” is the essential characteristic of the system for our study). Let x,(r) be the amount of tracer in the ith compartment. Assuming the kinetics of the system are basically linear, the mathematical model for the system dynamics is taken to be (a s dx/dt) ~2, = - k,,x,
- kS1xl - k,,x,
+ k12x2 + k,,x,
Jz =k2,x1
+ kz3x3 -k,zxz
-kszxz
z& = k,,x,
+ k3,x,
- kzJxg + b,
- k,,x,
+ b,
+b,
in which k, is the nonnegative fractional transfer coefficient from compartment j to compartment i, k,, is the positive excretion rate constant, and b,(t) is some input function [2, p. 481. This system can be written in the vector-matrix form i=Kx+b 317
DAVID H. ANDERSON
318
where
-621 + k,, + kc,,) K=
k 21
k
k,2 -&2
k 31
+
k::
k32)
k 32
-h
+
k23)
This matrix has properties (a)-(d) listed above. 2. PRELIMINARIES Following the general development adopted in [2, p. 481 or [3], let K z [k,] be an n x n real matrix in which k, 2 0 if i # j, and for all j = 1,2,. . . , n, l#j
where k, are nonnegative numbers and the summation is over all i = 1,2, . . . , n, except i = j. To be nontrivial, we assume that k,, < 0 for all j. We now delineate the important subclass of interest to our study. Suppose that in addition to the above properties, the quantities k, are taken to be zero for all j = 1,2,. . . , n with the single exception ofj = m, for which we have k,, > 0. Motivated by the introductory example, such a matrix K possessing these properties will be called a single-exit compartmental (SEC) matrix.
There are at least 2 interesting problems concerning such matrices. First, that of determining any unknown entries of the matrix from observation of certain compartments of the system in time (the “inverse” or “system identification” problem), and second, that of solving a linear algebraic system having a SEC matrix as the known coefficient matrix, or equivalently, inverting a SEC matrix (e.g. equilibrium or steady-state problems associated with single exit compartmental systems give rise to these types of algebraic equations). The present paper is concerned with the second problem, that of obtaining a general theory and method of solution of the inversion problem. Given any n x n SEC matrix K and n x 1 real vector y, our goal then is to investigate the theory of solution of the linear algebraic system Kx = y. The method of approach is to construct an approximate inverse of K, from this derive conditions under which K is invertible and then design an iterative scheme that will always take us from a good initial guess to the solution x. Knowledge of the proposed computational scheme is appropriate for a variety of reasons. The matrix K is often ill-conditioned and this complication, which causes difficulties in direct methods such as Gaussian elimination, may be overcome by using the iterative technique to be developed. Ill-conditioning of K is definitely indicated since the condition number 11K /I 11 K- ’ 1)may be quite large. Also K- ’ may have large values relative to K, and a small excretion rate k,, is not atypical in biomedical problems (det K --+0 as k,, + 0). Through the iterative approach, upper bounds on the error involved are possible. An important consequence of the approach is the resulting mathematical theory obtained about the structure of K and its inverse matrix. The proposed iterative method may be applied in order to improve an approximate solution already obtained by some other method. Many present modeling efforts yield sparse matrices, and for such matrices the proposed computational scheme is simple and convenient for use, and its storage requirements economical. Probably the only disadvantage of the successive approximation method is that convergence may be slow. However, we are able to develop a method of accelerating convergence which does not complicate the iterative formulas and in fact decreases the scheme’s storage requirements.
3. INVERTIBILITY
OF A SEC
MATRIX
In treating the theory of solution of Kx = y, we should first face the question of invertibility of the matrix and try to ascertain simple conditions governing nonsingularity of
Iterative inversion of single compartment matrices K.
319
There do exist singular SEC matrices ; an immediate example is -3
0 1
0
I
2
-4
1 4 -2 K= [ Let us start with a known invertibility result for genercrl compartmental matrices (matrices K in which k, may be zero or nonzero). Draw the connectivity dioyrum D for the specified compartmental matrix K = [kli] (kij > 0 for i # j means to draw the directed line segmentJ) and let the associated compartmental system be denoted by s E {1,2,...,n} where the integers denote the numbering of the compartments (or vertices) in D. Suppose (renumbering the compartments if necessary) we define the subset T={p,p+l,...,n)(p
of S. This subsystem T of the compartmental system S is called a rrap provided k, = 0 for all i,j such that j 2 p and i < p. In other words, there does not exist a path in D from any compartment in T to a compartment which is in S but outside T. See [4] and [5] for information on the following result. Theorem 1. The compurtmentul mutrix compurtmentul system S has no trups
K is invertible if and only if its ussociuted
Consider the following (SEC) matrix K and its connectivity diagram D (Fig. 1):
K=
1
0
1
1
0
-
0
Exit .
. (1)
0 -1
I
Fig. 1.
Observe that this matrix is reducible so that the standard inversion theory for irreducible diagonally dominant matrices [6, p. 1063 does not apply. However, as can immediately be seen from D, the matrix is nonsingular since the associated compartmental system contains no traps. For the particular subclass of SEC matrices, an alternative necessary and sufficient condition for invertibility can be given (Theorem 4) which is vitally important for the iterative process to be developed in this paper (Section 4). We make use of the result that a square matrix is invertible if and only if there exits an “approximate inverse” of A [7]. To be specific, if for a given n x n real matrix A there is an invertible n x n real matrix B and a matrix norm I] . 11such that (II - BA II < 1, then B is said to be an approximate inverse of A (I is the n x n identity matrix). Note that if B = A- ‘, then I - BA is the zero matrix so that we may consider the size of 1 - BA to be a rough estimate of the inaccuracy of B as an approximation to A-‘. The underlying theme of this paper centers on the construction, for a specified n x n SEC matrix K, of an invertible n x n real matrix J such that (II - JK II1 I 1. Here the l-norm or column norm [8, p. 2123 i IIAII1= max I
(ulil
DAVID H. ANDERSON
320
is a natural choice since much of the structure of the SEC matrix is geared to its column sums. Let a SEC matrix K be given. Since J should look something like K- ’ (if K- ’ exists), then it is desirable to know certain fundamental information about this inverse in order to construct J. After calculating a few such inverses, the following theorem is not difficult to conjecture. All theorem proofs appear in the appendix. Theorem
2. If the SEC matrix K is invertible, then each entry in the mth row of K- ’ is - 1 /k,,,
Judging from the literature, some of the most interesting and important cases of SEC systems are the mammillary or cutenary types [2, p. 531. Theorem 2 allows us to get a direct solution of Kx = y in some of these cases. The crucial observation is that since each entry in the mth row of K- ’ is known to be - l/k,, then x = K- ‘y yields the closed form calculation n
1 Yilkww
Xm = -
i=l
If, in the case of a mammillary system, the number of the central compartment is designated one, then the coefficient matrix K assumes the form where there are nonzero entries only in its first row, first column, and main diagonal. The mth component x, is now known from the above equation and thus x, can be computed from the mth equation of Kx = y. Then for i = 2,3,. * . , m - 1, m + 1,. . . , n, we can read xi = (yi - kiixi )lkii directly from Kx = y. If K represents a catenary system with the usual sequential numbering ofcompartments, then it is a tridiagonal matrix possessing nonzero elements only on its main diagonal, first superdiagonal, and first subdiagonal. Suppose that excretion takes place just from the terminal or nth compartment (m = n). Then x, is known from the closed form (2) and the system Kx = y can be solved by back substitution. For, the component x, _ 1 is found from the equation k,.,- ix.- 1 + k,,x, = Y,, and then the equations ki.i-ixi-i i = n - 1, n - 2,...,
+ kiixi + ki,i+rxi+r
2, can be successively
= Yi3
solved for x._~, x,_ 3,. . . , x1.
In terms ofdesigning J, the conclusion coming out of Theorem 2 is that a reasonable choice for each element in the mth row of J is - l/k,,. There is a second observation that aids in the formation of J, and that is the diagonal dominance of the K matrix suggests that generally the off-diagonal elements of K will be small in comparison to the diagonal elements. Thus in a loose sense, K z ditry{k,},
and so the diagonal
entries
of K- ’ might be close to l/k, ,r
V&,.
..,
ilk,,.
Unfortunately, if the diagonal elements of J are chosen as {l/k,,}, we do not always get )II 1 for any SEC matrix K selected. Instead of taking the option ofchoosing a different norm, we make a slight modification in the choice of the diagonal elements of J which will give us 111 - JKII 1 I 1 for every K. Define the parameter c to be any number satisfying JK 11l I
c 2 max Ikii(. I
Then let every diagonal position of J be filled by - l/c. This discussion motivates the actual definition of the J matrix. For a specified n x n SEC matrix K, let J [or, if we wish, J(c)] be the n x n matrix with - l/k,, as each element in its mth row, the value - l/c as each element on its main diagonal (in other than the (m, m)-position), and zeros elsewhere. For example, consider the 4 x 4 matrix displayed in (1). It is such that
Iterative inversion of single compartment
m =
matrices
321
1, k,, = 1 and c 2 2. If we take c = 2, then the J associated with (1) is -1
-1
-1
0-i J = J(2)=
0
-1 0
0-t
0 0
0 0 0 -9 [ I Let us observe a few things about the .Z matrix. It is always invertible since detJ = (-l)“/k,,c”-’
#O.
It is easy to verify that J-l has -k,, as its (m,m)entry, every other element on the main diagonal is -c, each nondiagonal entry of the mth row is c and zeros elsewhere. In fact, .Z- ’ is itself an SEC matrix. For any SEC matrix K form the associated J matrix according to the above rules. Then define the n x n matrix E = E(c) E I - J(c)K.
We can now prove that 11Z - .ZKII1 I especially the nonnegativity of E. Theorem
1
as well as get some interesting ancillary results,
3. Define the subset of indices H E {j: kmj = O}.
Then )IE (11 < 1 if H is the empty set, and 1)E 111 = 1 if H 1s . nonempty. Moreover, E is a singular reducible non-negative
matrix.
We conclude via Theorem 3 that J is a reasonable approximation of K-’ and thus in a rough sense K z 1-l. This observation provides an interesting “first order” interpretation of the interworkings of a single exit compartmental system. To within this approximation, the only non-zero fractional transfer coefficients are those coming into the mth compartment, k,,
= c,
j = 1,2,. . . , n, j # m,
and thus the system is acting like a mammillary system with excretion from the mth or central compartment [2]. This theorem allows us to derive a number of results. The first of these ‘is the basic invertibility condition, which is stated in terms of the spectrul radius of a square matrix A, p(A) = max{ Irl : r is an eigenvulue of A}. Theorem
4. The SEC matrix
K is invertible if and only if p(E) c 1.
The above theorem tells us that J & an approximate inverse of K whenever K is nonsingular. The result also reduces the invertibility of K to the algebraic problem of demonstrating that every eigenvalue df E has magnitude less than unity. Even though there are presently available excellent methods of computing eigenvalues, the equivalence of Theorem 4 is theoretical in nature and not readily applicable to a particular matrix E. However, based on this fundamental equivalence we can derive simple sufficient conditions on the entries of K to determine its invertibility. The next criterion is especially satisfying from the modeling point of view because it depends only on the disposition of zero and nonzero elements in the matrix K and thus relates to Theorem 1. Upon commencing the modeling process, this disposition is often all that is known or even conjectured about kii. Theorem 5. Recall that H = {j : kmj = 0). Zffor each j E H there is at leust one i 4 H such that kij is nonzero, then K is invertible
Applying Theorem 5 is very easy. Simply strike out all rows with row numbers in H and all columns with column numbers not in H. If in each of the remaining columns there exist at least one nonzero element not eliminated, then K is invertible.
DAVID H. ANDERSON
322
4. THE
FUNDAMENTAL
ITERATIVE
SCHEME
The problem to be considered in this section is that of finding the solution vector x to satisfy the linear system Kx = y where the SEC matrix K will be assumed nonsingular. We will develop an iterative scheme based on our knowledge of the approximate inverse J of K, a theory which has already been demonstrated to be rich in ideas. A method of successive approximation will be set up as follows. The equation Kx = y implies that x = (I - JK)x + Jy = Ex + Jy, a fixed point equation for the true solution x. Thus the computational X(P+“=EX(P’+Jy,
p=o,1,2
scheme
,...)
is suggested. Moreover, because E is small in the sense that p(E) < 1, the successive approximation is guaranteed to provide a sequence that converges regardless of the starting vector [6]. Since x = K-‘y z Jy, we choose to initiate the process at x(“’ E Jy. If the iterates are successively evaluated, then
xcl' =
Jy + EJy
xc” = Jy + EJy + E2Jy
is a sequence converging to EPJy
x= f
(E” = I).
p=o
When the scheme is written in component computation on most machines, we have x(P+ll I
=
jgl
form, which should be used for actual
eijXy'
+
X$O’,
i= 1,2 ,..., n. Since eij = 0 when i = m, we immediately get the mth component solution x in closed form:
(3) of the
n
xm = x’,“’= - i;l Yilkmw Thus the scheme (3) can be reduced in dimension by one by dropping the mth iterative equation and replacing all x2’ by x,. Let us at this point introduce the component-wise partial ordering I (or 2) for real vectors x = (x,,. . . ,x.) and y = (vi,... . ,y,) by x I-y
if and only if x, I yi,
i = 1,. . . ,n,
and similarly for n x n real matrices A = [Uij] and B = [bIj], AI B
if and only if Ui, I bu
i,j = 1,. . . ,n.
Now suppose y I 0, a common situation in biomedical work. Since J I 0 and E 2 0, then Jy 2 0 and EPJy 2 0 for all p = 1,2,... . Hence {x@‘}forms a monotonically increasing sequence that converges unilaterally to the solution x of Kx = y. If there are many y vectors to be used, there may be an advantage in finding K- ’ and then proceeding by the simple multiplication K- ‘y in each case of y to the solution x. One way of generating K- ’ is through an algorithm parallel to that above. The relationship E = I - JK yields the fixed point equation K-‘=EK-‘+J
Iterative inversion of single compartment
for the inverse matrix; thus the computational X p+i =EX,+J,
matrices
323
scheme p=o,1,2
)..‘)
with X, = J is motivated. Because this is a linear iterative process and p(E) < 1, the sequence X,=J 2,
=J+EJ
X2 = J + EJ + E’J
is guaranteed to converge to K-‘.
Moreover, since EPJ 5 0 for all p = 1,2,. . . , we have
This is an improvement of the known result that K-l < 0 [3]. These inequalities on the sequence elements are very useful not only because they yeild such information as x > 0 (in the case when y I 0) and K- ’ I 0, but also because if either iteration is terminated at any step p, then it is known in advance that x(P)is a lower bound for x and X, is an upper bound for K-l. The particular upper bound J of K - 1 is a refinement of the approximation K- ’ % J. This inequality K-l 5 J is also useful. It is better than K- ’ I 0 and can be used to get information about K and its inverse. For instance, if the eigenvalues {ni} of K are ordered such that
then the above inequality can be used to show that k,, 2 IA,1 : [A.[-’
= p(-K-l)
2 p(-J)
= k,,
[6].
Moreover, the inequality is sharp in that there are SEC matrices for which K- 1 = J( II, I). For example [9] and [lo] use a two compartment single exit model in which
K = [ -:z:
_::49]
If we choose II, I = 0.716, K-i
= [ 1;:;;
_;:;,I
= J@,I,.
The inequality K - ’ I 0 can also be improved to K - 1 < 0 in certain cases. If the SEC matrix K is irreducible, then K- ’ < 0 [3]. Another set of sufficient conditions is contained in the following result. Theorem 6. l’fthe mth column of the SEC matrix K consists entirely of nonzero elements, then
K-l<0 As an example, suppose K is the SEC matrix associated with a mammillary system with excretion from the mth or central compartment [2]. If we assume that k,j # 0 for allj, then by Theorem 5, K is invertible. Moreover, if as is common, k, # 0 for all i, then K- l < 0 via the last Theorem. 5. AN UPDATING
PROCEDURE
If n is large, there can be a very practical inefficiency with the iterative scheme (3) for computing x. The iteration requires us to retain all of the components of xtp’ until the computation of xtp+ ” is complete. A more natural idea which cuts storage is to update the iterates, that is, if the computations are done sequentially, to start using each component of the new xtp+ ” as soon as it is calculated. Thus xtp+ ” replaces xcp’a component at a time and
DAVID H. ANDERSON
324
so x(p)can be dumped as x (p+i’ is created. For actual computation, this modified process can be expressed as i-l Xi
(P+l)
=
1
f?ijX!f
+ 1’ +
i
e&P’
+
xp,
(4)
j=i
j=l
i= 12 , ,a.., n. Since this scheme makes use of more recent information than does the previous
scheme, it is natural to think this new method will converge and its asymptotic rate of convergence will probably be faster. However, this is not so obvious if we observe that the relationship between our original and updated schemes is analogous to the relationship between the Jacobi and the Gauss-Siedel methods for the solution of a linear algebraic system and that it is possible for either of these methods to converge while the other diverges. Also, for all systems of dimension greater than or equal to 3, when both of these methods converge, the rate of convergence of the Jacobi iteration may be less than, equal to, or greater than that of the Gauss-Seidel iteration [ 111. To analyze our updating procedure (4), we go to a matrix format. Let L be the strictly lower triangular part of E and U be the upper triangular part of E including its diagonal. Then (4) is equivalent to xtP+ll = Lx@+11+ Ux(P) + x(oJ. This equality may be rewritten to give the matrix-vector
form of the new algorithm as
x’p+l’ = (I - L)_‘Ux(P’ + (I - f,-‘Jy.
(5)
The matrix I - L is invertible since det (I - L) = 1. Theorem 7. The iterative solution x of Kx = y
scheme
(5) always
produces
a sequence
which converges
to the
The updating procedure does in fact accelerate the convergence of the solution process without complicating the iterative method. This is the result of the next theorem. Theorem 8. The iterative scheme (4) hus at least as fast a rate of convergence 6.
ERROR
us (3)
ANALYSIS
In what has already been described, the procedure x(P+
11 =
QX(P’
+
d,
d = (I - L)-‘Jy,
is given for approximating, with any preassigned degree of accuracy, the solution x of Kx = Y. If a matrix norm 11.(1is found (one always exists) for which 11Q II < 1, then a standard error estimate can be made if desired 1123. We start by noting the representations x(P)
=
i i=O
x = t
Q’4 Q’d.
i=O
If II.II is a vector norm which is compatible with the above matrix norm, then 11 x(P) - x/I = 11i-t+1
Q'dIi
=
IIQ""f. Q'dIl 5 IIQll4"Ildll
=
lIQ(IP+llldll(l
- lIQII)-‘3
izo ItQIl'
Iterative inversion of single compartment matrices
325
where the last expression can be computed from knowledge of 1)Q I]and 1)d )I. In particular, to find an index p such that IIx(P) - xl] < ( for prescribed c > 0, the inequality
IIQIIP+llldtIU - IIQIO-' -=E can be solved for the least such integer p. It should be noted that when performing practical computations on most digital computers the inherent peculiarities of the machine’s floating point arithmetic may require some alteration of the above-mentioned p value.
7. CONDITIONING
OF A SEC MATRIX
In this section we examine the vulnerability or sensitivity of a solution: If K is slightly perturbed, how great is the effect on x = K-‘y? A standard measure of sensitivity is the condition number for a matrix A, cond (A) = 1)A II II A-’ )I 2 1, using an appropriate matrix norm [13, p. 2711. If cond (A) zzz 1, the matrix A is wellconditioned. The larger cond (A) is, the more ill-conditioned is the system since the solution is more vulnerable to round-off errors. As we will see from the calculations below, a convenient matrix norm to use is the infinity norm. This matrix norm for an n x n real matrix A 3 [I,] can be calculated by
11 A 11 max i laijl[8]I
=
I
j=l
Since each element in the mth row of K- l is -k,-,‘, then
where the notation K- ’ = [ki; ‘I] h as b een used. Thus the condition number
Another lower bound on this condition number can be used when the maximal eigenvalue of K and the excretion rate constant k,, are known (as in the inverse problem). Retain the ordering of the eigenvalues {&} of K as 11112 IArl r ... 2 ]&I, and use the fact that 11K (Im 2 I& I. Then cond,(K)
8. A SIMPLE
2 II, Ink&,‘.
EXAMPLE
In our previous work the algorithm using updated components for finding x = K- ‘y has been developed. An illustration of the methodology is now presented. The example is to solve the system Kx = y where Y = (Yl,Y,,Y#
= (-LO,W
and ‘-0.805 K=
0.7 1
0.1
0.2 -0.7 0.5
0.4 0.7 - 1.1
1,
J
DAVID H. ANDERSON
326
The matrix K is seen to be invertible by an appeal to Theorem problem is x = (200, 400, 2OO)r. Note that k,, = 0.005,
5. The exact solution
of this
(1KIIw = 2.1,
so that cond,(K) 2 1260. Using the method developed in this paper, we get an immediate advantage over other methods since we have a closed form solution for the mth component in the solution x z (x,, x2, X3)=: 3 xl
Now the J matrix
=
-
C i=l
y,/k,,
= 200.
is (with c E 1.1)
[;
so that the starting
---l;l
_-;I,
vector x(O) E (x(p), x$0’, x:“‘)’ = Jy
is (2OO,O,O)r. The E matrix
e13 -0 -
e12 --0
-
ell --0
has entries
ezl = k2Jc = 7/11
ez2 = (c + k,,)/c
e31 = k3Jc = l/11
eJ2 = k32/c = 5111
ez3 = k,,Jc = 7111
= 4/11
e33 = (c + k33)/c = 0
(note that if c is defined to be max [kit] z k,,, then always eq4 = 0, which simplifies algorithm even further). Thus the computational scheme for this problem becomes
the
x’p+ l’ = (4x’p’ + 7x’9’ + 1400)/l 1 x(g+ 1’ = (Sxcp’ I’ + 200)/l 1 where the first equation in (4) has been disregarded since we know the exact component x1 = 200. These expressions can be easily programmed and the resulting sequences generated by them converge rapidly (- 30 iterations) and monotonically from the initial vector up to X) = 399.99999, using the stopping
criterion I(x(p+l’-
9.
xj = 199.99999,
COMPARISON
WITH
x(p’I~l/I]x(p+l’(J1 DIRECT
< 10-S.
NUMERICAL
METHODS
Thirty-two test systems Kx = y in the range 3 I n I 8 (the usual size of the problems we have encountered) and 1 I cond,(K) 5 lo6 with known exact solutions were devised. A subroutine LEQSEC (Linear EQuation solver for SEC systems), utilizing the above updating iterative scheme, was written specifically for linear systems involving SEC matrices and has been extensively tested and compared to the direct methods of Gaussian elimination, Crout reduction with iterative improvement, and singular value decomposition [ 141. Gaussian elimination failed to give “exact solutions” (answers satisfactorily close to the known exact solutions) in 4 of the 32 cases. Crout reduction failed to produce exact solutions for 6 cases; of these Crout reduction with iterative improvement was able to compute exact solutions for only 2 cases. In all of the trials that failed, the matrix was classified by the program as “algorithmically singular” because too small a pivot element was encountered. When the test matrices were ill-conditioned to this extent, the standard stabilization strategy of singular value decomposition was used. The LEQSEC algorithm is also numerically stable
Iterative
inversion
of single compartment
matrices
327
since p(Q) < 1 [ 121. On none of the remaining 4 test cases did singular value decomposition give more accurate answers than did LEQSEC, and in some cases the solution was much less accurate. In fact, no SEC system has been encountered thus far on which the LEQSEC algorithm failed to yield the exact solution and so LEQSEC would seem to be a preferable method for solving SEC systems on machines which have a comparatively small range of representable numbers or machines which are known to be numerically “dirty”. In this testing process, all computations were done on a CDC CYBER 70, model 72, or a DEC system- 10. The respective subroutines used were DECOMP-SOLVE [ 151, LEQTlF and LEQT2F in the International Mathematical and Statistical Libraries package, and SVD [ 163. The experimental comparisons did shed some light on the one shortcoming of LEQSEC, which was the relatively larger time required on a few problems to obtain a solution as compared to the direct methods. However, since these problems were exactly those most illconditioned, slow convergence of the iterative scheme was to be expected. Work is now in progress to modify the algorithm so as to further increase its rate of convergence. SUMMARY In the study of biological systems a frequently used model is an open compartmental system having only one compartment from which elimination takes place. Let K be the compartmental matrix associated with the mathematical model of such a system, i = Kz + b. This article treats the theory of solution of the simultaneous set of linear equations Kx = y which arises in the analysis of the model. Simple sufficient conditions for the invertibility of K are found, properties of the inverse matrix are discovered, and an iterative scheme is developed for which convergence to x is guaranteed. The algorithm’s progress is accelerated and its storage requirements are made more modest by an updating procedure. An example comparing this computational technique to other methods is given, and a general error analysis is discussed. Acknowledgements - In the formulation of this report, Drs. Myron Ginsberg and G. Milton Wing made many valuable comments. I would also like to express my appreciation to Sonja Sandberg who did substantial work in the programming and computing necessary for the research.
REFERENCES 1. P. B. Bright, The volumes of some compartment systems with sampling and loss from one compartment, Bull. Math. Biol. 35, 69-79 (1973). Anulysis in Biology und Medicine. Elsevier, Amsterdam (1972). 2. J. A. Jacquez, Compartmentul 3. J. 2. Hearon, Theorems on linear systems, N. Y. Acud. oJSci. APIIIA/S108, 36-68 (1963). 4. D. Fife, Which linear compartmental systems contain traps? Muth. Biosci. 14, 311-315 (1972). in single-exit compartmental systems, Applied Nonlineur Anulysis, 5. D. H. Anderson, The volume of distribution V. LAKSHMKANI-HAM Ed. Academic Press, New York (1979). 6. J. M. Ortega, Numerical Analysis. A Second Course. Academic Press, New York (1972). I. L. B. Rail, Computational Solution of Nonlineur Operator Equutions. John Wiley, New York (1969). Theory of Matrices. Academic Press, New York (1969). 8. P. Lancaster, 9. B. K. Shah, Data analysis problems in the area of pharmacokinetics research. Biomen+cs 32, 145% I57 (1976). C/in. Pharmucol. Therupeur. 8, 201-218 (1967). IO. J. G. Wagner, Use of computers in pharmacokinetics. 11. S. Venit, The convergence-of Jacobiand Gauss-Seidel iteration. Muth. Maya& 48, 163-167 (1975). Theory und Practice. Addison-Wesley, Reading, Ma12. E. K. Blum, Numerical Analysis und Computation. ssachusetts (1972). Academic Press, New York (1976). 13. G. Strang, Linear Algebra and ifs Applicufions. 14. C. H. Farnsworth, Comparing numerical methods for solving linear systems of equations involving single-exit compartmental systems, Unpublished report (1978). 15. G. E. Forsythe, M. A. Malcolm and C. B. Moler, Computer merhods for muthematicul computurions. Prentice-Hall, Englewood Clifis, N.J. (1977). 16. G. H. Golub and C. Reinsch, Singular value decomposition and linear systems solutions. Num. Muth. 14, 403-420 (1970). 17. G. Forsythe and C. B. Moler, Computer Solution ofLinear Algebraic Systems. Prentice-Hall, Englewood Cliffs, N.J. (1967). About the Author - DAVID H. ANDERSON received his B.S. degree in Mathematics from the University of Mississippi in 1961. From 1961 to 1965 he was an instructor of mathematics at the U.S. Naval Nuclear Power School, Vallejo, California, involved in training the crews of the nuclear ships of the
DAVID H. ANDERSON
328
U.S. Navy. After working one summer for Exxon in Baton Rouge, Louisiana, he took the M.A. degree in Mathematics from the University of California, Berkeley, followed by a year as Assistant Professor of Mathematics at Millsaps College, Jackson, Mississippi. Earning his Ph.D. in Mathematics from Duke University in 1969, he joined the faculty at SMU in Dallas. At present, Dr. Anderson is an Associate Professor in the Department of Mathematics at SMU and in the Department of Medical Computer Science at the University of Texas Health Science Center at Dallas, and is a Visiting Lecturer for the Mathematical Association of America. During the period January-May 1979, he had a visiting position with the Division of Applied Mathematics, Brown University, Providence, R.I. He has published several articles dealing with mathematical modeling, numerical analysis, and algorithm design in compartmental analysis and population mathematics.
APPENDIX Proof of Theorem 2 Let u be the I x n row vector of zeros except for a I in the mth position. each entry unity. Then the condition
Take w to be the I x n row vector with
leads to the equation WK = - k,,u. Since K is nonsingular,
then UK-’ = (-l/k&v.
But UK-’ is just the mth row (vector) of K-l. - l/k,.
(6) (6) says that each entry in the mth row of K-’
Hence equation
is
Proof of Theorem 3 The (i,j)-element
of E is 0 eij
-i
i=m
if
1 + k,Jc
if
i # m,
i=j.
k,k
if
if
i#j
m,
Since each entry in the mth row of E is zero, then E is a singular reducible matrix. Also since c 2 1kiilfor all i, and k, r 0 whenever
i #j, then each e,, 2 0. Thus E is a nonnegative
matrix.
Now the jth column
sum of E (j # m) is
= I - k,+, which is less than or equal to unity. The mth column
sum of E is
=(-km, + ILmlVc, a value strictly
less than one since k,
z 0. These sums tell us I
(in particuiar, nonempty.
observe
j#H;
C e,j < l i=,
if
,$, eu=
if jeH.
that m $ H.) This is equivalent
I
to saying
(/E
//, <: 1 ifHisempty,and
l/E/,
=
I ifHis
Proof of Theorem 4 Let K be invertible. We wish to prove that p(E) < 1. The magnitude of any matrix norm of E. Thus, by the last theorem,
of any eigenvalue
of E cannot exceed the value
p(E) 5 IIgIIt 5 1 By way of contradiction, let us assume that p(E) = I. Appealing to the PerronFrobenius theory for nonnegative matrices [8]. we can say that there is a nonnegative eigenvalue r of the matrix E such that r = p(E) = 1. Since rl - E = JK, then 0 = det(rf
- E) = det(JK)
= det J
det K.
Iterative
inversion
of single compartment
329
matrices
Since J is nonsingular, we must conclude that det K = 0. But this contradicts the invertibility of K ; hence p(E) < 1. On the other hand, suppose p(E) < 1. Given any 8 > 0, there is a matrix norm 11 t 11such that
11 El\ I p(E) + 463. Let I: = (1 - p(E))/2. Then III - JK )I s
IIE 11< 1 since p(E) < 1. Thus J is an approximate
inverse of K and so K is invertible.
Proojof nteorem 5 From the proof of Theorem
3, we know that C elj + 1
d H. Since this is a strict inequality
wheneverj
1
for every j, there exists an E > 1 which will serve for all j $ H such that
2 ei,+tz matrix
1 et,< 1. id/
r$fl
Let D s Cd,,] be the n x n diagonal
eij <
idl
i$ll
with d,, =
i
if jeH
;
if j$H
Define the matrix C = [c,,] 2 0 through the similarity transformation C s DED- ‘. Thejth row of C,j E H, is the jth row of E multiplied by ~zand the jth column of C, j E H, is the jth column of E multiplied by l/c. Otherwise the entries of C are the same as those of E. Form the Gerschgorin disks [S] based on the columns of C. This will give information about the eigenvalues of E since C is similar to E. For each j = 1,2,. , n, we have the disk of complex numbers 0,s
z:]z - c,,] I {
For the nonnegative
maximal
can be written r I
I
as
1 eilo + 1 eij& id”
k”
of the similarity transformation. Consider the last sum just given. Since in the summation only when i # j,. Recall that when i # j, ei, =
By hypothesis, Thus
“i,
r = p(E) of E, there exists at least one index j, for which r E bj0. Thus
eigenvalue
Suppose j0 E H. Then the last inequality
Via the construction quantity eilo occurs
1 i#i
0
if
i=m
k,/c
if
i#m
there is at least one i, 4 H for which k ,0,0 # 0 (sincej,, E H). This index i,, cannot
i $ H and& E ff, the
be m because k,,,,” = 0.
and thus the sum
is nonzero.
obtains
As a result, the strict inequality
and thus
r<
1 $I/
eij, +
1 eijo= j,
i+tl
eih =
1.
Now suppose j0 $ H. Then we can write
But since j, 4 H, the right hand side of this inequality must be less than unity. Thus for either possibility, j, d H or j, E H, we have p(E) = r < 1. Hence by Theorem proof is complete.
4, K is invertible
and the
DAVID H. ANDERSON
330 Prwf
of Theorem 6
The matrix
J has negative
entries everywhere
on its mth row and main diagonal.
Thus the matrix
x,=.J+JV, which is less than or equal to J, has negative elements in the same positions. Now concentrate on row numbers i # m and column numbers j # i. Introduce the matrix notation J = [Ji,]. For these i and j values, the (i, j)-element of EJ is i
p=*
ei,,Jpj = e,,,,J,,,, + e,lJ,, = - k,&k,,)-’
- kifi- L
ifj # tn. and is
%Jnn = - h,kkm) - 1 ifj = tn. In either case, the (i,j)-entry
of EJ is negative since, by hypothesis,
ki, > 0 for all
i = I, 2,.
, n, i # m. Thus
K-‘
(5) is of the general
linear iterative
form
xc&‘+11=
QX’P’
+
d
for the square matrix Q = (I - L)- ‘U. In order that the sequence of iterates converge, that p(Q) < 1 [6]. We prove that p(Q) < 1 by showing that the series
it is necessary
and sufficient
I + Q + Q’ + Q3 + . . converges. Now
Q = (I - L)-‘u
= (I + L + L* + . . . + L”-‘)U
is a non-negative matrix [6]. The series expansion triangular. Define A= Then E = I - JK imples B-’
-J
-K,
= (I - E)A-’
= (I - L)-‘(-J)
for (I Es
L)-’
terminates
with L’-’
since L is strictly
lower
-J-‘(I-L).
and so
= (I - L)-‘(I
- E)A-’
= (I - L)-‘((I
- L) - U)A-’
=A-‘-(I-L)-‘(IA-‘=(I-Q)A-‘. Since (I-L)-‘=I+L+L2+...+L’-t20 and -J
2 0, then 8-l
2 0. We also observe
that A-’
in [6, p. 1191.
= - K - ’ h 0. From here on, the proof is similar to that
Proof of Theorem 8. For the general
linear iterative
process
r
x@‘+ ’ G Mx@ + d,
p =
0, 1,2,.
,
(7)
the asymptotic convergence factor is p(M), provided p(M) < 1[6]. For (4), the matrix A4 z Q and p(Q) < 1 by Theorem 7. In the case of (3). M = E and via Theorem 4, we have p(E) 4 1. Thus the asymptotic convergence factors of the 2 iterative schemes are p(Q) and p(E), respectively. The smaller p(M) is, the faster the rate ofconvergence of the process (7) [6]. Clearly then, we need to prove that p(Q) 5 p(E). Let n = p(Q). Since Q is a nonnegative matrix, the Perron-Frobenius theory guarantees the existence of a nonzero vector IAsuch that (I - L)-
’ Uu z Qu = qu [S-J.
Then (?$ + U)u = vu, which says that q is a nonnegative
eigenvalue
of the matrix
qL + U. Thus
P(rlL + U) 2 n. SinceLzO,U20,andO~rr<
l,then O
Hence
the desired
inequality.
Compur. Biol. Med. Vol. 9, pp. 317-330. @ Pergmon Press Ltd. 1979. Printed in Great Britain.
ITERATIVE INVERSION COMPARTMENTAL DAVID
OF SINGLE EXIT MATRICES
H. ANDERSON
Department of Mathematics, Southern Methodist University, Dallas, Texas 75275 and Department of Medical Computer Science, The University of Texas Health Science Center, Dallas, Texas 75235, U.S.A. (Received
25
July 1978;
in
revised form 5 February 1979)
Abstract - This article develops an iterative method for solving a system of linear algebraic equations in which the coefficient matrix is a compartmental matrix of the single-exit type. Matrix properties, invertibility, ill-conditioning, algorithm acceleration and error analysis are discussed. Matrix Single-exit compartmental systems Simultaneous linear algebraic equations Numerical iterative methods Error analysis and invertibility Acceleration of convergence rate of convergence
1. INTRODUCTION In the mathematical modeling and simulation of phenomena in such areas as biology, medicine and chemistry, there naturally arises an important class of square matrices in which each matrix has the following properties: (a) each off-diagonal element of the matrix is nonnegative; (b) each diagonal element is negative ; (c) all columns except one add to zero ; (d) that exceptional column sums to a negative number. A simple example is -0.3 0.1
K=
0.1 -0.1
0.2
0.0
0.0
0.0
0.0
0.2
0.2
0.2
-0.4 0.0
0.1 -0.5
J
One way such matrices come about in practice is through linear compartmental models associated with radioactive tracer studies dealing with human physiology [ 11. For example, suppose there are 3 identifiable compartments in a certain clinical study and of these there is only a single compartment, the first, through which exit to the outside -environment is possible (this “single exit” is the essential characteristic of the system for our study). Let x,(r) be the amount of tracer in the ith compartment. Assuming the kinetics of the system are basically linear, the mathematical model for the system dynamics is taken to be (a s dx/dt) ~2, = - k,,x,
- kS1xl - k,,x,
+ k12x2 + k,,x,
Jz =k2,x1
+ kz3x3 -k,zxz
-kszxz
z& = k,,x,
+ k3,x,
- kzJxg + b,
- k,,x,
+ b,
+b,
in which k, is the nonnegative fractional transfer coefficient from compartment j to compartment i, k,, is the positive excretion rate constant, and b,(t) is some input function [2, p. 481. This system can be written in the vector-matrix form i=Kx+b 317
DAVID H. ANDERSON
318
where
-621 + k,, + kc,,) K=
k 21
k
k,2 -&2
k 31
+
k::
k32)
k 32
-h
+
k23)
This matrix has properties (a)-(d) listed above. 2. PRELIMINARIES Following the general development adopted in [2, p. 481 or [3], let K z [k,] be an n x n real matrix in which k, 2 0 if i # j, and for all j = 1,2,. . . , n, l#j
where k, are nonnegative numbers and the summation is over all i = 1,2, . . . , n, except i = j. To be nontrivial, we assume that k,, < 0 for all j. We now delineate the important subclass of interest to our study. Suppose that in addition to the above properties, the quantities k, are taken to be zero for all j = 1,2,. . . , n with the single exception ofj = m, for which we have k,, > 0. Motivated by the introductory example, such a matrix K possessing these properties will be called a single-exit compartmental (SEC) matrix.
There are at least 2 interesting problems concerning such matrices. First, that of determining any unknown entries of the matrix from observation of certain compartments of the system in time (the “inverse” or “system identification” problem), and second, that of solving a linear algebraic system having a SEC matrix as the known coefficient matrix, or equivalently, inverting a SEC matrix (e.g. equilibrium or steady-state problems associated with single exit compartmental systems give rise to these types of algebraic equations). The present paper is concerned with the second problem, that of obtaining a general theory and method of solution of the inversion problem. Given any n x n SEC matrix K and n x 1 real vector y, our goal then is to investigate the theory of solution of the linear algebraic system Kx = y. The method of approach is to construct an approximate inverse of K, from this derive conditions under which K is invertible and then design an iterative scheme that will always take us from a good initial guess to the solution x. Knowledge of the proposed computational scheme is appropriate for a variety of reasons. The matrix K is often ill-conditioned and this complication, which causes difficulties in direct methods such as Gaussian elimination, may be overcome by using the iterative technique to be developed. Ill-conditioning of K is definitely indicated since the condition number 11K /I 11 K- ’ 1)may be quite large. Also K- ’ may have large values relative to K, and a small excretion rate k,, is not atypical in biomedical problems (det K --+0 as k,, + 0). Through the iterative approach, upper bounds on the error involved are possible. An important consequence of the approach is the resulting mathematical theory obtained about the structure of K and its inverse matrix. The proposed iterative method may be applied in order to improve an approximate solution already obtained by some other method. Many present modeling efforts yield sparse matrices, and for such matrices the proposed computational scheme is simple and convenient for use, and its storage requirements economical. Probably the only disadvantage of the successive approximation method is that convergence may be slow. However, we are able to develop a method of accelerating convergence which does not complicate the iterative formulas and in fact decreases the scheme’s storage requirements.
3. INVERTIBILITY
OF A SEC
MATRIX
In treating the theory of solution of Kx = y, we should first face the question of invertibility of the matrix and try to ascertain simple conditions governing nonsingularity of
Iterative inversion of single compartment matrices K.
319
There do exist singular SEC matrices ; an immediate example is -3
0 1
0
I
2
-4
1 4 -2 K= [ Let us start with a known invertibility result for genercrl compartmental matrices (matrices K in which k, may be zero or nonzero). Draw the connectivity dioyrum D for the specified compartmental matrix K = [kli] (kij > 0 for i # j means to draw the directed line segmentJ) and let the associated compartmental system be denoted by s E {1,2,...,n} where the integers denote the numbering of the compartments (or vertices) in D. Suppose (renumbering the compartments if necessary) we define the subset T={p,p+l,...,n)(p
of S. This subsystem T of the compartmental system S is called a rrap provided k, = 0 for all i,j such that j 2 p and i < p. In other words, there does not exist a path in D from any compartment in T to a compartment which is in S but outside T. See [4] and [5] for information on the following result. Theorem 1. The compurtmentul mutrix compurtmentul system S has no trups
K is invertible if and only if its ussociuted
Consider the following (SEC) matrix K and its connectivity diagram D (Fig. 1):
K=
1
0
1
1
0
-
0
Exit .
. (1)
0 -1
I
Fig. 1.
Observe that this matrix is reducible so that the standard inversion theory for irreducible diagonally dominant matrices [6, p. 1063 does not apply. However, as can immediately be seen from D, the matrix is nonsingular since the associated compartmental system contains no traps. For the particular subclass of SEC matrices, an alternative necessary and sufficient condition for invertibility can be given (Theorem 4) which is vitally important for the iterative process to be developed in this paper (Section 4). We make use of the result that a square matrix is invertible if and only if there exits an “approximate inverse” of A [7]. To be specific, if for a given n x n real matrix A there is an invertible n x n real matrix B and a matrix norm I] . 11such that (II - BA II < 1, then B is said to be an approximate inverse of A (I is the n x n identity matrix). Note that if B = A- ‘, then I - BA is the zero matrix so that we may consider the size of 1 - BA to be a rough estimate of the inaccuracy of B as an approximation to A-‘. The underlying theme of this paper centers on the construction, for a specified n x n SEC matrix K, of an invertible n x n real matrix J such that (II - JK II1 I 1. Here the l-norm or column norm [8, p. 2123 i IIAII1= max I
(ulil
DAVID H. ANDERSON
320
is a natural choice since much of the structure of the SEC matrix is geared to its column sums. Let a SEC matrix K be given. Since J should look something like K- ’ (if K- ’ exists), then it is desirable to know certain fundamental information about this inverse in order to construct J. After calculating a few such inverses, the following theorem is not difficult to conjecture. All theorem proofs appear in the appendix. Theorem
2. If the SEC matrix K is invertible, then each entry in the mth row of K- ’ is - 1 /k,,,
Judging from the literature, some of the most interesting and important cases of SEC systems are the mammillary or cutenary types [2, p. 531. Theorem 2 allows us to get a direct solution of Kx = y in some of these cases. The crucial observation is that since each entry in the mth row of K- ’ is known to be - l/k,, then x = K- ‘y yields the closed form calculation n
1 Yilkww
Xm = -
i=l
If, in the case of a mammillary system, the number of the central compartment is designated one, then the coefficient matrix K assumes the form where there are nonzero entries only in its first row, first column, and main diagonal. The mth component x, is now known from the above equation and thus x, can be computed from the mth equation of Kx = y. Then for i = 2,3,. * . , m - 1, m + 1,. . . , n, we can read xi = (yi - kiixi )lkii directly from Kx = y. If K represents a catenary system with the usual sequential numbering ofcompartments, then it is a tridiagonal matrix possessing nonzero elements only on its main diagonal, first superdiagonal, and first subdiagonal. Suppose that excretion takes place just from the terminal or nth compartment (m = n). Then x, is known from the closed form (2) and the system Kx = y can be solved by back substitution. For, the component x, _ 1 is found from the equation k,.,- ix.- 1 + k,,x, = Y,, and then the equations ki.i-ixi-i i = n - 1, n - 2,...,
+ kiixi + ki,i+rxi+r
2, can be successively
= Yi3
solved for x._~, x,_ 3,. . . , x1.
In terms ofdesigning J, the conclusion coming out of Theorem 2 is that a reasonable choice for each element in the mth row of J is - l/k,,. There is a second observation that aids in the formation of J, and that is the diagonal dominance of the K matrix suggests that generally the off-diagonal elements of K will be small in comparison to the diagonal elements. Thus in a loose sense, K z ditry{k,},
and so the diagonal
entries
of K- ’ might be close to l/k, ,r
V&,.
..,
ilk,,.
Unfortunately, if the diagonal elements of J are chosen as {l/k,,}, we do not always get )II 1 for any SEC matrix K selected. Instead of taking the option ofchoosing a different norm, we make a slight modification in the choice of the diagonal elements of J which will give us 111 - JKII 1 I 1 for every K. Define the parameter c to be any number satisfying JK 11l I
c 2 max Ikii(. I
Then let every diagonal position of J be filled by - l/c. This discussion motivates the actual definition of the J matrix. For a specified n x n SEC matrix K, let J [or, if we wish, J(c)] be the n x n matrix with - l/k,, as each element in its mth row, the value - l/c as each element on its main diagonal (in other than the (m, m)-position), and zeros elsewhere. For example, consider the 4 x 4 matrix displayed in (1). It is such that
Iterative inversion of single compartment
m =
matrices
321
1, k,, = 1 and c 2 2. If we take c = 2, then the J associated with (1) is -1
-1
-1
0-i J = J(2)=
0
-1 0
0-t
0 0
0 0 0 -9 [ I Let us observe a few things about the .Z matrix. It is always invertible since detJ = (-l)“/k,,c”-’
#O.
It is easy to verify that J-l has -k,, as its (m,m)entry, every other element on the main diagonal is -c, each nondiagonal entry of the mth row is c and zeros elsewhere. In fact, .Z- ’ is itself an SEC matrix. For any SEC matrix K form the associated J matrix according to the above rules. Then define the n x n matrix E = E(c) E I - J(c)K.
We can now prove that 11Z - .ZKII1 I especially the nonnegativity of E. Theorem
1
as well as get some interesting ancillary results,
3. Define the subset of indices H E {j: kmj = O}.
Then )IE (11 < 1 if H is the empty set, and 1)E 111 = 1 if H 1s . nonempty. Moreover, E is a singular reducible non-negative
matrix.
We conclude via Theorem 3 that J is a reasonable approximation of K-’ and thus in a rough sense K z 1-l. This observation provides an interesting “first order” interpretation of the interworkings of a single exit compartmental system. To within this approximation, the only non-zero fractional transfer coefficients are those coming into the mth compartment, k,,
= c,
j = 1,2,. . . , n, j # m,
and thus the system is acting like a mammillary system with excretion from the mth or central compartment [2]. This theorem allows us to derive a number of results. The first of these ‘is the basic invertibility condition, which is stated in terms of the spectrul radius of a square matrix A, p(A) = max{ Irl : r is an eigenvulue of A}. Theorem
4. The SEC matrix
K is invertible if and only if p(E) c 1.
The above theorem tells us that J & an approximate inverse of K whenever K is nonsingular. The result also reduces the invertibility of K to the algebraic problem of demonstrating that every eigenvalue df E has magnitude less than unity. Even though there are presently available excellent methods of computing eigenvalues, the equivalence of Theorem 4 is theoretical in nature and not readily applicable to a particular matrix E. However, based on this fundamental equivalence we can derive simple sufficient conditions on the entries of K to determine its invertibility. The next criterion is especially satisfying from the modeling point of view because it depends only on the disposition of zero and nonzero elements in the matrix K and thus relates to Theorem 1. Upon commencing the modeling process, this disposition is often all that is known or even conjectured about kii. Theorem 5. Recall that H = {j : kmj = 0). Zffor each j E H there is at leust one i 4 H such that kij is nonzero, then K is invertible
Applying Theorem 5 is very easy. Simply strike out all rows with row numbers in H and all columns with column numbers not in H. If in each of the remaining columns there exist at least one nonzero element not eliminated, then K is invertible.
DAVID H. ANDERSON
322
4. THE
FUNDAMENTAL
ITERATIVE
SCHEME
The problem to be considered in this section is that of finding the solution vector x to satisfy the linear system Kx = y where the SEC matrix K will be assumed nonsingular. We will develop an iterative scheme based on our knowledge of the approximate inverse J of K, a theory which has already been demonstrated to be rich in ideas. A method of successive approximation will be set up as follows. The equation Kx = y implies that x = (I - JK)x + Jy = Ex + Jy, a fixed point equation for the true solution x. Thus the computational X(P+“=EX(P’+Jy,
p=o,1,2
scheme
,...)
is suggested. Moreover, because E is small in the sense that p(E) < 1, the successive approximation is guaranteed to provide a sequence that converges regardless of the starting vector [6]. Since x = K-‘y z Jy, we choose to initiate the process at x(“’ E Jy. If the iterates are successively evaluated, then
xcl' =
Jy + EJy
xc” = Jy + EJy + E2Jy
is a sequence converging to EPJy
x= f
(E” = I).
p=o
When the scheme is written in component computation on most machines, we have x(P+ll I
=
jgl
form, which should be used for actual
eijXy'
+
X$O’,
i= 1,2 ,..., n. Since eij = 0 when i = m, we immediately get the mth component solution x in closed form:
(3) of the
n
xm = x’,“’= - i;l Yilkmw Thus the scheme (3) can be reduced in dimension by one by dropping the mth iterative equation and replacing all x2’ by x,. Let us at this point introduce the component-wise partial ordering I (or 2) for real vectors x = (x,,. . . ,x.) and y = (vi,... . ,y,) by x I-y
if and only if x, I yi,
i = 1,. . . ,n,
and similarly for n x n real matrices A = [Uij] and B = [bIj], AI B
if and only if Ui, I bu
i,j = 1,. . . ,n.
Now suppose y I 0, a common situation in biomedical work. Since J I 0 and E 2 0, then Jy 2 0 and EPJy 2 0 for all p = 1,2,... . Hence {x@‘}forms a monotonically increasing sequence that converges unilaterally to the solution x of Kx = y. If there are many y vectors to be used, there may be an advantage in finding K- ’ and then proceeding by the simple multiplication K- ‘y in each case of y to the solution x. One way of generating K- ’ is through an algorithm parallel to that above. The relationship E = I - JK yields the fixed point equation K-‘=EK-‘+J
Iterative inversion of single compartment
for the inverse matrix; thus the computational X p+i =EX,+J,
matrices
323
scheme p=o,1,2
)..‘)
with X, = J is motivated. Because this is a linear iterative process and p(E) < 1, the sequence X,=J 2,
=J+EJ
X2 = J + EJ + E’J
is guaranteed to converge to K-‘.
Moreover, since EPJ 5 0 for all p = 1,2,. . . , we have
This is an improvement of the known result that K-l < 0 [3]. These inequalities on the sequence elements are very useful not only because they yeild such information as x > 0 (in the case when y I 0) and K- ’ I 0, but also because if either iteration is terminated at any step p, then it is known in advance that x(P)is a lower bound for x and X, is an upper bound for K-l. The particular upper bound J of K - 1 is a refinement of the approximation K- ’ % J. This inequality K-l 5 J is also useful. It is better than K- ’ I 0 and can be used to get information about K and its inverse. For instance, if the eigenvalues {ni} of K are ordered such that
then the above inequality can be used to show that k,, 2 IA,1 : [A.[-’
= p(-K-l)
2 p(-J)
= k,,
[6].
Moreover, the inequality is sharp in that there are SEC matrices for which K- 1 = J( II, I). For example [9] and [lo] use a two compartment single exit model in which
K = [ -:z:
_::49]
If we choose II, I = 0.716, K-i
= [ 1;:;;
_;:;,I
= J@,I,.
The inequality K - ’ I 0 can also be improved to K - 1 < 0 in certain cases. If the SEC matrix K is irreducible, then K- ’ < 0 [3]. Another set of sufficient conditions is contained in the following result. Theorem 6. l’fthe mth column of the SEC matrix K consists entirely of nonzero elements, then
K-l<0 As an example, suppose K is the SEC matrix associated with a mammillary system with excretion from the mth or central compartment [2]. If we assume that k,j # 0 for allj, then by Theorem 5, K is invertible. Moreover, if as is common, k, # 0 for all i, then K- l < 0 via the last Theorem. 5. AN UPDATING
PROCEDURE
If n is large, there can be a very practical inefficiency with the iterative scheme (3) for computing x. The iteration requires us to retain all of the components of xtp’ until the computation of xtp+ ” is complete. A more natural idea which cuts storage is to update the iterates, that is, if the computations are done sequentially, to start using each component of the new xtp+ ” as soon as it is calculated. Thus xtp+ ” replaces xcp’a component at a time and
DAVID H. ANDERSON
324
so x(p)can be dumped as x (p+i’ is created. For actual computation, this modified process can be expressed as i-l Xi
(P+l)
=
1
f?ijX!f
+ 1’ +
i
e&P’
+
xp,
(4)
j=i
j=l
i= 12 , ,a.., n. Since this scheme makes use of more recent information than does the previous
scheme, it is natural to think this new method will converge and its asymptotic rate of convergence will probably be faster. However, this is not so obvious if we observe that the relationship between our original and updated schemes is analogous to the relationship between the Jacobi and the Gauss-Siedel methods for the solution of a linear algebraic system and that it is possible for either of these methods to converge while the other diverges. Also, for all systems of dimension greater than or equal to 3, when both of these methods converge, the rate of convergence of the Jacobi iteration may be less than, equal to, or greater than that of the Gauss-Seidel iteration [ 111. To analyze our updating procedure (4), we go to a matrix format. Let L be the strictly lower triangular part of E and U be the upper triangular part of E including its diagonal. Then (4) is equivalent to xtP+ll = Lx@+11+ Ux(P) + x(oJ. This equality may be rewritten to give the matrix-vector
form of the new algorithm as
x’p+l’ = (I - L)_‘Ux(P’ + (I - f,-‘Jy.
(5)
The matrix I - L is invertible since det (I - L) = 1. Theorem 7. The iterative solution x of Kx = y
scheme
(5) always
produces
a sequence
which converges
to the
The updating procedure does in fact accelerate the convergence of the solution process without complicating the iterative method. This is the result of the next theorem. Theorem 8. The iterative scheme (4) hus at least as fast a rate of convergence 6.
ERROR
us (3)
ANALYSIS
In what has already been described, the procedure x(P+
11 =
QX(P’
+
d,
d = (I - L)-‘Jy,
is given for approximating, with any preassigned degree of accuracy, the solution x of Kx = Y. If a matrix norm 11.(1is found (one always exists) for which 11Q II < 1, then a standard error estimate can be made if desired 1123. We start by noting the representations x(P)
=
i i=O
x = t
Q’4 Q’d.
i=O
If II.II is a vector norm which is compatible with the above matrix norm, then 11 x(P) - x/I = 11i-t+1
Q'dIi
=
IIQ""f. Q'dIl 5 IIQll4"Ildll
=
lIQ(IP+llldll(l
- lIQII)-‘3
izo ItQIl'
Iterative inversion of single compartment matrices
325
where the last expression can be computed from knowledge of 1)Q I]and 1)d )I. In particular, to find an index p such that IIx(P) - xl] < ( for prescribed c > 0, the inequality
IIQIIP+llldtIU - IIQIO-' -=E can be solved for the least such integer p. It should be noted that when performing practical computations on most digital computers the inherent peculiarities of the machine’s floating point arithmetic may require some alteration of the above-mentioned p value.
7. CONDITIONING
OF A SEC MATRIX
In this section we examine the vulnerability or sensitivity of a solution: If K is slightly perturbed, how great is the effect on x = K-‘y? A standard measure of sensitivity is the condition number for a matrix A, cond (A) = 1)A II II A-’ )I 2 1, using an appropriate matrix norm [13, p. 2711. If cond (A) zzz 1, the matrix A is wellconditioned. The larger cond (A) is, the more ill-conditioned is the system since the solution is more vulnerable to round-off errors. As we will see from the calculations below, a convenient matrix norm to use is the infinity norm. This matrix norm for an n x n real matrix A 3 [I,] can be calculated by
11 A 11 max i laijl[8]I
=
I
j=l
Since each element in the mth row of K- l is -k,-,‘, then
where the notation K- ’ = [ki; ‘I] h as b een used. Thus the condition number
Another lower bound on this condition number can be used when the maximal eigenvalue of K and the excretion rate constant k,, are known (as in the inverse problem). Retain the ordering of the eigenvalues {&} of K as 11112 IArl r ... 2 ]&I, and use the fact that 11K (Im 2 I& I. Then cond,(K)
8. A SIMPLE
2 II, Ink&,‘.
EXAMPLE
In our previous work the algorithm using updated components for finding x = K- ‘y has been developed. An illustration of the methodology is now presented. The example is to solve the system Kx = y where Y = (Yl,Y,,Y#
= (-LO,W
and ‘-0.805 K=
0.7 1
0.1
0.2 -0.7 0.5
0.4 0.7 - 1.1
1,
J
DAVID H. ANDERSON
326
The matrix K is seen to be invertible by an appeal to Theorem problem is x = (200, 400, 2OO)r. Note that k,, = 0.005,
5. The exact solution
of this
(1KIIw = 2.1,
so that cond,(K) 2 1260. Using the method developed in this paper, we get an immediate advantage over other methods since we have a closed form solution for the mth component in the solution x z (x,, x2, X3)=: 3 xl
Now the J matrix
=
-
C i=l
y,/k,,
= 200.
is (with c E 1.1)
[;
so that the starting
---l;l
_-;I,
vector x(O) E (x(p), x$0’, x:“‘)’ = Jy
is (2OO,O,O)r. The E matrix
e13 -0 -
e12 --0
-
ell --0
has entries
ezl = k2Jc = 7/11
ez2 = (c + k,,)/c
e31 = k3Jc = l/11
eJ2 = k32/c = 5111
ez3 = k,,Jc = 7111
= 4/11
e33 = (c + k33)/c = 0
(note that if c is defined to be max [kit] z k,,, then always eq4 = 0, which simplifies algorithm even further). Thus the computational scheme for this problem becomes
the
x’p+ l’ = (4x’p’ + 7x’9’ + 1400)/l 1 x(g+ 1’ = (Sxcp’ I’ + 200)/l 1 where the first equation in (4) has been disregarded since we know the exact component x1 = 200. These expressions can be easily programmed and the resulting sequences generated by them converge rapidly (- 30 iterations) and monotonically from the initial vector up to X) = 399.99999, using the stopping
criterion I(x(p+l’-
9.
xj = 199.99999,
COMPARISON
WITH
x(p’I~l/I]x(p+l’(J1 DIRECT
< 10-S.
NUMERICAL
METHODS
Thirty-two test systems Kx = y in the range 3 I n I 8 (the usual size of the problems we have encountered) and 1 I cond,(K) 5 lo6 with known exact solutions were devised. A subroutine LEQSEC (Linear EQuation solver for SEC systems), utilizing the above updating iterative scheme, was written specifically for linear systems involving SEC matrices and has been extensively tested and compared to the direct methods of Gaussian elimination, Crout reduction with iterative improvement, and singular value decomposition [ 141. Gaussian elimination failed to give “exact solutions” (answers satisfactorily close to the known exact solutions) in 4 of the 32 cases. Crout reduction failed to produce exact solutions for 6 cases; of these Crout reduction with iterative improvement was able to compute exact solutions for only 2 cases. In all of the trials that failed, the matrix was classified by the program as “algorithmically singular” because too small a pivot element was encountered. When the test matrices were ill-conditioned to this extent, the standard stabilization strategy of singular value decomposition was used. The LEQSEC algorithm is also numerically stable
Iterative
inversion
of single compartment
matrices
327
since p(Q) < 1 [ 121. On none of the remaining 4 test cases did singular value decomposition give more accurate answers than did LEQSEC, and in some cases the solution was much less accurate. In fact, no SEC system has been encountered thus far on which the LEQSEC algorithm failed to yield the exact solution and so LEQSEC would seem to be a preferable method for solving SEC systems on machines which have a comparatively small range of representable numbers or machines which are known to be numerically “dirty”. In this testing process, all computations were done on a CDC CYBER 70, model 72, or a DEC system- 10. The respective subroutines used were DECOMP-SOLVE [ 151, LEQTlF and LEQT2F in the International Mathematical and Statistical Libraries package, and SVD [ 163. The experimental comparisons did shed some light on the one shortcoming of LEQSEC, which was the relatively larger time required on a few problems to obtain a solution as compared to the direct methods. However, since these problems were exactly those most illconditioned, slow convergence of the iterative scheme was to be expected. Work is now in progress to modify the algorithm so as to further increase its rate of convergence. SUMMARY In the study of biological systems a frequently used model is an open compartmental system having only one compartment from which elimination takes place. Let K be the compartmental matrix associated with the mathematical model of such a system, i = Kz + b. This article treats the theory of solution of the simultaneous set of linear equations Kx = y which arises in the analysis of the model. Simple sufficient conditions for the invertibility of K are found, properties of the inverse matrix are discovered, and an iterative scheme is developed for which convergence to x is guaranteed. The algorithm’s progress is accelerated and its storage requirements are made more modest by an updating procedure. An example comparing this computational technique to other methods is given, and a general error analysis is discussed. Acknowledgements - In the formulation of this report, Drs. Myron Ginsberg and G. Milton Wing made many valuable comments. I would also like to express my appreciation to Sonja Sandberg who did substantial work in the programming and computing necessary for the research.
REFERENCES 1. P. B. Bright, The volumes of some compartment systems with sampling and loss from one compartment, Bull. Math. Biol. 35, 69-79 (1973). Anulysis in Biology und Medicine. Elsevier, Amsterdam (1972). 2. J. A. Jacquez, Compartmentul 3. J. 2. Hearon, Theorems on linear systems, N. Y. Acud. oJSci. APIIIA/S108, 36-68 (1963). 4. D. Fife, Which linear compartmental systems contain traps? Muth. Biosci. 14, 311-315 (1972). in single-exit compartmental systems, Applied Nonlineur Anulysis, 5. D. H. Anderson, The volume of distribution V. LAKSHMKANI-HAM Ed. Academic Press, New York (1979). 6. J. M. Ortega, Numerical Analysis. A Second Course. Academic Press, New York (1972). I. L. B. Rail, Computational Solution of Nonlineur Operator Equutions. John Wiley, New York (1969). Theory of Matrices. Academic Press, New York (1969). 8. P. Lancaster, 9. B. K. Shah, Data analysis problems in the area of pharmacokinetics research. Biomen+cs 32, 145% I57 (1976). C/in. Pharmucol. Therupeur. 8, 201-218 (1967). IO. J. G. Wagner, Use of computers in pharmacokinetics. 11. S. Venit, The convergence-of Jacobiand Gauss-Seidel iteration. Muth. Maya& 48, 163-167 (1975). Theory und Practice. Addison-Wesley, Reading, Ma12. E. K. Blum, Numerical Analysis und Computation. ssachusetts (1972). Academic Press, New York (1976). 13. G. Strang, Linear Algebra and ifs Applicufions. 14. C. H. Farnsworth, Comparing numerical methods for solving linear systems of equations involving single-exit compartmental systems, Unpublished report (1978). 15. G. E. Forsythe, M. A. Malcolm and C. B. Moler, Computer merhods for muthematicul computurions. Prentice-Hall, Englewood Clifis, N.J. (1977). 16. G. H. Golub and C. Reinsch, Singular value decomposition and linear systems solutions. Num. Muth. 14, 403-420 (1970). 17. G. Forsythe and C. B. Moler, Computer Solution ofLinear Algebraic Systems. Prentice-Hall, Englewood Cliffs, N.J. (1967). About the Author - DAVID H. ANDERSON received his B.S. degree in Mathematics from the University of Mississippi in 1961. From 1961 to 1965 he was an instructor of mathematics at the U.S. Naval Nuclear Power School, Vallejo, California, involved in training the crews of the nuclear ships of the
DAVID H. ANDERSON
328
U.S. Navy. After working one summer for Exxon in Baton Rouge, Louisiana, he took the M.A. degree in Mathematics from the University of California, Berkeley, followed by a year as Assistant Professor of Mathematics at Millsaps College, Jackson, Mississippi. Earning his Ph.D. in Mathematics from Duke University in 1969, he joined the faculty at SMU in Dallas. At present, Dr. Anderson is an Associate Professor in the Department of Mathematics at SMU and in the Department of Medical Computer Science at the University of Texas Health Science Center at Dallas, and is a Visiting Lecturer for the Mathematical Association of America. During the period January-May 1979, he had a visiting position with the Division of Applied Mathematics, Brown University, Providence, R.I. He has published several articles dealing with mathematical modeling, numerical analysis, and algorithm design in compartmental analysis and population mathematics.
APPENDIX Proof of Theorem 2 Let u be the I x n row vector of zeros except for a I in the mth position. each entry unity. Then the condition
Take w to be the I x n row vector with
leads to the equation WK = - k,,u. Since K is nonsingular,
then UK-’ = (-l/k&v.
But UK-’ is just the mth row (vector) of K-l. - l/k,.
(6) (6) says that each entry in the mth row of K-’
Hence equation
is
Proof of Theorem 3 The (i,j)-element
of E is 0 eij
-i
i=m
if
1 + k,Jc
if
i # m,
i=j.
k,k
if
if
i#j
m,
Since each entry in the mth row of E is zero, then E is a singular reducible matrix. Also since c 2 1kiilfor all i, and k, r 0 whenever
i #j, then each e,, 2 0. Thus E is a nonnegative
matrix.
Now the jth column
sum of E (j # m) is
= I - k,+, which is less than or equal to unity. The mth column
sum of E is
=(-km, + ILmlVc, a value strictly
less than one since k,
z 0. These sums tell us I
(in particuiar, nonempty.
observe
j#H;
C e,j < l i=,
if
,$, eu=
if jeH.
that m $ H.) This is equivalent
I
to saying
(/E
//, <: 1 ifHisempty,and
l/E/,
=
I ifHis
Proof of Theorem 4 Let K be invertible. We wish to prove that p(E) < 1. The magnitude of any matrix norm of E. Thus, by the last theorem,
of any eigenvalue
of E cannot exceed the value
p(E) 5 IIgIIt 5 1 By way of contradiction, let us assume that p(E) = I. Appealing to the PerronFrobenius theory for nonnegative matrices [8]. we can say that there is a nonnegative eigenvalue r of the matrix E such that r = p(E) = 1. Since rl - E = JK, then 0 = det(rf
- E) = det(JK)
= det J
det K.
Iterative
inversion
of single compartment
329
matrices
Since J is nonsingular, we must conclude that det K = 0. But this contradicts the invertibility of K ; hence p(E) < 1. On the other hand, suppose p(E) < 1. Given any 8 > 0, there is a matrix norm 11 t 11such that
11 El\ I p(E) + 463. Let I: = (1 - p(E))/2. Then III - JK )I s
IIE 11< 1 since p(E) < 1. Thus J is an approximate
inverse of K and so K is invertible.
Proojof nteorem 5 From the proof of Theorem
3, we know that C elj + 1
d H. Since this is a strict inequality
wheneverj
1
for every j, there exists an E > 1 which will serve for all j $ H such that
2 ei,+tz matrix
1 et,< 1. id/
r$fl
Let D s Cd,,] be the n x n diagonal
eij <
idl
i$ll
with d,, =
i
if jeH
;
if j$H
Define the matrix C = [c,,] 2 0 through the similarity transformation C s DED- ‘. Thejth row of C,j E H, is the jth row of E multiplied by ~zand the jth column of C, j E H, is the jth column of E multiplied by l/c. Otherwise the entries of C are the same as those of E. Form the Gerschgorin disks [S] based on the columns of C. This will give information about the eigenvalues of E since C is similar to E. For each j = 1,2,. , n, we have the disk of complex numbers 0,s
z:]z - c,,] I {
For the nonnegative
maximal
can be written r I
I
as
1 eilo + 1 eij& id”
k”
of the similarity transformation. Consider the last sum just given. Since in the summation only when i # j,. Recall that when i # j, ei, =
By hypothesis, Thus
“i,
r = p(E) of E, there exists at least one index j, for which r E bj0. Thus
eigenvalue
Suppose j0 E H. Then the last inequality
Via the construction quantity eilo occurs
1 i#i
0
if
i=m
k,/c
if
i#m
there is at least one i, 4 H for which k ,0,0 # 0 (sincej,, E H). This index i,, cannot
i $ H and& E ff, the
be m because k,,,,” = 0.
and thus the sum
is nonzero.
obtains
As a result, the strict inequality
and thus
r<
1 $I/
eij, +
1 eijo= j,
i+tl
eih =
1.
Now suppose j0 $ H. Then we can write
But since j, 4 H, the right hand side of this inequality must be less than unity. Thus for either possibility, j, d H or j, E H, we have p(E) = r < 1. Hence by Theorem proof is complete.
4, K is invertible
and the
DAVID H. ANDERSON
330 Prwf
of Theorem 6
The matrix
J has negative
entries everywhere
on its mth row and main diagonal.
Thus the matrix
x,=.J+JV, which is less than or equal to J, has negative elements in the same positions. Now concentrate on row numbers i # m and column numbers j # i. Introduce the matrix notation J = [Ji,]. For these i and j values, the (i, j)-element of EJ is i
p=*
ei,,Jpj = e,,,,J,,,, + e,lJ,, = - k,&k,,)-’
- kifi- L
ifj # tn. and is
%Jnn = - h,kkm) - 1 ifj = tn. In either case, the (i,j)-entry
of EJ is negative since, by hypothesis,
ki, > 0 for all
i = I, 2,.
, n, i # m. Thus
K-‘
(5) is of the general
linear iterative
form
xc&‘+11=
QX’P’
+
d
for the square matrix Q = (I - L)- ‘U. In order that the sequence of iterates converge, that p(Q) < 1 [6]. We prove that p(Q) < 1 by showing that the series
it is necessary
and sufficient
I + Q + Q’ + Q3 + . . converges. Now
Q = (I - L)-‘u
= (I + L + L* + . . . + L”-‘)U
is a non-negative matrix [6]. The series expansion triangular. Define A= Then E = I - JK imples B-’
-J
-K,
= (I - E)A-’
= (I - L)-‘(-J)
for (I Es
L)-’
terminates
with L’-’
since L is strictly
lower
-J-‘(I-L).
and so
= (I - L)-‘(I
- E)A-’
= (I - L)-‘((I
- L) - U)A-’
=A-‘-(I-L)-‘(IA-‘=(I-Q)A-‘. Since (I-L)-‘=I+L+L2+...+L’-t20 and -J
2 0, then 8-l
2 0. We also observe
that A-’
in [6, p. 1191.
= - K - ’ h 0. From here on, the proof is similar to that
Proof of Theorem 8. For the general
linear iterative
process
r
x@‘+ ’ G Mx@ + d,
p =
0, 1,2,.
,
(7)
the asymptotic convergence factor is p(M), provided p(M) < 1[6]. For (4), the matrix A4 z Q and p(Q) < 1 by Theorem 7. In the case of (3). M = E and via Theorem 4, we have p(E) 4 1. Thus the asymptotic convergence factors of the 2 iterative schemes are p(Q) and p(E), respectively. The smaller p(M) is, the faster the rate ofconvergence of the process (7) [6]. Clearly then, we need to prove that p(Q) 5 p(E). Let n = p(Q). Since Q is a nonnegative matrix, the Perron-Frobenius theory guarantees the existence of a nonzero vector IAsuch that (I - L)-
’ Uu z Qu = qu [S-J.
Then (?$ + U)u = vu, which says that q is a nonnegative
eigenvalue
of the matrix
qL + U. Thus
P(rlL + U) 2 n. SinceLzO,U20,andO~rr<
l,then O
Hence
the desired
inequality.