The Moore–Penrose inverses of matrices over quaternion polynomial rings

Linear Algebra and its Applications 475 (2015) 45–61

Contents lists available at ScienceDirect

Linear Algebra and its Applications www.elsevier.com/locate/laa

The Moore–Penrose inverses of matrices over quaternion polynomial rings Liji Huang b , Qing-Wen Wang a,∗ , Yang Zhang b a

Department of Mathematics, Shanghai University, Shanghai, PR China Department of Mathematics, University of Manitoba, Winnipeg, MB R3T 2N2, Canada b

a r t i c l e

i n f o

Article history: Received 8 May 2014 Accepted 1 February 2015 Available online xxxx Submitted by W.B. Wu MSC: 15A09 68W30 Keywords: Quaternion polynomial matrix Moore–Penrose inverse Leverrier–Faddeev algorithm

a b s t r a c t In this paper, we define and discuss the Moore–Penrose inverses of matrices with quaternion polynomial entries. When the Moore–Penrose inverses exist, we prove that Leverrier– Faddeev algorithm works for these matrices by using generalized characteristic polynomials. Furthermore, after studying interpolations for quaternion polynomials, we give an efficient algorithm to compute the Moore–Penrose inverses. We developed a Maple package for quaternion polynomial matrices. All algorithms in this paper are implemented, and tested on examples. © 2015 Elsevier Inc. All rights reserved.

1. Introduction In 1843 Sir Rowan Hamilton discovered the algebra H of real quaternion, which is a four-dimensional non-commutative algebra over real number field R with canonical basis 1, i, j, k satisfying the conditions: i2 = j2 = k2 = ijk = −1, * Corresponding author. E-mail addresses: [email protected] (Q.-W. Wang), [email protected] (Y. Zhang). http://dx.doi.org/10.1016/j.laa.2015.02.004 0024-3795/© 2015 Elsevier Inc. All rights reserved.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

46

that implies ij = −ji = k, jk = −kj = i, and ki = −ik = j. The elements in H can be written in a unique way: α = a + bi + cj + dk, where a, b, c and d are real numbers, that is, H = {a + bi + cj + dk | a, b, c, d ∈ R}. The conjugate of √ α is defined as α ¯ = a − bi − cj − dk, and the norm |α| is |α| = αα ¯ . It is well-known −1 2 that H is a skew field, that is, for 0 = α ∈ H, α = α/|α| ¯ . The study of polynomials with quaternion coefficients may go back to Niven [17] and [18] in the early 1940’s. Due to the non-commutativity of H, there are several forms of quaternion polynomials depending on the positions of coefficients. In this paper, we will use the following Definition 1.1, which puts the coefficients on the left side of a variable x. Also they are called regular quaternion polynomials in [5] and quaternion simple polynomials in [19]. Some properties have been discussed (see, for example, [14] and [20]). Definition 1.1. A quaternion polynomial f (x) over H is defined as f (x) = an xn + · · · + a1 x + a0 ,

ai ∈ H, i = 0, . . . , n,

where x commutes element-wise with H. The set of all quaternion polynomials in x is denoted by H[x]. The degrees, leading terms and leading coefficients are defined in a common way. Also the degree of 0 is n m i j defined as −∞. For f (x) = i=0 ai x , g(x) = j=0 bj x ∈ H[x], the addition and multiplication are defined as: 

max{n,m}

f (x) + g(x) =

k=0

(ak + bk )xk ,

f (x)g(x) =

n+m 



ai bj xk ,

k=0 i+j=k

where some ai , bj may be assumed to be zero for simplicity. The conjugate of f (x) = an xn + · · · + a0 ∈ H[x] is defined as f (x) = a ¯ n xn + · · · + a ¯0 , and has the following properties: Lemma 1.2. (See [21].) Let f, g ∈ H[x]. Then (i) f g = g¯f¯ (ii) f f¯ = f¯f ∈ R[x], where R are reals, (iii) if f g ∈ R[x], then f g = gf . The quaternion polynomials and matrices with quaternion polynomial entries have been widely studied with many applications in past decades. For example, in [27], the Fast Fourier Transform for the product of two quaternion polynomial has been discussed with the complexity analysis. In [5], Damiano et al. studied Gröbner basis theory for the ring of quaternion polynomials and explored how to compute the module syzygy. Smith–McMilian forms of quaternion polynomial matrices are defined in [22] and some applications to dynamical systems are given.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

47

For matrices over commutative rings, it is well-known that the Moore–Penrose inverses have been defined and explored for many years (see, for example, [3,12,13,23]). This motivates us to consider the Moore–Penrose inverses of quaternion polynomial matrices. The structure of this paper is as follows. In Section 2, we give the definition of the Moore–Penrose inverse for quaternion polynomial matrices and discuss some basic properties including a sufficient and necessary condition for the existence of the Moore–Penrose inverse. Leverrier–Faddeev algorithm is extended to quaternion polynomial matrices in Section 3 by using generalized characteristic polynomials. In Section 4, we discuss the interpolation problems for quaternion polynomials and give an efficient algorithm to compute the Moore–Penrose inverses. Finally using algorithms in this paper and our Maple package some examples are given in Section 5. 2. Definitions and basic properties In this section, we first define the Moore–Penrose inverse for matrices over H[x], and then discuss some basic properties as well as the sufficient and necessary conditions for the existence of the Moore–Penrose inverses. Let H[x]m×n denote the set of all m × n matrices with entries from H[x]. For A ∈ H[x]m×n , the conjugate  A ofA is defined as A = (Aij ). If A = P + Qj with P, Q ∈

C[x]m×n , then χA =

P −Q

Q P

∈ C[x]2m×2n denotes the complex adjoint of A. Moreover,

AT , A∗ ∈ H[x]n×m denote the transpose and the conjugate transpose of A, respectively. In particular, for all f ∈ H[x], f ∗ = f . More properties of H[x]m×n can be found, for example, in [21,22]. Definition 2.1. A† ∈ H[x]n×m is called a Moore–Penrose inverse of A ∈ H[x]m×n if it is a solution of the following system of equations: ∗

∗

AXA = A, XAX = X, (AX) = AX, (XA) = XA. Note that we require that A† must be in H[x]n×m . Therefore unlike the matrices over fields, the Moore–Penrose inverses for some quaternion polynomial matrices might not exist. Recall that A ∈ Hm×m is unitary if AA∗ = A∗ A = Im . Based on the properties of quaternions, one can easily prove the following elementary properties that will be often used throughout this paper. Proposition 2.2. Let A ∈ H[x]m×n and B ∈ H[x]n×l . Then ∗

(i) (AB) = B ∗ A∗ and AA∗ = (AA∗ )∗ .  ∗  ∗ † (ii) If A has a Moore–Penrose inverse A† , then (A∗ ) = A† , A† A† A∗ = A† =  ∗ A∗ A† A† and A† AA∗ = A∗ = A∗ AA† . (iii) If A has a Moore–Penrose inverse A† , then A† is unique.

48

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

(iv) Let A have the Moore–Penrose inverse A† . If U ∈ Hm×m is a unitary matrix, then † (U A) = A† U ∗ . Next we will give conditions for quaternion polynomial matrices to have Moore– Penrose inverses. We need to prove some lemmas first. n×1 It is easy to see that A ∈ Hm×n [x] induces an additive homomorphism from H [x] m×1 n×1 m×1 to H [x] , that is, for all P, Q ∈ H [x] , A (P + Q) = AP + AQ ∈ H [x] . By the definition of Moore–Penrose inverses and Proposition 2.2, it is easy to prove the following lemma: Lemma 2.3. Let A ∈ H [x] have the Moore–Penrose inverse A† . Consider A as a hon×1 m×1 momorphism from H [x] to H [x] . Then Image(A) = Image(AA∗ ) = Image(AA† ) ∗ ∗ and Image(A ) = Image(A A) = Image(A† A). m×n

m×m

Lemma 2.4. If E ∈ H [x] E ∈ Hm×m .

is a symmetric projection, that is, E = E 2 = E ∗ , then

Proof. Let f1 , . . . , fm be the entries on the first row of E. From E = E ∗ , we may assume that f1 = f1 = 0. Then by E = E 2 , we have

f1 = f1 f1 +

m 

fi fi = f12 +

i=2

m 

fi fi .

i=2

Since f1 = f1 , the leading coefficient of f12 is a positive real number. Note that the m leading coefficient of i=2 fi fi is also a positive real number. Thus,

deg



f12



 ≥ deg f1 = deg

f12

+

m 

fi fi

i=2

m      2 fi fi ≥ deg f12 . = max deg f1 , deg

i=2

m  This shows that f1 ∈ H. Furthermore, 0 = deg f1 = deg i=2 fi fi and the leading coefficients of {fi f¯i } (fi = 0) are positive reals imply that fi ∈ H for all 1 ≤ i ≤ m. The same discussions can be done for the other rows of E. Therefore, E ∈ Hm×m . 2 Due to the non-commutativity of quaternions, there are two types of eigenvalues: right eigenvalues and left eigenvalues. Right eigenvalues have been studied extensively (see, for example, [1,4,15]). We will work with right eigenvalues towards our main result. For convenience, we will just use the term “eigenvalue” from now on. The following result is well-known and very useful.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

49

Lemma 2.5. (See [26].) A ∈ Hm×m is hermitian, that is, A = A∗ , if and only if there exists a unitary matrix U ∈ Hm×m such that U ∗ AU = diag(d1 , . . . , dm ), where di are the eigenvalues of A. Now we are ready to give conditions that quaternion polynomial matrices must satisfy in order to have Moore–Penrose inverses. The following theorem is well-known in some cases, see, for example, [3,24]. Here is an analog version for quaternion polynomial matrices. Theorem 2.6. Let A . Then A has the Moore–Penrose inverse A† if and only  ∈ H [x] A1 A2 r×r if A = U 0 0 with U ∈ Hm×m unitary and A1 A∗1 + A2 A∗2 a unit in H [x] with r ≤ min {m, n}. Moreover, m×n

A† =



−1

A∗1 (A1 A∗1 + A2 A∗2 ) −1 A∗2 (A1 A∗1 + A2 A∗2 )

0 0

U ∗.

Proof. (=⇒) If A has the Moore–Penrose inverse A† , then  2  ∗ AA† = AA† AA† = AA† = AA† . By Lemma 2.4, AA† ∈ Hm×m . AA† is hermitian and hence, by Lemma 2.5, there exists a unitary matrix U ∈ Hm×m such that U ∗ AA† U = D where D is diagonal. Since D2 = (U ∗ AA† U )(U ∗ AA† U ) = U ∗ AA† AA† U = U ∗ AA† U = D, the diagonal  entries  of D are either 1 or 0. Therefore, we can rearrange the rows of U so Ir 0 that D = 0 0 with r ≤ min {m, n}.

Set A = U∗ A. By Lemma 2.2, A has its own generalized inverse A† and  Ir 0 A1 A2 †  A A = 0 0 . Set A = A A , for arbitrary quaternion polynomial matrices 

3

r×r

4

r×(n−r)

(m−r)×r

(m−r)×(n−r)

A1 ∈ H [x] , A2∈ H [x] ,A and A4 ∈ H [x] .Since    3 ∈ H [x]  Ir 0 A1 A2 A1 A2 A1 A2   †   A = AA A = 0 0 = 0 0 , we must have A = 0 0 and A3 A4     ∗ ∗ A1 A1 + A2 A2 0 B1 0 † therefore A A∗ = Similarly, A = . for some B1 and B2 . By 0 0 B2 0 Lemma 2.3, 

∗







†

Image(A A ) = Image(A ) = Image(A A ) = Image This implies the surjectivity of A1 A∗1 + A2 A∗2 on H [x] r×r unit in H [x] and

r×1

A = U A = U



A1 0

A2 0

.

Ir 0

0 0

.

. Therefore, A1 A∗1 + A2 A∗2 is a

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

50

Next, we have that:  ∗ † † A† = A† A† A∗ = A† (A∗ ) A∗ = A∗ (A A∗ )



∗ ∗ −1 −1 (A1 A∗1 + A2 A∗2 ) 0 A1 (A1 A∗1 + A2 A∗2 ) A1 0 = = ∗ ∗ ∗ ∗ −1 A2 0 0 0 A2 (A1 A1 + A2 A2 )

0 0

,

which gives: A† =



−1

A∗1 (A1 A∗1 + A2 A∗2 ) −1 A∗2 (A1 A∗1 + A2 A∗2 )

0 0

U ∗.

(⇐=) The converse can be proved by direct computation. 2 3. Leverrier–Faddeev algorithm Leverrier–Faddeev algorithm has been used to compute the Moore–Penrose inverses for many years. We refer the reader to [2,6,11,25] for more details. In this section, we define the characteristic polynomial for quaternion polynomial matrix A by using AA∗ . In particular, we prove that the coefficients of this characteristic polynomial are reals. Then we show that Leverrier–Faddeev algorithm works very well for quaternion polynomial matrices. Given A ∈ H[x]m×m . An element λ ∈ H is called an eigenvalue of A if there exists a vector X ∈ Hm×1 [x] such that AX = Xλ. Lemma 3.1. Let A ∈ H[x]m×n . Then eigenvalues of AA∗ are real. Proof. Let B = AA∗ and λ ∈ H be an eigenvalue of B with corresponding eigenvector T X = ( x1 · · · xm ) = 0 such that BX = Xλ. Then X ∗ BX = X ∗ Xλ. Note that B = B ∗ . We have that X ∗ BX = λ∗ X ∗ X, and thus ∗

X ∗ Xλ = λ∗ X ∗ X = (X ∗ Xλ) , that is, ⎞ x1     ⎜ . ⎟ ∗ ∗ ¯i xi ) λ = (( x (¯ x1 , · · · x ¯m ) ⎝ .. ⎠ λ = ( x ¯i xi ) λ) = λ∗ ( x ¯i xi ) = λ∗ ( x ¯ i xi ) . xm ⎛

By Lemma 1.2, 0 = λ ∈ R. 2



x ¯i xi ∈ R [x]. The above equation gives λ = λ∗ , which implies

Cayley–Hamilton theorem for quaternion matrices has been extensively discussed. A survey can be found in [26]. Next we define the characteristic polynomial for a given quaternion polynomial matrix.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

51

Definition 3.2. For A ∈ H[x]m×n , let B = AA∗ and χB be its complex adjoint. Then fB (λ) = det (λI2m − χB ) is called the characteristic polynomial of A. Remark 3.3. By Lemma 3.1, λ can be assumed to be a real indeterminate that enjoys the following: λ = λ and λ commutes element-wise with H [x]. and B = AA∗ . Then fB (λ) = g (λ) where g (λ) ∈

m×n

2

Theorem 3.4. Let A ∈ H [x] (R [x]) [λ].

Proof. We first show that fB (λ) ∈ (R [x]) [λ]. Note that B = AA∗ , we have    T ∗ det (λI2m − χB ) = det (λI2m − χB ) = det (λI2m − χB ) , and thus det (λI2m − χB ) = det(λI2m − χB ) = det (λI2m − χB ). Therefore det (λI2m − χB ) = fB (λ) ∈ (R [x]) [λ] .

(1)

2

Next, we show that fB (λ) = g (λ) where g (λ) ∈ (C [x]) [λ]. Let B = P + Qj. For any fixed 1 ≤ i, j ≤ m, we have Bij = a + bi + cj + dk where a, b, c and d ∈ R [x]. Since B is hermitian, Bji = a − bi − cj − dk and therefore Pij = a + bi, Pji = a − bi and Qij = c + di, Qji = −c − di. So P T = P and Q = −QT . Therefore, χB =

P Q −Q P

=

P Q −Q P T

=⇒ λI2m − χB =

λIm − P −Q

Q λIm − P T

Next, we have:

Im 0



=

−Im Im



Q λIm − P

Im Im



−Im Im

Im 0

T P − λIm . Q 0 Im



λIm − P −Q

Q λIm − P T

Therefore, fB (λ) = det

λIm − P −Q

Q λIm − P T

= det

Q λIm − P

P T − λIm Q

Note that

Q λIm − P

P T − λIm Q

T

=−

Q λIm − P

P T − λIm Q

,

.

.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

52

  Q P T − λIm which implies that λI − P is skew-symmetric. By [16], the determinant of Q   m T Q P − λIm , also called its Pfaffian, can be written as the square of a polynoλI − P Q m

2

mial in its entries. Therefore, fB (λ) = g (λ) where g (λ) ∈ (C [x]) [λ]. Finally we show that g (λ) ∈ (R [x]) [λ]. Suppose otherwise. Then g (λ) = a (λ)+b (λ) i 2 where a)λ) and b(λ) ∈ (R [x]) [λ] with b (λ) = 0. By (1), g (λ) = a(λ)2 − b(λ)2 + 2 2 2 2a(λ)b(λ)i ∈ (R [x]) [λ]. Hence a (λ) = 0 and fB (λ) = g (λ) = (b (λ) i) = −b (λ) where b (λ) ∈ (R [x]) [λ]. For a fixed x ∈ R, let λ ∈ R be large enough such that λ I2m − χB ∈ C2m×2m is diagonally dominant with non-negative diagonal entries and that (b (x)) (λ ) = 0. Since λ I2m −χB is also hermitian, λ I2m −χB is positive definite [9]. 2 But det (λ I2m − χB ) = − ((b (x)) (λ )) < 0, a contradiction. Therefore, b = 0 and thus 2 fB (λ) = g (λ) where g (λ) ∈ (R [x]) [λ]. 2 Corollary 3.5. Let A ∈ H [x] , B = AA∗ and fB (λ) = g (λ) . Then g (B) = 0. We will call g (λ) the generalized characteristic polynomial of A. m×n

2

Proof. Note that g (λ) ∈ (R[x]) [λ] by Theorem 3.4. Then χg(B) = g (χB ). Next, fB (χB ) = 0 by the Cayley–Hamilton theorem for complex polynomial matrices [9]. Therefore g (χB ) = 0, and 0 = g (χB ) = χg(B) , that is, g (B) = 0. 2 From the definition, it is easy to check the following lemma, which have the analogues in the complex case. Lemma 3.6. Let A ∈ H[x]m×n have the Moore–Penrose inverse A† . Set B = AA∗ . Then (i) (ii) (iii)

†

B † = (A∗ ) A† and B † B = AA† . B † B = BB † and (B † B)2 = B † B.  † k  k † B = B and (B n−k )† B n−k = B † B, for any k ∈ N.

The following result is well-known for quaternion case and it is easy to check that the result also holds for quaternion polynomials. Lemma 3.7. Let A ∈ H [x] , B ∈ H [x] and C ∈ H [x] . If A† and B † both exist, n×p then the quaternion polynomial matrix equation AXB = C has a solution in H [x] if † † and only if AA CB B = C, in which case the general solution is m×n

p×q

m×q

X = A† CB † + Y − A† AY BB † , n×p

where Y ∈ H [x]

is arbitrary.

Theorem 3.8. Let A ∈ H [x] have the Moore–Penrose inverse A† and B = AA∗ . Suppose the generalized characteristic polynomial of A is m×n

g (λ) = λm + a1 λm−1 + · · · + ak λm−k + · · · + am−1 λ + am ,

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

53

where ai ∈ R [x]. If k is the largest integer such that ak = 0, then the generalized inverse of A is given by A† = −

 1 ∗  k−1 A B + a1 B k−2 + · · · + ak−1 I . ak

If ai = 0 for all 1 ≤ i ≤ m, then A† = 0. Proof. The proof is similar to the complex case in [6] by using Corollary 3.5, Lemmas 3.6, 3.7 and 2.2. 2 From the above theorem, we can find the Moore–Penrose inverse A† of A by computing its generalized characteristic polynomials. Faddeev [7] modified Leverrier’s method and gave one algorithm to compute {ai } without computing g(λ). Next we extend this algorithm to quaternion polynomial matrices. m×n

Lemma 3.9. Let A ∈ H [x] Then for 1 ≤ k ≤ m, tr



have the Moore–Penrose inverse A† and set B = AA∗ .

B k + a1 B k−1 + · · · + ak−1 B



= −kak ,

where the ai arise from the generalized characteristic polynomial of A: g (λ) = λm + a1 λm−1 + · · · + ak λm−k + · · · + am−1 λ + am . Proof. Let Y = yI where y ∈ R. We can write g (Y ) as: g (Y ) = g(Y ) − g(B)    = (Y − B) Y m−1 + (B + a1 I) Y m−2 + · · · + B m−1 + a1 B m−2 + · · · + am I . As long as y is not an eigenvalue of B, (yI − B) = Y − B is non-singular, so we can write: (Y − B)

−1

  g (Y ) = Y m−1 + (B + a1 I) Y m−2 + B 2 + a1 B + a2 I Y m−3 + · · ·   + B m−1 + a1 B m−2 + · · · + am I .

Taking the traces gives:   −1 tr (Y − B) g (Y ) = my m−1 + tr [(B + a1 I)] y m−2 + · · ·   + tr B m−1 + a1 B m−2 + · · · + am I . Let C = (Y − B) fore,

−1

g (Y ). Since g (Y ) = g (yI) = g (y) I, C = g (y) (Y − B)

−1

. There-

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

54

  −1 trC = g (y) tr (Y − B) .   −1 Let λ1 , · · · , λm , where m ≤ m, be all the non-zero eigenvalues of B. tr (Y − B) −1

is the sum of the eigenvalue of (Y − B) 1 fractions y−λ , · · · , y−λ1  . 1 m

−1

Let ζ be an eigenvalue of (Y − B)

. We will show that these eigenvalues are the

with corresponding eigenvector Z such that: −1

(Y − B)

Z = Zζ.

ζ is real by Lemma 3.1, and hence (Y − B) Z = Z Therefore, y −

1 ζ

= λi ⇒ ζ =

1 y−λi



1 1 =⇒ BZ = Z y − . ζ ζ

for some 1 ≤ i ≤ m .

Since g (y) = (y − λ1 ) (y − λ2 ) · · · (y − λm ), we have that g  (y) = g (y)  1 and trC = g  (y). The derivative of g is also equal to: y−λ 



1 y−λ1

+ ···+

m

g  (y) = my m−1 + a1 (m − 1) y m−2 + · · · + am−1 . Therefore, my m−1 + a1 (m − 1) y m−2 + · · · + am−1

  = my m−1 + tr (B + a1 I) y m−2 + · · · + tr B m−1 + a1 B m−2 + · · · + am I .

Comparing the coefficient of y m−k−1 on both sides, we obtain   ak (m − k) = tr B k + a1 B k−1 + · · · + ak−1 B + ak I   = tr B k + a1 B k−1 + · · · + ak−1 B + tr(ak I), and then   −kak = tr B k + a1 B k−1 + · · · + ak−1 B .

2

Next we give the Leverrier–Faddeev algorithm for finding Moore–Penrose inverses of quaternion polynomial matrices. have the generalized inverse A† and B = AA∗ . Proposition 3.10. Let A ∈ H [x] Suppose the generalized characteristic polynomial of A is m×n

g (λ) = λm + a1 λm−1 + · · · + ak λm−k + · · · + am−1 λ + am ,

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

55

where ai ∈ R [x]. Define a0 = 1. If p is the largest integer such that ap = 0 and we construct the sequence A0 , · · · , Ap as follows:

Ap−1

=0 −1 trA1 = AA∗ B0 1 .. . trAp−1 = AA∗ Bp−2 p−1

= q0 B0 = q1 B1 .. . = qp−1 Bp−1

=I = A1 − q1 I .. . = Ap−1 − qp−1 I

Ap

= AA∗ Bp−1

= qp

= A p − qp I

A0 A1

trAp p

Bp

then qi (x) = −ai (x), i = 0, · · · , p. Proof. We will show qi (x) = −ai (x) by mathematical induction. By the definition, clearly q0 = −a0 holds. Now we assume that qi (x) = −ai (x) holds for all 0 ≤ i ≤ k − 1. Then Ak = AA∗ Bk−1 = BBk−1 = B (Ak−1 − qk−1 I) = B (B (Ak−2 − qk−2 I) − qk−1 I) ··· = B k − q1 B k−1 − q2 B k−2 − · · · − qk−1 B = B k + a1 B k−1 + a2 B k−2 + · · · + ak−1 B, and thus   tr(Ak ) = tr B k + a1 B k−1 + · · · + ak−1 B , which, by Lemma 3.9, is equal to −kak . So qk = for all p ≥ i ≥ 0. 2

trAk k

= −ak . Therefore, qi (x) = −ai (x)

Now combining Theorem 3.8 and Proposition 3.10, we have the following algorithm to compute Moore–Penrose inverses. Algorithm 3.11. Leverrier–Faddeev algorithm for quaternion polynomial matrices Input: A ∈ H[x]m×n Output: The Moore–Penrose inverse A† of A in H[x]n×m if exists. 1. B0 ← Im , a0 ← 1 2. for i = 1, . . . , m do Ai ← AA∗ Bi−1 , ai ← − triAi , Bi ← Ai + ai Im

56

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

3. Find the maximal

index p such that ap = 0. − a1p A∗ Bp−1 , p > 0, 4. Return A† = 0, p = 0. Note that we have to compute many matrix products in Proposition 3.10, which means that Leverrier–Faddeev method is not efficient. In the next section, we will give an efficient way by combining Theorem 3.8 and interpolation methods. 4. Finding Moore–Penrose inverses by interpolation In this section we present an efficient method to obtain the Moore–Penrose inverse of a quaternion polynomial matrix by interpolation at data points of real numbers. Since H is a skew field, the roots of quaternion polynomials in H[x] and interpolations are quite complicated comparing with complex number case. First we recall some basic definitions and properties. An element r ∈ H is a root of a nonzero polynomial f ∈ H [x] iff x − r is a right divisor of f . The set of polynomials in H [x] having r as a root is the left ideal H [x] · (x − r). It is worth mentioning that the evaluations of quaternion polynomials are quite different from commutative case. Let f , g and h ∈ H [x], f = gh and r ∈ H. If h(r) = 0, then f (r) = 0. Otherwise, set β = h (r) = 0. Then the evaluation of f (x) at x = r is defined as   f (r) = g βrβ −1 h (r) .

(2)

In particular, if r is a root of f but not of h, then βrβ −1 is a root of g. We refer the reader to [14] for more details. In [8], Gordon and Motzkin proved that if f ∈ H [x] is of degree n, then the roots of f lie in at most n conjugacy classes of H. It is well-known that Newton’s interpolation and Lagrange’s interpolation play important roles in studying polynomials over fields. Unfortunately one cannot get similar nice formulas in quaternion case. But we still can find a quaternion polynomial from a given set of points. Lemma 4.1. Let c1 , . . . , cn be n pairwise non-conjugate elements of H. Then there is a unique monic polynomial gn ∈ H [x] of degree n such that gn (c1 ) = · · · = gn (cn ) = 0. Moreover, c1 , . . . , cn are the only roots (up to conjugacy classes) of gn in H. Proof. We first show the existence of gn for all n ≥ 1 by mathematical induction. For n = 1, it is trivially true as g = x − c1 . Suppose the claim holds for all 1 ≤ n ≤ k − 1. Let c1 , · · · , ck ∈ H be pairwise non-conjugate. Invoking the inductive hypothesis, there exists a monic polynomial gk−1 of degree k − 1 with c2 , · · · , ck as its only roots (up to conjugacy classes), that is, gk−1 (c1 ) = 0. Construct gk as follows,

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

57

  −1 gk (x) = x − gk−1 (c1 ) c1 gk−1 (c1 ) · gk−1 (x). By Eq. (2), gk (c1 ) = 0. Thus, the claim holds for k. Therefore this claim holds for all n ≥ 1. We next show that gn is unique. For a fixed n, let g  = gn be a monic polynomial of degree n such that g  (c1 ) = · · · = g  (cn ) = 0 too. Then deg (gn − g  ) ≤ n − 1 but gn −g  has roots c1 , · · · , cn which lie in n different conjugacy classes of H, a contradiction. Therefore, gn is unique for all n ≥ 1. 2 Proposition 4.2. Let c1 , · · · , cn+1 ∈ H be pairwise non-conjugate and let d1 , · · · , dn+1 ∈ H. Then there exists a unique lowest degree polynomial f ∈ H [x], of degree p ≤ n, such that f (ci ) = di for all 1 ≤ i ≤ n + 1. Proof. For any 1 ≤ s ≤ n + 1, let S = {1, · · · , n + 1} \ {s }. By Lemma 4.1, we can find a unique monic hS ∈ H [x] of degree n such that hS (cs ) = 0, s ∈ S and that {cs | s ∈ S} and thus we are the only roots (up to conjugacy class) of hS in H. Then hS (cs ) = 0,  0 α ∈ S, can construct a quaternion polynomial gS of degree n such that gS (cα ) = 1 α = s , as follows: −1

gS (x) = hS (cs )

hS (x).

Furthermore we construct a quaternion polynomial f of degree at most n such that f (ci ) = di for all 1 ≤ i ≤ n + 1 as follows:

f=

n+1 

d s  gS .

s =1

Finally we show that f is unique. Suppose we have f  ∈ H [x] of degree p ≤ n such that f  = f and that f  (ci ) = di for all 1 ≤ i ≤ n + 1 too. Then f − f  = 0 is of degree at most ≤ n. But f − f  has roots c1 , · · · , cn+1 which lie in n + 1 conjugacy classes of H, a contradiction. Therefore, f is unique. 2 Next we consider the interpolation for quaternion polynomial matrices. Recall that m×n the degree of a given A ∈ H [x] is equal to deg A = max {deg (Aij ) | 1 ≤ i ≤ m , 1 ≤ j ≤ n}. Then the upper bound of the degree of its Moore–Penrose inverse (if it exists) A† is given by the following lemma. m×n

Lemma 4.3. Let A ∈ H [x]

have the Moore–Penrose inverse A† . Then deg A† ≤ (2m − 1) deg A.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

58

Proof. By Theorem 3.8, deg A† (2m − 1) deg A. 2

≤

  m−2 deg A∗ (AA∗ )

≤

  deg A2m−1

≤

For A = (Aij ) ∈ H[x] and c ∈ H, the evaluation of A at c can be defined as entry-wise in a common sense, that is, A(c) = (Aij (c)). One has to pay an attention that the evaluations of quaternion polynomials have some special rules as we explained at the beginning of this section. Proposition 4.4. Let c1 , · · · , ck+1 ∈ H be pairwise non-conjugate and let A1 , · · · , Ak+1 ∈ n×m Hn×m . Then there is a unique lowest degree matrix A ∈ H [x] of degree p ≤ k, such that A (ci ) = Ai for all 1 ≤ i ≤ k + 1. Proof. For any 1 ≤ n1 ≤ n and 1 ≤ m1 ≤ m, by Proposition 4.2, there is a lowest degree polynomial An1 m1 (x) determined by the values c1 , · · · , ck+1 and (A1 )n1 m1 , · · · , (Ak+1 )n1 m1 . In fact, for any 1 ≤ s ≤ k + 1, let S = {1, · · · , k + 1} \ {s }. Then An1 m1 (x) =

k+1 

(As )n1 m1 gS (x),

s =1



0 α ∈ S, Since n1 and m1 are chosen randomly, a lowest degree ma1 α = s . k+1 trix A that satisfies A (ci ) = Ai for all 1 ≤ i ≤ k +1 is determined by A = ( s =1 As gS ). Next we show that A is unique. Suppose A = A of degree p ≤ p also satisfies A (ci ) = Ai for all 1 ≤ i ≤ k + 1. Then for some 1 ≤ n2 ≤ n and 1 ≤ m2 ≤ m, (A − A )n2 m2 = 0. But (A − A )n2 m2 , of degree at most p ≤ k, has roots c1 , · · · , ck+1 which lie in k + 1 conjugacy classes of H, a contradiction. Therefore, A is unique. 2 where gS (cα ) =

have the Moore–Penrose inverse A† , and set B = AA∗ . Let p be Let A ∈ H [x] the largest integer such that ap = 0. We can construct the sequence A0 , · · · , Ap as in Proposition 3.10. The next theorem gives the interpolation version of Leverrier–Faddeev algorithm. m×n

Theorem 4.5. In the above setting, let k = (2m − 1) deg A and c1 , · · · , ck+1 ∈ R be k + 1 distinct real numbers such that qp (cs ) = 0 for any 1 ≤ s ≤ k + 1. Let S = {1, · · · , k + 1} \ {s }. Then A† =

k+1 

†

A (cs ) gS

s =1

where †

A (cs ) =

  1 ∗ p−1 p−2 A (cs ) B (cs ) − q1 (cs ) B (cs ) − · · · − qp−1 (cs ) I qp (cs )

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

59

and  gS (cα ) =

0 1

α ∈ S, α = s .

Proof. It follows from Theorem 3.8, Proposition 3.10 and Proposition 4.4. 2 5. Implementation and examples The calculations of quaternions are very complicated and time-consuming. It is almost impossible to do some calculations for quaternion polynomials and quaternion polynomial matrices with a little big size by hand. There are only few quaternion packages in computer algebra system Maple. But no one has commends for quaternion polynomials and quaternion polynomial matrices. In [10], we developed a Maple package which includes all basic operations for quaternion polynomials and quaternion polynomial matrices. In particular, we implemented all algorithms in this paper. Next we will give some examples by using our Maple package. Example 5.1. Find, by interpolation, the generalized inverse of the quaternion polynomial matrix ⎛

14x + 14 + 76i + 70j + 56k

⎜ −2x − 2 − 43i − 10j − 8k A=⎝ −3x − 3 + 3i − 15j − 12k −4x − 4 + 4i − 20j − 16k

56 − 28i − 70j + 70k −8 + 4i + 10j − 10k −12 + 6i + 15j − 15k −16 + 8i + 20j − 20k

28j − 56k −4j + 8k −6j + 12k −8j + 16k

⎞

14x − 56 − 8i − 14j − 56k −2x + 8 − 31i + 2j + 8k ⎟ 4×3 ∈H [x] −3x + 12 + 21i + 3j + 12k ⎠ −4x + 16 + 28i + 4j + 16k

From Lemma 4.3, we know that the upper bound of the degree of A† is less than (2m − 1) deg A = (2 × 4 − 1) · 1 = 7. In practice, we don’t need to start from the upper bound. For this example, we may guess deg A† = 2, and choose c1 = 0 and c2 = 1. Then using our Maple package, it is easy to do the following calculations: ⎛

14 + 76i + 70j + 56k ⎜ −2 − 43i − 10j − 8k A (c1 ) = ⎝ −3 + 3i − 15j − 12k −4 + 4i − 20j − 16k

56 − 28i − 70j + 70k −8 + 4i + 10j − 10k −12 + 6i + 15j − 15k −16 + 8i + 20j − 20k

28j − 56k −4j + 8k −6j + 12k −8j + 16k

⎞ −56 − 8i − 14j − 56k 8 − 31i + 2j + 8k ⎟ 12 + 21i + 3j + 12k ⎠ 16 + 28i + 4j + 16k

56 − 28i − 70j + 70k −8 + 4i + 10j − 10k −12 + 6i + 15j − 15k −16 + 8i + 20j − 20k

28j − 56k −4j + 8k −6j + 12k −8j + 16k

⎞ −42 − 8i − 14j − 56k 6 − 31i + 2j + 8k ⎟ . 9 + 21i + 3j + 12k ⎠ 12 + 28i + 4j + 16k

and ⎛

28 + 76i + 70j + 56k ⎜ −4 − 43i − 10j − 8k A (c2 ) = ⎝ −6 + 3i − 15j − 12k −8 + 4i − 20j − 16k

By the algorithm stated in Theorem 4.5, we calculate and obtain: †

†

A (c1 ) = A (0) =



1 × 230175

140−560i−228j−342k 276+88i+426j−382k 32+16i−176j+292k −140−122i+228j+342k

355+1730i−96j+81k −255−870i+126j+54k −340−1160i+168j+72k 282+416i−93j−149k −252−276i−72j+204k −336−368i−96j+272k −176−88i+68j+194k 96+48i+12j−204k 128+64i+16j−272k −355+2021i+96j−81k 255−1176i−126j−54k 340−1568i−168j−72k



L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

60

and †

†

A (c2 ) = A (1) =



1 × 230175 289+1675i−8j+15k −219−840i+78j+90k −292−1120i+104j+120k 326+328i+17j−39k −276−228i−132j+144k −368−304i−176j+192k −176−88i−20j+150k 96+48i+60j−180k 128+64i+80j−240k −289+2076i+8j−15k 219−1206i−78j−90k 292−1608i−104j−120k

152−550i−244j−330k 268+104i+406j−402k 32+16i−160j+300k −152−132i+244j+330k

.

By Theorem 4.5, we have: †

A =

2 

†

†

†

A (cs ) gS = A (0) (1 − x) + A (1) x

s =1

=



1 × 230175 (12 + 10i − 16j + 12k) x + 140 − 560i − 228j − 342k (−66 − 55i + 88j − 66k) x + 355 + 1730i − 96j + 81k (−8 + 16i − 20j − 20k) x + 276 + 88i + 426j − 382k (44 − 88i + 110j + 110k) x + 282 + 416i − 93j − 149k (16j + 8k) x + 32 + 16i − 176j + 292k (−88j − 44k) x − 176 − 88i + 68j + 194k (−12 − 10i + 16j − 12k) x − 140 − 122i + 228j − 342k (66 + 55i − 88j + 66k) x − 355 + 2021i + 96j − 81k

(36 + 30i − 48j + 36k) x − 255 − 870i + 126j + 54k (48 + 40i − 64j + 48k) x − 340 − 1160i + 168j + 72k (−24 + 48i − 60j − 60k) x − 252 − 276i − 72j + 204k (−32 + 64i − 80j + 80k) x − 366 − 368i − 96j + 272k (48j + 24k) x + 96 + 48i + 12j − 204k (64j + 32k) x + 128 + 64i + 16j − 272k (−36 − 30i + 48j − 36k) x + 255 − 1176i − 126j − 54k (−48 − 40i + 64j − 48k) x + 340 − 1568i − 168j − 72k

.

Direct calculation shows that A† satisfies the four defining relations of the Moore– Penrose inverse. Therefore A† is the Moore–Penrose inverse of A. Example 5.2. As special cases, the algorithms in this paper and our Maple package can compute the Moore–Penrose inverses for (real) complex polynomial matrices and quaternion matrices.

1 i + 2k 3 . Then its Moore–Penrose inverse can be For instance, let A = i 6 + j 7 2×3 found by using our Maple package as follows: ⎛

47 347

⎜ 63 A† = ⎝ − 347 −

+

21 694 i

+

11 694 j

28 21 101 347 i + 694 j − 694 k 57 49 77 347 + 694 i + 694 k

⎞ 11 11 347 i − 694 k ⎟ 61 21 6 21 694 + 694 i − 347 j + 347 k ⎠ 21 21 33 3×2 347 − 694 i − 694 k 21 − 694 −

.

Clearly these examples also show that it could be very time-consuming to find the Moore– Penrose inverses by hand. Acknowledgments We would like to thank two anonymous referees for their valuable comments. This research was supported by the grants from the National Natural Science Foundation of China (11171205), the Natural Science Foundation of Shanghai (11ZR1412500), Shanghai Leading Academic Discipline Project (J50101), the National Sciences and Engineering Research Council of Canada (NSERC) and URGP from the University of Manitoba.

L. Huang et al. / Linear Algebra and its Applications 475 (2015) 45–61

61

References [1] A. Baker, Right eigenvalues for quaternionic matrices: a topological approach, Linear Algebra Appl. 286 (1–3) (1999) 303–309. [2] S. Barnett, Leverrier’s algorithm: a new proof and extensions, SIAM J. Matrix Anal. Appl. 10 (1989) 551–556. [3] A. Ben-Israel, T.N.E. Greville, Generalized Inverses: Theory and Applications, second ed., SpringerVerlag, New York, 2003. [4] J.L. Brenner, Matrices of quaternions, Pacific J. Math. 1 (1951) 329–335. [5] A. Damiano, G. Gentili, D. Sreuppa, Computations in the ring of quaternionic polynomials, J. Symbolic Comput. 45 (2010) 38–45. [6] H.P. Decell, An application of the Cayley–Hamilton theorem to generalized matrix inversion, SIAM Rev. 7 (1965) 526–528. [7] V.N. Faddeeva, Computational Methods of Linear Algebra, Dover Books on Advanced Mathematics, Dover Publications, 1959. [8] B. Gordon, T.S. Motzkin, On the zeros of polynomials over division rings, Trans. Amer. Math. Soc. 116 (1965) 218–226. [9] R.A. Horn, C.R. Johnson, Matrix Analysis, Cambridge University Press, 1990. [10] Liji Huang, Yang Zhang, Maple Package for Quaternion Matrices, The University of Manitoba, 2013. [11] J. Jones, N.P. Karampetakis, A.C. Pugh, The computation and application of the generalized inverse via Maple, J. Symbolic Comput. 25 (1998) 99–124. [12] I.I. Kyrchei, Determinantal representation of the Moore–Penrose inverse matrix over the quaternion skew field, J. Math. Sci. 180 (1) (2012) 23–33. [13] I.I. Kyrchei, Cramers rule for some quaternion matrix equations, Appl. Math. Comput. 217 (2010) 2024–2030. [14] T.Y. Lam, A First Course in Noncommutative Rings, Graduate Texts in Mathematics, vol. 131, Springer, 2001. [15] H.C. Lee, Eigenvalues and canonical forms of matrices with quaternion coefficients, Proc. R. Ir. Acad. A Math. Phys. Sci. 52 (1949) 253–260. [16] T. Muir, A Treatise on the Theory of Determinants, Dover Phoenix Editions, Dover Publications, 2003. [17] I. Niven, Equations in quaternions, Amer. Math. Monthly 48 (1941) 654–661. [18] I. Niven, The roots of a quaternion, Amer. Math. Monthly 49 (1942) 386–388. [19] G. Opfer, Polynomials and Vandermonde matrices over the field of quaternions, Electron. Trans. Numer. Anal. 36 (2009) 9–16. [20] R. Pereira, P. Rocha, On the determinant of quaternionic polynomial matrices and its application to system stability, Math. Methods Appl. Sci. 31 (2008) 99–122. [21] R. Pereira, P. Rocha, P. Vettori, Algebraic tools for the study of quaternionic behavioral systems, Linear Algebra Appl. 400 (2005) 121–140. [22] R. Pereira, P. Vettori, Stability of quaternionic linear systems, IEEE Trans. Automat. Control 51 (2006) 518–523. [23] R. Penrose, A generalized inverse for matrices, Math. Proc. Cambridge Philos. Soc. 51 (1955) 406–413. [24] R. Puystjens, Moore–Penrose inverses for matrices over some Noetherian rings, J. Pure Appl. Algebra 31 (1–3) (1984) 191–198. [25] P.S. Stanimirovic, M.D. Petkovic, Computing generalized inverse of polynomial matrices by interpolation, Appl. Math. Comput. 172 (2006) 508–523. [26] F. Zhang, Quaternions and matrices of quaternions, Linear Algebra Appl. 251 (0) (1997) 21–57. [27] M. Ziegler, Quasi-optimal arithmetic for quaternion polynomials, ISAAC 2003, pp. 705–715.