- Email: [email protected]
A recurrence among the elements of functions of triangular matrices
A Recurrence among the Elements of Functions of Triangular Matrices* B. N. Parlett
Departmmt of Mathematics and Computer Science Division Department of Electrical Engineering and Computer Sciences and the Electronics Research Laboratory University of California at Berkeley
Submitted by Hans Schneider
ABSTRACT A simple relation exists among the elements of G(T) when @I is an analytic function and T is triangular. This permits the rapid build up of @b(T) from its diagonal. Moreover, exp(B -‘A) can be formed without inverting B.
1.
A RELATION
AMONG
Any square complex
THE
ELEMENTS
matrix B can be factored
OF QJ(T)
into
B = PTP*, where P is unitary (P-l = P*) and T is upper triangular [2]. Here M* denotes the conjugate transpose of M. If B is real, then P can be taken as real (i.e., orthogonal) while T is real and block upper triangular, with diagonal blocks of orders 1 and 2.l The decomposition is not unique, but is robust in the face of rounding errors. Let + be a function analytic in some region in the complex plane containing B’s spectrum. Then
G(B)=PG(T)P*, *The author gratefully acknowledges partial support from Office of Naval Research Contract NOOO14-69-A-0200-1017. ‘The real Schur form is sometimes called the Wintner-Mumaghan form.
LINEAR ALGEBRA AND ITS APPLICATIONS 14, 117-121 (1976) Q American Elsevier Publishing Company, Inc., 1976
117
B. N. PARLETT
118
and so 9(B) can be obtained from c$(T) and P. A recurrence relation among the elements of G(T) can be derived from two simple lemmas whose proofs may be found in [4]. LEMMA 1. If T is block upper triangular, then c$(T) is also block upper triangular and has the same block structure. LEMMA 2. Let + be an analytic function and B be a square matrix. If G(B) is well defined, then it commutes with B.
Figure 1 illustrates Lemma 1. x x
T=
x x
x x
x x
x x
x x
x
x
xi
ooxxxx 0
0
0
F=+(T)I=
YYYYYY YYYYYY 0 0
Y
Y
Y
Y
0
Y
Y
Y
OOOYYY
oooxxx
0
00000x
oooooy
0
FIG. 1
THEOREM 1. Let T = (7J be block upper triangular. Then F E+(T) the same structure and, for r < s,
has
s-r-l
T-J,--CT,=
Proof.
2 Pr,r+J-r+k,s-Tr,s-kFs-k,s)*
k=O
By Lemma 2, FT-
TF-0.
By Lemma 1, F is also block upper triangular. On writing out the (r,s) block n of FT- TF and rearranging terms the recurrence relation is obtained. The diagonal blocks of F can be obtained from F,= +( T,). Then the upper triangular blocks of F may be computed from the formula, by diagonals, moving out from the main diagonal as indicated in Fig. 2. When the formula is used in this way, the right hand side will be known, and so the elements of F, are determined by a set of simultaneous linear
119
FUNCTIONS OF TRIANGULAR MATRICES x x x x--x-x-x-
x
x
x x
x x
x
x
x x
1
+
X
T
4
x
zi z x x
FIG. 2
equations T, F, - F, T, =
R.
A unique solution exists provided that T, and T, have no eigenvalues in common. An efficient algorithm exists for solving equations of this form: T, is reduced to lower triangular form and T, is reduced to upper triangular form by separate unitary similarity transformations, and the resulting system of equations is triangular. A program is given in [l]. So much for the general case. When I3 is complex, the Schur form T is triangular, and the scalar quantities F, are obtained with one division. When B is real, the F, may be 1 X 1, 2 X 1, 1 X 2, or 2 X 2, depending on the block structure of the real Schur form T. The first three cases are easy to program, but the 2 X 2 case involves four unknowns and the nuisance of having to allow for some form of pivoting. However, if the diagonal elements of T,, and of T, are made equal, then the solution can be expressed in a convenient closed form. Details are given in [5]. The formula breaks down when T, and T, have eigenvalues in common, in particular when T,= T,. Moreover implementations in finite precision arithmetic can be expected to give inaccurate answers when T, is very close to T,. When F is n X n, fewer than in” multiplications are required to form F from its diagonal using the formula, whereas 7n3 multiplications are needed to form T and P from B; see [4, 51. In the confluent case (T, = T,) there are alternative formulae which can be used to compute the corresponding F, before the general formula is applied. See [4] for more details.
120
2.
B. N. PARLETT FUNCTIONS
OF MATRIX PENCILS
An initial value problem which arises in many applications has the form: find u(t) for t > 0 such that Bti(t)=Au(t),
u (0) = ug>
where A and B are constant. When B is invertible, the solution is, formally, u
(t) = exp( B - ‘At) uo.
When working with floating point arithmetic of fixed precision, the explicit inversion of B may sometimes provoke unnecessarily large roundoff errors. The next theorem shows how exp(B -‘At) can be computed without inverting B. The storage requirements of this phase are doubled, but the “extra” arithmetic operations (for F) are balanced by avoiding the formation of B -‘At. The method is not confined to the exponential function. The first step is to reduce A and B simultaneously to upper triangular form by unitary transformations. This is accomplished by the QZ algorithm
E31. THEOREM 2. Let A and B - ’ be upper triangular and let + be a scalar function for which F = C#J( B -‘A) and F= +(AB L’) are well defined. The elements of the upper triangular matrices F and F are given by
and, for T< s, s-r-l
a~fVs-_&ass=G-Ma~
+
_
jZl
(fr,r+k’t+k,s-ar.s-kfs-k,s),
k-l
Proof.
Since AB -‘=
B (B -lA)B -l, F”= BFB -l. By Lemma 2,
B -‘AF=FB
-‘A,
AB-%=gAB-‘.
FUNCTIONS OF TRIANGULAR MATRICES
121
Thus AF-iA=O,
BF-
fB=O.
(1)
Moreover, F and 1; are upper triangular (by Lemma 1). On equating the (r,s) element on each side of (l), the given relations are obtained. a By tolerating the storage of the auxiliary matrix r”, the pair F, f can be built up together from the diagonal outwards, without ever inverting B. When the two linear equations do not have a unique solution, i.e., when a,.,b,, - a,, b, = det [
aI7 b ~
1
ass b,,
=Oy
the confluent form, involving derivatives of +, must be used.
3.
IMPLEMENTATION
A paper presenting programs based on these results and discussing their numerical properties will be submitted elsewhere. REFERENCES R. H. Bartels and G. W. Stewart, Algorithm 432, solution of the matrix equation AX + XB = C, Commun. Assoc. Comput. Mach. 15 (1972), 826826. C. C. MacDuffee, The Theory of Matrices, Chelsea, 1956. C. B. Moler and G. W. Stewart, An algorithm for generalized matrix eigenvalue problems, SIAM 1. Numer. And. 10 (1973), 241-256. B. N. Parlett, Computation of functions of triangular matrices, Memo. No. ERL-M481, Univ. of Cahf. at Berkeley, Nov. 1974. S. P. Chan, R. Feldman, and B. N. Parlett, A program to compute the condition numbers of matrix eigenvalues without computing eigenvectors, Memo. No. ERL-517, Univ. of Cahf. at Berkeley, May 1975. Received 27 May 1975
Departmmt of Mathematics and Computer Science Division Department of Electrical Engineering and Computer Sciences and the Electronics Research Laboratory University of California at Berkeley
Submitted by Hans Schneider
ABSTRACT A simple relation exists among the elements of G(T) when @I is an analytic function and T is triangular. This permits the rapid build up of @b(T) from its diagonal. Moreover, exp(B -‘A) can be formed without inverting B.
1.
A RELATION
AMONG
Any square complex
THE
ELEMENTS
matrix B can be factored
OF QJ(T)
into
B = PTP*, where P is unitary (P-l = P*) and T is upper triangular [2]. Here M* denotes the conjugate transpose of M. If B is real, then P can be taken as real (i.e., orthogonal) while T is real and block upper triangular, with diagonal blocks of orders 1 and 2.l The decomposition is not unique, but is robust in the face of rounding errors. Let + be a function analytic in some region in the complex plane containing B’s spectrum. Then
G(B)=PG(T)P*, *The author gratefully acknowledges partial support from Office of Naval Research Contract NOOO14-69-A-0200-1017. ‘The real Schur form is sometimes called the Wintner-Mumaghan form.
LINEAR ALGEBRA AND ITS APPLICATIONS 14, 117-121 (1976) Q American Elsevier Publishing Company, Inc., 1976
117
B. N. PARLETT
118
and so 9(B) can be obtained from c$(T) and P. A recurrence relation among the elements of G(T) can be derived from two simple lemmas whose proofs may be found in [4]. LEMMA 1. If T is block upper triangular, then c$(T) is also block upper triangular and has the same block structure. LEMMA 2. Let + be an analytic function and B be a square matrix. If G(B) is well defined, then it commutes with B.
Figure 1 illustrates Lemma 1. x x
T=
x x
x x
x x
x x
x x
x
x
xi
ooxxxx 0
0
0
F=+(T)I=
YYYYYY YYYYYY 0 0
Y
Y
Y
Y
0
Y
Y
Y
OOOYYY
oooxxx
0
00000x
oooooy
0
FIG. 1
THEOREM 1. Let T = (7J be block upper triangular. Then F E+(T) the same structure and, for r < s,
has
s-r-l
T-J,--CT,=
Proof.
2 Pr,r+J-r+k,s-Tr,s-kFs-k,s)*
k=O
By Lemma 2, FT-
TF-0.
By Lemma 1, F is also block upper triangular. On writing out the (r,s) block n of FT- TF and rearranging terms the recurrence relation is obtained. The diagonal blocks of F can be obtained from F,= +( T,). Then the upper triangular blocks of F may be computed from the formula, by diagonals, moving out from the main diagonal as indicated in Fig. 2. When the formula is used in this way, the right hand side will be known, and so the elements of F, are determined by a set of simultaneous linear
119
FUNCTIONS OF TRIANGULAR MATRICES x x x x--x-x-x-
x
x
x x
x x
x
x
x x
1
+
X
T
4
x
zi z x x
FIG. 2
equations T, F, - F, T, =
R.
A unique solution exists provided that T, and T, have no eigenvalues in common. An efficient algorithm exists for solving equations of this form: T, is reduced to lower triangular form and T, is reduced to upper triangular form by separate unitary similarity transformations, and the resulting system of equations is triangular. A program is given in [l]. So much for the general case. When I3 is complex, the Schur form T is triangular, and the scalar quantities F, are obtained with one division. When B is real, the F, may be 1 X 1, 2 X 1, 1 X 2, or 2 X 2, depending on the block structure of the real Schur form T. The first three cases are easy to program, but the 2 X 2 case involves four unknowns and the nuisance of having to allow for some form of pivoting. However, if the diagonal elements of T,, and of T, are made equal, then the solution can be expressed in a convenient closed form. Details are given in [5]. The formula breaks down when T, and T, have eigenvalues in common, in particular when T,= T,. Moreover implementations in finite precision arithmetic can be expected to give inaccurate answers when T, is very close to T,. When F is n X n, fewer than in” multiplications are required to form F from its diagonal using the formula, whereas 7n3 multiplications are needed to form T and P from B; see [4, 51. In the confluent case (T, = T,) there are alternative formulae which can be used to compute the corresponding F, before the general formula is applied. See [4] for more details.
120
2.
B. N. PARLETT FUNCTIONS
OF MATRIX PENCILS
An initial value problem which arises in many applications has the form: find u(t) for t > 0 such that Bti(t)=Au(t),
u (0) = ug>
where A and B are constant. When B is invertible, the solution is, formally, u
(t) = exp( B - ‘At) uo.
When working with floating point arithmetic of fixed precision, the explicit inversion of B may sometimes provoke unnecessarily large roundoff errors. The next theorem shows how exp(B -‘At) can be computed without inverting B. The storage requirements of this phase are doubled, but the “extra” arithmetic operations (for F) are balanced by avoiding the formation of B -‘At. The method is not confined to the exponential function. The first step is to reduce A and B simultaneously to upper triangular form by unitary transformations. This is accomplished by the QZ algorithm
E31. THEOREM 2. Let A and B - ’ be upper triangular and let + be a scalar function for which F = C#J( B -‘A) and F= +(AB L’) are well defined. The elements of the upper triangular matrices F and F are given by
and, for T< s, s-r-l
a~fVs-_&ass=G-Ma~
+
_
jZl
(fr,r+k’t+k,s-ar.s-kfs-k,s),
k-l
Proof.
Since AB -‘=
B (B -lA)B -l, F”= BFB -l. By Lemma 2,
B -‘AF=FB
-‘A,
AB-%=gAB-‘.
FUNCTIONS OF TRIANGULAR MATRICES
121
Thus AF-iA=O,
BF-
fB=O.
(1)
Moreover, F and 1; are upper triangular (by Lemma 1). On equating the (r,s) element on each side of (l), the given relations are obtained. a By tolerating the storage of the auxiliary matrix r”, the pair F, f can be built up together from the diagonal outwards, without ever inverting B. When the two linear equations do not have a unique solution, i.e., when a,.,b,, - a,, b, = det [
aI7 b ~
1
ass b,,
=Oy
the confluent form, involving derivatives of +, must be used.
3.
IMPLEMENTATION
A paper presenting programs based on these results and discussing their numerical properties will be submitted elsewhere. REFERENCES R. H. Bartels and G. W. Stewart, Algorithm 432, solution of the matrix equation AX + XB = C, Commun. Assoc. Comput. Mach. 15 (1972), 826826. C. C. MacDuffee, The Theory of Matrices, Chelsea, 1956. C. B. Moler and G. W. Stewart, An algorithm for generalized matrix eigenvalue problems, SIAM 1. Numer. And. 10 (1973), 241-256. B. N. Parlett, Computation of functions of triangular matrices, Memo. No. ERL-M481, Univ. of Cahf. at Berkeley, Nov. 1974. S. P. Chan, R. Feldman, and B. N. Parlett, A program to compute the condition numbers of matrix eigenvalues without computing eigenvectors, Memo. No. ERL-517, Univ. of Cahf. at Berkeley, May 1975. Received 27 May 1975