# Alternating Minimization Methods for Solving Multilinear Systems

Authors: Maolin Liang; Lifang Dai
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Recent works on the multilinear systemAxm−1=b with an order-m and dimension-n tensor A and a vector b of dimension-n have been motivated by their applications in data mining, numerical PDEs, tensor complementary problems, and so on. In this paper, we propose an alternating minimization method for the solution of the system mentioned above and present several randomized versions of this algorithm in order to improve its performance. The provided numerical experiments show that our methods are feasible for any tensor A and outperform some existing ones in the same case.


---

## Body

## 1. Introduction

Many higher-dimensional problems in scientific computing come down to solve the following system of equations:(1)Axm−1=b,where A=ai1i2,…,im is a tensor in ℝm,n, the set of all mth-order n-dimensional tensors over the real field, and the right-hand side b a vector in ℝn, and Axm−1 is a vector in ℝn, whose ith-entry is defined by(2)Axm−1i=∑i2,…,imaii2,…,imxi2xi3,…,xim.We call (1) multilinear system since it is obviously linear over x for each mode. Similar to linear systems case, we also call (1) a nonhomogeneous multilinear system if the right-hand side is nonzero; otherwise, it is a homogeneous one.Many problems in scientific computing including data mining [1], the numerical solution of nonlinear PDEs [2–6], tensor complementary problems [7], and higher-dimensional statistics [8], can be characterized via the multilinear systems mentioned above. To be precise, by using the finite difference method, the discretization of nonlinear PDE,(3)uxm−2⋅∑k=1d∂2∂xk2ux=−fx,inΩ=0,1m,ux=gx,on∂Ω,can be reformulated as the form of (1) [3]; see Section 4 for details. Particularly, if A is a Z-tensor, it has been shown in [7] that the sparsest solution of the tensor complementarity problem,(4)minx0s.t.Axm−1−b≥0,x≥0,xTAxm−1−b=0,can be obtained by solving the following relaxed form:(5)minx1s.t.Axm−1=b,x≥0.In addition, using the “optimize then discretize” approach [8, 9], the differential equation,(6)−maxγ,λ∈Γ,ΛLλU−ηU−12αγ2U+βγ=0,onΩ,U=g,on∂Ω,with the elliptic operator LλUx=1/2σx,λ2U″x+μx,λU′x, Ω=0,1, Γ=0,∞ and compact metric space Λ can be discretized as the well-known 3rd-order Bellman equation, which has the form of (1); see Section 4 for details.Recently, several works concerning the multilinear system (1) have been obtained when the coefficient tensor A is a structured tensor. In [3], Ding and Wei have shown that this system has a unique positive solution when the coefficient tensor A is a nonsingular M-tensor with a positive vector b and then generalized the classical iteration methods for linear systems and the Newton method for nonlinear equations to the multilinear case. Specifically, a homotopy method was also presented by Han [10] for finding the unique positive solution. Very recently, for symmetric M-tensor system, Li et al. [11] and Xie et al. [12], respectively, proposed splitting methods based on the ones in matrix theory and tensor methods for (1). Besides, an algorithm based on the fast Fourier transform was proposed by Xie et al. [13] for the tensor system (1) when the tensor A is a circulant one.In this paper, we continue to study the multilinear system (1) in the case that the contained tensor A is a general one. In contrast to the structured case mentioned above, there is little work on (1) with a general tensor A as far as we are informed. Generally, it is difficult to derive solvability conditions for a general multilinear system as in [3], so we alternatively study the least-squares problem associated with (1), that is,(7)minxAxm−1−b.Throughout this paper, the symbol⋅ denotes the 2-norm of a vector or matrix. We call (7) the tensor least-square (TLS) problem corresponding to the multilinear system (1). The approaches presented in this paper are quite different from the ones mentioned above. By the definition of tensor-vector multiplication, we rewrite the TLS problem (7) as the following optimization problem (see Section 3 for details):(8)minimizefx2,…,xm=12A⋅2x2⋅3x3⋯⋅mxm−b2s.t.x2=⋯=xmxp∈ℝn,p=2,3,…,m.This is a multiblock optimization problem with consensus constrained conditions.Alternating minimization iteration is popular for solving large-scale optimization problems; see, for example, [14, 15], which realizes the reduction of the problem from a nonlinear case to some linear subproblems via updating only one variable while the others are fixed. We apply the alternating minimization strategy to (8), followed by a projection concerning with the consensus constraints. We call it the alternating minimization method for solving the TLS problem, denoted by AMM (i.e., Algorithm 1). Although the alternating minimization scheme and its various variants are widely exploited, their convergence analyses are commonly difficult [16]. We leave the convergence analysis of the AMM an open problem. It should be pointed out that we proposed an alternating iterative method for (1) [17]; this kind of iterative method is ADMM-type one, so the method proposed here is very different from the ADMM-type method.Algorithm 1: AMM: alternating minimization method for the TLS problem (7).
Input: A∈ℝm,n, b∈ℝn.Initialize: xp0∈ℝn with p=2,3,…,m.Roundkk=0,1,2,…,:Step 1. Forp=2,3,…,m, update xpk+1 by (16), and thenStep 2. Computex¯k+1=∑p=2mxpk+1/m−1Notably, a revived interest in the design and analysis of algorithms for solving large-scale optimization problems is to add randomness to the update order of the multiblock updates, which is beneficial to alleviate the iteration complexity and improve the convergence; see, for example, [18–20] and the references therein. As variants of the AMM, we present two random iterations for solving the TLS problem (7); see Algorithms 2 and 3 for details. Besides, Algorithm 4 is another modified version of the AMM. At present, we have not found a similar iterative method for solving multilinear systems.Algorithm 2: Random AMM (R-AMM-I) for problem (8).
Input: A∈ℝm,n, b∈ℝn.Initialize: xp0∈ℝn with p=2,3,…,m.Roundkk=0,1,2,…,:Pick a permutationσ of 2,3,…,m uniformly at random.Step 1. Fori=1,2,…,m−1, update xσik+1 by (16), and thenStep 2. Computex¯k+1=∑p=2mxpk+1/m−1.Algorithm 3: Random AMM (R-AMM-II) for problem (8).
Input:A∈ℝm,n, b∈ℝn.Initialize: xp0∈ℝn with p=2,3,…,m.Roundkk=0,1,2,…,:LetS be a nonempty set with m−1 repeatable elements derived from 2,…,m.Step 1. For everyi∈S, update xik+1 by (16), and thenStep 2. Computex¯k+1=∑p=2mxpk+1/m−1.Algorithm 4: A modified AMM (AMM-III) for the TLS problem (7).
Input: A∈ℝm,n, b∈ℝn.Initialize: z10∈ℝn.Roundkk=0,1,2,…,:Forp=2,…,m, computexpk+1=argminxpfx2k+1,…,xp−1k+1,xp,zp−1k,…,zp−1k,zpk=∑j=2pxjk+1+zp−1k/p.The contribution of this paper is as follows: (1) we considered the least-square problem of the multilinear system (1) by transforming it as an optimization problem. (2) We proposed an alternative iterative method for the optimization problem mentioned above, which is easy to carry out in a parallel version. (3) As modifications of the alternative iterative method, we apply the randomized technique to accelerate it, that is, Algorithms 2 and 3. The performed numerical results show that the randomized methods may have better performance under suitable initial iterative vectors.The outline of this paper is as follows: in Section2, we recall some notations and operations and then give a short review of the classical alternating minimization method. In Section 3, we establish some alternating minimization methods for the TLS problem (7), which are the variants of the classical one. In Section 4, some numerical experiments are performed to illustrate the efficiency and feasibility of the proposed algorithms. Finally, we conclude this paper.


## 2. Preliminaries

### 2.1. Notation and Operations

We will briefly describe notation for tensors and tensor operations used in this paper. Scalars are denoted by lowercase letters, for example,a; vectors are denoted by boldface lowercase letters, for example, a; matrices are denoted by boldface capital letters, for example, I, identity matrix with proper size; tensors are denoted by calligraphic script letters, for example, A. The set of all N-order I1×I2×⋯×IN-dimensional tensors over real field is denoted by ℝI1×I2×⋯×IN; particularly, the set of all mth-order n-dimensional tensors over real field is denoted by ℝm,n.We need the following concepts related to tensors; see, for example, [21] for details.Definition 1.
(tensor-vector multiplication). Let a tensorT∈ℝI1×I2×⋯×IN and a vector v∈ℝIn; then, the n-mode (vector) product of the tensor T with the vector v, denoted by T•nv, is an I1×,⋯,In−1×In+1×⋯×IN tensor, elementwise,(9)T⋅nvi1,…,in−1in+1,…,iN=∑in=1Inti1,…,in,…,iNvin.Particularly, ifT∈ℝI1×I2×I3 and v∈ℝI2, then T⋅2v=∑i2=1I2T:,i2,:vi2, where T:,i2,: stands for the i2th lateral slice of this tensor. In other words, the 2-mode (vector) product between a tensor of order 3 and a vector with the appropriate size is a matrix, the linear combination of the lateral slices taking the associated elements of the vector as the coefficient.By the definition of tensor-by-vector multiplication, one can rewrite (1) as the form of(10)A⋅2x⋅3x⋯⋅mx=b.Definition 2.
(tensor-matrix multiplication). For the tensorT in Definition 1 and the matrix A∈ℝJn×In, the n-mode (matrix) product, denoted by T×nA, is an I1×,⋯,In−1×Jn×In+1×⋯×IN tensor, elementwise,(11)T×nAi1,…,in−1jnin+1,…,iN=∑in=1Inti1,…,in,…,iNajnin.Definition 3.
(matricization/unfolding). For the tensorT in Definition 1, the mode-n matricization (unfolding) of the tensor T, denoted by Tn, is an In×I1I2,…,IN matrix, whose in,j−entry is defined by Tnin,j=ti1i2,…,iN,withj=1+∑k=1,k≠nNik−1∏s=1,s≠nk−1Is.


### 2.2. The Classical Alternating Minimization Method

Consider the following multiblock minimization problem:(12)minimizeFx1,x2,…,xds.t.xp∈Vp⊆ℝnp,p=1,2,…,d,where F:ℝn1×ℝn2×⋯×ℝnd⟶ℝ is a multivariate function. This problem can be solved by the so-called block nonlinear Gauss–Seidel (B-NGS) iteration [15, 22], namely,(13)xpk+1=argminxpFx1k+1,…,xp−1k+1,xp,xp+1k,…,xdk,p=1,2,…,d,which reduces to the ALS method if each subproblem embedded is a least-squares problem; see [16, 21, 23, 24] for some recent developments and applications.


## 2.1. Notation and Operations

We will briefly describe notation for tensors and tensor operations used in this paper. Scalars are denoted by lowercase letters, for example,a; vectors are denoted by boldface lowercase letters, for example, a; matrices are denoted by boldface capital letters, for example, I, identity matrix with proper size; tensors are denoted by calligraphic script letters, for example, A. The set of all N-order I1×I2×⋯×IN-dimensional tensors over real field is denoted by ℝI1×I2×⋯×IN; particularly, the set of all mth-order n-dimensional tensors over real field is denoted by ℝm,n.We need the following concepts related to tensors; see, for example, [21] for details.Definition 1.
(tensor-vector multiplication). Let a tensorT∈ℝI1×I2×⋯×IN and a vector v∈ℝIn; then, the n-mode (vector) product of the tensor T with the vector v, denoted by T•nv, is an I1×,⋯,In−1×In+1×⋯×IN tensor, elementwise,(9)T⋅nvi1,…,in−1in+1,…,iN=∑in=1Inti1,…,in,…,iNvin.Particularly, ifT∈ℝI1×I2×I3 and v∈ℝI2, then T⋅2v=∑i2=1I2T:,i2,:vi2, where T:,i2,: stands for the i2th lateral slice of this tensor. In other words, the 2-mode (vector) product between a tensor of order 3 and a vector with the appropriate size is a matrix, the linear combination of the lateral slices taking the associated elements of the vector as the coefficient.By the definition of tensor-by-vector multiplication, one can rewrite (1) as the form of(10)A⋅2x⋅3x⋯⋅mx=b.Definition 2.
(tensor-matrix multiplication). For the tensorT in Definition 1 and the matrix A∈ℝJn×In, the n-mode (matrix) product, denoted by T×nA, is an I1×,⋯,In−1×Jn×In+1×⋯×IN tensor, elementwise,(11)T×nAi1,…,in−1jnin+1,…,iN=∑in=1Inti1,…,in,…,iNajnin.Definition 3.
(matricization/unfolding). For the tensorT in Definition 1, the mode-n matricization (unfolding) of the tensor T, denoted by Tn, is an In×I1I2,…,IN matrix, whose in,j−entry is defined by Tnin,j=ti1i2,…,iN,withj=1+∑k=1,k≠nNik−1∏s=1,s≠nk−1Is.


## 2.2. The Classical Alternating Minimization Method

Consider the following multiblock minimization problem:(12)minimizeFx1,x2,…,xds.t.xp∈Vp⊆ℝnp,p=1,2,…,d,where F:ℝn1×ℝn2×⋯×ℝnd⟶ℝ is a multivariate function. This problem can be solved by the so-called block nonlinear Gauss–Seidel (B-NGS) iteration [15, 22], namely,(13)xpk+1=argminxpFx1k+1,…,xp−1k+1,xp,xp+1k,…,xdk,p=1,2,…,d,which reduces to the ALS method if each subproblem embedded is a least-squares problem; see [16, 21, 23, 24] for some recent developments and applications.


## 3. Alternating Minimization Methods

This section is devoted to the derivation of the alternating minimization method (AMM) as well as some variants for the solution of the system of multilinear equation (1).To facilitate implementation of the algorithms, we convert equivalently (10) into the following form:(14)A⋅2x2⋅3x3⋯⋅mxm=b,x2=x3=⋯=xm,And, correspondingly, the optimization problem (8) can be written as(15)minimizefx2,…,xm≕12A≠pxp−b2s.t.x2=⋯=xm,where A≠p=A⋅2x2,…,⋅p−1xp−1⋅p+1xp+1,…,⋅mxm∈ℝn×n. For any p∈2,3,…,m, updating one block xp while fixing the remaining ones, we have(16)xpk+1=argminxpfk+1xp≔12A≠pkxp−b2,where(17)fk+1xp≔fx2k+1,…,xp−1k+1,xp,x¯k,…,x¯k,and A≠pk=A⋅2x2k+1,…,⋅p−1xp−1k+1⋅p+1x¯k,…,⋅mx¯k with x¯k=∑p=2mxpk/m−1.The above procedure reformulates the alternating minimization method for solving the TLS problem (7), which can be precisely represented in the following algorithm:We have the following comments in order.Remark 1.
(a)
This algorithm enables that the large-scale TLS problem (7) is solved by several small-size subproblems, which may be computed in an efficient and perhaps highly parallel manner. In other words, Algorithm 1 is more convenient to be parallelized than the existing methods proposed in, for instance, [3, 11, 13].(b)
Obviously, every subproblem contained in Algorithm1 is a linear least-squares problem, and the least-norm solution to the subproblem in (16) has closed form xpk+1=A≠pk†b, where the subscript † stands for the Moore–Penrose inverse of a matrix [25]. Additionally, for the large-scale linear (least-squares) subproblems contained in this algorithm, one can apply the iterative methods (e.g., LSQR, GMRES [25]) to them.(c)
The update in step 2 corresponds to the consensus restrictions in (15). In other words, x¯k+1 is the minimizer of the function gy≔∑p=2my−xpk+12.(d)
From Algorithm1, one can observe that the main amount of calculations hinges on forming the matrices A≠pk in (17) for p=2,3,…,m. So, it is easy to know that the computation complexity of this algorithm is of order Onm−1. Notably, if the tensor A in (1) has the CP decomposition with CP-rank r, the computation complexity of the above algorithm decreases to Om−12nr.Just as mentioned in the Introduction, a revived interest in the design and analysis of algorithms for solving large-scale optimization problems is to add randomness to the updating order of the multiblock variants. The above randomly updating technique can be utilized in Algorithm1. In fact, updating xpk+1 for stochastically chosen p in 2,3,…,m, we obtain the first randomized AMM.In Algorithm2, we randomly update xp by (16) for each p∈2,3,…,m. In this procedure, all the block coordinates are updated after one loop. Furthermore, if we randomly update xp in a subset of 2,3,…,m, we obtain the second randomized AMM.It follows from Algorithm3 that some variants in this algorithm may not be updated within one loop, which is the difference between the above two algorithms. In numerical experiments, we will see that the randomized algorithms presented above indeed have better performance than the original AMM, that is, Algorithm 1.In addition, we note that, in Algorithm1, x¯k+1 is the average value of all the coordinates xpk+1 with p=2,3,…,m. That is to say, the average value is obtained only after all the coordinates are updated. On the other hand, if we can take the average value at once after one coordinate is updated until all the coordinates are contained, we then obtain another modification of Algorithm 1.Certainly, we can also add the randomness in Algorithm4 as in Algorithm 2.We should mention that the alternating minimization methods presented above are numerically feasible and efficient for solving the multilinear system (1) or the associated least-squares problem, which are very different from some existing ones; see, for example, [3, 11, 12], but we have not found a feasible approach to complete their convergence analyses, which will be further studied in future work.As is well-known, the Newton method is an efficient approach in solving nonlinear equations [22]. In [3], the authors established the Newton iteration for solving the multilinear system (1) in the case that the coefficient tensor A is a symmetrical nonsingular M-tensor. When the tensor A in (1) is not symmetric, we can also derive a similar Newton-type iteration method (say NS-Newton for convenience) as presented in [3].Suppose that the tensor system (1) is solvability, we rewrite it as(18)gx≔Axm−1−b=0.This is a nonlinear equation. Notice that the first-order gradient ofgx is ∇gx=∑p=2mA≠p∈ℝn×n, where A≠p=A⋅2x,…,⋅p−1x⋅p+1x,…,⋅mx. Then, the iterations are as follows.If the tensorA is symmetric over the last m−1 modes, we obtain ∇gx=m−1A≠p for any p∈2,3,…,m. In particular, the above algorithm reduces to the Newton iteration proposed in [3] when the tensor A was a symmetrically nonsingular M-tensor. In this case, the matrix ∇gxk is nonsingular as well. It should be emphasized that the matrix ∇gxk may not be nonsingular for a general tensor; in this case, one can replace the inverse with its Moore–Penrose inverse.Our numerical tests show that this is an efficient algorithm. In contrast to Algorithm5. The first three algorithms proposed previously can be implemented in parallel when updating xpk+1, which indicates that these methods are promising for solving large-scale multilinear problems.Algorithm 5: Newton-type method (NS-Newton) for the multilinear system (1).
Input: A∈ℝm,n, b∈ℝn.Initialize: x0∈ℝn.Roundkk=0,1,2,…,:Step 1. Compute∇gxkStep 2. Computerk=b−Axkm−1Step 3. Computedk=∇gxk−1rkStep 4. Updatexk+1=xk+1+dk


## 4. Numerical Experiments

In this section, we provide some numerical experiments to illustrate the feasibility and efficiency of our algorithms proposed in the previous section. All codes were implemented in MATLAB (version 2016a) with a machine precision1.0e−16 on a personal computer with Inter(R) Core(TM) i5-4200M@2.5GHz, and 4.00 G memory. All tensor-vector multiplications were computed turning to the Tensor Toolbox Version 2.6 (http://www.sandia.gov/tgkolda/TensorToolbox/index-2.6.html).For the convenience of expression, we use “IT” and “CPU” to denote the number of iterations and the elapsed CPU time (in seconds), respectively. Let “ERR” and “RES” represent, respectively, two symbols defined by(19)ERR≔x¯k+1−x¯k,RES≔Ax¯km−1−b,where x¯k stands for the kth iteration value. For consistent multilinear systems, the stopping rule is that RES<ɛ≔1.0e−03, or that the number of iterations exceeds the maximum kmax=10000. On the other hand, for the inconsistent case, the terminated criterion is that the error ERR<ɛ or that the number of iterations exceeds the above maximum kmax. In addition, the number of iteration steps, the CPU time, and the residual contained in the tables below are, respectively, the average of 5 runs from different initial iteration vectors unless otherwise stated. In particular, the symbol “– –” in the tables means that the iteration did not derive the desired solution before the terminated rule is satisfied.Example 1.
Let the tensorA in (1) be given by(20)A=sℐ−ℬ,ℬ=bi1i2,…,im∈ℝm,nwithbi1i2,…,im=tani1+i2+⋯+im,where s=1+α⋅max1≤i≤nℬem−1i with e=onesn,1 and α=0.01, and the right-hand side b is chosen such that Ax∗m−1=b, where x∗=8∗e.It is obvious that the tensorA in this example is a symmetrically nonsingular M-tensor, which follows from Theorem 3.2 of [3] that the multilinear system (1) has the unique solution x∗. We compared our algorithms (i.e., Algorithms 1–5) with the Newton method presented in [3] for the different choices of the order m and dimension n of the tensor A. Notably, the latter is a feasible method only for symmetrically nonsingular M-tensor systems and denoted by S-Newton for simplicity.For randomly chosen initial vectorsx0=z10=randn,1, the corresponding numerical results including the number of iterations, the elapsed CPU time, and the residuals are listed in Tables 1–3. From Table 1, one can observe that the AMM-III takes more iteration steps and CPU time than the others, whose numbers of iteration steps exceed the maximum when m=5:12. As two modifications of the AMM method, the R-AMM-I and R-AMM-II methods generally have fewer numbers of iterations than the original one; the paid price is that both of them cost more CPU time. For the S-Newton method, the number of iteration steps is not necessarily less than those of the AMM method, the R-AMM-I method, and the R-AMM-II method, but the clasped CPU time of which is always less than the others. The numerical results in Table 2 illustrate that one can improve the performance of the R-AMM-III method by adding randomized acceleration. Nevertheless, the improved version still needs more iteration steps and CPU time. Hence, in the sequel, we only discuss the performance of Algorithms 1–3.Table 1
Numerical results for the multilinear systems in Example1.
AlgorithmsAMMR-AMM-IR-AMM-IIAMM-IIIS-Newton3,2IT1010101010CPU0.05950.05870.06910.05880.0098RES1.3859e − 061.3859e − 061.3859e − 061.3859e − 061.3859e − 064,2IT1011609069158CPU1.37402.19481.35680.94040.2814RES1.2806e − 041.2828e − 041.2785e − 052.4082e − 046.6883e − 065,2IT353535—35CPU0.81290.81110.8509—0.1036RES5.5163e − 075.5165e − 075.5167e − 07—5.5165e − 076,2IT149199204—195CPU5.43327.26267.6269—0.7267RES7.0794e − 053.6912e − 062.4473e − 05—7.2131e − 067,2IT757575—75CPU3.84043.84133.8947—0.3609RES2.1842e − 072.6061e − 071.7296e − 07—1.7407e − 078,2IT172216404—263CPU12.282915.567929.4014—1.7491RES2.5257e − 047.2879e − 057.0054e − 05—6.3380e − 059,2IT136136136—136CPU12.485612.346412.4625—1.0673RES5.1420e − 055.2800e − 054.6484e − 05—2.7729e − 0510,2IT458299374—288CPU52.889834.545943.4561—2.6043RES5.0163e − 045.9289e − 044.2793e − 04—2.4541e − 0411,2IT238238238—238CPU33.905133.927734.1376—2.4006RES1.3987e − 042.5207e − 044.9489e − 04—1.2302e − 0412,2IT222913733103——CPU403253558——RES0.0000e + 000.0000e + 000.0000e + 00——Table 2
Numerical results for the multilinear systems in Example1.
Algorithms[m, n]4,104,204,404,60AMMIT115192228725CPU1.11972.827917.4153290.0682RES6.4585E − 061.1342E − 076.4105E − 059.7471E − 05R-AMM-IIT123401227494CPU1.21665.962717.4263196.4768RES3.1021E − 053.7321E − 053.3234E − 053.5282E − 05R-AMM-IIIT603074741027CPU0.76066.236437.8383433.1855RES4.7718E − 051.9313E − 041.7141E − 045.0952E − 05AMM-IIIIT– –– –– –– –CPU– –– –– –– –RES– –– –– –– –AMM-III-IIT1583761896553CPU1.62605.5448145.3878258.5192RES5.6223E − 066.7505E − 051.1723E − 041.9261E − 06AMM-III-IIIT989441431929CPU1.179719.9625114.3256386.6988RES4.0573E − 053.9176E − 042.8023E − 067.1267E − 06S-NewtonIT95412490653CPU0.22101.42159.772168.4803RES1.9016E − 043.1456E − 053.3576E − 071.9708E − 04Table 3
Numerical results for the multilinear systems in Example1.
AlgorithmsAMMR-AMM-IR-AMM-IIS-NewtonSORADMM-type3,10IT1010101060—CPU0.07980.08370.11090.06530.977663.2002RES2.0627e − 102.0627e − 101.9949e − 102.0249e − 109.1981e − 046.5422e − 033,30IT1414141443—CPU0.05430.05810.10470.05280.567055.3239RES1.4074e − 091.4074e − 091.3932e − 091.3854e − 097.6229e − 042.9197e − 033,50IT1414141445—CPU0.06880.06960.11180.06311.178959.9732RES5.5115e − 065.5115e − 065.5115e − 065.5114e − 069.9423e − 042.8297e − 034,10IT656060619758CPU0.42380.36710.45550.33840.45630.3600RES6.4552e − 091.0080e − 064.0921e − 065.8740e − 079.4153e − 041.9697e − 044,30IT1022438348387—CPU2.14584.99191.92809.61931.4587219.0979RES1.4629e − 076.6867e − 041.2782e − 078.0313e − 088.2593e − 041.7089e − 024,50IT411141815375190478CPU3.828110.1056165.991831.29889.780640.5234RES5.7258e − 072.2322e − 056.8631e − 078.2008e − 079.1916e − 042.3125e − 055,10IT373737377536CPU0.47770.47370.54190.44260.56110.4755RES5.9909e − 076.0299e − 075.3983e − 075.5564e − 079.9197e − 046.9649e − 065,20IT4949494912948CPU3.32403.32233.38053.31774.16113.5067RES5.8495e − 076.7395e − 077.9920e − 076.1547e − 079.7748e − 041.4438e − 065,30IT4141414116140CPU17.573617.247017.370817.302230.444917.4773RES3.8089e − 064.1881e − 063.9757e − 062.3715e − 069.4562e − 042.9135e − 066,10IT59489558956244205652CPU28.533342.4125303.612229.63897.7784273.8786RES3.1993e − 055.2957e − 043.0382e − 052.4495e − 059.5355e − 045.0279e − 05Since the ADMM-type method proposed by Liang et al. in [17] is an alternating iterative method and the SOR method presented in [26] is an extension of the well-known SOR method for solving linear systems, we compared these two methods with the proposed ones here; see Table 3 for details. This table shows that the AMM and its two randomized versions as well as the S-Newton method have better performance than the SOR method in the number of iteration steps and the elapsed CPU time, while the ADMM-type method takes less iteration steps sometimes for larger m.In order to show the feature of Algorithms1–3, and 5, we, respectively, plotted the figure regarding the residual versus the iteration step k in Figure 1 for m=4,8 and n=2 with initial iterative vectors x0=1,2T. From these figures, one can see that the convergence of these algorithms is nonlinear.Figure 1
Convergence behavior of five algorithms for the multilinear system from Example1.
(a)(b)Next, the proposed iterative methods of this paper will be applied to two numerical examples deriving from realistic applications.Example 2.
Consider the numerical solution of the Klein–Gordon equation [2, 27]:(21)uxm−2⋅△ux=−fx,inΩ,ux=gx,on∂Ω,where △=∑k=0d∂2/∂xk2, Ω=0,1d, and m=3,4,….Whend=1 and the above Klein–Gordon equation can be discretized as the following multilinear system [3]:(22)ℒhum−1=f,in which h=1/n−1 and ℒh=ℒhi1i2,…,im∈ℝm,n with the entries defined by(23)ℒh11,…,1=ℒhnn,…,n=1h2,ℒhii,…,i=2h2,i=2,3,…,n−1,ℒhii−1i,…,i=ℒhiii−1i,…,i=⋯=ℒhi,…,ii−1=−1h2m−1,i=2,3,…,n−1,ℒhii+1i,…,i=ℒhiii+1i,…,i=⋯=ℒhi,…,ii+1=−1h2m−1,i=2,3,…,n−1.Letf be the vector such that ℒhu∗m−1=f with u∗=8∗onesn,1.Since the coefficient tensorA (22) in this example is a nonsymmetric and nonsingular M-tensor, we compared the three proposed iterative algorithms with the NS-Newton method, the SOR method, and the ADMM-type method for initial iterative vector x0 chosen randomly and reported the numerical results in Table 4. This table reflects that the numbers of iterations steps of the three AMM-type methods proposed in this article and the NS-Newton method are identical and the CPU time is almost the same. Relatively, the SOR method has better performance than the ADMM-type method.Table 4
Numerical results for the multilinear systems in Example2.
AlgorithmsAMMR-AMM-IR-AMM-IINS-NewtonSORADMM-type3,10IT8888423—CPU0.15171.30870.19810.08691.205228.4212RES1.9308e − 041.9308e − 041.9308e − 041.9308e − 049.7660e − 041.1703e − 023,30IT15151 5152955—CPU0.06110.05940.07730.052412.379442.0934RES5.9708e − 045.9708e − 046.1631e − 046.0411e − 049.9562e − 041.2772e − 033,50IT16161616——CPU0.07140.06870.09970.064050.058498.4823RES5.3146e − 045.3146e − 045.2894e − 045.3783e − 042.6927e − 021.0036e − 034,10IT19191919705—CPU0.16410.12490.15300.12043.2752147.3018RES7.7367e − 047.7367e − 047.7367e − 047.7367e − 049.8150e − 043.9236e − 024,30IT191919194635335CPU0.36440.35750.37480.341675.407210.5496RES6.7312e − 046.7312e − 046.7312e − 046.7312e − 049.9853e − 049.9831e − 044,50IT252525257856885CPU1.20061.14781.18881.1115293.373154.6492RES8.8715e − 048.8715e − 048.8715e − 048.8715e − 049.9992e − 049.9910e − 045,10IT1919191911553709CPU0.22830.22000.25050.19997.225750.8193RES4.0658e − 044.0658e − 044.0658e − 044.0658e − 049.8803e − 049.9987e − 045,20IT202020203314—CPU1.32221.30611.39441.2267116.2169768.1777RES8.7162e − 048.7162e − 048.7162e − 048.7162e − 049.9932e − 042.2657e − 035,30IT2222222288973175CPU9.06249.43469.66499.33471549.56681406.9804RES4.7169e − 044.7169e − 044.7169e − 044.7169e − 049.9987e − 049.9994e − 046,10IT2222222215137835CPU0.85280.84720.85410.822924.7529361.4588RES6.8396e − 046.8396e − 046.8396e − 046.8396e − 049.9128e − 049.9992e − 04Example 3.
Consider the differential equation(24)−maxγ,λ∈Γ,ΛLλU−ηU−12αγ2U+βγ=0,onΩ,U=g,on∂Ω,where LλUx=1/2σx,λ2U″x+μx,λU′x, Ω=0,1, Γ=0,∞ and Λ is a compact metric space.By applying the “optimize then discretize” approach to the differential equation (24), it can be discretized as a 3rd-order Bellman equation [8]:(25)maxλ∈ΛΔxAλu2=b,in which u=ui∈ℝn, ui≈UiΔx with Δx=1/n for i=0,1,2,…,n, Aλ is a 3rd-order and n+1-dimensional parameterized tensor, and the entries of them are defined as follows:(26)Ai,i−1,iλ=Ai,i,i−1λ,2Ai,i,i−1λ=−12σi2λi1Δx2+μiλi1Δx1−∞,0μiλi,Ai,i,iλ=12σi2λi2Δx2+μiλi1Δx+ηi,2Ai,i,i+1λ=−12σi2λi1Δx2−μiλi1Δx10,+∞μiλi,Ai,i+1,iλ=Ai,i,i+1λ,i=1,2,…,n−1,and Ai,i,iλ=1, if i=0,n. The vector b=bi is defined by(27)bi=12βi2αi,ifi=1,2,…,n−1,gi2,ifi=0,n.In this example, we letσx,λ=0.2, αx=2−x, ηx=0.04, ux,λ=0.04λ, βx=1+x, gx=1, and λ=−1. The coefficient tensor A here is a nonsymmetric and nonsingular M-tensor. We compared the three proposed iterative algorithms with the NS-Newton method, the SOR method, and the ADMM-type method for the initial iterative vector x0 chosen randomly and reported the numerical results in Table 5. From this table, one can see that the three proposed iteration methods in the present paper as well as the NS-Newton method have almost the same convergence; they take fewer iteration steps and CPU time than the SOR method and the ADMM-type method. Notably, the SOR method needs more iteration steps than the ADMM-type method.Table 5
Numerical results for the multilinear systems in Example3.
Algorithms[m, n]3,103,303,503,1003,200AMMIT78101411CPU0.01940.02890.04150.12810.5046RES2.8217e − 054.1736e − 055.7067e − 063.6582e − 061.8151e − 04R-AMM-IIT78101411CPU0.01970.02500.03810.14650.4947RES2.8217e − 054.1736e − 055.7067e − 063.6582e − 061.8151e − 04R-AMM-IIIT78101411CPU0.02740.03500.05040.15140.5053RES2.8217e − 054.1736e − 055.7067e − 063.6582e − 061.8151e − 04NS-NewtonIT78101411CPU0.01520.02210.03380.13320.4939RES2.8217e − 054.1736e − 055.7067e − 063.6582e − 061.8151e − 04SORIT7484324339643– –CPU1.43853.172710.7665118.9747– –RES9.1750e − 049.9189e − 049.9983e − 049.9915e − 04– –ADMM-typeIT301186277400– –CPU0.09080.43983.012789.5569– –RES8.6000e − 049.9739e − 049.9547e − 049.9987e − 04– –Example 4.
Let the tensorA in (1) be given by A=tenrandv∈ℝm,n with v=ones1,m∗n, and let the right-hand side b be chosen such that Ax∗m−1=b, where tenrandv returns an mth-order n-dimensional tensor containing pseudorandom values drawn from the uniform distribution and x∗ is the same as that defined in Example 1.In this example, the tensorA is a randomly generated tensor; then, the solution of the multilinear system (1) may not be unique. We compared the proposed algorithms (i.e., Algorithms 1–3) with the NS-Newton method and the ADMM-type method for randomly produced initial vectors x0=z10=randn,1. The obtained results were displayed in Table 6, which indicates that all the tested methods converge before the number of iteration steps arrives at the maximum kmax. Particularly, the randomized technique may improve the performance of the AMM method. Besides, in Figure 2, we also displayed the convergence behavior of Algorithms 1–3 and 5 when m=4,8 with n=2 and initial iterative vector x0=1,2T.Table 6
Numerical results for the multilinear systems in Example4.
AlgorithmsAMMR-AMM-IR-AMM-IINS-Newton3,2IT1111108CPU0.06740.06360.06940.0231RES1.2805e − 043.8136e − 045.2323e − 041.7398e − 044,2IT34748832CPU0.45611.01081.31560.2187RES7.0644e − 046.6789e − 046.9935e − 043.0967e − 055,2IT75617627CPU1.73561.40591.83330.3656RES4.8703e − 044.7376e − 046.6511e − 041.0511e − 056,2IT10714920550CPU3.84855.34457.55211.0975RES8.9606e − 044.9813e − 045.5719e − 046.6814e − 067,2IT7911611567CPU4.11146.05896.14782.1800RES3.5437e − 045.9663e − 043.7370e − 041.4862e − 048,2IT192175215108CPU13.660212.406815.53094.8704RES5.9745e − 047.8879e − 046.9943e − 048.1033e − 059,2IT331494467150CPU31.442148.226945.93298.9267RES8.2688e − 046.4713e − 048.5859e − 047.9815e − 0810,2IT247446430157CPU29.216452.711551.343511.9918RES8.6464e − 047.6783e − 047.1800e − 045.6479e − 0511,2IT223478509181CPU33.290271.521976.916417.1811RES6.0747e − 046.3177e − 044.8267e − 041.2165e − 0512,2IT617181710102234CPU112.0114333.9998188.5232261.1745RES4.8158e − 046.4279e − 045.7918e − 046.9879e − 04Figure 2
Convergence behavior of four algorithms for the multilinear systems in Example4.
(a)(b)From the above numerical tests, we claim that Algorithm1 is an efficient method for solving some consistent multilinear system (1); what is more, the first two modified versions of this algorithm are also feasible. In particular, one of the advantages of these algorithms is that they can be carried out in a parallel version, which indicates that our algorithms may be implementable for large-scale tensor systems, while the Newton type methods are not easy to run in a parallel version.Additionally, we may always meet the inconsistent multilinear system because of errors. In the next example, we consider the least-squares problem associated with the multilinear system (1); there, the coefficient tensor A and the right-hand side b are chosen randomly.Example 5.
Let the tensorA in (1) be given by A=tenrandv∈ℝm,n with v=ones1,m∗n, and the right-hand side b=randnn,1∈ℝn, where tenrandnv returns an mth-order n-dimensional tensor containing pseudorandom values drawn from the standard normal distribution.We consider the least-squares problem of the multilinear system (1), since the tensor A and the right-hand side b are random, which may lead to an inconsistent multilinear system. In view of the bad performance of Algorithm 4, we only compare Algorithms 1–3 and 5 with the SOR method [26] and the ADMM-type method [17].For the randomly chosen initial vectorx0=randn,1, we still run 5 trials of the algorithms in the cases that m=3:3:18 with n=2, respectively; the numerical results are listed in Table 7. In the tests, the SOR method is not convergent before the number of iteration steps arrives at kmax, and so it is not contained in Table 7. From this table, we can observe that all the tested algorithms are convergent. Relatively, the NS-Newton method often possesses the least number of iterations for the lower m, while the ADMM-type method has better performance than the others for the higher case m. Moreover, the ADMM-type method takes less CPU time than the former when their iteration steps are identical. On the whole, the original AMM method performs better than the two random versions.Table 7
Numerical results for the multilinear systems with fixedn=2 from Example 5.
AlgorithmsAMMR-AMM-IR-AMM-IINS-NewtonADMM-type3,2IT917185135CPU0.05950.11140.14140.02360.5301ERR1.7753e − 049.2725e − 051.4149e − 046.0588e − 049.9810e − 046,2IT4055541578CPU1.47792.02482.03270.53751.3020ERR2.5086e − 021.0281e − 021.0412e − 022.7398e − 039.9766e − 049,2IT791882153435CPU7.448217.599520.98233.12331.3917ERR1.6793e − 035.0516e − 031.5765e − 031.1789e − 039.9366e − 0412,2IT1591543134647CPU28.464427.737356.43538.25073.3506ERR3.3652e − 021.6893e − 029.4988e − 034.5323e − 039.5194e − 0415,2IT3578215527842CPU10.8379241.7929480.695923.74976.7673ERR1.9348e − 024.4889e − 023.5673e − 022.8820e − 032.5703e − 0418,2IT2122919610459CPU143.4707195.262064.584970.319713.3903ERR1.8530e − 026.7065e − 032.2555e − 022.3766e − 039.9950e − 04Next, we present an example in which the entries of the coefficient tensor vary greatly, which means that it leads to ill-posed multilinear systems.Example 6.
Let the tensorA=tensor1:nm,v with v=ones1,m∗n, and let the vector b be chosen in the same way as in Example 2, where tensor1:nm,v reshapes the numbers 1,2,…,nm as an mth-order n-dimensional tensor.Form=5 and n=2, we, respectively, ran Algorithms 2 and 3 for randomly chosen 1000 initial vectors x0=randn,1 and displayed the numerical results in Figure 3, in which the x-axis refers to the kth initial vector and the y-axis represents the difference of the numbers of iteration steps (left) or the elapsed CPU time (right) of the two algorithms. Specifically, the red symbol “+” means that the number of iteration steps or the CPU time of Algorithm 2 is more than that of Algorithm 3, and the green symbol “∗” means that their iteration steps are identical. Otherwise, we use the blue cycle “°.” In Figure 3, the red symbols account for 49.6% of the total (left) and 39.8% of the total (right). Moreover, for the blue symbols “°”, there are 48.1% of cases (left) and 60.2% of cases (right). We should mention that the number of iteration steps of Algorithm 1 always exceeds the kmax for the chosen initial vectors.Figure 3
Comparisons of the iteration steps (a) and the elapsed CPU time (b) between AMM-I and AMM-II for the multilinear system from Example6.
(a)(b)The obtained numerical results indicate that all the algorithms introduced in this paper are feasible and efficient for the tested multilinear systems. In particular, the randomized versions always have better performance than the original one. Notably, many others trials not listed here also reflect that our algorithms are commonly comparable to the existing ones when solving some structured multilinear systems. However, there are still many problems to be studied; for instance, although the proposed algorithms are numerically effective, unfortunately, we could not give their proofs of convergence at present. Moreover, for very largem, our algorithms may be intractable due to the curse-of-dimensionality [28]; thus, how to break the curse is an interesting but important work.


## 5. Conclusions

In this paper, we introduced the alternating minimal method as well as its randomized versions for solving the least-squares problem corresponding to the multilinear systems (1), that is, Algorithms 3.1–3.4. The performed numerical experiments reflect that the proposed algorithms are feasible and efficient. However, we cannot provide the convergence proofs of the two iterative algorithms and leave it as an open problem, which will be considered in the future work.



---

*Source: 6629243-2021-07-13.xml*