Ergodicity coefficients for higher-order stochastic processes

TL;DR

This paper introduces ergodicity coefficients for stochastic tensors to analyze the long-term behavior of higher-order stochastic processes, ensuring convergence and uniqueness of stationary distributions in complex models.

Contribution

It extends ergodicity coefficients from matrices to tensors, providing explicit formulas for uniqueness, sensitivity bounds, and convergence guarantees for higher-order stochastic processes.

Findings

Guarantees uniqueness of Perron Z-eigenvectors for stochastic tensors

Provides bounds on eigenvector sensitivity to tensor changes

Ensures convergence of higher-order stochastic processes

Abstract

The use of higher-order stochastic processes such as nonlinear Markov chains or vertex-reinforced random walks is significantly growing in recent years as they are much better at modeling high dimensional data and nonlinear dynamics in numerous application settings. In many cases of practical interest, these processes are identified with a stochastic tensor and their stationary distribution is a tensor $Z$-eigenvector. However, fundamental questions such as the convergence of the process towards a limiting distribution and the uniqueness of such a limit are still not well understood and are the subject of rich recent literature. Ergodicity coefficients for stochastic matrices provide a valuable and widely used tool to analyze the long-term behavior of standard, first-order, Markov processes. In this work, we extend an important class of ergodicity coefficients to the setting of stochastic tensors. We show that the proposed higher-order ergodicity coefficients provide new explicit formulas that (a) guarantee the uniqueness of Perron $Z$-eigenvectors of stochastic tensors, (b) provide bounds on the sensitivity of such eigenvectors with respect to changes in the tensor and (c) ensure the convergence of different types of higher-order stochastic processes governed by cubical stochastic tensors. Moreover, we illustrate the advantages of the proposed ergodicity coefficients on several example application settings, including the analysis of PageRank vectors for triangle-based random walks and the convergence of lazy higher-order random walks.

Equations

$$(\mathbf{P}xy{)}_i=\sum _{j,k=1}^n{\mathbf{P}}_{ijk}{x}_j{y}_k$$

$$(\mathbf{P}x{)}_{ij}=\sum _{k=1}^n{\mathbf{P}}_{ikj}{x}_k.$$

$${\mathbf{P}}^S={\mathbf{P}}^{⟨[132]⟩},({\mathbf{P}}^S{)}_{ijk}={\mathbf{P}}_{ikj}.$$

$${\mathbf{E}}_{ijk}^L={\delta }_{ij}\text{and}{\mathbf{E}}_{ijk}^R={\delta }_{ik},\text{for all }i,j,k=1,\dots ,n.$$

$${\mathbf{P}}_{ijk}=\mathbb{P}(X(t+1)=i|X(t)=j,X(t-1)=k),$$

$$\mathbb{P}(X(t+1)=i)=\sum _{j,k}{\mathbf{P}}_{ijk}\mathbb{P}(X(t)=j,X(t-1)=k).$$

$$\{\begin{array}{ll}(Y_{t+1}{)}_{ij}={\sum }_k{\mathbf{P}}_{ijk}(Y_t{)}_{jk}& \\ (x_{t+1}{)}_i={\sum }_j(Y_{t+1}{)}_{ij}.& \end{array}$$

$$Y_{ij}=\sum _{k=1}^n{\mathbf{P}}_{ijk}{Y}_{jk}.$$

$$x_{t+1}=\mathbf{P}{x}_t{x}_{t-1},$$

$$x=\mathbf{P}xx,$$

$$(1^{T}Y{)}_j=\sum _i{Y}_{ij}=\sum _k{Y}_{jk}\sum _i{\mathbf{P}}_{ijk}=\sum _k{Y}_{jk}=(Y1{)}_j.$$

$${\tau }_d(f)=\underset{\begin{array}{c}x,y\in {\mathcal{S}}_1\\ x\ne y\end{array}}{\sup }\frac{d(f(x),f(y))}{d(x,y)}.$$

$$d_{H}(x,y)=\log\Big{(}\max_{i}\frac{x_{i}}{y_{i}}\max_{i}\frac{y_{i}}{x_{i}}\Big{)}$$

$$\tau_{H}(A)=\tanh\left(\frac{1}{4}\log\Big{(}\max_{ijhk}\frac{A_{ij}A_{hk}}{A_{ik}A_{hj}}\Big{)}\right),$$

$${\tau }_p(P)=\underset{\begin{array}{c}x,y\in {\mathcal{S}}_1\\ x\ne y\end{array}}{\sup }\frac{\Vert Px-Py{\Vert }_p}{\Vert x-y{\Vert }_p}.$$

$${\tau }_p(P)=\underset{\begin{array}{c}\Vert x{\Vert }_p=1\\ x^{T}1=0\end{array}}{\max }\Vert Px{\Vert }_p.$$

$$\Vert x-x^{\prime }{\Vert }_p\le \frac{\Vert P-P^{\prime }{\Vert }_p}{1-{\tau }_p(P)}.$$

$$\Vert x_t-x{\Vert }_p\le {\tau }_p(P{)}^t\Vert x_0-x{\Vert }_p.$$

$${\tau }_1(P)=\frac{1}{2}\underset{j,k}{\max }\sum _i|P_{ij}-P_{ik}|.$$

$${\tau }_1(P)$$

$$\begin{array}{rl}{\mathcal{T}}_L(\mathbf{P})& =\underset{x\in \text{SS}}{\max }\underset{y\in {\mathcal{Z}}_1}{\max }\Vert \mathbf{P}xy{\Vert }_1\\ {\mathcal{T}}_R(\mathbf{P})& =\underset{x\in \text{SS}}{\max }\underset{y\in {\mathcal{Z}}_1}{\max }\Vert \mathbf{P}yx{\Vert }_1\\ \mathcal{T}(\mathbf{P})& =\underset{x\in \text{SS}}{\max }\underset{y\in {\mathcal{Z}}_1}{\max }\Vert \mathbf{P}xy+\mathbf{P}yx{\Vert }_1.\end{array}$$

$${\mathcal{T}}_L(\mathbf{P})={\mathcal{T}}_R({\mathbf{P}}^S),\mathcal{T}(\mathbf{P})=2{\mathcal{T}}_L(\frac{1}{2}\mathbf{P}+\frac{1}{2}{\mathbf{P}}^S)=2{\mathcal{T}}_R(\frac{1}{2}\mathbf{P}+\frac{1}{2}{\mathbf{P}}^S).$$

$${\mathcal{T}}_L(\mathbf{P})=\underset{x\in \text{SS}}{\max }{\tau }_1(\mathbf{P}x),{\mathcal{T}}_R(\mathbf{P})=\underset{x\in \text{SS}}{\max }{\tau }_1({\mathbf{P}}^{S}x),\mathcal{T}(\mathbf{P})=\underset{x\in \text{SS}}{\max }{\tau }_1(\mathbf{P}x+{\mathbf{P}}^{S}x).$$

$${\mathcal{T}}_L(\mathbf{P})=\frac{1}{2}\underset{j,k_1,k_2}{\max }\sum _i|{\mathbf{P}}_{ij{k}_1}-{\mathbf{P}}_{ij{k}_2}|.$$

$${\mathcal{T}}_L(\mathbf{P})$$

$$\displaystyle=1-\min_{\mathcal{I}\subset[n]}\min_{j}\Big{(}\min_{k_{1}}\sum_{i\notin\mathcal{I}}\boldsymbol{P}_{ijk_{1}}+\min_{k_{2}}\sum_{i\in\mathcal{I}}\boldsymbol{P}_{ijk_{2}}\Big{)}.$$

$$\|\boldsymbol{P}xy\|_{1}=\Big{\|}\sum_{i}x_{i}P^{(i)}y\Big{\|}_{1}\leq\sum_{i}x_{i}\|P^{(i)}y\|_{1}\leq\max_{i}\tau_{1}(P^{(i)}).$$

$${\mathcal{T}}_R(\mathbf{P})$$

$$\mathcal{T}(\mathbf{P})$$

$${\mathcal{T}}_R(\mathbf{P})$$

$$\displaystyle=1-\min_{\mathcal{I}\subset[n]}\min_{j_{1},j_{2},k}\Big{(}\sum_{i\notin\mathcal{I}}\boldsymbol{P}_{ij_{1}k}+\sum_{i\in\mathcal{I}}\boldsymbol{P}_{ij_{2}k}\Big{)}$$

$$\mathcal{T}(\mathbf{P})$$

Full text

{AMS}

15B51, 65F35, 60J10, 65C40

1 Introduction

Markov processes are among the best known and most popular stochastic processes in computational mathematics and mathematics of data science. For these types of processes the state transitions only depend on the last state. This is modeled by a stochastic matrix $P$ whose entries $P_{ij}$ quantify the probability of the process of transitioning from state $j$ to state $i$. Any such a matrix leaves the simplex $\SS=\{x\geq 0:x^{T}\mathbb{1}=1\}$ invariant and the classical Brouwer’s fixed point theorem thus implies that there exists at least one stationary distribution $x=Px$ for the Markov chain described by $P$. In other words, there exists at least one eigenvector $x$ of $P$ corresponding to the eigenvalue $1$, such that $x$ has nonnegative entries that sum up to one. While Brouwer’s theorem holds in general for mappings leaving a closed convex set invariant, much more can be said for the specific case of stochastic matrices. In particular, if the Markov chain described by $P$ is ergodic, then $P$ has a unique nonnegative eigenvector $x$ in $\mathcal{S}_{1}$ which corresponds to the eigenvalue $1$, the magnitude of any other eigenvalue of $P$ is strictly smaller than one and the power method $x_{t+1}=Px_{t}$ converges to $x$, for any choice of $x_{0}\in\mathcal{S}_{1}$, with a convergence rate that depends on the largest sub-dominant eigenvalue. In a way, these properties characterize the concept of ergodic chain and the so-called ergodicity coefficients were introduced to estimate whether or not a Markov chain is ergodic without resorting to spectral properties [29, 51].

A natural extension of a Markov process is to have the state transitions depend on several previous states, rather than just the last one. While the study and application of this kind of higher-order stochastic processes has a relatively long history, see e.g., [5, 45], their interest has grown significantly in more recent years due to their ability to improve the mathematical modeling and understanding of numerous problems in data and network sciences, such as detecting communities and analyzing spreading dynamics in networks [49, 55], understanding the behavior of web browsers and drivers trajectories [6, 15], defining new clustering algorithms in data mining that exploit motifs and non-backtracking walks [8, 33], improving centrality and ranking models for networks and hypergraphs [2, 25, 38, 41] and forecasting the appearance of new links or finding missing links in networks [3, 40].

Many higher-order stochastic processes of practical interest can be modeled by hypermatrices, or tensors, with $m$ modes $\boldsymbol{P}=(\boldsymbol{P}_{i_{1},\dots,i_{m}})$. When $m=3$, for example, we have a second-order Markov chain if $\boldsymbol{P}_{ijk}$ quantifies the probability of transitioning to state $i$, given that the last state was $j$ and the previous one $k$. Another example is given by linear vertex-reinforced random walks, where the transition probability from state $j$ to state $i$ is defined as $\sum_{k}\boldsymbol{P}_{ijk}y_{k}$, for a vector $y$ which depends on the history of the states that have been visited. The stationary distribution equation for these higher-order stochastic processes boils down to the $Z$-eigenvector of the corresponding stochastic tensor. While the existence of such stationary distribution is ensured also in this setting by Brouwer’s fixed point theorem, the ergodicity of those processes is much less understood than that of matrix-based processes. By extending the wide and influential literature on ergodicity coefficients for matrices, in this work we introduce a family of higher-order ergodicity coefficients for stochastic cubical tensors and discuss how these allow us to derive new conditions on the existence, uniqueness and computability of stationary distributions for different type of higher-order stochastic processes described by tensors. In particular, second-order Markov chains and a new class of linear vertex-reinforced random walks for which, to the best of our knowledge, we provide the first convergence result for both the occupation vector and the density distribution. This class includes previously considered vertex-reinforced stochastic processes such as the spacey random walk [6].

From the linear algebraic perspective, our new conditions allow us to prove guarantees for existence, uniqueness and computability of the Perron $Z$-eigenvector of a stochastic tensor of order three. Dominant $Z$-eigenvectors of nonnegative tensors appear in a large variety of applications, including diffusion kurtosis imaging in medical engineering [44], low-rank factorization and signal processing [1, 18], quantum processing, quantum geometry and data mining [4, 7, 28]. Even though computing a prescribed $Z$-eigenvector is in general NP-hard [26], the use of higher-order ergodicity coefficients allows us to identify a class of nonnegative tensors for which the Perron eigenvector can be approximated efficiently to an arbitrary precision. While we focus here on stochastic tensors of order three, we believe the results here presented can be further extended to more general eigenvector problems for nonnegative tensors.

The remainder of the paper is structured as follows: We fix the relevant notation in the next section. In Section 3 we review the concept of higher-order Markov chain, its associated $Z$-eigenvector stationary distribution and the issues related to the ergodicity of this type of higher-order stochastic process. In Section 4 we recall the concept of ergodicity coefficient for a stochastic matrix and some of its properties. Then, in Section 5, we introduce our new higher-order ergodicity coefficients for stochastic cubic tensors and we prove some of our main results. In Section 6 we show how these apply to the ergodicity of higher-order stochastic processes. In particular, after recalling the definition of vertex-reinforced and spacey random walks, we introduce in Subsection 6.2 a general family of Markov processes with memory that includes the spacey random walk as particular case and we prove a new convergence result for this general stochastic process. In Section 7 we compare the ergodicity coefficients introduced in Section 5 with analogous coefficients found in the recent literature. Finally, in Section 8, we discuss a number of application examples that showcase the advantages of the proposed results. In particular, we consider the computation of the multilinear PageRank and its application to triangle-based random walks in networks, and the convergence of the shifted higher-order power method.

2 Notation

Let $e_{i}$ be the $i$-th canonical basis vector in $\mathbb{R}^{n}$ and let $\mathbb{1}$ be the all-ones vector. Define the sets $\SS=\{x\in\mathbb{R}^{n}:x\geq 0,\|x\|_{1}=1\}$ and $\mathcal{Z}_{1}=\{x\in\mathbb{R}^{n}:\mathbb{1}^{T}x=0,\|x\|_{1}=1\}$. A real cubical tensor $\boldsymbol{P}$ of order $3$ (or, equivalently, with $3$ modes) is a three-way array with real entries of size $n\times n\times n$. We denote by $\mathbb{R}^{[3,n]}$ the set of such tensors and use capital bold letters to denote its elements. The $(i,j,k)$-entry of $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ is denoted by $\boldsymbol{P}_{ijk}$. Matrices are tensors with only $2$ modes and are denoted with standard capital letters.

Given a tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$, we write $\boldsymbol{P}xy$ to denote the tensor-vector multiplication over the second and third modes:

[TABLE]

for $i=1,\dots,n$. Moreover, the product $\boldsymbol{P}x$ denotes the matrix associated with the linear map $y\mapsto\boldsymbol{P}xy$, that is,

[TABLE]

With this notation, it holds $(\boldsymbol{P}x)y=\boldsymbol{P}xy$. A $Z$-eigenvalue of a tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ is a real number $\lambda$ such that there exists a nonzero vector $x\in\mathbb{R}^{n}$ such that $\boldsymbol{P}xx=\lambda x$. Such vector $x$ is a $Z$-eigenvector associated with $\lambda$, see [43].

There are $6=3!$ possible transpositions of a tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$, each corresponding to a different permutation $\pi$ of the set $\{1,2,3\}$. Using the notation proposed in [47], the transposed tensor corresponding to the permutation $\pi$ can be denoted by $\boldsymbol{P}^{\left<\pi\right>}$, namely, $(\boldsymbol{P}^{\left<\pi\right>})_{ijk}=\boldsymbol{P}_{\pi(i),\pi(j),\pi(k)}.$ As it will be of particular importance to us, we devote the special notation $\boldsymbol{P}^{S}$ to denote the tensor obtained by transposing the entries of $\boldsymbol{P}$ over the second and third modes, namely

[TABLE]

Moreover, we say that a tensor $\boldsymbol{P}$ is $S$-symmetric whenever $\boldsymbol{P}=\boldsymbol{P}^{S}$.

All inequalities in this work are meant entry-wise. In particular, we write $\boldsymbol{P}\geq 0$ (resp., $\boldsymbol{P}>0$) to denote a tensor such that $\boldsymbol{P}_{ijk}\geq 0$ (resp., $\boldsymbol{P}_{ijk}>0$) for all indices $i,j,k=1,\dots,n$. A tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ is said to be column stochastic or simply stochastic, if $\boldsymbol{P}\geq 0$ and its first mode entries all sum up to one, i.e., $\sum_{i=1}^{n}\boldsymbol{P}_{ijk}=1$, $\forall j,k=1,\dots,n$. A tensor acting as the identity on the unit sphere $x^{T}x=1$ can be defined in the case of tensors with an even number of modes, see [31]. For tensors with three modes we define the following two left $\boldsymbol{E}^{L}$ and right $\boldsymbol{E}^{R}$ “one-sided identity tensors”:

[TABLE]

Both $\boldsymbol{E}^{L}$ and $\boldsymbol{E}^{R}$ are stochastic tensors and for all $x\in\SS$ and $v\in\mathbb{R}^{n}$ one has $\boldsymbol{E}^{L}vx=v$ and $\boldsymbol{E}^{R}xv=v$. Indeed, $(\boldsymbol{E}^{L}vx)_{i}=\sum_{jk}\boldsymbol{E}^{L}_{ijk}v_{j}x_{k}=\sum_{jk}\delta_{ij}v_{j}x_{k}=v_{i}\sum_{k}x_{k}=v_{i}$ and similarly for $\boldsymbol{E}^{R}$. Note that, letting $\boldsymbol{E}=\alpha\boldsymbol{E}^{L}+(1-\alpha)\boldsymbol{E}^{R}$ for any $\alpha\in[0,1]$, it holds $\boldsymbol{E}xx=x$ for all $x\in\SS$.

3 Higher-order Markov chains

Higher-order Markov chains are a natural extension of Markov chains, where the transition probabilities depend on the past few states, rather than just the last one. For a plain introduction, see e.g., [6, 56]. For example, a discrete-time second-order Markov chain is defined by a third-order tensor $\boldsymbol{P}=(\boldsymbol{P}_{ijk})$ where $\boldsymbol{P}_{ijk}$ is the conditional probability of transitioning to state $i$, given that the last state was $j$ and the second last state was $k$. More precisely, if $X(t)$ is the random variable describing the status of the chain on the set $\{1,\ldots,n\}$ at time $t=0,1,\ldots$, then

[TABLE]

where $\mathbb{P}$ denotes probability. Hence, the sequence $\{X(t)\}$ obeys the rule

[TABLE]

Obviously it must hold $\sum_{i}\boldsymbol{P}_{ijk}=1$ for $j,k=1,\ldots,n$, i.e., the tensor $\boldsymbol{P}$ is stochastic.

Let $x_{t}\in\SS$ be the probability vector of the random variable $X(t)$, i.e., the vector with entries $(x_{t})_{i}=\mathbb{P}(X(t)=i)$. Let $Y_{t}$ denote the joint probability function $(Y_{t})_{ij}=\mathbb{P}(X(t)=i,X(t-1)=j)$. Then, $x_{t}$ is the marginal probability $Y_{t}\mathbb{1}$, i.e., the vector with entries $(x_{t})_{i}=\sum_{j}(Y_{t})_{ij}$. Hence, the dynamic of the second-order Markov chain (3) is described by the two-phase process

[TABLE]

Note that both steps in (4) are linear and thus their convergence can be analyzed using standard ergodicity arguments. In fact, the second-order Markov chain over the state set $\{1,\dots,n\}$ can be easily reduced to a first-order Markov chain with state space $\{1,\dots,n\}\times\{1,\dots,n\}$, see e.g., [6, 56]. Thus, under appropriate hypotheses on $\boldsymbol{P}$, the iteration (4) has a unique limit $Y\geq 0$ such that

[TABLE]

However, this approach has a computational drawback: the size of the joint probability function of a second-order Markov chain is the square of the number of states. The situation gets even worse for an $m$-th order Markov chain due to the “curse of dimensionality” effect: the memory space required by the joint density grows exponentially with the states space size, requiring $n^{m}$ entries. Moreover, the convergence analysis of the iteration (4) and its natural extension to the $m>2$ setting becomes cumbersome.

In order to circumvent these issues, Raftery [45] proposed a technique to approximate higher-order Markov chains by means of a linear combination of first-order ones, by assuming that the joint probability distribution of the lagged random variables $X(t),\ldots,X(t-m+1)$ can be replaced by a mixture of its marginals. In the second-oder ($m=2$) case, that assumption reduces to replacing the conditional probabilities $\boldsymbol{P}_{ijk}$ by an expression of the form $\lambda Q_{ij}+(1-\lambda)Q_{ik}$, where $Q$ is a stochastic matrix and $\lambda\in[0,1]$. This technique, known as the Mixture Transition Distribution model, has been widely used to fit stochastic models with far fewer parameters than the fully parameterized model to multi-dimensional data in a variety of applications [9, 46].

A more recent approach, which maintains all the information contained in the transition tensor $\boldsymbol{P}$, is the one proposed by Li and Ng in [37]. Here, still in the $m=2$ case, one assumes that the joint probability distribution of the higher-order Markov chain is the tensor product of its marginal distributions, that is, $Y_{t}=x_{t}\,x_{t-1}^{T}$. This hypothesis, which is equivalent to assuming that the random variables $X(t)$ and $X(t-1)$ are independent, is a conceptual simplification of the Markov chain formalism that is introduced in order to obtain a computationally tractable extension to the second-order case. The resulting process is the quadratic version of a nonlinear Markov process [32] and it is still called a second-order Markov chain by many authors, see e.g., [27, 34, 37]. In this work, we will follow this well established convention.

Using our tensor-vector product notation, this “reduced” higher-order Markov process boils down to the iteration

[TABLE]

which replaces (4) and is the higher-order counterpart of the usual Markov process for a stochastic matrix in the classical (first-order) Markov chain setting. The limit of this sequence, if it exists, is a nonnegative vector $x\in\mathcal{S}_{1}$ such that

[TABLE]

that is, $x$ is a $Z$-eigenvector of $\boldsymbol{P}$ associated with the $Z$-eigenvalue $1$. Thus, it is natural to consider any such vector as a stationary density of the Markov chain (3).

Note that the limit matrix $Y$ of (4) is such that $Y\mathbb{1}=Y^{T}\mathbb{1}$. Indeed, from (5) we have

[TABLE]

But that row-column sum is generally different from the vector $x$ in (7). In fact, that solution corresponds to a case where $Y$ has rank one, namely, $Y=xx^{T}$. Indeed, if $Y$ in (5) is such that $\mathrm{rank}(Y)=1$, then $x=Y\mathbb{1}$ must solve (7).

On the other hand, the converse implication is false in general; that is, if $x$ solves (7) then the matrix $Y=xx^{T}$ may not be a solution of (5). Indeed, extensive numerical experiments reported in [56] show that the vector $x$ is strongly correlated with the row-column sum vector of $Y$, but the matrix $Y$ has full rank in general and $x\neq Y\mathbb{1}$.

3.1 Ergodicity of higher-order Markov chains

In the matrix case, a Markov chain is called ergodic whenever it has a unique stationary vector and, for any initial probability distribution, that vector is the limiting distribution of the chain. Necessary and sufficient conditions for ergodicity of Markov chains are well known, and are essentially related to spectral and structural properties, e.g., irreducibility and aperiodicity, of the transition matrix [51].

The situation complicates significantly when moving from matrices to tensors and, more generally, from linear to nonlinear cases [32]. In fact, even though the existence of a solution $x\in\SS$ to (7) is a direct consequence of the Brouwer’s fixed point theorem, the properties that characterize uniqueness and convergence of the process to the stationary distribution do not extend straightforwardly from the matrix case. For example, unlike the matrix case, the irreducibility of $\boldsymbol{P}$ is not enough to ensure the uniqueness of $x$ and additional assumptions are required. In fact, it is not too difficult to produce examples of entrywise positive stochastic tensors for which the equation (7) has multiple solutions, or the solution of (7) is unique but the iteration (6) fails to converge to that solution. For instance, a $\boldsymbol{P}\in\mathbb{R}^{[4,2]}$ example is provided by Chang and Zhang in [13], while several $\boldsymbol{P}\in\mathbb{R}^{[3,3]}$ examples are provided by Saburov in [50].

A sufficient condition that ensures ergodicity is the existence of a metric with respect to which the system is contractive. Even though, as in the linear case, this is a sufficient but not a necessary requirement in general, suitable choices of the metric can provide valuable conditions for the ergodicity of higher-order stochastic processes that can be given in terms of the entries of the tensor $\boldsymbol{P}$.

By considering the $\ell^{1}$ and the Hilbert metrics on $\SS$, in the following we introduce a family of ergodicity coefficients for stochastic cubic tensors of order three and we show, in Section 6, how they allow us to prove new conditions for the ergodicity of various higher-order stochastic processes. The conditions we obtain in this way can be easily computed and are, to the best of our knowledge, among the weakest conditions available in the literature so far.

4 Coefficients of ergoditicy

Let $d:\SS\times\SS\to\mathbb{R}_{+}$ be a metric on $\mathcal{S}_{1}=\{x\geq 0:x^{T}\mathbb{1}=1\}$ and consider a mapping $f:\mathcal{S}_{1}\to\SS$. Although other notions of ergodicity coefficient are available in the literature, see e.g., [29], for the purpose of this work a coefficient of ergodicity for $f$ is the best Lipschitz constant of $f$ with respect to $d$, that is

[TABLE]

Different choices of the metric $d$ give rise to different notions of ergodicity coefficients. For example, if $d$ is the Hilbert projective distance

[TABLE]

then (8) is the so-called Birkhoff contraction ratio [10], which we denote by $\tau_{H}(f)$. This choice of metric is particularly interesting because it extends very naturally to the case of mappings $f$ that leave the slice of a generic proper cone invariant (not just $\mathcal{S}_{1}$). Moreover, when $f$ is a linear map described by the matrix $A$, the Birkhoff–Hopf theorem [20] provides an explicit formula for $\tau_{H}(f)=\tau_{H}(A)$, which we recall below:

[TABLE]

where $\tanh(\lambda)=(e^{\lambda}-e^{-\lambda})/(e^{\lambda}+e^{-\lambda})$ denotes the hyperbolic tangent. An equivalent formula can be found also in [51, §3.4]. More recently, in [22], an analogous explicit formula has been proved for the case where $f$ is a (weakly) multilinear mapping induced by a nonnegative tensor. In particular, this formula holds for the case of $Z$-eigenvectors of cubic stochastic tensors and we will review it in this setting in Subsection 7.1.

Another popular and successful choice for the distance in (8) is $d(x,y)=\|x-y\|_{p}$, where $\|\cdot\|_{p}$ is the $p$-norm on $\mathbb{R}^{n}$. Norm-based coefficients were introduced by Dobrushin in 1956 [19] for the case of linear mappings and have been the subject of numerous investigations afterwards, see e.g., [29, 50, 51].

In Section 5 we analyze properties of norm-based coefficients for mappings defined by a stochastic tensor $\boldsymbol{P}$. To this end, we first review some relevant properties of these coefficients for the case of linear maps.

4.1 Norm-based ergodicity coefficients for matrices

Let $P$ be a stochastic matrix and $p\geq 1$. The $p$-norm ergodic coefficient of $P$ is

[TABLE]

This definition extends obviously to any matrix $P\in\mathbb{R}^{n\times n}$, when appropriate. The linearity of $P$, the continuity of $\|\cdot\|_{p}$ and the fact that the set $\{z\in\mathbb{R}^{n}:z=(x-y)/\|x-y\|_{p},\,x,y\in\mathcal{S}_{1}\}$ coincides with $\mathcal{Z}_{p}=\{z\in\mathbb{R}^{n}:\|z\|_{p}=1,\,\mathbb{1}^{T}z=0\}$, which is compact, yield the equivalent formula

[TABLE]

Norm-based ergodicity coefficients for stochastic matrices $P$ are particularly useful for three reasons: they provide sufficient conditions for the ergodicity of the Markov chain associated with $P$; they can be used to derive bounds on the variation of the stationary distribution of the Markov chain, when the transition probabilities face a small perturbation; and, in the case $p=1$ they yield easily computable upper bounds on the convergence rate of the Markov process $x_{t+1}=Px_{t}$. We review these properties in the next Theorems 4.2, 4.3 and 4.4. Then, in Sections 5 and 6, we will use the $1$-norm ergodicity coefficients for matrices to introduce what we call higher-order ergodicity coefficients for stochastic tensors of order three and we will show that the above three fundamental properties carry over to the tensor case. We refer to [29, 51, 54] for more details on $\tau_{p}(P)$.

The following properties follow directly from the definition of $p$-norm ergodic coefficient.

Theorem 4.1**.**

If $P,Q\in\mathbb{R}^{n\times n}$ are stochastic, then

$0\leq\tau_{p}(P)\leq\|P\|_{p}$** 2. 2.

$|\tau_{p}(P)-\tau_{p}(Q)|\leq\tau_{p}(P-Q)$** 3. 3.

$\tau_{p}(P)=0$* if and only if $\mathrm{rank}(P)=1$.*

Moreover, the following perturbation bound holds (see e.g. [52] or [29, Thm. 3.14]).

Theorem 4.2**.**

Let $P,P^{\prime}$ be two stochastic irreducible matrices, and let $x,x^{\prime}$ be their corresponding stationary probability vectors. Then

[TABLE]

As an immediate consequence of the definition (8), the inequality $\tau_{p}(P)<1$ implies that the map $f:\SS\mapsto\SS$ defined by $f(x)=Px$ is a contraction. This observation implies the following result.

Theorem 4.3**.**

If $P$ is a stochastic matrix with $\tau_{p}(P)<1$ for some $p\geq 1$ then $P$ is ergodic, i.e., there exists a unique eigenvector $x\in\mathcal{S}_{1}$ such that $Px=x$. Moreover, the power method $x_{t+1}=Px_{t}$ converges to $x$ for any $x_{0}\in\mathcal{S}_{1}$, and

[TABLE]

The theorem above gives a sufficient condition for the ergodicity of $P$ which is very useful in practice when combined with a number of explicit formulas that allow to compute $\tau_{p}(P)$ using only the entries of $P$. Here we recall those for the particular case $p=1$ [19].

Theorem 4.4**.**

Let $P\in\mathbb{R}^{n\times n}$. Then

[TABLE]

Moreover, if $P$ is stochastic then

[TABLE]

We will devote Sections 5 and 6 to extend the ergodicity coefficient $\tau_{1}(P)$ to three-mode tensors, to prove analogous theorems to the preceding ones and to discuss further properties and applications.

5 Ergodicity coefficients for third-order tensors

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a cubic stochastic tensor. We define the following higher-order ergodicity coefficients:

[TABLE]

The preceding definitions are extended obviously to any tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$, when appropriate. We remark the following immediate identities:

[TABLE]

In particular, for an $S$-symmetric tensor $\boldsymbol{P}$ we have $\mathcal{T}_{L}(\boldsymbol{P})=\mathcal{T}_{R}(\boldsymbol{P})=\frac{1}{2}\mathcal{T}(\boldsymbol{P})$.

The relationship between the preceding definitions and the norm-based ergodicity coefficients considered in Section 4.1 can be revealed by considering the matrices associated with the tensor-vector products $\boldsymbol{P}x$ and $\boldsymbol{P}^{S}x$ defined as in (1). In fact, it is not difficult to see that the following identities hold,

[TABLE]

The above formulas yield a characterization of the three coefficients in (10) which, for example, was used in [50] to define $\mathcal{T}(\boldsymbol{P})$ in the case of $S$-symmetric tensors. Hereafter, we exploit these formulas to derive explicit expressions for computing the coefficients above from the knowledge of the tensor entries and provide the tensor equivalent of Theorem 4.4.

Theorem 5.1**.**

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$. Then,

[TABLE]

Moreover, if $\boldsymbol{P}$ is stochastic then

[TABLE]

Proof 5.2**.**

For $i=1,\ldots,n$ let $P^{(i)}$ be the stochastic matrix $P^{(i)}_{jk}=\boldsymbol{P}_{jik}$. Hence, for any $x\in\SS$ and $y\in\mathcal{Z}_{1}$ we have $\boldsymbol{P}xy=\sum_{i}x_{i}P^{(i)}y$. By the triangle inequality,

[TABLE]

*On the other hand, if $x=e_{i}$ where $i$ is an index such that $\tau_{1}(P^{(i)})\geq\tau_{1}(P^{(j)})$ for $j=1,\ldots,n$ and $y\in\mathcal{Z}_{1}$ is a vector such that $\|P^{(i)}y\|_{1}=\tau_{1}(P^{(i)})$ then the inequalities above are actually equalities. Thus all the claims follow at once from Theorem 4.4. *

The analogous formulas for the other higher-order coefficients in (10) are derived hereafter.

Corollary 5.3**.**

*Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$. The following identities hold: *

[TABLE]

Moreover, if $\boldsymbol{P}$ is stochastic then

[TABLE]

Proof 5.4**.**

Equations (14), (16) and (17) derive from the identity $\mathcal{T}_{R}(\boldsymbol{P})=\mathcal{T}_{L}(\boldsymbol{P}^{S})$ and equations (11), (12) and (13), respectively. Now, define $\boldsymbol{Q}=\frac{1}{2}(\boldsymbol{P}+\boldsymbol{P}^{S})$. Note that if $\boldsymbol{P}$ is stochastic then also $\boldsymbol{Q}$ is stochastic. Since

[TABLE]

*we have $\mathcal{T}(\boldsymbol{P})=2\mathcal{T}_{L}(\boldsymbol{Q})$. Hence, equations (15), (18) and (19) derive from the latter identity and equations (11), (12) and (13), respectively. *

By (12), (16) and (18), it is immediate to observe that for a stochastic tensor $\boldsymbol{P}$ it holds $0\leq\mathcal{T}_{L}(\boldsymbol{P}),\mathcal{T}_{R}(\boldsymbol{P})\leq 1$ and

[TABLE]

Stronger inequalities can be easily obtained for positive tensors, as shown in the next result.

Corollary 5.5**.**

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a stochastic tensor. If there exists a positive number $\alpha>0$ such that $\boldsymbol{P}_{ijk}\geq\alpha$ for all $i,j,k$ then

[TABLE]

Proof 5.6**.**

*The three inequalities in the claim follow immediately from equations (12), (16) and (18), respectively. *

Remark 5.7**.**

A close look at Theorem 5.1 reveals that, for any tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ we have $\mathcal{T}_{L}(\boldsymbol{P})=0$ if and only if $\boldsymbol{P}_{ijk}=A_{ij}$ for some matrix $A\in\mathbb{R}^{n\times n}$. In particular, $\boldsymbol{P}$ is stochastic if and only if $A$ is stochastic. Analogously, from Corollary 5.3 we derive that $\mathcal{T}_{R}(\boldsymbol{P})=0$ if and only if $\boldsymbol{P}_{ijk}=A_{ik}$ for some matrix $A$. Consequently, $\mathcal{T}_{L}(\boldsymbol{P})+\mathcal{T}_{R}(\boldsymbol{P})=0$ if and only if $\boldsymbol{P}_{ijk}=v_{i}$ for some vector $v$. It is not difficult to prove that the latter is also equivalent to $\mathcal{T}(\boldsymbol{P})=0$. Hence, if $\boldsymbol{P}$ is nonzero, we have

[TABLE]

*In the matrix case, a coefficient of ergodicity $\tau$ is called proper when the identity $\tau(P)=0$ for a stochastic matrix $P$ is equivalent to the condition $\mathrm{rank}(P)=1$, see [29, 51]. For example, both the Birkhoff coefficient $\tau_{H}$ and all the norm-based ergodicity coefficients $\tau_{p}$ are proper. By extending that definition to the tensor case, we can say that $\mathcal{T}$ is proper, while $\mathcal{T}_{L}$ and $\mathcal{T}_{R}$ are not proper. *

The remark above shows that Property 3 of Theorem 4.1 carries over to the higher-order setting. In the next Subsection 5.1 we show that also Properties 1 and 2 of that theorem enjoy a tensor counterpart. In Subsection 6.1, instead, we show how the perturbation result of Theorem 4.2 transfers to stochastic tensors.

5.1 Bounding the variation of higher-order coefficients

When working with stochastic tensors, it is quite natural to endow $\mathbb{R}^{[3,n]}$ with the norm

[TABLE]

so that, if $\boldsymbol{P}$ is stochastic, we have $\|\boldsymbol{P}\|_{1}=1$.

With the next theorem we prove a Lipschitz-continuity condition for the higher-order ergodicity coefficients with respect to the tensor $1$-norm above.

Theorem 5.8**.**

For arbitrary $\boldsymbol{P},\boldsymbol{Q}\in\mathbb{R}^{[3,n]}$ we have

[TABLE]

where $\mathcal{T}_{*}$ is any of $\mathcal{T}_{L}$ or $\mathcal{T}_{R}$. Moreover,

[TABLE]

Proof 5.9**.**

Let $\mathcal{T}_{*}=\mathcal{T}_{L}$, the other case being completely analogous. Suppose that $\mathcal{T}_{L}(\boldsymbol{P})\geq\mathcal{T}_{L}(\boldsymbol{Q})$. Hence, for some $x\in\SS$ and $y\in\mathcal{Z}_{1}$ we have

[TABLE]

Hence, $\mathcal{T}_{L}(\boldsymbol{P})-\mathcal{T}_{L}(\boldsymbol{Q})\leq\mathcal{T}_{L}(\boldsymbol{P}-\boldsymbol{Q})$. By reversing the roles of $\boldsymbol{P}$ and $\boldsymbol{Q}$ we obtain $\mathcal{T}_{L}(\boldsymbol{Q})-\mathcal{T}_{L}(\boldsymbol{P})\leq\mathcal{T}_{L}(\boldsymbol{P}-\boldsymbol{Q})$ and we arrive at the first claim. Analogously, for some $x\in\SS$ and $y\in\mathcal{Z}_{1}$ we have

[TABLE]

*The inequality $\mathcal{T}(\boldsymbol{Q})-\mathcal{T}(\boldsymbol{P})\leq\mathcal{T}(\boldsymbol{P}-\boldsymbol{Q})$ follows from the preceding one by exchanging $\boldsymbol{P}$ and $\boldsymbol{Q}$, and the second claim follows. The rightmost inequalities follow immediately from the definition of the ergodicity coefficients. *

6 Second-order stochastic processes and $Z$-eigenvectors

In this section we prove an analogous of Theorem 4.3 for tensor $Z$-eigenvectors. Precisely, given $\boldsymbol{P}$ stochastic, we provide a new condition that ensures the existence and uniqueness of a positive vector $x\in\SS$ such that $x=\boldsymbol{P}xx$. Moreover, we show that under the same condition the higher-order power method $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$, which is the prototypical nonlinear Markov process [32], always converges to $x$ and we provide an analogous, but stronger, condition that guarantees the global convergence of the alternate scheme $x_{t+1}=\boldsymbol{P}x_{t}x_{t-1}$.

The next theorem provides the tensor analogous of Theorem 4.3.

Theorem 6.1**.**

If $\boldsymbol{P}$ is stochastic, then $\mathcal{T}(\boldsymbol{P})$ is the best Lipschitz constant of the quadratic map $f:\SS\mapsto\SS$ given by $f(x)=\boldsymbol{P}xx$, that is,

[TABLE]

Therefore, if $\mathcal{T}(\boldsymbol{P})<1$ then there exists a unique $Z$-eigenvector $x\in\SS$ such that $\boldsymbol{P}xx=x$. Moreover, the higher-order power method $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$ converges to $x$ for any $x_{0}\in\SS$, and

[TABLE]

Proof 6.2**.**

Let $f:\SS\mapsto\SS$ be given by $f(x)=\boldsymbol{P}xx$. Let $\boldsymbol{Q}=\frac{1}{2}(\boldsymbol{P}+\boldsymbol{P}^{S})$. Note that $\boldsymbol{Q}$ is a stochastic tensor such that $\boldsymbol{Q}=\boldsymbol{Q}^{S}$. Moreover, the equation $f(x)=x$ is equivalent to $\boldsymbol{Q}xx=x$. Then, for all $x,y\in\SS$ we have

[TABLE]

Hence,

[TABLE]

*Since $2\mathcal{T}_{L}(\boldsymbol{Q})=\mathcal{T}(\boldsymbol{P})$, we obtain the first part of the claim. In particular, we get $\|f(x)-f(y)\|_{1}\leq\mathcal{T}(\boldsymbol{P})\|x-y\|_{1}$ for any $x,y\in\SS$. Hence, if $\mathcal{T}(\boldsymbol{P})<1$ then $f$ is contractive with respect to the $1$-norm. By the Banach fixed point theorem, there exists a unique fixed point $x\in\SS$ such that $x=f(x)$. Moreover, the iteration $x_{t+1}=f(x_{t})$ converges to $x$ with $\|x_{t}-x\|_{1}\leq\mathcal{T}(\boldsymbol{P})\|x_{t-1}-x\|$ for any $x_{0}\in\SS$ and the proof is complete. *

We note in passing that the following result, which has been derived from a well-known uniqueness result in the fixed point theory several times by different authors [17, 36, 37], is a direct consequence of the theorem above.

Corollary 6.3**.**

*If $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ is a stochastic tensor such that $\boldsymbol{P}_{ijk}>1/(2n)$ for all $i,j,k$, then there exists a unique $Z$-eigenvector $x\in\SS$ such that $\boldsymbol{P}xx=x$ and the higher-order power method $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$ converges to $x$ for any $x_{0}\in\SS$. *

Proof 6.4**.**

*In the stated hypotheses we have $\mathcal{T}(\boldsymbol{P})<1$ by virtue of Corollary 5.5. Hence, the claim is a direct consequence of Theorem 6.1. *

Given a stochastic tensor $\boldsymbol{P}$ and two initial points $x_{0},x_{-1}\in\SS$, the following alternate higher-order power method has been considered in [27]:

[TABLE]

Note that this coincides with the second-order stochastic process described in (6). In [27] the convergence of (21) has been proven when $\boldsymbol{P}_{ijk}>1/(2n)$ and under restrictive hypotheses on the choice of $x_{0}$ and $x_{-1}$. The following theorem provides a condition in terms of $\mathcal{T}_{L}(\boldsymbol{P})$ and $\mathcal{T}_{R}(\boldsymbol{P})$ that ensures that (21) converges globally to the unique fixed point of $\boldsymbol{P}$.

Theorem 6.5**.**

Let $\boldsymbol{P}$ be a stochastic tensor and let $s=\mathcal{T}_{L}(\boldsymbol{P})+\mathcal{T}_{R}(\boldsymbol{P})$. If $s<1$ then the iteration (21) converges to the unique $Z$-eigenvector $x\in\SS$ such that $x=\boldsymbol{P}xx$. In fact, for all $x_{0},x_{-1}\in\SS$ and $t=0,1,\ldots$ it holds

[TABLE]

Proof 6.6**.**

First notice that the assumption $\mathcal{T}_{L}(\boldsymbol{P})+\mathcal{T}_{R}(\boldsymbol{P})<1$ implies $\mathcal{T}(\boldsymbol{P})<1$, thus, by Theorem 6.1, there exists a unique positive $x\in\SS$ such that $x=\boldsymbol{P}xx$. We have

[TABLE]

Thus, for any $t\geq 0$,

[TABLE]

In particular, the claim is true for $t=0$. The proof is completed by a simple inductive argument. Indeed, let $m=\max\{\|x_{0}-x\|_{1},\|x_{-1}-x\|_{1}\}$ and $\epsilon_{t}=\|x_{t}-x\|_{1}$ to simplify notations. For $t>0$, suppose the claim true up to $t-1$. Then,

[TABLE]

*and the theorem is proved. *

6.1 A perturbation result for the stochastic $Z$-eigenvector

A fundamental perturbation analysis problem is to obtain quality bounds on the variation of the ergodic distribution of the non-negative stochastic tensor $\boldsymbol{P}$, when $\boldsymbol{P}$ is perturbed. The following result provides a bound in terms of the higher-order norm-based ergodicity coefficients, and represents a tensor counterpart of Theorem 4.2.

Theorem 6.7**.**

Let $\boldsymbol{P}$ and its perturbation $\boldsymbol{P}^{\prime}$ be two stochastic tensors in $\mathbb{R}^{[3,m]}$. If $\mathcal{T}(\boldsymbol{P})<1$ then the stochastic solution of $x=\boldsymbol{P}xx$ is unique, and for any stochastic vector $x^{\prime}$ such that $x^{\prime}=\boldsymbol{P}^{\prime}x^{\prime}x^{\prime}$ it holds

[TABLE]

Proof 6.8**.**

Suppose first that both $\boldsymbol{P}$ and $\boldsymbol{P}^{\prime}$ are $S$-symmetric. By adding and subtracting $\boldsymbol{P}x^{\prime}x^{\prime}$ we have

[TABLE]

Rearranging terms we find $\|x-x^{\prime}\|_{1}(1-\mathcal{T}(\boldsymbol{P}))\leq\|\boldsymbol{P}-\boldsymbol{P}^{\prime}\|_{1}$ and the claim follows.

*In the general case, define $\boldsymbol{Q}=\frac{1}{2}(\boldsymbol{P}+\boldsymbol{P}^{S})$ and $\boldsymbol{Q}^{\prime}=\frac{1}{2}(\boldsymbol{P}^{\prime}+\boldsymbol{P}^{\prime S})$ and repeat the previous arguments. Finally, note that $\mathcal{T}(\boldsymbol{Q})=\mathcal{T}(\boldsymbol{P})$ and $\|\boldsymbol{Q}-\boldsymbol{Q}^{\prime}\|_{1}\leq\|\boldsymbol{P}-\boldsymbol{P}^{\prime}\|_{1}$. *

6.2 Convergence of a class of vertex reinforced random walks

Vertex reinforced random walks are another important example of higher-order discrete-time stochastic process $\{X(t)\}_{t}$ on the state space $\{1,\dots,n\}$, where the state transitions at time $t$ depend on the whole history $X(0),\ldots,X(t-1)$ [5, 42]. Starting from an initial state $X(0)\in\{1,\ldots,n\}$, the process evolves according to the formulas

[TABLE]

where $\mathcal{F}_{t}$ is the $\sigma$-field generated by $X(1),\dots,X(t)$, and $\mathcal{M}$ is a map from $\SS$ to the set of stochastic $n\times n$ matrices. The vector $y_{t}$, which is called the occupation vector, is an auxiliary stochastic vector that is introduced in order to record the history of the process. Indeed, the $i$-th entry of $y_{t}$ is proportional to the number of times the process visited state $i$ up to the $t$-th time step, plus one. Now, let $x_{t}$ be the probability vector of $X(t)$, that is, the $n$-vector whose $i$-th entry is $\mathbb{P}(X(t)=i)$. Then, the process (22) can be equivalently described via the coupled equations

[TABLE]

When $\mathcal{M}$ is linear, there exists a stochastic tensor $\boldsymbol{P}$ such that $\mathcal{M}(v)_{ij}=\sum_{k}\boldsymbol{P}_{ijk}v_{k}$ and the corresponding stochastic process is the so-called spacey random walk, introduced in [6]. In this case, with minor notation changes with respect to the original version, the previous iteration can be recast as

[TABLE]

On the basis of key results by Benaïm [5], Benson, Gleich and Lim established the convergence of the spacey random walk in terms of the convergence of a certain ordinary differential equation to a stable equilibrium, and one auxiliary condition placed on $\boldsymbol{P}$ [6, Thm. 9]. However, only the convergence of the occupation vectors $\{y_{t}\}$ (which corresponds to the convergence in the Cesàro average sense of the random variables $X(t)$) can be derived from the results in [5, 6]. In fact, the second equation in (23) yields

[TABLE]

Hence, even if the sequence $\{y_{t}\}$ has a limit and the left hand side converges to zero, that does not imply the convergence of the sequence $\{x_{t}\}$.

In what follows, we consider the following generalization of (23),

[TABLE]

with $c_{t}\in[0,1]$ and we show in the next theorem that, if the higher-order ergodicity coefficients $\mathcal{T}_{L}(\boldsymbol{P})$ and $\mathcal{T}_{R}(\boldsymbol{P})$ are small enough, then the stochastic process (24) is globally convergent, provided that the sequence $\{c_{t}\}$ is not too small. This requirement on $\{c_{t}\}$ can be seen as a condition that avoids the process from freezing along the way on a limit point that is far away from the $Z$-eigenvector of $\boldsymbol{P}$. In fact, the possibility of such a behavior has been shown in [11] for a stochastic process closely related to (24).

Theorem 6.9**.**

Let the sequence $\{c_{t}\}$ in (24) be non-increasing and such that

[TABLE]

If $\mathcal{T}_{L}(\boldsymbol{P})+\mathcal{T}_{R}(\boldsymbol{P})<1$, then the vertex reinforced random walk (24) converges globally, i.e., for any starting points $x_{0},y_{0}\in\SS$ we have

[TABLE]

*where $x$ is the unique stochastic solution of $x=\boldsymbol{P}xx$. Moreover, if there exists a positive constant $\alpha$ such that $c_{t}\geq\alpha$ then the convergence is linear. *

Proof 6.10**.**

Firstly, note that, in the stated hypotheses, the vector $x$ exists and is unique owing to Theorem 6.1. Subtracting the identity $x=\boldsymbol{P}xx$ from (24) we obtain

[TABLE]

Let $\alpha_{t}=\|x_{t}-x\|_{1}$ and $\beta_{t}=\|y_{t}-x\|_{1}$. Using vector inequalities, we have

[TABLE]

Let $\gamma_{t}=\|(\alpha_{t},\beta_{t})^{T}\|_{\infty}=\max\{\alpha_{t},\beta_{t}\}$. For notational simplicity, let $\ell=\mathcal{T}_{L}(\boldsymbol{P})$, $r=\mathcal{T}_{R}(\boldsymbol{P})$, and define

[TABLE]

Hence, for $t=1,2\ldots$ we have

[TABLE]

Moreover, since $\gamma_{t+1}\leq\|A_{t}\|_{\infty}\gamma_{t}$ and $\|A_{t}\|_{\infty}=1$, we have $\gamma_{t+1}\leq\gamma_{t}$, that is, the sequence $\{\gamma_{t}\}$ is non-increasing. Now, for $t\geq 1$ consider the product $A_{t}A_{t-1}$. Simple computations show that

[TABLE]

In particular, if $\lim_{t\to\infty}c_{t}=0$ then there exists an integer $t_{*}$ such that for $t\geq t_{*}$ it holds $\|A_{2t}A_{2t-1}\|_{\infty}=1-c_{2t}(\ell+r)$. Consequently, we have

[TABLE]

where $C=\|A_{2t_{*}-2}\cdots A_{1}A_{0}\|_{\infty}\gamma_{0}$. In order to prove that $\lim_{t\to\infty}\gamma_{t}=0$ it is sufficient to discuss the limit

[TABLE]

which exists and is nonnegative since all factors belong to $(0,1)$. By a known result on the convergence of infinite products, see e.g., [30, p. 223], the preceding limit is positive if and only if the series

[TABLE]

is convergent. Hence, if (25) holds then $\lim_{t\to\infty}\gamma_{t}=0$ and we are done.

On the other hand, if $c_{t}\geq\alpha>0$ then there exists a number $s\in(0,1)$ such that $\|A_{t}A_{t-1}\|_{\infty}\leq s$. Hence,

[TABLE]

*and the last claim follows. *

Note that both the spacey random walk (23) and the second-order Markov chain (21) are particular cases of the stochastic processes (24), corresponding to the choices $c_{t}=\frac{1}{t+1}$ and $c_{t}=1$, respectively. Observe that both these choices satisfy the assumption (25). Thus, the convergence condition for the second-order Markov chain of Theorem 6.5 also follows as a consequence of Theorem 6.9. Moreover, we obtain the following convergence result for the spacey random walk which, to the best of our knowledge, is the first result that gives explicit conditions that guarantee the convergence of both the occupation vector and the density distribution for this stochastic process.

Corollary 6.11**.**

If $\mathcal{T}_{L}(\boldsymbol{P})+\mathcal{T}_{R}(\boldsymbol{P})<1$ then the spacey random walk (23) converges globally, i.e., for any starting points $x_{0},y_{0}\in\SS$ we have

[TABLE]

*where $x$ is the unique stochastic solution of $x=\boldsymbol{P}xx$. *

Proof 6.12**.**

*It suffices to observe that the coefficient sequence $\{c_{t}\}$ of the spacey random walk is a trailing sub-sequence of the harmonic sequence, hence the hypothesis (25) is fulfilled. *

7 Comparison with previous works

In this section we discuss how the newly proposed higher-order ergodicity coefficient $\mathcal{T}(\boldsymbol{P})$, based on the $1$-norm, compares with previous works. In particular, we compare it with the contraction ratios proposed by Gautier and Tudisco in [22], where the Hilbert metric is used to quantify the contractivity of multilinear operators, and with the coefficients introduced by Li and Ng in [37] in order to characterize the uniqueness of stationary distributions of stochastic tensors.

7.1 Higher-order Birkhoff coefficients

When $d$ is the Hilbert projective metric $d_{H}$ defined in (9), the ergodicity coefficient (8) is known as Birkhoff contraction ratio and the renowned Birkhoff–Hopf theorem provides an explicit formula for such coefficient when $f$ is a linear map. Recently, the Birkhoff–Hopf theorem has been extended to the case of multilinear mappings [22]. We review that theorem in the following, for the case of a bilinear map $f:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ described by a cubic tensor $\boldsymbol{P}$ as $f(x,y)=\boldsymbol{P}xy$.

Theorem 7.1**.**

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a nonnegative tensor, let

[TABLE]

and let $\kappa(\boldsymbol{P})=\tanh(\frac{1}{4}\log\triangle(\boldsymbol{P}))$. Then

[TABLE]

From Theorem 7.1 we immediately derive a formula for the higher-order Birkhoff ergodicity coefficient for stochastic tensors, and the corresponding analogous of Theorem 6.1. Precisely, we have the following result.

Corollary 7.2**.**

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a stochastic tensor and let

[TABLE]

where

[TABLE]

*If $\mathcal{T}_{H}(\boldsymbol{P})<1$ then there exists a unique $Z$-eigenvector $x\in\mathcal{S}_{1}$ such that $\boldsymbol{P}xx=x$ and the higher-order power iteration $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$ converges to $x$ for any starting point $x_{0}\in S_{1}$. *

Proof 7.3**.**

Consider the $S$-symmetric tensor $\boldsymbol{Q}=\frac{1}{2}(\boldsymbol{P}+\boldsymbol{P}^{S})$. Note that $\triangle(\boldsymbol{Q})=\widehat{\triangle}(\boldsymbol{P})$ and thus $\kappa(\boldsymbol{Q})=\frac{1}{2}\mathcal{T}_{H}(\boldsymbol{P})$. Therefore, using the identity $\boldsymbol{P}xx=\boldsymbol{Q}xx$, which holds for all $x\in\mathbb{R}^{n}$, the triangle inequality for $d_{H}$ and Theorem 7.1, we have

[TABLE]

*This shows that $x\mapsto\boldsymbol{P}xx$ is a contraction with respect to the Hilbert metric. As $(\mathcal{S}_{1},d_{H})$ is a complete metric space, the proof continues as that of Theorem 6.1. *

Note that, similarly to the $1$-norm case, $\mathcal{T}_{H}(\boldsymbol{P})=0$ if and only if $\boldsymbol{P}$ has rank one, that is, $\mathcal{T}_{H}$ is proper. However, while $\mathcal{T}_{H}(\boldsymbol{P})=2$ for any tensor $\boldsymbol{P}$ not of rank one and having at least one zero entry, $\mathcal{T}(\boldsymbol{P})$ can be smaller than one even for sparse tensors. For example, if $\boldsymbol{P}$ is the tensor

[TABLE]

then one easily verifies that $\mathcal{T}_{H}(\boldsymbol{P})=2$, while $\mathcal{T}(\boldsymbol{P})=1/2$.

The left panel of Figure 1 scatter plots these two coefficients computed on a set of ten thousand random stochastic $n\times n\times n$ tensors with size $n$ between $2$ and $10$. In the matrix case it is well known that, for any stochastic matrix $P$ it holds $\tau_{1}(P)\leq\tau_{H}(P)$, see [51, §3.4]. While the numerical comparison shown in Figure 1 suggests the inequality $\mathcal{T}(\boldsymbol{P})\leq\mathcal{T}_{H}(\boldsymbol{P})$, an explicit comparison between the $1$-norm and the Birkhoff higher-order coefficients $\mathcal{T}(\boldsymbol{P})$ and $\mathcal{T}_{H}(\boldsymbol{P})$, for general tensors, is out of scope and is left open to future work.

7.2 Li and Ng’s coefficients

Given a stochastic tensor $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$, consider the following quantities introduced in [35, 37]:

[TABLE]

Li and Ng proved in [37] two conditions for the uniqueness of the stationary distribution and the convergence of the iteration $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$ in terms of the entries of $\boldsymbol{P}$, that we review in the following.

Theorem 7.4** ([37]).**

*Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a stochastic tensor. If $\gamma(\boldsymbol{P})>1$ then there exists an unique solution $x\in\SS$ of the equation $x=\boldsymbol{P}xx$. Moreover, the iteration $x_{t+1}=\boldsymbol{P}x_{t}x_{t}$ converges to $x$. *

As $\gamma(\boldsymbol{P})\geq 2\delta(\boldsymbol{P})$, the following consequence is immediate.

Corollary 7.5**.**

*Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a stochastic tensor. If $\delta(\boldsymbol{P})>1/2$ then all the claims in the preceding theorem are true. *

Moreover, we recall from [35, Thm. 4] the three-mode case of a perturbation bound for the stationary probability vector of a stochastic tensor of order $m>2$.

Theorem 7.6** ([35]).**

Let $\boldsymbol{P}$ and its perturbation $\boldsymbol{P}^{\prime}$ be two stochastic tensors in $\mathbb{R}^{[3,m]}$. If $\delta(\boldsymbol{P})>1/2$ then the stochastic solution of $x=\boldsymbol{P}xx$ is unique, and for any stochastic vector $x^{\prime}$ such that $x^{\prime}=\boldsymbol{P}^{\prime}x^{\prime}x^{\prime}$ it holds

[TABLE]

In the sequel, we aim to compare the above results with the ones we proved in the previous sections. First, we prove a special characterization of $\delta(\boldsymbol{P})$ in (27), which provides an explicit formula for $\delta(\boldsymbol{P})$ in terms of the entries of $\boldsymbol{P}$.

Lemma 7.7**.**

Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be a stochastic tensor. Then

[TABLE]

Proof 7.8**.**

First, note that for any zero-sum vector $y\in\mathbb{R}^{n}$ it holds \sum_{i}|y_{i}|=2\max\big{\{}\sum_{i\in\mathcal{I}}y_{i}:\mathcal{I}\subseteq\{1,\dots,n\}\big{\}}. Let $j_{1},j_{2},k_{1},k_{2}$ be fixed. Then we have

[TABLE]

Therefore

[TABLE]

*which coincides with (27), after rearranging terms. *

Using the characterization of $\delta(\boldsymbol{P})$ in the preceding lemma, the following theorem compares $\delta(\boldsymbol{P})$ and $\gamma(\boldsymbol{P})$ with the higher-order ergodic coefficient $\mathcal{T}(\boldsymbol{P})$:

Theorem 7.9**.**

*Let $\boldsymbol{P}\in\mathbb{R}^{[3,n]}$ be stochastic. Then $\mathcal{T}(\boldsymbol{P})\leq 2-2\delta(\boldsymbol{P})$. Moreover, if $\boldsymbol{P}=\boldsymbol{P}^{S}$ then $2-\gamma(\boldsymbol{P})\leq\mathcal{T}(\boldsymbol{P})$. *

Proof 7.10**.**

The formulas (11) and (14) can be rewritten as

[TABLE]

respectively. Using the preceding formulas and (20) it is immediate to obtain

[TABLE]

From Lemma 7.7 we conclude $1-\delta(\boldsymbol{P})\geq\frac{1}{2}\mathcal{T}(\boldsymbol{P})$ and this proves the first part of the claim. Furthermore, using the symmetry $\boldsymbol{P}=\boldsymbol{P}^{S}$, the formulas (26) and (19) simplify to

[TABLE]

*The inequality $2-\gamma(\boldsymbol{P})\leq\mathcal{T}(\boldsymbol{P})$ follows, and the proof is complete. *

We conclude with several important remarks that we obtain as a consequence of the preceding results.

First, notice that the requirement $\delta(\boldsymbol{P})>\frac{1}{2}$ appearing in Theorem 7.6 is stronger than the one of Theorem 6.7, namely, if $\delta(\boldsymbol{P})>\frac{1}{2}$ holds then $\mathcal{T}(\boldsymbol{P})<1$ must hold as well. Moreover, $2\delta(\boldsymbol{P})-1\geq 1-\mathcal{T}(\boldsymbol{P})$. Thus the right hand side of Theorem 7.6 is larger than the one of Theorem 6.7. This shows that Theorem 6.7 is an improvement over Theorem 7.6.

On the other hand, the condition $\gamma(\boldsymbol{P})>1$ is weaker than $\mathcal{T}(\boldsymbol{P})<1$. Hence, the hypothesis in Theorem 7.4 ensuring uniqueness of the solution of $x=\boldsymbol{P}xx$ and convergence of the higher-order power method can be more general than the one in Theorem 6.1, at least when $\boldsymbol{P}=\boldsymbol{P}^{S}$. Additionally, it is important to point out that the inequality $\mathcal{T}(\boldsymbol{P})<1$ can be checked using $\mathcal{O}(n^{4})$ arithmetic operations, while the computation of $\gamma(\boldsymbol{P})$ requires the solution of a nontrivial combinatorial optimization problem, which is in general significantly more expensive.

The central and the rightmost panels of Figure 1 compare numerically, via scatter plots, the condition $2-2\delta(\boldsymbol{P})<1$ and the ergodicity conditions $\mathcal{T}(\boldsymbol{P})<1$ and $\mathcal{T}_{H}(\boldsymbol{P})<1$ obtained via the higher-order ergoditicity coefficients, on a test set of $10,000$ randomly generated tensors with varying size.

8 Examples

We conclude with a number of example applications of Theorem 6.1. The examples here below further demonstrate the usefulness of the newly introduced higher-order ergodicity coefficients in a variety of contexts.

8.1 Multilinear PageRank

Given a stochastic tensor $\boldsymbol{P}$, a scalar $0<\alpha<1$ and a probability vector $v\in\SS$, the multilinear PageRank is a solution of the equation

[TABLE]

This definition has been introduced by Gleich, Lim, and Yu [25] in analogy to the renowned Google’s PageRank vector, defined as the solution of $\alpha Px+(1-\alpha)v=x$ where $P$ is a stochastic transition probability matrix. Pursuing that analogy, the solution of (28) gives the stationary probability of a stochastic process that with probability $\alpha$ behaves like the second-order Markov chain (21) and with probability $1-\alpha$ teleports to a random state chosen according to the discrete density $v$.

A detailed analysis of the possibly multiple nonnegative solutions to (28) is provided by Meini and Poloni in [39]. They also discuss various first- and second-order iterative methods to compute a solution to (28). In particular, fixed-point type methods are often a choice of preference, due to their inexpensive iterations and simple implementation. Also, these types of methods can be easily extrapolated achieving fast converge rates, see [16]. However, in practice one is interested in values of $\alpha$ not too far from $1$ but, unlike the matrix case, requiring $\alpha<1$ is not enough to ensure the uniqueness of the multilinear PageRank nor the convergence of the fixed-point iterates. In the original paper [25], the condition $\alpha<1/2$ is proved to be sufficient to ensure both these properties (28). More recently, a tighter sufficient condition for the uniqueness of the multilinear PageRank has been proved by Li et al. [36], in terms of the following quantity,

[TABLE]

where $\sigma$ is any real vector. Precisely, Theorems 1 and 2 in [36] show that if there exists $\sigma\in\mathbb{R}^{n}$ such that $\alpha\,\theta(\boldsymbol{P},\sigma)<1$, then (28) has a unique nonnegative solution and the fixed-point iteration for (28) converges to such a solution.

Theorem 6.1 provides a new condition that improves the range of values of $\alpha$ for which we can guarantee both the uniqueness of a nonnegative solution of (28) and the convergence of the associated fixed point iteration, as shown by the following result.

Corollary 8.1**.**

*If $\alpha\mathcal{T}(\boldsymbol{P})<1$ then (28) has a unique solution $x\in\SS$. Moreover, the fixed point iteration $x_{t+1}=\alpha\boldsymbol{P}x_{t}x_{t}+(1-\alpha)v$ converges linearly to $x$, with a convergence rate of at least $\alpha\mathcal{T}(\boldsymbol{P})$. Finally, it holds $\|x-v\|_{1}\leq 2\alpha$. *

Proof 8.2**.**

Let $\boldsymbol{V}\in\mathbb{R}^{[3,n]}$ be the rank-one tensor $\boldsymbol{V}_{ijk}=v_{i}$. Since $\boldsymbol{V}xx=v$ for any $x\in\SS$, the equation (28) can rewritten as $x=\boldsymbol{P}_{\alpha}xx$ where

[TABLE]

By Theorem 6.1, the condition $\mathcal{T}(\boldsymbol{P}_{\alpha})<1$ guarantees uniqueness of the solution and convergence of the fixed point iteration. However, $\mathcal{T}(\boldsymbol{P}_{\alpha})\leq\alpha\mathcal{T}(\boldsymbol{P})+(1-\alpha)\mathcal{T}(\boldsymbol{V})=\alpha\mathcal{T}(\boldsymbol{P}),$ due to the fact that $\mathcal{T}(\boldsymbol{V})=0$, as noted in Remark 5.7.

Finally, note that the vector $v$ is characterized by the identity $v=\boldsymbol{V}vv=\boldsymbol{P}_{0}vv$. Hence, by considering $\boldsymbol{P}_{\alpha}$ as a perturbation of $\boldsymbol{V}=\boldsymbol{P}_{0}$, from Theorem 6.7 we get

[TABLE]

*since $\mathcal{T}(\boldsymbol{V})=0$, and the proof is complete. *

Note that the condition for the uniqueness given by Corollary 8.1 is always an improvement with respect to the one of [36]. In fact, using the formula (15) for $\mathcal{T}(\boldsymbol{P})$, for any $\sigma\in\mathbb{R}^{n}$ we have

[TABLE]

In order to illustrate how the various conditions differ in practice, we consider two small example tensors borrowed from [25]

[TABLE]

Figure 2 compares the range of values of $\alpha$ that guarantee uniqueness of the multilinear PageRank and convergence of the corresponding fixed-point iteration for the two tensors $\boldsymbol{P}_{1}$ and $\boldsymbol{P}_{2}$, according to the original Theorem 5.1 in [25], Corollary 7.2, Theorem 6.1 and Theorem 1 in [36]. For the latter result, we show the value of the quantities $\alpha\,\theta(\boldsymbol{P},\sigma_{k})$, $k=1,2,3$ obtained with the three choices of vectors

[TABLE]

as proposed in Corollaries 1 and 2 in the same paper. The interesting ranges are those where the corresponding graphs stay below the dashed line.

8.2 Triangle-based PageRank on networks

Random walks are an important tool for exploratory network analysis. For example, they are at the basis of widely used methods for local clustering, link prediction and network centrality. The typical random walk on a network is a Markov process where the probability to move from a node $i$ to a node $j$ is proportional to the number of outgoing edges leaving from $i$. This classical first-order process only takes into account pairwise node-node relationships. However, recent work has highlighted that many important network features arise by exploiting the interaction of larger groups of nodes acting together, see e.g., [3, 8].

In order to account for this type of second-order node interaction, we can consider a second-order stochastic process on the network, where the probability to move to a node $i$ depends on the number of triangles that point towards $i$. We show in this section how the higher-order ergodicity coefficients for stochastic tensors help dealing with triangle-based random walks.

Let $G=(V,E)$ be an undirected graph, with $V=\{1,\dots,n\}$, and consider the tensor

[TABLE]

where $\triangle(j,k)$ is the number of triangles that contain both nodes $j$ and $k$. This tensor is the triangle-based version of the transition matrix of the standard random walk in $G$,

[TABLE]

where $d(i)=\sum_{j}A_{ij}$ is the degree of node $i$. Clearly, $\boldsymbol{T}$ has many vanishing columns as in general two nodes $j,k\in V$ may not participate in any triangle in $G$. In that case, we set $\boldsymbol{T}_{ijk}=1/n$ for all $i=1,\ldots,n$. Similarly, we set $A_{ij}=1/n$ for all $i=1,\dots,n$ if $j$ is an isolated node in $G$ (i.e. if the $j$-th column of $A$ is zero).

Now, define the tensor $\boldsymbol{A}$ as $\boldsymbol{A}_{ijk}=A_{ij}$ and, for $\beta\in[0,1]$, consider the stochastic tensor

[TABLE]

This construction has been considered for example in [3, 25], within the multilinear PageRank equation (28), in order to combine the standard and the triangle-based random walks on real-world networks. The next result specializes Theorem 6.1 to the multilinear PageRank problem associated with the tensor in (32) and, additionally, provides a bound on the distance between the solution to that problem and the standard PageRank vector.

Corollary 8.3**.**

Let $\boldsymbol{P}$ be defined as in (32) and let $\gamma=\alpha(1+\beta)$. If $\gamma<1$ then (28) has a unique solution $x\in\SS$, and the fixed point iteration $x_{t+1}=\alpha\boldsymbol{P}x_{t}x_{t}+(1-\alpha)v$ converges linearly to $x$, with a convergence rate of at least $\gamma$.

Moreover, let $z\in\SS$ be the solution of the ordinary PageRank problem corresponding to the transition matrix $A$ and the teleportation vector $v$,

[TABLE]

Then,

[TABLE]

Proof 8.4**.**

Observe that, as $\mathcal{T}(\boldsymbol{A})=\tau_{1}(A)\leq 1$, we have the trivial upper bound $\mathcal{T}(\boldsymbol{P})\leq\beta+1$, hence the first part of the claim follows from Corollary 8.1. Moreover, simple passages allow us to recast $x$ as the stochastic solution of the equation $x=\boldsymbol{P}_{\alpha,\beta}xx$ where

[TABLE]

and $\boldsymbol{V}$ is the rank-one tensor $\boldsymbol{V}_{ijk}=v_{i}$. Analogously, the vector $z$ can also be considered as the stochastic solution of $z=(\alpha\boldsymbol{A}+(1-\alpha)\boldsymbol{V})zz$, that is $z=\boldsymbol{P}_{\alpha,0}zz$. We have $\mathcal{T}(\boldsymbol{P}_{\alpha,0})\leq\alpha$, hence from Theorem 6.7 we get

[TABLE]

*which completes the proof. *

Together with Corollary 8.1, the previous result shows that small values of $\alpha$ and $\beta$ produce a multilinear PageRank vector that does not differ sensibly from the ordinary PageRank vector. Note that the uniqueness and convergence condition $\gamma<1$ in Corollary 8.3 can be fulfilled by any value $\alpha<1/(1+\beta)$. This condition is evidently less restrictive than the better known inequality $\alpha<1/2$. This is one of several implications of Corollary 8.3. Below we consider an example real-wold network to further showcase the advantages of that corollary.

The socfb-Carnegie49 network is a Facebook graph considered for example in the study of the social structure of Facebook users [53], and available online on NetworkRepository [48]. The graph has $6637$ nodes and $249967$ undirected edges. The largest connected component consists of $6621$ vertices, and the triangle tensor $\boldsymbol{T}$ has $13860318$ nonzero entries. Figure 3 shows the results of a number of multilinear PageRank problems with coefficients $\alpha$ and $\beta$ varying in $[0,1]$. The equation $x=\alpha\boldsymbol{P}xx+(1-\alpha)v$, with $\boldsymbol{P}=\beta\boldsymbol{T}+(1-\beta)\boldsymbol{A}$ as in (32) and with uniform teleportation vector $v=\mathbb{1}/n$ has been solved via the fixed point iteration $x_{t+1}=\alpha\boldsymbol{P}x_{t}x_{t}+(1-\alpha)v$ endowed with the stopping criterion $\|x_{t+1}-x_{t}\|_{1}<10^{-8}$.

The leftmost panel in Figure 3 shows the distance $\|x-v\|_{1}$, whereas the central panel shows $\|x-z\|_{1}$ where $z$ is the usual PageRank vector, defined as in (33), with the same $\alpha$ value chosen for the multilinear version. The iteration number to convergence is shown in the rightmost panel.

While the overall behavior of $\|x-z\|_{1}$ reflects the estimate in Corollary 8.3, the panel on the left shows that $x$ approaches $v$ not only when $\alpha\approx 0$ but also when $\beta$ is large. This is due to the fact that the triangle tensor $\boldsymbol{T}$ has many zero columns. Consequently, the vast majority of the columns of $\boldsymbol{P}$ coincide with the uniform vector $v$ and, for an arbitrary vector $x\in\SS$, the product $\boldsymbol{P}xx$ is in general very close to $v$, in the sense that the $1$-norm of the vector $r(x)=\boldsymbol{P}xx-v$ is rather small. With this notation we obtain

[TABLE]

Broadly speaking, when $\boldsymbol{T}$ is very sparse $\|r(x)\|_{1}$ is usually negligible, and we can adopt the estimate $\|x-v\|_{1}\approx\alpha(1-\beta)\|Ax-v\|_{1}$. This approximation justifies the small error $\|x-v\|_{1}$ observed when $\beta\approx 1$.

In conclusion, the most informative results are obtained when both $\alpha$ and $\beta$ are neither too small nor too close to $1$. Extensive numerical experiments we performed on several real-world networks suggest the “reference” choice $\alpha,\beta\approx 0.6$. These values yield a good balance between first- and second-order information, fulfill the condition $\alpha(1+\beta)<1$ in Corollary 8.3 and ensure a fast convergence of the fixed point iteration. As an illustration, in Figure 4 we compare via scatter plots the reference solution for $\alpha=\beta=0.6$ with other multilinear PageRank vectors for different values of $\alpha$ and $\beta$ chosen as follows: in the leftmost panel we compare the reference solution against the standard PageRank vector with $\alpha=0.6$; in the central panel the solution for $\alpha=\beta=0.6$ is compared against the purely triangle-based case $\alpha=0.6$ and $\beta=1$; in the rightmost panel the vector corresponding to $\alpha=\beta=0.6$ is scatter plotted against the solution for the three choices $\alpha=\beta\in\{0.7,0.8,0.9\}$. The first two panels show the sensitivity of the solution with respect to different choices of $\beta$ highlighting, in particular, the importance of both edge- and triangle-based walks in the graph. The last panel on the right, instead, shows that the reference solution for $\alpha=\beta=0.6$ highly correlates with other numerical solutions obtained with larger values of the coefficients $\alpha=\beta$. This illustrates that larger choices of the coefficients $\alpha=\beta$ essentially do not alter the information on the nodes, but require a much larger iteration count.

8.3 Higher-order shifted power method and lazy random walk

Let $\boldsymbol{P}$ be symmetric, that is, $\boldsymbol{P}=\boldsymbol{P}^{\langle\pi\rangle}$ for every permutation $\pi$ of $\{1,2,3\}$. In [31], Kolda and Mayo analyzed the convergence of the “shifted symmetric higher-order power method”

[TABLE]

Their starting point is the optimization of the cubic form $f(x)=x^{T}(\boldsymbol{P}xx)$ over the sphere $x^{T}x=1$, whose stationary points are, for symmetric tensors, $Z$-eigenvectors of $\boldsymbol{P}$ related to the best symmetric rank-one approximation of $\boldsymbol{P}$. The coefficient $\alpha$ can be chosen positive or negative, in order to make the modified function $f(x)+\alpha x^{T}x$ convex or concave, respectively. Using fixed point theory, the authors of [31] prove that, given an appropriate shift $\alpha$ the iterates in (34) generically converge to some $Z$-eigenvector. The shifting technique has been considered also for tensors that are not symmetric. For example, it has been considered in the framework of the multilinear PageRank [25] or in the case of $\ell^{p}$-eigenvalue computation [24].

Let the coefficient $\beta(\boldsymbol{P})$ be defined as $\beta(\boldsymbol{P})=2\max_{\|x\|_{2}=1}\rho(\boldsymbol{P}x)$, where $\rho(\boldsymbol{P}x)$ denotes the spectral radius of the matrix $\boldsymbol{P}x$. One of the main results from [31] is that, if $\boldsymbol{P}$ is symmetric and $\alpha>\beta(\boldsymbol{P})$, then the method (34) converges to some stationary point of $f$, which is a $Z$-eigenvector of $\boldsymbol{P}$.

If $\boldsymbol{P}$ is stochastic (but not necessarily symmetric), then it is natural to replace the sphere $x^{T}x=1$ with the simplex $\SS$ and the vector $2$-norm with the $1$-norm. With these replacements, and other minor notation changes, the iteration (34) boils down to

[TABLE]

which, for an initial stochastic vector $x_{0}$, will remain in $\SS$ throughout. This iteration coincides with the higher-order power method $x_{t+1}=\boldsymbol{P}_{\sigma}x_{t}x_{t}$ for the “shifted tensor”

[TABLE]

where $\boldsymbol{E}$ is any tensor such that $\boldsymbol{E}xx=x$, for all $x\in\SS$. For example, $\boldsymbol{E}$ can be chosen as a convex combination of the left and right identities $\boldsymbol{E}^{L}$ and $\boldsymbol{E}^{R}$, defined in (2).

Note that $\boldsymbol{P}_{\sigma}$ is stochastic, for any choice of $\sigma\in(0,1)$ and thus the iteration (35) can be interpreted as a form of higher-order lazy random walk. In fact, recall that if $P\in\mathbb{R}^{n\times n}$ is a stochastic matrix, then the Markov chain associated with $\sigma P+(1-\sigma)I$ is called lazy random walk, as it describes a walker that, with probability $\sigma$ performs a transition according to $P$, and remains in its current state otherwise. Hence, we can use Theorem 6.1 to provide a condition on $\sigma$, in terms of the entries of $\boldsymbol{P}$, that guarantees global convergence of the shifted power method (35).

Even though the higher-order ergodicity coefficient $\mathcal{T}(\boldsymbol{P})$ of the original tensor may be larger than one, suitable values of $\sigma$ can ensure that Theorem 6.1 holds for $\boldsymbol{P}_{\sigma}$. In fact, it is interesting to note that the function $\sigma\mapsto\mathcal{T}(\boldsymbol{P}_{\sigma})$ is continuous, piecewise linear and convex, with $\mathcal{T}(\boldsymbol{P}_{0})=\mathcal{T}(\boldsymbol{E})=1$. As $x=\boldsymbol{P}_{\sigma}xx$ if and only if $x=\boldsymbol{P}xx$, we deduce that

Corollary 8.5**.**

*If $\boldsymbol{P}$ is stochastic and $\mathcal{T}(\boldsymbol{P}_{\sigma})<1$ for some $\sigma\geq 0$, then $\boldsymbol{P}$ has a unique positive $Z$-eigenvector $x\in\SS$ and the method (34) converges to $x$, for any starting point $x_{0}\in\SS$, with a convergence rate of at least $\mathcal{T}(\boldsymbol{P}_{\sigma})$. *

Proof 8.6**.**

*The claim follows straightforwardly from Theorem 6.1 applied to $\boldsymbol{P}_{\sigma}$. *

In Figure 5 we show the value of $\mathcal{T}(\boldsymbol{P}_{\sigma})$ as a function of $\sigma$, for the two example tensors (30) and for the choice $\boldsymbol{E}=\frac{1}{2}(\boldsymbol{E}^{L}+\boldsymbol{E}^{R})$. Notice that for both the examples shown there exists an optimal $\sigma_{*}$ such that $\min_{\sigma}\mathcal{T}(\boldsymbol{P}_{\sigma})=\mathcal{T}(\boldsymbol{P}_{\sigma_{*}})<1$. Thus, although the higher-order ergodicity coefficient $\mathcal{T}(\boldsymbol{P})$ of the original tensor is larger than one, by Corollary 8.5 there exists a unique positive $x\in\SS$ such that $x=\boldsymbol{P}xx$ and we can compute it with a method that converges as $\|x_{t+1}-x\|_{1}\leq\mathcal{T}(\boldsymbol{P}_{\sigma^{*}})^{t}\|x_{0}-x\|_{1}$, for an arbitrary $x_{0}\in\SS$.

9 Conclusions

This work adds to the long and continuing history of applications of tensor methods to data science by providing a novel analysis of the long-term behaviour of higher-order stochastic processes governed by stochastic tensors. These types of processes are used in a large number of network science and data mining applications due to their ability to improve the underlining models and offer additional valuable insights. Even though stationary distributions of these processes are often required, fundamental mathematical questions such as the uniqueness of the distribution and the convergence behavior of the stochastic process remain unanswered. In fact, this is a relatively newly born and actively growing research area, with many open questions.

Following a natural extension of the widely used ergodicity coefficients for Markov chains, we have introduced a new family of higher-order ergodicity coefficients for higher-order processes that provides new and easily computable conditions to ensure existence, uniqueness and convergence towards the corresponding stationary distribution. The proposed analysis adds to previous work on uniqueness of $Z$-eigenvectors of stochastic tensors [12, 35, 37] and non-negative tensors in general [14, 21, 23, 24] by providing new conditions that are either less restrictive or are computationally easier to verify, or both.

Acknowledgments

The main results of this work have been developed during a visiting period that F.T. has spent at the Department of Mathematics, Computer Science and Physics of the University of Udine, Italy. He would like to thank the department and D.F. for the warm hospitality he received during that period.