A Power Method for Computing Square Roots of Complex Matrices

A Power Method for Computing Square Roots of Complex Matrices

JOURNAL OF MATHEMATICAL ANALYSIS AND APPLICATIONS ARTICLE NO. 213, 393]405 Ž1997. AY975517 A Power Method for Computing Square Roots of Complex Mat...

188KB Sizes 3 Downloads 67 Views

JOURNAL OF MATHEMATICAL ANALYSIS AND APPLICATIONS ARTICLE NO.

213, 393]405 Ž1997.

AY975517

A Power Method for Computing Square Roots of Complex Matrices Mohammed A. Hasan Department of Electrical Engineering, Colorado State Uni¨ ersity, Fort Collins, Colorado 80523 Submitted by Harlan W. Stech Received August 22, 1995

In this paper higher order convergent methods for computing square roots of nonsingular complex matrices are derived. These methods are globally convergent and are based on eigenvalue shifting and powering. Specifically, it is shown for each positive integer r G 2, a convergent method of order r can be developed. These algorithms can be used to compute square roots of general nonsingular complex matrices such as computing square roots of matrices with negative Q 1997 Academic Press eigenvalues.

1. INTRODUCTION A square root of a complex matrix A g C m= m is defined to be any matrix B g C m= m such that B 2 s A, where C is the field of complex numbers. If all eigenvalues of an m = m matrix A are distinct, then the matrix equation X 2 s A generally has exactly 2 m solutions. This follows from the fact that A is diagonalizable, i.e., there exists a similarity matrix U such that A s UDUy1 , where D s diagŽ l1 , . . . , l m . and thus B s Uy1 'D U, where 'D s diagŽŽy1. i1 l1 , . . . , Žy1. i m l m ., and i k s 0 or 1 for k s 1, 2, . . . , m. However, if A has multiple eigenvalues, the number of solutions will be different from 2 m as shown next. Let m s 2, then without 2 2 loss of generality, we can assume that A s l 02 or A s l 12 .

'

'

0

Assume that A s l I, then the family l 2

½

cos Ž u . sin Ž u .

l

sin Ž u . ycos Ž u .

0

5

0 F u - 2p

infinite set of square roots of A. On the other hand, if A s

l

2

0

l

forms an 1

l2

, then

l A has only two square roots given by "0l "1r2 provided that l / 0. "l Unlike the square roots of complex numbers, square roots of complex

393 0022-247Xr97 $25.00 Copyright Q 1997 by Academic Press All rights of reproduction in any form reserved.

394

MOHAMMED A. HASAN

matrices may not exist. For example, when l s 0 in the last matrix no square root exists. From this observation, it is obvious that for 2 = 2 matrices, the equation B 2 s A / 0 has a solution if and only if A has a nonzero eigenvalue. To understand the structure of solutions of the equation B 2 s A for b b m s 2, let B s b11 b12 , and A s aa1121 aa1222 such that B 2 s A. Then we 21

22

have the following four equations bi1 b1 j q bi2 b 2 j s a i j for i, j s 1, 2 which are equivalent to F s 0, where

¡b

q b12 b 21 y a11 b11 b12 q b12 b 22 y a12 F Ž b11 , b12 , b 21 , b 22 . s b 21 b11 q b 22 b 21 y a21 b 21 b12 q b 22 b 22 y a22 .

~

11 b 11

¢

Ž 1.

The Jacobian of this system can be shown to be 2 b11 ­ F Ž b11 , b12 , b 21 , b 22 . b12 Js s b 21 ­ Ž b11 , b12 , b 21 , b 22 . 0

b 21 b11 q b 22 0 b 21

b12 0 b11 q b 22 b12

0 b12 . b 21 2 b 22

It can be verified that < J < s 4Ž b 11 q b 22 . 2 Ž b 11 b 22 y b 12 b 21 . s 4 Trace 2 Ž B .< B <. Here the notation < J < denotes the determinant of J, and < < TraceŽ B . s Ý m is1 bii . Since A is nonsingular, it follows that B / 0 and therefore J is nonsingular if and only if b11 q b 22 / 0. Now assume that A is nonsingular and let B0 s

0 b11

0 b12

0 b 21

0 b 22

0 0 be a solution of the equation B 2 s A such that b11 q b 22 / 0. Since J is 0 0 0 0 . nonsingular at Ž b11 , b12 , b 21 , b 22 , it follows from the implicit function theorem that B0 is the only solution. From the eigendecomposition of A indicated before, one can see that there are at least four square roots of A. The implicit function theorem guarantees exactly four square roots with nonzero traces. These square roots of A which have nonzero traces are referred to as functions of A w1x. Essentially, B is a function of A if B can be expressed as a polynomial in A. Now if b11 q b 22 s 0, then it follows from Ž1. that a11 q a22 s l2 , a12 s a21 s 0, i.e., A is diagonal of the form l2 I. In this case the equation

SQUARE ROOTS OF COMPLEX MATRICES

395

B 2 s A has a two-dimensional family of solutions given by

½

"'l2 y rs

r

s

.'l y rs 2

5

r, s g C.

Note that when r s s s l sinŽ u ., we get the one-parameter family

½l

cos Ž u . sin Ž u .

sin Ž u . ycos Ž u .

5 0 F u - 2p which is described before.

The following result provides conditions on the eigenstructure of A which ensure the existence of square roots which are functions of A. PROPOSITION 1 w1x. Let A be nonsingular and its p elementary di¨ isors be coprime, that is, each eigen¨ alue appears in only one Jordan block. Then A has precisely 2 p square roots, each of which is a function of A. Several computational methods of square roots of complex matrices have been reported in the literature. In w2x, the Newton-Raphson method was used for computing the principal square root of a complex matrix. An accelerated algorithm for computing the positive definite square root of a positive definite matrix was presented in w3x. A matrix continued fraction method was presented in w4x. The matrix sign algorithm was developed in w5x. A Schur method for computing square roots was developed in w6x. Fast stable methods for computing square roots were also presented in w7, 8x. It is noted in almost all of the above methods either a linear or quadratic convergence can be obtained. In this paper, higher order convergent methods of order r G 2 will be derived. The essence of these methods is a process whereby a sequence of matrices which in the limit converges to a square root of A is generated. This process involves creating gaps between the magnitudes of eigenvalues of different square roots of A so that for sufficiently high powers the eigenvalues will become decoupled. This is similar in principle to well-known methods such as those of Graeffe, Bernoulli, and the qd algorithm for solving polynomial equations in that these methods are based on eigenvalue powering. For a survey of some of these methods the reader is referred to w9, 10x and the references therein. Let S be a set of commutative and thus simultaneously diagonalizable matrices. In the sequel, the notation l i Ž X . denotes the ith eigenvalue of the square matrix X g S relative to a fixed similarity matrix which diagonalizes the set S . The notation s Ž A. denotes the set of eigenvalues of A. The symbol R is used to denote the set of real numbers and 5 A 5 denotes any vector norm of the matrix A.

396

MOHAMMED A. HASAN

2. DERIVATION OF THE MAIN RESULTS In the next theorem we will generate a sequence which converges to a square root of a square matrix. THEOREM 2. Let A g C m= m be a nonsingular matrix. Let r be a positi¨ e integer such that r G 2 and define A k and Bk recursi¨ ely as follow. Let r

A kq 1 s

r l 2l l Ary2 Bk A , k 2l

Ž 2.

r ly1 2 lq1 l Ary2 Bk A . k 2l q 1

Ž 3.

Ý ls0

ž /

and r

Bkq 1 s

Ý ls0

ž

/

Then there exists an a g C such that Bk is nonsingular for all sufficiently large k. Set X k s Bky1A k , then the initial guess A 0 s aIm and B0 s Im the sequence X k con¨ erges to a square root W of A. Moreo¨ er, r

r X kq 1 " W s By1 kq1 Bk Ž X k " W . ,

Ž 4.

i.e., if the sequence X k con¨ erges, it is rth order con¨ ergent to W. Additionally, r lim k ª` Ay1 kq1 A k s I

y1 lim k ª` Bkq1 Bkr s I.

and

Proof. Let W be any square root of A, i.e., W 2 s A and show by induction that k

A k " BkW s Ž aI " W . r .

Ž 5.

Clearly Ž5. holds for k s 0. Assume that Ž5. holds for the positive integer k. Then

Ž aI " W . r

kq 1

r

s Ž A k " BkW . s

r

r l 2l l A ry2 Bk A k 2l

ž / Ýž

Ý

ls0

r

"

ls0

r ly1 2 lq1 l A ry2 Bk A W k 2l q 1

/

s A kq 1 " Bkq1W , where the last equality follows from Ž2. and Ž3.. Hence Ž5. is true for the integer k q 1. This shows that Ž5. is true for each nonnegative integer k.

397

SQUARE ROOTS OF COMPLEX MATRICES

The nonsingularity of A implies that there exists an a g C such that < l j Ž aI q W .< ) < l j Ž aI y W .<, for j s 1, . . . , m. From Ž5. we have Ž aI q k k W . r s A k q BkW and Ž aI y W . r s A k y BkW. Solving the last two equations for A k and Bk yields Ak s

1

 Ž aI q W . r 2

k

k

q Ž aI y W . r 4 ,

Ž 6.

and Bk s

1

 Ž aI q W . r 2

k

k

y Ž aI y W . r 4 Wy1 .

Ž 7.

Note that the order of matrix multiplication is immaterial in this case since the choice X 0 s aI implies that the elements of the sequence  X k 4`ks1 commute with each other. It should also be noted that for sufficiently large k both A k and Bk are invertible since < l j Ž aI q W .< ) < l j Ž aI y W .<, j s 1, . . . , m. Thus multiplying both sides of Ž5. by By1 yields k k

r y1 By1 s 2W I y  Ž aI q W . k A k q W s Bk Ž aI q W .

s 2W I q Ž aI q W .

½

y1

½

y1

Ž aI y W . r

rk

y1 Ž aI y W . 4 q  Ž aI q W . Ž aI y W . 4

2rk

k

4

y1

q??? .

5

Since < l i Ž aI q W .< ) < l i Ž aI q W .<, for i s 1, . . . , m, it follows from the last equation that k

r y1 lim By1 s 2W. k A k q W s lim Bk Ž aI q W .

kª`

kª`

lim k ª` Bky1A k

Therefore, s W. To prove Ž4. we have from the relation A kq 1 " Bkq1W s Ž A k " BkW . r that r

y1 r y1 By1 kq 1 A kq1 " W s Bkq1 Bk Ž Bk A k " W . ,

or equivalently r

r X kq 1 " W s By1 kq1 Bk Ž X k " W . .

Finally, the last conclusion follows from r r By1 kq 1 Bk s  Ž aI q W .

s

as k ª `.

kq 1

qŽ aI y W . r

½ ž I y Ž Ž aI q W . ½ ž I y Ž Ž aI q W .

kq 1

y1

4  Ž aI q W . r

y1

Ž aI y W . .

y1

Ž aI y W . .

r kq 1

rk

k

q Ž aI y W . r

k

r

4

y1

5

r

/ ªI Q.E.D.

398

MOHAMMED A. HASAN

An analogous algorithm to that of Theorem 2 can be obtained by setting X k s Bky1A k as shown next. THEOREM 3. Let A be a nonsingular matrix and let r be a positi¨ e integer such that r G 2. Let r

Ck s

r X kry2 lAl , 2l

Ž 8.

r X kry2 ly1Al , 2l q 1

Ž 9.

Ý ls0

r

Dk s

Ý ls0

ž

ž /

/

m= m and set X kq 1 s Dy1 , k C k . Then there exists an initial matrix X 0 s aI g C ` a g C for which the sequence  X k 4ks1 con¨ erges to a square root W of A. Moreo¨ er, r

X kq 1 " W s Dy1 k Ž Xk " W . ,

Ž 10 .

i.e., if the sequence X k con¨ erges it is rth order con¨ ergent to W. Proof. Theorem 3 follows directly from Theorem 2 by setting X k s Bky1A k . Q.E.D. In Theorems 2 and 3 and if r s 2, quadratically convergent algorithms are obtained as follows. COROLLARY 4. Let A be an m = m nonsingular complex matrix and define the following sequence A kq 1 s A2k q Bk2 A, and Bkq1 s 2 A k Bk with A 0 s aI and B0 s I. Then for some a g C , the sequence Bky1A k con¨ erges quadratically to 'A . Moreo¨ er, if we set X k s By1 k A k , then the iteration Ž X k2 y A. with X 0 s aI con¨ erges quadratically to some X kq 1 s X k y 12 Xy1 k 'A and X kq 1 y 'A s 12 Xy1 Ž X k y 'A . 2 . k Ž X k2 y A. of the above Remark. The iteration X kq 1 s X k y 12 Xy1 k corollary is exactly the Newton iteration for solving X 2 y A s 0. Several variants of Newton’s method and their implementations were presented in w2x, where the matrix A is assumed diagonalizable. Note that the algorithm of Corollary 4 and the Newton method are equivalent. However, the derivation of the above algorithm lends itself to the development of other methods whose convergence is of any prescribed order. A cubically convergent algorithm can be derived by setting r s 3 in Theorems 2 and 3 as in the next result. COROLLARY 5. Let A be an m = m nonsingular complex matrix and define the following sequence A kq 1 s A3k q 3 Bk2 A k A and Bkq1 s 3 A2k Bk q Bk3 A, with A 0 s aI and B0 s I. Then for some a g C , the sequence Bky1A k

399

SQUARE ROOTS OF COMPLEX MATRICES

con¨ erges cubically to a square root of A. Moreo¨ er, if we set X k s Bky1A k then the sequence satisfies the recursion X kq 1 s X k y 2Ž3 X k2 q A.y1 Ž X k2 y A. which with the initial guess X 0 s aI con¨ erges cubically to 'A and X kq 1 y 'A s 2Ž3 X k2 q A.y1 Ž X k y 'A . 3. 3. ANALYSIS OF CONVERGENCE To analyze the convergence of the sequence  X k 4`ks1 we formulate the iteration in Theorem 2 as a fixed point iteration as established in the next result. Assume there exists an a g C such that < l i Ž aI y for i s 1, . . . , m. Then the sequence X k defined in Theorem 3 is rth order con¨ ergent, i.e., THEOREM 6.

'A .rli Ž aI q 'A .< - 1,

r

X kq 1 y 'A s O Ž X k y 'A . , or equi¨ alently, there exists a constant S G 0 such that 5 X kq 1 y 'A 5 F S 5 X k y 'A 5 r . Proof. Let Z g C m= m be any square matrix of A and define C Ž Z . s r Ý rls0 2rl Z ry2 lAl and DŽ Z . s Ý rls0 2 l q Z ry2 ly1Al. Then 1

ž /

ž

/

r

C Ž Z . " D Ž Z . 'A s Ž Z " 'A . . Define F Ž Z . s DŽ Z .y1 C Ž Z .. Then F: C m= m ª C m=m and r

F Ž "'A . s

½Ýž ls0

r 2l q 1

/ r

Ž "'A .

s 2 ry1 Ž "'A . Ž "'A .

½

ry2 ly1

y1

Al

r

ls0

y1 y1

5 ½2

r 2l

5 ½ Ý ž /Ž

ry1

"'A .

ry2 l

Al

5

r Ž "'A . 5 s "'A .

Therefore 'A and y 'A are fixed points for the function F. Now, the function F can be written as F Ž Z . s .'A q D Ž Z .

y1

Ž C Ž Z . " D Ž Z . 'A . r y1 s .'A q D Ž Z . Ž Z " 'A . , F Ž Z . " 'A s DŽ Z .y1 Ž Z " 'A . r. Note that when Z is a scalar it

hence can be shown

d dZ

F Ž Z . s rD Ž Z .

y1

Ž Z 2 y A.

ry1

.

400

MOHAMMED A. HASAN

Therefore if j 2 y A s 0, then F Ž j . s j and Ž d lrdZ l . F Ž j . s 0, for l s 1, 2, . . . , r y 1. It follows that there is an R ) 0 and a neighborhood, SR Ž j . s  Z g C m= m < 5 Z y j 5 - R4 , of j and a constant K, 0 F K - 1 such that FŽ Z. y FŽ j . F K 5 Z y j 5 for all Z g SR Ž j . [  Z < g C m= m 5 Z y j 5 F R4 Ži.e., F is a contractive mapping in SR Ž j ... Hence for any X 0 g SR Ž j . the generated sequence X kq 1 s F Ž X k ., k s 1, 2, ??? has the properties X k g SR Ž j . and 5 X k y j 5 F K k 5 X 0 y z 5 for k s 1, 2, ??? . This shows that the iteration X kq1 s F Ž X k . converges to j for any initial guess X 0 in SR Ž j . for which X 0 A s AX 0 . To prove the last conclusion of the theorem, we have r

Ck " D k'A s Ž X k " 'A . , where

C k [ C Ž X k . s Ý rls 0

r 2l

ž /X

ry 2 l l A, k

and

Ž 11 . Dk [ DŽ Xk . s

r Ý rls0 2 l q X kry2 ly1Al. Since DŽ'A . s Ž2'A . ry1 is nonsingular, D k is 1 nonsingular for all sufficiently large k. It follows that there exists a 5 F S for all sufficiently large k. Therefore constant S G 0 such that 5 Dy1 k 5 X kq 1 y 'A 5 F S 5 X k y 'A 5 r. This proves the r th order convergence of Xk . Q.E.D.

ž

/

4. CONTINUED FRACTION EXPANSION There are many formulas in the literature that approximate square roots of numbers and matrices using continued fraction expansion w4x. In this section, we express the results of Theorem 3 in continued fraction form. THEOREM 7. Let A be a nonsingular matrix and let W be a square root of A. Then for each a g C for which < l i Ž aI q W .< ) < l i Ž aI y W .<, i s 1, . . . , m, W has the representation W s aI q Ž A y a2 I . 2 aI q Ž A y a2 I .

½

 2 aI q Ž A y a2 I .  2 aI q ??? 4 y1 4

y1 y1

5

. Ž 12 .

401

SQUARE ROOTS OF COMPLEX MATRICES

It can easily be verified that the r th order truncation of the continued fraction of Ž12. yields DŽ aI .y1 C Ž aI ., where C Ž aI . s Ý rls0 2rl a ry2 lAl , and

ž /

DŽ aI . s

r Ý rls0 2 l q 1

ž

/a

ry2 ly1

l

A . When A is a scalar, Formula Ž12. reduces

to

'A

saq

A y a2 2 a q Ž A y a2 . r Ž 2 a q Ž A y a2 . r Ž 2 a q ??? . .

,

Ž 13 .

which in the case a s 1 becomes that of w4x. Note that when A is positive, Ž13. converges for any a g C with nonzero real part. The free parameter a provides some flexibility in choosing the initial guess in that the closer a is to 'A the more accelerated is the convergence. Remark on the Choice of the Initial Guesses. In this remark we provide special cases where conditions on the initial matrix X 0 can be imposed to ensure convergence. Assume that all eigenvalues of A are not on the negative real line. Then each eigenvalue of 'A has nonzero real part. The initial guess A 0 s aI with a ) 0 Ž a - 0. forces the sequence X k to converge to a square root of A with eigenvalues having positive Žnegative. real part. However, if some of the eigenvalues of 'A are on the imaginary axis, then the choice X 0 s Ž a q ib . I, where a ) 0 and b ) 0, will lead to a sequence which converges to a square root of A having eigenvalues with nonnegative real parts. These observations show in particular that if A is positive Žor negative definite. then the methods of Theorems 2 and 3 converge for any initial guess X 0 of the form X 0 s aI, a g R Ž ia g R . , where i s 'y 1 . This is an improvement over the algorithm in w2x which is applicable only to matrices which have no negative eigenvalue.

5. COMPUTATIONAL RESULTS To demonstrate the performance of the algorithms proposed in this paper we consider in this section some examples to show the behavior of some of these methods in finite-precision arithmetic. These computations were carried out on approximately seven decimal digit accuracy. For the purpose of comparison we will apply Theorem 3 to four examples. The notation X k, r will denote the computed square root using k iterations and r th order method. The error is measured in the Frobenious norm 5 X k,2 r y A5 F .

402

MOHAMMED A. HASAN

EXAMPLE 1. Consider the complex 3 = 3 matrix

As

8 q 15i

y1 y 3i

y4 y 9i

3 y1 y 3i

3 5 q 9i

3 y1 y 3i

3 y4 y 9i

3 y1 y 3i

3 8 q 15i

3

3

3

.

Applying five iterations of the algorithm of Theorem 3 with r s 2, X 0 s Ž1 q i . I, yields 1.938066 q 1.123146i X5, 2 s y0.2334079 q y 0.2188985i y0.6059740 q y 0.4491569i

y0.2334078 q y 0.2188987i 1.565499 q 0.8928874i y0.2334077 q y 0.2188986i y0.6059739 q y 0.4491573i y0.2334077 q y 0.2188989i , 1.938066 q 1.123145i

which agrees with the exact square root to five decimal places, i.e., 5 X5,2 2 y A 5 F s O Ž10y6 .. Comparable accuracy can also be achieved using only four iterations with a third order method with the same initial guess in which case we obtain 5 X 4,2 3 y A 5 F s O Ž10y6 .. EXAMPLE 2. Consider the 4 = 4 matrix w2x 5 4 As 1 1

4 5 1 1

1 1 4 2

1 1 . 2 4

The eigenvalues of this matrix are lŽ A. s  1, 2, 5, 104 and 2-norm condition number k 2 Ž A. s 10. When X 0 s I and 4 iterations are used with r s 3, we obtain 1.988520 0.9885170 X 4, 3 s 0.1852420 0.1852418

0.9885167 1.988519 0.1852420 0.1852418

0.1852414 0.1852422 1.917762 0.5035481

0.1852417 0.1852421 , 0.5035480 1.917762

and 5 X 4,2 3 y A 5 F s O Ž10y6 .. We also note that 5 X 4,2 2 y A 5 F s O Ž10y3 . while 5 X5,2 2 y A 5 F s O Ž10y7 ..

SQUARE ROOTS OF COMPLEX MATRICES

403

EXAMPLE 3. Consider the 4 = 4 near singular matrix 0 1.31 As 1.06 y2.64

.07 y0.36 2.86 y1.84

.27 1.21 1.49 y0.24

y.33 0.4 . y1.34 y2.01

The eigenvalues of this matrix are lŽ A. s  0.03, 3.03, y1.97 " i4 . Applying the iteration of Theorem 3 with r s 3 and X 0 s Ž1 q i . I we obtain 0.2441971 1.317404 X 7, 3 s 1.0659751 E y 02 y0.6756389

y9.0401828 E y 02 1.181676 0.1509080 y1.982304

0.1997281 0.2569374 1.370772 0.3442460 y8.5186012 E y 02 0.8455525 , y1.248934 y0.1964209

with error 5 X 7,2 3 y A 5 F s O Ž10y4 .. However, when r s 2, similar accuracy can be obtained with the same X 0 as 0.2437900 1.317573 X6, 2 s 1.0955708 E y 02 y0.6754091

y9.0409003E y 02 1.181780 0.1507968 y1.982364

0.1997136 0.2569094 1.370820 0.3442680 y8.5148215E y 02 0.8455746 , y1.249036 y0.1964408

with error measured in Frobenious norm as 5 X6,2 2 y A 5 F s O Ž10y4 .. EXAMPLE 4. Consider the 3 = 3 matrix w2x 1 A s y1 0

1 2 1

y2 1 , y1

which has the eigenvalues lŽ A. s  y1, 1, 24 , i.e., this matrix has a negative eigenvalue. Thus the initial matrix X 0 s aI should be chosen so that

404

MOHAMMED A. HASAN

the imaginary part of a is nonzero. Using X 0 s Ž1 q i . I, we obtain 1.0285987 y 0.1666667i X 9, 3 s y0.4142135 y 0.0000000i y0.0285954 q 1.1666666i

0.4714045 y 0.3333333i 1.4142135 q 0.0000000i 0.4714045 y 0.3333333i y1.0285954 q 1.16666667i 0.4142135 q 0.0000000i , y0.0285954 q 1.1666666i

with error 5 X5,2 2 y A 5 F s O Ž10y1 0 .. If a method of order three is used, then 5 X 3,2 3 y A 5 F s O Ž10y8 .. In summary, some square roots of nonsingular complex matrices can be computed using only initial matrices of the form X 0 s Ž a q ib . I. The case b / 0 is required for matrices having negative eigenvalues.

6. CONCLUSION A new set of algorithms for computing square roots of nonsingular complex matrices was developed. Given any positive integer r G 2 we presented a systematic way of deriving an r th order convergent square root algorithm. The convergence and its speed are largely affected by the ratios R i s < l i  X 0 y 'A 4
ACKNOWLEDGMENT The author thanks the referees for their helpful remarks and suggestions which improved the quality of this work.

SQUARE ROOTS OF COMPLEX MATRICES

405

REFERENCES 1. N. J. Higham, Computing real square roots of a real matrix, Linear Algebra Appl. 88r89 Ž1987., 405]430. 2. N. J. Higham, Newton’s method for the matrix square root, Math. Comp. 46, No. 174 Ž1986., 537]549. 3. E. D. Denman, Roots of real matrices, Linear Algebra Appl. 36 Ž1981., 133]139. 4. L. S. Shieh and N. Chahin, A computer-aided method for the factorization of matrix polynomials, Appl. Math. Comput. 2 Ž1976., 63]94. 5. E. D. Denman and A. N. Beavers, The matrix sign function and computation of systems, Appl. Math. Comput. 2 Ž1976., 63]94. 6. A. Bjorck and S. Hammarling, A Schur method for the square root of a matrix, Linear Algebra Appl. 52r53 Ž1983., 127]140. 7. W. D. Hoskins and D. J. Walton, A fast method of computing the square root of a matrix, IEEE Trans. Automat. Control AC-23, No. 3 Ž1978., 494]495. 8. W. D. Hoskins and D. J. Walton, A fast, more stable method for computing the pth root of positive definite matrices, Linear Algebra Appl. 26 Ž1979., 139]163. 9. P. Henrici, ‘‘Applied and Computational Complex Analysis,’’ Vol. 1, Wiley, New York, 1974. 10. A. S. Householder, ‘‘The Numerical Treatment of a Single Nonlinear Equation,’’ McGraw]Hill, New York, 1970. 11. M. A. Hasan, Higher order convergent algorithms for computing nth roots of complex matrices, submitted for publication.