Solving Multi-linear Systems with $\mathcal {M}$-Tensors

Abstract

This paper is concerned with solving some structured multi-linear systems, especially focusing on the equations whose coefficient tensors are $\mathcal {M}$-tensors, or called $\mathcal {M}$-equations for short. We prove that a nonsingular $\mathcal {M}$-equation with a positive right-hand side always has a unique positive solution. Several iterative algorithms are proposed for solving multi-linear nonsingular $\mathcal {M}$-equations, generalizing the classical iterative methods and the Newton method for linear systems. Furthermore, we apply the $\mathcal {M}$-equations to some nonlinear differential equations and the inverse iteration for spectral radii of nonnegative tensors.

Appendix

Consider the following the ordinary differential equation with Dirichlets boundary condition

$$\begin{aligned} \left\{ \begin{array}{l} u(x)^{m-2} \cdot u''(x) = -f(x),\ x \in (0,1), \\ u(0) = g_0,\ u(1) = g_1, \end{array}\right. \end{aligned}$$

where $f(x) > 0$ in (0, 1) and $g_0,g_1>0$.

Partition the interval [0, 1] into $n-1$ small intervals with the same length $h=1/(n-1)$, and denote

$$\begin{aligned} \begin{aligned} \mathbf{u}_h&= \Big [u(0),u(h),\dots ,u\big ((n-1)h\big )\Big ]^\top , \\ \mathbf{v}_h&= \Big [-u^{(4)}(0),-u^{(4)}(h),\dots ,-u^{(4)}\big ((n-1)h\big )\Big ]^\top , \\ \mathbf{f}_h&= \Big [g_0/h^2,f(h),\dots ,f\big ((n-2)h\big ),g_1/h^2\Big ]^\top , \end{aligned} \end{aligned}$$

where u(x) is the exact solution of the above boundary-value problem.

The discretization tensor $\mathcal {L}_h$ of the operator $u \mapsto u^{m-2} \cdot u''$ is introduced in the first section, and our numerical solution $\widehat{\mathbf{u}}_h$ is obtained by solving the unique positive solution of an $\mathcal {M}$-equation $\mathcal {L}_h \widehat{\mathbf{u}}_h^{m-1} = \mathbf{f}_h$. It is well-known that the truncated error of the discretization

$$\begin{aligned} \textstyle \frac{u(x-h) - 2u(x) + u(x+h)}{h^2} - u''(x) =\frac{h^2}{12} u^{(4)}(x) + O(h^4). \end{aligned}$$

Thus we have

$$\begin{aligned} \mathcal {L}_h \mathbf{u}_h^{m-1} - \mathcal {L}_h \widehat{\mathbf{u}}_h^{m-1} = \mathcal {L}_h \mathbf{u}_h^{m-1} - \mathbf{f}_h = \textstyle \frac{h^2}{12} \cdot \mathbf{u}_h^{[m-2]} \circ \mathbf{v}_h + O(h^4), \end{aligned}$$

which further implies that

$$\begin{aligned} d_h(\mathbf{u}_h,\widehat{\mathbf{u}}_h) := \big \Vert \mathcal {L}_h \mathbf{u}_h^{m-1} - \mathcal {L}_h \widehat{\mathbf{u}}_h^{m-1}\big \Vert _\infty \le \textstyle \frac{h^2}{12} \cdot \Vert u\Vert _{L^\infty }^{m-2} \cdot \Vert u^{(4)}\Vert _{L^\infty } + O(h^4). \end{aligned}$$

It can be verified that $d_h(\cdot ,\cdot )$ is a metric in the cone $\{\mathbf{x}>\mathbf{0}:\, \mathcal {L}_h \mathbf{x}^{m-1} > \mathbf{0}\}$. Then we can say that the numerical solution $\widehat{\mathbf{u}}_h$ is very close to the exact solution $\mathbf{u}_h$ when the parameter h is small enough. Next, we shall estimate the convergence of the discretization scheme.

Note that $\mathcal {L}_h \mathbf{u}_h^{m-1}$ is also a positive vector when h is small enough, then the matrix $\mathcal {L}_h \mathbf{u}_h^{m-2}$ is a nonsingular M-matrix as discussed in Sect. 4. Hence we have the first order approximation

$$\begin{aligned} \mathbf{u}_h-\widehat{\mathbf{u}}_h \approx (\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1} (\mathcal {L}_h \mathbf{u}_h^{m-1} - \mathcal {L}_h \widehat{\mathbf{u}}_h^{m-1}) \end{aligned}$$

when h is small enough, and thus

$$\begin{aligned} \Vert \mathbf{u}_h-\widehat{\mathbf{u}}_h\Vert _\infty \lesssim \big \Vert (\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1}\big \Vert _\infty \cdot \big \Vert \mathcal {L}_h \mathbf{u}_h^{m-1} - \mathcal {L}_h \widehat{\mathbf{u}}_h^{m-1}\big \Vert _\infty . \end{aligned}$$

We thus need to bound the $\infty$-norm of the inverse of the M-matrix $\mathcal {L}_h \mathbf{u}_h^{m-2}$. First denote $\mathcal {L}_h = s_h \mathcal {I} - \mathcal {A}_h$, where $s_h = 2/h^2$ and $\mathcal {A}_h$ is nonnegative. Then we can write

$$\begin{aligned} \begin{aligned} (\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1}&= \big ( (s_h \mathcal {I} - \mathcal {A}_h) \mathbf{u}_h^{m-2}\big )^{-1} \\&= \big (s_h U_h^{m-2} - \mathcal {A}_h \mathbf{u}_h^{m-2}\big )^{-1} \\&= s_h^{-1} U_h \left[ I - s_h^{-1} U_h^{-(m-1)} (\mathcal {A}_h \mathbf{u}_h^{m-2}) U_h\right] ^{-1} U_h^{-(m-1)}, \end{aligned} \end{aligned}$$

where $U_h = \mathrm{diag}\big ((\mathbf{u}_h)_1,(\mathbf{u}_h)_2,\dots ,(\mathbf{u}_h)_n\big )$.

Denote $W_h = U_h^{-(m-1)} (\mathcal {A}_h \mathbf{u}_h^{m-2}) U_h$, which is a nonnegative matrix. Note that $(\mathcal {A}_h \mathbf{u}_h^{m-2}) U_h \mathbf{1} = \mathcal {A}_h \mathbf{u}_h^{m-1}$, thus the summations of all the rows of $W_h$ are

$$\begin{aligned} \begin{aligned} W_h \mathbf{1}&= U_h^{-(m-1)} \mathcal {A}_h \mathbf{u}_h^{m-1} \le \mathbf{1} \cdot \max _{i=1:n} \frac{(\mathcal {A}_h \mathbf{u}_h^{m-1})_i }{(\mathbf{u}_h)_i^{m-1}} \\&= \mathbf{1} \cdot \max _{i=1:n} \frac{s_h (\mathbf{u}_h)_i^{m-1} - (\mathcal {L}_h \mathbf{u}_h^{m-1})_i }{(\mathbf{u}_h)_i^{m-1}} \\&= \mathbf{1} \cdot \bigg [s_h - \min _{i=1:n} \frac{(\mathcal {L}_h \mathbf{u}_h^{m-1})_i}{(\mathbf{u}_h)_i^{m-1}}\bigg ]\\&=: \mathbf{1} \cdot (s_h - \gamma _h). \end{aligned} \end{aligned}$$

Similarly, we have $W_h^k \mathbf{1} \le W_h^{k-1} \mathbf{1} \cdot (s_h - \gamma _h) \le \dots \le \mathbf{1} \cdot (s_h - \gamma _h)^k$. Also, because $W_h \mathbf{1} \le \mathbf{1} \cdot (s_h - \gamma _h) < \mathbf{1} \cdot s_h$ and $W_h$ is an irreducible nonnegative matrix, we have $\rho (s_h^{-1} W_h) < 1$. Employ the Taylor expansion of the matrix $(I-X)^{-1}=I+X+X^2+\dots$ for $\rho (X) < 1$, and we can obtain that

$$\begin{aligned} \big (I - s_h^{-1} W_h\big )^{-1} \mathbf{1} = \sum _{k=0}^\infty \big (s_h^{-1} W_h\big )^k \mathbf{1} \le \sum _{k=0}^\infty \mathbf{1} \cdot (1 - \gamma _h/s_h)^k = \mathbf{1} \cdot (s_h/\gamma _h). \end{aligned}$$

Finally, we get a upper bound of the $\infty$-norm of the nonnegative matrix $(\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1}$ that

$$\begin{aligned} \begin{aligned} \Vert (\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1}\Vert _\infty&= \max _{i=1:n} \big ((\mathcal {L}_h \mathbf{u}_h^{m-2})^{-1} \mathbf{1}\big )_i \\&= \max _{i=1:n} \big (s_h^{-1} U_h \big (I - s_h^{-1} W_h\big )^{-1} U_h^{-(m-1)} \mathbf{1}\big )_i \\&\le \Big (\min _{i=1:n} \mathbf{u}_h\Big )^{-(m-1)} \cdot \max _{i=1:n} \frac{(\mathbf{u}_h)_i^{m-1}}{(\mathcal {L}_h \mathbf{u}_h^{m-1})_i} \cdot \max _{i=1:n} \mathbf{u}_h \\&\approx \Big (\min _{i=1:n} \mathbf{u}_h\Big )^{-(m-1)} \cdot \max _{i=1:n} \frac{(\mathbf{u}_h)_i^{m-1}}{(\mathbf{f}_h)_i} \cdot \max _{i=1:n} \mathbf{u}_h \\&\le \frac{\max _x u(x)^m}{\min _x u(x)^{m-1} \cdot \min _x f(x)}. \end{aligned} \end{aligned}$$

Note that u(x) can also be regarded as the solution of the elliptic problem

$$\begin{aligned} \left\{ \begin{array}{l} u''(x) = f_1(x) := -f(x)/u(x)^{m-2},\ x \in (0,1), \\ u(0) = g_0,\ u(1) = g_1. \end{array}\right. \end{aligned}$$

So we know that $u(x) \ge \min \{g_0,g_1\}$ since $f(x) > 0$ and $g_0,g_1 > 0$. Then

$$\begin{aligned} \Vert \mathbf{u}_h-\widehat{\mathbf{u}}_h\Vert _\infty \lesssim \textstyle \frac{\Vert u\Vert _{L^\infty }^m}{\min \{g_0,g_1\}^{m-1} \cdot \min _x f(x)} \cdot \textstyle \frac{h^2}{12} \cdot \Vert u\Vert _{L^\infty }^{m-2} \cdot \Vert u^{(4)}\Vert _{L^\infty } =: K h^2, \end{aligned}$$

where the constant K is independent with the parameter h.

Therefore, the sequence $\{\widehat{\mathbf{u}}_h\}$ converges to the exact solution when $h \rightarrow 0$.