Perron-based algorithms for the multilinear pagerank

TL;DR

This paper introduces new Perron-based algorithms for solving the multilinear PageRank problem, improving reliability and ability to handle more complex instances by leveraging quadratic vector equation theory and homotopy continuation.

Contribution

It develops novel fixed-point algorithms and a homotopy strategy for multilinear PageRank, addressing solution existence and computational challenges.

Findings

Proves existence of a minimal solution distinct from the stochastic one.

Develops more reliable algorithms capable of solving larger problem sets.

Demonstrates improved performance over existing methods.

Abstract

We consider the multilinear pagerank problem studied in [Gleich, Lim and Yu, Multilinear Pagerank, 2015], which is a system of quadratic equations with stochasticity and nonnegativity constraints. We use the theory of quadratic vector equations to prove several properties of its solutions and suggest new numerical algorithms. In particular, we prove the existence of a certain minimal solution, which does not always coincide with the stochastic one that is required by the problem. We use an interpretation of the solution as a Perron eigenvector to devise new fixed-point algorithms for its computation, and pair them with a homotopy continuation strategy. The resulting numerical method is more reliable than the existing alternatives, being able to solve a larger number of problems.

Tables

Table 1. Table 1: Iteration counts and times for the 29 benchmark tensors and α = 0.90 𝛼 0.90 \alpha=0.90

F	IO	N	O	ON	ON-T1	ON-IMP	ON-T2	ON-EXT	N-T1	N-IMP	N-T2	N-EXT
$184$	$45$	$4$	$19$	$7$	$30$	$23$	$23$	$47$	$8$	$6$	$23$	$13$
$18$	$40$	$5$	$22$	$9$	$40$	$32$	$25$	$64$	$11$	$9$	$25$	$21$
$18$	$40$	$5$	$14$	$8$	$40$	$32$	$25$	$64$	$11$	$9$	$25$	$21$
$67$	$28$	$4$	$17$	$7$	$23$	$23$	$23$	$40$	$6$	$6$	$23$	$11$
$54$	$83$	$6$	$29$	$10$	$62$	$30$	$26$	$94$	$19$	$13$	$26$	$30$
$227$	$169$	$7$	$59$	$9$	$55$	$32$	$27$	$82$	$19$	$10$	$27$	$28$
$680$	$142$	$5$	$63$	$8$	$39$	$25$	$25$	$57$	$11$	$7$	$25$	$17$
$117$	$153$	$6$	$63$	$8$	$32$	$25$	$25$	$51$	$9$	$7$	$25$	$15$
$512$	$166$	$5$	$63$	$8$	$37$	$25$	$25$	$51$	$11$	$7$	$25$	$15$
$407$	$160$	$6$	$95$	$9$	$47$	$25$	$26$	$73$	$15$	$7$	$26$	$23$
$654$	$140$	$5$	$97$	$8$	$39$	$24$	$25$	$58$	$12$	$6$	$25$	$17$
$186$	$141$	$6$	$51$	$9$	$47$	$31$	$25$	$72$	$15$	$8$	$25$	$22$
NaN	$118$	$6$	$73$	$8$	$46$	$31$	$25$	$64$	$15$	$9$	$25$	$21$
$248$	$134$	$6$	$58$	$8$	$53$	$32$	$25$	$79$	$17$	$11$	$25$	$27$
$2,259$	$151$	$5$	$79$	$8$	$31$	$24$	$24$	$49$	$10$	$7$	$24$	$15$
$171$	$186$	$5$	$92$	$7$	$31$	$23$	$23$	$42$	$9$	$6$	$23$	$12$
$114$	$155$	$6$	$43$	$8$	$32$	$31$	$25$	$51$	$11$	$10$	$25$	$17$
$240$	$175$	$6$	$68$	$8$	$47$	$25$	$26$	$72$	$16$	$8$	$26$	$24$
$204$	$124$	$6$	$52$	$8$	$47$	$32$	$26$	$66$	$14$	$10$	$26$	$22$
$594$	$140$	$5$	$57$	$8$	$31$	$24$	$25$	$50$	$9$	$6$	$25$	$14$
$227$	$167$	$5$	$79$	$8$	$47$	$25$	$26$	$66$	$14$	$7$	$26$	$21$
$1,514$	$207$	$8$	$103$	$7$	$31$	$25$	$24$	$50$	$9$	$8$	$24$	$16$
$119$	$127$	$8$	$55$	$10$	$40$	$27$	$28$	$59$	$14$	$10$	$28$	$20$
$164$	$109$	$5$	$51$	$7$	$40$	$24$	$25$	$58$	$12$	$7$	$25$	$17$
$170$	$225$	$7$	$46$	$8$	$48$	$33$	$25$	$65$	$15$	$11$	$25$	$23$
$91$	$110$	$6$	$40$	$9$	$48$	$26$	$26$	$71$	$14$	$8$	$26$	$22$
$69$	$99$	$5$	$38$	$7$	$31$	$24$	$24$	$43$	$8$	$6$	$24$	$12$
$112$	$131$	$6$	$58$	$8$	$38$	$31$	$25$	$65$	$13$	$10$	$25$	$22$
$70$	$106$	$6$	$38$	$7$	$30$	$24$	$24$	$48$	$10$	$8$	$24$	$16$

Table 2. Table 2: Iteration counts and times for the 29 benchmark tensors and α = 0.99 𝛼 0.99 \alpha=0.99

F	IO	N	O	ON	ON-T1	ON-IMP	ON-T2	ON-EXT	N-T1	N-IMP	N-T2	N-EXT
NaN	$42$	$5$	$21$	$6$	$35$	$22$	$22$	$46$	$10$	$6$	$22$	$13$
$23$	$48$	$6$	$26$	$7$	$50$	$36$	$24$	$80$	$17$	$13$	$24$	$29$
$23$	$48$	$6$	$10$	$7$	$50$	$36$	$24$	$80$	$17$	$13$	$24$	$29$
$114$	$24$	$4$	$19$	$6$	$28$	$22$	$22$	$39$	$8$	$6$	$22$	$11$
$22$	$41$	$17$	$15$	$7$	$74$	$26$	$24$	$110$	$23$	$29$	$24$	$37$
NaN	$639$	$7$	NaN	NaN	$73$	$39$	$25$	$109$	$26$	$14$	$25$	$38$
NaN	$234$	$5$	$369$	$8$	$44$	$24$	$24$	$63$	$13$	$7$	$24$	$20$
NaN	$1,221$	$6$	$1,845$	$9$	$46$	$26$	$26$	$69$	$16$	$9$	$26$	$24$
NaN	$532$	$5$	NaN	$7$	$44$	$24$	$25$	$63$	$15$	$7$	$25$	$20$
NaN	$339$	$6$	NaN	$7$	$59$	$30$	$25$	$85$	$20$	$9$	$25$	$29$
NaN	$228$	$5$	NaN	$7$	$45$	$23$	$25$	$64$	$15$	$7$	$25$	$20$
NaN	$357$	$6$	$362$	NaN	$62$	$31$	$26$	$90$	$21$	$9$	$26$	$30$
NaN	$156$	$6$	$7,509$	$6$	$52$	$31$	$25$	$77$	$18$	$10$	$25$	$27$
NaN	$329$	$6$	NaN	NaN	$65$	$32$	$25$	$92$	$22$	$11$	$25$	$33$
NaN	$550$	$6$	NaN	$8$	$38$	$23$	$25$	$56$	$13$	$7$	$25$	$19$
NaN	NaN	$5$	NaN	$7$	$31$	$23$	$23$	$47$	$10$	$7$	$23$	$14$
NaN	$1,173$	NaN	NaN	$71$	$99$	$32$	NaN	$87$	NaN	$11$	NaN	NaN
NaN	$541$	$6$	$1,595$	NaN	$60$	$31$	$25$	$89$	$21$	$10$	$25$	$32$
NaN	$273$	NaN	$248$	$8$	$53$	$37$	NaN	$78$	$17$	$12$	NaN	$27$
NaN	$230$	$7$	$196$	$7$	$37$	$24$	$24$	$56$	$12$	$7$	$24$	$17$
NaN	$533$	$6$	NaN	$9$	$59$	$31$	$26$	$90$	$19$	$10$	$26$	$29$
NaN	NaN	$7$	$205$	$6$	$31$	$32$	$24$	$50$	$10$	$16$	$24$	$17$
NaN	$435$	NaN	$504$	NaN	$59$	$35$	NaN	$89$	$23$	NaN	NaN	$31$
NaN	$147$	$5$	$137$	$7$	$46$	$29$	$25$	$64$	$15$	$8$	$25$	$20$
NaN	$835$	NaN	NaN	$44$	$76$	$59$	$41$	$121$	$40$	$36$	$41$	$44$
NaN	$654$	NaN	NaN	NaN	$72$	$33$	NaN	$103$	$24$	$12$	NaN	$37$
NaN	NaN	NaN	$8,025$	$1,326$	NaN	$903$	$652$	$2,878$	NaN	NaN	$652$	NaN
NaN	$334$	$6$	$307$	NaN	$57$	$37$	$24$	$91$	$21$	$13$	$24$	$33$
$2,005$	$2,100$	$9$	$909$	$10$	$34$	$33$	$27$	$52$	$14$	$14$	$27$	$20$

Table 3. Table 3: Iteration counts and times for the 29 benchmark tensors and α = 0.999 𝛼 0.999 \alpha=0.999

F	IO	N	O	ON	ON-T1	ON-IMP	ON-T2	ON-EXT	N-T1	N-IMP	N-T2	N-EXT
NaN	$41$	$5$	$21$	$6$	$35$	$22$	$22$	$46$	$10$	$6$	$22$	$13$
$29$	$51$	$6$	$28$	$7$	$50$	$36$	$24$	$80$	$18$	$13$	$24$	$29$
$29$	$51$	$6$	$9$	$6$	$50$	$36$	$24$	$80$	$18$	$13$	$24$	$29$
$122$	$24$	$4$	$19$	$6$	$28$	$22$	$22$	$39$	$8$	$6$	$22$	$11$
$18$	$38$	$336$	$13$	$8$	$74$	$26$	$23$	$111$	$23$	$71$	$23$	$37$
NaN	NaN	$7$	NaN	NaN	$73$	$39$	$25$	$110$	$26$	$14$	$25$	$39$
NaN	$251$	$5$	$723$	$8$	$45$	$24$	$25$	$63$	$14$	$7$	$25$	$20$
NaN	NaN	$6$	NaN	$9$	$46$	$26$	$26$	$69$	$16$	$9$	$26$	$24$
NaN	NaN	$5$	NaN	$7$	$44$	$24$	$25$	$63$	$15$	$7$	$25$	$20$
NaN	$383$	$6$	NaN	$7$	$59$	$31$	$25$	$85$	$20$	$9$	$25$	$29$
NaN	$243$	$5$	NaN	$7$	$45$	$23$	$25$	$64$	$15$	$7$	$25$	$20$
NaN	$442$	$6$	$824$	NaN	$63$	$31$	$26$	$95$	$22$	$9$	$26$	$32$
NaN	$164$	$6$	NaN	$6$	$52$	$31$	$25$	$77$	$18$	$10$	$25$	$27$
NaN	$396$	$6$	NaN	NaN	$66$	$36$	$26$	$97$	$23$	$12$	$26$	$35$
NaN	$812$	$6$	NaN	$8$	$38$	$23$	$25$	$60$	$13$	$7$	$25$	$21$
NaN	NaN	$5$	NaN	$7$	$31$	$23$	$23$	$47$	$10$	$7$	$23$	$14$
NaN	NaN	NaN	NaN	NaN	$100$	$33$	NaN	$92$	NaN	$12$	NaN	NaN
NaN	$848$	$6$	NaN	NaN	$65$	$31$	$25$	$93$	$23$	$10$	$25$	$33$
NaN	$323$	NaN	$388$	$8$	$54$	$37$	NaN	$78$	$18$	$12$	NaN	$27$
NaN	$253$	$7$	$258$	$7$	$37$	$24$	$24$	$56$	$12$	$7$	$24$	$17$
NaN	$857$	$6$	NaN	$9$	$59$	$31$	$26$	$90$	$19$	$10$	$26$	$29$
NaN	$334$	$7$	$165$	$6$	$30$	$32$	$24$	$50$	$10$	$16$	$24$	$17$
NaN	$577$	NaN	$1,185$	NaN	$59$	$35$	NaN	$94$	$23$	NaN	NaN	$33$
NaN	$150$	$5$	$160$	$7$	$46$	$29$	$25$	$69$	$15$	$8$	$25$	$22$
NaN	NaN	NaN	NaN	NaN	$81$	$60$	NaN	$126$	$42$	$37$	NaN	$46$
NaN	NaN	NaN	NaN	NaN	$72$	$38$	NaN	$108$	$25$	$14$	NaN	$39$
$348$	$358$	$50$	$276$	$192$	$108$	$64$	$53$	$71$	$232$	$54$	$53$	$131$
NaN	$406$	$7$	$522$	NaN	$57$	$38$	$24$	$91$	$21$	$13$	$24$	$33$
$205$	$224$	$12$	$108$	NaN	$48$	NaN	$30$	$56$	$31$	$17$	$30$	NaN

Equations

$$1_n^{\top }\mathbf{v}=1,1_n^{\top }R=1_{n^2}^{\top },$$

$$\mathbf{x}=\alpha R(\mathbf{x}\otimes \mathbf{x})+(1-\alpha )\mathbf{v}.$$

$$G(\mathbf{x}):=\alpha R(\mathbf{x}\otimes \mathbf{x})+(1-\alpha )\mathbf{v},$$

$$G_{\mathbf{x}}^{\prime }:=\alpha R(\mathbf{x}\otimes I_n)+\alpha R(I_n\otimes \mathbf{x}).$$

$${\mathbf{x}}_{k+1}=G({\mathbf{x}}_k),k=0,1,\dots ,$$

$$1^{\top }G(\mathbf{x})$$

$$=\alpha 1^{\top }(\mathbf{x}\otimes \mathbf{x})+(1-\alpha )=\alpha (1^{\top }\mathbf{x}{)}^2+(1-\alpha ).\text{qed}$$

$$(I-G_{{\mathbf{x}}_k}^{\prime }){\mathbf{x}}_{k+1}=(1-\alpha )\mathbf{v}-\alpha R({\mathbf{x}}_k\otimes {\mathbf{x}}_k),k=0,1,\dots .$$

$$\mathbf{y}=\alpha R(\mathbf{y}\otimes \mathbf{y})+\alpha R(\mathbf{y}\otimes \mathbf{m})+\alpha R(\mathbf{m}\otimes \mathbf{y})=\alpha R(\mathbf{y}\otimes \mathbf{y})+G_{\mathbf{m}}^{\prime }\mathbf{y}=(\alpha R(\mathbf{y}\otimes I_n)+G_{\mathbf{m}}^{\prime })\mathbf{y}.$$

$$1^{\top }\mathbf{y}=1^{\top }(\mathbf{s}-\mathbf{m})=1-\frac{1-\alpha }{\alpha }=\frac{2\alpha -1}{\alpha }.$$

$$1^{\top }{P}_{\mathbf{y}}$$

$$=\alpha 1^{\top }(\mathbf{y}\otimes I)+\alpha 1^{\top }(I\otimes \mathbf{m})+\alpha 1^{\top }(\mathbf{m}\otimes I)$$

$$=\alpha \frac{2\alpha -1}{\alpha }{1}^{\top }+\alpha \frac{1-\alpha }{\alpha }{1}^{\top }+\alpha \frac{1-\alpha }{\alpha }{1}^{\top }=1^{\top },$$

$$\mathbf{y}=\frac{2\alpha -1}{\alpha }\mathcal{P}\mathcal{V}(P_{\mathbf{y}}),$$

$${\mathbf{y}}_{k+1}=\frac{2\alpha -1}{\alpha }\mathcal{P}\mathcal{V}(P_{{\mathbf{y}}_k}),$$

$${\mathbf{w}}^{\prime }(\mathbf{y})=\alpha \left((I-P_{\mathbf{y}}+\mathbf{w}(\mathbf{y})1^{\top }{)}^{-1}R(I\otimes \mathbf{w}(\mathbf{y}))-\mathbf{w}(\mathbf{y})1^{\top }\right).$$

$$\alpha R(\mathbf{z}\otimes I)\mathbf{v}(0)+P_{\mathbf{y}}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=\alpha (1^{\top }\mathbf{z})\mathbf{v}(0)+\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0),$$

$$(I-P_{\mathbf{y}})\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=\alpha \left(R(I\otimes \mathbf{v}(0))-\mathbf{v}(0)1^{\top }\right)\mathbf{z},$$

$$(I-P_{\mathbf{y}}+\mathbf{v}(0)1^{\top })\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=\left(R(I\otimes \mathbf{v}(0))-\mathbf{v}(0)1^{\top }\right)\mathbf{z}.$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=\alpha (I-P_{\mathbf{y}}+\mathbf{v}(0)1^{\top }{)}^{-1}\left(R(I\otimes \mathbf{v}(0))-\mathbf{v}(0)1^{\top }\right)\mathbf{z}.$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=\alpha \left((I-P_{\mathbf{y}}+\mathbf{v}(0)1^{\top }{)}^{-1}R(I\otimes \mathbf{v}(0))-\mathbf{v}(0)1^{\top }\right)\mathbf{z},$$

$$H_{\mathbf{y}}^{\prime }=I_n+(2\alpha -1)\mathbf{w}(\mathbf{y})1^{\top }-(2\alpha -1)(I-P_{\mathbf{y}}+\mathbf{w}(\mathbf{y})1^{\top }{)}^{-1}R(I\otimes \mathbf{w}(\mathbf{y}))$$

$${\mathbf{y}}_{k+1}={\mathbf{y}}_k-(H_{{\mathbf{y}}_k}^{\prime }{)}^{-1}H({\mathbf{y}}_k)$$

$${\mathbf{s}}_{\hat{\alpha }}=\mathbf{s}+{\mathbf{s}}^{(1)}(\hat{\alpha }-\alpha )+{\mathbf{s}}^{(2)}(\hat{\alpha }-\alpha {)}^2+O((\hat{\alpha }-\alpha {)}^3),$$

$$0=F_{\alpha }({\mathbf{s}}_{\alpha })={\mathbf{s}}_{\alpha }-\alpha R({\mathbf{s}}_{\alpha }\otimes {\mathbf{s}}_{\alpha })-(1-\alpha )\mathbf{v},$$

$$\frac{\partial F_{\alpha }}{\partial {\mathbf{s}}_{\alpha }}$$

$$\frac{\partial F_{\alpha }}{\partial \alpha }$$

$${\mathbf{s}}^{(1)}=\frac{d}{d\alpha }{\mathbf{s}}_{\alpha }$$

$$\frac{d}{d\alpha }\frac{\partial F_{\alpha }}{\partial {\mathbf{s}}_{\alpha }}$$

$$\frac{d}{d\alpha }\frac{\partial F_{\alpha }}{\partial \alpha }$$

$${\mathbf{s}}^{(2)}=\frac{d^2{\mathbf{s}}_{\alpha }}{d{\alpha }^2}$$

$$=-{\left(\frac{\partial F_{\alpha }}{\partial {\mathbf{s}}_{\alpha }}\right)}^{-1}\left(\frac{d}{d\alpha }\frac{\partial F_{\alpha }}{\partial {\mathbf{s}}_{\alpha }}\right){\mathbf{s}}^{(1)}-{\left(\frac{\partial F_{\alpha }}{\partial {\mathbf{s}}_{\alpha }}\right)}^{-1}\frac{d}{d\alpha }\frac{\partial F_{\alpha }}{\partial \alpha }$$

$$=(I-G_{\alpha ,{\mathbf{s}}_{\alpha }}^{\prime }{)}^{-1}\left(R\left({\mathbf{t}}_{\alpha }\otimes {\mathbf{s}}^{(1)}\right)+R\left({\mathbf{s}}^{(1)}\otimes {\mathbf{t}}_{\alpha }\right)\right).$$

$${\mathbf{s}}^{(1)}=(I-{\alpha }_{h+1}R({\mathbf{s}}_{{\alpha }_h}\otimes I)-{\alpha }_{h+1}R(I\otimes {\mathbf{s}}_{{\alpha }_h}){)}^{-1}R({\mathbf{s}}_{{\alpha }_h}\otimes {\mathbf{s}}_{{\alpha }_h}-\mathbf{v}),$$

$$\frac{{\mathbf{s}}_{{\alpha }_h}-{\mathbf{s}}_{{\alpha }_{h-1}}}{({\alpha }_h-{\alpha }_{h-1}{)}^2}\approx \frac{{\mathbf{s}}_{{\alpha }_{h+1}}-{\mathbf{s}}_{{\alpha }_h}}{({\alpha }_{h+1}-{\alpha }_h{)}^2}.$$

Full text

1 Introduction

Gleich, Lim and Yu [1] consider the following problem, arising as an approximate computation of the stationary measure of an order-2 Markov chain: given $\mathbf{v}\in\mathbb{R}^{n},R\in\mathbb{R}^{n\times n^{2}}$, $\alpha\in\mathbb{R}$ with $\mathbf{v}\geq 0,R\geq 0,\alpha\in(0,1)$ and

[TABLE]

solve for $\mathbf{x}$ the equation

[TABLE]

The solution of interest $\mathbf{s}$ is stochastic, i.e., $\mathbf{s}\geq 0$ and $\mathbf{1}_{n}^{\top}\mathbf{s}=1$. Here $\mathbf{1}$ denotes a column vector of all ones (with an optional subscript to specify its length), and inequalities between vectors and matrices are intended in the componentwise sense. In the paper [1], they prove some theoretical properties, consider several solution algorithms, and evaluate their performance. In the more recent paper [2], the authors improve some results concerning the uniqueness of the solution.

This problem originally appeared in [3], and is a variation of problems related to tensor eigenvalue problems and Perron–Frobenius theory for tensors; see, e.g., [4, 5, 6]. However, it also fits in the framework of quadratic vector equations derived from Markovian binary tree models introduced in [7] and later considered in [8, 9, 10]. Indeed, the paper [10] considers a more general problem, which is essentially (2) without the hypotheses (1). Hence, all of its results apply here, and they can be used in the context of multilinear PageRank. In particular, [10] considers the minimal nonnegative solution of (2) (in the componentwise sense), which is not necessarily stochastic as the one sought in [1].

In this paper, we use the theory of quadratic vector equations in [8, 9, 10] to better understand the behavior of the solutions of (2) and suggest new algorithms for computing the stochastic solution. More specifically, we show that if one considers the minimal nonnegative solution of (2) as well, the theoretical properties of (2) become clearer, even if one is only interested in stochastic solutions. Indeed we prove that there always exists a minimal nonnegative solution, which is the unique stochastic solution when $\alpha\leq 1/2$. When $\alpha>\frac{1}{2}$, the minimal nonnegative solution $\mathbf{m}$ is not stochastic and there is at least one stochastic solution $\mathbf{s}\geq\mathbf{m}$. Note that [1, Theorem 4.3] already proves that when $\alpha\leq\frac{1}{2}$ the stochastic solution is unique and [2, Theorem 1] slightly improves this bound; our results give a broader characterization. All this is in Section 2.

When $\alpha\leq 1/2$, as already pointed out in [1], computing the stochastic solution of (2) is easy. Indeed, this is also due to the fact that the stochastic solution is the minimal solution, and for instance the numerical methods proposed in [7, 8, 9] perform very well. The most difficult case is when $\alpha>1/2$, in particular when $\alpha\approx 1$. Since the minimal solution $\mathbf{m}$ of (2) can be easily computed, the idea is to compute and deflate it, with a similar strategy to the one developed in [8, 9], hence allowing us to compute stochastic solutions even when they are not minimal. The main tool in this approach is rearranging (2) to show that (after a change of variables) a solution $\mathbf{x}$ corresponds to the Perron eigenvector of a certain matrix that depends on $\mathbf{x}$ itself. This interpretation in terms of Perron vector allows to devise new algorithms based on fixed point iteration and on Newton’s method. Sections 3 and 4 describe this deflation technique and the algorithms based on the Perron vector computation.

Finally, we propose in Section 5 a continuation strategy based on a perturbative approach that allows one to solve the problem for values $\hat{\alpha}<\alpha$ in order to obtain better starting values for the more challenging cases when $\alpha\approx 1$.

We report several numerical experiments in Section 6, to show the effectiveness of these new techniques for the set of small-scale benchmark problems introduced in [1], and draw some conclusions in Section 7.

2 Properties of the nonnegative solutions

In this section, we show properties of the nonnegative solutions of the equation (2). In particular, we prove that there always exists a minimal nonnegative solution, which is stochastic when $\alpha\leq 1/2$. These properties can be derived by specializing the results of [10], which apply to more general vector equations defined by bilinear forms.

We introduce the map

[TABLE]

and its Jacobian

[TABLE]

We have the following result.

Lemma 1**.**

Consider the fixed-point iteration

[TABLE]

*started from $\mathbf{x}_{0}=0$. Then the sequence of vectors $\{\mathbf{x}_{k}\}$ is such that $0\leq\mathbf{x}_{k}\leq\mathbf{x}_{k+1}\leq\mathbf{1}$, there exists $\lim_{k\to\infty}\mathbf{x}_{k}=\mathbf{m}$ and $\mathbf{m}$ is the minimal nonnegative solution of (2), i.e., equation (2) has a (unique) solution $\mathbf{m}\geq 0$ such that $\mathbf{m}\leq\mathbf{x}$ for any other possible solution $\mathbf{x}\geq 0$. *

Proof.

The map $G(\mathbf{x})$ is weakly positive, i.e., $G(\mathbf{x})\geq 0$, $G(\mathbf{x})\neq 0$ whenever $\mathbf{x}\geq 0$, $\mathbf{x}\neq 0$. Moreover, if $0\leq\mathbf{x}\leq\mathbf{1}$ then $0\leq G(\mathbf{x})\leq\mathbf{1}$. Therefore Condition A1 of [10] is satisfied which, according to Theorem 4 of [10], implies that the sequence of vectors $\{\mathbf{x}_{k}\}$ is bounded and converges monotonically to a vector $\mathbf{m}$, which is the minimal nonnegative solution of (2). ∎

2.1 Sum of entries and criticality

In this specific problem, the hypotheses (1) enforce a stronger structure on the iterates of (3): the sum of the entries of $G(\mathbf{x})$ is a function of the sum of the entries of $\mathbf{x}$ only.

Lemma 2**.**

Let $g(u):=\alpha u^{2}+(1-\alpha)$. Then, $\mathbf{1}^{\top}G(\mathbf{x})=g(\mathbf{1}^{\top}\mathbf{x})$ for any $\mathbf{x}\in\mathbb{R}^{n}$.

Proof.

[TABLE]

This fact has important consequences for the sum of the entries of the solutions of (2).

Corollary 3**.**

For each solution $\mathbf{x}$ of (2), $\mathbf{1}^{\top}\mathbf{x}$ is one of the two solutions of the quadratic $u=g(u)$, i.e., $u=1$ or $u=\frac{1-\alpha}{\alpha}$.

Let $u$ be one of the solutions of $u=g(u)$ and define the level set $\ell_{u}=\{\mathbf{x}:\mathbf{1}^{\top}\mathbf{x}=u,\mathbf{x}\geq 0\}$. Since $\ell_{u}$ is convex and compact, and since $G(\mathbf{x})$ maps $\ell_{u}$ to itself by Lemma 2, then the Brouwer fixed-point theorem implies the following result.

Corollary 4**.**

There exists at least a solution $\mathbf{x}\geq 0$ to (2) with $\mathbf{1}^{\top}\mathbf{x}=1$ and a solution $\mathbf{x}\geq 0$ with $\mathbf{1}^{\top}\mathbf{x}=\frac{1-\alpha}{\alpha}$.

Hence we can have two different settings, for which we borrow the terminology from [10].

Subcritical case

$\alpha\leq\frac{1}{2}$, hence the minimal nonnegative solution $\mathbf{m}=\mathbf{s}$ is the unique stochastic solution.

Supercritical case

$\alpha>\frac{1}{2}$, hence the minimal nonnegative solution $\mathbf{m}$ satisfies $\mathbf{1}^{\top}\mathbf{m}=\frac{1-\alpha}{\alpha}<1$ and there is at least one stochastic solution $\mathbf{s}\geq\mathbf{m}$.

Note that [1, Theorem 4.3] already proves that when $\alpha\leq\frac{1}{2}$ the stochastic solution is unique; these results give a broader characterization.

The tools that we have introduced can already be used to determine the behavior of simple iterations such as (3).

Theorem 5**.**

Consider the fixed-point iteration (3), with a certain initial value $\mathbf{x}_{0}\geq 0$, for the problem (2) with $\alpha>\frac{1}{2}$. Define $z_{k}:=\mathbf{1}^{\top}\mathbf{x}_{k}$. Then,

• •

If $z_{0}\in(0,\frac{1-\alpha}{\alpha}]$, then $\lim_{k\to\infty}z_{k}=\frac{1-\alpha}{\alpha}$, and the iteration (3) can converge only to the minimal solution $\mathbf{m}$ (if it converges).

• •

If $z_{0}\in(\frac{1-\alpha}{\alpha},1]$, then $\lim_{k\to\infty}z_{k}=1$, hence the iteration (3) can converge only to a stochastic solution (if it converges).

• •

If $z_{0}\in(1,+\infty)$, then $\lim_{k\to\infty}z_{k}=+\infty$, hence the iteration (3) diverges.

Proof.

Thanks to Lemma 2, the quantity $z_{k}:=\mathbf{1}^{\top}\mathbf{x}_{k}$ evolves according to $z_{k+1}=g(z_{k})$. So the result follows by the theory of scalar fixed-point iterations, since this iteration converges to $\frac{1-\alpha}{\alpha}$ for $z_{0}\in(0,\frac{1-\alpha}{\alpha}]$, to $1$ for $z_{0}\in(\frac{1-\alpha}{\alpha},1]$, and diverges for $z_{0}\in(1,+\infty)$. ∎

An analogous result holds for the subcritical case.

The papers [8, 9, 10] describe several methods to compute the minimal solution $\mathbf{m}$. In particular, all the ones described in [10] exhibit monotonic convergence, that is, $0=\mathbf{x}_{0}\leq\mathbf{x}_{1}\leq\mathbf{x}_{2}\leq\cdots\leq\mathbf{x}_{k}\leq\cdots\leq\mathbf{m}$. Due to the uniqueness and the monotonic convergence properties, computing the minimal solution $\mathbf{m}$ is typically simple, fast, and free of numerical issues. Hence in the subcritical case the multilinear PageRank problem is easy to solve. The supercritical case is more problematic.

Among all available algorithms to compute the minimal solution $\mathbf{m}$, we recall Newton’s method, which is one of the most efficient ones. The Newton–Raphson method applied to the function $F(\mathbf{x})=\mathbf{x}-G(\mathbf{x})$ generates the sequence

[TABLE]

The following result holds [10, Theorem 13].

Lemma 6**.**

Suppose that $\mathbf{m}>0$, and that $G^{\prime}_{\mathbf{m}}$ is irreducible. Then, Newton’s method (4) starting from $\mathbf{x}_{0}=0$ is well defined and converges monotonically to $\mathbf{m}$ (i.e., $0=\mathbf{x}_{0}\leq\mathbf{x}_{1}\leq\mathbf{x}_{2}\leq\cdots\leq\mathbf{x}_{k}\leq\cdots\leq\mathbf{m}$).

Algorithm (1) shows a straightforward implementation of Newton’s method as described above.

Note that the theory in [10, Section 9] shows how one can predict where zeros appear in $\mathbf{m}$ and eliminate them reducing the problem to a smaller one. Indeed, in view of the probabilistic interpretation of multilinear PageRank, zero entries in $\mathbf{m}$ can appear only when the second-order Markov chain associated to $R$ is not irreducible. So we can assume in the following that $\mathbf{m}>0$, and that the nonnegative matrix $G^{\prime}_{\mathbf{m}}$ is irreducible. In particular, in this case $G^{\prime}_{\mathbf{m}}$ also substochastic (as proved in [10, Theorem 6]).

3 Deflating the minimal solution

From now on, we assume to be in the supercritical case, i.e., $\alpha>1/2$, and that $\mathbf{m}$ has been already computed and is explicitly available.

We wish adapt to this setting the deflation strategy introduced in [9]. Since all solutions $\mathbf{s}$ to (2) satisfy $\mathbf{s}\geq\mathbf{m}$, it makes sense to change variable to obtain an equation in $\mathbf{y}:=\mathbf{s}-\mathbf{m}\geq 0$. After a few manipulations, using bilinearity of $(\mathbf{m}+\mathbf{y})\otimes(\mathbf{m}+\mathbf{y})$ and the fact that $\mathbf{m}=\alpha R(\mathbf{m}\otimes\mathbf{m})+(1-\alpha)\mathbf{v}$, one gets

[TABLE]

Moreover,

[TABLE]

We set $P_{\mathbf{y}}:=\alpha R(\mathbf{y}\otimes I_{n})+G^{\prime}_{\mathbf{m}}$. Note that $P_{\mathbf{y}}\geq G^{\prime}_{\mathbf{m}}$ for each $\mathbf{y}\geq 0$, hence it is irreducible. In addition, if $\mathbf{y}$ is chosen such that $\mathbf{1}^{\top}\mathbf{y}=\frac{2\alpha-1}{\alpha}$ as in (6), then

[TABLE]

so $P_{\mathbf{y}}$ is a stochastic matrix.

Let us introduce the map $\mathcal{PV}(A)$ that gives the Perron vector $\mathbf{w}$ of an irreducible matrix $A\geq 0$, normalized such that $\mathbf{1}^{\top}\mathbf{w}=1$. Then, since $P_{\mathbf{y}}$ is irreducible and stochastic, (5) and (6) show that

[TABLE]

i.e., the sought vector $\mathbf{y}$ is the Perron vector of the matrix $P_{\mathbf{y}}$, multiplied by the constant $\frac{2\alpha-1}{\alpha}$.

4 Perron-based algorithms

Equation (8) suggests a new fixed-point iteration for computing $\mathbf{y}$, which is analogous to the one appearing in [9],

[TABLE]

starting from a given nonnegative vector $\mathbf{y}_{0}$ such that $\mathbf{1}^{\top}\mathbf{y}_{0}=\frac{2\alpha-1}{\alpha}$. This iteration is implemented in Algorithm 2

We may also apply Newton’s method to find a solution to (8), following [8]. To this end, we first compute the Jacobian of the function $\mathbf{w}(\mathbf{y}):=\mathcal{PV}(P_{\mathbf{y}})$.

Lemma 7**.**

Let $\mathbf{w}(\mathbf{y}):=\mathcal{PV}(P_{\mathbf{y}})$, with $\mathbf{y}\geq 0$ such that $\mathbf{1}^{\top}\mathbf{y}=\frac{2\alpha-1}{\alpha}$. Then, its Jacobian is given by

[TABLE]

Proof.

Let us compute its directional derivative along the direction $\mathbf{z}$. We set $\mathbf{y}(h)=\mathbf{y}+h\mathbf{z}$; hence, $\frac{\mathrm{d}}{\mathrm{d}h}P_{\mathbf{y}(h)}=\alpha R(\mathbf{z}\otimes I)$. Since $P_{\mathbf{y}}$ is irreducible, its Perron eigenvalue is a simple eigenvalue, and hence we can define locally smooth functions $\lambda(h)$ as the Perron eigenvalue of $P_{\mathbf{y}(h)}$ and $\mathbf{v}(h)=\mathcal{PV}(P_{\mathbf{y}(h)})$. A computation analogous to (7) shows that $\mathbf{1}^{\top}P_{\mathbf{y}(h)}=(1+h\alpha\mathbf{1}^{\top}\mathbf{z})\mathbf{1}^{\top}$, hence $\lambda(h)=1+h\alpha\mathbf{1}^{\top}\mathbf{z}$ and $\frac{\mathrm{d}}{\mathrm{d}h}\lambda(h)=\alpha(\mathbf{1}^{\top}\mathbf{z})$.

By differentiating $P_{\mathbf{y}(h)}\mathbf{v}(h)=\lambda(h)\mathbf{v}(h)$ and evaluating at $h=0$, we get

[TABLE]

or, rearranging terms,

[TABLE]

where $\mathbf{v}(0)=\mathbf{w}(\mathbf{y})$. Since we defined $\mathbf{v}(h)$ so that $\mathbf{1}^{\top}\mathbf{v}(h)=1$, we have $\mathbf{1}^{\top}\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}h}(0)=0$, and hence also

[TABLE]

The matrix in the left-hand side is nonsingular, since it can be obtained by replacing the eigenvalue [math] with $1$ in the eigendecomposition of the singular irreducible M-matrix $I-P_{\mathbf{y}}$. Thus we can write

[TABLE]

Since $(I-P_{\mathbf{y}}+\mathbf{v}(0)\mathbf{1}^{\top})\mathbf{v}(0)=\mathbf{v}(0)$, we can simplify this expression further to

[TABLE]

from which we get (10). ∎

From the above result, the Jacobian of the function $H(\mathbf{y})=\mathbf{y}-\frac{2\alpha-1}{\alpha}\mathcal{PV}(P_{\mathbf{y}})$ is given by

[TABLE]

and Newton’s method consists of generating the sequence of vectors

[TABLE]

The Perron–Newton method is described in Algorithm 3.

The standard theorems [11] on local convergence of Newton’s method imply the following result.

Theorem 8**.**

Suppose that the matrix $H^{\prime}_{\mathbf{s}-\mathbf{m}}$ is nonsingular. Then the Perron–Newton method is locally convergent to a stochastic solution $\mathbf{s}$ of (2).

Remark 9**.**

Since $\mathbf{1}^{\top}\left(\mathbf{w}(\mathbf{y})\mathbf{1}^{\top}-(I-P_{\mathbf{y}}+\mathbf{w}(\mathbf{y})\mathbf{1}^{\top})^{-1}R(I\otimes\mathbf{w}(\mathbf{y}))\right)=0$, the matrix $H^{\prime}_{\mathbf{y}}$ has an eigenvalue 1 with left eigenvector $\mathbf{1}^{\top}$, for each value of $\mathbf{y}$.

5 Continuation techniques

The above algorithms, as well as those in [1], are sufficient to solve most of the test problems that are explored in [1]. However, especially when $\alpha\approx 1$, the algorithms may converge very slowly or stagnate far away from the minimal solution. For this reason, we explore an additional technique based on a perturbation expansion of the solution as a function of $\alpha$, inspired loosely by the ideas of homotopy continuation [12], which is a well-known strategy to derive approximate solutions for parameter-dependent equations.

The main result is the following.

Lemma 10**.**

Let $\mathbf{s}$ be a stochastic solution of problem (2) (for a certain $\alpha>\frac{1}{2}$), and suppose that $I-G^{\prime}_{\mathbf{s}}$ is nonsingular. Then, there is a stochastic solution $\mathbf{s}_{\hat{\alpha}}$ to (2) with the parameter $\alpha$ replaced by $\hat{\alpha}$ such that

[TABLE]

with $\mathbf{s}^{(1)}=(I-G^{\prime}_{\alpha,\mathbf{s}_{\alpha}})^{-1}(R(\mathbf{s}_{\alpha}\otimes\mathbf{s}_{\alpha})-\mathbf{v})$, $\mathbf{s}^{(2)}=(I-G^{\prime}_{\alpha,\mathbf{s}_{\alpha}})^{-1}R\left(\mathbf{t}_{\alpha}\otimes\mathbf{s}^{(1)}+\mathbf{s}^{(1)}\otimes\mathbf{t}_{\alpha}\right)$, where $\mathbf{t}_{\alpha}:=2\mathbf{s}_{\alpha}+\alpha\mathbf{s}^{(1)}$.

Proof.

Let us make the dependence of the various quantities on the parameter $\alpha$ explicit, i.e., we write $\mathbf{s}_{\alpha}$ to denote a stochastic solution of (2) for a certain value of the parameter $\alpha$, and similarly $G_{\alpha},G^{\prime}_{\alpha,\mathbf{x}}$ and $F_{\alpha}$.

We apply the implicit function theorem [13, Theorem 9.28] to

[TABLE]

obtaining

[TABLE]

Differentiating (13), we get

[TABLE]

The function $\mathbf{s}_{\hat{\alpha}}$ obtained by the theorem must satisfy $\mathbf{1}^{\top}\mathbf{s}_{\hat{\alpha}}=1$ for each $\hat{\alpha}>\frac{1}{2}$: indeed, by Corollary 3 for a solution $\mathbf{s}_{\hat{\alpha}}$ of the equation $\mathbf{x}=G_{\hat{\alpha}}(\mathbf{x})$ it must hold either $\mathbf{1}^{\top}\mathbf{s}_{\hat{\alpha}}=1$ or $\mathbf{1}^{\top}\mathbf{s}_{\hat{\alpha}}=\frac{1-\hat{\alpha}}{\hat{\alpha}}<1$, and the continuous function $\mathbf{1}^{\top}\mathbf{s}_{\hat{\alpha}}$ cannot jump from $1$ to $\frac{1-\hat{\alpha}}{\hat{\alpha}}$ without taking any intermediate value.

The formula (12) now follows by Taylor expansion. ∎

This result suggests a strategy to improve convergence for $\alpha\approx 1$: we start solving the problem with a small value of $\alpha=\alpha_{0}$, e.g., $0.6$, then we solve it repeatedly for increasing values $\alpha_{0}<\alpha_{1}<\alpha_{2}<\dots<\alpha_{k}=\alpha$; at each step $h$ we use as a starting point for our iterative algorithms an approximation of $\mathbf{s}_{\alpha_{h+1}}$ constructed from $\mathbf{s}_{\alpha_{h}}$.

We tried several different strategies to construct this approximation.

T2

Second-order Taylor expansion $\mathbf{s}_{\alpha_{h+1}}\approx\mathbf{s}_{\alpha_{h}}+\frac{d\mathbf{s}_{\alpha}}{d\alpha}|_{\alpha_{h}}(\alpha_{h+1}-\alpha_{h})+\frac{1}{2}\frac{d^{2}\mathbf{s}_{\alpha}}{d\alpha^{2}}|_{\alpha_{h}}(\alpha_{h+1}-\alpha_{h})^{2}$;

T1

First-order Taylor expansion $\mathbf{s}_{\alpha_{h+1}}\approx\mathbf{s}_{\alpha_{h}}+\frac{d\mathbf{s}_{\alpha}}{d\alpha}|_{\alpha_{h}}(\alpha_{h+1}-\alpha_{h})$;

EXT

Linear extrapolation from the previous two values $\mathbf{s}_{\alpha_{h+1}}\approx\mathbf{s}_{\alpha_{h}}+\frac{\mathbf{s}_{\alpha_{h}}-\mathbf{s}_{\alpha_{h-1}}}{\alpha_{h}-\alpha_{h-1}}(\alpha_{h+1}-\alpha_{h})$;

IMP

An attempt at building a ‘semi-implicit method’ with only the values available: ideally, one would like to set $\mathbf{s}_{\alpha_{h+1}}\approx\mathbf{s}_{\alpha_{h}}+\frac{d\mathbf{s}_{\alpha}}{d\alpha}|_{\alpha_{h+1}}(\alpha_{h+1}-\alpha_{h})$ (which would correspond to expressing $\mathbf{s}_{\alpha_{h}}$ with a Taylor expansion around $\alpha_{h+1}$); however, we do not know $\frac{d\mathbf{s}_{\alpha}}{d\alpha}|_{\alpha_{h+1}}$. Instead, we take $\mathbf{s}_{\alpha_{h+1}}\approx\mathbf{s}_{\alpha_{h}}+\mathbf{s}^{(1)}(\alpha_{h+1}-\alpha_{h})$ with

[TABLE]

i.e., we use the new value of $\alpha$ and the old value of $\mathbf{s}$ in the expression for the first derivative.

The only remaining choice is designing an effective strategy to choose the next point $\alpha_{h+1}$. We adopt the following one. Note that for T1,EXT,IMP we expect $\mathbf{s}_{\alpha_{h+1}}-\mathbf{s}_{\alpha_{h}}=O((\alpha_{h+1}-\alpha_{h})^{2})$. Hence one can expect

[TABLE]

Thus we choose a ‘target’ value $\tau$ (e.g., $\tau=0.01$) for $\mathbf{s}_{\alpha_{h+1}}-\mathbf{s}_{\alpha_{h}}$, and then we choose $\alpha_{h+1}$ by solving the equation

[TABLE]

The resulting continuation algorithm is described in Algorithm 4. Note that we start from $\alpha_{0}=0.6$, to steer clear of the double solution for $\alpha=0.5$. Convergence for such a value of $\alpha_{0}$ is typically not problematic.

Note that for strategy T2 one would expect $\mathbf{s}_{\alpha_{h+1}}-\mathbf{s}_{\alpha_{h}}=O((\alpha_{h+1}-\alpha_{h})^{3})$ instead, but for simplicity (and to ensure a fair comparison between the methods) we use the same step size selection strategy.

6 Numerical experiments

We have implemented the described methods using Matlab, and compared them to some of the algorithms in [14] on the same benchmark set of small-size models ($n\in\{3,4,6\}$) used in [1]. We have used tolerance $\varepsilon=\sqrt{\mathbf{u}}$, where $\mathbf{u}$ is the machine precision, $\mathbf{x}_{0}=\mathbf{v}$ as an initial value unless differently specified, and $\tau=0.01$ for Algorithm 4. To compute Perron vectors, we have used the output of Matlab’s eig. For problems with larger $n$, different methods (eigs, inverse iteration, Gaussian elimination for kernel computation, etc.) can be considered [15].

The results, for various values of $\alpha$, are reported in Tables 1 to 3. Each of the 29 rows represents a different experiment in the benchmark set, and the columns stand for the different algorithms, according to the following legend.

F

The fixed-point iteration (3), from an initial value such that $\mathbf{1}^{\top}\mathbf{x}_{0}=1$ (and with renormalization to $\mathbf{1}^{\top}\mathbf{x}_{k}=1$ at each step).

IO

The inner-outer iteration method, as described in [1].

N

Newton’s method with normalization to $\mathbf{1}^{\top}\mathbf{x}_{k}=1$ at each step, as described in [1]. Note that this is not Algorithm 1, which would converge instead to the minimal solution $\mathbf{m}$.

P

The Perron method (Algorithm 2).

PN

The Perron–Newton method (Algorithm 3).

N-xxx, PN-xxx

N or PN methods combined with continuation with one of the four strategies described in Section 5.

For each experiment, we report the number of iterations required, and the CPU times in seconds (obtained on Matlab R2017a on a computer equipped with a 64-bit Intel core i5-6200U processor). The value NaN inside a table represents failure to converge after 10,000 iterations. For P and PN, the number of iterations is defined as the sum of the number of Newton iterations required to compute $\mathbf{m}$ with Algorithm 1 and the number of iterations in the main algorithm. For the continuation algorithms, we report the number of inner iterations, that is, the total number of iterations performed by PN or N (respectively), summing along all calls to the subroutine.

In addition to the tables, performance profiles [16, Section 22.4] are reported in Figures 1–3.

Note that neither the number of iterations nor the CPU time are particularly indicative of the true performance of the algorithms: indeed, the iterations in the different methods amount to different quantity of work, since some require merely linear system solutions and some require eigenvalue computations. Moreover, Matlab introduces overheads which are difficult to predict for several instructions (such as variable access and function calls). In order to make the time comparisons more fair, for methods IO and N we did not use the code in [14], but rather rewrote it without making use of Matlab’s object-oriented features, which made the original functions significantly slower. The only change introduced with respect to their implementations is that we changed the stopping criterion to $\lVert G(\mathbf{x})-\mathbf{x}\rVert_{1}\leq\sqrt{\mathbf{u}}$, so that the stopping criterion is the same one for all tested methods.

In any case, the performance comparison between the various methods should be taken with a pinch of salt.

We comment briefly on the results. Newton-based methods typically require a constant number of iterations to converge, but on some of the benchmarks they fail to converge or require a much larger number of iterations. From the point of view of reliability, combining the Perron–Newton algorithm with continuation strategies gives a definite advantage: the resulting methods are the only ones (among the ones considered) that can solve successfully all the problems in the original benchmark set in [1], which contained the experiments with all values of $\alpha$ up to $0.99$. Raising $\alpha$ to the more extreme value of $0.999$ reveals failure cases for these methods as well. (However, if one reduces the step-size selection parameter $\tau$ to $0.001$, methods PN-EXT and PN-T1 succeed on all examples also for $\alpha=0.999$.)

Our results confirm the findings of [1] that problem R6_3 (the third to last one) for $\alpha=0.99$ is the hardest problem of the benchmark set; Perron-based algorithms can solve it successfully in the end, but (like the algorithms in [1]) they stagnate for a large number of iterations around a point which is far away from the true solution. The analysis in [1, End of Section 6.3 and caption of Figure 9] also attributes this difficulty to the presence of a ‘pseudo-solution’ that is an attraction basin for a large set of starting points.

To understand the specific difficulties encountered by continuation methods on problem R6_3, it is instructive to consider the plot in Figure 4.

This figure, which has been generated with the help of the software Bertini [17], shows the first entry of each valid stochastic solution $\mathbf{s}$ to (2) (recall that there may be more than one solution), for varying values of $\alpha$. It turns out that the problem has a single solution for all values of $\alpha$ up to $\alpha\approx 0.975$, then two more solutions appear, initially coincident and then slowly separating one from the other. Two of these three solutions merge at about $\alpha\approx 0.9899$ and then disappear (i.e., they are replaced by a pair of complex conjugate solutions), so for $\alpha=0.99$ there is only one real stochastic solution remaining. Note that the points where the number of solutions changes are points in which the Jacobian $I-G^{\prime}_{\mathbf{x}}$ is singular, in accordance with the implicit function theorem.

Continuation methods ‘track’ the upper solution, and hence fail to provide a good starting point for $\alpha=0.99$, when it disappears suddenly and is replaced by a completely different solution. Moreover, for $\alpha=0.99$, the numerical behavior of these iterations is affected by the presence of a point with a very small residual. For an intuition, the reader may consider the behavior of Newton’s method (on a function of a single variable) in presence of a local minimum very close to [math], e.g., on $f(x)=x^{2}+\varepsilon$ for a very small $\varepsilon$.

7 Conclusions

We have used the theory of quadratic vector equations in [8, 9, 10] to attack the multilinear PageRank problem described in [1]. The tools that we have introduced improve the theoretical understanding of this equation and of its solutions. In particular, considering the minimal solution in addition to the stochastic ones allows one to use a broader array of algorithms, with computational advantage. The new algorithms achieve better results when $\alpha\approx 1$, which is the most interesting and computationally challenging case.

Possible future improvements to these algorithms include experimenting with further step-size selection strategies, and more importantly integrating them with inexact eigenvalue solvers that are more suitable to computation with large sparse matrices.