# Computing the extremal nonnegative solutions of the $M$-tensor equation with
a nonnegative right side vector

Chun-Hua Guo  Department of Mathematics and Statistics, University of Regina,
Regina, Saskatchewan S4S 0A2, Canada<EMAIL_ADDRESS>
###### Abstract.

We consider the tensor equation whose coefficient tensor is a nonsingular
$M$-tensor and whose right side vector is nonnegative. Such a tensor equation
may have a large number of nonnegative solutions. It is already known that the
tensor equation has a maximal nonnegative solution and a minimal nonnegative
solution (called extremal solutions collectively). However, the existing
proofs do not show how the extremal solutions can be computed. The existing
numerical methods can find one of the nonnegative solutions, without knowing
whether the computed solution is an extremal solution. In this paper, we
present new proofs for the existence of extremal solutions. Our proofs are
much shorter than existing ones and more importantly they give numerical
methods that can compute the extremal solutions. Linear convergence of these
numerical methods is also proved under mild assumptions. Some of our
discussions also allow the coefficient tensor to be a $Z$-tensor or allow the
right side vector to have some negative elements.

###### Key words and phrases:

$Z$-tensor; $M$-tensor; Maximal nonnegative solution; Minimal nonnegative
solution; Iterative methods; Linear convergence.

###### 1991 Mathematics Subject Classification:

Primary 65F10; Secondary 15A69, 65H10

##  1\. Introduction

We consider the tensor equation

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

where $x,b\in{\mathbb{R}}^{n}$ and ${\mathcal{A}}$ is a real-valued $m$th-
order $n$-dimensional tensor that has the form

${\mathcal{A}}=(a_{i_{1}i_{2}\ldots i_{m}}),\quad a_{i_{1}i_{2}\ldots
i_{m}}\in{\mathbb{R}},\quad 1\leq i_{1},i_{2},\ldots,i_{m}\leq n,$

and ${\mathcal{A}}x^{m-1}\in{\mathbb{R}}^{n}$ has elements

$({\mathcal{A}}x^{m-1})_{i}=\sum_{i_{2},\ldots,i_{m}=1}^{n}a_{ii_{2}\ldots
i_{m}}x_{i_{2}}\cdots x_{i_{m}},i=1,2,\ldots,n.$

Later on, we will also use $\\{x_{k}\\}$ to denote a sequence of vectors. The
$i$th element of a vector $x$ will then be denoted by $(x)_{i}$ (instead of
$x_{i}$) to avoid confusion. The element $a_{i_{1}i_{2}\ldots i_{m}}$ of the
tensor ${\mathcal{A}}$ will also be denoted by
${\mathcal{A}}_{i_{1}i_{2}\ldots i_{m}}$. The tensor equation (1) is also
called a multilinear system of equations. It appears in many applications
including data mining, numerical partial differential equations, and tensor
complementarity problems; see [11] and the references cited therein.

We denote the set of all real-valued $m$th-order $n$-dimensional tensors by
${\mathbb{R}}^{[m,n]}$. A tensor ${\mathcal{A}}=(a_{i_{1}i_{2}\ldots
i_{m}})\in{\mathbb{R}}^{[m,n]}$ is called semi-symmetric [16] if
$a_{ij_{2}\ldots j_{m}}=a_{ii_{2}\ldots i_{m}}$, $1\leq i\leq n$, $j_{2}\ldots
j_{m}$ is any permutation of $i_{2}\ldots i_{m}$, $1\leq
i_{2},\ldots,i_{m}\leq n$. For any tensor ${\mathcal{A}}=(a_{i_{1}i_{2}\ldots
i_{m}})$, we have a semi-symmetric tensor
$\bar{{\mathcal{A}}}=(\bar{a}_{i_{1}i_{2}\ldots i_{m}})$ defined by

$\bar{a}_{i_{1}i_{2}\ldots i_{m}}=\frac{1}{(m-1)!}\sum_{j_{2}\ldots
j_{m}}a_{i_{1}j_{2}\ldots j_{m}},$

where $j_{2}\ldots j_{m}$ is any permutation of $i_{2}\ldots i_{m}.$ Then
${\mathcal{A}}x^{m-1}=\bar{{\mathcal{A}}}x^{m-1}$ for all
$x\in{\mathbb{R}}^{n}$, and
$({\mathcal{A}}x^{m-1})^{\prime}=(m-1)\bar{{\mathcal{A}}}x^{m-2}$ (see [12,
16]), where for any ${\mathcal{A}}\in{\mathbb{R}}^{[m,n]}$, the $(i,j)$
element of ${\mathcal{A}}x^{m-2}\in{\mathbb{R}}^{n\times n}$ is defined by

$({\mathcal{A}}x^{m-2})_{ij}=\sum_{i_{3},\ldots,i_{m}=1}^{n}a_{iji_{3}\ldots
i_{m}}{x}_{i_{3}}\cdots{x}_{i_{m}}.$

We call $\lambda\in\mathbb{R}$ an eigenvalue and
$x\in{\mathbb{R}}^{n}\setminus\\{0\\}$ a corresponding eigenvector of
${\mathcal{A}}$ if [18]

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

where for any positive real number $r$, $x^{[r]}\in{\mathbb{R}}^{n}$ is given
by $(x^{[r]})_{i}=(x_{i})^{r}$, $i=1,2,\ldots,n$. The spectral radius of
${\mathcal{A}}$ is the maximum modulus of its eigenvalues, and is denoted by
$\rho({\mathcal{A}})$.

Let $[n]=\\{1,2,\ldots,n\\}$. For $x,y\in{\mathbb{R}}^{n}$, we write $x\geq y$
if $x_{i}\geq y_{i}$ for all $i\in[n]$, and write $x>y$ if $x_{i}>y_{i}$ for
all $i\in[n]$. If $x\geq 0$, we say $x$ is nonnegative; if $x>0$ we say $x$ is
positive. The set ${\mathbb{R}}_{+}=\\{x\in{\mathbb{R}}^{n}\ |\ x\geq 0\\}$
will be used frequently. A solution $x$ of (1) is said to be a maximal
nonnegative solution if $x\geq y$ for every solution $y\geq 0$; a solution $x$
of (1) is said to be a minimal nonnegative solution if $0\leq x\leq y$ for
every solution $y\geq 0$. A tensor ${\mathcal{A}}$ is said to be nonnegative,
denoted by ${\mathcal{A}}\geq 0$, if all its elements are nonnegative. The
identity tensor ${\mathcal{I}}\in{\mathbb{R}}^{[m,n]}$ is such that its
diagonal elements are all ones and its off-diagonal elements are all zeros,
i.e., ${\mathcal{I}}_{ii\ldots i}=1$ for $i\in[n]$ and
${\mathcal{I}}_{i_{1}i_{2}\ldots i_{m}}=0$ elsewhere.

######  Definition 1.

[3] A tensor ${\mathcal{A}}\in{\mathbb{R}}^{[m,n]}$ is called a $Z$-tensor if
its off-diagonal elements are nonpositive. If ${\mathcal{A}}$ can be written
as ${\mathcal{A}}=s{\mathcal{I}}-{\mathcal{B}}$ with ${\mathcal{B}}\geq 0$ and
$s>\rho({\mathcal{B}})$, then the tensor ${\mathcal{A}}$ is called a
nonsingular $M$-tensor.

For ${\mathcal{B}}\geq 0$, $\rho({\mathcal{B}})$ can be found by the power
method in [15] and the Newton–Noda iteration in [12]. Thus, one can determine
whether a $Z$-tensor is a nonsingular $M$-tensor using the definition. It is
known [3] that a $Z$-tensor ${\mathcal{A}}$ is a nonsingular
${\mathcal{M}}$-tensor if and only if ${\mathcal{A}}x^{m-1}>0$ for some $x>0$.
It follows that the diagonal elements of a nonsingular ${\mathcal{M}}$-tensor
are all positive.

Suppose ${\mathcal{A}}$ is a nonsingular $M$-tensor. When $b$ is positive,
equation (1) has a unique positive solution [4]. When $b$ is nonnegative,
equation (1) has nonnegative solutions and moreover it has a minimal
nonnegative solution [14] and has a maximal nonnegative solution [9].

The following simple example is an extension of Example 1.1 in [1]. It shows
that the number of nonnegative solutions could be huge when $b\geq 0$.

######  Example 1.

Let $k\geq 1$ be an integer and let
${\mathcal{A}}=(a_{i_{1}i_{2}i_{3}i_{4}})\in{\mathbb{R}}^{[4,2k]}$, where
$a_{iiii}=1$ for $i\in[2k]$, $a_{2i-1,2i-1,2i-1,2i}=-2$ for $i\in[k]$, and all
other elements are zero. It is clear that ${\mathcal{A}}$ is a $Z$-tensor with
${\mathcal{A}}x^{3}>0$ for $x=[3,1,\ldots,3,1]^{T}$. So ${\mathcal{A}}$ is a
nonsingular $M$-tensor. For $b=[0,1,\ldots,0,1]^{T}$, we find that the
equation ${\mathcal{A}}x^{3}=b$ has $2^{k}$ solutions
$x=[x_{1},x_{2},\ldots,x_{2k-1},x_{2k}]^{T}$, where for $i\in[k]$,
$[x_{2i-1},x_{2i}]=[0,1]^{T}$ or $[2,1]^{T}$. The minimal nonnegative solution
is $[0,1,\ldots,0,1]^{T}$ and the maximal nonnegative solution is
$[2,1,\ldots,2,1]^{T}$.

Since equation (1) may have a large number of nonnegative solutions, it is
unlikely that each of these solutions is of practical interest. Intuitively,
the extremal nonnegative solutions are of particular interest. By computing
the maximal nonnegative solution, we can answer the question whether the
equation has a positive solution. By computing the extremal nonnegative
solutions, we get bounds for all other nonnegative solutions.

In [14], the tensor complementarity problem (TCP):

$x\geq 0,\quad{\mathcal{A}}x^{m-1}-b\geq 0,\quad
x^{T}({\mathcal{A}}x^{m-1}-b)=0$

is considered. When ${\mathcal{A}}$ is a $Z$-tensor and
$b\in{\mathbb{R}}_{+}$, it is shown in [14] that the solution set of the TCP
is the same as the set of all nonnegative solutions of equation (1). So a
sparsest solution to the TCP is the minimal nonnegative solution of equation
(1). The problem of finding a sparsest solution to the TCP is of practical
interest [14]. Note that, for this application, the case $b\geq 0$ with some
zero elements is much more interesting than the case $b>0$ (where the equation
has a unique positive solution and there are no other nonnegative solutions).

When ${\mathcal{A}}$ is a nonsingular $M$-tensor and $b>0$, the unique
positive solution $x({\mathcal{A}},b)$ of ${\mathcal{A}}x^{m-1}=b$ is the
solution of particular interest. But some elements of $b$ may be very tiny.
What happens to $x({\mathcal{A}},b)$ if $b$ decreases monotonically towards a
vector $b^{(0)}$ with some zero elements? We will see in the next section that
$x({\mathcal{A}},b)$ converges to the maximal nonnegative solution of the
equation ${\mathcal{A}}x^{m-1}=b^{(0)}$.

In [9], the authors state in the introduction that their purpose is to find
the largest (maximal) nonnegative solution. But their algorithms there only
find a nonnegative solution, which is usually not maximal. Numerical methods
have also been presented in [1, 10, 11]. Those methods can find a nonnegative
solution, but the solution is usually not maximal and cannot be guaranteed to
be minimal.

While our main interest is on equation (1) with ${\mathcal{A}}$ being a
nonsingular $M$-tensor and $b\geq 0$, some of our discussions also allow
${\mathcal{A}}$ to be a $Z$-tensor or allow $b$ to have some negative
elements. In Section 2, we present new proofs for the existence of extremal
nonnegative solutions by using simple iterative methods. These iterative
methods can actually compute the extremal solutions. Some other results also
follow directly from the existence theorems and their proofs. In Section 3, we
show how the simple iterations used in Section 2 can be generalized to have
faster convergence. Linear convergence of these iterative methods is also
proved under mild assumptions. Some concluding remarks are given in Section 4.

##  2\. Existence of extremal nonnegative solutions

The following theorem is a main result in [14] (stated differently in [14];
see [14, Theorem 3] and its proof). The proof there is based on several other
results and does not show how the minimal nonnegative solution can be
computed. Our new proof here is very short and it provides a way to compute
the minimal solution. The approach we take here is similar to the approach we
used many years ago for determining the existence of the minimal nonnegative
solution of $M$-matrix algebraic Riccati equations [5, 6].

######  Theorem 1.

Let ${\mathcal{A}}$ be a $Z$-tensor and $b\in{\mathbb{R}}_{+}$. Suppose that

${\mathcal{S}}_{g}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\geq
b\\}\neq\emptyset.$

Then ${\mathcal{S}}_{g}$ has a minimal element that is also the minimal
nonnegative solution of ${\mathcal{A}}x^{m-1}=b$.

###### Proof.

We write ${\mathcal{A}}={\mathcal{D}}-{\mathcal{B}}$, where ${\mathcal{D}}$ is
a diagonal tensor with positive diagonal elements and ${\mathcal{B}}\geq 0$.
Equation (1) becomes ${\mathcal{D}}x^{m-1}={\mathcal{B}}x^{m-1}+b$. Then we
have the fixed-point iteration, given implicitly as follows:

(2) ${\mathcal{D}}x_{k+1}^{m-1}={\mathcal{B}}x_{k}^{m-1}+b.$

(This is the Jacobi iteration when the diagonal elements of ${\mathcal{A}}$
are positive and ${\mathcal{D}}$ is the diagonal part of ${\mathcal{A}}$.)
Note that ${\mathcal{D}}x_{k+1}^{m-1}=Dx_{k+1}^{[m-1]}$, where $D$ is the
diagonal matrix having the diagonal elements of the tensor ${\mathcal{D}}$ on
the diagonal. When ${\mathcal{B}}x_{k}^{m-1}+b\geq 0$, the iteration can be
given explicitly as

$x_{k+1}=\left(D^{-1}({\mathcal{B}}x_{k}^{m-1}+b)\right)^{[1/(m-1)]},$

but the implicit form will be more convenient for discussions.

Let $x$ be any element in ${\mathcal{S}}_{g}$. Then $x\geq 0$ and
${\mathcal{B}}x^{m-1}+b\leq{\mathcal{D}}x^{m-1}$. Take $x_{0}=0$. We can
generate a sequence $\\{x_{k}\\}$ by iteration (2). It is clear that
$x_{0}\leq x_{1}$. Suppose $x_{k-1}\leq x_{k}$ ($k\geq 1$). Then
${\mathcal{D}}x_{k}^{m-1}={\mathcal{B}}x_{k-1}^{m-1}+b\leq{\mathcal{B}}x_{k}^{m-1}+b={\mathcal{D}}x_{k+1}^{m-1}$.
Thus $x_{k}\leq x_{k+1}$. Therefore, $x_{k}\leq x_{k+1}$ for all $k\geq 0$.
Also, $x_{0}\leq x$. Suppose $x_{k}\leq x$ ($k\geq 0$). Then
${\mathcal{D}}x_{k+1}^{m-1}={\mathcal{B}}x_{k}^{m-1}+b\leq{\mathcal{B}}x^{m-1}+b\leq{\mathcal{D}}x^{m-1}$.
Thus $x_{k+1}\leq x$. Therefore, $x_{k}\leq x$ for all $k\geq 0$. Now,
$\\{x_{k}\\}$ is monotonically increasing and bounded above by $x$. Thus,
$\lim_{k\to\infty}x_{k}=x_{*}$ exists and $x_{*}\leq x$. Letting $k\to\infty$
in (2), we see that $x_{*}$ is a nonnegative solution of (1) and $x_{*}\leq x$
for every $x$ in ${\mathcal{S}}_{g}$. In particular, $x_{*}\leq x$ for every
nonnegative solution $x$ of (1), so $x_{*}$ is the minimal nonnegative
solution of (1), and also the minimal element in ${\mathcal{S}}_{g}$. ∎

######  Corollary 2.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and $b\geq 0$. Then equation
(1) has a minimal nonnegative solution.

###### Proof.

Since ${\mathcal{A}}$ be a nonsingular $M$-tensor, we have
${\mathcal{A}}\hat{x}^{m-1}>0$ for some $\hat{x}>0$. Then for scalar $t>0$
sufficiently large, $\tilde{x}=t\hat{x}>0$ is such that
${\mathcal{A}}\tilde{x}^{m-1}=t^{m-1}{\mathcal{A}}\hat{x}^{m-1}\geq b$. ∎

The next result is already known in [10]. It also follows quickly from our
proof of Theorem 1.

######  Corollary 3.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and $b>0$. Then the unique
positive solution of (1) is the only nonnegative solution.

###### Proof.

From our proof of Theorem 1, we see that the minimal nonnegative solution is
positive in this case. ∎

######  Corollary 4.

Let ${\mathcal{A}}$ be a $Z$-tensor and $b\in{\mathbb{R}}_{+}$. Suppose that

${\mathcal{S}}_{g}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\geq
b\\}\neq\emptyset$

and let $x_{\min}({\mathcal{A}},b)$ be the minimal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$. If any element of $b$ decreases but remains
nonnegative, or if any diagonal element of ${\mathcal{A}}$ increases, or if
any off-diagonal element of ${\mathcal{A}}$ increases but remains nonpositive,
then the new equation $\tilde{{\mathcal{A}}}x^{m-1}=\tilde{b}$ also has a
minimal nonnegative solution $x_{\min}(\tilde{{\mathcal{A}}},\tilde{b})$.
Moreover, $x_{\min}(\tilde{{\mathcal{A}}},\tilde{b})\leq
x_{\min}({\mathcal{A}},b)$.

###### Proof.

Under any of those changes, $\tilde{{\mathcal{A}}}$ is a $Z$-tensor and
$\tilde{b}\in{\mathbb{R}}_{+}$. Let
$\tilde{{\mathcal{S}}}_{g}=\\{x\in{\mathbb{R}}_{+}\ |\
\tilde{{\mathcal{A}}}x^{m-1}\geq\tilde{b}\\}$. Then
${\mathcal{S}}_{g}\subseteq\tilde{{\mathcal{S}}}_{g}$ and the conclusions
follow immediately. ∎

The following theorem has been proved in [9]. The proof there is somewhat
complicated and does not show how the maximal solution can be computed. Here
we present a short proof and also a way to compute the maximal solution.

######  Theorem 5.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and suppose that
${\mathcal{S}}_{l}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\leq
b\\}\neq\emptyset$. Then ${\mathcal{S}}_{l}$ has a maximal element that is
also the maximal nonnegative solution of ${\mathcal{A}}x^{m-1}=b$.

###### Proof.

Let ${\mathcal{D}}$ be the diagonal part of ${\mathcal{A}}$ and write
${\mathcal{A}}={\mathcal{D}}-{\mathcal{B}}$. Then the diagonal elements of
${\mathcal{D}}$ are positive and ${\mathcal{B}}\geq 0$. Equation (1) becomes
${\mathcal{D}}x^{m-1}={\mathcal{B}}x^{m-1}+b$. Then we have the fixed-point
iteration (Jacobi iteration), given implicitly as follows:

(3) ${\mathcal{D}}x_{k+1}^{m-1}={\mathcal{B}}x_{k}^{m-1}+b.$

Let $x$ be any element in ${\mathcal{S}}_{l}$. Then $x\geq 0$ and
${\mathcal{B}}x^{m-1}+b\geq{\mathcal{D}}x^{m-1}$. Take $x_{0}\geq x$ such that
${\mathcal{A}}x_{0}^{m-1}\geq b$, so
${\mathcal{B}}x_{0}^{m-1}+b\leq{\mathcal{D}}x_{0}^{m-1}$. Note that

${\mathcal{B}}x_{0}^{m-1}+b\geq{\mathcal{B}}x^{m-1}+b\geq{\mathcal{D}}x^{m-1}\geq
0.$

Thus $x_{1}$ is determined by iteration (3) and
${\mathcal{D}}x_{1}^{m-1}\geq{\mathcal{D}}x^{m-1}$, so $x_{1}\geq x$. Also,
${\mathcal{D}}x_{1}^{m-1}\leq{\mathcal{D}}x_{0}^{m-1}$, so $x_{1}\leq x_{0}$.
By induction, we can show that $x_{k+1}\leq x_{k}$ and $x_{k}\geq x$ for all
$k\geq 0$. Therefore, $\lim_{k\to\infty}x_{k}=x_{*}$ exists and $x_{*}\geq x$.
Then $x_{*}$ is a nonnegative solution of (1).

We still need to show that we can choose one fixed $x_{0}$ such that
${\mathcal{A}}x_{0}^{m-1}\geq b$ and $x_{0}\geq x$ for all
$x\in{\mathcal{S}}_{l}$. To this end, we take any $\hat{b}>0$ such that
$\hat{b}\geq b$, and apply the above argument to the equation
${\mathcal{A}}x^{m-1}=\hat{b}$. We can conclude that
${\mathcal{A}}x^{m-1}=\hat{b}$ has a nonnegative solution $\hat{x}_{*}$
(actually the unique positive solution) with $\hat{x}_{*}\geq x$ for every
element $x$ in $\\{x\in{\mathbb{R}}_{+}\ |\
{\mathcal{A}}x^{m-1}\leq\hat{b}\\}$ and thus for every element $x$ in
${\mathcal{S}}_{l}$.

We now return to the equation ${\mathcal{A}}x^{m-1}=b$ and take
$x_{0}\geq\hat{x}_{*}$ such that ${\mathcal{A}}x_{0}^{m-1}\geq b$. Then the
sequence $\\{x_{k}\\}$ from the Jacobi iteration converges to a nonnegative
solution $x_{*}$ of (1) and $x_{*}\geq x$ for all elements $x$ in
${\mathcal{S}}_{l}$. This $x_{*}$ is then the maximal element in
${\mathcal{S}}_{l}$ and also the maximal nonnegative solution of (1). ∎

######  Remark 1.

From the proof, we see that the maximal nonnegative solution (when exists) can
be found by iteration (3) using any $x_{0}$ such that

$x_{0}>0,\quad{\mathcal{A}}x_{0}^{m-1}>0,\quad{\mathcal{A}}x_{0}^{m-1}\geq b.$

In other words, we can take $x_{0}$ to be the unique positive solution of
${\mathcal{A}}x_{0}^{m-1}=\hat{b}$, where $\hat{b}$ is any vector such that
$\hat{b}>0$ and $\hat{b}\geq b$.

######  Remark 2.

If ${\mathcal{A}}$ is a $Z$-tensor, but not a nonsingular $M$-tensor, then
equation (1) may not have a maximal nonnegative solution when it has
nonnegative solutions. One example is the equation (1) with
${\mathcal{A}}\in{\mathbb{R}}^{[4,3]}$ given by:

$a_{1111}=0,a_{2222}=a_{3333}=1,a_{1112}=a_{3111}=-1,\mbox{ and
}a_{i_{1}i_{2}i_{3}i_{4}}=0\mbox{ elsewhere,}$

and $b=[0,0,1]^{T}$. The equation has infinitely many nonnegative solutions,
given by $[c,0,(1+c^{3})^{1/3}]^{T}$ for any $c\geq 0$. The minimal
nonnegative solution is $[0,0,1]^{T}$, but the maximal nonnegative solution
does not exist.

######  Corollary 6.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and $b\geq 0$. Then equation
(1) has a maximal nonnegative solution.

######  Remark 3.

When ${\mathcal{A}}$ is a nonsingular $M$-tensor and $b\geq 0$, we can use
Jacobi iteration with $x_{0}=0$ to get the minimal solution, and use $x_{0}>0$
with ${\mathcal{A}}x_{0}^{m-1}>0$ and ${\mathcal{A}}x_{0}^{m-1}\geq b$ to get
the maximal solution. When $b>0$, we can use either of these two choices to
get the unique positive solution. For the case $b>0$, the Jacobi iteration has
been studied in [4] with $x_{0}>0$ satisfying $0<{\mathcal{A}}x_{0}^{m-1}\leq
b$. More general tensor splitting methods have been studied in [13], again
with this requirement on $x_{0}$ (see [13, Theorem 5.4]). It is interesting to
note that both our choices are excluded by this requirement (unless $x_{0}$ is
already the solution).

######  Corollary 7.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor. Suppose that

${\mathcal{S}}_{l}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\leq
b\\}\neq\emptyset$

and let $x_{\max}({\mathcal{A}},b)$ be the maximal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$. If any element of $b$ increases, or if any entry of
${\mathcal{A}}$ decreases, then the new equation
$\tilde{{\mathcal{A}}}x^{m-1}=\tilde{b}$ also has a maximal nonnegative
solution $x_{\max}(\tilde{{\mathcal{A}}},\tilde{b})$ provided that
$\tilde{{\mathcal{A}}}$ is still a nonsingular $M$-tensor. Moreover,
$x_{\max}(\tilde{{\mathcal{A}}},\tilde{b})\geq x_{\max}({\mathcal{A}},b)$.

###### Proof.

Let $\tilde{{\mathcal{S}}}_{l}=\\{x\in{\mathbb{R}}_{+}\ |\
\tilde{{\mathcal{A}}}x^{m-1}\leq\tilde{b}\\}$. Then
${\mathcal{S}}_{g}\subseteq\tilde{{\mathcal{S}}}_{g}$ and the conclusions
follow immediately. ∎

The next result partially explain why the maximal nonnegative solution is of
particular interest for equation (1) with $b\geq 0$.

######  Corollary 8.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and $b\geq 0$. Suppose
$x^{(k)}$ is the unique positive solution of ${\mathcal{A}}x^{m-1}=b^{(k)}$,
where $b^{(k)}>0$ ($k=1,2,\ldots$) and $b^{(k)}$ is monotonically decreasing
and converges to $b$ as $k\to\infty$. Then, as $k\to\infty$, $x^{(k)}$
converges to the maximal nonnegative solution $x_{\max}({\mathcal{A}},b)$ of
${\mathcal{A}}x^{m-1}=b$.

###### Proof.

By Corollary 7, $x^{(k)}\geq x^{(k+1)}\geq x_{\max}({\mathcal{A}},b)$ for all
$k\geq 0$, so $\lim_{k\to\infty}x^{(k)}=x_{*}$ exists. Letting $k\to\infty$ in
${\mathcal{A}}(x^{(k)})^{m-1}=b^{(k)}$, we get ${\mathcal{A}}x_{*}^{m-1}=b$
and $x_{*}\geq x_{\max}({\mathcal{A}},b)$. Thus
$x_{*}=x_{\max}({\mathcal{A}},b)$. ∎

##  3\. Iterative methods for extremal nonnegative solutions

In the previous section, we have presented new proofs for the existence of
extremal nonnegative solutions of the tensor equation
${\mathcal{A}}x^{m-1}=b$, where ${\mathcal{A}}$ is a nonsingular $M$-tensor
and $b$ is a nonnegative vector, by using the Jacobi iteration with suitable
initial guesses. We have also presented some results when ${\mathcal{A}}$ is a
$Z$-tensor or $b$ is a general real vector. Now, we would like to present some
iterative methods that may be more efficient than the Jacobi iteration for
actual computation of the extremal solutions, determine and compare the rates
of convergence of these methods.

The Jacobi iteration is just a very simple fixed-point iteration. In the
Jacobi iteration, we keep the term involving $x_{i}^{m-1}$ in the $i$th
equation ($i\in[n]$) on the left and move all other terms to the right. But in
the $i$th equation, we may also have terms $x_{j}^{m-1}$ with $j\neq i$. We
may consider keeping all unmixed terms ($x_{1}^{m-1},\ldots,x_{n}^{m-1}$) on
the left, regardless which equation they are from. In other words, we have a
splitting of the tensor ${\mathcal{A}}$:
${\mathcal{A}}={\mathcal{M}}-{\mathcal{N}}$, where
${\mathcal{M}}=(m_{i_{1}i_{2}\ldots i_{m}})$ with $m_{ij\ldots j}=a_{ij\ldots
j}$ for $i,j\in[n]$ and $m_{i_{1}i_{2}\ldots i_{m}}=0$ elsewhere. We may call
this splitting a level-1 splitting. If we keep all unmixed terms on the left,
we then have the equation

${\mathcal{M}}x^{m-1}={\mathcal{N}}x^{m-1}+b.$

Note that ${\mathcal{M}}x^{m-1}=Mx^{[m-1]}$, where the $n\times n$ matrix $M$
has $(i,j)$ element $a_{ij...j}$ for $i,j\in[n]$, and is called [17] the
majorization matrix associated with ${\mathcal{A}}$. It is easy to see [13]
that $M$ is a nonsingular $M$-matrix when ${\mathcal{A}}$ is a nonsingular
$M$-tensor. When ${\mathcal{A}}$ is a $Z$-tensor, $M$ is obviously a
$Z$-matrix.

We then have a splitting of the matrix $M$: $M=P-Q$, where $P$ is a
nonsingular $M$-matrix and $Q\geq 0$. This may be called a level-2 splitting.
For example, the following splitting of a $Z$-matrix is permitted, although
not a good one.

$\left[\begin{array}[]{ccc}-1&0&0\\\ -3&2&-2\\\
0&-3&4\end{array}\right]=\left[\begin{array}[]{ccc}2&0&0\\\ -1&3&-1\\\
0&-2&5\end{array}\right]-\left[\begin{array}[]{ccc}3&0&0\\\ 2&1&1\\\
0&1&1\end{array}\right].$

We can now rewrite ${\mathcal{A}}x^{m-1}=b$ as

$Px^{[m-1]}=Qx^{[m-1]}+{\mathcal{N}}x^{m-1}+b$

and get fixed-point iteration in implicit form

(4) $Px_{k+1}^{[m-1]}=Qx_{k}^{[m-1]}+{\mathcal{N}}x_{k}^{m-1}+b$

or in explicit form

$x_{k+1}=\left(P^{-1}\left(Qx_{k}^{[m-1]}+{\mathcal{N}}x_{k}^{m-1}+b\right)\right)^{[1/(m-1)]},$

assuming $P^{-1}\left(Qx_{k}^{[m-1]}+{\mathcal{N}}x_{k}^{m-1}+b\right)\geq 0$
for each $k$. In [13], fixed-point iterations like this are called tensor
splitting iterative methods and studied for finding the unique positive
solution of the tensor equation ${\mathcal{A}}x^{m-1}=b$ with ${\mathcal{A}}$
being a nonsingular $M$-tensor and $b$ a positive vector.

We first study iteration (4) for finding the minimal nonnegative solution.

######  Theorem 9.

Let ${\mathcal{A}}$ be a $Z$-tensor and $b\in{\mathbb{R}}_{+}$. Suppose that

${\mathcal{S}}_{g}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\geq
b\\}\neq\emptyset.$

Then for $x_{0}=0$, the sequence $\\{x_{k}\\}$ from iteration (4) is
monotonically increasing and converges to the minimal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$.

###### Proof.

The proof is almost the same as the proof of Theorem 1. In that proof, we use
the obvious fact that ${\mathcal{D}}x^{m-1}\geq{\mathcal{D}}y^{m-1}$ (with
$x,y\geq 0$) implies $x\geq y$. We now need ${P}x^{[m-1]}\geq{P}y^{[m-1]}$
implies $x\geq y$, which is true since $P^{-1}\geq 0$ for the nonsingular
$M$-matrix $P$. ∎

######  Remark 4.

When ${\mathcal{A}}$ is a $Z$-tensor but not a nonsingular $M$-tensor, it may
not be easy to determine whether ${\mathcal{S}}_{g}\neq\emptyset$. We may
apply iteration (4) without checking this condition. In this case, the
sequence $\\{x_{k}\\}$ is still well-defined and monotonically increasing. The
sequence is bounded above if and only if ${\mathcal{S}}_{g}\neq\emptyset$.

In iteration (4), ${\mathcal{N}}$ is uniquely determined by ${\mathcal{A}}$,
but we have the freedom to choose $Q$ ($P$ is uniquely determined by $Q$). The
next result gives a comparison of convergence rate by examining the iterates
right from the beginning (so we are not talking about asymptotic rate of
convergence here).

######  Theorem 10.

Under the conditions of Theorem 9, let $\\{x_{k}\\}$ and $\\{\hat{x}_{k}\\}$
be the sequences generated by iteration (4) with two splittings of $M$:
$M=P-Q$ and $M=\hat{P}-\hat{Q}$, respectively, and with $x_{0}=\hat{x}_{0}=0$.
If $\hat{Q}\leq Q$, then $x_{k}\leq\hat{x}_{k}$ for all $k\geq 0$. In other
words, a smaller matrix $Q$ gives faster termwise convergence for finding the
minimal nonnegative solution.

###### Proof.

We need to show that $x_{k}\leq\hat{x}_{k}$ implies $x_{k+1}\leq\hat{x}_{k+1}$
for each $k\geq 0$. We have (4) and

(5)
$\hat{P}\hat{x}_{k+1}^{[m-1]}=\hat{Q}\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b.$

From $P-Q=\hat{P}-\hat{Q}$, we get $\hat{P}={P}-(Q-\hat{Q})$ and get from (5)
that

$\displaystyle P\hat{x}_{k+1}^{[m-1]}$ $\displaystyle=$
$\displaystyle(Q-\hat{Q})\hat{x}_{k+1}^{[m-1]}+\hat{Q}\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$
$\displaystyle\geq$
$\displaystyle(Q-\hat{Q})\hat{x}_{k}^{[m-1]}+\hat{Q}\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$
$\displaystyle=$ $\displaystyle
Q\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$ $\displaystyle\geq$
$\displaystyle Qx_{k}^{[m-1]}+{\mathcal{N}}{x}_{k}^{m-1}+b$ $\displaystyle=$
$\displaystyle Px_{k+1}^{[m-1]}.$

Thus $x_{k+1}\leq\hat{x}_{k+1}$. ∎

We now show that iteration (4) has linear convergence under suitable
assumptions. By linear convergence we actually mean at least linear
convergence. For example, $x_{0}=0$ is already a solution when $b=0$. Let
${\mathcal{L}}$ be the nonnegative tensor such that
${\mathcal{L}}x^{m-1}=Qx^{[m-1]}+{\mathcal{N}}x^{m-1}$. Then iteration (4)
becomes

(6) $Px_{k+1}^{[m-1]}={\mathcal{L}}x_{k}^{m-1}+b.$

For $x\in{\mathbb{R}}^{n}$ and index set $I\subseteq[n]$, we denote by $x_{I}$
the subvector of $x$, whose elements are $x_{i}$, $i\in I$. For tensor
${\mathcal{A}}=(a_{i_{1}\ldots i_{m}})\in{\mathbb{R}}^{[m,n]}$, we denote by
${\mathcal{A}}_{I}$ the subtensor of ${\mathcal{A}}$ with elements
$a_{i_{1}\ldots i_{m}}$, $i_{1},\ldots,i_{m}\in I$.

######  Theorem 11.

Let ${\mathcal{A}}$ be a $Z$-tensor and $b\in{\mathbb{R}}_{+}$. Let
$I_{0}=\\{i\ |\ b_{i}=0\\}$ and suppose that
${\mathcal{S}}_{g}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\geq
b\\}\neq\emptyset$. Let $\\{x_{k}\\}$ be the sequence from iteration (6) with
$x_{0}=0$.

  1. (1)

If $I_{0}=[n]$, then $x_{0}=0$ is already the minimal nonnegative solution.

  2. (2)

If $I_{0}=\emptyset$, then $\\{x_{k}\\}$ converges linearly to the minimal
nonnegative solution (which is actually the unique positive solution) of
${\mathcal{A}}x^{m-1}=b$.

  3. (3)

If $I_{0}$ is a proper subset of $[n]$ and $k_{0}$ is the smallest integer
such that $x_{k_{0}}$ and $x_{k_{0}+1}$ have the same zero pattern, then
$1\leq k_{0}\leq n$ and $x_{k}$ have the same zero pattern for all $k\geq
k_{0}$, which is also the zero pattern of the minimal nonnegative solution
$x_{\min}$. Let $I=\\{i\ |\ (x_{k_{0}})_{i}=0\\}$. Then $I\subseteq I_{0}$.
Let $I_{c}=[n]\setminus I$. Then the iteration (6) for $k\geq k_{0}$ is
reduced to an iteration for the lower-dimensional tensor equation:

(7) ${\mathcal{A}}_{I_{c}}\hat{x}^{m-1}=b_{I_{c}}.$

For $k\geq k_{0}$, $\hat{x}_{k}$ from the reduced iteration is the same as
$(x_{k})_{I_{c}}$. The minimal nonnegative solution $\hat{x}_{\min}$ of (7) is
positive and is the same as $(x_{\min})_{I_{c}}$. Thus $x_{k}$ converges to
$x_{\min}$ linearly if and only if $\hat{x}_{k}$ converges to $\hat{x}_{\min}$
linearly. Assume ${\mathcal{A}}{x}^{m-1}=b$ is already the reduced equation
for notation convenience (all diagonal elements of ${\mathcal{A}}$ are now
positive). Then $x_{k}$ converges to $x_{\min}$ linearly under any of the
following four conditions:

    1. (a)

$P^{-1}b>0$.

    2. (b)

$P^{-1}\bar{\mathcal{L}}e^{m-2}$ is irreducible, where $e$ is the vector of
ones and $\bar{\mathcal{L}}$ is the semi-symmetric tensor from
${\mathcal{L}}$.

    3. (c)

With all diagonal elements of ${\mathcal{A}}$ included entirely in $P$ and
${\mathcal{B}}={\mathcal{D}}-{\mathcal{A}}$ (${\mathcal{D}}$ is the diagonal
part of ${\mathcal{A}}$), the matrix $\bar{\mathcal{B}}e^{m-2}$ is strictly
upper or lower triangular.

    4. (d)

For each $i\in I_{0}$, there is an element $a_{ii_{2}\ldots i_{m}}\neq 0$,
where $i_{j}\notin I_{0}$ for at least one $j$ ($2\leq j\leq m$).

###### Proof.

Conclusion 1 in the theorem is obvious. The proof of conclusion 2 is the same
as the proof of conclusion 3 for the reduced equation under condition (a). We
now prove conclusion 3.

With $x_{0}=0$ for iteration (6), we have $x_{k}\leq x_{k+1}$ for each $k\geq
0$. Let $k_{0}$ be the smallest integer such that $x_{k_{0}}$ and
$x_{k_{0}+1}$ have the same zero pattern. It is clear that $1\leq k_{0}\leq n$
and that the zero pattern will not change afterwards and the minimal solution
$x_{\min}$ has the same zero pattern. Let $I=\\{i\ |\ (x_{k_{0}})_{i}=0\\}$.
Since $b_{i}>0$ implies $(x_{k_{0}})_{i}\geq(x_{1})_{i}>0$ for a diagonal $P$,
we have $I\subseteq I_{0}$. For $k\geq k_{0}$, we have $(x_{k})_{I}=0$. Let
$I_{c}=[n]\setminus I$. Since
$(x_{k})_{i_{2}}(x_{k})_{i_{3}}\cdots(x_{k})_{i_{m}}=0$ for all $k\geq k_{0}$
when $i_{j}\in I$ for at least one $j$ ($j=2,\ldots,m$), we only need the
elements in ${\mathcal{A}}_{I_{c}}$ to continue the iteration (6) for $k\geq
k_{0}$. In other words, the iteration (6) for $k\geq k_{0}$ is reduced to an
iteration for the lower-dimensional tensor equation with indices running
through $I_{c}$ only:

(8) ${\mathcal{A}}_{I_{c}}\hat{x}^{m-1}=b_{I_{c}},$

with $\hat{x}_{k_{0}}$ obtained from $x_{k_{0}}$ by deleting all zero
elements. Note that the level-1 and level-2 splittings of
${\mathcal{A}}_{I_{c}}$ are the ones inherited from those for ${\mathcal{A}}$.
The minimal nonnegative solution $\hat{x}_{\min}$ of (8) is positive (obtained
from $x_{\min}$ by deleting all zero elements), but $b_{I_{c}}$ will have some
zero elements when $I\neq I_{0}$. For $k\geq k_{0}$, $\hat{x}_{k}$ from the
reduced iteration can also be obtained from $x_{k}$ by deleting all zero
elements. Thus $x_{k}$ converges to $x_{\min}$ linearly if and only if
$\hat{x}_{k}$ converges to $\hat{x}_{\min}$ linearly. We now assume that
${\mathcal{A}}{x}^{m-1}=b$ is already the reduced equation for notation
convenience, and denote its positive minimal solution by $\bar{x}$ for short.
Note that all diagonal elements of the reduced tensor are positive.

By Theorem 10, we may assume $P$ is a diagonal matrix (where slower
convergence happens). So we will assume $P$ is diagonal when needed.

The iteration (6) can be written explicitly as $x_{k+1}=\phi(x_{k})$, with

$\phi(x)=\left(P^{-1}({\mathcal{L}}x^{m-1}+b)\right)^{[1/(m-1)]}.$

From

$P\phi(x)^{[m-1]}={\mathcal{L}}x^{m-1}+b=\bar{\mathcal{L}}x^{m-1}+b,$

where $\bar{\mathcal{L}}$ is the semi-symmetric tensor obtained from
${\mathcal{L}}$, we take derivative on both sides to obtain

$P(m-1){\rm
diag}(\phi(x)^{[m-2]})\phi^{\prime}(x)=(m-1)\bar{\mathcal{L}}x^{m-2},$

where for $x\in{\mathbb{R}}^{n}$, ${\rm diag}(x)$ is the diagonal matrix with
the elements of $x$ on the diagonal. We need to show
$\rho(\phi^{\prime}(\bar{x}))<1$ for

(9) $\phi^{\prime}(\bar{x})={\rm
diag}(\bar{x}^{[-(m-2)]})P^{-1}\bar{\mathcal{L}}{\bar{x}}^{m-2}.$

Note that the nonnegative matrix $\phi^{\prime}(\bar{x})$ is such that

$\displaystyle\phi^{\prime}(\bar{x})\bar{x}$ $\displaystyle=$
$\displaystyle{\rm
diag}(\bar{x}^{[-(m-2)]})P^{-1}\bar{\mathcal{L}}{\bar{x}}^{m-2}\bar{x}$
$\displaystyle=$ $\displaystyle{\rm
diag}(\bar{x}^{[-(m-2)]})P^{-1}\bar{\mathcal{L}}{\bar{x}}^{m-1}$
$\displaystyle=$ $\displaystyle{\rm
diag}(\bar{x}^{[-(m-2)]})P^{-1}{\mathcal{L}}{\bar{x}}^{m-1}$ $\displaystyle=$
$\displaystyle{\rm
diag}(\bar{x}^{[-(m-2)]})P^{-1}(P{\bar{x}}^{[m-1]}-{\mathcal{A}}{\bar{x}}^{m-1})$
$\displaystyle=$ $\displaystyle\bar{x}-{\rm diag}(\bar{x}^{[-(m-2)]})P^{-1}b.$

If $P^{-1}b>0$, then we have $\phi^{\prime}(\bar{x})\bar{x}<\bar{x}$ and thus
[2] $\rho(\phi^{\prime}(\bar{x}))<1$. This shows that condition (a) is
sufficient for linear convergence. Conditions (b), (c) and (d) are needed only
when $I_{0}\neq\emptyset$ for the reduced equation.

Since $P^{-1}b\geq 0$ and $P^{-1}b\neq 0$, we have
$\phi^{\prime}(\bar{x})\bar{x}\leq\bar{x}$ and
$\phi^{\prime}(\bar{x})\bar{x}\neq\bar{x}$. Thus we still have [2]
$\rho(\phi^{\prime}(\bar{x}))<1$ if $\phi^{\prime}(\bar{x})$ is irreducible.
From (9), we see that $\phi^{\prime}(\bar{x})$ is irreducible if and only if
$P^{-1}\bar{\mathcal{L}}e^{m-2}$ is irreducible. Thus, condition (b) is also
sufficient for linear convergence.

Under condition (c), we just need to show linear convergence for the Jacobi
iteration. Now,

$\phi^{\prime}(\bar{x})={\rm
diag}(\bar{x}^{[-(m-2)]})D^{-1}\bar{\mathcal{B}}{\bar{x}}^{m-2}.$

Since $\bar{\mathcal{B}}e^{m-2}$ is strictly upper or lower triangular, so is
$\phi^{\prime}(\bar{x})$. Thus $\rho(\phi^{\prime}(\bar{x}))=0$ and linear
convergence follows.

Finally, we show that condition (d) is also sufficient for linear convergence.
We now assume $P$ is diagonal. We modify $\bar{x}$ to $\hat{x}$ by changing
$(\bar{x})_{i}$ to $(\bar{x})_{i}-\epsilon>0$ for each $i\notin I_{0}$, where
$0<\epsilon<\min_{i\notin I_{0}}(\bar{x})_{i}^{-(m-2)}(P^{-1}b)_{i}$. Now for
$i\notin I_{0}$,

$(\phi^{\prime}(\bar{x})\hat{x})_{i}\leq(\phi^{\prime}(\bar{x})\bar{x})_{i}=(\bar{x})_{i}-(\bar{x})_{i}^{-(m-2)}(P^{-1}b)_{i}<(\hat{x})_{i}.$

For $i\in I_{0}$,

$(\phi^{\prime}(\bar{x})\hat{x})_{i}\leq(\phi^{\prime}(\bar{x})\bar{x})_{i}=(\bar{x})_{i}=(\hat{x})_{i}.$

Since $P$ is diagonal,
$(\phi^{\prime}(\bar{x})\hat{x})_{i}=(\phi^{\prime}(\bar{x})\bar{x})_{i}$ if
and only if
$(\bar{\mathcal{L}}{\bar{x}}^{m-2}\hat{x})_{i}=(\bar{\mathcal{L}}{\bar{x}}^{m-2}\bar{x})_{i}$,
i.e.,

$\displaystyle\sum_{j=1}^{n}\left(\sum_{i_{3},\ldots,i_{m}=1}^{n}\bar{\mathcal{L}}_{iji_{3}\ldots
i_{m}}\bar{x}_{i_{3}}\cdots\bar{x}_{i_{m}}\right)\hat{x}_{j}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}\left(\sum_{i_{3},\ldots,i_{m}=1}^{n}\bar{\mathcal{L}}_{iji_{3}\ldots
i_{m}}\bar{x}_{i_{3}}\cdots\bar{x}_{i_{m}}\right)\bar{x}_{j}.$

This holds if and only if $\bar{\mathcal{L}}_{iji_{3}\ldots i_{m}}=0$ for all
$j\notin I_{0}$ and for all $i_{3},\ldots,i_{m}\in[n]$, i.e.,
${\mathcal{L}}_{ii_{2}i_{3}\ldots i_{m}}=0$ for all $i_{2},\ldots,i_{m}\in[n]$
such that $i_{j}\notin I_{0}$ for at least one $j$ ($2\leq j\leq m$).
Therefore, when condition (d) holds,
$(\phi^{\prime}(\bar{x})\hat{x})_{i}<(\phi^{\prime}(\bar{x})\bar{x})_{i}=(\hat{x})_{i}$
for each $i\in I_{0}$. Linear convergence follows since we again have
$\rho(\phi^{\prime}(\bar{x}))<1$. ∎

We now study iteration (4) for finding the maximal nonnegative solution.

######  Theorem 12.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and suppose that
${\mathcal{S}}_{l}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\leq
b\\}\neq\emptyset$. Then for any $x_{0}>0$ such that
${\mathcal{A}}x_{0}^{m-1}>0$ and ${\mathcal{A}}x_{0}^{m-1}\geq b$. The
sequence from iteration (4) is monotonically decreasing and converges to the
maximal nonnegative solution of ${\mathcal{A}}x^{m-1}=b$.

###### Proof.

The proof is almost the same as the proof of Theorem 5. In that proof, we use
the fact that ${\mathcal{D}}x^{m-1}\geq{\mathcal{D}}y^{m-1}$ (with $x,y\geq
0$) implies $x\geq y$. We now use the fact ${P}x^{[m-1]}\geq{P}y^{[m-1]}$
implies $x\geq y$. ∎

######  Remark 5.

When $b\notin{\mathbb{R}}^{+}$, it may not be easy to determine whether
${\mathcal{S}}_{l}\neq\emptyset$. We may apply iteration (4) without checking
this condition. In this case, ${\mathcal{S}}_{l}\neq\emptyset$ if and only if
$P^{-1}\left(Qx_{k}^{[m-1]}+{\mathcal{N}}x_{k}^{m-1}+b\right)\geq 0$ for each
$k\geq 0$. Indeed, if
$P^{-1}\left(Qx_{k}^{[m-1]}+{\mathcal{N}}x_{k}^{m-1}+b\right)\geq 0$ for each
$k\geq 0$, then we see from a proof similar to that of Theorem 5 that
$x_{k}\geq x_{k+1}\geq 0$ for all $k\geq 0$. In this case,
$\lim_{k\to\infty}x_{k}=x_{*}$ exists and $x_{*}$ is a solution of
${\mathcal{A}}x^{m-1}=b$ and thus ${\mathcal{S}}_{l}\neq\emptyset$.

######  Remark 6.

To obtain a suitable $x_{0}$ in Theorem 12, we can always take a vector
$\hat{b}>0$ with $\hat{b}\geq b$, and use the methods in [4, 7, 8] to get the
unique positive solution $x_{*}$ of the equation
${\mathcal{A}}x^{m-1}=\hat{b}$ and then take $x_{0}=x_{*}$. When we need to
solve the tensor equation with the same ${\mathcal{A}}$ and many different
right side vectors $b_{1},\ldots,b_{p}$, we can take a vector $\hat{b}>0$ with
$\hat{b}\geq\max\\{b_{1},\ldots,b_{p}\\}$ (where the maximum is taken
elementwise) and get the unique positive solution $x_{*}$ of the equation
${\mathcal{A}}x^{m-1}=\hat{b}$ and then take $x_{0}=x_{*}$ for all equations
${\mathcal{A}}x^{m-1}={b}_{i}$, $i=1,\ldots,p$.

######  Theorem 13.

Under the conditions of Theorem 12, let $\\{x_{k}\\}$ and $\\{\hat{x}_{k}\\}$
be the sequences generated by iteration (4) with two splittings of $M$:
$M=P-Q$ and $M=\hat{P}-\hat{Q}$, respectively, and with
$x_{0}\geq\hat{x}_{0}>0$ satisfying

${\mathcal{A}}x_{0}^{m-1}>0,\quad{\mathcal{A}}x_{0}^{m-1}\geq
b,\quad{\mathcal{A}}\hat{x}_{0}^{m-1}>0,\quad{\mathcal{A}}\hat{x}_{0}^{m-1}\geq
b.$

If $\hat{Q}\leq Q$, then $x_{k}\geq\hat{x}_{k}$ for all $k\geq 0$. In other
words, a smaller matrix $Q$ gives faster termwise convergence for finding the
maximal solution.

###### Proof.

We need to show that $x_{k}\geq\hat{x}_{k}$ implies $x_{k+1}\geq\hat{x}_{k+1}$
for each $k\geq 0$. We have (4) and (5). Since $\hat{P}={P}-(Q-\hat{Q})$, we
get from (5) that

$\displaystyle P\hat{x}_{k+1}^{[m-1]}$ $\displaystyle=$
$\displaystyle(Q-\hat{Q})\hat{x}_{k+1}^{[m-1]}+\hat{Q}\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$
$\displaystyle\leq$
$\displaystyle(Q-\hat{Q})\hat{x}_{k}^{[m-1]}+\hat{Q}\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$
$\displaystyle=$ $\displaystyle
Q\hat{x}_{k}^{[m-1]}+{\mathcal{N}}\hat{x}_{k}^{m-1}+b$ $\displaystyle\leq$
$\displaystyle Qx_{k}^{[m-1]}+{\mathcal{N}}{x}_{k}^{m-1}+b$ $\displaystyle=$
$\displaystyle Px_{k+1}^{[m-1]}.$

Thus $x_{k+1}\geq\hat{x}_{k+1}$. ∎

We now study the convergence rate of iteration (4) for finding the maximal
nonnegative solution of (1), under the conditions of Theorem 12. In our first
result we assume that the maximal nonnegative solution is positive. We know
from [10, Theorem 2.4] that every nonnegative solution of (1) is positive if
${\mathcal{A}}$ is irreducible and $b\geq 0$ is nonzero. We also know from
Corollary 7 that if the maximal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$ is positive, then the maximal nonnegative solution of
$\tilde{{\mathcal{A}}}x^{m-1}=\tilde{b}$ will also be positive if
$\tilde{{\mathcal{A}}}$ and $\tilde{b}$ are obtained from ${\mathcal{A}}$ and
$b$ in ways described there.

######  Theorem 14.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and suppose that
${\mathcal{S}}_{l}=\\{x\in{\mathbb{R}}_{+}\ |\ {\mathcal{A}}x^{m-1}\leq
b\\}\neq\emptyset$. Let $x_{\max}$ be the maximal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$ and assume that $x_{\max}>0$. Let $I_{0}=\\{i\ |\
b_{i}=0\\}$. Then for any $x_{0}>0$ such that ${\mathcal{A}}x_{0}^{m-1}>0$ and
${\mathcal{A}}x_{0}^{m-1}\geq b$, the sequence $\\{x_{k}\\}$ from iteration
(4) converges to $x_{max}$ linearly under any of the following four
conditions:

  1. (1)

$P^{-1}b>0$.

  2. (2)

$P^{-1}\bar{\mathcal{L}}e^{m-2}$ is irreducible, where $e$ is the vector of
ones and $\bar{\mathcal{L}}$ is the semi-symmetric tensor from
${\mathcal{L}}$.

  3. (3)

With all diagonal elements of ${\mathcal{A}}$ included entirely in $P$ and
${\mathcal{B}}={\mathcal{D}}-{\mathcal{A}}$, the matrix
$\bar{\mathcal{B}}e^{m-2}$ is strictly upper or lower triangular.

  4. (4)

For each $i\in I_{0}$, there is an element $a_{ii_{2}\ldots i_{m}}\neq 0$,
where $i_{j}\notin I_{0}$ for at least one $j$ ($2\leq j\leq m$).

###### Proof.

When $b=0$, ${\mathcal{A}}x^{m-1}=b$ has a unique solution $x=0$. Thus $b\neq
0$ when $x_{\max}>0$. Exactly as in the proof of Theorem 11 (we do not need
the reduction process there since we assume $x_{\max}>0$), we can show that
the iteration map $\phi$ is such that $\rho(\phi^{\prime}(x_{\max}))<1$ under
any of the four conditions in the theorem. In proving linear convergence under
condition (c) or (d), we may assume that $P$ is diagonal. ∎

We now assume $b\geq 0$, but allow $x_{\max}$ to have some zero elements.

######  Theorem 15.

Let ${\mathcal{A}}$ be a nonsingular $M$-tensor and $b\in{\mathbb{R}}_{+}$ be
nonzero. Let $x_{\max}$ be the maximal nonnegative solution of
${\mathcal{A}}x^{m-1}=b$. Let $I_{0}=\\{i\ |\ b_{i}=0\\}$. Then for iteration
(4) with any $x_{0}>0$ such that ${\mathcal{A}}x_{0}^{m-1}>0$ and
${\mathcal{A}}x_{0}^{m-1}\geq b$, there is a smallest integer $k_{0}$ ($0\leq
k_{0}\leq n-1$) such that $x_{k_{0}}$ and $x_{k_{0}+1}$ have the same zero
pattern (including the case with no zero elements), and $x_{k}$ have the same
zero pattern for all $k\geq k_{0}$. Let $I=\\{i\ |\ (x_{k_{0}})_{i}=0\\}$
(which may be empty). Then $I\subseteq I_{0}$. Let $I_{c}=[n]\setminus I$.
Then the iteration (4) for $k\geq k_{0}$ is reduced to an iteration for the
lower-dimensional tensor equation:

(10) ${\mathcal{A}}_{I_{c}}\hat{x}^{m-1}=b_{I_{c}}.$

For $k\geq k_{0}$, $\hat{x}_{k}$ from the reduced iteration is the same as
$(x_{k})_{I_{c}}$. The maximal nonnegative solution $\hat{x}_{\max}$ of (10)
is the same as $(x_{\max})_{I_{c}}$ (which is not necessarily positive). Thus
$x_{k}$ converges to $x_{\max}$ linearly if and only if $\hat{x}_{k}$
converges to $\hat{x}_{\max}$ linearly. Assume ${\mathcal{A}}{x}^{m-1}=b$ is
already the reduced equation and assume that its maximal solution $x_{\max}$
is positive. Then $x_{k}$ converges to $x_{\max}$ linearly under any of the
following four conditions:

  1. (1)

$P^{-1}b>0$.

  2. (2)

$P^{-1}\bar{\mathcal{L}}e^{m-2}$ is irreducible, where $e$ is the vector of
ones and $\bar{\mathcal{L}}$ is the semi-symmetric tensor from
${\mathcal{L}}$.

  3. (3)

With all diagonal elements of ${\mathcal{A}}$ included entirely in $P$ and
${\mathcal{B}}={\mathcal{D}}-{\mathcal{A}}$, the matrix
$\bar{\mathcal{B}}e^{m-2}$ is strictly upper or lower triangular.

  4. (4)

For each $i\in I_{0}$, there is an element $a_{ii_{2}\ldots i_{m}}\neq 0$,
where $i_{j}\notin I_{0}$ for at least one $j$ ($2\leq j\leq m$).

###### Proof.

We have $x_{0}>0$ and $x_{k}\geq x_{k+1}$ for all $k\geq 0$. Each $x_{k}$ has
at least one nonzero entry since $x_{\max}\neq 0$. Thus, there is a smallest
integer $k_{0}$ ($0\leq k_{0}\leq n-1$) such that $x_{k_{0}}$ and
$x_{k_{0}+1}$ have the same zero pattern. Since $b\geq 0$, $x_{k}$ have the
same zero pattern for all $k\geq k_{0}$. It is clear that $I\subseteq I_{0}$.
The reduced equation (10) is then obtained as in the proof of Theorem 11.
However, its maximal solution $x_{\max}$ (as the limit of a positive sequence)
may have some zero elements. With the additional assumption that $x_{\max}>0$
for the reduced equation, the proof of linear convergence is exactly the same
as in the proof of Theorem 11. ∎

######  Remark 7.

The reduction in Theorem 11 is obtained from the iteration for finding the
minimal solution and the reduction in Theorem 15 is obtained from the
iteration for finding the maximal solution. A reduction process has also been
described in [10], without mentioning maximal and minimal solutions. We will
explain that one can only find a nonnegative solution with the same number of
nonzero elements as the minimal solution by finding a positive solution of the
reduced equation in [10].

To describe the approach in [10], we recall the following definition.

######  Definition 2.

A tensor ${\mathcal{A}}\in{\mathbb{R}}^{[m,n]}$ is called reducible with
respect to $I\subset[n]$ if its elements satisfy

$a_{i_{1}i_{2}\dots i_{m}}=0,\forall i_{1}\in I,\forall
i_{2},\ldots,i_{m}\notin I.$

The next result has been presented in [10] (see Corollary 2.8 there).

######  Theorem 16.

Suppose that ${\mathcal{A}}$ is a nonsingular $M$-tensor and
$b\in{\mathbb{R}}^{n}_{+}$ is nonzero. Then there is an index set $I\subseteq
I_{0}$ (which could be empty) such that every nonnegative solution to the
following lower dimensional tensor equation with $I_{c}=[n]\setminus I$

${\mathcal{A}}_{I_{c}}x^{m-1}_{I_{c}}=b_{I_{c}}$

is positive. Moreover, every positive solution $x_{I_{c}}$ of the last
equation together with $x_{I}=0$ forms a nonnegative solution to equation (1).

By comparing the discussions in [10] and in this paper, we can see that the
index set $I\subseteq I_{0}$ in Theorem 16 is the largest set such that
${\mathcal{A}}$ is reducible with respect to $I$. We also see that the index
set $I$ can be determined automatically by iteration (4) for computing the
minimal nonnegative solution, with $x_{0}=0$. Indeed, $I=\\{i\ |\
(x_{k})_{i}=0\\}$, where $k$ is the smallest integer such that the vectors
$x_{k}$ and $x_{k+1}$ from iteration (4) have the same zero pattern (see
Theorem 11; $I=\emptyset$ when $b>0$). This is a very easy way to determine
the set $I$. It would be much more expensive to determine $I$ by using the
definition of reducibility to find the largest index set $I\subseteq I_{0}$
such that ${\mathcal{A}}$ is reducible with respect to $I$.

The numerical methods in [10] and [11] for computing a nonnegative solution of
(1) are based on Theorem 16. They have quadratic convergence under suitable
assumptions. For example, assuming equation (1) has already been reduced, the
assumption needed for a Newton method in [11] is that for each $i\in I_{0}$,
there is an element $a_{ii_{2}\ldots i_{m}}\neq 0$ with all $i_{j}\notin
I_{0}$ ($j=2,\ldots,m$). So the assumption is stronger than our condition 3(d)
in Theorem 11 for the linear convergence of our simple iteration for finding
the minimal solution. From our comments on Theorem 16, we know that the
methods in [10] and [11] can only find one of the nonnegative solutions that
has the same number of zero elements as the minimal solution. In particular,
they will never find the maximal solution if it have more nonzero elements
than the minimal solution. Moreover, to use the methods in [10] and [11], the
tensor ${\mathcal{A}}$ should be semi-symmetrized first, which increases
computational work and makes the tensor much less sparse.

######  Example 2.

We consider Example 1 with $k=1$. Then equation (1) has two solutions:
$[0,1]^{T}$ and $[2,1]^{T}$. Note that $I_{0}=\\{1\\}$. We have $I=\\{1\\}$
for Theorem 16 since ${\mathcal{A}}$ is reducible with respect to $\\{1\\}$.
The reduced equation is $x_{2}^{3}=1$ with nonnegative solution $x_{2}=1$,
which is positive. A nonnegative solution of the original equation is then
$[0,1]^{T}$ by Theorem 16, but the other solution is lost in the reduction. We
now apply Theorem 11. With $x_{0}=[0,0]^{T}$, we get $x_{1}=[0,1]^{T}$ and
$x_{2}=[0,1]^{T}$. So $k_{0}=1$ and $I=\\{1\\}$ in Theorem 11. The reduced
equation is $x_{2}^{3}=1$ with minimal nonnegative solution $x_{2}=1$, which
is positive, and the minimal nonnegative solution of the original equation is
then $[0,1]^{T}$ by Theorem 11 (we got this solution after just one
iteration). We then apply Theorem 15. With $x_{0}=[3,1]^{T}$, we get
$x_{1}=[18^{1/3},1]^{T}$. So $k_{0}=0$ and $I=\\{1,2\\}$ in Theorem 15. The
equation is thus not reduced. We have $x_{k}=[s_{k},1]^{T}$, where $s_{k}$ is
determined by the iteration: $t_{0}=3$, $t_{k+1}=(2t_{k}^{2})^{1/3}$. So
$t_{k}$ converges to $2$ linearly with rate $2/3$, and $x_{k}$ converges to
$[2,1]^{T}$ at the same rate. The linear convergence of $\\{x_{k}\\}$ is also
guaranteed by Theorem 14 since condition 4 there is satisfied ($a_{1112}\neq
0$).

There are examples for which iteration (4) converges linearly with a rate very
close to $1$. We consider another extension of Example 1.1 in [1].

######  Example 3.

We consider equation (1) with ${\mathcal{A}}\in{\mathbb{R}}^{[m,2]}$ given by
$a_{11\ldots 1}=a_{22\ldots 2}=1$ and $a_{11\ldots 12}=-2$ and with
$b=[0,1]^{T}$. Then equation (1) has two solutions: $[0,1]^{T}$ and
$[2,1]^{T}$. Note that $I_{0}=\\{1\\}$. With $x_{0}=[3,1]^{T}$, we get
$x_{1}=[(2\cdot 3^{m-2})^{1/(m-1)},1]^{T}$. So $k_{0}=0$ and $I=\\{1,2\\}$ in
Theorem 15. The equation is thus not reduced. We have $x_{k}=[s_{k},1]^{T}$,
where $s_{k}$ is determined by the iteration: $t_{0}=3$,
$t_{k+1}=(2t_{k}^{m-2})^{1/(m-1)}$. So $t_{k}$ converges to $2$ linearly with
rate $(m-2)/(m-1)$, and $x_{k}$ converges to $[2,1]^{T}$ at the same rate. The
linear convergence of $\\{x_{k}\\}$ is also guaranteed by Theorem 14 since
$a_{11\ldots 12}\neq 0$. However, when $m$ is large, the rate is very close to
$1$, so $x_{k}$ converges to $[2,1]^{T}$ very slowly.

######  Remark 8.

We have the freedom to choose the splitting $M=P-Q$ for iteration (4), when
${\mathcal{A}}$ is a nonsingular $M$-tensor. Recall that a smaller $Q$ is
going to give faster termwise convergence (see Theorems 10 and 13). For a
dense tensor with order $m\geq 4$, computing ${\mathcal{N}}x^{m-1}$ will
require $O(n^{m})$ flops (for large $n$ and fixed $m$) and solving the linear
system $My=c$ requires $O(n^{3})$ flops. So it is advisable to use $P=M$ and
$Q=0$ in the splitting $M=P-Q$ for a dense tensor with order $m\geq 4$. For a
dense tensor with order $3$, it is advisable to take $P$ to be the lower
triangular part of $M$ or the upper triangular part of $M$, since solving the
linear system $Py=c$ requires $O(n^{2})$ flops in this case.

##  4\. Conclusion

We have presented new proofs for the existence of extremal nonnegative
solutions of the $M$-tensor equation with a nonnegative right side vector by
using some simple fixed-point iterations. We have studied these iterative
methods and their generalizations. With a suitable starting point, each of
these methods has monotonic convergence (to the maximal nonnegative solution
or to the minimal nonnegative solution) and the rate of convergence is (at
least) linear under some mild assumptions. These methods are currently the
only methods that are guaranteed to compute the maximal nonnegative solution
or the minimal nonnegative solution with suitable initial guesses. There are
examples for which the presented iterations have linear convergence at a rate
very close to $1$. Iterative methods with faster convergence (without
increasing computational work each iteration by too much) are still desirable.
Those methods should be ones that can compute the maximal nonnegative solution
and/or the minimal nonnegative solution, not just any one of the nonnegative
solutions.

## Acknowledgments

The author was supported in part by an NSERC Discovery Grant.

## References

  * [1] Bai, X., He, H., Ling, C., Zhou, G.: A nonnegativity preserving algorithm for multilinear systems with nonsingular $M$-tensors. Numer. Algorithms. 87, 1301–1320 (2021)
  * [2] Berman, A., Plemmons, R.: Nonnegative Matrices in the Mathematical Sciences. SIAM, Philadephia (1994)
  * [3] Ding, W., Qi, L., Wei, Y.: $M$-tensors and nonsingular $M$-tensors. Linear Algebra Appl. 439, 3264–3278 (2013)
  * [4] Ding, W., Wei, Y.: Solving multi-linear systems with $M$-tensors. J. Sci. Comput. 68, 683–715 (2016)
  * [5] Guo, C.-H.: Nonsymmetric algebraic Riccati equations and Wiener–Hopf factorization for $M$-matrices. SIAM J. Matrix Anal. Appl. 23, 225–242 (2001)
  * [6] Guo, C.-H., Laub, A.J.: On the iterative solution of a class of nonsymmetric algebraic Riccati equations. SIAM J. Matrix Anal. Appl. 22, 376–391 (2000)
  * [7] Han, L.: A homotopy method for solving multilinear systems with $M$-tensors. Appl. Math. Lett. 69, 49–54 (2017)
  * [8] He, H., Ling, C., Qi, L., Zhou, G.: A globally and quadratically convergent algorithm for solving multilinear systems with $M$-tensors. J. Sci. Comput. 76, 1718–1741 (2018)
  * [9] Li, D.H., Guan, H.B., Wang, X.Z.: Finding a nonnegative solution to an $M$-tensor equation. Pac. J. Optim. 16, 419–440 (2020)
  * [10] Li, D.H., Guan, H.B., Xu, J.F.: Inexact Newton method for $M$-tensor equations. Pac. J. Optim. 17, 617–643 (2021)
  * [11] Li, D.H., Xu, J.F., Guan, H.B.: Newton’s method for $M$-tensor equations. J. Optim. Theory Appl. 190, 628–649 (2021)
  * [12] Liu, C.-S., Guo, C.-H., Lin, W.-W.: Newton–Noda iteration for finding the Perron pair of a weakly irreducible nonnegative tensor. Numer. Math. 137, 63–90 (2017)
  * [13] Liu, D., Li, W., Vong, S.-W.: The tensor splitting with application to solve multi-linear systems. J. Comput. Appl. Math. 330, 75–94 (2018)
  * [14] Luo, Z., Qi, L., Xiu, N.: The sparsest solutions to $Z$-tensor complementarity problems. Optim. Lett. 11, 71–482 (2017)
  * [15] Ng, M., Qi, L., Zhou, G.: Finding the largest eigenvalue of a nonnegative tensor. SIAM J. Matrix Anal. Appl. 31, 1090–1099 (2009)
  * [16] Ni, Q., Qi, L.:A quadratically convergent algorithm for finding the largest eigenvalue of a nonnegative homogeneous polynomial map. J. Global Optim. 61, 627–641 (2015)
  * [17] Pearson, K.J.: Essentially positive tensors. Int. J. Algebra 4, 421–427 (2010)
  * [18] Qi, L., Luo, Z.: Tensor Analysis: Spectral Theory and Special Tensors. SIAM, Philadephia (2017)