Algorithms for solving a class of real quasi-symmetric Toeplitz linear systems and its applications

1.   Introduction

There have been extensive studies and applications of Toeplitz and quasi-Toeplitz matrices recently. Toeplitz or quasi-Toeplitz matrices are adopted as a core tool in the quantum magnetic field ^[1], generation of random bits ^[2], quantum key distribution ^[3,4,5], the evaluation of the short-time rate of change in the trend of CO$_2$ data ^[6], the deconvolution of hemodynamic responses ^[7] and the scattering and radiation on thin wire models ^[8]. Some studies have utilized Toeplitz or quasi-Toeplitz matrices to solve the ordinary and partial differential equations ^[9,10], also applying them to solve the fractional differential equations ^{[11,12,13,14,15,16]}.

A large number of linear systems with Toeplitz coefficient matrices need to be solved. Recently, Liu et al. developed fast solvers for Toeplitz linear systems ^[17,18,19]. The inversion of the Toeplitz matrix is represented as a combination of circulant and skew-circulant matrices in ^[20,21] so that the solution can be gained directly by using the fast Fourier transform (FFT) and the inverse FFT (IFFT). Here, we turn to express the Toeplitz inversion by using the sum products of the skew-imaginary circulant ^[22] and skew-circulant matrices, where the skew-imaginary circulant matrix is an $\omega$-circulant matrix with $\omega = -i$, $i = \sqrt{-1}$.

Regarding the design of fast algorithms for a quasi-Toeplitz linear system, the authors of ^[23] presented new efficient algorithms to solve the CUPL-Toeplitz linear system ^[24,25,26] by first splitting the CUPL-Toeplitz matrix into a Toeplitz matrix and a low-rank matrix. In ^[27], Zhang et al. studied order-reduction algorithms to further reduce the computation complexity. The tridiagonal quasi-Toeplitz linear system is solved in ^[28].

We are concerned with a class of quasi-symmetric Toeplitz linear systems. The authors of ^[18] showed the trigonometric transformation splitting method to solve the real symmetric Toeplitz linear system. Different from the iterative method, we adopt the direct method based on the symmetric Toeplitz inversion representation, which is simplified due to the symmetry. We also show the numerical stability of the splitting symmetric Toeplitz inversion inspired by previous studies ^[26,29,30].

In this paper, we study an $n \times n$ real symmetric Toeplitz matrix with perturbations in the first and last columns $P = (p_{j, k})^n_{j, k = 1}\in \mathbb{R}^{n\times n}$:

Obviously,

where $\alpha = [0, \varsigma_1, 0, \dots, 0]^T$, $e_1 = [1, 0, \dots, 0]^T$ is the first unit vector, $\beta = [0, \dots, 0, \varsigma_2, 0]^T$, $e_n = [0, \dots, 0, 1]^T$ is the last unit vector and $A = (t_{|j-k|})^n_{j, k = 1}$ is a nonsingular symmetric Toeplitz matrix.

An outline of this paper is as follows. Section 2 solves a sequence of linear systems with the same symmetric Toeplitz coefficient matrix. Section 3 presents the solvers for the quasi-symmetric Toeplitz linear system. A discussion of the quasi-symmetric Toeplitz matrix-vector multiplication is presented in Section 4. In Section 5, the stability analysis of the splitting inverse of the symmetric Toeplitz matrix is shown. Different numerical examples are given in Section 6. In Section 7, images are encrypted and decrypted by using the proposed quasi-symmetric Toeplitz matrix. Finally, we express our conclusions in Section 8.

2.   Solving symmetric Toeplitz linear systems with the same coefficient matrix

The starting point of the algorithms derived in this paper is the following inversion formula of Toeplitz matrices:

Lemma 1. [31,p. 738] Let $A = (t_{j-k})_{j, k = 1}^{n }\in \mathbb{R}^{n\times n}$ be a Toeplitz matrix. If the vectors $x = [x_1, x_2, \dots, x_n]^T$, $x_1\neq0$ and $y = [y_1, y_2, \dots, y_n])^T$ are the solutions of the linear systems

Then $A$ is invertible and

where $S_{I1}$ and $S_{I2}$ are skew-imaginary circulant matrices ^[22] holding the first columns as $x = [x_1, x_2, \dots, x_n]^{T}$ and $\tilde{y} = [i y_n, y_1, y_2, \dots, y_{n-1}]^{T}$ respectively; $S_1$ and $S_2$ are skew-circulant matrices with the first columns $\hat{y} = [-y_n, y_1, y_2, \dots, y_{n-1}]^{T}$ and $x = [x_1, x_2, \dots, x_n]^{T}$ respectively.

Note that we are concerned with a symmetric Toeplitz matrix $A$; then, the special structure means that the solutions $x$ and $y$ in Eq (2.1) satisfy $y = Jx$, where the matrix $J$ is an anti-identity matrix of order $n$.

If $Ax = e_1$ has the solution $x = [x_1, x_2, \dots, x_n]^{T}$ and $x_1\neq0$, then the symmetric Toeplitz matrix $A$ is invertible and Eq (2.2) can be rewritten as

where $S_I$ and $S$ have the same first column $x = [x_1, x_2, \dots, x_n]^{T}$ and the symbol $S_I^{\ast}$ denotes the conjugate transpose of $S_I$.

Performing the diagonalization scheme on the skew-imaginary circulant matrix, that is, $S_I = \Omega_n ^{\ast}F_n^{\ast}\Lambda_{S_I} F_n\Omega_n$, we can obtain

where $F_n = (F_{j, k})_{j, k = 1}^n$, $F_{j, k} = \frac{1}{\sqrt{n}}e^{\frac{2\pi i(j-1)(k-1)}{n}}$, $1\leq j, k\leq n$, $\Omega_n = \mathrm{diag}(1, e^{\frac{-3i\pi}{2n}}, \dots, e^{\frac{-3i(n-1)\pi}{2n}})$ and $\Lambda_{S_I}$ is a diagonal matrix containing the eigenvalues of $S_I$.

Equation (2.4) shows a new decomposition of the symmetric Toeplitz inversion. By taking full account of the structural characteristics of the symmetric matrix, we solve only one system $Ax = e_1$ instead of two linear systems.

In the case of solving a sequence of linear systems with the same coefficient matrix, it is effective to solve $Ax = e_1$ first and then use the new decomposition of the inverse of the symmetric Toeplitz matrix for calculation. Therefore, Algorithm 1 is proposed for the computation of $x$ and the eigenvalues of $S_I$. By applying Eq (2.4), an order-reduction algorithm ^[27] and FFT and IFFT operations, Algorithm 2 realizes the product of $A^{-1}$ and $v$ in the real number fields.

In the proposed algorithms, the diagonal matrix can be represented by a vector containing the diagonal entries; then, the multiplication of the diagonal matrix and the vector is implemented by using a dot product, thus reducing the storage and computational complexity.

To solve a sequence of linear systems with a constant symmetric Toeplitz coefficient matrix, Algorithm 1 only needs to be computed once. The first step in Algorithm 1 is to solve the symmetric Toeplitz linear system $Ax = e_1$. Any symmetric Toeplitz linear solver can be used, such as the generalized minimal residual (GMRES) method, simplified quasi-minimal residual method, conjugate gradient (CG) method, etc. Since $A$ is set as a symmetric positive definite matrix in our numerical experiments in this paper, we choose a the preconditioned CG (PCG) method with the Strang circulant preconditioner ^[32] to solve it. The workload of Algorithm 1 is considered to be $O(n\log_2 n)$.

One FFT or IFFT needs $5n\log_2 n + O(n)$ real arithmetic operations [33,p. 75], and one run of Algorithm 1 in ^[27] needs $\frac{15n}{2}\log_2 n + O(n)$. Algorithm 2 mainly consists of two FFTs, one IFFT and two runs of Algorithm 1 in ^[27]. In the first step of Algorithm 2, it can save two FFTs with the length $\frac{n}{2}$ when calculating $S^T v$ and $Sv$ simultaneously. Finally, we regard the workload of Algorithm 2 as $25n\log_2 n+ O(n)$. The application of the order-reduction algorithm reduces the entire computational complexity.

A sequence of Toeplitz linear systems with the same coefficient matrix is solved in ^[20,21,34]. The method of using the splitting Toeplitz matrix inversion greatly saves the computational time. This reminds us that in a mathematical or engineering problem, we may be faced with thousands of linear systems with the same symmetric Toeplitz matrix. By applying Algorithm 2, it will be more efficient to solve these symmetric Toeplitz linear systems.

3.   Quasi-symmetric Toeplitz linear solvers

In this section, we solve a real quasi-symmetric Toeplitz linear system

where $a = [a_1, a_2, \dots, a_n]^T$ is an unknown vector and $b = [b_1, b_2, \dots, b_n]^T$ is known. The coefficient matrix has a decomposition shown in Eq (1.2). We substitute Eq (1.2) into Eq (3.1) and multiply both sides of the equation by $A^{-1}$ from the left; given that $I_n$ is an identity matrix, we have

Let $A^{-1} \alpha = \mu$, $A^{-1}\beta = \nu$ and $A^{-1}b = \eta$, where $\mu = [\mu_1, \mu_2, \dots, \mu_n]^T$, $\nu = [\nu_1, \nu_2, \dots, \nu_n]^T$ and $\eta = [\eta_1, \eta_2, \dots, \eta_n]^T$. The solution of Eq (3.2) can be

Therefore, solving $Pa = b$ is the same as calculating Eq (3.3). Based on the Sherman-Morrison-Woodbury (SMW) formula ^[35], $(I_n+\mu e_1^T+\nu e_n^T)^{-1}$ can be represented as

where $C = (c_{j, k})^n_{j, k = 1}$ and

According to Eqs (3.3)–(3.5), we have

To solve the linear system $Pa = b$, we must solve the symmetric Toeplitz systems $A \mu = \alpha$, $A \nu = \beta$ and $A\eta = b$ at first. By using Algorithms 1 and 2 and Eqs (3.5) and (3.6), we obtain the following:

Similarly, if we multiply Eq (1.2) by $A^{-1}$ from the right, we can get $P = (I_n+ \alpha e_1^T A^{-1}+ \beta e_n^T A^{-1})A$ and substitute it into Eq (3.1); then,

Denote $e_1^T A^{-1} = \rho^T$, $e_n^T A^{-1} = \sigma^T$ and $Aa = \psi$, where $\rho = [\rho_1, \rho_2, \dots, \rho_n]^T$, $\sigma = [\sigma_1, \sigma_2, \dots, \sigma_n]^T$ and $\psi = [\psi_1, \psi_2, \dots, \psi_n]^T$. Then, Eq (3.7) can be rewritten as

The SMW formula is applied to Eq (3.8); we can write

where $D = (d_{j, k})^n_{j, k = 1}$ and

From Eqs (3.8)–(3.10), we can get

Since $A$ is a symmetric Toeplitz matrix, if $e_1^T A^{-1} = \rho^T$ and $e_n^T A^{-1} = \sigma^T$, then $A\rho = e_1$ and $A\sigma = e_n$. According to Eqs (3.10) and (3.11) and the proposed algorithms in Section 2, we propose another algorithm for solving $Pa = b$.

From the above analysis, it is quite evident that the main work of Algorithms 3 and 4 is solving three symmetric Toeplitz linear systems. Moreover, we can omit one calculation of $A\rho = e_1$ in Algorithm 4.

A greater number of symmetric Toeplitz linear systems will need to be solved if there are more perturbations in the coefficient matrix of the quasi-symmetric Toeplitz linear system. It is therefore meaningful to decompose the inverse of the symmetric Toeplitz matrix, as it will save computational time, especially in the case of high-order systems.

4.   Fast algorithm for the product of the quasi-symmetric Toeplitz matrix and vector

In terms of quasi-symmetric Toeplitz matrix-vector multiplication, the symmetric Toeplitz matrix $A$ can be seen as the sum of the circulant and skew-circulant matrices ^[36], that is,

where $v = (v_1, v_2, \dots, v_n)$ is the vector and the circulant matrix $C$ and the skew-circulant matrix $S$ are all symmetric.

The first columns of circulant and skew-circulant matrices need to be computed in the splitting step. Considering the symmetry of $n$-th order $A$, the first columns of $C$ and $S$ can be obtained by applying the first $\frac{n}{2}+1$ ($n$ is even) or $\frac{n+1}{2}$ ($n$ is odd) elements instead of $n$ elements, which can reduce the computation complexity.

Based on the diagonalization scheme of the circulant matrix, $C$ has the spectral decomposition $C = F_n^{\ast}\Delta_C F_n$, where $F_n = (F_{j, k})_{j, k = 1}^n$, $F_{j, k} = \frac{1}{\sqrt{n}}e^{\frac{2\pi i(j-1)(k-1)}{n}}$, $1\leq j, k\leq n$, and $\Delta_C$ is a diagonal matrix containing the eigenvalues of $C$. In ^[27], an order-reduction algorithm for the multiplication of the real skew-circulant matrix and vector is proposed. So the product of the quasi-symmetric Toeplitz matrix and vector can be expressed as

By applying Eq (4.2), the order-reduction algorithm for real skew-circulant matrix-vector multiplication, the FFT and IFFT operations, we give Algorithm 5 for $Pv$ in the real number field as follows.

Because the complexity of one FFT or IFFT is $5n\log_2 n + O(n)$ [33,p. 75], Algorithm 1 in ^[27] needs $\frac{15n}{2}\log_2 n + O(n)$ real arithmetic operations. The complexity of Algorithm 5 is $\frac{45n}{2}\log_2 n + O(n)$, and it consists of one order-reduction algorithm, two FFTs and one IFFT.

5.   Stability analysis

We find that the inverse factorization of the symmetric Toeplitz matrix is critical for the solution of quasi-symmetric Toeplitz linear equations. The following is the error analysis of Eq (2.3) in terms of the 1-norm, $\infty$-norm and 2-norm, respectively. Assume that $\hat{x} = [\hat{x}_1, \hat{x}_2, \dots, \hat{x}_n]^{T}$ is the numerical solution of $Ax = e_1$. If $\hat{x}_1 \neq0$, we denote

as a perturbation of $A^{-1}$, where the forms of $\hat{S}_I$ and $\hat{S}$ are the same as those in Eq (2.3) and have corresponding perturbations.

Theorem 1. Let $\epsilon > 0$; if $x_1 \neq0$ and $\hat{x}_1 \neq0$, suppose that the relative error of $\frac{1}{x_1}$ is $\hat{\epsilon} = \frac{|1/x_1-1/\hat{x}_1|}{|1/x_1|}$ and

we get

and

Proof. From the representation of $A^{-1}$ and $\hat{A}^{-1}$ and Eqs (2.3) and (5.1), we can get

On the one hand, we can get

According to the structural characteristics of $S_I$ and $S$ and Eq (5.2), we note that

$\hat{\epsilon} = \frac{|1/x_1-1/\hat{x}_1|}{|1/x_1|}$ is the relative error; then,

Combining Eqs (5.6)–(5.8), we thus obtain

Similarly, on the other hand,

Based on Eqs (5.5), (5.9) and (5.10), we can write

By $Ax = e_1$, it is proven that $\Vert x \Vert_1 \leq \Vert A^{-1} \Vert_1$, so

Theorem 2. Let $\epsilon > 0$; if $x_1 \neq0$, $\hat{x}_1 \neq0$ and

then

and

where $\hat{\epsilon} = \frac{|1/x_1-1/\hat{x}_1|}{|1/x_1|}$ is the relative error of $\frac{1}{x_1}$.

Proof. Similar to Theorem 1, this proof is based on the same conditions:

and

Theorem 3. Under the assumptions and concepts of the previous two theorems, the upper bound of the $2$-norm is said to be

and

Proof. As we all know, $\Vert A^{-1}-\hat{A}^{-1} \Vert_2 ^2 \leq \Vert A^{-1}-\hat{A}^{-1} \Vert_1 \cdot \Vert A^{-1}-\hat{A}^{-1} \Vert_\infty$; from Eqs (5.3) and (5.12), we have

Recall that

From Eqs (5.14) and (5.16), we can get

We show the stability analysis of inverse factorization of the symmetric Toeplitz matrix in this section. The absolute and relative perturbation upper bounds are shown in Eqs (5.3) and (5.4), (5.12) and (5.13), (5.14) and (5.15), as far as $Ax = e_1$ is solvable.

6.   Numerical simulations

In this section, we give some examples for the comparison of different algorithms. The first two examples present a comparison of different quasi-symmetric Toeplitz linear solvers. Examples 3 and 4 are devoted to the fast algorithm for the quasi-symmetric Toeplitz matrix-vector multiplication.

These experiments were done by using MATLAB (R2022a) on a laptop with the following specifications: 16 GB RAM, AMD Ryzen 7 5800H CPU 3.20 GHz. In the following tables, the calculation time is in seconds, "$n$" denotes the matrix order and "—" indicates that it is out of MATLAB's memory.

Example 1. An $n \times n$ quasi-symmetric Toeplitz linear system $Pa = b$ is considered in this example. According to Eq (1.2), the first column of $A$ is $(t_{i})^n_{i = 1} = \frac{1}{i}$, and $\varsigma_1$ and $\varsigma_2$ are random values in the range of $(0, 1)$. Assume that $a_{\text{exact}} = [1, 1, ..., 1]^T$ is the exact solution, so the vector $b$ is $b = Pa_{\text{exact}}$.

A comparison of the error and time required to solve a quasi-symmetric Toeplitz linear system is presented in Table 1. The derived two algorithms for solving $Pa = b$ are shown in the table as Algorithms 3 and 4. Based on Algorithm 3, Alg$_{\text{Huang}}$ Ⅰ replaces Steps 1 and 2 with Huang's algorithm ^[20,21] for solving three symmetric Toeplitz systems. Similarly, Alg$_{\text{Huang}}$ Ⅱ was executed based on Algorithm 4. The back-slash method means solving $Pa = b$ directly by using the back-slash operator in MATLAB. $\text{"Error"}$ was set to be $\| \frac{a-a_{\text{exact}}}{a_{\text{exact}}}\|_\infty$, and the stopping criterion for those algorithms, besides Back-slash, was set to be $1\times10^{-7}$.

Algorithms 3 and 4 have significate superiority over the other methods. They can solve high-order quasi-symmetric Toeplitz linear systems. However, Back-slash cannot work when the matrix order is above $2^{14}$. The proposed Algorithms 3 and 4 have the shortest computational time. They perform better than Alg$_{\text{Huang}}$ Ⅰ and Alg$_{\text{Huang}}$ Ⅱ in high-order linear systems due to the application of the order-reduction algorithm.

Example 2. Consider another quasi-symmetric Toeplitz linear system $Pa = b$, the first column of $A$ and that $\varsigma_1$ and $\varsigma_2$ in Eq (1.2) are generated from an open interval $(0, 1)$. The sum of the elements in the first column was added to the diagonal entries of $A$ to keep $A$ as a diagonally dominant matrix. The exact solution of the system was set to be $a_{\text{exact}} = [1, 1, ..., 1]^T$.

Table 2 shows a comparison of different methods for solving $Pa = b$ for Example 2. Alg$_{\text{PGMRES}}$ refers to solving the linear system by using the GMRES method with the preconditioner. Since $P$ is a nonsymmetric matrix, the GMRES method can be utilized to solve $Pa = b$. In order to accelerate the convergence rate, the Strang circulant preconditioner ^[32] of the symmetric Toeplitz matrix $A$ was established.

According to theoretical analysis, the complexity of Algorithm 4 is lower than that of Algorithm 3. The results show that Algorithm 3 takes more time than Algorithm 4, which is consistent with the theoretical analysis. Alg$_{\text{PGMRES}}$ performs better than the other methods. However, compared with the algorithm for solving symmetric systems, the convergence analysis of the GMRES algorithm is very difficult, and there is no clear conclusion at present.

Example 3. Consider the multiplication of $Pv$. The vector $v$ is Gaussian distributed. The quasi-symmetric Toeplitz matrix $P$ is decomposed as Eq (1.2), where the first column of $A$ is $(t_{i})^n_{i = 1} = \frac{1}{i}$ and $\varsigma_1$ and $\varsigma_2$ are random values in the range $(0, 1)$.

In Table 3, $Pv$-direct denotes calculation of the quasi-symmetric Toeplitz matrix-vector multiplication via MATLAB. $Pv$-splitting refers to splitting $P$ according to Eq (4.1) first, and then utilizing the diagonalization scheme of the circulant matrix and skew-circulant matrix for fast calculation. Algorithm 5 is the proposed quasi-symmetric Toeplitz matrix-vector multiplication. Obviously, Algorithm 5 consumed the least amount of time.

Example 4. Consider the multiplication of $Pv$, the first column of $A$, $\varsigma_1$ and $\varsigma_2$ in Eq (1.2), and that the vector $v$ are random values in the range $(0, 1)$. The sum of the elements in the first column was added to the diagonal entries of $A$.

Table 4 shows another example for the quasi-symmetric Toeplitz matrix-vector multiplication. As the matrix order increases, the efficiency of Algorithm 5 becomes increases. Our proposed algorithm is suitable for high-order matrix-vector multiplication. However, the $Pv$-direct method not only takes the longest time, but it also cannot work when the matrix order is above $2^{14}$.

7.   Applications

In the past decades, image encryption and decryption have been widely researched in the area of information security. The results show that the proposed algorithms can be used to encrypt and decrypt images.

Example 1. Based on Eq (1.2), we consider a quasi-symmetric Toeplitz matrix $P$, where the entries in the first column of $A$ and $\varsigma_1$ and $\varsigma_2$ are random values in the range $(0, 1)$. To keep $P$ as a diagonally dominant matrix, we added a parameter to the diagonal entries of $P$. Images were encrypted and decrypted by left-multiplying by $P^2$ and $P^{-2}$, respectively.

Figures 1–3 present the effects of image encryption and decryption. The famous figure "Lenna" is taken as an example. Let the original image matrix be $X = [x_1, x_2, \dots, x_{256}]$; the encrypted image matrix can be obtained by matrix-vector multiplication via Algorithm 5, that is, $Y = [y_1, y_2, \dots, y_{256}] = P(P[x_1, x_2, \dots, x_{256}])$. The decrypted image matrix can be computed by $\hat{X} = [\hat{x}_1, \hat{x}_2, \dots, \hat{x}_{256}] = P^{-1}(P^{-1}[y_1, y_2, \dots, y_{256}])$. Algorithm 4 was performed for image decryption because of its lower complexity than Algorithm 3.

In the process of image encryption and decryption, it is necessary to do multiple matrix-vector multiplications with the same quasi-symmetric Toeplitz matrix and solve linear systems with the same coefficient matrix. Some calculations can be done only once, like Steps 1 and 2 of Algorithm 5 in image encryption, as well as Steps 1 and 2 of Algorithm 4 in image decryption. So, the computational time will be saved.

It is obvious that our proposed algorithms can realize efficient encryption and decryption for different pixel-sized images.

8.   Conclusions

Algorithms for solving a class of real quasi-symmetric Toeplitz linear systems are shown in this paper. We have proposed methods for solving a sequence of linear systems with the constant symmetric Toeplitz matrix based on the decomposition of Toeplitz matrix inversion. By splitting the coefficient matrix into a symmetric Toeplitz matrix plus two low-rank matrices and combining it with the SMW formula, two fast algorithms with $O(n\log_2 n)$ complexity have been given to solve the real quasi-symmetric Toeplitz linear system. Quasi-symmetric Toeplitz matrix-vector multiplication in the real number field has been presented. Structural perturbation analysis of the inverse factorization of symmetric Toeplitz was also analyzed. Numerical simulations and applications have shown that the proposed algorithms are accurate and efficient.

Acknowledgements

The research was supported by the National Natural Science Foundation of China (Grant No. 12101284) and Natural Science Foundation of Shandong Province (Grant No. ZR2020MA051).

Conflict of interest

The authors declare that there is no conflict of interest.