A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems Ru-Ru Ma a , Zhi-Gang Jia b, *, Zheng-Jian Bai a a b

School of Mathematical Sciences, Xiamen University, Xiamen 361005, PR China School of Mathematics and Statistics, Jiangsu Normal University, Xuzhou 221116, PR China

article

info

Article history: Received 26 May 2017 Received in revised form 4 September 2017 Accepted 14 October 2017 Available online xxxx Keywords: Real counterpart Quaternion Hermitian eigenvalue problem Jacobi rotation Structure-preserving algorithm

a b s t r a c t A new real structure-preserving Jacobi algorithm is proposed for solving the eigenvalue problem of quaternion Hermitian matrix. By employing the generalized JRS-symplectic Jacobi rotations, the new quaternion Jacobi algorithm can preserve the symmetry and JRS-symmetry of the real counterpart of quaternion Hermitian matrix. Moreover, the proposed algorithm only includes real operations without dimension-expanding and is generally superior to the state-of-the-art algorithm. Numerical experiments are reported to indicate its efficiency and accuracy. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Quaternions, proposed by William Rowan Hamilton in 1843 [1], and quaternion matrices [2] have extensive applications in many research fields, consisting of quantum mechanics [3], field theory [4] and image processing [5]. Among the topics about quaternion matrices, the Hermitian eigenvalue problem stands on the indispensable position. High-performance algorithms for solving such problem are still greatly in need. In this paper, we present an efficient and stable Jacobi algorithm for calculating the eigenvalues and corresponding eigenvectors of a quaternion Hermitian matrix, with only using real arithmetic and without dimension expanding. Jacobi’s method is one of the classical methods to compute full eigenvalues and eigenvectors of a real symmetric matrix, which was first put forward by Jacobi in 1846 [6]. Utilizing a series of Jacobi rotations, any real symmetric matrix can be reduced to a diagonal form. The whole processing can be decomposed into sweeps of Jacobi rotations on 2-by-2 principle submatrices, based on the fact that any 2-by-2 symmetric matrix can be diagonalized through a rotation matrix which is denotative in terms of A by means of a simple formula [7]. Demmel et al. [8] elucidated that Jacobi’s method is more accurate than QR method. Later, Mathias [9] strengthened this result. Not only the preferable accuracy of Jacobi method but also the parallelism [10,11] gives researchers enthusiasm to further study Jacobi method. To solve the quaternion Hermitian eigenvalue problem, one may extend classical algorithms, which are designed for real or complex matrices, to quaternion Hermitian matrices. However, we have to overcome difficulties caused by the non-commutativity of quaternions. Bunse-Gerstner, Byers and Mehrmann [12] firstly proposed the quaternion QR Algorithm, which reduces a general quaternion matrix to a Schur-like triangular form by a sequence of implicit QR steps. Bihan and Sangwine [13] raised a quaternion Jacobi algorithm for quaternion Hermitian matrices, in which a quaternion Jacobi rotation performs on a 2-by-2 Hermitian matrix straightly and obtains a real 2-by-2 diagonal matrix. These two algorithms use complex operations, which are more tiresome than real operations. The eigendata of a quaternion matrix

*

Corresponding author. E-mail addresses: [email protected] (R.-R. Ma), [email protected] (Z.-G. Jia), [email protected] (Z.-J. Bai).

https://doi.org/10.1016/j.camwa.2017.10.009 0898-1221/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

2

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

can be generated by the eigendata of its real counterpart. Instead of the quaternion Hermitian matrix, this paper focuses on its real counterpart without dimension-expanding. The real counterpart of an n × n quaternion Hermitian matrix has double structures: symmetry and JRS-symmetry. To improve the efficiency of an algorithm, these two inherent structures should be preserved. Jia, Wei and Lin [14] designed a real structure-preserving method for reducing a quaternion Hermitian matrix to the real tridiagonal form. Li et al. [15] generalized the real Householder-based transformation in [14] and proposed other three Householder-based transformations with preserving structures of the real counterpart of quaternion matrix. Sparked by these results, we propose a structure-preserving Jacobi algorithm for solving the quaternion Hermitian eigenvalue problem, which directly reduces the real counterpart to a diagonal form. To this end, we firstly construct a generalized JRS-symplectic Jacobi rotation for any quaternion Hermitian matrix. Secondly, we show that the real counterpart of a 2-by-2 quaternion Hermitian matrix can be reduced to a diagonal form by using a generalized JRS-symplectic Jacobi rotation. Thirdly, we show that the real counterpart of an n-by-n quaternion Hermitian matrix can be diagonalized by performing a series of generalized JRSsymplectic Jacobi rotations. Finally, we propose a structure-preserving Jacobi algorithm for solving the quaternion Hermitian eigenvalue problem, where all operations are in real arithmetic and the dimension-expanding problem, caused by the real counterpart method, is avoided by taking full use of the symmetry and JRS-symmetry of the real counterpart. The rest of this paper is organized as follows. In Section 2, we present preliminary results. In Section 3, we propose a structure-preserving Jacobi algorithm for solving the quaternion Hermitian eigenvalue problem. In Section 4, we report some numerical experiments to indicate the effectiveness of the proposed algorithm. Finally, we give some concluding remarks in Section 5. 2. Preliminaries In this section, we briefly review some basic properties of quaternions and the real counterpart of a quaternion matrix. Interested reader may find more information in [1,2,16] and references therein. We need the following notations. R and H denote the real field and the quaternion field, respectively. Fm×n denotes the set of all m × n matrices on F. For any matrix A ∈ Fm×n , AT , A¯ and AH refer to the transpose, conjugate and conjugate transpose of A accordingly. Let In and On be the n × n unit matrix and zero matrix separately. ∥ · ∥F denotes the Frobenius matrix norm. A quaternion a ∈ H is defined by a = a0 + a1 i + a2 j + a3 k, where a0 , a1 , a2 , a3 ∈ R and the imaginary units are determined by the rules i2 = j2 = k2 = ijk = −1, which imply jk = −kj = i, ki = −ik = j, ij = −ji = k. The conjugate of a ∈ H is given by a¯ = a0 − a1 i − a2 j − a3 k and the module |a| is defined by

|a|2 = a20 + a21 + a22 + a23 = aa¯ = a¯ a. Thus the multiplicative inverse of any nonzero quaternion a is given by a−1 = a¯ /|a|2 . The quaternion field H is an associative but noncommutative division algebra over R with the standard basis {1, i, j, k}. A quaternion matrix A ∈ Hn×n is Hermitian if A = AH . For a quaternion matrix A = A0 + A1 i + A2 j + A3 k ∈ Hn×n where A0 , A1 , A2 , A3 ∈ Rn×n , the corresponding real counterpart is defined by

⎡

⎤

A0 A2 A1 A3 ⎢−A2 A0 A3 −A1 ⎥ ΓA = ⎣ ∈ R4n×4n . −A1 −A3 A0 A2 ⎦ −A3 A1 −A2 A0

(2.1)

Now we introduce the structural properties of ΓA . Define three orthogonal matrices:

⎡

0 ⎢0 Jn = ⎣ In 0

⎤

⎡

0 −In 0 0 −In 0 0 0 − In ⎥ ⎢In 0 0 ,R = ⎣ 0 0 0 ⎦ n 0 0 0 In 0 0 0 0 − In

⎤

⎡

0 0 0 0⎥ ⎢0 0 ,S = ⎣ In ⎦ n 0 − In 0 In 0

⎤

0 − In In 0 ⎥ . 0 0 ⎦ 0 0

A matrix Ω ∈ R4n×4n is called JRS-symmetric if Jn Ω JnT = Ω , Rn Ω RTn = Ω and Sn Ω SnT = Ω . A matrix Ω ∈ R4n×4n is called JRS-symplectic if Ω Jn Ω T = Jn , Ω Rn Ω T = Rn and Ω Sn Ω T = Sn . A matrix Ω ∈ R4n×4n is called orthogonal JRS-symplectic if it is orthogonal and JRS-symplectic. The quaternion matrix A is Hermitian if and only if A0 is symmetric, A1 , A2 , A3 are skew-symmetric. The real counterpart ΓA of Hermitian quaternion matrix A is symmetric and JRS-symmetric. We now recall some results on the relationship between a quaternion matrix and its real counterpart. Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

3

Lemma 2.1 ([14]). Let A = X + Yj be a quaternion matrix, where X = A0 + A1 i and Y = A2 + A3 i with A0 , A1 , A2 , A3 ∈ Rn×n . Then there exists a unitary quaternion matrix

⎡

⎤ −jIn −iIn −kIn jIn −iIn kIn ⎥ −jIn iIn kIn ⎦ jIn iIn −kIn

In 1 ⎢In U = ⎣ 2 In In such that

X + Yj 0 0 0 X − Yj 0 H ⎢ ΓA = U ⎣ 0 0 X¯ + Y¯ j 0 0 0 X¯

⎡

0 0 0

⎤ ⎥ ⎦ U,

(2.2)

− Y¯ j

where ΓA is the real counterpart of A defined as in (2.1). From Lemma 2.1 we know that the eigenvalue problem of a quaternion matrix is equivalent to the eigenvalue problem of its real counterpart. Lemma 2.2 ([14]). Let Q and H be two n × n quaternion matrices, α ∈ R. Then we have the following results. (1) (2) (3) (4)

ΓQ +H = ΓQ + ΓH ; ΓαQ = α ΓQ ; ΓQH = ΓQ ΓH . ΓQ is JRS-symmetric. Q is unitary if ΓQ is orthogonal. If ΓQ is orthogonal then it must be orthogonal JRS-symplectic.

Lemma 2.3 ([17]). Let Q ∈ Hn×n . Then Q is a diagonalizable quaternion matrix if and only if its real counterpart ΓQ is a diagonalizable matrix. Recall that λ ∈ C is an eigenvalue of a quaternion matrix A = A0 + A1 i + A2 j + A3 k ∈ Hn×n with nonzero eigenvector x ∈ Hn if Ax = xλ. The eigenvalues of A are precisely those eigenvalues of ΓA , and consequently the real eigenvalues appear in fours and the purely imaginary eigenvalues appear in pairs and in conjugate pairs. Our main work is calculating the eigenvalues of the real counterpart of a quaternion Hermitian matrix by employing a structured-preserving algorithm. 3. Structure-preserving Jacobi algorithm In this section, we propose a structure-preserving Jacobi algorithm for solving the quaternion Hermitian eigenvalue problem. The structure-preserving algorithm works by performing a series of orthogonal JRS-symplectic transformations ΓA ← ΓGT ΓA ΓG with the property that each updated ΓA is full but more approximating to a diagonal form than the previous one. Eventually, the off-diagonal entries of the updated ΓA are small enough and the diagonal entries can be seen as the eigenvalues of ΓA numerically. Suppose that A = A0 + A1 i + A2 j + A3 k ∈ Hn×n is Hermitian, where A0 , A1 , A2 , A3 ∈ Rn×n with AT0 = A0 and ATs = −As for s = 1, 2, 3. Let ΓA ∈ R4n×4n be the real counterpart of A as defined by (2.1). The idea behind the structure-preserving Jacobi method is to gradually reduce the quantity 4n 4n 4n ( ) 21 (∑ ) 12 ∑ ∑ off(ΓA ) = ∥ΓA ∥2F − ΓA 2ii = ΓA 2ij

(3.1)

i=1 j=1,j̸ =i

i=1

to zero. To do so, we mainly adopt the following so-called generalized orthogonal JRS-symplectic Jacobi rotation: G0 (p, q; p, q; θ ) G2 (p, q; p, q; θ ) G1 (p, q; p, q; θ ) G3 (p, q; p, q; θ ) ⎢−G2 (p, q; p, q; θ ) G0 (p, q; p, q; θ ) G3 (p, q; p, q; θ ) −G1 (p, q; p, q; θ )⎥ := ⎣ , −G1 (p, q; p, q; θ ) −G3 (p, q; p, q; θ ) G0 (p, q; p, q; θ ) G2 (p, q; p, q; θ )⎦ −G3 (p, q; p, q; θ ) G1 (p, q; p, q; θ ) −G2 (p, q; p, q; θ ) G0 (p, q; p, q; θ )

⎡

ΓG(p,q;p,q;θ )

⎤

where G(p, q; p, q; θ ) = In + ep eq

[

[ ][ T] ] cr − 1 s ep ∈ Hn×n −¯s cr − 1 eTq

(3.2)

is an n × n unitary quaternion matrix. Here, ep is the pth unit vector, cr = cos(θ ) ∈ R, s = s0 + s1 i + s2 j + s3 k ∈ H with cr2 + |s|2 = 1, and G0 (p, q; p, q; θ ) = In + ep eq

[

][ T] [ ] cr − 1 s0 ep , −s0 cr − 1 eTq

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

4

R.-R. Ma et al. / Computers and Mathematics with Applications (

G1 (p, q; p, q; θ ) = ep eq

[

G2 (p, q; p, q; θ ) = ep eq

[

G3 (p, q; p, q; θ ) = ep eq

[

[ ][ ] ] 0 s1 eTp eTq

s1 0

[ ][ ] ] 0 s2 eTp eTq

s2 0

[ ][ ] ] 0 s3 eTp eTq

s3 0

)

–

, , .

Then we can easily verify that ΓG is orthogonal JRS-symplectic, i.e., ΓGT ΓG = I4n , ΓG Jn ΓGT = Jn , ΓG Rn ΓGT = Rn , and

ΓG Sn ΓGT = Sn .

The basic step in our structure-preserving Jacobi algorithm includes: (a) choosing an index pair (p, q) such that 1 ≤ p < q ≤ n; (b) calculating a rotation angle θ such that B(p, q; p, q) ≡

[

bpp bpq bqp bqq

]

[

c s = r −¯s cr

]H [

app apq aqp aqq

][

cr s −¯s cr

]

≡ G(p, q; p, q; θ )H A(p, q; p, q)G(p, q; p, q; θ )

(3.3)

is diagonal in the sense that T ΓB(p,q;p,q) = ΓG(p ,q;p,q;θ ) ΓA(p,q;p,q) ΓG(p,q;p,q;θ )

= B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q), where |bpq | = |bqp | = 0 and cr = cos(θ ); (c) overwriting ΓA with ΓB = ΓGT ΓA ΓG , where G = G(p, q; p, q; θ ) and B = [bij ] ∈ Hn×n . Now we conclude the diagonalization of a 2-by-2 quaternion Hermitian matrix in the following theorem. Theorem 3.1. Suppose

[

a a A(p, q; p, q) = pp pq aqp aqq

]

is a 2-by-2 principal submatrix of quaternion Hermitian matrix A, where app , aqq ∈ R, apq = apq0 + apq1 i + apq2 j + apq3 k ∈ H, and aqp = apq . Then there exists a 2-by-2 unitary quaternion matrix given by

[

cr s G(p, q; p, q; θ ) = −¯s cr

]

with cr and s defined in (3.2) such that T ΓB(p,q;p,q) = ΓG(p ,q;p,q;θ ) ΓA(p,q;p,q) ΓG(p,q;p,q;θ )

= B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q), for some B(p, q; p, q) =

[

]

bpp 0 . 0 bqq

Proof. By the hypothesis, we have T ΓG(p ,q;p,q;θ ) ΓA(p,q;p,q) ΓG(p,q;p,q;θ ) ⎡

cr

s0 ⎢−s0 cr ⎢0 −s 2 ⎢ ⎢−s2 0 ⎢ =⎢ ⎢0 −s1 ⎢−s1 0 ⎣ 0 −s3 −s3 0

0 s2 cr

s2 0 s0 −s0 cr 0 − s3 −s3 0 0 s1 s1 0

0 s1 0 s3 cr

s1 0 s3 0 s0 −s0 cr 0 − s2 − s2 0

0 s3 0

⎤T ⎡

s3 0 ⎥ −s1 ⎥ ⎥ − s1 0 ⎥ ⎥ 0 s2 ⎥ ⎥ s2 0 ⎥ ⎦ cr s0 −s0 cr

app ⎢apq0 ⎢0 ⎢ ⎢apq2 ⎢ ⎢0 ⎢ ⎢apq1 ⎣ 0 apq3

apq0 aqq −apq2 0 −apq1 0 −apq3 0

0

−apq2 app apq0 0 apq3 0 −apq1

apq2 0 apq0 aqq −apq3 0 apq1 0

0

−apq1 0

−apq3 app apq0 0 apq2

apq1 0 apq3 0 apq0 aqq −apq2 0

0

−apq3 0 apq1 0 −apq2 app apq0

apq3 0 ⎥ −apq1 ⎥ ⎥ ⎥ 0 ⎥ apq2 ⎥ ⎥ ⎥ 0 ⎦ apq0 aqq

⎤

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

R.-R. Ma et al. / Computers and Mathematics with Applications (

cr

0 s2 −s2 cr 0 −s0 −s1 0 0 −s3 0 −s3 0 −s3 0 s1

⎡

⎢−s0 ⎢0 ⎢ ⎢−s2 ⎢ ⎢0 ⎢ ⎢−s1 ⎣ b11 ⎢b12 ⎢0 ⎢ ⎢b14 =⎢ ⎢0 ⎢ ⎢b16 ⎣ 0 b18

⎡

s0 cr

b12 b22 b23 0 b25 0 b27 0

0 b23 b33 b34 0 b36 0 b38

s2 0 s0 cr

0 s1 0 s3 −s3 cr 0 − s0 s1 0 0 − s2

b14 0 b34 b44 b45 0 b47 0

0 b25 0 b45 b55 b56 0 b58

b16 0 b36 0 b56 b66 b67 0

s1 0 s3 0 s0 cr

0 s3 0

0 b27 0 b47 0 b67 b77 b78

b18 0 ⎥ b38 ⎥ ⎥ 0 ⎥ ⎥, b58 ⎥ ⎥ 0 ⎥ ⎦ b78 b88

)

–

5

s3 0 ⎥ − s1 ⎥ ⎥ − s1 0 ⎥ ⎥ 0 s2 ⎥ ⎥ s2 0 ⎥ ⎦ −s2 cr s0 0 −s0 cr

⎤

⎤

(3.4)

where b12 = b34 = b56 = b78

b16

b14

= cr s0 (app − aqq ) + (cr2 − s20 + s21 + s22 + s23 )apq0 − 2s0 s1 apq1 − 2s0 s2 apq2 − 2s0 s3 apq3 , = −b25 = b47 = −b38 = cr s1 (app − aqq ) + (cr2 + s20 − s21 + s22 + s23 )apq1 − 2s0 s1 apq0 − 2s1 s2 apq2 − 2s1 s3 apq3 , = −b23 = b58 = −b67 = cr s2 (app − aqq ) + (cr2 + s20 + s21 − s22 + s23 )apq2 − 2s0 s2 apq0 − 2s1 s2 apq1 − 2s3 s2 apq3 ,

b18 = −b27 = b36 = −b45

= cr s3 (app − aqq ) + (cr2 + s20 + s21 + s22 − s23 )apq3 − 2s0 s3 apq0 − 2s1 s3 apq1 − 2s2 s3 apq2 . Let b12 = b34 = b56 = b78 = b16 = −b25 = b47 = −b38 = b14 = −b23 = b58 = −b67 = b18 = −b27 = b36 = −b45 = 0. By simple calculation, we get 1

cr = cos(θ ) = √ s0 = √ s2 = √

1 + t2 apq0

a2pq0

+

a2pq1

+

a2pq1

+

,

a2pq2

+

a2pq3

+

a2pq3

apq2 a2pq0

+

a2pq2

sin θ,

apq1 s1 = √ sin θ, 2 2 apq0 + apq1 + a2pq2 + a2pq3

sin θ,

apq3 s3 = √ sin θ, 2 2 apq0 + apq1 + a2pq2 + a2pq3

where sin θ = tcr ,

t=

sign(τ ) , √ 1 + τ2

|τ | +

τ= √ 2

a2pq0

aqq − app

+ a2pq1 + a2pq2 + a2pq3

.

Thus T ΓG(p ,q;p,q;θ ) ΓA(p,q;p,q) ΓG(p,q;p,q;θ ) = B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q) ⊕ B(p, q; p, q),

where B(p, q; p, q) =

[

bpp 0 0 bqq

]

with bpp = b11 = b33 = b55 = b77 Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

6

R.-R. Ma et al. / Computers and Mathematics with Applications (

bqq

)

–

= cr2 app + (s20 + s21 + s22 + s23 )aqq − 2cr s0 apq0 − 2cr s1 apq1 − 2cr s2 apq2 − 2cr s3 apq3 , = b22 = b44 = b66 = b88 = cr2 aqq + (s20 + s21 + s22 + s23 )app + 2cr s0 apq0 + 2cr s1 apq1 + 2cr s2 apq2 + 2cr s3 apq3 .

This completes the proof. □ We observe from Theorem 3.1 that the updated matrix B agree with A except the pth row (column) and the qth row (column). We also note that the Frobenius norm is unitary invariant. We get by (3.3),

|app |2 + |aqq |2 + 2|apq0 |2 + 2|apq1 |2 + 2|apq2 |2 + 2|apq3 |2 = |bpp |2 + |bqq |2 . Thus

off(B)2 = ∥B∥2F −

n ∑

b2ii

i=1

= ∥A∥2F −

n (∑

a2ii − |app |2 − |aqq |2 + |bpp |2 + |bqq |2

)

i=1

= ∥A∥2F −

n ∑

a2ii + |app |2 + |aqq |2 − |bpp |2 − |bqq |2

(

)

i=1

= off(A)2 − 2|apq0 |2 − 2|apq1 |2 − 2|apq2 |2 − 2|apq3 |2 .

(3.5)

We see that ∥ ∥ = 4∥A∥ Hence, off(ΓA )2 = 4off(A)2 . Our goal is to minimize off(B). The best choice of p, q should be |apq | = maxi̸=j |aij |. In this sense, ΓA moves closer to a diagonal form under generalized orthogonal JRS-symplectic rotations. Based on Theorem 3.1, we present the following algorithm for generating a generalized orthogonal JRS-symplectic Jacobi rotation for a 2-by-2 quaternion Hermitian matrix. This algorithm costs 31 flops.

ΓA 2F

2 F.

Algorithm 3.2. Given a 2-by-2 principal submatrix of quaternion Hermitian matrix A: A(p, q; p, q) =

]

[

app aqp

apq , aqq

where app , aqq ∈ R, apq = apq0 + apq1 i + apq2 j + apq3 k ∈ H, and aqp = apq , this algorithm computes a cosine–sine group T (cr , s0 , s1 , s2 , s3 ), i.e., a 2-by-2 unitary quaternion matrix G(p, q; p, q; θ ), such that ΓG(p ,q;p,q;θ ) ΓA(p,q;p,q) ΓG(p,q;p,q;θ ) is a diagonal matrix. function [cr , s0 , s1 , s2 , s3 ] = RQJR(app , aqq , apq0 , apq1 , apq2 , apq3 ) if apq0 = apq1 = apq2 = apq3 = 0 Γ G = I8 ; else √ a1 = a2pq0 + a2pq1 + a2pq2 + a2pq3 ;

τ = (app − aqq )/(2a1 ); if |τ |= 0 √ t = 1/(|τ |+ 1 + τ 2 ); else √ t = sign(τ )/(|τ |+ 1 + τ 2 ); end if √ a2 = 1/ 1 + t 2 ; a3 = ta2 ; t˜ = a3 /a1 ; cr = a2 ; [s0 s1 s2 s3 ] = t˜[apq0 apq1 apq2 apq3 ];

We now present a structure-preserving Jacobi algorithm for solving the eigenvalue problem for the quaternion Hermitian matrix A = A0 + A1 i + A2 j + A3 k ∈ Hn×n . The basic iterative scheme is given by (k)

A1

⎢ (k) (k) ⎢−A2 A0 ΓA(k) = ⎢ ⎢ (k) (k) ⎣−A1 −A3

A3

⎡

(k)

A0

−A(k) 3

A2

(k)

A1

(k) (k) (k)

A0

−A(k) 2

(k)

A3

⎤

⎥ ⎥ −A(k) 1 ⎥ = ΓG(k) T ΓA(k−1) ΓG(k) , k = 1, 2 . . . , (k) ⎥ A2 ⎦

(3.6)

(k)

A0

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

7

where A(0) = A, ΓG(k) is a generalized orthogonal JRS-symplectic Jacobi rotation obtained by using Algorithm 3.2 to A(k−1) . By Lemma 2.2 and (3.6), we have

ΓA(k) = ΓG(k) T ΓA(k−1) ΓG(k) = Γ(G(k) )H ΓA(k−1) ΓG(k) = Γ(G(k) )H A(k−1) G(k) , which implies that A(k) = (G(k) )H A(k−1) G(k) .

(3.7)

On the convergence of the structure-preserving Jacobi algorithm, we have the following result. Theorem 3.3. Let the eigenvalues of quaternion Hermitian matrix A ∈ Hn×n be λ1 , . . . , λn . Then lim ΓA(k) = diag(λσ1 , . . ., λσn , λσ1 , . . ., λσn , λσ1 , . . ., λσn , λσ1 , . . ., λσn ),

k→∞

where {σ1 , . . . , σn } is a permutation of {1, 2, . . . , n}. Proof. Firstly, we prove that off(ΓA(k) ) tends to zero as k → ∞. From (3.5),

off(ΓA(k) )2 = 4off(A(k) )2 ) ( 2 −1) 2 = 4 off(A(k−1) ) − 2|a(k | , pq where

(k−1) apq

(3.8)

is the largest (in modulus) off-diagonal entry of A 2

(k−1)

. Notice that

2

−1) off(A(k−1) ) ≤ n(n − 1)|a(k | . pq

(3.9)

Substituting (3.9) into (3.8) yields

( ) 2 2 2 off(ΓA(k) )2 ≤ 4 off(A(k−1) ) − off(A(k−1) ) n(n − 1) ( 1) (k−1) 2 ≤ 4 1− off(A ) N ) ( 1 off(ΓA(k−1) )2 = 1− N ( 1 )k ≤ 1− off(ΓA(0) )2 , N

where N = − 1). Since 1 − N1 < 1, limk→∞ off(ΓA(k) ) = 0. Secondly, we show that there exists a permutation σ = {σ1 , . . . , σn } such that 1 n(n 2

(k)

lim aii = λσi , i = 1, 2, . . ., n,

(3.10)

k→∞ (k)

where aii denotes the (i, i)th entry of A(k) . Suppose that the smallest distance between distinct eigenvalues of A is δ , i.e.,

δ = min{|λi − λj | | λi , λj ∈ λ(A), λi ̸= λj }. For any positive number ϵ with ϵ < δ/4, limk→∞ off(ΓA(k) ) = 0 implies that there exists k0 > 0 such that for all k ≥ k0 ,

off(ΓA(k) ) < ϵ <

δ 4

.

Recall that A(k0 ) and A have the same eigenvalues. Define the diagonal matrix ΓD(k0 ) as a diagonal matrix generated by the diagonal elements of ΓA(k0 ) and thus

( ) (k) (k) ΓD(k0 ) = diag a11 , a22 , . . . , a(k) nn . By using the Wielandt–Hoffman theorem (see [7, Theorem 8.1.4]), we have (k )

|λσi − aii 0 | ≤ ∥ΓA(k0 ) − ΓD(k0 ) ∥2 ≤ off(ΓA(k) ) < ϵ < 0

δ 4

, i = 1, 2, . . . , n.

(3.11)

Next we prove that for any k ≥ k0 + 1,

|λσi − a(k) ii | < ϵ, i = 1, 2, . . . , n,

(3.12) (k)

under the assumption that (3.12) holds for k − 1. Since only two diagonal entries of A are different from that of A(k−1) (k) (k) (i.e., app and aqq ), it is sufficient to show that (3.12) holds for i = p, q. By Theorem 3.1, we have 2 (k−1) −1) −1) −1) a(k) − cr sa(k − cr a(k s¯ + |s|2 a(k pp = cr app qp pq qq

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

8

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

−1) −1) −1) −1) −1) = a(k − |s|2 a(k s¯ + |s|2 a(k − cr a(k − cr sa(k qq pq qp pp pp −1) −1) −1) −1) −1) = a(k − t 2 cr2 a(k − cr2 tma(k − cr2 ta(k m + t 2 cr2 a(k pp pp qp pq qq −1) −1) −1) −1) −1) − ta(k − tma(k − cr2 (t 2 a(k = a(k ) m + t 2 a(k pq qp pp pp qq −1) −1) −1) −1) −1) − ta(k ) − tma(k − a(k − cr2 (t 2 (a(k = a(k m) pq qp pp qq pp −1) −1) −1) −1) −1) = a(k m − t 3 ma(k − tma(k m) − cr2 (ta(k − ta(k pq pp pq qp qp −1) −1) − cr2 t(1 + t 2 )a(k = a(k qp pp −1) −1) − ta(k = a(k qp pp

(3.13)

−1) (k−1) , + ta(k a(k) pq qq = aqq

(3.14)

and

where m = m0 + m1 i + m2 j + m3 k with ms = √ For any λσl ̸ = λσp , we have

apqs a2pq0 +a2pq1 +a2pq2 +a2pq3

, s = 0, 1, 2, 3.

(k−1) −1) |a(k) − λσp − ta(k | pp − λσl | = |λσp − λσl + app qp

−1) −1) ≥ |λσp − λσl | − |a(k − λσp | − |t ||a(k | pp qp −1) ≥ |λσp − λσl | − |a(k − λσp | − off(A(k−1) ) pp

≥ δ−ϵ−ϵ ≥ 2ϵ.

(3.15)

We note that A and A have the same eigenvalues and off(A ) < ϵ . By using the Wielandt–Hoffman theorem again, we (k) (k) have |app − λσp | < ϵ . The proof is complete under the fact that the quaternion matrix A(k) and its real part A0 have the same diagonal entries. □ (k)

(k)

By Theorem 3.3, we present the structure-preserving Jacobi method for any n-by-n quaternion Hermitian matrix A. Algorithm 3.4 (Cyclic Structure-preserving Jacobi). Given an n-by-n quaternion Hermitian matrix A = A0 + A1 i + A2 j + A3 k ∈ Hn×n and a tolerance tol > 0, this algorithm overlaps the real representation matrix ΓA by ΓGT ΓA ΓG , where ΓG is a generalized JRS-symplectic rotation and off(ΓGT ΓA ΓG ) ≤ tol · ∥ΓA ∥F . V = ΓIn , ζ = tol · ∥ΓA ∥F while off(ΓA ) > tol for p = 1 : n − 1 for q = p + 1 : n app = A0 (p, p); aqq = A0 (q, q); apq0 = A0 (p, q); apq1 = A1 (p, q); apq2 = A2 (p, q); apq3 = A3 (p, q); [cr , s0 , s1 , s2 , s3 ] = RQJR(app , aqq , apq0 , apq1 , apq2 , apq3 ); T ΓA = ΓG(p ,q,θ ) ΓA ΓG(p,q,θ ) ;

V = V ΓG(p,q,θ ) ; end for end for end while

Remark 3.1. Algorithm 3.4 is cyclic; maximum element structure-preserving Jacobi algorithm costs much more CPU time than the cyclic structure-preserving Jacobi algorithm. In Algorithm 3.4, it only needs to store A0 A2 A1 A3 of ΓA in each sweep of n(n − 1)/2 iterations of structurepreserving Jacobi rotations. In the numerical examples later, one can see that Algorithm 3.4 costs less CPU time than the classical complex Jacobi algorithm in [13] under the same accuracy. Parallelism is one of the most attractive advantages of the structure-preserving Jacobi algorithm. The main idea of the parallel structure-preserving Jacobi algorithm is same as the parallel Jacobi algorithm in [7] except the cosine–sine groups we get from Algorithm 3.2. Now we present a simple example for illustration. Let A ∈ H8×8 be Hermitian. By a parallel computer with 4 processors, we can group the 28 coordinate planes into 7 groups as follows:

[

• • • •

]

group(1): (1,2), (3,4), (5,6), (7,8); group(2): (1,3), (2,4), (5,7), (6,8); group(3): (1,4), (2,3), (5,8), (6,7); group(4): (1,5), (2,6), (3,7), (4,8);

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

9

• group(5): (1,6), (2,5), (3,8), (4,7); • group(6): (1,7), (2,8), (3,5), (4,6); • group(7): (1,8), (2,7), (3,6), (4,5). In each group, we compute four generalized JRS-symplectic Jacobi rotations in parallel by 4 processors, which are not conflicting with each other. Such parallel algorithm only needs 1/4 computing time of the computer with a single processor. 4. Numerical experiment In this section, we demonstrate the efficiency of the newly proposed structure-preserving algorithm by three numerical examples. All these experiments are performed on an Inter(R)Core(TM)i5–3230M 2.60GHz/4.00GB computer using MATLAB R2010b, and the tolerance is taken to be tol = 10−14 . Example 4.1. Let A be a 5-by-5 quaternion Hermitian matrix:

⎡

⎤

3 −2i −j −2k 0 ⎢ 2i 3 −2i −j −2k⎥ ⎢ ⎥ 3 −2i −j⎥ . A = ⎢ j 2i ⎣2k j 2i 3 −2i⎦ 0 2k j 2i 3 In this experiment, we apply the cyclic structure-preserving Jacobi method √ √ (Algorithm 3.4) to compute the eigenvalues of A ˜ = {3 − 2 2, 0, 3, 6, 3 + 2 2}. This is a toy example to show the validity and and compare them with the explicit ones Λ reliability of our algorithm. Recall that the real counterpart ΓA of A has the form in (2.1), where

⎤

⎡

⎤

⎡

0 3 0 0 0

⎡

0 0 −1 0 0 ⎢0 0 0 −1 0⎥ ⎢ ⎥ 0 0 0 −1⎥ , A3 = ⎢0 ⎣2 1 0 0 0⎦ 0 0 1 0 0

3 ⎢0 ⎢ A0 = ⎢0 ⎣0 0 0 ⎢0 ⎢ A2 = ⎢1 ⎣0 0

0 0 3 0 0

0 0 0 3 0

0 −2 0 0 0 0 ⎢2 0 −2 0 0⎥ 0⎥ ⎥ ⎢ ⎥ 0⎥ , A1 = ⎢0 2 0 −2 0⎥ , ⎣ ⎦ 0 0 2 0 −2⎦ 0 0 0 0 2 0 3

⎤

⎡

0 0 0 0 2

⎤

0 −2 0 0 0 −2⎥ ⎥ 0 0 0⎥ . 0 0 0⎦ 0 0 0 (5)

Using Algorithm 3.4, we get the numerical results ΓA

(5)

⎡

−1.2426

⎢ ⎢ ⎣

0 0 0 0

A0 = ⎢

with

⎤

⎡

0 0 0 0 0 ⎢0 7.2426 0 0 0 ⎥ ⎢ ⎥ (5) (5) (5) 0 0.0000 0 0 ⎥ , A1 = A2 = A3 = ⎢0 ⎣0 0 0 6.0000 0 ⎦ 0 0 0 0 3.0000

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

⎤

0 0⎥ ⎥ 0⎥ . 0⎦ 0

(5)

The diagonal entries of A0 are computed eigenvalues of A, collected in the set Λ. Since

˜ = 7.1054e − 015, min{min|λ − λ|} ˜ = 1.72201e − 016, max{min|λ − λ|} ˜ Λ ˜ λ∈Λ λ∈

˜ Λ ˜ λ∈Λ λ∈

the errors of the computed eigenvalues have a upper bound 7.1054e − 015, which is smaller than the given tol = 10−14 . Example 4.2. In this experiment, we deal with the diagonalization of 150 random quaternion Hermitian matrices, whose dimension is n = 11 : 160 and figure out the CPU times by utilizing Algorithm 3.4 and the classical complex Jacobi algorithm in [13]. The CPU times are shown in Fig. 1, where the star line and the dot line show the CPU times costed by the classical complex Jacobi algorithm in [13] and Algorithm 3.4, respectively. From Fig. 1, one can see that Algorithm 3.4 is faster than the classical complex Jacobi method when n is large, and this advantage becomes more obvious as the dimensions of quaternion matrices become larger. Meantime, Algorithm 3.4 and the classical complex Jacobi algorithm almost have same calculation accuracy. Taking n = 100 for an example, the computed eigenvalues of Algorithm 3.4 and the classical complex Jacobi algorithm, shown in Fig. 2, are almost equal to each other. These results indicate that Algorithm 3.4 is superior to the classical complex Jacobi algorithm under the same high accuracy. Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

10

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

Fig. 1. CPU times for computing eigenvalues and eigenvectors.

Fig. 2. Eigenvalues in descending order.

Example 4.3. In face recognition, color information can help us to get a higher rate of recognition. Each color face image can be represented by a quaternion matrix, and in this way the covariance matrix of training samples is a quaternion Hermitian matrix. To compute the eigenfaces is to calculate the eigenvectors of the covariance matrix. In this experiment, we take training color face images from Faces300, which contains three hundred 128 × 128 color face images of young students from JSNU, and calculate their eigenfaces by Algorithm 3.4 and the classical complex Jacobi method in [13]. The first twelve eigenfaces in color are shown in Fig. 3. In Fig. 4 we present the CPU times and the residual errors for computing eigenfaces with tolerances changing from 1.0e−1 to 1.0e−11. The results of two methods are almost of the same residual error, but Algorithm 3.4 is faster than the classical complex Jacobi method. 5. Conclusions A new structure-preserving Jacobi method is proposed for solving the eigenvalue problems of quaternion Hermitian matrices. By taking full use of the double structures throughout whole iteration processing, we construct a generalized JRS-symplectic Jacobi rotation to diagonalize the real counterpart. What is more, the 4n-dimensional eigenvalue problem is reduced to an n-dimensional problem. Numerical experiments have shown that our algorithm is faster than the classical complex Jacobi algorithm and the larger is the dimension of the problem, the more obvious is the advantage. Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

11

Fig. 3. Eigenfaces in color of Faces300 (scaled). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. CPU times and residual errors for computing eigenfaces with tol = 1.0e−1:1.0e−11.

Acknowledgments We are grateful to the editor and the anonymous referees for their useful comments and suggestions, which help us improve the original presentation. The research of second author is partially supported by National Natural Science Foundation of China under grants 11771188 and 11201193, TAPP (PPZY2015A013) and PAPD of Jiangsu Higher Education Institutions. The research of third author is partially supported by the National Natural Science Foundation of China (No. 11671337), the Natural Science Foundation of Fujian Province of China (No. 2016J01035), and the Fundamental Research Funds for the Central Universities (No. 20720150001). Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.

12

R.-R. Ma et al. / Computers and Mathematics with Applications (

)

–

References [1] W.R. Hamilton, On quaternions, Proc. R. Ir. Acad. (1844) November 11. [2] F.Z. Zhang, Quaternions and matrices of quaternions, Linear Algebra Appl. 251 (1997) 21–57. [3] D. Finkelstein, J.M. Jauch, D. Speiser, Notes on quaternion quantum mechanics (CERN, Report 59-7), Logico-Algebraic Approach To Quantum Mechanics II, 1979, pp. 367–421. [4] S.D. Leo, P. Rotelli, Quaternion scalar field, Phys. Rev. D 45 (1992) 575–579. [5] N.L. Bihan, S.J. Sangwine, Quaternion principal component analysis of color images, in: IEEE International Conference on Image Processing, Vol. 1, 2003, pp. 809–812. [6] C.G.J. Jacobi, Über ein leichtes verfahren die in der theorie der säcularstörungen vorkommenden gleichungen numerisch aufzulösen, J. Reine Angew. Math. 30 (1846) 51–94. [7] G.H. Golub, C. Van Loan, Matrix Computation, fourth ed., Johns Hopkins University Press, Baltimore, 2013. [8] J. Demmel, K. Veselić, Jacobi’s method is more accurate than QR, SIAM J. Matrix Anal. Appl. 13 (1992) 1204–1245. [9] R. Mathias, Accurate eigensystem computations by Jacobi methods, SIAM J. Matrix Anal. Appl. 16 (1995) 977–1003. [10] R.P. Brent, F.T. Luk, The solution of singular-value and symmetric eigenvalue problems on multiprocessor arrays, SIAM J. Sci. Stat. Comput. 6 (1985) 69–84. [11] P.J. Eberlein, On one-sided Jacobi methods for parallel computation, SIAM J. Algebr. Discrete Methods 8 (1987) 790–796. [12] A. Bunse-Gerstner, R. Byers, V. Mehrmann, A quaternion QR algorithm, Numer. Math. 55 (1989) 83–95. [13] N.L. Bihan, S.J. Sangwine, Jacobi method for quaternion matrix singular value decomposition, Appl. Math. Comput. 187 (2007) 1265–1271. [14] Z.G. Jia, M.S. Wei, S.T. Ling, A new structure-preserving method for quaternion Hermitian eigenvalue problems, J. Comput. Appl. Math. 239 (2013) 12–24. [15] Y. Li, M.S. Wei, F.X. Zhang, J.L. Zhao, Real structure-preserving algorithms of Householder based transformations for quaternion matrices, J. Comput. Appl. Math. 305 (2016) 82–91. [16] I.L. Kantor, A.S. Solodovnikov, Hypercomplex Numbers, An Elementary Introduction To Algebras, Springer, Berlin, 1989. [17] T. Jiang, Algebraic methods for diagonalization of a quaternion matrix in quaternionic quantum theory, J. Math. Phys. 46 (2005) 052106–052108.

Please cite this article in press as: R.-R. Ma, et al., A structure-preserving Jacobi algorithm for quaternion Hermitian eigenvalue problems, Computers and Mathematics with Applications (2017), https://doi.org/10.1016/j.camwa.2017.10.009.