Preconditioned TBiCOR and TCORS Algorithms for Solving the Sylvester Tensor Equation

In this paper, the preconditioned TBiCOR and TCORS methods are presented for solving the Sylvester tensor equation. A tensor Lanczos $\mathcal{L}$-Biorthogonalization algorithm (TLB) is derived for solving the Sylvester tensor equation. Two improved TLB methods are presented. One is the biconjugate $\mathcal{L}$-orthogonal residual algorithm in tensor form (TBiCOR), which implements the $LU$ decomposition for the triangular coefficient matrix derived by the TLB method. The other is the conjugate $\mathcal{L}$-orthogonal residual squared algorithm in tensor form (TCORS), which introduces a square operator to the residual of the TBiCOR algorithm. A preconditioner based on the nearest Kronecker product is used to accelerate the TBiCOR and TCORS algorithms, and we obtain the preconditioned TBiCOR algorithm (PTBiCOR) and preconditioned TCORS algorithm (PTCORS). The proposed algorithms are proved to be convergent within finite steps of iteration without roundoff errors. Several examples illustrate that the preconditioned TBiCOR and TCORS algorithms present excellent convergence.


Introduction
This paper is concerned of the computation of the Sylvester tensor equation of the form where matrices A n ∈ R In×In (n = 1, 2, ..., N ) and the tensor D ∈ R I1×I2×...×IN are given, and the tensor X ∈ R I1×I2×...×IN is unknown. The Sylvester tensor equation (1.1) plays vital roles in many fields such as image processing [8], blind source separation [22] and the situation when we describes a chain of spin particles [2]. If N = 2, then (1.1) can be reduced to the Sylvester matrix equation which has many applications in system and control theory [12][13][14]. When N = 3, (1.1) becomes which often arises from the finite element [20], finite difference [3] and spectral methods [21]. Many approaches are constructed to solve the Sylvester tensor equation (1.1) in recent years. Chen and Lu [9] proposed the GMRES method based on a tensor format for solving (1.1) and presented the gradient based iterative algorithms [10] for solving (1.3). Beik et al. [24] also presented some global iterative schemes based on Hessenberg process to solve the Sylvester tensor equation (1.3). Beik et al. [5] solved the Sylvester tensor equation (1.1) with severely ill-conditioned coefficient matrices and considered its application in color image restoration. Heyouni et al. [15] proposed a general framework by using tensor Krylov projection techniques to solve high order the Sylvester tensor equation (1.1). Beik et al. [1] proposed the Arnoldi process and full orthogonalization method in tensor form, and the conjugate gradient and nested conjugate gradient algorithms in tensor form to solve the Sylvester tensor equation (1.1). When D in (1.1) is a tensor with low rank, Bentbib et al. [6] proposed Arnoldi-based block and global methods. Kressner and Tobler [19] developed Krylov subspace methods based on extended Arnoldi process for solving the system of equation (1.1). The perturbation bounds and backward error are presented in [27] for solving (1.1). For more methods on other linear systems in tensor form we refer to [4,16,18,23]. Using the nearest Kronecker product (NKP) in [28], Chen and Lu [9] presented an efficient preconditioner for solving Eq.(1.1) based on GMRES in tensor form. Very recently Zhang and Wang in [30] gave a preconditioned BiCG (PBiCG) and a preconditioned BiCR (PBiCR) based on the nearest Kronecker product (NKP) in [28].
Inspired by the Lanczos biorthogonalization (LB) algorithm in [26], BICOR and CORS methods in [11] for non-symmetric linear equation, in this paper, we present two improved Lanczos L-orthogonal algorithms in tensor form for solving the Sylvester tensor equation (1.1). We further present preconditioned TBiCOR (PTLB) and preconditioned TCORS (TCORS) algorithms by using the NKP preconditioner in [9] for solving Eq (1.1). The preconditioned LB in tensor form (PTLB) is also considered.
The rest of this paper is organized as follows. Section 2 reviews some related symbols, concepts and lemmas that will be used in the contexture. Section 3 presents a tensor Lanczos L-biorthogonalization algorithm (TLB) and two improved TLB methods are shown in section 4. The tensor biconjugate L-orthogonal residual(TBiCOR) and tensor conjugate L-orthogonal residual squared(TCORS) algorithms for solving the tensor equation (1.1) are presented in subsections 4.1 and 4.2, respectively. The convergence of the TBiCOR and TCORS methods are proved. Section 5 presents the preconditioned TLB, TBi-COR and TCORS algorithms and the convergence of the preconditioned TBi-COR and TCORS algorithms. Section 6 presents several examples and some conclusions are drawn in section 7.

Preliminaries
The notations and definitions as follows are needed. For a positive integer N , an N -way or N th-order tensor X = (x i1 ... iN ) is a multidimensional array with I 1 I 2 ...I N entries, where 1 ≤ i j ≤ I j , j = 1, ..., N . R I1×...×IN denotes the set of the N th-order I 1 × ... × I N dimension tensors over the real field R, while C I1×...×IN defines the set of the N th-order I 1 × ... × I N dimension tensors over the complex field C.
Let X × n A define the n-mode (matrix) product of a tensor X ∈ R I1×I2×...×IN with a matrix A ∈ R J×In , i.e., The n-mode (vector) product of a tensor X ∈ R I1×I2×...×IN and a vector v ∈ R In is denoted by X × n v, i.e., The inner product of X , Y ∈ R I1×I2×...×IN is defined by x i1i2...iN y i1i2...iN , and the norm of X is denoted by Furthermore, it derives from [9] that We refer more notions and definitions in [17].
The following results from [1] will be used later.
Lemma 1. Suppose A ∈ R Jn×In , y ∈ R Jn and X ∈ R I1×I2×...×IN , we have X × n A× n e = X × n (A T y).
where e j is the j-th column of the I N -order identity matrix E IN , and X j denotes the j-th frontal slice of X .
and (2.4) can be represented as We also need the following results. (2.8)

A Tensor Lanczos L-Biorthogonalization Algorithm
Define the linear operator of the form Let L T define the dual linear operator of L, i.e., then it holds that L(X ), Y = X , L T (Y) for any X , Y ∈ R I1×I2×...×IN . Define the Krylov subspaces in tensor form as follows: Algorithm 1 lists the Lanczos L-Biorthogonalization procedure in tensor form that will be used to produce two series of biorthogonalization tensors.
Algorithm 1 A Lanczos L-biorthogonalization procedure in tensor form.
We have the following results for Algorithm 1. The proofs of these results are similar to the proof of Proposition 1 in [11] by using the definitions of the inner product (2.4), (2.5) and linear operator L in (3.1) and are omitted. Proposition 1. If Alogrithm 1 stops at the m-th step, then the tensors V j and W i (i, j = 1, 2, ..., m) produced Algorithm 1 are L-biorthogonal, i.e., Proposition 2. Suppose that V m is the (N + 1)-order tensor with columns V 1 , V 2 , ..., V m , and W m is the (N +1)-order tensor with the columns W 1 , W 2 , ..., W m , H m and G m are the (N + 1)-order tensors with the columns H j := L(V j ) and G j := L T (W j ) (j = 1, 2, ..., m), respectively. Then we have and being a triangular matrix with its elements generated by Algorithm 1. Moreover, it holds that and where E m denotes the identity matrix with m order.
We remark that (3.8) and (3.9) can be represented as and

1). For any initial tensor
It is easy to verify that a series of tensors {V 1 , V 2 , ..., V m } produced via Algorithm 1 form a basis of K m (L, V 1 ). Thus we have where y m ∈ R m . By Eq (3.17) and (3.18), we have A further computation results in Through a simple inner product operation and according to Eq.(2.6) and Lemma Submitting Eq.(3.13) into Eq.(3.21) results in the tridiagonal system on y m : Once we compute y m by (3.22), we get the solution X m of (1.1) by (3.18). We summarize this method in Algorithm 2, which is called Tensor Lanczos L-Biorthogonalization algorithm (TLB).
We remark Algorithm 2 have to compute the inverse of T m . When T m is of large size, it needs much computation. We present two improved algorithms for Algorithm 2 in the next section.

The TBiCOR Algorithm
In this subsection, we develop an improved algorithm by introducing the LU decomposition to T m in Algorithm 2.
Let the LU decomposition of T m be  .22), the solution of the dual system L T (X * ) = D * can be represented as .
Proposition 4. Let P i and P * i (i = 1, ..., k) are the i-th column tensor of P k and P * k , respectively. It holds that Proof. According to Lemma 4 and (3.13), we have which implies that (4.8) holds.
For a given initial guess X 0 , let R 0 = D − L(X 0 ) and P 0 = R 0 . Set Similarly, for the dual linear system L T (X * ) = D * , we set Now we determine α j and β j in (4.11)-(4.14). According to (4.11) we have that then by Propositions 3, 4 and (4.11)-(4.14) it holds that . (4.16) Similarly according to we have Algorithm 3 summarizes the biconjugate L-orthogonal residual algorithm in Tensor form for solving (1.1), which is abbreviated as TBiCOR.
We have the following convergence properties on Algorithm 3.

Algorithm 3
TBiCOR: A tensor biconjugate L-orthogonal residual algorithm for solving (1.1) According to Proposition 3, we have

The TCORS Algorithm
This subsection presents an improved method on Algorithm 3 by introducing a squared operator of the residual of X n produced by Algorithm 3. The proposed method is called the conjugate L-orthogonal residual squared algorithm in tensor form, which is abbreviated as TCORS.
Algorithm 4 TCORS: A tensor conjugate L-orthogonal residual squared algorithm for solving (1.1) The residual tensor of X n produced by Algorithm 3 can be represented as R n = a 0 R 0 + a 1 L(R 0 ) + a 2 L 2 (R 0 ) + ... + a n L n (R 0 ) = (a 0 L 0 + a 1 L + a 2 L 2 + ... + a n L n )R 0 , where a i is determined by Algorithm 3. Denote ϕ n (L) = a 0 L 0 + a 1 L + a 2 L 2 + ... + a n L n , then (4.19) can be represented as Similarly, we have where φ n (L) = b 0 L 0 +b 1 L+b 2 L 2 +...+b n L n , and b i can be derived by Algorithm 3. For the directions R * n and P * n in Algorithm (3), replacing L in (4.20) and (4.21) with L T results in Thus α n in (4.16) and β n in (4.18) can be represented as According to (4.11)-(4.12), ϕ j and φ j can be expressed as respectively. Squaring on both sides of (4.24) and (4.25) results in Furthermore, we have Taking (4.28) into (4.26) results in then (4.34)-(4.36) can be represented as  Proof. The proof of Theorem 2 is similar to that of Theorem 1 by replacing R k with U k , thus is omitted.

Preconditioned BiCOR and TCORSs Algorithms
This section presents two preconditioned methods based on Algorithms 3-4 for solving Eq.(5.1).
Using the definition of the Kronecker product in [4], one can transform Eq.(1.1) to its equivalent linear system . We refer to [17] for more details.
Algorithm 5 PTBiCOR: A preconditioned tensor biconjugateL-orthogonal residual algorithm for solving (1.1)  [28], Chen and Lu [9] presented an efficient preconditioner for solving Eq.(5.1) based on GMRES in tensor form, which is abbre-viated as preconditioned GMRES (PGMRES) later. Zhang and Wang in [30] gave a preconditioned BiCG (PBiCG) and a preconditioned BiCR (PBiCR) based on NKP in [28]. The preconditioner based on NKP approximates where the optimal parameters a ij in (5.3) can be computed by using the nonlinear optimization software, such as fminsearch in MATLAB.
Proof. The proof is very similar to those of Propositions (3) and (4) with the operator L being replaced byL in Algorithm 5, and is omitted. Proof. The proof is similar to that of Theorem 1 by replacing L withL and is omitted.
We can also obtain a preconditioned TLB (PTLB) by introducing the NKP preconditioner in [9] to Algorithm 2, which is listed in Algorithm 7.

Numerical Experiments
In this section, we show several numerical examples to illustrate Algorithms 2-7 and compare them with CGLS in [16], MCG in [23], preconditioned GMRES (PGMRES) in [9], preconditioned BiCG (PBiCG) and preconditioned BiCR (PBiCR) in [30]. All experiments are implemented on a computer with macOS Big Sur 11.1 and 8G memory. The MATLAB R2018a (9.4.0) is used to run all examples. All algorithms are stopped when the relative error r k = X k − X * / X * < 10 −10 , where X * is assumed to be an exact solution of (1.1). Example 6.1. In this example, we consider the Poisson equation in d-dimensional space [4] −△u = f, in Ω = (0, 1) d , u = 0, on ∂Ω.
A finite difference discretization leads to the Sylvester tensor equation (1.1), with the mesh-width h = 1 11 . We set d = 3 and let the initial tensor X 0 = O. The right-hand side D of (1.1) is constructed by (1.1) with the exact solution X * of (1.3) derived by the MATLAB command tenones(10, 10, 10) in [7]. Algorithms 2-6 are used to solve (1.1) with the matrices A i in (6.1). These methods are compared with CGLS in [16] and MCG in [23], respectively. Figure 1 shows the convergence of the relative error r k versus the number of iterations for all methods. From Figure 1, we can see that our preconditioned Algorithms 5-7 present better convergence than Algorithms 2-4 without preconditioning, PGMRES in [9], PBiCG and PBiCR in [30]. While Algorithms 2-4 can compare with PGMRES, PBiCG and PBiCR, and are better than CGLS [16] and MCG [23]. Algorithm 6 converges fastest among all algorithms. Algorithm 5 converges the second fastest among all algorithms.
A standard finite difference discretization on equidistant nodes combined with the second order convergent scheme [18,20] for the convection term leads to the linear system (1.1) with 2) where n = 1, 2, ..., N and the mesh-size h = 1 p+1 . We consider the case when N = 3 and p = 10. The right-hand side D is constructed by (1.1) with the exact solution X * of (1.3) produced by the MATLAB commend tenones(10, 10, 10) in [7]. Let the initial solution X 0 be a tensor with each element being zero. Algorithms 2-7 are used to solve (1.1) with A i given in (6.2). These methods are compared with CGLS [16], MCG [23], PGMRES in [9], PBiCG and PBiCR in [30]. Table 1 displays the running time, total iteration number (TIN) and relative error of different method with different parameters v = 1, 0.1, 0.001 and c i . Figure 2 shows the convergence of the relative error r k for each method with the parameters v = 0.01, c 1 = 1, c 2 = 2 and c 3 = 3. Table 1 shows that, when the stop criterion is satisfied, preconditioned Algorithms 5-6 require less CPU time and iterations than Algorithms 2-4, PGMRES, PBiCG and PBiCR. In most cases Algorithms 2-4 requires much less CPU time but more iterations than PGMRES, PBiCG and PBiCR, and are better than CGLS [16] and MCG [23] both in CPU time and the number of iterations. Algorithm 6 requires the minimal CPU time and iterations among all methods. Figure 2 shows similar results to that in Figure 1. on the unit square [0, 1]× [0, 1] with homogeneous Dirichlet boundary conditions. We use the MATLAB command fdm 2d matrix in the Lyapack package [25] to generate matrices A i : A i = f dm 2d matrix(I i , e xy , sin(xy), y 2 − x 2 ), ,
The initial solution X 0 is selected as zero tensor. Algorithms 2-6 are used to solve (1.1) with A i given in (6.4). These methods are compared with CGLS [16], MCG [23], PGMRES in [9], PBiCG and PBiCR in [30]. Figure 3 shows that Algorithm 6 converges fastest among all methods and Algorithm 6 converges the second fastest among all methods, which are very similar to those in Figures 1 and 2.

Conclusion
This paper first presents a tensor Lanczos L-Biorthogonalization (TLB) algorithm for solving the Sylvester tensor equation (1.1) based on the Lanczos L-Biorthogonalization procedure. Then two improved methods based on the TLB algorithm are developed. The one is the biconjugate L-orthogonal residual algorithm in tensor form (TBiCOR). The other is the conjugate L-orthogonal residual squared algorithm in tensor form (TCORS). The preconditioner based on the nearest Kronecker product (NKP) are used to accelerate the TBiCOR and TCORS algorithms, thus we present preconditioned a preconditioned TBi-COR method and a preconditioned TCORS method. The convergence of these proposed algorithms are proved. Numerical examples show the advantage of the preconditioned TBiCOR and TCORS methods. 8 Acknowledgments