A NOTE ON THE EVALUATION
OF MATRIX POLYNOMIALS
by
Charles Van Loan
Department of Computer Science
Cornell University
Ithaca, New York 14853
TR78--H353
A NOTE ON THE EVALUATION OF NATRIX POLYNONIALS
Charles Van Loan
Abstract:
The problem of evaluating a polynomial p(x) in a matrix
A arises in many applications, e.g. the Taylor approximation of
e .			O(? n3)
A			The			algorithm of Paterson and Stoc?meyer has the
drawbac? that it requires O( ? n2) storage, where q is the degree
of p and n is the dimension of A. An algorithm which greatly
reduces this storage requirement without undue loss of speed is
presented.
Author?s Current Address:
Department of Computer science
Cornoll Univorsity
Ithaca, Now York l4?53
A NOTE ON THE EVALUATION OF NATRIX POLYNOMIALS
Charles Van Loan
One method for evaluating the matrix pelynomial
p(A) = b01 + b1A + ... + b
q
is to use Homer's rule:
R , A  RflXfl
- 1
= bql
For k = 1,2g
5k = ASk?l + bg?kI
whereupon 5 =?p(A). Ignoring low order terms in n, this algorithm re-
q
quires 2n2 storage and (q-l)n3 multiplications.
In an interesting paper by Paterson and Stockmeyer 21 , it
is shown how this computation can be rearranged so that only O[?q n31
multiplications are needed. The main idea behind their algorithm is
to apply Homer's rule to p(A) regarded as a polynomial in AS
where a . For exanple, if q = 9 , then
p(A) = (b9I) (A33 + (b8A2+ b7A + b6I) (A32 +
2 + + (A3) + (b2A2 + b1A + b0I)
(b5A b4A b3I)
Onco A2 and A3 aro computod, then only 2n3 multiplications are
2-
required to evaluate p(A) an evidenced by the followin(j "l?rner
regrouping"
p(A) = A3 A3[ b9A3 + (b8A2 + b7A + b61)] + (b5A2 + b4A + b31)]
+ (b2A + b1A + b01)
3
Hence, by this technique, only 4n nultiplications are required-
to evaluate p(A). This is in contrast to the count of Sn3 as-
sociated with the "ordinary" application of Homer's rule.
In general, if 5 is any integer satisfying
l<s?<q
and
r = [q?5?			= the integral part of
then
(1)
where
p(A) =			f			Bk (AS)k
k
Bk = bsk+sl A5?1 +			+ b5)c+1 A + b5)cI			(k = 0r-l)
B			= b			+ ... + bsr+1 A + b - I
r			q			sr
and
--H3--H
Applying Horner?s rulo to (1) we obtain:
F0			B
r
For k			1r
Fk = (A5)Fk?l + Br?k
where F = p(A). The fact that there can be a significant reduction
r
in work as a result of these manipulations is made clear when this
iteration is spelled out in greater detail:
Algorithm I
1.			?k = Ak			(k = 0,15)
2.			F0 = b Yq?? + b			Y + ... + bsr+1Y1 + b ?
q			q-l?q-sr-l			sr 0
3. For k = 1,
= YFk?l +b			Y +			+b			Y
Fk			5			s(r-k)+s-l 5-1			s(r-k) 0
where F = p(A).
r
This is essentially Algorithm "B" in [21 applied to the ev?alu-
ation of matrix polynomials. Noting that (s-l)n3 multiplications
are required to compute the yk and that there is one matrix--Hmatrix
multiplication per iteration in step 3, we conclude that for a given
3
a, Algorithm I requires w1ts?n???rmultipli?ations where
w1(s)
--H4--H
s +. [q,5?			-			2
a +			-			1
if 5 divides q
otherwise
Setting a = 1 in Algorithm I renders the ordinary version of
Homer's rule. Setting 5 ? 1 approxinately nininizes w1(s)
and we then find w1( 1 ? ? 2 ?
A discouraging feature of Algorithn I is that it requires
(s+l)n3 storage. (Arrays are required for A,A2A5 , and the
current" Fk o+) However, a slightly less efficient algorithn
which requires only 3n2 storage results if p(A) is conputed one
column at a tine. To derive this algorithn, let a denote the
i-tn column of the identity natrix. Fron (I) we have
p(A)e. =			r			(A5k
)			(Bker)
k=O
which nay be?conputed as follows:
= B e,
0			r ?
For k - 1r
- (A5)fk?1l? + ? kel
--H5--H
where ?r??? .P(A)e? , the j-th column.of p(A). Notice that once
the vectors Ae. , A2e3 A9?1e. are computed, then the
3 3
vectors Skej			k = 0r are "free"
Bk?j = bsk+s?l (???ie.) + + bsk+l. (Ae ) + (k = 0,
3			j			bs?e3
B e.			=			b			???S??			+			+ bsr+i Ae.			+			b e.
r 3			q			3			3			sr 3
We therefore obtain the following algorithm:
Al orithin II
1.			Y = AS
2.			For j = 1n
y0			=
For k =			1s-l
= A Yk?-?t
= bq ?(?)			+			+ + bsr y0(?)
0 q-sr bsr+i y1(i)
For k = 1r
+			+ b
= v ?k(-?l? + bs(r?k)+s?l Ys?l s(r--Hk) C
whereupon p(A) = ( ?r?1)			fr(?) 1
--H6--H
Since arrays are needed for only A , A8 , and p(A), Algor-
ithin II requires 3n2 storage. A careful operation count for 8 ? 1
shows that w2(s)n			multiplications are required where
2(.s --H 1) + ???5? --H 2 if a divides g
w2(s) =
2(5-- 1) - + [q?5? - 1 otherwise
This work count is approximately minimized by setting a			? 3
whereupon w2 (s) cc 2? - about .1.5 times the work associated with
Algorithm I
At the expense of one extra nxn array, A5 can be computed
via its binary expansion
A5 = (A?:20 ) (A?l21)-. (A?m 2??			? =			m?			?k2k
in tilog2(s)n multiplications where a ? [1,2] depends upon the
number of nonzero - (If a is an exact power of two, the extra
storage is not required.) When this feature is incorporated in Al-
gorithm 11, the work count becomes w2(s)n3 where
= a log2(s) --H 1 + w1(s)
These ideas have been incorporated in a subroutine for approx--H
A with a
xnatimg the matrix exponential e			truncated Taylor series
eA			?			?k/k:			= T (A)
k=O			q
For reasons based on the rnachine precision used (lo?l5) the norm
of A and the desired accuracy? q was chosen to be 16. (See [11
A
for an analysis of Taylor approximation to e .) Setting s 4 and
using binary squaring to compute A4 Aigorithn II requires 7n3
multiplications and 3n2 storage to evaluate T16(A). This is less
than half the work associated with traditional Homer evaluation
and only marginally slower than the 6n3 multiplications needed by
Algorithm I which requires 5n2 storage for this problem.
Acknowledgenents
The author greatly profited from discussions he had with
Dr. Len Berman who pointed out reference f21.
References
tl? C,3.Moler?and C.Van Loan, Nineteen Dubious Way to Conpute
the Matrix Exponential, to appear, SlAM Review.
M.S.Paterson and I,.Stockmeyer, On the Nunber of Nonscalar
Multiplications Necessary to Evaluate Polynornials,
SlAM J.Cornputing, 2(l973),60-66.
E2?