1 Introduction↩︎

A quaternion number is a number of form \(a_0+a_1{{\tt i}}+a_2{{\tt j}}+a_3{{\tt k}}\), where \(a_s\)’s are all real numbers; \({{\tt i}},{{\tt j}},{{\tt k}}\) are imaginary units such that \[{{\tt i}}^2={{\tt j}}^2={{\tt k}}^2=-1,\quad{{\tt i}}{{\tt j}}=-{{\tt j}}{{\tt i}}={{\tt k}},\quad {{\tt j}}{{\tt k}}=-{{\tt k}}{{\tt j}}={{\tt i}},\quad {{\tt k}}{{\tt i}}=-{{\tt i}}{{\tt k}}={{\tt j}};\] We call a quaternion as pure quaternion if \(a_0=0\). In this work, we are interested in solving Hermitian quaternion Toeplitz systems. \[\label{hpdqtsystem} T_n \mathbf{u}=\mathbf{b},\tag{1}\] where \(b\) is a given quaternion vector of size \(n\times 1\); \(u\) is the unknown to be solved; \(T_n\) is an \(n\times n\) Hermitian quaternion Toeplitz matrix of the following form \[T_n=\left[ \begin{array} [c]{ccccc} t_0 & t_{-1} &\ldots & t_{2-n} & t_{1-n} \\ t_1 & t_0 & t_{-1} &\ddots & t_{2-n} \\ \vdots&\ddots &\ddots&\ddots&\vdots\\ t_{n-2} & \ddots& t_1 & t_0 & t_{-1} \\ t_{n-1} & t_{n-2} & \ldots & t_1 & t_{0} \end{array} \right],\] with \(t_s\)’s are quaternion numbers; \(t_{-s}\) is conjugate of \(t_s\) (see Section 2 for the definition) for each \(s>0\); \(t_0\) is a real number. Quaternion Toeplitz systems have attracted a growing attention in recently years, thanks to its various applications including signal and image processing [1]–[7].

The main aim of this paper is to show that the eigenvalues of Hermitian quaternion Toeplitz matrices can be characterized by their quaternion-valued generating function. This refers to a quaternion-version of Grenander-Szegö theorem. We note that such generating function must be the sum of a real-valued function and an odd function with imaginary component. This setting is different from the case of Hermitian complex Toeplitz matrices generated by real-valued functions only. As an application, we investigate circulant preconditioned conjugate gradient methods for solving Hermitian quaternion Toeplitz systems arising from quaternion signal processing. Here we employ Strang’s circulant preconditioners (as examples) and show that the eigenvalues of Strang’s circulant matrices are computed by the values evaluated at the partial sum of generating functions. We remark that a quaternion circulant matrix in general is not diagonalized by discrete Fourier transform matrix. Also we verify when the circulant preconditioned conjugate gradient method is applied to solving such Hermitian quaternion Toeplitz systems, its superlinear convergence can be achieved.

The outline of this paper is given as follows. In Section 2, we study the spectra spectra of Hermitian quaternion Toeplitz matrices constructed by quaternion-valued generating function. In Section 3, we study the spectra of circulant preconditioned Toeplitz matrices. Numerical results for Hermitian quaternion Toeplitz matrices arising quaternion signal processing are reported in Section 4. Finally, some concluding remarks are given in Section 5.

1.1 Preliminaries↩︎

The real number field, the complex number field and the quaternion algebra are denoted by \(\mathbb{R}\), \(\mathbb{C}\) and \(\mathbb{Q}\), respectively. The set of all positive real numbers is denoted by \(\mathbb{R}^{+}\). The spaces of all \(m\times n\) real matrices, complex matrices and quaternion matrices are denoted by \(\mathbb{R}^{m\times n}\), \(\mathbb{C}^{m\times n}\) and \(\mathbb{Q}^{m\times n}\), respectively. If there is no ambiguity, we treat \(\mathbb{R}^{1\times 1}\), \(\mathbb{C}^{1\times 1}\) and \(\mathbb{Q}^{1\times 1}\) as \(\mathbb{R}\), \(\mathbb{C}\) and \(\mathbb{Q}\), respectively. The set of all positive integers and the of all nonnegative integers are denoted by \(\mathbb{N}^{+}\) and \(\mathbb{N}\), respectively.

The dot product of two quaternion numbers \(x=x_0+x_1{{\tt i}}+x_2{{\tt j}}+x_3{{\tt k}}\) and \(y=y_0+y_1{{\tt i}}+y_2{{\tt j}}+y_3{{\tt k}}\) with \(x_s,y_s\in\mathbb{R}\) for \(s=0,1,2,3\), is defined as \[x \cdot y:=\sum\limits_{s=0}^{3}x_sy_s.\] The modulus of \(x\) is defined by \(|x|:=\sqrt{x\cdot x}=\sqrt{|x_0|^2+|x_1|^2+|x_2|^2+|x_3|^2}\). We say \(x\) is a unit if and only if \(|x|=1\). Two pure quaternions \(x,y\) are said to be orthogonal if and only if \(x\cdot y=0\).

It can be seen that \(({{\tt i}},{{\tt j}},{{\tt k}})\) is orthonormal. Any orthonormal triple of pure quaternions \(({{\tt p}},{{\tt q}},{{\tt r}})\) can be treated as a three-axis imaginary system like \(({{\tt i}},{{\tt j}},{{\tt k}})\), see Property 1 in [8].

For any \(A=\tilde{A}_0+\tilde{A}_1{{\tt p}}+\tilde{A}_2{{\tt q}}+\tilde{A}_3{{\tt r}}\in\mathbb{Q}^{m_1\times m_2}\) with \(\tilde{A}_s\in\mathbb{R}^{m_1\times m_2}~(s=0,1,2,3)\), its transpose \(A^{{\rm T}}\) and its conjugate \(\bar{A}\) are defined as \[\begin{align} A^{{{\rm T}}}:=\tilde{A}_0^{{\rm T}}+\tilde{A}_1^{{\rm T}}{{\tt p}}+\tilde{A}_2^{{\rm T}}{{\tt q}}+\tilde{A}_3^{{\rm T}}{{\tt r}},\quad \bar{A}:=\tilde{A}_0-\tilde{A}_1{{\tt p}}-\tilde{A}_2{{\tt q}}-\tilde{A}_3{{\tt r}}; \end{align}\] For any \(A\in\mathbb{Q}^{m_1\times m_2}\), the conjugate transpose of \(A\), denoted by \(A^{*}\), is \(A^{*}:=\bar{A}^{\rm T}\). \(A\) is said to be Hermitian if and only if \(A=A^{*}\).

For any \({\boldsymbol{x}}=(x_1,x_2,...,x_m)^{\rm T}\in\mathbb{Q}^{m\times 1}\), its vector 2-norm is defined as \(||{\boldsymbol{x}}||_2=\sqrt{\sum\limits_{i=1}^{m}|x_i|^2}\). Clearly, \(||{\boldsymbol{x}}||_2=\sqrt{{\boldsymbol{x}}^{*}{\boldsymbol{x}}}\) for any quaternion-valued vector \({\boldsymbol{x}}\).

Denote \[\mathbb{Q}_{(0,1)}:=\{x+y{{\tt p}}|x,y\in\mathbb{R}\},\quad \mathbb{Q}_{(2,3)}:=\{x{{\tt q}}+y{{\tt r}}|x,y\in\mathbb{R}\}.\] The set of all \(m\times n\) matrices whose entries belonging to \(\mathbb{Q}_{(0,1)}\) (or \(\mathbb{Q}_{(2,3)}\)) is denoted as \(\mathbb{Q}_{(0,1)}^{m\times n}\) (or \(\mathbb{Q}_{(2,3)}^{m\times n}\)), We remark that \(\mathbb{Q}_{(0,1)}\) is a field isomorphic to the complex field. Hence, results on numbers in \(\mathbb{C}\) (or matrices in \(\mathbb{C}^{m\times n}\)) can be applied to the numbers in \(\mathbb{Q}_{(0,1)}\) (or matrices in \(\mathbb{Q}_{(0,1)}^{m\times n}\)) under isomorphism. For a square matrix \(A\in\mathbb{K}^{m\times m}\) (\(\mathbb{K}=\mathbb{Q}\) or \(\mathbb{Q}_{(0,1)}\)), we denote \[\sigma(A):=\{\lambda\in\mathbb{K}|A{\mathbf{z}}=\mathbf{ z}\lambda{\rm~for~some~nonzero~vector~}\mathbf{z}\in\mathbb{K}^{m\times 1}\}.\] Here we refer \((\lambda,\mathbf{z})\) be the right eigenvalue and corresponding eigenvector of \(A\), and \(\sigma(A)\) be the set of all right eigenvalues of \(A\).

For any \(A=A_0+A_1{{\tt p}}+A_2{{\tt q}}+A_3{{\tt r}}\in\mathbb{Q}^{m_1\times m_2}\) with \(A_s\in\mathbb{R}^{m_1\times m_2}\) for \(s=0,1,2,3\), \(A\) can be equivalently expressed as \[A=(A_0+A_1{{\tt p}})+(A_2+A_3{{\tt p}}){{\tt q}},\] and such representation is unique.

From the above definition, \(A=\phi_1(A)+\phi_2(A){{\tt q}}\) holds for any positive integers, \(m_1\) and \(m_2\) and any \(A\in\mathbb{Q}^{m_1\times m_2}\). If \(A\) is a scalar, \(\phi_1(A)\) and \(\phi_2(A)\) can be defined in a similar way to ?? .

With \(\phi_1(\cdot)\) and \(\phi_2(\cdot)\), we define two mappings \(\mathcal{M}(\cdot)\) and \(\mathcal{V}(\cdot)\) as follows: \[\begin{align} &\mathcal{M}(\cdot):A\in\mathbb{Q}^{n\times n}\mapsto\mathcal{M}{(A)}:=\left[\begin{array}[c]{cc} \phi_1(A)&-\phi_2(A)\\ &\\ \overline{\phi_2(A)}&\overline{\phi_1(A)} \end{array}\right]\in\mathbb{Q}_{(0,1)}^{2n\times 2n};\\ &\mathcal{V}(\cdot):\mathbf{x}\in\mathbb{Q}^{n\times 1}\mapsto\mathcal{V}{\left(\mathbf{x}\right)}:=\left[\begin{array}[c]{c} \phi_1(\mathbf{x})\\ ~\\ \overline{\phi_2(\mathbf{x})} \end{array}\right]\in\mathbb{Q}_{(0,1)}^{2n\times 1}. \end{align}\] Note that if \(A \in \mathbb{Q}^{n \times n}\) is Hermitian, then \(\mathbf{x}^* A \mathbf{x} = \mathcal{V}(\mathbf{x})^* \mathcal{M}(A) \mathcal{V}(\mathbf{x})\) and \(\sigma(A) =\sigma( \mathcal{M}(A))\).

2 Quaternion Toeplitz Matrices↩︎

2.1 Generating Functions↩︎

For a quaternion-valued function \(f\), we define two \(\mathbb{Q}_{(0,1)}\)-valued functions \(\Phi_1[f]\) and \(\Phi_2[f]\) as follows: \[\Phi_1[f](x):=\phi_1(f(x)),\quad \Phi_2[f](x):=\phi_2(f(x)).\] With the notations introduced above, we define a pair of functions for generating a sequence of quaternion Hermitian Toeplitz matrices.

According to Lemma 2.2, we remark that a sequence of Hermitian quaternion Toeplitz matrices is generated by a quaternion-valued function which is the sum of a real-valued function and an odd function with imaginary component. This setting is different from the case of Hermitian complex Toeplitz matrices generated by real-valued functions only.

For any vector \(\mathbf{y}\in\mathbb{Q}^{n\times 1}\), denote \[\mathcal{F}\left(\mathbf{ y},x\right):=\sum\limits_{s=1}^{n}\mathbf{y}(s)\exp({{\tt p}}s x),\quad x \in[-\pi,\pi].\] Next we characterize the eigenvalues of quaternion Hermitian Toeplitz matrices.

Because we are interested in the spectra of Hermitian quaternion Toeplitz matrices \(T_n\), we give the conditions on \(f\) such that \(T_n\) are positive definite.

Another important theorem on spectral distribution of Hermitian Toeplitz matrices is the Grenander-Szegö Theorem. Here we give a quaternion-version Grenander-Szegö theorem based on quaternion-valued generating functions. The proof is analogous to that of the Grendander-Szego theorem in the complex field. We present the details in Appendix A.

2.2 Quaternion Signal Processing↩︎

In this subsection, we give examples of generating functions arising from quaternion signal processing.

We consider a discrete-time quaternion signal \(\{ x(t) = x_0(t) + x_1(t){{\tt p}}+x_2(t){{\tt q}}+x_3(t){{\tt r}}\}\) with a finite 2nd-moment \(\mathcal{E}(|x(t)|^2)<+\infty\) for each \(t\), where \(\mathcal{E}(\cdot)\) denotes the expectation operator, see [11], [12]. The linear prediction is to find a set of linear combination coefficients \(\{\alpha_s\}_{s=1}^{n}\subset\mathbb{Q}\) such that the linear combination \(\sum\limits_{s=1}^{n}\alpha_{s}x(t-s)\) best fitting \(x(t)\) under the expected value: \[\begin{align} &\mathcal{E}\left(\left(x(t)-\sum\limits_{s=1}^{n}\alpha_{s}x(t-s)\right)\overline{\left(x(t)-\sum\limits_{s=1}^{n}\alpha_{s}x(t-s)\right)}\right)\notag\\ &=\min\limits_{\{w_s\}_{s=1}^{n}\subset\mathbb{Q}}\mathcal{E}\left(\left(x(t)-\sum\limits_{s=1}^{n}w_{s}x(t-s)\right)\overline{\left(x(t)-\sum\limits_{s=1}^{n}w_{s}x(t-s)\right)}\right).\label{alphasoptimiz} \end{align}\tag{5}\] Once \(\{\alpha_s\}_{s=1}^{n}\) is computed, one can then predict the value of the signal at next time step by computing the linear combination of the most recent \(n\) many observations with \(\{\alpha_s\}_{s=1}^{n}\) as the left-hand side coefficients. It is straightforward to verify that 5 is equivalent to the systems of equations \[\mathcal{E}\left(\left(x(t)-\sum\limits_{s=1}^{n}\alpha_{s}x(t-s)\right)\overline{x(t-s')}\right)=0, \;s'=1,2,...,n,\] which means \[\begin{align} \label{linearsys} \sum\limits_{s=1}^{n}\mathcal{E}\left(x(t-s')\overline{x(t-s)}\right)\overline{\alpha_s}&=\overline{\sum\limits_{s=1}^{n}\alpha_s\mathcal{E}\left(x(t-s)\overline{x(t-s')}\right)}\notag\\ &=\overline{\mathcal{E}\left(x(t)\overline{x(t-s')}\right)}=\mathcal{E}\left(x(t-s')\overline{x(t)}\right), \;s'=1,2,...,n. \end{align}\tag{6}\] Suppose \(\mathcal{E}\left(x(t)\overline{x(t\textcolor{black}{+}s)}\right)\) is a function only of the time-lag \(s\), i.e., there exists a covariance function of \(s\), \(\eta(s)\) such that \[\label{covarianceshiftinvarianteq} \mathcal{E}\left(x(t)\overline{x(t+s)}\right)=\eta(s),\quad s=0,1, 2,\cdots.\tag{7}\] With 7 , \(\mathcal{E}\left(x(t)\overline{x(t+s)}\right)\) for negative \(s\) is also a function only of time-lag \(s\). This is due to the following facts \[\begin{align} \mathcal{E}\left(x(t)\overline{x(t+s)}\right)&=\overline{\mathcal{E}\left(x(t+s)\overline{x(t+s+(-s))}\right)}=\overline{\eta(-s)},\quad s=-1,-2,\cdots. \end{align}\] Then (6 ) indicates a Hermitian quaternion Toeplitz linear system as follows: \[\label{hermqtsystem} T_n \boldsymbol{\alpha} = \mathbf{w} = \left[\eta(1),\eta(2),...,\eta(n)\right]^{\rm T},\tag{8}\] where \[[ T_n ]_{s',s} =\begin{cases} \mathcal{E}\left(x(t-s')\overline{x(t-s'+(s'-s))}\right)=\eta(|s'-s|),\quad s'\geq s,\\ \mathcal{E}\left(x(t-s')\overline{x(t-s'+(s'-s))}\right)=\overline{\eta(|s'-s|)},\quad s'<s, \end{cases}\] and the components \(\overline{\alpha_{s'}}\) \((s'=1,2,...,n)\) of the unknown vector \(\alpha\). When \(\{\eta(s)\}_{s=0}^{\infty}\) is absolutely summable, it is straightforward to verify that \[\label{genf} f(\theta)=\eta(0)+\sum\limits_{s=1}^{+\infty} \eta(s) \exp({{\tt p}}s \theta)+ \overline{\eta(s)} \exp(-{{\tt p}}s \theta), \;\theta \in [-\pi,\pi],\tag{9}\] meets Definition 4 as a generating function for \(T_n\) in 8 .

3 Circulant Preconditioned Conjugate Gradient Method↩︎

Theory of circulant preconditioners for preconditioning real or complex Toeplitz linear systems has been well studied in the literature; see, e.g., [13]–[22]. However, circulant preconditioners for quaternion Toeplitz systems are rarely to see. The main purpose of this section is to demonstrate how the classical theory of circulant preconditioner for complex HPD Toeplitz matrix is extended to that for quaternion HPD Toeplitz matrix. To know more about the classical preconditioning theories for block Toeplitz matrices, one may refer to [23]–[28] and the references therein.

3.1 Circulant Preconditioners↩︎

A Toeplitz matrix \(C_n\in\mathbb{Q}^{n\times n}\) of the following form is called a circulant matrix. \[\label{circmatform} C_n=\left[ \begin{array} [c]{cccc} c_0 & c_{n-1}&\ldots & c_1 \\ c_1 & c_0 & \ddots &\vdots \\ \vdots&\ddots &\ddots&c_{n-1}\\ c_{n-1} & \ldots & c_1 & c_{0} \end{array} \right]\tag{10}\] which its rows are composed of the same elements and each row is rotated one element to the right relative to the preceding row. Let \(\mathbf{v}:=(c_0,c_1,\cdots,c_{n-1})^T\), the first column of \(C\), we notice that \(C\) can be identified by \(\mathbf{v}\). For simplicity, we denote \(C_n\) by \(\mathcal{C}(\mathbf{v})\). Given a general \(n\times n\) Toeplitz matrix \(T_n=[t_{s-l}]_{s,l=1}^{n}\), let \[\begin{align} &{\boldsymbol{v}}_S\left(T_n\right):=(t_0,t_1,...,t_{\lfloor (n-1)/2\rfloor},t_{-\lfloor (n-1)/2\rfloor},...,t_{-2},t_{-1})^{\rm T}, ~ \text{if n is odd} ;\\ &{\boldsymbol{v}}_S\left(T_n\right):=(t_0,t_1,...,t_{\lfloor (n-1)/2\rfloor},0,t_{-\lfloor (n-1)/2\rfloor},...,t_{-2},t_{-1})^{\rm T}, ~ \text{if n is even.} \end{align}\] Then, the Strang circulant preconditioner denoted by \(c(T_n)\) for a Toeplitz matrix is defined by \[c(T_n):=\mathcal{C}({\boldsymbol{v}}_S\left(T_n\right)).\] Apparently, \(c(T_n)\) is Hermitian whenever \(T_n\) is a Hermitian Toeplitz matrix. Here we employ the Strang circulant preconditioner as an example, other circulant preconditioners can be constructed and studied similarly.

For any \(m\in\mathbb{N}^{+}\), denote \[\textcolor{black}{F_{{{\tt p}},m}:=\frac{1}{\sqrt{m}}\left[\exp\left(\frac{2\pi{{\tt p}}(s-1)(l-1)}{m}\right)\right]_{s,l=1}^{m}}.\] \(F_{{{\tt p}},m}\) is called a quaternion discrete Fourier matrix, which is unitary; see, e.g., [8]. As shown in [8], for general \(\mathbf{v}\in\mathbb{Q}^{m\times 1}\), \(F_{{{\tt p}},m}\mathcal{C}(\mathbf{v})F_{{{\tt p}},m}^{*}\) has an \(X\)-shape sparse pattern shown as follows: \[F_{{{\tt p}},m}\mathcal{C}(\mathbf{v})F_{{{\tt p}},m}^{*}=\left[\begin{array}[c]{cccccc} *&0&\ldots&\ldots&\ldots&0\\ 0&*&0&\ldots&0&*\\ \vdots&0&\ddots&& \iddots&0\\ \vdots&\vdots&& &&\vdots\\ \vdots&0&\iddots&&\ddots&0\\ 0&*&0&\ldots&0&* \end{array}\right].\] from which it is clear to see that for general \(\mathbf{v}\in\mathbb{Q}^{m\times 1}\), applying Fourier transform to the first column of \(\mathcal{C}(\mathbf{v})\) would no longer obtain a vector with components in \(\sigma(\mathcal{C}(\mathbf{v}))\), the right spectrum of \(\mathcal{C}(\mathbf{v})\). Interestingly, the next lemma states that \(c(T_n)\) can be diagonalized by \(F_{{{\tt p}},m}\) and its eigenvalues are the eigenvalues of 2-by-2 block matrices where their entries are the values evaluated at the partial sum of quaternion-valued generating functions.

3.2 The Spectra of Preconditioned Matrices↩︎

In this subsection, we show the spectra of \(c(T_n)^{-1} T_n\) is clustered around 1 under the assumption that the sequence \(\{t_k\}_{k\in\mathbb{N}^{+}}\) is absolutely summable.

Referring to the analysis in [13], [14], one can prove show the following lemma.

3.3 The Circulant Preconditioned Conjugate Gradient Method↩︎

The PCG method is an efficient and powerful iterative algorithm to solve HPD linear systems, that improves convergence rates compared to the standard Conjugate Gradient (CG) method by incorporating a preconditioner. The PCG solver using an HPD preconditioner \(P\in\mathbb{Q}^{n\times n}\) for solving a general \(n\times n\) HPD quaternion linear system \[\label{generalhpdqsystem} A\mathbf{w}=\mathbf{z},\tag{12}\] follows the same framework as that used in solving HPD linear systems, which is described in detail in many relevant textbooks, such as [13]. Hence we omit the details here.

When applying PCG algorithm to solve 1 , the dominant operation cost is on computing some matrix-vector multiplications of form \(P^{-1}\mathbf{z}_1\) and \(T_n\mathbf{z}_2\) for some given vectors \(\mathbf{z}_1,\mathbf{z}_2\in\mathbb{Q}^{n\times 1}\) during each iteration. It has been proven in [8] that any circulant quaternion matrix is fast block diagonalizable by means of fast Fourier transforms (FFTs), with each eigen-block being of size at most \(2\times 2\). As a result, for any invertible circulant matrix \(C\in\mathbb{Q}^{n\times n}\), the matrix vector product \(C^{-1}\mathbf{x}\) can be fast computed within \(\mathcal{O}(n\log n)\) operations for any given vector \(\mathbf{x}\in\mathbb{Q}^{n\times 1}\). Moreover, it is well-known that the matrix-vector \(T_n\mathbf{y}\) for any given \(\mathbf{y}\in\mathbb{Q}^{n\times 1}\) can be fast computed within \(\mathcal{O}{(n\log n)}\) by embedding the \(n\times n\) Toeplitz matrix \(T_n\) into a \(2n\times 2n\) circulant matrix (see, e.g., [13], [14]). Therefore, when applying PCG algorithm with an HPD circulant matrix as preconditioner to solve 1 , it requires \(\mathcal{O}{(n\log n)}\) operations in each iteration.

The complexity of PCG depends not only on the operation cost at each iteration but also on the number of iterations required for achieving a specific stopping criterion. Hence, the convergence rate of PCG solver is another important property that has to be investigated.

We will discuss it in the following theorems. To start the discussion, we first need some notations. For any vector \(\mathbf{x}\), we denote the \(i\)th component of \(\mathbf{x}\) by \(x{\boldsymbol{(}i)}\). For any \(\mathbf{x},\mathbf{y}\in\mathbb{K}^{n\times 1}\) (\(\mathbb{K}=\mathbb{Q}\) or \(\mathbb{Q}_{(0,1)}\)), define \[\left\langle{\mathbf{x}},\mathbf{y}\right\rangle:=\mathbf{y}^{*}\mathbf{x}=\sum\limits_{s=1}^{n}\overline{y(s)}x(s).\] For any HPD matrix \(A\in\mathbb{K}^{n\times 1}\), define \[\label{bnormdef} ||\mathbf{x}||_{A}:=\left\langle{A\mathbf{x}},\mathbf{x}\right\rangle^{\frac{1}{2}},\quad \forall \mathbf{x}\in\mathbb{K}^{n\times 1}.\tag{13}\] Clearly, the above defined \(||\cdot||_{B}\) is a vector norm on \(\mathbb{K}^{n\times 1}\).

Let \(\mathbf{w}^{(k)}~(k\geq 1)\) being the \(k\)th iterative solution to 12 , \(\mathbf{w}^{(0)}\) being the initial guess and \(\mathbf{w}\) being the exact solution to 12 . It is known that when applying the PCG solver to linear system of equations in complex field, the step parameters are generated by the inner product of vectors. Similarly, the step parameters of the PCG in quaternion algebra are also given by the inner product of quaternion vectors. Similar to the algorithm in complex filed, the solution of \(k\)-th iteration \({\boldsymbol{w}}^{(k)}\) is given by \[\mathbf{w}^{(k)}=\mathbf{w}^{(0)}+\sum\limits^{k-1}_{s=0} \zeta_s( {P^{-1}A})^sP^{-1}(\mathbf{z}-A\mathbf{w}^{(0)}),\] where the coefficients \(\zeta_s\) (\(s=0,1,...\)) are generated by Hermitian quadratic forms / inner products (e.g., \({\boldsymbol{r}}_k^*P^{-1}{\boldsymbol{r}}_k\) and \({\boldsymbol{p}}_k^*A{\boldsymbol{p}}_k\)) with \({\boldsymbol{r}}_k\) and \({\boldsymbol{p}}_k\) being some vectors computed in the iteration. As \(P^{-1}\) and \(A\) are both quaternionic Hermitian matrices, \(\zeta_s\) (\(s=0,1,...\)) must be real numbers. Therefore, we can estimate the error of the PCG solver in the following theorem.

4 Numerical Experiments↩︎

In this section, we test examples 1 and 2 arising from quaternion signal processing as outlined in Subsection 2.2 by the PCG algorithm. In practice, \(\eta(\cdot)\) is usually unknown, which means the evaluation of \(T_n\) and \(\mathbf{w}\) given in 8 may not be practical. In general, samplings of a discrete-time signal are available, which is useful in approximate evaluation of \(T_n\) and \(\mathbf{w}\). There are several types of well-known windowing methods with samplings of the signal, such as, correlation, covariance, pre-windowed and post-windowed methods, see, e.g., [13], [30]. Let \(\mathbf{x}_1,\mathbf{x}_2,...,\mathbf{x}_n,\mathbf{x}_{n+1},...,\mathbf{x}_{M}\) be \(M\) (\(M>n\)) many successive samplings of the signal \(x(t)\). Now that value of \(\mathcal{E}\left(x(t)\overline{x(t+s)}\right)=\eta(s)\) is independent of \(t\). Following the correlation windowing method, we use the following sample mean \(\tilde{\eta}(s)\) to approximate the expectation value \(\eta(s)\): \[\eta(s)\approx\tilde{\eta}(s):=\frac{1}{M}\sum\limits_{l=1}^{M-s}\mathbf{x}_l\bar{\mathbf{x}}_{l+s},\quad s=0,1,...,n.\] With the computed \(\tilde{\eta}(s)\) \((s=0,1,...,n)\), we obtain another Hermitian quaternion Toeplitz linear system 14 to approximate 8 , \[\label{apphermqtsystem} \tilde{H}\tilde{\mathbf{v}}=\tilde{\mathbf{w}},\tag{14}\] where \[\begin{align} &\tilde{H}:=\left[ \begin{array} [c]{ccccc} \tilde{\eta}(0) & \overline{\tilde{\eta}(1)} &\ldots & \overline{\tilde{\eta}(n-2)} & \overline{\tilde{\eta}(n-1)} \\ \tilde{\eta}(1 ) & \tilde{\eta}(0) & \overline{\tilde{\eta}(1)} &\ddots & \overline{\tilde{\eta}(n-2)} \\ \vdots&\ddots &\ddots&\ddots&\vdots\\ \tilde{\eta}(n-2) & \ddots& \tilde{\eta}(1) & \tilde{\eta}(0) & \overline{\tilde{\eta}(1)} \\ \tilde{\eta}(n-1) & \tilde{\eta}(n-2) & \ldots & \tilde{\eta}(1) & \tilde{\eta}(0) \end{array} \right],\\ &\tilde{\mathbf{w}}:=\left[\tilde{\eta}(1),\tilde{\eta}(2),..., \tilde{\eta}(n)\right]^{\rm T}. \end{align}\] On the other hand, \(\tilde{H}\) can be rewritten as \(\tilde{H}=\frac{1}{M}\tilde{T}^{*}\tilde{T}\) with \[\tilde{T}:=\left[\begin{array}[c]{ccc} \bar{x}_{1}&&\\ \vdots&\ddots&\\ \bar{x}_{n}&\ldots&\bar{x}_{1}\\ \vdots&\ddots&\vdots\\ \vdots&\ddots&\vdots\\ \bar{x}_{M}&\ldots&\bar{x}_{M-n+1}\\ &\ddots&\vdots\\ &&\bar{x}_{M} \end{array}\right]\in\mathbb{Q}^{(M+n-1)\times n}.\] Hence, if \(\tilde{T}\) is a full rank matrix, then \(\tilde{H}\) is an HPD quaternion matrix. One can then apply the PCG solver with the circulant preconditioner to solve the linear system 14 .

Settings: In the experiments, zero initial guess is taken for the PCG algorithm; the stopping criterion of PCG algorithm is set as \(||\hat{\mathbf{r}}^{(k)}||_2\leq 10^{-7}||\hat{\mathbf{r}}^{(0)}||_2\), where \(\hat{\mathbf{r}}^{(k)}\) denotes the residual vector at \(k\)-th PCG iteration and \(\hat{\mathbf{r}}_0\) denotes the initial residual vector. To demonstrate the effectiveness of Strang’s circulant preconditioner on quaternionic Toeplitz systems, we test the PCG algorithm with the Strang’s circulant preconditioner and unpreconditioned conjugate gradient algorithm in the lateral experiments and compare their performance. For ease of statement, we use PCG-\({\boldsymbol{C}}\) (CG, resp.) to represent the PCG algorithm with Strang’s circulant preconditioner (unpreconditioned conjugate gradient solver, resp.). We shall use ‘CPU’ to represent the computational time; and ‘Iter’ for the iteration numbers of PCG solvers. To quantify the accuracy of PCG solvers, we define error measure of iterative solution as \[{\rm Error}:=||\mathbf{b}-T_n\mathbf{u}^{(k)}||_{2},\] where \(\mathbf{u}^{(k)}\) denotes some iterative solution to 1 or 14 .

In the following, we present the results of Examples 1 and 2 by applying the PCG solvers and list the results in Tables 1-2. From these tables, we see that (i) PCG-\({\boldsymbol{C}}\) is more efficient than CG in terms of iteration number and computational time; (ii) PCG-\({\boldsymbol{C}}\) is more accurate than CG in terms of Error. The better performance of PCG-\({\boldsymbol{C}}\) compared with that of CG demonstrates the effectiveness of the circulant preconditioner. Such performance of circulant preconditioner for quaternion Hermitian Toeplitz system is consistent with what we have usually observed on circulant preconditioning for complex Hermitian Toeplitz system.

0.7em

0.7em

0.7em

0.7em

5 Concluding Remarks↩︎

In this paper, we have studied a sequence of Hermitian quaternion Toeplitz matrices generated by a quaternion-valued function which is the sum of a real-valued function and an odd function with imaginary component. This is different from the case of Hermitian complex Toeplitz matrices generated by real-valued functions only. As an example, we studied the Strang circulant preconditioner for quaternion Hermitian Toeplitz matrix. We have shown that the Strang circulant preconditioner can be diagonalized by discrete Fourier transform matrix whereas a general quaternion circulant matrix cannot be diagonized. Both theoretical and numerical results are presented to demonstrate the effectiveness of quaternion circulant preconditioner for solving quaternion Toeplitz system.

It is interesting to note that the spectra of complex Toeplitz matrices are characterized by their complex generating functions on the unit circle in the complex plane and their winding numbers [31]. As a future research work, it is worth how the spectra of quaternion Toeplitz matrices are related to their quaternion generating functions and the possible generalization of winding number in the hypercomplex domain.

Acknowledgments↩︎

The work of Xue-Lei Lin was partially supported by research grants: 12301480 from NSFC, 2025A1515010945 from Natural Science Foundation of Guangdong Province.

References↩︎