A fast structure-preserving method for computing the singular value decomposition of quaternion matrices

Abstract

In this paper we propose a fast structure-preserving algorithm for computing the singular value decomposition of quaternion matrices. The algorithm is based on the structure-preserving bidiagonalization of the real counterpart for quaternion matrices by applying orthogonal JRS-symplectic matrices. The algorithm is efficient and numerically stable.

Introduction

Quaternion and quaternion matrices have wide applications in applied science, such as special relativity and non-relativistic [11], [13], group representations [9], [10], [12], relativistic dynamics [14], [15], field theory [16], Lagrangian formalism [17], electro weak model [18] and grand unification theories [21]. With the rapid development of the above disciplines, it is getting more and more necessary for us to further study the theoretical properties and numerical computations of quaternion and quaternion matrices.

There are several ways to design an algorithm for numerical computations involving quaternion matrices. The first way is to directly use the methods for real matrix computations, then the algorithm will need quaternion arithmetic, and so sometimes is inefficient. The second way is based on the real presentations of quaternion matrices to design an algorithm. Notice that, both rows and columns of the real presentation of a quaternion matrix expand four times of that of the original quaternion matrix. Therefore, in this case we need to pay extra attention to design an efficient algorithm. When quaternion matrices have special properties, then one may design efficient algorithms by applying the properties of the matrices.

As a special case, consider the following linear complex symmetric system

$$\mathit{Ax}=b\text{,}A\in {\mathbb{C}}^{n\times n}\text{,}\text{and}x\text{,}b\in {\mathbb{C}}^n\text{,}$$

where $A=W+\mathit{iT}$, in which W is real symmetric positive definite, and $T\ne 0$ is real symmetric positive semi-definite, $x=y+\mathit{iz}\text{,}b=c+\mathit{id}$ with $y\text{,}z\text{,}c\text{,}d$ n-dimensional real vectors. Complex symmetric linear systems of this kind arise in many problems in scientific computing and engineering applications, including diffuse optical tomography, FFT-based solution of certain time-dependent PDEs, quantum mechanics, molecular scattering, structural dynamics, and lattice quantum chromo-dynamics. If we design any iterative method direct from (1.1), we would use complex arithmetic throughout the code and it would be wasteful and inefficient.

During the recent years, several efficient iterative algorithms for solving (1.1) were proposed.

• Benzi et al. [4] proposed the following procedure. From the real presentation of a complex matrix, $f(W+\mathit{iT})=\left(\begin{array}{cc}\hfill W\hfill & \hfill -T\hfill \\ \hfill T\hfill & \hfill W\hfill \end{array}\right)$, the above complex symmetrical Eq. (1.1) can be rewritten as

  $$\left(\begin{array}{cc}\hfill W\hfill & \hfill -T\hfill \\ \hfill T\hfill & \hfill W\hfill \end{array}\right)\left(\begin{array}{c}\hfill y\hfill \\ \hfill z\hfill \end{array}\right)=\left(\begin{array}{c}\hfill c\hfill \\ \hfill d\hfill \end{array}\right)\text{.}$$

  Because the left side matrix is real, then one can propose any appropriate iterative method and precondition, such as Hermitian/skew-Hermitian splitting (HSS) [3] and preconditioned HSS (PHSS).

• Bai et al. [1] observed that, the complex symmetric system (1.1) can be rewritten as

  $$\left\{\begin{array}{c}(\alpha I+W)x=(\alpha I-\mathit{iT})x+b\text{,}\hfill \\ (\alpha I+T)x=(\alpha I+\mathit{iW})x-\mathit{ib}\text{,}\hfill \end{array}\right.$$

  in which $\alpha >0$ is a parameter. From the above equations the authors of [1] propose the modified HSS (MHSS) iterative method.

• Bai et al. [2] observed that, the complex symmetric system (1.1) can also be rewritten as

  $$\left\{\begin{array}{c}(\alpha V+W)x=(\alpha V-\mathit{iT})+b\text{,}\hfill \\ (\alpha V+T)x=(\alpha V+\mathit{iW})x-\mathit{ib}\text{,}\hfill \end{array}\right.$$

  in which $\alpha >0$ is a parameter, and V is a real symmetric positive definite matrix. From the above equations the authors of [2] propose the preconditioned MHSS (PMHSS) iterative method.

For detailed analysis, examples and additional references, we refer to [4], [1], [2].

The singular value decomposition (SVD) [19], [20] is one of the most powerful tools in matrix computations. The SVD of a quaternion matrix was theoretically derived in 1997 by Zhang [28]. Now the quaternion SVD has been widely used in reduced-rank signal processing where the idea is to extract the significant parts of a signal [6], algebraic feature extraction method for any colour image in image pattern recognition [26], [7], [25].

We now describe application background of the Quaternion SVD.

• reduced-rank signal processing. A color image can be represented by a quaternion matrix $Q\in {\mathbb{H}}^{m\times n}$. Let the SVD of Q be

  $$Q=U\mathrm{\Sigma }{V}^H\text{,}$$

  in which $U\in {\mathbb{H}}^{m\times m}\text{,}V\in {\mathbb{H}}^{n\times n}$ are quaternion unitary matrices, respectively, $V^H$ is the conjugate transpose matrix of V, $\mathrm{\Sigma }=\mathrm{diag}({\sigma }_1\text{,}\dots \text{,}{\sigma }_l)$ with $l=\min \{m\text{,}n\}$, and ${\sigma }_1⩾{\sigma }_2⩾\cdots ⩾{\sigma }_l$ are the singular values of the matrix Q.

  Then the reduced-rank signal processing model can be described as: find a quaternion matrix ${\hat{Q}}_0\in {\mathbb{H}}^{m\times n}$, such that

  $$\mathrm{rank}({\widehat{Q}}_0)=\min \{\mathrm{rank}(\hat{Q}):\hat{Q}\in {\mathbb{H}}^{m\times n}\}\text{,}\text{s}.\text{t}.\Vert Q-{\hat{Q}}_0{\Vert }_2<∊\text{,}$$

  where $∊>0$ is a problem related parameter and $\Vert ·|{|}_2$ is the spectral matrix norm. If the singular values of Q satisfy ${\sigma }_1⩾\cdots ⩾{\sigma }_k⩾∊>{\sigma }_{k+1}⩾\cdots ⩾{\sigma }_l$, then a solution to (1.6) has the form

  $${\hat{Q}}_0=U_1{\mathrm{\Sigma }}_1{V}_1^H\text{,}$$

  where ${\mathrm{\Sigma }}_1=\mathrm{diag}({\sigma }_1\text{,}\dots \text{,}{\sigma }_k)\text{,}{U}_1\text{,}{V}_1$ are respectively the first k columns of $U\text{,}V$.

• The Karhunen–Lóeve transform. The Karhunen–Lóeve transform is a well-known technique in image and signal processing, and is based on the EVD of the covariance matrix of the different lines or columns. Given an image $Q\in {\mathbb{H}}^{m\times n}$, one can build the covariance matrix of $q_l$, the lth column of Q as:

$${\mathrm{\Gamma }}_l={\tilde{q}}_l{\tilde{q}}_l^H\text{,}$$

where ${\tilde{q}}_l=q_l-m[q_l]$ is the lth centered (corrected form its mean color value) column of Q. Then, the mean covariance matrix is given as:

$$\mathrm{\Gamma }=\frac{1}{n}{\displaystyle \sum _{l=1}^n}{\mathrm{\Gamma }}_l=\frac{1}{n}\tilde{Q}{\tilde{Q}}^H\text{,}$$

in which $\tilde{Q}=({\tilde{q}}_1\text{,}\dots \text{,}{\tilde{q}}_n)$, and $\mathrm{\Gamma }$ is a positive semi-definite quaternion matrix. Then $\mathrm{\Gamma }$ has the following quaternion eigenvalue decomposition:

$$\mathrm{\Gamma }=U\mathrm{\Lambda }{U}^H\text{,}$$

in which U is a quaternion unitary matrix and $\mathrm{\Lambda }=\mathrm{diag}({\lambda }_1\text{,}\dots \text{,}{\lambda }_n)$ is a positive semi-definite diagonal matrix. From the above transformation we obtain a new image

$$Y=U^{T}Q\text{,}$$

here $U^T$ denotes the transpose of the matrix U. On the other hand, the SVD of $\tilde{Q}$ can be denoted as

$$\tilde{Q}=U\mathrm{\Sigma }{V}^H\text{,}$$

where U is the same quaternion unitary matrix as in the quaternion eigenvalue decomposition of $\mathrm{\Gamma }$, therefore we can compute the SVD of $\tilde{Q}$ instead the QEVD of $\mathrm{\Gamma }$ to obtain U.

Bihan [6], [5] introduced the computations of quaternion SVD using the complex adjoint matrix. [23] presented a method to compute quaternion SVD based on transformation of the quaternion matrix to bidiagonal form using quaternion Householder transformations. Bihan and Sangwine [8] proposed the Jacobi SVD algorithm for quaternion matrices. Recently, the authors of [22] proposed a structure-preserving algorithm for solving the right eigenvalue problem of Hermitian quaternion matrices, using the special structures and properties of the real representation of Hermitian quaternion matrices. By applying the same idea, the authors of [27] proposed a structure-preserving algorithm to compute the Cholesky decomposition of the Hermitian positive definite quaternion matrix.

In this paper we will propose the SVD of a quaternion matrix Q using the same idea of [22]. Let ${\mathrm{\Upsilon }}_Q$ be the real representation of Q. Then ${\mathrm{\Upsilon }}_Q$ has a nice structure: generalized JRS-symmetry. This structure is unchanged under orthogonal JRS-symplectic transformations (for definitions of generalized JRS-symmetry and orthogonal JRS-symplectic transformation, see §2). That makes it possible to derive a structure-preserving method for bidiagonalization of ${\mathrm{\Upsilon }}_Q$. This structure-preserving method only cost about a quarter of arithmetic operations of the Householder bidiagonalization algorithm directly applied to ${\mathrm{\Upsilon }}_Q$. What is more important is that the bidiagonal matrix obtained from $m\times n$ quaternion matrix by the structure-preserving method is still generalized JRS-symmetric, and has the following form

$$\left[\begin{array}{cccc}\hfill D\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill D\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill D\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill D\hfill \end{array}\right]$$

where D is an $m\times n$ real bidiagonal matrix. Therefore to evaluate the singular values of the quaternion matrix, we only need to compute the singular values of D.

This paper is organized as follows. In Section 2, we will recall some preliminary results used in the paper. In Section 3, we will propose a structure-preserving algorithm for computing the SVD for quaternion matrices. In Section 4, we will provide experiments to compare this algorithm with two other standard algorithms to demonstrate the efficiency of our algorithm. Finally in Section 5 we will make some concluding remarks.