The general tensor regular splitting iterative method for multilinear PageRank problem

1.   Introduction

As the most crucial algorithm for web ranking in the web search engines, Google's PageRank has been extensively researched over the years ^{[1,2,3,4,5,6,7,8]}. Gleich et al. ^[9] extended PageRank to higher-order Markov chains and proposed the following multilinear PageRank problem:

where $\alpha\in(0, 1)$ is a parameter, ${\bf{v}}\in\mathbb{R}^{n}$ is a stochastic vector, i.e., ${\bf{v}}\geq 0, \left \|{\bf{v}}\right \|_{1} = 1$, $\mathbb{R}^{n}$ is the set of all $n$-dimensional real vectors. The stochastic solution ${\bf{x}}\in\mathbb{R}^{n}$ is called the multilinear PageRank vector. Here $\mathcal{P} = (p_{i_{1}, i_{2}, \cdots, i_{m}})$ $(\forall i_{2}, \cdots, i_{m} \in \langle n\rangle)$ is an $m^{th}$-order $n$-dimensional stochastic tensor representing an $(m-1)^{th}$-order Markov chain, namely,

where $\langle n\rangle = \{1, 2, \cdots, n\}$. The tensor-vector product $\mathcal{P}{\bf{x}}^{m-1}$ is a vector in $\mathbb{R}^{n}$ and its definition is given in Section 2.

In particular, the multilinear PageRank model (1.1) can be transformed into higher-order Markov chain when we take $\alpha = 1$ (see ^[10] for the general case). When $m = 2$, (1.1) simplifies to the classical PageRank as described in ^[11].

The model (1.1) can be transformed into

where $\bar{\mathcal{P}} = \alpha \mathcal{P}+(1-\alpha){\mathcal{V}}$ and $\mathcal{V}$ is an $m^{th}$-order $n$-dimensional tensor with $(\mathcal{V})_{i, i_2, \cdots, i_m} = v_{i}, {\forall}i, i_{2}, \cdots, i_{m}\in\langle n \rangle$. Clearly, $\bar{\mathcal{P}}$ is also a transition probability tensor.

In recent years, the multilinear PageRank problem has attracted much attention. The theory of existence and uniqueness of solutions for the multilinear PageRank problem is a fundamental basis for analyzing the convergence of algorithms. Gleich et al. ^[9] pointed to the fact that the multilinear PageRank vector is unique in (1.1) for $\alpha < \frac{1}{m-1}$. The above condition can be easily applied to compute the stochastic solution of (1.1) when $\alpha < \frac{1}{2}$. However, when $\frac{1}{2} < \alpha < 1$, especially when $\alpha \approx 1$, it becomes increasingly challenging to solve (1.1). This challenge inspired researchers to search for more general but stricter conditions, such as those discussed in studies by Li et al. ^[12,13], Huang and Wu ^[14], Fasino and Tudisco ^[15] and Liu et al. ^[16]. It should be noted that we cannot theoretically compare the merits of these conditions and there is no optimal condition in numerical experiments.

Gleich et al. first proposed Newton's method to solve this problem in ^[9]. Subsequently, Meini and Poloni ^[17] and Guo et al. ^[18] proposed the Perron Newton iterative algorithm and the multi-step Newton method, respectively. These algorithms mentioned above involve gradient computation, which have shown excellent performance. But when solving asymmetric tensor models, the difficulty of gradient computation can lead to expensive calculations. Unlike gradient algorithms, non-gradient algorithms typically do not require the gradient information of the objective function. Instead, they estimate the optimal solution of the objective function through other methods. Non-gradient algorithms have proven to be effective in solving this problem and have gained a significant amount of attention in recent years. For example, Gleich et al. ^[9] proposed the fixed-point method, the shifted fixed-point method, the inner-outer method and the inverse method. Liu et al. ^[19] introduced several relaxation methods. Other novel methods have been proposed by various scholars in an effort to achieve better convergence performances (see References ^{[20,21,22,23]} and their respective citations). Hence, the problem about how to design high-speed, robust and flexible non-gradient algorithms has become the key challenge in solving the multilinear PageRank problem.

As an important class of non-gradient algorithms, the splitting algorithm has significant applications in tensor computations. For example, Liu et al. ^[24] proposed a tensor splitting algorithm for solving tensor equations and high-order Markov chain models. Cui et al. ^[25] proposed an iterative refinement method by using higher-order singular value decomposition to solve general multilinear systems. Jiang and Li ^[26] proposed a new preconditioner for improving the AOR-type method. Cui and Zhang ^[27] proposed a new constructive way to estimate bounds on the H-eigenvalues of two kinds of interval tensors. This paper will explore a tensor splitting algorithm for solving the multilinear PageRank problem.

The main contributions of this paper can be outlined as follows:

(1) The general tensor regular splitting (GTRS) iterative method is proposed to address the multilinear PageRank problem.

(2) Additionally, we offer five typical splitting methods of the GTRS iteration.

(3) We give the convergence analysis for the proposed method based on the uniqueness condition provided by Li et al. ^[13].

(4) Several numerical examples are provided and the outcomes significantly outperform the existing ones in most cases.

The reminder of this work is organized as follows. In Section $2$, we present some essential definitions and key properties. By making use of the (weak) regular splitting of the coefficient matrix $I-\alpha\mathcal{P}{\bf{x}}^{m-2}$, we propose the GTRS iterative method in Section $3$. Section $4$ is devoted to the convergence analysis of the GTRS method. In Section $5$, some application experiments are presented to demonstrate the effectiveness of our method. Finally, the conclusion is given in Section $6$.

2.   Preliminaries

In this section, we present a brief overview of the notations and definitions that are necessary for this paper.

An $m^{th}$-order $n$-dimensional tensor $\mathcal{A}$ with $n^{m}$ entries is defined as

where $\mathbb{R}$ is the real field. $\mathcal{A}$ is called non-negative (positive) if $a_{i_{1}, \cdots, i_{m}}\geq 0$ ($a_{i_{1}, \cdots, i_{m}} > 0$).

For any two matrices $A = (a_{ij})_{n\times n}, B = (b_{ij})_{n\times n}\in\mathbb{R}^{n}$, $A\geq B$ means that $a_{ij}\geq b_{ij}$ for $i, j\in\langle n\rangle$.

Definition 1. ^[13] Let $\mathcal {A}$ be an $m^{th}$-order $n$-dimensional tensor, ${\bf{x}}$ be an $n$-dimensional vector and ${\bf{x}}_{i}$ represent the $i^{th}$ entry of ${\bf{x}}$. $\mathcal{A}{\bf{x}}^{m-r}$ is an $r^{th}$-order $n$-dimensional tensor is given by

Definition 2. ^[13] Let $\mathcal {A}$ be an $m^{th}$-order $n$-dimensional tensor and both ${\bf{x}}$ and ${\bf{y}}$ be $n$-dimensional tensors. Then we get the following definition

Particularly, when $r = 1$ and $r = 2$, we obtain $\mathcal{A}{\bf{x}}^{m-1}$ and $\mathcal{A}{\bf{x}}^{m-2}$ from (2.1). Noticed that $\mathcal{A}{\bf{x}}^{m-1} = \mathcal{A}{\bf{x}}^{m-2}{\bf{x}}$. Then, (2.2) can be expressed as

Definition 3. ^[13] Let $\mathcal {A}$ be an $m^{th}$-order $n$-dimensional tensor and both ${\bf{x}}$ and ${\bf{y}}$ be $n$-dimensional tensors. $\mathcal {A}_{xy}^{(k)}$ is an $n\times n$ matrix whose $(i, j)^{th}$- entry is defined as

It should be noted that (1.1) can be rewritten as

Definition 4. A real square matrix $\hat{A}$ is defined as an M-matrix when it satisfies the following conditions: $\hat{A}$ is a Z-matrix (i.e., all of its off-diagonal elements $\hat{a}_{ij}\leq 0$) and $\hat{A}$ can be expressed as $hI-\hat{B}$, where $\hat{B}$ is a nonnegative matrix with $\hat{b}_{ij}\geq 0$, and its spectral radius $\rho(\hat{B}) < h$.

Definition 5. ^[28] Let $O$ be a null $n\times n$ matrix and $\check{A}$, $\check{M}$ and $\check{N}$ be $n\times n$ real matrices and $\check{A} = \check{M}-\check{N}$ is a regular splitting of matrix $\check{A}$ if $\check{M}$ is nonsingular with $\check{M}^{-1}\geq O$ and $\check{N} \geq O$. Similarly, $\check{A} = \check{M}-\check{N}$ is a weak regular splitting of matrix $\check{A}$ if $\check{M}$ is nonsingular with $\check{M}^{-1}\geq O$ and $\check{M}^{-1}\check{N} \geq O$. $\check{A} = \check{M}-\check{N}$ is a convergent splitting of matrix $\check{A}$ if $\check{M}$ is nonsingular with $\rho(\check{M}^{-1}\check{N}) < 1$.

Definition 6. Let ${\bf{z}}\in \mathbb{R}^{n}$ be a stochastic vector, ${\bf{0}}\in\mathbb{R}^{n}$ and ${\bf{z}}_{+} = \max({\bf{z}}, {\bf{0}})$. We define the projection ${\rm{proj}}(\cdot)$ as follows:

It is clear that ${\rm{proj}}(\cdot)$ is a stochastic vector.

3.   The GTRS iterative method

Building upon the success of the general inner-outer iterative method (see ^[29]) for solving the typical PageRank problem, we introduce the GTRS iterative method as an efficient solution for computing the multilinear PageRank vector.

Let

be a (weak) regular splitting. Based on (3.1), the iterative method, i.e., GTRS, for (1.1) is constructed as follows:

with $0 < \phi < 1$ and ${\bf{b}} = (1-\alpha){\bf{v}}$. The scheme (3.2) is derived from the matrix splitting

where $\check{M}_{{\bf{x}}_{k}} = \bar{M}_{{\bf{x}}_{k}}-\phi \bar{N}_{{\bf{x}}_{k}}$ and $\check{N}_{{\bf{x}}_{k}} = (1-\phi) \bar{N}_{{\bf{x}}_{k}}$.

Let $\check{U}_{\bf{x}}\geq 0$ be the strictly upper-triangular matrix of the matrix $\mathcal{P}{\bf{x}}^{m-2}$, $\check{D}_{\bf{x}}\geq 0$ be a diagonal matrix and $\check{L}_{\bf{x}}\geq 0$ be a lower-triangular matrix such that

Next, we give several typical practical choices of the matrices $\bar{M}_{\bf{x}}$ and $\bar{N}_{\bf{x}}$. Derived from the above matrix splitting (3.3), those well-known iterative methods for solving (1.1) are presented below:

(1) Let $\bar{M}_{{\bf{x}}_{k}} = I$ and $\bar{N}_{{\bf{x}}_{k}} = \alpha\mathcal{P}{\bf{x}}_{k}^{m-2}$; then, we can obtain the fixed point iteration (denote by GTRS-FP):

(2) Let $\bar{M}_{{\bf{x}}_{k}} = I-\beta\mathcal{P}{\bf{x}}^{m-2}_{k}$ and $\bar{N}_{{\bf{x}}_{k}} = (\alpha-\beta)\mathcal{P}{\bf{x}}^{m-2}_{k}$; then, we get another method, i.e., the inner-outer iteration (denote by GTRS-IO):

where $0 < \beta < \alpha$.

(3) Let $\bar{M}_{{\bf{x}}_{k}} = I-\alpha \check{D}_{{\bf{x}}_{k}}$ and $\bar{N}_{{\bf{x}}_{k}} = \alpha(\check{L}_{{\bf{x}}_{k}}+\check{U}_{{\bf{x}}_{k}})$; we derive the Jacobi iteration (denote by GTRS-JCB):

(4) Let $\bar{M}_{{\bf{x}}_{k}} = I-\alpha \check{D}_{{\bf{x}}_{k}}-\alpha \check{U}_{{\bf{x}}_{k}}$ and $\bar{N}_{{\bf{x}}_{k}} = \alpha \check{L}_{{\bf{x}}_{k}}$; then, we can get the Gauss-Seidel iteration (denote by GTRS-GS):

(5) Let $\bar{M}_{{\bf{x}}_{k}} = \frac{1}{\omega}\left(I-\alpha \check{D}_{{\bf{x}}_{k}}-\omega\alpha \check{L}_{{\bf{x}}_{k}} \right) \text { and } \bar{N}_{{\bf{x}}_{k}} = \frac{1}{\omega}\left[(1-\omega)(I-\alpha \check{D}_{{\bf{x}}_{k}})+\omega\alpha \check{U}_{{\bf{x}}_{k}}\right]$; then, we have the successive overrelaxation iteration (denote by GTRS-SOR):

Finally, the whole algorithm of the GTRS iterative method can be outlined as follows:

4.   Convergence analysis

In this section, we analyze the convergence of the GTRS iterative method. We establish the convergence of the GTRS iteration by using the uniqueness condition of the solution to the multilinear PageRank problem in ^[13]. At first, we give several lemmas which will be used later.

Lemma 1. ^[13] Let ${\bf{x}}$ and ${\bf{y}}$ be $n$-dimensional stochastic vectors, and let $\mathcal{J}^{(k)}(k = 2, 3, \cdots, m)$ be $m^{th}$-order $n$-dimensional tensors defined as

where $\sigma_{i_{1}, i_{2}, \cdots, i_{k-1}, i_{k+1}, \cdots, i_{m}}^{(k)}\in\mathbb{R}$ for any $i_{l}\in\langle n\rangle$, $l = 1, 2, \cdots, k-1, k+1, \cdots, i_{m}$ and $k = 2, 3, \cdots, m$; then,

where $\Delta {\bf x} = {\bf x}-{\bf y}$.

Lemma 2. ^[13] Let $\mathcal{P}$ be an $m^{th}$-order stochastic tensor and ${\bf{v}}$ be a stochastic vector; then, the multilinear PageRank model $(1.1)$ has the unique solution if

where $\mu = \min\limits_{\mathcal{J}^{(k)}, k = 2, 3, \cdots, m} \mu\left(\mathcal{J}^{(2)}, \cdots, \mathcal{J}^{(m)}\right) \mathit{\text{and}}\nu = \min\limits_{\mathcal{J}^{(k)}, k = 2, 3, \cdots, m} \nu\left(\mathcal{J}^{(2)}, \cdots, \mathcal{J}^{(m)}\right).$

Lemma 3. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, and let ${\bf{y}}$ and ${\bf{z}}$ be $n$-dimensional stochastic vectors; then, we have

where $\gamma=\min \{\bar{\gamma}, \breve{\gamma}\}$ with $\bar{\gamma}=1-\min\limits_{i_{3}, \cdots, i_{m} \in\langle n\rangle} \sum\limits_{i=1}^{n} \min \limits_{i_{2} \in\langle n\rangle} p_{i, i_{2}, \cdots, i_{m}}$ and $\breve{\gamma}=\max \limits_{i_{3}, \cdots, i_{m} \in\langle n\rangle} \sum\limits_{i=1}^{n} \max \limits_{i_{2} \in\langle n\rangle} p_{i, i_{2}, \cdots, i_{m}}-1$.

Proof. Let $\Delta{\bf y} = {\bf y}-{\bf z}$ and $\Delta y_{i}$ be the $i$-th entry of $\Delta{\bf y}$, and let $\mathcal{J}^{(2)}$ be defined in Lemma 1; we have

Substituting $\sigma_{i_{1}, i_{3}, \cdots, i_{m}}^{(2)} = \min\limits_{i_{2}} p_{i_{1}, i_{2}, \cdots, i_{m}}$ and $\sigma_{i_{1}, i_{3} \cdots, i_{m}}^{(2)} = \max\limits_{i_{2}} p_{i_{1}, i_{2}, \cdots, i_{m}}$ into (4.1) respectively, it follows that

and

where

Let $\gamma = \min \{\bar{\gamma}, \breve{\gamma}\}$; we have

This completes the proof. □

Lemma 4. ^[13] Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, and let ${\bf{x}}$, ${\bf{y}}$ and ${\bf{z}}$ be $n$-dimensional stochastic vectors. For any set $\mathcal{J}^{(k)}(k = 3, \cdots, m)$, we have

where

Lemma 5. ^[19] Let $\hat{{\bf{y}}} = (\hat{y}_{i})\in\mathbb{R}^{n}$ with $\sum_{i = 1}^{n}\hat{y}_{i} = 1;\; {\bf{z}} = (z_{i})$ is a stochastic vector. If ${\bf{y}} = {\rm{proj}}(\hat{{\bf{y}}})$, then $\left \|{\bf{y}}-{\bf{z}}\right\|_{1}\leq \left \|\hat{{\bf{y}}}-{\bf{z}}\right\|_{1}$.

Lemma 6. Let $\hat{{\bf{y}}} = (\hat{y}_{i})\in\mathbb{R}^{n}; \; {\bf{z}} = (z_{i})$ is a stochastic vector. If ${\bf{y}} = {\rm{proj}}(\hat{{\bf{y}}})$, then $\left \|{\bf{y}}-{\bf{z}}\right\|_{1}\leq 2\left \|\hat{{\bf{y}}}-{\bf{z}}\right\|_{1}$.

Proof. Since $\textbf{y} = {\rm{proj}}(\hat{\textbf{y}}) = \frac{\hat{{\bf{y}}}_{+}}{\left \|\hat{{\bf{y}}}_{+}\right \|_{1}}$, it follows that

The proof is completed. □

Theorem 1. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, $\alpha < \frac{1}{\min \{\mu, \nu\}}$ and ${\bf{x}}$ be the exact solution of (1.1). Then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.2) is convergent and

where $\mu_{{\bf{x}}_{k}} = \frac{(1-\phi)\left\|\bar{M}_{{\bf{x}}_{k}}^{-1}\bar{N}_{{\bf{x}}_{k}}\right\|_{1} +\alpha\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)\left\|\bar{M}^{-1}_{{\bf{x}}_{k}}\right\|_{1}} {\left (1-\phi\left\|\bar{M}_{{\bf{x}}_{k}}^{-1}\bar{N}_{{\bf{x}}_{k}}\right\|_{1}\right)}$ and $\delta_{{\bf{x}}_{k}} = 2\mu_{{\bf{x}}_{k}}$.

Proof. By (3.2), we have

Note that ${\bf{x}}$ is the solution of (1.1); hence,

Let $\hat{{\bf{e}}}_{k} = {\bf{y}}_{k}-{\bf{x}}$ and ${\bf{e}}_{k} = {{\bf{x}}}_{k}-{\bf{x}}$. Subtracting (4.3) from (4.2) yields

Hence

By taking the $1$-norm of both sides of (4.4) and Lemma 4 we get

Thus

From (4.5) and Lemma 6, we obtain

Then the GTRS iteration converges linearly with the convergence rate $\delta_{{\bf{x}}_{k}}$ when $\mu_{{\bf{x}}_{k}} < \frac{1}{2}$. □

Remark 1. In Theorem 1, the value of $\delta_{{\bf{x}}_{k}}$ is associated with ${\bf{x}}_{k}$, which is difficult to prove in practical scenarios. Therefore, we will provide convergence analysis for specific splittings of the GTRS iterative method, such as (3.4)–(3.8).

Theorem 2. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, ${\bf{x}}$ be a solution of (1.1) , and $\alpha < \min\{\frac{1}{\phi\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+\gamma}, \frac{1}{\min\{\mu, \nu\}}\}$; then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.4) is convergent and

where $\varsigma = \frac{ \phi\alpha\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+(\alpha-\phi\alpha)\gamma}{1-\phi\alpha\gamma}.$

Proof. It is easy to verify that $\alpha < \frac{1}{\min \{\mu, \nu\}}$ when $\alpha < \min\{\frac{1}{\phi\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+\gamma}, \frac{1}{\min\{\mu, \nu\}}\}$. According to Lemma 2, the solution ${\bf{x}}$ is unique.

Set $\hat{{\bf{e}}}_{k} = {\bf{y}}_{k}-{\bf{x}}$ and ${\bf{e}}_{k} = {{\bf{x}}}_{k}-{\bf{x}}$. From (3.4), it follows that

Note that ${\bf{x}}$ is the solution of (1.1) yields

Subtracting (4.8) from (4.7), we have

By (4.7), we have

On the other hand

Therefore, we get

Combining (4.9) and Lemmas 3 and 4 together gives

It is noted that

Then by (4.10) and Lemma 5, we get

By the above proof, we have that $\varsigma < 1$; then, the proof is completed. □

Theorem 3. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, ${\bf{x}}$ be a solution of (1.1) and $\alpha < \min\{\frac{1-(1-\phi)\beta\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}{\phi\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+\gamma}, \frac{1}{\min\{\mu, \nu\}}\}$; then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.5) is convergent and

where $\zeta = \frac{ (\beta+\phi\alpha-\phi\beta)\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+(1-\phi)(\alpha-\beta)\gamma}{1-(\beta+\phi\alpha-\phi\beta)\gamma}.$

Proof. Clearly, we have the unique solution ${\bf{x}}$ if $\alpha < \min\{\frac{1-(1-\phi)\beta\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}{\phi\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+\gamma}, \frac{1}{\min\{\mu, \nu\}}\}$ from {Lemma 2}.

By (3.5), we have

Let $\hat{{\bf{e}}}_{k} = {\bf{y}}_{k}-{\bf{x}}$ and ${\bf{e}}_{k} = {{\bf{x}}}_{k}-{\bf{x}}$. Then, we can get the following equation:

By the same proof as in Theorem 2, we get

Similarly, we have

Hence

By combining (4.11) with Lemmas 3 and 4, we get

which implies that

By (4.12) and Lemma 5, we have

It is easy to check that $\zeta < 1$ and the proof is completed. □

Theorem 4. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, ${\bf{x}}$ be a solution of (1.1) and $\alpha < \min\{\frac{1}{2\gamma(1-\phi)+3\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}, \frac{1}{\min\{\mu, \nu\}}\}$; then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.6) is convergent and

where $\xi = 2\frac{(1-\phi)\alpha\gamma+\alpha\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}{1-\alpha\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}.$

Proof. By {Lemma 2} and the condition that $\alpha < \min\{\frac{1}{2\gamma(1-\phi)+3\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}, \frac{1}{\min\{\mu, \nu\}}\}$, we know that x is unique.

Let $\hat{{\bf{e}}}_{k} = {\bf{y}}_{k}-{\bf{x}}$ and ${\bf{e}}_{k} = {{\bf{x}}}_{k}-{\bf{x}}$. By (3.6), we obtain

then we obtain

Combining this with Lemmas 3 and 4 leads to

Thus

By (4.13) and Lemma 6, we have

Notice that $\xi < 1$. This completes the proof of the theorem. □

Theorem 5. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, ${\bf{x}}$ be a solution of (1.1) and $\alpha < \min\{\frac{1}{2\gamma(1-\phi)+3\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}, \frac{1}{\min\{\mu, \nu\}}\}$; then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.7) is convergent and

Proof. It is similar to the proof of Theorem 4. □

Theorem 6. Let $\mathcal{P}$ be an $m^{th}$-order $n$-dimensional stochastic tensor, ${\bf{x}}$ be a solution of (1.1) , $\alpha < \min\{\frac{2\omega+(1-\omega)\phi}{3(1-\phi)+\omega\gamma+2\omega\bar{\mu}+\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}, \frac{1}{\min\{\mu, \nu\}}\}$ and $\omega\in(\frac{1-\phi}{2-\phi}, 1)$; then, for an arbitrary initial stochastic guess ${\bf{x}}_{0}$, the iterative sequence $\left\{{\bf{x}}_{k}\right\}_{k = 0}^{\infty}$ generated by (3.8) is convergent and

where $\vartheta = 2\frac{\omega\alpha\bar{\mu}\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)+(1-\phi)(1-\omega)(1+\alpha)+\omega\alpha(1-\phi)}{1-\phi(1-\omega)-\alpha(1-\phi)-\omega\alpha\gamma}.$

Proof. According to Lemma 2, it is obvious that the solution ${\bf{x}}$ is unique when $\alpha < \min\{\frac{2\omega+(1-\omega)\phi}{3(1-\phi)+\omega\gamma+2\omega\bar{\mu}+\left(\mathcal{J}^{(3)}, \cdots, \mathcal{J}^{(m)}\right)}, \frac{1}{\min\{\mu, \nu\}}\}$ and $\omega\in(\frac{1-\phi}{2-\phi}, 1)$.

By (3.8), we have

Taking $\hat{{\bf{e}}}_{k} = {\bf{y}}_{k}-{\bf{x}}$ and ${\bf{e}}_{k} = {{\bf{x}}}_{k}-{\bf{x}}$, by a proof that is analogous to that for the above theorems, we get

Due to Lemma 4, the above inequality has the following estimation

Hence

According to (4.15) and Lemma 6, we have

and $\vartheta < 1$. This completes the proof. □

Remark 2. For the uniqueness conditions stated in other papers, such as those presented in ^[12,13], we can also offer corresponding convergence analysis for our proposed algorithm. Here, we omit the detailed description.

5.   Numerical experiments

In this section, we present numerical experiments to exhibit the asset of the algorithm we are proposing.

All experiments were performed by using MATLAB R2016a on a Windows 10 64 bit computer equipped with a 1.00 GHz $\text { Intel }{ }^{\circledR} \text { Core }^{\mathrm{TM}}$ i4-29210M CPU processor and 8.00 GB of RAM. We chose to use three parameters to evaluate the effectiveness of the proposed method, which are the iteration number (denoted by IT), the computing time in seconds (denoted by CPU) and the error rate (denoted by $err$), respectively. The $err$ parameter is defined by

We have opted to employ large damping factors during the computation process, and in this case, iterative methods typically face major convergence problems when solving multilinear PageRank problems. Therefore, we tested the following values for the damping parameter $\alpha$ = ${0.95, 0.99, 0.995, 0.999}$, respectively. We selected that the maximum iterative number to be 1,000, the initial value ${\bf{x}}_{0} = {\bf{v}} = \frac{1}{n}{\bf{e}}$, where ${\bf{e}}$ is an $n$-dimensional vector with all entries being 1, and the termination tolerance as $\varepsilon = 10^{-10}$. The stopping criterion for all test methods was set as $\left\|{\bf{x}}_{k}-\alpha\mathcal{P}{\bf{x}}_{k}^{m-1}-{\bf{b}}\right\|_{1} < \varepsilon$. We searched the parameter $\phi$ (or $\omega$) from 0.1 to 1 with the step size 0.01 and $\beta\in(0, \alpha)$ from 0.1 to $\alpha$ with the interval of 0.01; also, we set

Next, we provide the numerical analysis for the GTRS iteration discussed in Section 3, using four examples. We compare our approach with the relaxation methods presented by Liu et al. ^[19] and the TARS method. The TARS iterative method is described in Algorithm $2$. To clarify, we denote Algorithms 1–4 in the work of Liu et al. as $Al_{1}, Al_{2}, Al_{3}$ and $Al_{4}$ respectively.

Example 1. In this example, we considered three tensors from the practical problems in ^[30,31]. Example (i) involves DNA sequence data, while Example (ii) applies interpersonal relationship data, and Example (iii) involves physicists' occupational mobility data. The numerical results are listed in Tables 1–3, where the minimum CPU times in each row are indicated in bold font.

$({\rm{i}})$

$({\rm{ii}})$

$({\rm{iv}})$

From Tables 1–3, we have the following observations and remarks:

(1) Although our algorithms require more iterative steps, the computational times are typically shorter than those of $Al_{1}, Al_{2}, Al_{3}$ and $Al_{4}$ and TARS, provided that the parameters are chosen appropriately.

(2) Considering the least CPU times obtained from GTRS-FP, GTRS-IO, GTRS-JCB, GTRS-GS and GTRS-SOR, we observe that they account for approximately $0\%$, $91.7\%$, $0\%$, $0\%$ and $8.3\%$, respectively. Therefore, in this example, GTRS-IO appears to be the most efficient algorithm among the five tested ones.

(3) Among the five types of splitting in the GTRS iteration, GTRS-JCB may result in higher CPU time usage compared to other methods, but it exhibits fewer iterative steps. This indicates that in cases with fewer iterations, GTRS-JCB can converge faster, while other methods may require more iterations to achieve the same convergence performance.

Next, we demonstrate the relationship between the iteration count and different values of the parameters $\alpha$ and $\phi$ in the GTRS iteration based on the inner-outer splitting approach, using Example 1(ⅱ). By choosing $\beta = \{0.2, 0.4, 0.6, 0.8\}$ and $\phi$ from 0.1 to 1 in the interval of 0.01, we obtain Figure 1. It is shown that the GTRS iteration requires fewer iteration steps for a larger $\phi$, which is more evident for a larger $\alpha$, such as when $\alpha = 0.999$. It is clear that the iteration number changes significantly depending on the chosen $\beta$, especially with $\beta = 0.2$.

Example 2. ^[12,19] Let $\Gamma$ be a directed graph with node set $\mathbb{V} = \langle n\rangle = \mathbb{D}\cup\mathbb{P}$, where $\mathbb{D}$ and $\mathbb{P}$ denote the sets of dangling nodes and pairwise connected ones, respectively. Denote $n_{p}\; (n_{p}\geq 2)$ as the number of nodes in $\mathbb{P}$, that is, $\mathbb{P} = \{1, 2, ..., n_{p}\}$. A nonnegative tensor $\mathcal{A}$ is constructed as follows:

where $a_{i_{1} i_{2} \cdots i_{l}}$ is taken randomly in (0, 1). Normalizing the entries with $p_{i_{1} i_{2} \cdots i_{l}} = \frac{a_{i_{1} i_{2} \cdots i_{l}}}{\sum_{i_{1} = 1}^{n} a_{i_{1} i_{2} \cdots i_{l}}}$ yields a stochastic tensor $\mathcal{P} = p_{i_{1} i_{2} \cdots i_{l}}$.

Based on the data presented in Table 4, there are observations that can be made as follows:

(1) The GTRS iteration demonstrates significantly shorter CPU time than those for $Al_{1}$–$Al_{4}$ and TARS. This suggests that the GTRS iterative method demonstrates higher efficiency in terms of computational time when appropriate parameters are chosen.

(2) In the case of all tested algorithms for a larger $n_{p}$, such as the case of $n_{p} = 120$, GTRS-SOR exhibits shorter CPU times than other methods.

(3) For $n = 60$, in terms of the CPU times, GTRS-IO achieves the best performance when $\alpha$ = $0.9$ or 0.95, while GTRS-FP, GTRS-JCB, GTRS-GS and GTRS-SOR outperform GTRS-IO when $\alpha$ = $0.99$ or 0.999. This implies that the behavior and reliability of these methods do not exhibit a monotonic relationship with $\alpha$.

Furthermore, we display the test results of the relationship between CPU time and $err$ at each step in Figure 2 with $n = 150$ and $n_{p} = 120$. The values of $err$ are taken as the logarithm of base-10. It is known from Figure 2 that the GTRS iterative method performs better than $Al_{1}$–$Al_{4}$ and TARS in terms of both CPU time and error measures.

Example 3. ^[13] We generated third-order stochastic tensors with dimension $n$ of 100,150 and 200 using the function rand of MATLAB. The results are reported in Table 5.

From Table 5, we observe the following:

(1) If the parameters are chosen suitably, the GTRS iterative method can achieve shorter CPU times than $Al_{1}$–$Al_{4}$ and TARS.

(2) For GTRS-JCB, the optimal performance is achieved when $\alpha = 0.99$ or 0.999, and $n = 100$ or 150. Similarly, GTRS-GS performs best when $n = 150$, while GTRS-SOR achieves the best results when $n = 200$. Notably, GTRS-IO shows the best performance when $\alpha = 0.95$. GTRS-FP requires the least iterative steps.

(3) In terms of CPU times, GTRS-IO outperforms GTRS-JCB when $\alpha = 0.9$ or 0.95. However, GTRS-JCB, GTRS-GS and GTRS-SOR outperform GTRS-IO when $\alpha = 0.99$ or 0.999.

Similarly, we present the test results for the correlation between CPU times and $err$ at each iterative step for $n = 200$ in Figure 3.

As depicted in Figure 3, it is evident that the GTRS iterative method consistently outperforms $Al_{1}$–$Al_{4}$ and TARS in terms of both CPU time and error measurement.

Example 4. ^[13] Let $t$ denote the density of zero entries in a tensor. $\mathcal{P}$ is a positive tensor or a zero tensor when $t = 0$ or $t = 1$ respectively. Test sparse tensors with different $t$ are generated by the function $randsample$ of MATLAB.

In Example 4, we generated third-order tensors of dimension 200 with densities ranging from 0.1 to 0.9 in intervals of 0.1. The numerical results are presented in Table 6. The conclusions drawn from Table 6 are as follows:

(1) When suitable parameters are found, our proposed algorithm consistently outperforms the methods presented in ^[19] and the TARS method.

(2) Out of the 36 test cases, in terms of the minimum CPU times, GTRS-SOR, GTRS-GS, GTRS-IO, GTRS-JCB and GTRS-FP occupy about $36.1\%$, $0\%$, $44.5\%$, $19.4\%$ and $0\%$, respectively. Based on these results, GTRS-SOR appears to be the most competitive algorithm, particularly when $t$ = 0.4 or 0.6.

(3) When $t$ is less than 0.5, GTRS-JCB has better convergence performance than GTRS-GS and GTRS-FP. In most cases, GTRS-IO demonstrates a shorter CPU time than the other splitting methods of the GTRS iterative type, especially when $t$ = 0.8.

Next, we discuss the relationship between the iterative steps and $\phi$ with different values of $\alpha$. Given the superior convergence properties of the GTRS iterative method with successive overrelaxation splitting, and the significant variation in the number of iterations in Example 4, we selected this splitting for the experiment. Specifically, we conducted the experiment by using $\omega = \{0.3, 0.5, 0.7, 0.9\}$, $t = \{0.4, 0.6\}$ and varying $\phi$ from 0.1 to 1 in increments of 0.01. We have plotted the associated comparison results in Figures 4 and 5. The results indicate that increasing $\phi$ leads to better performance of the GTRS iterative method, as evidenced by the reduction in the number of iterations required.

Figures 6–9 show the convergence history of the GTRS iteration under different densities of tensors and values of $\alpha$. From Figures 6–9, we find that the proposed algorithms converge faster than the existing methods in the most cases, while GTRS-IO performs the best. It is noticed that GTRS-FP, GTRS-IO and GTRS-JCB always require less computing time than the existing methods. This is due to the fact that they compute smaller errors, which reduce the number of iterations required for convergence.

Different parameter settings can lead to optimal performance for each GTRS iteration. By selecting the appropriate values for $\alpha$ and $n$, the computational efficiency and accuracy of the GTRS iteration can be improved. It is important to note that these optimal parameter settings may vary depending on the specific problem and dataset being analyzed. In general, the GTRS iteration can achieve better convergence performance than $Al_{1}$–$Al_{4}$ and TARS in terms of solving the multilinear PageRank problem.

6.   Conclusions

In this paper, based on the (weak) regular splitting of the matrix $I-\alpha\mathcal{P}{\bf{x}}^{m-2}$, we presented the GTRS iterative method for solving (1.1) from the perspective of the multilinear PageRank problem, and we have proved its overall convergence. In addition, we provided several splits of the GTRS iteration and constructed the corresponding convergence theorems. Numerical experiments indicated that the GTRS iterative method is superior to the existing ones in ^[19] and the TARS method as a result of choosing appropriate parameters. Overall, this research contributes to the advancement of non-gradient algorithms and their practical implications in tensor computations. In future work, we plan to investigate the latest uniqueness conditions to establish better convergence properties.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence tools in the creation of this paper.

Acknowledgments

This work was supported in part by the Postgraduate Educational Innovation Program of Guangdong (No. 2021SFKC030) and Guangzhou Basic and Applied Basic Research Project (No. 2023A04J1322).

The authors gratefully acknowledge the anonymous reviewers for their constructive and valuable suggestions, which greatly improved this paper.

Conflict of interest

The authors declare no conflict of interest.