A Gradient Based Iterative Solutions for Sylvester Tensor Equations

Abstract

This paper is concerned with the numerical solution of the Sylvester tensor equation, which includes the Sylvester matrix equation as special case. By applying hierarchical identification principle proposed by Ding and Chen, 2005, and by using tensor arithmetic concepts, an iterative algorithm and its modification are established to solve the Sylvester tensor equation. Convergence analysis indicates that the iterative solutions always converge to the exact solution for arbitrary initial value. Finally, some examples are provided to show that the proposed algorithms are effective.

1. Introduction

This paper is concerned with the solution of the following equation: 𝒳×₁𝐴+𝒳×₂𝐵+𝒳×₃𝐶=𝒟,(1) where matrices 𝐴∈𝑅^𝑚×𝑚, 𝐵∈𝑅^𝑛×𝑛, and 𝐶∈𝑅^𝑙×𝑙, and tensor 𝒟∈𝑅^{𝑚×𝑛×𝑙}. The properties of the operators ×₁, ×₂, and ×₃ are described in detail in Section 2. When 𝒳 is a 2-mode tensor, that is, a matrix, (1) reduces to 𝐴𝑋+𝑋𝐵^𝑇=𝐷.(2) It is just a Sylvester matrix equation, which has received much attention, often motivated by applications in control theory, model reduction, signal processing, filtering, and system identification [1–7]. So we call (1) a Sylvester tensor equation in this paper.

The simplest direct approach for solving the matrix equation (2) is Kronecker product method [8] in which the huge memory is needed, for example, if 𝑚=𝑛=100, the dimension of the new coefficient matrix is 10000. Alternative direct method is to firstly transform coefficient matrices into special forms, such as Schur canonical form [9], Hessenberg-Schur form [1], and [10–12], and then to solve another matrix equation which may be readily computed. Such kinds of methods have been widely adopted. For example, 3D Schur decomposition was developed to solve the Sylvester tensor equation (1) (for more details see [13]). However, these methods require the transformations of coefficient matrices, which need considerable cost.

Iterative algorithms for solving matrix equations have been studied extensively, for example, [14–21]. For example, Ding and Chen proposed a few simple iterative methods for matrix equation in [22–24]. The methods, resembling the classical Jacobi and Gaussian iterations for linear systems, are easy to implement and cost little per step. Gradient based iterative algorithms and least squares based iterative algorithms have also been proposed in [25–27] for (coupled) matrix equations, which include the Sylvester matrix equations and Lyapunov matrix equations as a special form. However, to the best knowledge of the authors, there is hardly any iterative method for the Sylvester tensor equation in references up to now.

In this paper, we first extend the gradient based iterative algorithm for Sylvester matrix equation (2) to solve the Sylvester matrix equation (1) based on previous work in [22]. The basic idea of such an algorithm is to regard the unknown tensor as the system parameter tensor by applying hierarchical identification principle. Then, we present a more efficient modified iterative method. The convergence properties of the two algorithms are proved. Numerical examples are given to show that the two iterations are effective.

The rest of this paper is organized as follows. In Section 2, we describe tensor notations and common operations used throughout the paper. An extension and modification of the gradient based iterative algorithms to the Sylvester tensor equation (1) are presented in Section 3. The convergence of the algorithms are analyzed in Section 4. Section 5 contains some numerical experiments to illustrate the theoretical results in this paper. Finally, some conclusions are outlined in Section 6.

2. Tensor Notations and Basic Operations

In this section, we briefly review some concepts and notations that are used throughout the paper. Matrices are denoted by capital letter, for example, 𝐴, 𝐵, and 𝐶. Tensors are denoted by Euler script letter, for example, 𝒳,𝒴, and 𝒵. The 𝑛×𝑛 identity matrix is denoted by 𝐼_𝑛. The order of a tensor is the number of dimensions, also known as ways or modes. As a special case, a vector is 1-mode tensor and a matrix is 2-mode tensor. The operator vec(⋅) stacks the columns of a matrix to form a vector. ‖⋅‖ denotes the Frobenius norm unless otherwise stated.

2.1. Tensor-Matrix Multiplication

An important operation for a tensor is the tensor-matrix multiplication [28–30]. The 1-mode product of a tensor 𝒳∈𝑅^{𝐼×𝐽×𝐾} by a matrix 𝐴∈𝑅^𝑃×𝐼, denoted by 𝒳×₁𝐴, is a 𝑃×𝐽×𝐾 tensor in which the entries are given by 𝒳×₁𝐴𝑝,𝑗,𝑘=^𝐼_𝑖=1𝑥_𝑖𝑗𝑘𝑎_𝑝𝑖.(3) Similarly, 2-mode multiplication of a tensor by a matrix 𝐵∈𝑅^𝑄×𝐽𝒳×₂𝐵𝑖,𝑞,𝑘=^𝐽_𝑗=1𝑥_𝑖𝑗𝑘𝑏_𝑞𝑗.(4) 3-mode multiplication is analogous and is omitted here.

It follows immediately from the definition that 𝑛-mode and 𝑚-mode multiplication commutes if 𝑚≠𝑛: 𝒳×_𝑚𝐴×_𝑛𝐵=𝒳×_𝑛𝐵×_𝑚𝐴.(5) If the modes are the same, then 𝒳×_𝑛𝐴×_𝑛𝐵=𝒳×_𝑛(𝐵𝐴).(6)

2.2. Inner Product and Tensor Norm

The inner product of two tensors 𝒳,𝒴∈𝑅^{𝐼×𝐽×𝐾} is defined by ⟨𝒳,𝒴⟩=^𝐼𝐽_𝑖=1^𝐾_𝑗=1𝑘=1𝑥_𝑖𝑗𝑘𝑦_𝑖𝑗𝑘,(7) and the norm induced by this inner product is ‖𝒳‖=⟨𝒳,𝒳⟩.(8)

For two tensors 𝒳 and 𝒴, the computation of ‖𝒳−𝒴‖ can be simplified as follows: ‖𝒳−𝒴‖²=‖𝒳‖²−2⟨𝒳,𝒴⟩+‖𝒴‖².(9)

Moreover, 𝑛-mode multiplication commutes with respect to the inner product, that is, 𝒳,𝒴×_𝑛𝐴=𝒳×_𝑛𝐴^𝑇,𝒴.(10)

Finally, the tensor norm of 𝑛-mode multiplication is mutually consistent, that is, 𝒳×_𝑛𝐴≤‖𝐴‖₂‖𝒳‖.(11)

3. Iterative Algorithms

In this section, we apply the tensor arithmetic concepts and the hierarchical identification principle to solve the Sylvester tensor equation (1).

Firstly, we consider the existence and uniqueness of the solution of (1). By using the Kronecker product, we have 𝐼_𝑙⊗𝐼_𝑛⊗𝐴+𝐼_𝑙⊗𝐵⊗𝐼_𝑚+𝐶⊗𝐼_𝑛⊗𝐼_𝑚(𝒳)(𝒟),vec=vec(12) which is equivalent to (1), This vectorized representation immediately leads to the following result.

Lemma 1. Equation (1) has a unique solution if and only if 𝜆_𝑖(𝐴)+𝜆_𝑗(𝐵)+𝜆_𝑘(𝐶)≠0  for any 𝑖, 𝑗, and 𝑘, where 𝜆_𝑖(𝐴), 𝜆_𝑗(𝐵) and 𝜆_𝑘(𝐶) represent the eigenvalues of 𝐴,𝐵, and 𝐶, respectively.

In the case of the Sylvester matrix equation (2), this corresponds to the well-known condition 𝜆_𝑖(𝐴)+𝜆_𝑗(𝐵)≠0, for any 𝑖 and 𝑗.

In order to drive the iterative algorithms for solving (1), we need to introduce the following residual tensors: ℰ₁∶=𝒟−𝒳×₂𝐵−𝒳×₃ℰ𝐶,₂∶=𝒟−𝒳×₃𝐶−𝒳×₁ℰ𝐴,₃∶=𝒟−𝒳×₁𝐴−𝒳×₂𝐵.(13) According to the hierarchical identification principle [22], (1) is decomposed into three subsystems: 𝒳×₁𝐴=ℰ₁,𝒳×₂𝐵=ℰ₂,𝒳×₃𝐶=ℰ₃.(14) Consider the following least squares problem: 𝐽₁(𝒳)=12𝒳×₁𝐴−ℰ₁²,𝐽₂(𝒳)=12𝒳×₂𝐵−ℰ₂²,𝐽₃(𝒳)=12𝒳×₃𝐶−ℰ₃².(15) Note that the derivatives of 𝐽_𝑖(𝑖=1,2,3) with respect to 𝒳 are derived in [30] and are given by 𝐽₁(𝒳)=𝜕𝒳𝒳×₁𝐴−ℰ₁×₁𝐴^𝑇,𝐽₂(𝒳)=𝜕𝒳𝒳×₂𝐵−ℰ₂×₂𝐵^𝑇,𝐽₃(𝒳)=𝜕𝒳𝒳×₃𝐶−ℰ₃×₃𝐶^𝑇.(16) Minimizing 𝐽_𝑖(𝒳)(𝑖=1,2,3) leads to the following recursive equations: 𝒳^𝑘₁=𝒳₁^𝑘−1+𝜇ℰ₁−𝒳₁^𝑘−1×₁𝐴×₁𝐴^𝑇,𝒳^𝑘₂=𝒳₂^𝑘−1+𝜇ℰ₂−𝒳₂^𝑘−1×₂𝐵×₂𝐵^𝑇,𝒳^𝑘₃=𝒳₃^𝑘−1+𝜇ℰ₃−𝒳₃^𝑘−1×₃𝐶×₃𝐶^𝑇,(17) where 𝜇 is iterative step length and is given in the following section.

Substituting (13) into (17) and replacing the unknown variable 𝒳 in (13) with their estimate at time (𝑘−1), we have 𝒳^𝑘₁=𝒳₁^𝑘−1+𝜇ℛ₁^𝑘−1×₁𝐴^𝑇,𝒳^𝑘₂=𝒳₂^𝑘−1+𝜇ℛ₂^𝑘−1×₂𝐵^𝑇,𝒳^𝑘₃=𝒳₃^𝑘−1+𝜇ℛ₃^𝑘−1×₃𝐶^𝑇,(18) where ℛ^𝑘_𝑖=𝒟−𝒳^𝑘_𝑖×₁𝐴−𝒳^𝑘_𝑖×₂𝐵−𝒳^𝑘_𝑖×₃𝐶  (𝑖=1,2,3). Taking the average of 𝒳^𝑘₁, 𝒳^𝑘₂, and 𝒳^𝑘₃, we obtain Algorithm 1.

Input:𝐴,𝐵,𝐶,𝒟, and 𝜇. Output:𝒳. Initialize 𝒳0. For 𝑘=1,2,…,𝑘max  Compute ℛ𝑘−1=𝒟−𝒳𝑘−1×1𝐴−𝒳𝑘−1×2𝐵−𝒳𝑘−1×3𝐶,  𝒳𝑘1=𝒳𝑘−1+𝜇ℛ𝑘−1×1𝐴𝑇,  𝒳𝑘2=𝒳𝑘−1+𝜇ℛ𝑘−1×2𝐵𝑇,  𝒳𝑘3=𝒳𝑘−1+𝜇ℛ𝑘−1×3𝐶𝑇,  𝒳𝑘=(𝒳𝑘1+𝒳𝑘2+𝒳𝑘3)/3. End

Algorithm 1 

The gradient based iterative algorithm (GI).

Remark 2. Note that when we compute 𝒳^𝑘₂ in Algorithm 1, the 𝒳^𝑘₁ has been computed. Similarly, when we compute 𝒳^𝑘₃, the 𝒳^𝑘₁ and 𝒳^𝑘₂ have been computed. Hence, we can fully take advantage of the information of 𝒳^𝑘_𝑖(𝑖=1,2,3) to update the 𝒳^𝑘−1 and get a more efficient algorithm as in Algorithm 2.

Input:𝐴,𝐵,𝐶,𝒟, and 𝜅. Output:𝒳. Initialize 𝒳0. For 𝑘=1,2,…,𝑘max  𝒳𝑘1=𝒳𝑘−1+𝜅(𝒟−𝒳𝑘−1×1𝐴−𝒳𝑘−1×2𝐵−𝒳𝑘−1×3𝐶)×1𝐴𝑇,  𝒳𝑘−1=(𝒳𝑘1+𝒳2𝑘−1+𝒳3𝑘−1)/3,  𝒳𝑘2=𝒳𝑘−1+𝜅(𝒟−𝒳𝑘−1×1𝐴−𝒳𝑘−1×2𝐵−𝒳𝑘−1×3𝐶)×2𝐵𝑇,  𝒳𝑘−1=(𝒳𝑘1+𝒳𝑘2+𝒳3𝑘−1)/3,  𝒳𝑘3=𝒳𝑘−1+𝜅(𝒟−𝒳𝑘−1×1𝐴−𝒳𝑘−1×2𝐵−𝒳𝑘−1×3𝐶)×3𝐶𝑇,  𝒳𝑘=(𝒳𝑘1+𝒳𝑘2+𝒳𝑘3)/3. End

Algorithm 2 

The modified gradient based iterative algorithm (MGI).

4. Convergence Analysis

In this section, we will discuss convergence properties of the algorithms in the previous section. The following theorem provides a sufficient condition to guarantee convergence of the Algorithm 1.

Theorem 3. If the Sylvester tensor equation (1) has a unique solution 𝒳^∗ and 0<𝜇<2/(‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂), then the iterative solution 𝒳^𝑘 given by the Algorithm 1 converges to 𝒳^∗ for arbitrary initial value 𝒳⁰.

Proof. Define the error tensors 𝒳^𝑘∶=𝒳^𝑘−𝒳,𝒳^𝑘_𝑖∶=𝒳^𝑘_𝑖(),𝒰−𝒳,𝑖=1,2,3(19)_𝑘∶=𝒳^𝑘−1×₁𝐴,𝒱^𝑘∶=𝒳^𝑘−1×₂𝒲𝐵,_𝑘∶=𝒳^𝑘−1×₃𝐶.(20) Using (1) and Algorithm 1, the following equalities are trivial: 𝒳^𝑘₁=𝒳^𝑘−1−𝜇ℒ^𝑘×₁𝐴^𝑇,𝒳^𝑘₂=𝒳^𝑘−1−𝜇ℒ^𝑘×₂𝐵^𝑇,𝒳^𝑘₃=𝒳^𝑘−1−𝜇ℒ^𝑘×₃𝐶^𝑇,(21) where ℒ^𝑘=𝒰^𝑘+𝒱^𝑘+𝒲^𝑘. Taking norm on both sides of (21), and using (9)–(11), we have 𝒳^𝑘₁²≤𝒳^𝑘−12+𝜇²‖𝐴‖²₂ℒ^𝑘2−2𝜇𝒰^𝑘,ℒ^𝑘,𝒳^𝑘₂²≤𝒳^𝑘−12+𝜇²‖𝐵‖²₂ℒ^𝑘2−2𝜇𝒱^𝑘,ℒ^𝑘,𝒳^𝑘₃²≤𝒳^𝑘−12+𝜇²‖𝐶‖²₂ℒ^𝑘2−2𝜇𝒲^𝑘,ℒ^𝑘.(22) According to (19) and (22), we have 𝒳^𝑘2≤𝒳^𝑘₁²+𝒳^𝑘₂²+𝒳^𝑘₃²3≤𝒳^𝑘−12+𝜇²3ℒ^𝑘2‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂−2𝜇3‖ℒ‖²=𝒳^𝑘−12−𝜇32−𝜇‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂ℒ^𝑘2≤𝒳⁰²−𝜇32−𝜇‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂^𝑘_𝑖=1ℒ^𝑖2,(23) that is, 𝒳^𝑘2+𝜇32−𝜇‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂^𝑘_𝑖=1ℒ^𝑖2≤𝒳⁰².(24) Since 0<𝜇<2/(‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂), we have ^∞_𝑖=1ℒ^𝑖2≤𝒳⁰²<∞.(25) This implies that as 𝑘→∞, ℒ^𝑘→0, that is, 𝒳^𝑘−1×₁𝐴+𝒳^𝑘−1×₂𝐵+𝒳^𝑘−1×₃𝐶⟶0.(26) According to Lemma 1, 𝒳^𝑘→0 as 𝑘→∞. Thus, the theorem holds.

Based on Theorem 3, we have the following theorem.

Theorem 4. If the Sylvester tensor equation (1) has a unique solution 𝒳^∗ and 0<𝜅<min{1/‖𝐴‖²₂,1/‖𝐵‖²₂,1/‖𝐶‖²₂}, then the iterative solution 𝒳^𝑘 given by the Algorithm 2 converges to 𝒳^∗ for any initial value 𝒳⁰.

Proof. The proof is similar to Theorem 3 and hence omitted.

Now we consider the convergence rate of the Algorithm 1. From (19) and (21), it is not difficult to get the following residual error equation: vec𝒳^𝑘=𝐼_𝑚𝑛𝑙−𝜇3𝑃^𝑘vec𝒳⁰,(27) where 𝑃=𝐼_𝑙⊗𝐼_𝑛⊗𝐴^𝑇𝐴+𝐼_𝑙⊗𝐵⊗𝐴^𝑇+𝐶⊗𝐼_𝑛⊗𝐴^𝑇+𝐼_𝑙⊗𝐵^𝑇⊗𝐴+𝐼_𝑙⊗𝐵^𝑇𝐵⊗𝐼_𝑚+𝐶^𝑇⊗𝐼_𝑛⊗𝐴+𝐶^𝑇𝐼⊗𝐵⊗_𝑚+𝐶⊗𝐵^𝑇⊗𝐼_𝑚+𝐶^𝑇𝐶⊗𝐼_𝑛⊗𝐼_𝑚. We observe that the closer the eigenvalues of 𝜇𝑃 are to 3, the closer the eigenvalues of 𝐼_𝑚𝑛𝑙−(𝜇/3)𝑃 tend to zero, and hence, the faster the residual error 𝒳^𝑘 converges to zero. So the convergence rate of Algorithm 1 depends on the condition number of the associated system [22]; that is, Algorithm 1 has a fast convergence rate for a small condition number of 𝑃. The analysis of Algorithm 2 is similar; hence, it is omitted here.

5. Numerical Experiments

This section gives some examples to illustrate the performance of the proposed algorithms. All the implementations are based on the codes from the MATLAB Tensor Toolbox developed by Bader and Kolda [31, 32].

Example 5. Let ,,𝒟()=()=.𝐴=31−12,𝐵=11−11𝐶=101−2∶,∶,110131511,𝒟∶,∶,214330(28) Then, the solution of 𝒳 is 𝒳(∶,∶,1)=()=1234,𝒳∶,∶,24231.(29) Taking 𝒳⁰=𝟎, 𝜇=2.8∗1/(‖𝐴‖²₂+‖𝐵‖²₂+‖𝐶‖²₂)≈0.1067, 𝜅=2.8∗min{1/‖𝐴‖²₂,1/‖𝐵‖²₂,1/‖𝐶‖²₂}≈0.2747, we apply Algorithms 1 and 2 to compute approximate solution in (1). Since the exact solution is known in this example, we plot the values ||𝒳^𝑘−𝒳||/||𝒳|| to indicate the convergence of the sequence {𝒳^𝑘}. These values are plotted on logarithmic scales, if the curve is approximate decreasing line, then it is indicated that the convergence of the sequence is at least linear.

From Figure 1, we can see that two algorithms are convergent linearly and MGI algorithm converges much faster than GI algorithm.

Figure 1 

Comparison of convergence curves using GI and MGI algorithms.

Remark 6. The effect of changing the step length 𝜇 and 𝜅 is illustrated in Figure 2. We observe that the larger the step length is, the faster the convergence of the algorithm. However, if the step length is too large, the algorithm may diverge. How to choose a best step length is still a issue to be studied [22].

(a)

(b)

Figure 2 

The effect of changing the step length.

Example 7. This example comes from [22], except for 𝐶 and 𝒟 which will be set to a random matrix (and tensor). The MATLAB codes which generate these matrices and compute the step lengths are as follow:  m=30; n=30; l=30; rand(’state’,0); alpha = 3;  A = triu(rand(m,m),1) + diag(alpha + diag(rand(m)));  B = triu(rand(n,n),1) + diag(alpha + diag (rand(n)));  C = triu(rand(l,l),1) + diag(alpha + diag (rand(l)));  D = tenrand([m,n,l]);  mu = 2.2/(norm(A)∧2 + norm(B)∧2 + norm (C)∧2); %Compute 𝜇  kappa = 1.2 * min([1/norm(A)^∧2,1/norm (B)^∧2,1/norm(C)∧2]); %Compute 𝜅.Since the exact solution in this example is unknown, we plot the relative residual errors 𝒳^𝑘×₁𝐴+𝒳^𝑘×₂𝐵+𝒳^𝑘×₃𝐶−𝒟‖𝒟‖.(30)

From Figure 3, we can see that the convergence of GI algorithm tends to slow down after some iterative steps. The MGI algorithm converges linearly, and much faster than GI method.

Figure 3 

Comparison of convergence curves using GI and MGI algorithms.

Example 8. In this example, we consider the effect of the condition number of matrix 𝑃 on two algorithms. The stoping criterion is set to be 10⁻⁶. The MATLAB codes which generate these matrices are as follows:  m=20; n=10; l=5; rand(’state’,0);  A = rand(m,m) + diag(ones(m,1)) * alpha;  B = rand(n,n) + diag(ones(n,1)) * alpha;  C = rand(l,l) + diag(ones(l,1)) * alpha;  D = tenrand([m,n,l]).

For different 𝛼, the number of iterations and the corresponding condition number (𝜅(𝑃)) of the matrix 𝑃 are outlined in Table 1.

Table 1 

The effect of the condition number of matrix 𝑃 on two algorithms.

The results of Table 1 affirm the analysis at the end of Section 4. In particular, we see that the convergence rate becomes faster as the condition number of 𝑃 decreases.

6. Conclusions

We have presented the gradient based iterative algorithms for solving Sylvester tensor equation, Sylvester matrix equation as a special case included. They are based on the hierarchical identification principle and tensor arithmetic concepts. The convergence analysis indicates that the iterative solutions given by the proposed algorithms converge to exact solution for any initial values. Such statement is also confirmed by the given numerical examples. In addition, the idea of this paper can be easily extended to study iterative solutions of high order (𝒳 is 𝑁-mode tensor) Sylvester tensor equation.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 11201092, no. 11261012) and the Scientific Research Starting Foundation for Doctors of Guizhou Normal University, China.