Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

Applied Mathematics and Computation xxx (xxxx) xxx

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equationsR Baohua Huang a,b, Changfeng Ma a,∗ a b

College of Mathematics and Informatics & FJKLMAA, Fujian Normal University, Fuzhou 350117, PR China School of Mathematical Sciences, South China Normal University, Guangzhou 510631, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 29 August 2018 Revised 24 June 2019 Accepted 27 October 2019 Available online xxx

This paper is concerned with some of well-known iterative methods in their tensor forms to solve a class of tensor equations via the Einstein product and the associated with least squares problem. Especially, the tensor forms of the L SQR and L SMR methods are presented. The proposed methods use tensor computations with no matricizations involved. We prove that the norm of residual is monotonically decreasing for the tensor form of the LSQR method. The norm of residual of normal equation is also monotonically decreasing for the tensor form of the LSMR method. We also show that the minimum-norm solution (or the minimum-norm least squares solution) of the tensor equation can be obtained by the proposed methods. Numerical examples are provided to illustrate the efficiency of the proposed methods and testify the conclusions suggested in this paper.

Keywords: Sylvester tensor equations Einstein product LSQR method LSMR method Minimum-norm solution

© 2019 Elsevier Inc. All rights reserved.

1. Introduction For a positive integer N, an order N tensor A = (ai1 ···iN )(1 ≤ i j ≤ I j , j = 1, . . . , N ) is a multidimensional array with I1 I2 . . . IN entries. Let RI1 ×···×IN and CI1 ×···×IN be the sets of the order N dimension I1 ×  × IN tensors over real field R and complex field C, respectively. It is easy to see that an order 2 tensor is a matrix. There has been a lot of research works on tensors, which can be seen in [6,21,22]. Let A ∈ CI1 ×···×IN ×K1 ×···×KN , B ∈ CK1 ×···×KN ×J1 ×···×JM ; the Einstein product of tensors A and B is defined by the operation ∗ N via

(A ∗N B )i1 ···iN j1 ··· jM =



ai1 ···iN k1 ···kN bk1 ···kN j1 ··· jM .

(1.1)

k1 ···kN

It is obvious that A ∗N B ∈ CI1 ×···×IN ×J1 ×···×JM and the associative law of this tensor product is true. By (1.1), if B ∈ CK1 ×···×KN , then

(A ∗N B )i1 ···iN =



ai1 ···iN k1 ···kN bk1 ···kN ,

k1 ···kN

where A ∗N B ∈ CI1 ×···×IN . R Supported by National Key Research and Development Program of China (2018YFC1504200 and 2018YFC0603500) and National Natural Science Foundation of China (41725017). ∗ Corresponding author. E-mail address: [email protected] (C. Ma).

https://doi.org/10.1016/j.amc.2019.124892 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 2

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

The high order Sylvester tensor equation via the Tucker product of tensors is as follows

X ×1 A(1) + X ×2 A(2) + · · · + X ×N A(N ) = D,

(1.2)

where the matrices A(i ) ∈ RIi ×Ii (i = 1, 2, . . . , N ) and the right hand side tensor D ∈ RI1 ×···×IN are known, and X ∈ RI1 ×···×IN is the unknown tensor to be determined. The operator X ×n A denotes the n-mode product of a tensor X ∈ RI1 ×···×IN with a matrix A ∈ RJ×In and the size is I1 × · · · × In−1 × J × In+1 × · · · × IN and its elements are defined by

(X ×n A )i1 ···in−1 jin+1 ···iN =

In 

xi1 i2 ···iN a jin .

in =1

In the special case that X ∈ RI1 ×I2 , (1.2) becomes the Sylvester matrix equation

AX + X BT = D,

(1.3)

which plays vital roles in many applications such as control theory and system theory [9,10,13]. In the special case that X ∈ RI1 ×I2 ×I3 , (1.2) reduces to the following form

X ×1 A + X ×2 B + X ×3 C = D,

(1.4)

which arises from the finite difference, finite element or spectral methods [14,17,18] and plays very important roles in discretization of the linear partial differential equation in high dimension. For example, in the discretization of the three dimensional radiative transform equation using Chebyshev collocation spectral method, the resultant algebraic system has the form (1.2) (see [18]). Consider the partial differential equation [2,27]



−u + cT ∇ u = f,

in  = [0, 1]N ,

u = 0,

on

∂ .

By using the finite-difference discretization, which together with a second order convergent scheme for the convection term, leads to a linear system of the form of Eq. (1.2). Chen and Lu [7] established the projection method to solve this tensor equation. They also applied Kronecker product preconditioner to accelerate the convergence of the iterative method. Later, Beik et al. [4] derived the Krylov subspace methods to solve the Sylvester tensor Eq. (1.4). Shi et al. [23] investigated the backward error and perturbation bounds for the tensor Eq. (1.4). The high order Sylvester tensor equation via the Einstein product is defined in [24], which has the following form

A ∗N X + X ∗M D = C, RI1 ×···×IN ×I1 ×···×IN ,

(1.5) RI1 ×···×IN ×J1 ×···×JM ,

RJ1 ×···×JM ×J1 ×···×JM

RI1 ×···×IN ×J1 ×···×JM ,

where A ∈ X ∈ D∈ and C ∈ comes from the discretization of the linear partial differential equation in high dimension. The special form of the higher order Sylvester tensor Eq. (1.5) is

A ∗N X = C,

(1.6)

where A ∈ RI1 ×···×IN ×J1 ×···×JN , X ∈ RJ1 ×···×JN ×K1 ×···×KM and C ∈ RI1 ×···×IN ×K1 ×···×KM . If N = 2 and X is a matrix, the tensor Eq. (1.6) has the following simple form

A ∗2 X = C, where A ∈ RI1 ×I2 ×I1 ×I2 , X ∈ RI1 ×I2 and C ∈ RI1 ×I2 , which arises from many applications, such as engineering, continuum physics, isotropic and anisotropic elastic models [19]. Recently, Sun et al. [24] investigated the generalized inverses of tensors via the Einstein product. By means of the generalized inverses of tensors, they also gave the general solutions of the tensor Eq. (1.6). Behera and Mishra [3] derived further results on generalized inverses of tensors via the Einstein product. Later, Wang and Xu [26] considered the iterative algorithms to solve the tensor Eq. (1.6). In this paper, we are concerned with the following high order generalized Sylvester tensor equation via the Einstein product

A ∗N X ∗M B + C ∗N X ∗M D = F, RI1 ×···×IN ×J1 ×···×JN ,

(1.7) RK1 ×···×KM ×L1 ×···×LM

RI1 ×···×IN ×L1 ×···×LM

where A, C ∈ B, D ∈ and F ∈ are known tensors, and X ∈ RJ1 ×···×JN ×K1 ×···×KM is an unknown tensor to be determined. For a tensor A ∈ RI1 ×···×IN ×J1 ×···×JN , the Frobenius norm  ·  of  1 A is defined as A = ( i1 ...iN j1 ... jN |ai1 ...iN j1 ... jN |2 ) 2 (see [5]). We say the tensor Eq. (1.7) is consistent, if there exists a ten satisfying (1.7). And, we also say that X  is a solution of (1.7). If the tensor Eq. (1.7) has no solution, we say the tensor sor X  which minimizes F − A ∗N X ∗M B − C ∗N X ∗M D is said to be the least Eq. (1.7) is inconsistent, in addition, a tensor X squares solution of the tensor Eq. (1.7). As we all know, the LSQR [19,20] and LSMR [11] methods are all efficient numerical methods for finding the solution x of linear equation Ax = b and the associated with least squares problem minx∈Rn Ax − b, where A ∈ Rm×n and b ∈ Rm . Karimi and Dehghan [15] applied the global least squares based on tensor form method (GLS-BTF) to approximate the solution of the Sylvester tensor Eq. (1.2). Motivated by the authors in Refs. [11,15,19,20], we are dedicated to derive the tensor form of Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

3

the LSQR and LSMR methods to find the solution or the least squares solution of the tensor Eq. (1.7). That is, we consider the following two cases. Problem 1.1. Given tensors A, C ∈ RI1 ×···×IN ×J1 ×···×JN , B, D ∈ RK1 ×···×KM ×L1 ×···×LM and F ∈ RI1 ×···×IN ×L1 ×···×LM . We find a tensor  ∈ RJ1 ×···×JN ×K1 ×···×KM such that X

 ∗M B + C ∗N X  ∗M D = F. A ∗N X Problem 1.2. Given tensors A, C ∈ RI1 ×···×IN ×J1 ×···×JN , B, D ∈ RK1 ×···×KM ×L1 ×···×LM and F ∈ RI1 ×···×IN ×L1 ×···×LM . We find a tensor  ∈ RJ1 ×···×JN ×K1 ×···×KM such that X

F − A ∗N X ∗M B − C ∗N X ∗M D =

min

X ∈RJ1 ×···×JN ×K1 ×···×KM

F − A ∗N X ∗M B − C ∗N X ∗M D.

The reminder of this paper is organized as follows. Section 2 presents some notations and definitions and lemmas which will be used throughout this paper. Section 3 derives the tensor form of the LSQR method (denoted by LSQR-BTF) for finding the solution of Problems 1.1 and 1.2. The tensor form of the LSMR method (denoted by LSMR-BTF) to solve Problems 1.1 and 1.2 is also given in Section 3. In Section 4, we prove that the minimum-norm solution or the minimum-norm least squares solution of the tensor Eq. (1.7) can be derived by the LSQR-BTF and LSMR-BTF methods. Numerical examples of the proposed methods are shown and analyzed in Section 5. The paper ends up with some conclusions in Section 6. 2. Preliminaries For convenience, we use the following notations throughout this paper. For a tensor A = (ai1 ···iN j1 ··· jM ) ∈ CI1 ×···×IN ×J1 ×···×JM , let B = (bi1 ···iM j1 ··· jN ) ∈ CJ1 ×···×JM ×I1 ×···×IN be the conjugate transpose of A, where bi1 ···iM j1 ··· jN = a j1 ··· jN i1 ···iM . The tensor B is denoted by A∗ . When bi1 ···iM j1 ··· jN = a j1 ··· jN i1 ···iM , the tensor B is called the transpose of A, denoted by AT . We say a tensor D = (di1 ···iN j1 ··· jN ) ∈ CI1 ×···×IN ×I1 ×···×IN is a diagonal tensor if di1 ···iN j1 ··· jN = 0 in the case that the indices i1 iN are different from j1 . . . jN . In this case, if all the diagonal entries di1 ···iN i1 ···iN = 1, then D is called a unit tensor, denoted by I. If all the entries of tensor are zero, we say this tensor is zero tensor, denoted by O. The trace of a tensor A ∈ CI1 ×···×IN ×I1 ×···×IN  is tr(A ) = i1 ···iN ai1 ···iN i1 ···iN [5]. Let I be the product of I1 , I2 , . . . , IN , i.e., I = I1 I2 · · · IN . Similarly, we denote J = J1 J2 · · · JN , K = K1 K2 · · · KM and L = L1 L2 · · · LM . Now, we recall the definition of inner product of two tensors (see [5]). Let A, B ∈ RI1 ×···×IN × J1 ×···×JM ; the inner product of two tensors A, B is defined as

A, B = tr(BT ∗N A ).

(2.1)

Then, for a tensor A ∈ the tensor norm induced by this inner product is Frobenius norm A =  1 ( i1 ···iN j1 ··· jN |ai1 ···iN j1 ··· jN |2 ) 2 . When A, B = 0, we say that A and B are orthogonal. Wang and Xu [26] derived the symmetric property of this inner product. That is, they established the following results. RI1 ×···×IN ×I1 ×···×IN ,

Lemma 2.1 [26]. Let A, B ∈ RI1 ×···×IN ×J1 ×···×JM . Then

A, B = tr(BT ∗N A ) = tr(A ∗M BT ) = tr(B ∗M AT ) = tr(AT ∗N B ) = B, A. In addition, it follows from the definition of inner product and the properties of tensor trace that for any A, B ∈ RI1 ×···×IN ×J1 ×···×JM and any scalar α ∈ R, (1) linearity in the first argument:

α A, B = αA, B, A + C, B = A, B + C, B, (2) positive definiteness: A, A > 0 for all nonzero tensor A. In what follows, we present some basic definitions and useful propositions which will be used later. Proposition 2.1 [24]. Let A ∈ CI1 ×···×IN ×K1 ×···×KN and B ∈ RK1 ×···×KN ×J1 ×···×JM . Then (1) (A ∗N B )∗ = B ∗ ∗N A∗ ; (2) (IN ∗N B ) = B and (B ∗M IM ) = B, where the unit tensor IN ∈ CK1 ×···×KN ×K1 ×···×KN and IM ∈ CJ1 ×···×JM ×J1 ×···×JM . Definition 2.1 [24]. For a tensor A = (ai1 ···iN j1 ··· jM ) ∈ CI1 ×···×IN ×J1 ×···×JM , A(: | j1 · · · jM ) = (a : · · · : j1 · · · jM ) ∈ CI1 ×···×IN is a subblock of A. Vec(A ) is obtained by lining up all the subtensors  in a column, and t-th subblock of Vec(A ) is A(: | j1 · · · jM ) =  M (a : · · · : j1 · · · jM ), where t = jM + M−1 ( j − 1 ) J . p q q= p+1 p=1 Definition 2.2 [24]. The Kronecker product of A ∈ RI1 ×···×IN ×J1 ×···×JN and B ∈ RK1 ×···× KM ×L1 ×···×LM , denoted by A  B = (ai1 ···iN j1 ··· jN B ), is a ‘Kr-block tensor’ whose (s, t) subblock is (ai1 ···iN j1 ··· jN B ) obtained via multiplying all the entries of B     N−1 N−1  N−1 by a constant ai1 ···iN j1 ··· jN , where s = iN + N−1 p=1 (i p − 1 ) q= p+1 Iq and t = jN + p=1 ( j p − 1 ) q= p+1 Jq . Proposition 2.2 [24]. Let A ∈ RI1 ×···×IN ×J1 ×···×JN , B ∈ RK1 ×···×KM ×L1 ×···×LM , C ∈ RK1 ×··· Then (1) (A  B )∗ = A∗  B ∗ ;

×KM ×L1 ×···×LM

and D ∈ RJ1 ×···×JN ×L1 ×···×LM .

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

(2) A  (B  C ) = (A  B )  C; (3) A  (B + C ) = A  B + A  C and (B + C )  A = B  A + C  A; (4) Vec(A ∗N D ∗M B T ) = (B  A ) ∗N Vec(D ). Definition 2.3 [26]. Define the transformation IJ : RI1 ×···×IN ×J1 ×···×JN → RI×J with IJ (A ) = A defined component-wise as

(A )i1 ···iN j1 ··· jN → (A )st where A ∈ RI1 ×···×IN ×J1 ×···×JN , A ∈ RI×J , s = iN +

N−1  p=1

(i p − 1 )

N−1

q= p+1 Iq



and t = jN +

N−1  p=1

( j p − 1)

N−1

q= p+1 Jq



(2.2)

.

The following proposition, which can be seen in [26], serves as a ‘bridge’ between the Einstein multiplication and the usual matrix multiplication. Proposition 2.3 [26]. Let IJ be defined as Definition 2.3. Then

A ∗N X = C ⇔ IJ (A ) · JK (X ) = IK (C ), where · refers to the usual matrix multiplication, A ∈

(2.3) RI1 ×···×IN ×J1 ×···×JN ,

X ∈

RJ1 ×···×JN ×K1 × ···×KM

and C ∈

RI1 ×···×IN ×K1 ×···×KM .

Definition 2.4. Let Qk ∈ RI1 ×···×IN ×J1 ×···×JM (k = 1, 2, . . . ) be a nonzero tensor sequence. We say the tensor sequence {Qk } is an orthogonal tensor sequence if Qi , Q j  = 0 for all i = j and i, j = 1, 2, . . .. 3. The proposed methods By Proposition 2.2, the tensor Eq. (1.7) can be reformulated as

(BT  A + DT  C ) ∗N Vec(X ) = Vec(F ),

(3.1)

where Vec( · ) is defined as in Definition 2.1. By the iterative algorithm proposed in [26], we can solve Problems 1.1 and 1.2. That is, we can find the solution of the tensor Eq. (1.7) when we apply Algorithm 3.1 [26] to the equivalent tensor Eq. (3.1); we can derive the least squares solution of the tensor Eq. (1.7) when we apply Algorithm 4.1 [26] to the equivalent tensor Eq. (3.1). However, the sizes of the coefficient tensors become very large. So, it is necessary to find iterative methods which make full use of the structure of the tensor Eq. (1.7). Now consider the tensor Eq. (1.7). This system in tensor format can be rewritten as



L : RJ1 ×···×JN ×K1 ×···×KM → RI1 ×···×IN ×L1 ×···×LM

(3.2)

L (X ) = F with

L(X ) = A ∗N X ∗M B + C ∗N X ∗M D, RI1 ×···×IN ×J1 ×···×JN

(3.3)

RK1 ×···×KM ×L1 ×···×LM .

where A, C ∈ and B, D ∈ For convenience, we introduce the adjoint operator of L. To this end, we start with the definition of adjoint operator. Definition 3.1. Let L : RJ1 ×···×JN ×K1 ×···×KM → RI1 ×···×IN ×L1 ×···×LM be a linear operator. Then the linear operator

L∗ : RI1 ×···×IN ×L1 ×···×LM → RJ1 ×···×JN ×K1 ×···×KM , which satisfies

L(X ), Y  = X , L∗ (Y ) for all X ∈ RJ1 ×···×JN ×K1 ×···×KM and Y ∈ RI1 ×···×IN ×L1 ×···×LM , is called the adjoint operator of L. Based on Proposition 2.3 and Definitions 2.3 and 3.1, we can establish the following result. Theorem 3.1. Let L be the linear operator defined as in (3.2) and (3.3). Then the adjoint operator L∗ : RI1 ×...×IN ×L1 ×...×LM → RJ1 ×···×JN ×K1 ×···×KM is

L∗ (Y ) = AT ∗N Y ∗M B T + C T ∗N Y ∗M D T . Proof. For arbitrary X ∈ RJ1 ×···×JN ×K1 ×···×KM and Y ∈ RI1 ×···×IN ×L1 ×···×LM , by (3.2) and (3.3), we have

L(X ), Y  = A ∗N X ∗M B + C ∗N X ∗M D, Y  = A ∗N X ∗M B, Y  + C ∗N X ∗M D, Y . RI1 ×···×IN ×J1 ×···×JN ,

RK1 ×···×KM ×L1 ×···×LM ,

RI1 ×···×IN ×L1 ×···×LM

(3.4)

RJ1 ×···×JN ×K1 ×···×KM ,

Since A ∈ B∈ Y∈ and X ∈ by Definition 2.3, we assume that IJ (A ) = A, KL (B ) = B, IL (Y ) = Y and JK (X ) = X. Then, it is obvious that A ∈ RIJ , B ∈ RKL , Y ∈ RIL and X ∈ RJK . According to Lemma 2.1, the definitions of Einstein products ∗ N and ∗ M and Proposition 2.3, we obtain

A ∗N X ∗M B, Y  = tr(Y T ∗N (A ∗N X ∗M B ))  = (Y T ∗N (A ∗N X ∗M B ))l1 l2 ...lM l1 l2 ...lM l1 l2 ...lM

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

=



(Y T (AX B ))ll =

l∈L



5

(Y T AX B )ll

l∈L

= tr(Y T AX B ) = tr(BT X T AT Y )  = tr(X T (AT Y BT ) = (X T (AT Y BT ))kk 

=

k∈K

(X T ∗N (AT ∗N Y ∗M BT ))k1 k2 ...kM k1 k2 ...kM

k1 k2 ...kM

= tr(X T ∗N (AT ∗N Y ∗M B T )) = AT ∗N Y ∗M B T , X  = X , AT ∗N Y ∗M B T .

(3.5)

Similarly, we have

C ∗N X ∗M D, Y  = X , C T ∗N Y ∗M DT , which, together with (3.4) and (3.5), yields

L(X ), Y  = A ∗N X ∗M B, Y  + C ∗N X ∗M D, Y  = X , AT ∗N Y ∗M BT + C T ∗N Y ∗M DT . It then follows from Definition 3.1 that

L∗ (Y ) = AT ∗N Y ∗M B T + C T ∗N Y ∗M D T , which completes the proof.



3.1. Operator bidiagonalization process For linear equation Ax = b and the least squares problem minx∈Rn b − Ax, where A ∈ Rm×n and b ∈ Rm , the least squares QR-factorization (LSQR) method [19,20] and the LSMR method [11] are based on the bidiagonalization procedure of Golub and Kahan [12] listed as Algorithm 3.1. Algorithm 3.1 Bidiagonalization process [12]. 1: 2: 3: 4: 5: 6: 7: 8: 9:

Given matrix A ∈ Rm×n (m ≥ n ), initial vector u1 ∈ Rm satisfying u1 2 = 1, v0 = 0 ∈ Rn , β1 = 0, i := 1. for i = 1, 2, . . ., k do vi = AT ui − βi vi−1 ; αi =  vi ; vi = vi /αi ; ui+1 = Avi − αi ui ; βi+1 =  ui+1 ; ui+1 = ui+1 /βi+1 ; end for It is easy to verify that

vi ∈ span{AT u1 , (AT A )AT u1 , . . . , (AT A )i−1 AT u1 }  Ki (AT A, AT u1 ) and

ui ∈ Ki (AAT , u1 ), vi ∈ Ki (AT A, v1 ). The following properties presented in [25] illustrate that Algorithm 3.1 has finite termination property. Proposition 3.1 [25]. Suppose that k steps of Algorithm 3.1 have been taken. Then the vectors v1 , v2 , . . . , vk and u1 , u2 , . . . , uk+1 are the orthonormal basis of the Krylov space Kk (AT A, v1 ) and Kk+1 (AAT , u1 ). Proposition 3.2 [25]. Algorithm 3.1 will stop at step m if and only if min{l, k} = m, where l is the grade of v1 with respect to AT A and k is the grade of u1 with respect to AAT . In order to derive the tensor form of the LSQR method to solve Problems 1.1 and 1.2, we first give the operator bidiagonalization procedure, which can be seen in Algorithm 3.2. It is easy to verify that

Vi ∈ span{L∗ (U1 ), (L∗ L )L∗ (U1 ), . . . , (L∗ L )i−1 L∗ (U1 )}  Ki (L∗ L, L∗ (U1 )) and

Ui ∈ Ki (LL∗ , U1 ), Vi ∈ Ki (L∗ L, V1 ). Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

Algorithm 3.2 Operator-bidiagonalization process. 1: 2: 3: 4: 5: 6: 7: 8: 9:

Given operator L : RJ1 ×···×JN ×K1 ×···×KM → RI1 ×···×IN ×L1 ×···×LM , initial tensor U1 ∈ RI1 ×···×IN ×L1 ×···×LM satisfying U1  = 1, V0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM , β1 = 0, i := 1. for i = 1, 2, . . ., k do i = L∗ (Ui ) − βi Vi−1 ; V i ; αi =  V i /αi ; Vi = V i+1 = L(Vi ) − αi Ui ; U βi+1 = U i+1 ; i+1 /βi+1 ; Ui+1 = U end for

In order to derive some other important properties of the operator-bidiagonalization process, we first recall the following definitions, which can be seen in [16]. For a given tensor X ∈ RI1 ×···×IN , the notation

X::···:k , k = 1, 2, . . . , IN denotes a tensor in RI1 ×I2 ×···×IN−1 , which is obtained by fixing the last index and is called frontal slice, for more details, one can refer to [16]. Definition 3.2 [16]. The n-mode product of a tensor X = (xi1 i2 ···iN ) ∈ RI1 ×I2 ×···×IN with a vector v = (vin ) ∈ RIn is denoted by ¯ n v, which is an (N − 1 )-mode tensor and X×

(X ׯ n v )i1 ···in−1 in+1 ···iN =

In 

xi1 i2 ···iN vin .

in =1

The following lemmas are also needed in our analysis. Lemma 3.1 [8]. Given tensor X ∈ RI1 ×I2 ×···×IN . (1) Given matrices A ∈ RJm ×Im , B ∈ RJn ×In (m = n), one has

X ×n A ×m B = X ×m B ×n A. (2) Given matrices A ∈ RJ×In , B ∈ RIn ×J , one has

X ×n A ×n B = X ×n (BA ). Lemma 3.2 [4]. If X ∈ RI1 ×I2 ×···×IN , A ∈ RJn ×In , y ∈ RJn , then

¯ ny = X × ¯ n ( AT y ). X ×n A × Lemma 3.3 [4]. If X , Y ∈

¯ Ny = X j, X×

(3.6)

RI1 ×I2 ×···×IN

and y = e j such that ej is the jth column of IN × IN identity matrix, then

j = 1, 2, . . . , IN ,

(3.7)

where X j is the jth frontal slice of X . Let Vk be the (M + N + 1 )-mode tensor with frontal slices V1 , V2 , . . . , Vk and Uk+1 be the (M + N + 1 )-mode tensor with frontal slices U1 , U2 , . . . , Uk , Uk+1 , respectively. That is,

Vk := [V1 , V2 , . . . , Vk ],

(3.8)

Uk+1 := [U1 , U2 , . . . , Uk , Uk+1 ]. Let L(Vk ) be the (M + N + 1 )-mode tensor with frontal slices L(V1 ), L(V2 ), . . . , L(Vk ) and mode tensor with frontal slices L∗ (U1 ), L∗ (U2 ), . . . , L∗ (Uk ), L∗ (Uk+1 ), respectively. That is,

(3.9) L∗ ( U

k+1 )

be the (M + N + 1 )-

L(Vk ) := [L(V1 ), L(V2 ), . . . , L(Vk )],

(3.10)

L∗ (Uk+1 ) := [L∗ (U1 ), L∗ (U2 ), . . . , L∗ (Uk ), L∗ (Uk+1 )].

(3.11)

Denote

⎡

α1 ⎢β2 ⎢ Bk = ⎢ ⎢ ⎣

⎤ α2 ..

.

..

.

βk

⎥ ⎥ ⎥. ⎥ αk ⎦ βk+1

(3.12)

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

7

By the procedure of operator-bidiagonalization, we immediately have the following result. Theorem 3.2. Suppose that k steps of Algorithm 3.2 have been taken. Then the following statements hold

L(Vk ) = Uk+1 ×(M+N+1) BTk ,

(3.13)

L∗ (Uk+1 ) = Vk ×(M+N+1) Bk + αk+1 Z ×(M+N+1) Ek ,

(3.14)

L∗ ( U

where Vk , Uk+1 , L(Vk ), k+1 ) and Bk are defined as in (3.8)–(3.12), respectively; Z is an (M + N + 1 )-mode tensor with (M + N )-mode column tensors O, . . . , O, Vk+1 ; Ek is an k × k matrix of the form Ek = [0, . . . , 0, ek ] with ek being the kth column of the identity matrix of order k. Proof. First, we prove (3.13). From the relation (3.9) and the definition of (M + N + 1 )-mode product of a tensor with a matrix, we have

(Uk+1 ×(M+N+1) BTk )i1 i2 ···iN l1 l2 ···lM i =

k+1 

(Uk+1 )i1 i2 ···iN l1 l2 ···lM j (Bk ) ji

j=1

= (Uk+1 )i1 i2 ···iN l1 l2 ···lM i (Bk )ii + (Uk+1 )i1 i2 ···iN l1 l2 ···lM ,i+1 (Bk )i+1,i = αi (Ui )i1 i2 ···iN l1 l2 ···lM + βi+1 (Ui+1 )i1 i2 ···iN l1 l2 ···lM .

(3.15)

By the relation (3.10) and the definition of frontal slice tensor, one has



L ( Vk )





i1 i2 ···iN l1 l2 ···lM i

= L ( Vi )



i1 i2 ···iN l1 l2 ···lM .

(3.16)

According to lines 6–8 of Algorithm 3.2, we obtain

L(Vi ) = αi Ui + βi+1 Ui+1 , which, together with (3.16), implies that



L ( Vk )



i1 i2 ···iN l1 l2 ···lM i=



L ( Vi )



i1 i2 ···iN l1 l2 ···lM =

αi (Ui )i1 i2 ···iN l1 l2 ···lM+ βi+1 (Ui+1 )i1 i2 ···iN l1 l2 ···lM .

(3.17)

Hence, combining (3.15) and (3.17), we conclude that (3.13) holds. Second, we claim that (3.14) is true. From the relation (3.11), the definition of frontal slice tensor and lines 3–6 of Algorithm 3.2, it follows that



L∗ (Uk+1 )





j1 j2 ··· jN k1 k2 ···kM i

= L∗ ( Ui )



j1 j2 ··· jN k1 k2 ···kM

= βi (Vi−1 ) j1 j2 ··· jN k1 k2 ···kM + αi (Vi ) j1 j2 ··· jN k1 k2 ···kM , ∀i = 1, 2, . . . , k, k + 1.

(3.18)

By (3.8), after some calculations, it is easy to see that



Vk ×(M+N+1) Bk + αk+1 Z ×(M+N+1) Ek



j1 j2 ··· jN k1 k2 ···kM i

= (Vk ×(M+N+1) Bk ) j1 j2 ··· jN k1 k2 ···kM i =

k 

(Vk ) j1 j2 ··· jN k1 k2 ···kM j (Bk )i j

j=1

= (Vk ) j1 j2 ··· jN k1 k2 ···kM i (Bk )ii + (Vk ) j1 j2 ··· jN k1 k2 ···kM ,i−1 (Bk )i,i−1

= αi (Vi ) j1 j2 ··· jN k1 k2 ···kM + βi (Vi−1 ) j1 j2 ··· jN k1 k2 ···kM , ∀i = 1, 2, . . . , k and



Vk ×(M+N+1) Bk + αk+1 Z ×(M+N+1) Ek

 j1 j2 ··· jN k1 k2 ···kM ,k+1



= (Vk ×(M+N+1) Bk ) j1 j2 ···iN k1 k2 ···kM ,k+1 + αk+1 Z ×(M+N+1) Ek =

k 

(3.19)

 j1 j2 ··· jN k1 k2 ···kM ,k+1

(Vk ) j1 j2 ··· jN k1 k2 ···kM j (Bk )k+1, j + αk+1 (Vk+1 ) j1 j2 ··· jN k1 k2 ···kM

j=1

= βk+1 (Vk ) j1 j2 ··· jN k1 k2 ···kM + αk+1 (Vk+1 ) j1 j2 ··· jN k1 k2 ···kM . It then follows from (3.18)–(3.20) that



L∗ (Uk+1 )





j1 j2 ··· jN k1 k2 ···kM i

= Vk ×(M+N+1) Bk + αk+1 Z ×(M+N+1) Ek

(3.20)



for i = 1, 2, . . . , k, k + 1. Hence, we obtain (3.14). The proof is completed.

j1 j2 ··· jN k1 k2 ···kM i



Moreover, similar to Proposition 3.1, we immediately have the following result, which illustrates that Algorithm 3.2 has finite termination property. Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

Proposition 3.3. Suppose that k steps of Algorithm 3.2 have been taken. Then the tensors V1 , V2 , . . . , Vk and U1 , U2 , . . . , Uk+1 are the orthonormal basis of the Krylov space Kk (L∗ L, V1 ) and Kk+1 (LL∗ , U1 ). Proof. This will be shown by induction on k. For k = 1, by lines 4–5 of Algorithm 3.2, we immediately have V1  = 1. By lines 7–8 of Algorithm 3.2, U2  = 1. According to Algorithm 3.2 and Definition 3.1, we have

  L∗ (U1 ), V1  − α1 U1 2 α1 V1 , V1  − α1 U1 2 U1 , U2  = U1 , L(V1 ) − α1 U1 /β2  = = = 0. β2 β2 Assume now that the result is true for some k. Then for k + 1, by Algorithm 3.2, Definition 3.1 and the induction principle, one has





L(Vi ), Uk+1  − βk+1 Vi , Vk  αk+1 αi Ui + βi+1 Ui+1 , Uk+1  − βk+1 Vi , Vk  = αk+1 ⎧  α U + β U i+1 i+1 , Uk+1  − βk+1 Vi , Vk  ⎪ , 1≤i≤k−1 ⎨ i i αk+1 = ⎪ ⎩ αk Uk + βk+1 Uk+1 , Uk+1  − βk+1 Vk , Vk  , i = k αk+1

Vi , Vk+1  = Vi ,

L∗ (Uk+1 ) − βk+1 Vk

αk+1

=

= 0. Similarly, we have

  L(Vk+1 ) − αk+1 Uk+1 L∗ (Ui ), Vk+1  − αk+1 Ui , Uk+1  Ui , Uk+2  = Ui , = βk+2 βk+2 αi Vi + βi Vi−1 , Vk+1  − αk+1 Ui , Uk+1  = βk+2 ⎧ α V + β V , V  − α U , U  i i i i−1 i k+1 k+1 k+1 ⎪ , 1≤i≤k ⎨ βk+2 = ⎪ ⎩ αk+1 Vk+1 + βk+1 Vk , Vk+1  − αk+1 Uk+1 , Uk+1  , i = k + 1 βk+2 = 0.

And, it is easy to see that Vk+1  = 1 and Uk+2  = 1 by lines 4–5 and lines 7–8 of Algorithm 3.2, respectively. The proof is completed.  3.2. LSQR-BTF: LSQR method based on tensor form In this subsection, we describe how the LSQR method can be extended based on tensor format. Let X0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM be the initial tensor and R0 denotes its corresponding residual. Take the initial tensor U1 in operator bidiagonalization process as U1 = R0 /β1 (β1 is redefined and β1 = R0 ). Based on Algorithm 3.2, the LSQR-BTF method constructs the approximate solution Xk of the form

¯ (M+N+1) yk = Vk × ¯ (M+N+1) yk , yk ∈ Rk . Xk = X0 + Vk ×

(3.21)

Then the residual at kth iterative step is defined as

Rk = F − L ( Xk ). According to the linearity of the operator L, it is obvious that

¯ (M+N+1) yk , Rk = R0 − L ( Vk )× where L(Vk ) is defined as in (3.10). By the fact R0 = β1 U1 , Lemma 3.2 and Theorem 3.2, we get

¯ (M+N+1) yk Rk = R0 − L ( Vk )× =

β1 U1 − L(Vk )ׯ (M+N+1) yk

¯ (M+N+1) (β1 e1(k+1) − Uk+1 ×(M+N+1) BTk × ¯ (M+N+1) yk = Uk+1 × ¯ (M+N+1) (β1 e1(k+1) ) − Uk+1 × ¯ (M+N+1) (Bk yk ) = Uk+1 × ¯ (M+N+1) (β1 e1(k+1) − Bk yk ), = Uk+1 ×

(3.22)

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

9

where e1(k+1 ) is the first column of identity matrix of order k + 1. Similar to the LSQR method, for the LSQR-BTF method, we

1) 2) k+1 ) T wish to minimize Rk . Let tk+1 := β1 e1(k+1 ) − Bk yk and tk+1 := (tk(+1 , tk(+1 , . . . , tk(+1 ) . Then, by (3.9) and the definition of (M + N + 1 )-mode product of a tensor with a vector, (3.22) can be rewritten as 1) 2) k+1 ) Rk = tk(+1 U1 + tk(+1 U2 + · · · + tk(+1 Uk+1 .

From the definition of tensor inner product and tensor norm, we have 1) 2) k+1 ) 1) 2) k+1 ) Rk 2 = tk(+1 U1 + tk(+1 U2 + · · · + tk(+1 Uk+1 , tk(+1 U1 + tk(+1 U2 + · · · + tk(+1 Uk+1 .

It then follows from Proposition 3.3 that

  k+1  j )   Rk  =  (tk(+1 )2 = tk+1  = β1 e1(k+1) − Bk yk . j=1

Hence, we have the following subproblem

    Rk  = β1 e1(k+1) − Bk yk  = mink β1 e1(k+1) − Bk y.

(3.23)

y∈R

Efficient solution of the least squares subproblem (3.23) is the heart of the LSQR-BTF method, which is similar to the LSQR method, one can refer to [19,20]. For completeness, we describe the detailed process. As in LSQR method, we form the QR decomposition

⎡



Qk Bk



 Rk β1 e1(k+1) = 0

zk ζ



k+1

⎢ ⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎢ ⎣

ρ1

θ2 ρ2

θ3 ..

..

.

.

ρk−1

θk ρk 0

⎤ ζ1 ζ2 ⎥ ⎥ .. ⎥ ⎥ . ⎥, ⎥ ζk−1 ⎥ ⎥ ζk ⎦ ζk+1



where Qk = Qk,k+1 · · · Q2,3 Q1,2 and Qk,k+1 is the kth plane rotation which operates on rows k and k + 1 of Bk β1 e1(k+1 )

 in

order to eliminate the subdiagonal element βk+1 . This gives the simple recurrence relation



ck

sk

sk

−ck



ρk βk+1

0

αk+1

ζk 0





=

ρk 0

θk+1 ρk+1

 ζk , ζk+1

1 = α1 , ζ1 = β1 , and the scalars ck and sk are the nontrivial elements of Qk,k+1 . The quantities ρ k and ζk are interwhere ρ mediate scalars that are subsequently replaced by ρ k and ζ k . Let Fk := Vk ×(M+N+1 ) R−T and Fk = [F1 , F2 , . . . , Fk ](Fk is the (M + N + 1 )-mode tensor with frontal slices F j for j = k 1, 2, . . . , k). By Lemmas 3.2 and 3.3, we have

¯ (M+N+1) e1(k ) , Vk × ¯ (M+N+1) e2(k ) , . . . , Vk × ¯ (M+N+1) ek(k ) ] Vk = [V1 , V2 , . . . , Vk ] = [Vk × = Vk ×(M+N+1) I (k ) = Vk ×(M+N+1) (RTk R−T ) k = Vk ×(M+N+1) R−T ×(M+N+1) RTk = Fk ×(M+N+1) RTk , k

(3.24)

where I(k) is the k × k identity matrix. Meanwhile, from the definition of (M + N + 1 )-mode product of a tensor with a matrix, for j = 1, . . . , k, we have

(Fk ×(M+N+1) RTk )i1 i2 ···iN j =

k 

(Fk )i1 i2 ···iN j (Rk )i j

i=1

= (Fk )i1 i2 ···iN , j−1 (Rk ) j−1, j + (Fk )i1 i2 ···iN j (Rk ) j j =

θ j (F j−1 )i1 i2 ···iN + ρ j (F j )i1 i2 ···iN ,

(3.25)

where θ1 = 0. Together (3.24) with (3.25) yields

F1 = (1/ρ1 )V1 , F j+1 = (1/ρ j+1 )(V j+1 − θ j+1 F j ), j = 1, 2, . . . , k − 1.

(3.26)

According to Lemma 3.2, (3.21) and (3.26), we have

¯ (M+N+1) yk = Vk × ¯ (M+N+1) (R−1 ¯ (M+N+1) zk Xk = Vk × zk ) = Vk ×(M+N+1) R−T × k k ¯ (M+N+1) zk = Xk−1 + ζk Fk . = Fk ×

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

JID: AMC 10

ARTICLE IN PRESS

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

Algorithm 3.3 LSQR-BTF method for solving (1.7). Given operator RI1 ×···×IN ×L1 ×···×LM . 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21:

L : RJ1 ×···×JN ×K1 ×···×KM → RI1 ×···×IN ×L1 ×···×LM ,

Set initial tensor X0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM ; 1 = L∗ (U1 ), α1 = V 1 , V1 = V 1 /α1 ; β1 = F , U1 = F/β1 , V  1 = α1 ; Set Z1 = V1 , ζ1 = β1 , ρ For k = 1, 2, . . . until convergence, Do: k+1 = L(Vk ) − αk Uk ; U βk+1 = U k+1 ; k+1 /βk+1 Uk+1 = U k+1 = L∗ (Uk+1 ) − βk+1 Vk ; V k+1 ; αk+1 = V k+1 /αk+1 ; Vk+1 = V

ρk =

ρk2 + βk2+1 ; k /ρk ; ck = ρ sk = βk+1 /ρk ; θk+1 = sk αk+1 ; ρk+1 = −ck αk+1 ; ζk = ck ζk ; ζk+1 = sk ζk ; Xk+1 = Xk + (ζk /ρk )Zk ; Zk+1 = Vk+1 − (θk+1 /ρk )Zk ; If |ζk+1 | is small enough then stop

the

right

hand

side

tensor

F∈

// Initialize operator bidiagonalization

//Continue operator bidiagonalization

//Construct and apply the next orthogonal transformation

//Update X and Z //Test convergence

End Do

The update of Fk+1 appears to require knowledge of ρk+1 , which is not available during the kth iteration. Following the idea of [19], we can define Zk := ρk Fk in order to circumvent this drawback. We then initialize Z1 := V1 and update according to Zk+1 = Vk+1 − (θk+1 /ρk )Zk . And, we can update Xk by the formula Xk = Xk−1 + ζk /ρk Zk . The main steps of the LSQR-BTF method are summarized as Algorithm 3.3. Remark 3.1. The stopping criterion can be chosen as Rk  = |ζk+1 |, where Rk is the kth residual. Other stopping criteria can also be used and are not listed here. Readers can see [19,20] for details. Remark 3.2. As ζk+1 = sk ζk , we conclude that

Rk  ≤ Rk−1 , where Rk = D − L(Xk ) is the kth residual. That is, Rk  is monotonically decreasing. Moreover, if there exists some constant η < 1 such that |sk | ≤ η < 1 for k = 1, 2, . . . , then limk→∞ Rk  = 0. ¯ (M+N+1 ) yk ∈ Kk (L∗ L, V1 ) generated by the LSQR–BTF method converges to the unique Remark 3.3. The sequence Xk = Vk × minimum-norm solution or the unique minimum-norm least squares solution of the tensor Eq. (1.7). (The detailed proof is presented in Section 4.) 3.3. LSMR-BTF: LSMR method based on tensor form In this section, we are devoted to deriving the tensor form of the LSMR method to solve the tensor Eq. (1.7). Let X0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM be the initial tensor and R0 denotes its corresponding residual. Take the initial tensor U1 in operator bidiagonalization process as U1 = R0 /β1 (β1 is redefined and β1 = R0 ). Similar to the LSMR method [11], we utilize the operator bidiagonalization process to solve the normal equation of (1.7), i.e., L∗ (L(X )) = L∗ (F ). Assume that the LSMR-BTF method forms approximate solution

¯ (M+N+1) yk = Vk × ¯ (M+N+1) yk , yk ∈ Rk . Xk = X0 + Vk ×

(3.27)

It then follows from the linearity of the operator L that the residual at kth iterative step is

¯ (M+N+1) yk , Rk = D − L ( Xk ) = R0 − L ( Vk )×

(3.28)

where L(Vk ) is defined as in (3.13). By (3.22), one has

¯ (M+N+1) (β1 e1(k+1) − Bk yk ), Rk = Uk+1 × Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

11

where e1(k+1 ) is the first column of identity matrix of order k + 1. Let β¯ k = αk βk for all k. Then, by Theorem 3.2 and Lemma 3.2, we have



¯ (M+N+1) (β1 e1(k+1) − Bk yk ) L∗ (Rk ) = L∗ Uk+1 ×



¯ (M+N+1) (β1 e1(k+1) − Bk yk ) = L∗ (Uk+1 )× ¯ (M+N+1) (β1 e1(k+1) − Bk yk ) = (Vk ×(M+N+1) Bk + αk+1 Z ×(M+N+1) Ek )× ¯ (M+N+1) (β1 e1(k+1) − Bk yk ) + αk+1 Z ×(M+N+1) Ek × ¯ (M+N+1) (β1 e1(k+1) − Bk yk ) = Vk ×(M+N+1) Bk ×







¯ (M+N+1) BTk (β1 e1(k+1) −Bk yk ) + αk+1 Z × ¯ (M+N+1) EkT (β1 e1(k+1) −Bk yk ) = Vk ×



¯ (M+N+1) (α1 β1 e1(k+1) −BTk Bk yk ) − αk+1 βk+1 (eTk yk )Vk+1 = Vk ×

! ! " " BTk Bk y α1 β1 e1(k+1) − αk+1 βk+1 eTk k ! ! T " " B B ¯ (M+N+1) β¯ 1 e1(k+1) − ¯ k k T yk . = Vk+1 × βk+1 ek ¯ (M+N+1) = Vk+1 ×

Denote

pk+1 := β¯ 1 e1(k+1) −

!

(3.29)

"

BTk Bk ) +1 ) T y and pk+1 := ( p(k1+1 , . . . , p(kk+1 ) . ¯ βk+1 eT k

(3.30)

k

Then, we can relax (3.29) as ) +1 ) L∗ (Rk ) = p(k1+1 V1 + · · · + p(kk+1 Vk+1 .

(3.31)

According to Proposition 3.3, we have

L∗ (Rk )2 = L∗ (Rk ), L∗ (Rk ) ) +1 ) ) +1 ) =  p(k1+1 V1 + · · · + p(kk+1 Vk+1 , p(k1+1 V1 + · · · + p(kk+1 Vk+1  k+1 

=

j) 2 ( p(k+1 ) =  pk+1 2

j=1

 ! T " 2  (k+1)  B B ¯  = β1 e1 − ¯ k k T yk  . βk+1 e k

Hence, we obtain the subproblem

  ! T "  ! T "      Bk Bk Bk Bk ¯ 1 e1 − ¯   L ∗ ( R k )  =  β y = min β e − y . 1 1 k T T ¯ ¯   βk+1 ek βk+1 ek  y∈R k

(3.32)

Efficient solution of the least squares subproblem (3.32) is the heart of the LSMR-BTF method, which is similar to the LSMR method, one can refer to [11]. For completeness, we describe the detailed process. As in the LSQR method, we form the first QR decomposition, i.e.,

⎛ρ

⎜ ⎜ , Rk = ⎜ ⎜ 0 ⎝

! "

ρ2

Rk

Qk+1 Bk =

⎞

θ2

1

..

.

..

.

⎟ ⎟ ⎟. ⎟ θk ⎠ ρk

(3.33)

The first QR decomposition is proceeded by constructing the k-th plane rotation Qk,k+1 to operate on rows k and k + 1 of Bk to annihilate βk+1 . That is

!

ck

sk

−sk

ck

"!

α¯ k βk+1

"

0

αk+1

=

! ρk 0

" θk+1 , α¯ k+1

where α¯ 1 = α1 . Let t := Rk y (tk := Rk yk ). Define also qk := R−T (β¯ k+1 ek ) = γk ek , where γk := β¯ k+1 /ρk . We now perform the k second QR decomposition

⎛

! Q¯ k+1

RTk

γk eTk

β¯ 1 e1 0

"

! =

R¯ k 0

"

⎜

⎜ zk , R¯ k = ⎜ ⎜ ζ¯k+1 ⎝

ρ¯ 1

⎞

θ¯2 ρ¯ 2

..

.

..

.

⎟ ⎟ ⎟. ⎟ θ¯k ⎠ ρ¯ k

(3.34)

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

The second QR decomposition is determined by constructing the kth plane rotation Q¯ k,k+1 in order to eliminate the subdiagonal of RTk , i.e., θk+1 . This gives the simple recurrence relation

!

c¯k

s¯k

−s¯k

c¯k

"!

c¯k−1 ρk

θk+1

0

ρk+1

ζ¯k

"

0

!

=

ρ¯ k 0

θ¯k+1 c¯k ρk+1

" ζk , ζ¯k+1

where ζ¯1 = α1 β1 . Moreover, together (3.33) with (3.34), yields

  ! T "  ! T "      Rk Rk B B k k ¯  β¯ 1 e1 − L ( R k )  =  y β1 e1 − β¯ k+1 eT yk  = min  T k ¯ y∈R βk+1 q Rk  k  ! ! T " " ! k"     zk Rk R¯ k  = min  β¯ 1 e1 − t = min  − t . T    k k ¯ t∈R γk ek t∈R 0  ζk+1 ∗

(3.35)

Hence, the subproblem (3.32) has been reformulated as the subproblem (3.35). The solution is obtained by solving R¯ k tk = zk and the value of the residual is |ζ¯k+1 |, i.e., L∗ (Rk ) = |ζ¯k+1 |. As in the LSQR-BTF method, let Fk := Vk ×(M+N+1 ) R−T and Fk = [F1 , F2 , . . . , Fk ] (Fk is the (M + N + 1 )-mode tensor with k frontal slices F j for j = 1, 2, . . . , k). By (3.26), we have

F1 = (1/ρ1 )V1 , F j+1 = (1/ρ j+1 )(V j+1 − θ j+1 F j ), j = 1, 2, . . . , k − 1.

(3.36)

We also assume that Gk = Fk ×(M+N+1 ) R¯ −T and Gk = [G1 , G2 , . . . , Gk ] (Gk is the (M + N + 1 )-mode tensor with frontal slices k G j for j = 1, 2, . . . , k). Similar to (3.24) and (3.25), we obtain

G1 = (1/ρ¯ 1 )F1 , G j+1 = (1/ρ¯ j+1 )(F j+1 − θ j+1 G j ), j = 1, 2, . . . , k − 1.

(3.37)

It then follows from (3.27), (3.36), (3.37) and Lemma 3.2 that

¯ (M+N+1) yk = Vk × ¯ (M+N+1) (R−1 ¯ (M+N+1) tk Xk = Vk × t ) = Vk ×(M+N+1) R−T × k k k ¯ (M+N+1) tk = Fk ×(M+N+1) R¯ −T ¯ (M+N+1) (R¯ k tk ) = Fk × × k ¯ (M+N+1) zk = Xk−1 + ζk Gk . = Gk ×

(3.38)

The update of Fk+1 appears to require knowledge of ρk+1 , which is not available during the kth iteration. And, the update of Gk+1 appears to require ρ¯ k+1 , which is available only during the next iteration. However, as Fong and Saunders [11] point out, it is possible to bypass this requirement by defining

Hk := ρk Fk , H¯ k := ρk ρ¯ k Gk and update

H¯ k = Hk − [(θ¯k ρk )/(ρk−1 ρ¯ k−1 )]H¯ k−1 , Xk = Xk−1 + [ζk /(ρk ρ¯ k ]H¯ k , Hk+1 = Vk+1 − [θk+1 /ρk ]Hk , where the relations (3.36)–(3.38) are used. The main steps of the LSMR-BTF method can be summarized as shown in Algorithm 3.4. Remark 3.4. The stopping criterion can be chosen as L∗ (Rk ) = L∗ (D − L(Xk )) = |ζ¯k+1 |. Other stopping criteria can also be used and are not listed here. Readers can see [11] for details. Remark 3.5. As ζ¯k+1 = −s¯k ζ¯k , we conclude that

L∗ (Rk ) ≤ L∗ (Rk−1 ), where Rk = D − L(Xk ) is the kth residual. That is, L∗ (Rk ) is monotonically decreasing. Moreover, if there exists some constant η < 1 such that |sk | ≤ η < 1 for k = 1, 2, . . . , then limk→∞ L∗ (Rk ) = 0. ¯ (M+N+1 ) yk ∈ Kk (L∗ L, V1 ) generated by the LSMR–BTF method converges to the unique Remark 3.6. The sequence Xk = Vk × minimum-norm solution or the unique minimum-norm least squares solution of the tensor Eq. (1.7). (The detailed proof is presented in Section 4.) 4. The minimum-norm solution In this section, we illustrate that the minimum norm solution or the minimum-norm least squares solution can be derived by our proposed methods. Theorem 4.1. Suppose that the tensorEq. (1.7) is inconsistent, then Algorithm 3.4 (or Algorithm 3.3) derives the unique minimumnorm least squares solution of the tensor Eq. (1.7). ¯ (M+N+1 ) yk ∈ Kk (L∗ L, V1 ). It then follows from Remark 3.5 that the solution generated Proof. By Algorithm 3.4, Xk = Vk × ∗ by Algorithm 3.4 has the form X = L∗ (M ). And, Remark 3.5 also implies that L∗ (R∗ ) = 0, where R∗ = D − L(X ∗ ). Now  is arbitrary solution of Problem 1.2, first we define tensor assume that X

 − X ∗. Y=X

(4.1)

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

13

Algorithm 3.4 LSMR-BTF method for solving (1.7). Given operator RI1 ×···×IN ×L1 ×···×LM . 1: 2: 3: 4: 5: 6: 7: 8: 9: 10:

16: 17:

ρ¯ k =

12: 13: 14: 15:

18: 19: 20: 21: 22: 23: 24: 25: 26:

the

right

hand

side

tensor

F∈

Set initial tensor X0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM ; // Initialize operator bidiagonalization 1 = L∗ (U1 ), α1 = V 1 , V1 = V 1 /α1 ; β1 = F , U1 = F/β1 , V Set H1 = V1 , H¯ 0 = O ∈ RJ1 ×···×JN ×K1 ×···×KM , α¯ 1 = α1 , ζ¯1 = α1 β1 , ρ0 = 1,ρ¯ 0 = 1, c¯0 = 1, s¯0 = 0; For k = 1, 2, . . . until convergence, Do: //Continue operator bidiagonalization k+1 = L(Vk ) − αk Uk ; U βk+1 = U k+1 ; k+1 /βk+1 Uk+1 = U k+1 = L∗ (Uk+1 ) − βk+1 Vk ; V k+1 ; αk+1 = V k+1 /αk+1 ; Vk+1 = V

α¯ k2 + βk2+1 ; ck = α¯ k /ρk ; sk = βk+1 /ρk ; θk+1 = sk αk+1 ; α¯ k+1 = ck αk+1 θ¯k = s¯k−1 ρk ;

11:

L : RJ1 ×···×JN ×K1 ×···×KM → RI1 ×···×IN ×L1 ×···×LM ,

ρk =

//Construct and apply the next orthogonal transformation

(c¯k−1 ρk )2 + θk2+1 ; c¯k = c¯k−1 ρk /ρ¯ k ; s¯k = θk+1 /ρ¯ k ; ζk = c¯k ζ¯k ; ζ¯k+1 = −s¯k ζ¯k ; H¯ k = Hk − [θ¯k ρk /(ρk−1 ρ¯ k−1 )]H¯ k−1 ; Xk = Xk−1 + [ζk /(ρk ρ¯ k )]H¯ k ; Hk+1 = Vk+1 − [θk+1 /ρk ]Hk If |ζ¯k+1 | is small enough then stop

//Update H, H¯ and X

//Test convergence

End Do

 be the residual of the tensor Eq. (1.7) corresponding to tensor X . That is, Let R

 = D − L (X ). R

(4.2)

 and X are solutions of Problem 1.2, together (4.1) with (4.2), yields Since X

2 = D − L(X )2 = D − L(X ∗ + Y )2 R ∗  2 =  R = R∗ − L(Y )2 = R∗ 2 − 2R∗ , L(Y ) + L(Y )2 = R∗ 2 − 2L∗ (R∗ ), Y  + L(Y )2 = R ∗  2 + L ( Y )  2 . Hence, we conclude that L(Y ) = O. By some calculations, we immediately have

X2 = X ∗ + Y 2 = X ∗ 2 + 2X ∗ , Y  + Y 2 = X ∗ 2 + 2L∗ (M ), Y  + Y 2 = X ∗ 2 + 2M, L(Y ) + Y 2 = X ∗  2 + Y  2 ≥ X ∗  2 . This implies that the solution X ∗ is the minimum-norm solution of Problem 1.2. In other words, X ∗ is the unique minimumnorm least squares solution of the tensor Eq. (1.7). The proof is completed.  Theorem 4.2. Suppose that the tensor Eq. (1.7) is consistent, then Algorithms 3.3 and 3.4 derive the unique minimum-norm solution of the tensor(1.7). Proof. By using the same proof method as that of Theorem 4.1, one can prove this result.



5. Numerical experiments In this section, we report some numerical examples to support our Algorithms 3.3 and 3.4. All of the tests were run on the Intel (R) Core (TM), where the CPU is 2.40 GHz and the memory is 8.0 GB. The programming language was Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC 14

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx Table 1 Numerical results for Example 5.1. Method

Iter

Time

Res

LSQR-BTF LSMR-BTF

90 90

0.6328 0.8848

4.9503e−12 6.6798e−12

MATLAB R2015a. We comment here that the tensor toolbox [1] are utilized for solving the succeeding discussed problems. We also illustrate how to use the Einstein product in Matlab tensor toolbox. For tensors A ∈ RI1 ×···×IN ×K1 ×···×KN and B ∈ RK1 ×···×KN ×J1 ×···×JM , the Einstein product A ∗N B can use the tensor command ‘ttt’ to implement. That is,

A = tenrand([I1 I2 · · · IN K1 K2 · · · KN ] );

//Input the tensor A

B = tenrand([K1 K2 · · · KN J1 J2 · · · JM ] );

//Input the tensor B

A ∗N B := ttt(A, B, [N + 1 N + 2 · · · 2N], [1 2 · · · N] ). //Compute the Einstein product Example 5.1. Consider tensor equation

A ∗N X ∗M B + C ∗N X ∗M D = F RI1 ×I2 ×J1 ×J2 ,

RK1 ×K2 ×L1 ×L2 ,

(5.1) RI1 ×I2 ×L1 ×L2

RJ1 ×J2 ×K1 ×K2 .

with A, C ∈ B, D ∈ F∈ and X ∈ We use the command ‘tenrand’ to generate tensors A, B, C, D. Especially, A = tenrand([6 5 3 3] ), B = tenrand([3 2 4 3] ), C = tenrand([6 5 3 3] ) and D = tenrand([3 2 4 3] ). The exact solution X ∗ of the tensor Eq. (5.1) is a tensor with X ∗ = tenrand([3 3 3 2] ). Then by F = A ∗N X ∗M B + C ∗N X ∗M D, we generate the right hand side tensor F. Take the initial tensor X0 = O. The stopping criteria of the LSQR-BTF method is |ζk+1 | = Rk  < 10−11 , and the LSMRBTF method is Rk  < 10−11 . By Algorithms 3.3 (LSQR-BTF method) and 3.4 (LSMR-BTF method), after 90 iterative steps, we obtain the minimum-norm solution of the tensor Eq. (5.1) as follows.

⎛

⎛

0.8893 0.8909

0.7586⎠, X90 (:, :, 2, 1 ) = ⎝ 0.1672

0.7367

0.6282

0.6401

0.2529

0.0086

0.3923

X90 (:, :, 1, 1 ) = ⎝0.4584

⎛

X90 (:, :, 3, 1 ) = ⎝ 0.1635

0.6641

⎞

0.0081

⎞

0.4831

0.0529⎠,

0.8219

0.0928

0.0108

0.9210

0.7650

⎛

0.4337

0.8185

0.9549

0.6424

0.2388

0.6925

0.9142

0.0640⎠, X90 (:, :, 3, 2 ) = ⎝0.9646

0.3915

0.4356

0.8611

X90 (:, :, 2, 2 ) = ⎝0.7227

⎞

0.1534

0.1081 ⎠, X90 (:, :, 1, 2 ) = ⎝0.4894

⎞

0.4745

0.5580

0.1008

⎛

0.6523

⎛

⎞

0.0413

0.2764⎠,

0.1741

0.8266

0.8453

0.8937

0.1890

0.3059

0.1780

0.0124 ⎠.

0.7835

0.2925

0.8840

⎞

And, we list the iterative steps (denoted by ‘Iter’), elapsed CPU time in seconds (denoted by ‘Time’) and the norm of residual (denoted by ‘Res’). Here, ‘Res’ is defined as

Res(Xk ) := Rk  = F − L(Xk ), where Xk is the kth approximate solution. The corresponding numerical results are listed in Table 1. Table 1 shows that the iteration number of the LSQR-BTF and LSMR-BTF methods is the same and the CPU time of the LSMR-BTF method is a little more than that of the LSQR-BTF method. This example confirms that both the LSQR-BTF method and the LSMR-BTF method are efficient in solving consistent tensor equation. Example 5.2. Consider tensor equation

A ∗N X ∗M B + C ∗N X ∗M D = F

(5.2)

with A = tenrand([8 7 5 4] ), B = tenrand[4 3 9 5], C = tenrand([8 7 5 4] ), D = tenrand([4 3 9 5] ) and F = 10 ∗ tenrand([8 7 9 5] ). In this case, the tensor Eq. (5.2) is inconsistent. We apply the LSQR-BTF and LSMR-BTF methods to find the least squares solution. Take the initial tensor X0 = O for all these methods. The stopping criteria of the LSMR-BTF method is |ζ¯k+1 | = L∗ (Rk ) < 10−11 , and the LSQR-BTF method is φ¯ k+1 αk+1 |ck | = L∗ (Rk ) < 10−11 as suggested for standard LSQR method (see [20]). We demonstrate the efficiency of the above mentioned iterative methods from the aspects of iterative steps (denoted by ‘Iter’), elapsed CPU time in seconds (denoted by ‘Time’) and the norm of absolute residual (denoted by ‘RES’). Here, ‘RES’ is defined as

RES(Xk ) := L∗ (Rk ) = L∗ (F − L(Xk )), Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

15

Table 2 Numerical results for Example 5.2. Method

Iter

Time

RES

LSQR-BTF LSMR-BTF

232 227

1.4249 1.4818

7.5242e−12 7.6657e−12

Table 3 The iteration number and the CPU time for the Poisson equation in three dimensions. Grid n×n×n

10 × 10 × 10 Iter/time

15 × 15 × 15 Iter/time

20 × 20 × 20 Iter/time

25 × 25 × 25 Iter/time

LSQR-BTF LSMR-BTF BiCG-BTF [5] Algorithm 3.1 [26]

41/0.1877 42/0.2628 43/0.2300 42/0.0958

99/1.7794 99/2.7080 101/0.8938 100/3.0370

166/13.1073 167/19.5284 169/24.1663 168/7.4748

263/74.5025 265/118.1222 266/40.6010 266/124.7202

where Xk is the kth approximate least squares solution of the tensor Eq. (5.2). The corresponding numerical results are listed in Table 2. As seen in Table 2, we find that the CPU time of the LSQR-BTF and LSMR-BTF methods is almost the same. However, the iteration number of the LSMR-BTF method is less than that of LSQR-BTF method. By multiple tests, we also find that the LSQR-BTF and LSMR-BTF methods ultimately converge to similar points, it is safer to use LSMR-BTF method in situations where the solver must be terminated early. This example further confirms that the LSQR-BTF and LSMR-BTF methods are quite efficient in finding the least squares solution of the tensor Eq. (1.7). Example 5.3. Consider the three-dimensional (3D) Poisson problem



−∇ 2 v = f,

in  = {(x, y, z ), 0 < x, y, z < 1},

v = 0,

on

(5.3)

∂ ,

where f is a given function, and

∇ 2v =

∂ 2v ∂ 2v ∂ 2v + + . ∂ x2 ∂ y2 ∂ z 2

We compute an approximation of the unknown function v(x, y, z) in (5.3). Take the uniform mesh step size. That is, the step sizes, x in the x-direction, y in the y-direction and z in the z-direction, satisfy x = y = z = h = 1/(n + 1 ). By the standard central difference approximations, we obtain the difference formula

6vi jk − vi−1 jk − vi+1 jk − vi j−1k − vi j+1k − vi jk−1 − vi jk+1 = h3 fi jk .

(5.4)

Hence, the higher order tensor representation of the 3D discretized Poisson problem (5.3) is (see [5])

A¯ ∗3 V = F,

(5.5) Rn×n×n×n×n×n

Rn×n×n .

where the Laplacian tensor A¯ ∈ and V, F ∈ Both V and F are discretized on the unit cube. The entries ( 2,4,6 ) on the tensor block (A¯ )l,m,p of A¯ would follow a seven-point stencil, i.e.,

((A¯ )α(2,,β4,,γ6) )α,β ,γ =

6 , h3

1 , h3 1 ,6 ) ,6 ) ((A¯ )α(2,,β4−1 ) = ((A¯ )α(2,,β4+1 ) = − 3, ,γ α ,β −1,γ ,γ α ,β +1,γ h 1 ((A¯ )α(2,,β4,,γ6)−1 )α,β ,γ −1 = ((A¯ )α(2,,β4,,γ6)+1 )α,β ,γ +1 = − 3 , h ,4,6 ) ,4,6 ) ((A¯ )α(2−1 ) = ((A¯ )α(2+1 ) =− ,β ,γ α −1,β ,γ ,β ,γ α +1,β ,γ

( 2,4,6 ) where (A¯ )l,m,p = A¯ (:, l, :, m, :, p) ∈ Rn×n×n .

We apply the LSQR-BTF and LSMR-BTF methods to solve the tensor Eq. (1.7). We also compare our proposed methods with the BiCG-BTF method developed by Brazell et al. [5] and Algorithm 3.1 developed by Wang and Xu [26]. The numerical results are reported in Table 3, where ‘Time’ denotes the running time when algorithm terminates, ‘Iter’ denotes the iteration steps and ‘Grid’ denotes the grid size. From Table 3, we find that it needs more calculation time and iterative steps to find the solution of the tensor Eq. (5.5) as the increasing of grid size. As the increasing of grid size, the CPU time and iteration numbers of the LSQR-BTF and LSMR-BTF methods are fewer than those of Algorithm 3.1 in [26]. However, the CPU time of the LSQR-BTF and LSMR-BTF methods is Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892

JID: AMC 16

ARTICLE IN PRESS

[m3Gsc;November 15, 2019;19:12]

B. Huang and C. Ma / Applied Mathematics and Computation xxx (xxxx) xxx

more than that of BiCG-BTF method in [5] with the increasing of grid size. This example also shows that the LSQR-BTF and LSMR-BTF methods are quite effective to solve Poisson problem in three dimensions. From this example, we also find some advantages in computing in a tensor structured domain. First, without vectorization, the solution and the computational grid have a one-to-one correspondence so that sophisticated boundary conditions are easily integrated in the multilinear system. Second, the Laplacian tensor preserves the low bandwidth since the main nodal points sit on the tensor diagonal entries and the rest of the stencil points lie on the off-diagonal positions. Third, the Laplacian tensor A¯ has a low bandwidth, hence the tensor-based methods can reduce the number of operations and storage locations. The same things occur in the biconjugate gradient method (BiCG-BTF method) and Jacobi method of Brazell et al. [5]. 6. Conclusions By means of the LSQR and LSMR methods, we have constructed the tensor form of the LSQR and LSMR methods to solve the tensor Eq. (1.7) and the associated with least squares problem. The proposed methods solely use tensor computations. That is, there is no matrix computations in Algorithms 3.3 and 3.4. The minimum-norm solution (or the minimum-norm least squares solution) of the tensor equation can be derived by Algorithms 3.3 and 3.4. Our text examples show that our proposed methods can solve the consistent or inconsistent tensor Eq. (1.7). Furthermore, we apply our proposed iterative methods to solve the tensor equation in three-dimensional Poisson problem. We can conclude that solving the tensor equation with the tensor-based iterative methods has several advantages: (1) higher order tensor representation of PDEs preserve low bandwidth thereby keeping the computational cost and memory requirement low and (2) a one to one correspondence between the solution and the computational grid facilitate the integration of complicated boundary conditions, which has been illustrated by Brazell et al. [5]. Acknowledgments The authors deeply thank the anonymous referees for helping to improve the original manuscript by valuable suggestions. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

B.W. Bader, T.G. Kolda, et al., MATLAB tensor toolbox version 2.5, 2012, http://www.sandia.gov/tgkolda/TensorToolbox/. J. Ballani, L. Grasedyck, A projection method to solve linear systems in tensor form, Numer. Linear Algebra Appl. 20 (2013) 27–43. R. Behera, D. Mishra, Further results on generalized inverses of tensors via the einstein product, Linear Multilinear Alg. 65 (2017) 1662–1682. A.F.P. Beik, F.S. Movahed, S. Ahmadi-Asl, On the Krylov subspace methods based on tensor format fo positive definite Sylvester tensor equations, Numer. Linear Alg. Appl. 23 (2016) 444–466. M. Brazell, N. Li, C. Navasca, C. Tamon, Solving multilinear systems via tensor inversion, SIAM J. Matrix Anal. Appl. 34 (2013) 542–570. K.C. Chang, K. Pearson, T. Zhang, Perron frobenius theorem for nonnegative tensors, Commun. Math. Sci. 6 (2008) 507–520. Z. Chen, L.Z. Lu, A projection method and Kronecker product preconditioner for solving sylvester tensor equations, Sci. China Math. 55 (2012) 1281–1292. L. De Lathauwer, B. De Moor, J. Vandewalle, A multilinear singular value decomposition, SIAM J. Matrix Anal. Appl. 21 (20 0 0) 1253–1278. F. Ding, T. Chen, Gradient based iterative algorithms for solving a class of matrix equations, IEEE Trans. Autom. Control 50 (2005) 1216–1221. F. Ding, T. Chen, Iterative least squares solutions of coupled sylvester matrix equations, Syst. Control Lett. 54 (2005) 95–107. D.C. Fong, M.A. Saunders, LSMR: an iterative algorithm for sparse least-squares problems, SIAM J. Sci. Comput. 33 (2011) 2950–2971. G.H. Golub, W. Kahan, Calculating the singular values and pseudoinverse of a matrix, SIAM J. Numer. Anal. 2 (1965) 205–224. G.H. Golub, S. Nash, C.F. Van Loan, A Hessenberg–Schur method for the problem AX + XB = c, IEEE Trans. Auto Control 24 (1979) 909–913. L. Grasedyck, Existence and computation of low Kronecker–Rank approximations for large linear systems of tensor product structure, Computing 72 (2004) 247–265. S. Karimi, M. Dehghan, Global least squares method based on tensor form to solve linear systems in Kronecker format, T. I. Meas. Control 40 (2018) 2378–2386. T.G. Kolda, B.W. Bader, Tensor decompositions and applications, SIAM Rev. 3 (2009) 455–500. B.W. Li, Y.S. Sun, D.W. Zhang, Chebyshev collocation spectral methods for coupled radiation and conduction in a concentric spherical participating medium, ASME J. Heat Transf. 131 (2009) 062701–062709. B.W. Li, S. Tian, Y.S. Sun, Z. Mao, Schur-decomposition for 3dmatrix equations and its application in solving radiative discrete ordinates equations discretized by chebyshev collocation spectral method, J. Comput. Phys. 229 (2010) 1198–1212. C.C. Paige, M.A. Saunders, LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw. 8 (1982) 43–71. C.C. Paige, M.A. Saunders, LSQR: Sparse linear equations and least squares problems, ACM Trans. Math. Softw. 8 (1982) 195–209. L. Qi, Eigenvalues of a real supersymmetric tensor, J. Symbolic Comput. 40 (2005) 1302–1324. J. Shao, A general product of tensors with applications, Linear Alg. Appl. 439 (2013) 2350–2366. X. Shi, Y. Wei, S. Ling, Backward error and perturbation bounds for high order Sylvester tensor equation, Linear Multilinear Alg. 61 (2013) 1436–1446. L. Sun, B. Zheng, C. Bu, Y. Wei, Moore–Penrose inverse of tensors via einstein product, Linear Multilinear Alg. 64 (2016) 686–698. F. Toutounian, S. Karimi, Global least squares method (gl-LSQR) for solving general linear systems with several right-hand sides, Appl. Math. Comput. 178 (2006) 452–460. Q.W. Wang, X. Xu, Iterative algorithms for solving some tensor equations, Linear Multilinear Alg. (2018), doi:10.1080/03081087.2018.1452889. H. Xiang, L. Grigori, Kronecker product approximation preconditioners for convection–diffusion model problems, Numer. Linear Alg. Appl. 17 (2010) 723–752.

Please cite this article as: B. Huang and C. Ma, Global least squares methods based on tensor form to solve a class of generalized Sylvester tensor equations, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124892