A homotopy method for solving multilinear systems with M-tensors

TL;DR

This paper introduces a homotopy method designed to efficiently find the unique positive solution of multilinear systems involving nonsingular M-tensors, with proven convergence and demonstrated numerical effectiveness.

Contribution

The paper presents a novel homotopy approach specifically tailored for solving multilinear systems with M-tensors, including convergence analysis and numerical validation.

Findings

The method converges to the unique positive solution.

Numerical experiments confirm the effectiveness of the approach.

The approach is applicable to systems arising in various applications.

Abstract

Multilinear systems of equations arise in various applications, such as numerical partial differential equations, data mining, and tensor complementarity problems. In this paper, we propose a homotopy method for finding the unique positive solution to a multilinear system with a nonsingular M-tensor and a positive right side vector. We analyze the method and prove its convergence to the desired solution. We report some numerical results based on an implementation of the proposed method using a prediction-correction approach for path following.

Tables

Table 1. Table 1 : Numerical Results

$(m,n)$	euitr	nwitr	time	residue
(3,10)	5	11	0.098	$3.8531\times {10}^{-15}$
(3,50)	5	10	0.126	$3.9837\times {10}^{-14}$
(3,100)	5	9	0.289	$5.8008\times {10}^{-13}$
(3,200)	5	8	0.929	$2.9054\times {10}^{-11}$
(3,400)	5	7	8.099	$1.3199\times {10}^{-8}$
(4,10)	5	10	0.134	$2.2176\times {10}^{-12}$
(4,50)	5	8	1.019	$1.6437\times {10}^{-8}$
(4,80)	5	8	8.902	$3.3791\times {10}^{-11}$
(4,100)	5	8	19.423	$1.3962\times {10}^{-10}$
(5,10)	5	9	0.165	$7.0535\times {10}^{-13}$
(5,20)	5	10	1.646	$2.4450\times {10}^{-13}$
(5,40)	5	9	55.656	$5.4560\times {10}^{-10}$
(6,5)	5	10	0.242	$1.3524\times {10}^{-12}$
(6,10)	5	9	1.483	$3.1020\times {10}^{-8}$
(6,15)	5	10	31.232	$6.2806\times {10}^{-11}$

Equations

$$\mathcal{A}{x}^{m-1}=b,$$

$$(\mathcal{A}{x}^{m-1}{)}_i=\sum _{i_2,\cdots ,i_m=1}^n{A}_{i{i}_2\cdots i_m}{x}_{i_2}\cdots x_{i_m},$$

$$\mathcal{A}{x}^{m-1}=\lambda x^{[m-1]},$$

$$\rho (\mathcal{A})=\max \{|\lambda |:\lambda \mathrm{is}\mathrm{an}\mathrm{eigenvalue}\mathrm{of}\mathcal{A}.\}.$$

$$P(x)=\mathcal{A}{x}^{m-1}-b=0.$$

$$Q(x)=\mathcal{I}{x}^{m-1}-b=0,$$

$$H(x,t)=(1-t)Q(x)+tP(x)=0,t\in [0,1].$$

$$x_0=[b_1^{1/(m-1)},b_2^{1/(m-1)},\dots ,b_n^{1/(m-1)}{]}^T.$$

$$H(x,t)=(t\mathcal{A}+(1-t)\mathcal{I})x^{m-1}-b.$$

$${\hat{A}}_{i_1{i}_2\dots i_m}=\frac{1}{(m-1)!}\sum _{\pi }{A}_{i_1\pi (i_2\dots i_m)},$$

$$D_x\mathcal{A}{x}^{m-1}=(m-1)\hat{\mathcal{A}}{x}^{m-2}.$$

$$D_{x}H(x,t)=(m-1)(t\hat{\mathcal{A}}+(1-t)\mathcal{I})x^{m-2},$$

$$D_{t}H(x,t)=(\mathcal{A}-\mathcal{I})x^{m-1},$$

$$D_{x}H(x(t),t)=(m-1)(t\hat{\mathcal{A}}+(1-t)\mathcal{I})x(t{)}^{m-2}$$

$$t\mathcal{A}+(1-t)\mathcal{I}=(st+1-t)\mathcal{I}-t\mathcal{B}.$$

$$\tau_{0}=\left\{\begin{array}[]{l}\frac{s-\rho({\mathcal{B}})}{\rho({\mathcal{B}})-s+2},\ \ {\rm if}\ \ \rho({\mathcal{B}})-s+2>0,\\ 1,\hskip 43.36243pt{\rm if}\ \ \rho({\mathcal{B}})-s+2\leq 0.\end{array}\right.$$

$$t\mathcal{A}+(1-t)\mathcal{I}=t\left(\frac{st+1-t}{t}\mathcal{I}-\mathcal{B}\right)$$

$$[(t\hat{\mathcal{A}}+(1-t)\mathcal{I})x(t{)}^{m-2}]x(t)=(t\hat{\mathcal{A}}+(1-t)\mathcal{I})x(t{)}^{m-1}=(t\mathcal{A}+(1-t)\mathcal{I})x(t{)}^{m-1}=b$$

$$D_{x}H(x(t),t)\cdot \frac{dx}{dt}=-D_{t}H(x(t),t).$$

$$\bar{\mathcal{A}}{x}^{m-1}=\bar{b},$$

$$\Vert \bar{\mathcal{A}}x(1{)}^{m-1}-\bar{b}{\Vert }_2\le 10^{-12}.$$

$$s=(1+ϵ)\cdot \underset{1\le i\le n}{\max }\sum _{i_2,\dots ,i_m}{\mathcal{B}}_{i,i_2,\dots ,i_m},$$

$$\rho (\mathcal{B})\le \underset{1\le i\le n}{\max }\sum _{i_2,\dots ,i_m}{\mathcal{B}}_{i,i_2,\dots ,i_m}.$$

Full text

Key words. M-tensor, multilinear system, homotopy method.

AMS subject classification (2010). 65H10, 65H20.

1 Introduction

Let $\mathbb{R}$ and $\mathbb{C}$ be the real field and complex field, respectively. We denote the set of all $m$th-order, $n$-dimensional real tensors by $\mathbb{R}^{[m,n]}$. For a tensor ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ and a vector $b\in\mathbb{R}^{n}$, a multilinear system is defined as

[TABLE]

where $x=[x_{1},x_{2},\ldots,x_{n}]^{T}$ is the unknown vector, and ${\mathcal{A}}x^{m-1}$ denotes the column vector whose $i$th entry is

[TABLE]

for $i=1,2,\ldots,n$. Multilinear systems of the form (1.1) arise in a number of applications, such as numerical partial differential equations, data mining, and tensor complementarity problems (see for example, [6, 7, 10]).

In their pioneering works, Qi [11] and Lim [9] independently introduced the concept of tensor eigenvalues. We say that $(\lambda,x)\in\mathbb{C}\times\mathbb{C}^{n}\backslash\{0\}$ is an eigenpair of a tensor ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ if

[TABLE]

where $x^{[m-1]}=[x_{1}^{m-1},x_{2}^{m-1},\ldots,x_{n}^{m-1}]^{T}$. Let $\rho({\mathcal{A}})$ denote the spectral radius of tensor ${\mathcal{A}}$, that is,

[TABLE]

Recently, the notion of M-tensors has been proposed and their properties have been studied in [5, 13]. A tensor ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ is called an M-tensor if it can be written as ${\mathcal{A}}=s{\mathcal{I}}-B$, in which ${\mathcal{I}}$ is the $m$th-order, $n$-dimensional identity tensor, ${\mathcal{B}}$ is a nonnegative tensor (that is, each entry of ${\mathcal{B}}$ is nonnegative), and $s\geq\rho({\mathcal{B}})$. Furthermore, ${\mathcal{A}}$ is called a nonsingular M-tensor if $s>\rho({\mathcal{B}})$.

In [6], Ding and Wei investigated the solutions of the multilinear system (1.1) when the coefficient tensor ${\mathcal{A}}$ is an M-tensor. In particular, they show that the system (1.1) has a unique positive solution if ${\mathcal{A}}$ a nonsingular M-tensor and $b$ is a positive vector (see [6, Theorem 3.2]). They generalized the Jacobi and Gauss-Seidel methods for linear systems to find the unique positive solution of the multilinear system (1.1). They also proposed to use Newton’s method when the nonsingular M-tensor ${\mathcal{A}}$ is symmetric and numerically showed that Newton’s method is much faster than the other methods. However, it is unclear whether or not the Newton method proposed in [6] always works when ${\mathcal{A}}$ is not symmetric.

In this paper, we propose a homotopy method for finding the unique positive solution of the multilinear system (1.1) and prove its convergence. The homotopy method is implemented using an Euler-Newton prediction-correction approach for path tracking. Numerical experiments show the efficiency of our method.

The paper is organized as follows. We introduce our homotopy method and prove its convergence in Section 2. Then we give some numerical results in Section 3.

2 A Homotopy Method

We are to design a homotopy method for finding the unique positive solution of the system (1.1) when ${\mathcal{A}}$ is a nonsingular M-tensor and $b$ is a positive vector. For this purpose, we will solve the following polynomial system:

[TABLE]

We choose the starting system

[TABLE]

and construct the following homotopy

[TABLE]

Note that the starting system (2.2) trivially has a unique positive solution

[TABLE]

Moreover, the homotopy $H(x,t)$ can be expressed as

[TABLE]

The partial derivatives matrix $D_{x}H(x,t)$ of the homotopy $H(x,t)$ plays an important role in our method. To compute this matrix, we need to partially symmetrize tensor ${\mathcal{A}}=(A_{i_{1},i_{2},\ldots,i_{m}})$ with respect to the indices $i_{2},\ldots,i_{m}$. Specifically, we define the partially symmetrized tensor $\hat{{\mathcal{A}}}=(\hat{A}_{i_{1},i_{2},\ldots,i_{m}})$ by

[TABLE]

where the sum is over all the permutations $\pi(i_{2}\ldots i_{m})$. The following lemma shows that this partial symmetrization preserves the nonsingular M-tensor structure.

LEMMA 2.1

If ${\mathcal{A}}\in{\mathcal{R}}^{[m,n]}$ is a nonsingular M-tensor, so is $\hat{{\mathcal{A}}}$.

Proof**: Since ${\mathcal{A}}$ is an M-tensor, there is a nonnegative tensor ${\mathcal{B}}$ such that ${\mathcal{A}}=s{\mathcal{I}}-{\mathcal{B}}$. Then $\hat{{\mathcal{B}}}$ is nonnegative and moreover, $\hat{{\mathcal{A}}}=s{\mathcal{I}}-\hat{{\mathcal{B}}}$. According to [5, Page 3277, Conditions D1 and D4], ${\mathcal{A}}$ is a nonsingular M-tensor if and only if there is a positive vector $y\in\mathbb{R}$ such that ${\mathcal{A}}y^{m-1}$ is a positive vector. Note that ${\mathcal{A}}x^{m-1}=\hat{{\mathcal{A}}}x^{m-1}$ for all $x\in\mathbb{R}^{n}$. Thus, $\hat{{\mathcal{A}}}y^{m-1}$ is a positive vector. Using the results in [5] again, $\hat{{\mathcal{A}}}$ is a nonsingular M-tensor. $\Box$

**

The partial derivatives matrix of ${\mathcal{A}}x^{m-1}$ with respect to $x$ is

[TABLE]

Therefore, the partial derivatives of $H$ with respect to $x$ and $t$ are:

[TABLE]

and

[TABLE]

respectively.

THEOREM 2.1

*Suppose that ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ is a nonsingular M-tensor and $b\in\mathbb{R}^{n}$ is a positive vector. Then there exists a number $\tau_{0}>0$ such that, for each $t\in[0,1+\tau_{0})$,

(a) $H(x,t)=0$ has a unique positive solution $x(t)$;

(b) the partial derivatives matrix*

[TABLE]

is nonsingular.

Proof**: Let ${\mathcal{A}}=s{\mathcal{I}}-{\mathcal{B}}$ such that ${\mathcal{B}}$ is a nonnegative tensor and $s>\rho({\mathcal{B}})$. Note that**

[TABLE]

Clearly, this tensor is a nonsingular M-tensor for $0\leq t\leq 1$. Choose

[TABLE]

Then $\tau_{0}>0$ and $\displaystyle{\frac{st+1-t}{t}>\frac{s+\rho({\mathcal{B}})}{2}>\rho({\mathcal{B}})}$ for $1\leq t<1+\tau_{0}$. This implies that

[TABLE]

is a nonsingular M-tensor when $1\leq t<1+\tau_{0}$. Therefore, $t{\mathcal{A}}+(1-t){\mathcal{I}}$ is a nonsingular M-tensor for each $t\in[0,1+\tau_{0})$. It follows that $H(x,t)=(t{\mathcal{A}}+(1-t){\mathcal{I}})x^{m-1}-b=0$ has a unique positive solution $x(t)$ for each $t\in[0,1+\tau_{0})$ by [6, Theorem 3.2]. Moreover, $t\hat{{\mathcal{A}}}+(1-t){\mathcal{I}}$ is a nonsingular M-tensor for each $t\in[0,1+\tau_{0})$ by Lemma 2.1. Note that the matrix $(t\hat{{\mathcal{A}}}+(1-t){\mathcal{I}})x(t)^{m-2}$ is Z-matrix, and

[TABLE]

**is a positive vector, we must have that $D_{x}H(x(t),t)=(m-1)(t\hat{{\mathcal{A}}}+(1-t){\mathcal{I}})x(t)^{m-2}$ is nonsingular by [3, Chapter 6]. $\Box$

**

COROLLARY 2.1

Suppose that ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ is a nonsingular M-tensor and $b\in\mathbb{R}^{n}$ is a positive vector. Then the positive solutions $x(t)$ of $H(x,t)=0$ for $t\in[0,1+\tau_{0})$ form a smooth curve in $\mathbb{R}_{++}^{n}$, where $\mathbb{R}_{++}^{n}$ is the set of positive $n$-vectors.

Proof**: By Theorem 2.1, $H(x,t)=0$ has a unique positive solution $x(t)$ for each $t\in[0,1+\tau_{0})$. The conclusion follows by using the Implicit Function Theorem and a continuation argument ([8]). $\Box$

**

Our next theorem shows that the homotopy (2.3) works.

THEOREM 2.2

Suppose that ${\mathcal{A}}\in\mathbb{R}^{[m,n]}$ is a nonsingular M-tensor and $b\in\mathbb{R}^{n}$ is a positive vector. Starting from the initial $x(0)=x_{0}=[b_{1}^{1/(m-1)},b_{2}^{1/(m-1)},\ldots,b_{n}^{1/(m-1)}]^{T}$, let $x(t)$ be the solution curve obtained by solving the homotopy $H(x,t)=0$ in $\mathbb{R}_{++}^{n}\times[0,1]$. Then $x(1)$ is the unique positive solution of the system (1.1).

Proof**: According to Corollary 2.1, the positive solutions $x(t)$ of $H(x,t)=0$ for $t\in[0,1+\tau_{0})$ form a smooth curve in $\mathbb{R}_{++}^{n}$. Differentiating $H(x(t),t)=0$ with respect to $t$ gives**

[TABLE]

**Since $D_{x}H(x(t),t)$ is nonsingular for all $t\in[0,1+\tau_{0})$, this system of differential equations is well defined for $t\in[0,1+\tau_{0})$. We can follow the curve by solving this system with the initial condition $x(0)=x_{0}$. Clearly, $x(1)$ is the unique positive solution of the system (1.1). $\Box$

**

We now present our homotopy method for finding the unique positive solution of (1.1).

ALGORITHM 2.1

Finding the positive solution for (1.1) when ${\mathcal{A}}$ is a nonsingular M-tensor and $b$ a positive vector.

Initialization.* Choose initial $x_{0}=[b_{1}^{1/(m-1)},b_{2}^{1/(m-1)},\ldots,b_{n}^{1/(m-1)}]^{T}$.

 

Path following. Solve the differential system (2.8) with the initial condition $x(0)=x_{0}$ in $\mathbb{R}_{++}^{n}$. $x(1)$ is the desired solution for system (1.1).*

3 Numerical Results

We have implemented Algorithm 2.1 in Matlab. The code can be downloaded from:

** http://homepages.umflint.edu/~lxhan/software.html **

In our implementation, an Euler-Newton type predication-correction approach (****[1, 12]****) with an adaptive stepsize for solving the system of differential equations (2.8) with initial condition $x(0)=x_{0}$ is used. Moreover, it solves the scaled system

[TABLE]

where $\bar{{\mathcal{A}}}={\mathcal{A}}/\omega$, $\bar{b}=b/\omega$, and $\omega$ is the largest value among the entries of $b$ and the absolute values of entries of ${\mathcal{A}}$. The code terminates if the residue of the scaled system

[TABLE]

To test the effectiveness of Algorithm 2.1, we did some numerical experiments. All the experiments were done using MATLAB 2014b on a laptop computer with Intel Core i7-4600U at 2.10 GHz and 8 GB memory running Microsoft Windows 7. The tensor toolbox of ****[2]**** was used to compute tensor-vector products and to compute partially symmetrized tensor $\hat{{\mathcal{A}}}$.

The examples we tested are generated by the following method. We chose nonnegative tensor ${\mathcal{B}}\in\mathbb{R}^{[m,n]}$ whose entires are uniformly distributed in $(0,1)$. Set

[TABLE]

for some $\epsilon>0$. Let ${\mathcal{A}}=s{\mathcal{I}}-{\mathcal{B}}$. According to ****[11, 4]****,

[TABLE]

Thus, $s>\rho({\mathcal{B}})$ and ${\mathcal{A}}$ is a nonsingular M-tensor. We chose the right side vector $b\in\mathbb{R}^{n}$ with entires uniformly distributed in $(0,1)$.

We tested the algorithm on tensors of various sizes by choosing different values of $m$ and $n$. We used $\epsilon=0.01$ as in ****[6]****. We now summarize the numerical results in Table 1. In this table, euitr and nwitr denote the number of Euler prediction steps and the total number of Newton correction steps were used, time denotes the CPU time used (in seconds) when the algorithm terminated, and residue denotes the residue $\|{\mathcal{A}}x(1)^{m-1}-b\|_{2}$ of the original system (1.1) at termination.

For the examples we tested, we observe that Algorithm 2.1 can find the positive solution of the multilinear system (1.1) when ${\mathcal{A}}$ is a nonsingular M-tensor and $b$ is a positive vector. It is efficient in terms of both euitr and nwitr. We remark that for the $(4,100),(5,40),(60,15)$ cases, the relatively large CPU time used by Algorithm 2.1 is mainly due to the procedure of partially symmetrizing tensor ${\mathcal{A}}$. A more efficient symmetrization method can help save the CPU time for such cases. Of course, the partial symmetrization is not needed if the tensor ${\mathcal{A}}$ is symmetric. Overall, the numerical results show that Algorithm 2.1 is quite promising.