Information Retrieval Under Network Uncertainty: Robust Internet Ranking

1. Introduction and Literature Review

PageRank (PR) is an automated information retrieval algorithm that was designed by Sergey Brin and Larry Page (Brin and Page 1998, Page et al. 1998) during the rise of the World Wide Web in order to improve the quality of search engines existing at that time. The algorithm relies not only on keyword matching but also on the additional structure present in the hypertext, resulting, therefore, in a higher-quality ranking.

Nowadays, different adaptations of the PageRank algorithm are successfully applied in a huge variety of areas (Agarwal and Venkatesh 2002, Aladwani and Palvia 2002, Pant and Sheng 2015, Lappas et al. 2016), starting with straightforward applications such as scientometrics (Liao et al. 2019, Xu et al. 2011), where ranking algorithms are used for measuring the importance of scientific literature, and social media, to seemingly unrelated fields such as financial risk management, where ranking is employed for the measurement of systemic risk in the banking sector (Sadoghi 2015, Chen et al. 2016, Capponi et al. 2020). In our article, we present the results in the context of web ranking; however, our findings can be easily applied in a wide variety of other areas.

The initial idea of the PageRank algorithm (Brin and Page 1998, Page et al. 1998) is based on the interchange of ranks between web pages in line with probability theory rules: a web page gives a part of its rank (or score) only to those web pages that it cites, that is, to which it refers with direct links. The amount of the transferred score is proportional to the probability for a random user to click on the link. In the initial PageRank formulation (Page et al. 1998, Bryan and Leise 2006, Langville and Meyer 2006, Franceschet 2011), this probability is assumed to be an inverse of the total number of links outgoing from a page. This implies equally probable transitions for a random user from some page to the pages cited on it. Similarly, transitions to referred pages are assumed to be uniform for huge-dimensional web ranking in practice. In general, though, this might not be true, as the user is more likely to move to the page with a higher rank than to the page with a lower rank. Therefore, for a better ranking quality, the transition probabilities should be assumed to be rank-dependent. One step forward is made in scientometrics, where the transition probabilities between journals are not uniform but are proportional to the number of cited articles belonging to the same journal (Timonina 2009, Xu et al. 2011). In our article, we make a further step forward, constructing a robust reformulation of the PageRank problem allowing transition probabilities to vary in a given ambiguity set.

The PageRank problem can be formulated as the problem of finding the principal eigenvector of a nonnegative transition matrix, where each column consists of transition probabilities from the corresponding page to all other pages. The principal eigenvector is not necessarily unique for nonnegative stochastic matrices (Horn and Johnson 1990). To impose uniqueness, the original PageRank problem is replaced with finding the principal eigenvector of a modified matrix with only positive components in the works (Brin and Page 1998, Page et al. 1998, Franceschet 2011). This matrix, known as the Google matrix, is used in practice and assumes a possibility of random transitions between pages. In the pioneering paper of Brin and Page (1998, section 2.1.2), the following intuitive justification is given to the matrix: “We assume there is a ‘random surfer’ who is given a web page at random and keeps clicking on links, never ‘hitting back’ but eventually gets bored and starts on another random page. The probability that the random surfer visits a page is its PageRank. And, the damping factor is the probability at each page the ‘random surfer’ will get bored and request another random page.” In the work of Polyak and Timonina (2011), it is shown that the principal eigenvector may be significantly influenced by this difference in matrices, and, importantly, for specific network structures, it may happen that the eigenvector of the Google matrix stops distinguishing between high-ranked pages that build the core of search engine results: if a large amount of high-ranked pages obtains the same score, then the ranking becomes ineffective, making the difference between “high-ranked” and “low-ranked” pages to be almost binary. In our article, we show that the Google matrix provides an approximate solution for the robust case with uniform probabilities. However, it accounts for neither rank-dependent transition probabilities nor a rapid network growth.

Generally, the approach with the modified matrix works well in the presence of a small number of high-ranked pages, which was true when the initial algorithm (Brin and Page 1998) was developed but which currently needs reconsideration as the amount of information on the web continues to grow and the number of experienced users in web research is increasing. In the work of Polyak and Timonina (2011), the authors propose an l₁-regularization method that avoids low-ranked pages. This optimization problem can be solved simply via the coordinate-wise descent method analogous to the Gauss-Seidel algorithm for solving linear equations (Polyak 1987, Tibshirani 1996, Bryan and Leise 2006, Langville and Meyer 2006, Polyak and Timonina 2011), enhancing efficiency of numerical ranking and improving accuracy by differentiating between high-importance pages. However, this optimization problem has a drawback: it is not robust to variations in the internet structure, while the World Wide Web is, clearly, changing over time; the amount of information on the web is growing, new pages appear, and some websites (old or spam sites) disappear (Zhu et al. 2014, Lappas et al. 2016, De los Santos and Koulayev 2017, Negahban et al. 2017). This is particularly important due to the fact that PageRank computation takes such a significant amount of time (Google reevaluation of PageRank takes weeks; Nazin and Polyak 2009) and around 1%–3% of new web pages open during this time with new links influencing the transportation matrix.

In the work of Juditsky and Polyak (2012), the robust PageRank model is proposed for networks, whose links may vary, influencing the transition matrix of a fixed dimension. Thus, by solving this problem, one obtains robust ranking for the case when the connections change. As for the internet, (i) the number of pages grows, and (ii) one does not explicitly know how the structure changes. In this article, we study networks which vary not only in links but also in the number of nodes, and we therefore consider growing structures. This is one of the main contributions of our work: allowing for the network growth.

Considering the growth rate to be fixed for the network, we formulate the robust PageRank in terms of a nonconvex optimization problem and bound this formulation from above by a strictly convex but nonsmooth optimization problem. Given some conditions on the parameters suggesting high uncertainty in new pages, the resulting problem appears to be independent of the growth rate and, thus, is also robust to its variations. In middle- to large-dimensional cases, convex nonsmooth optimization problems can be solved numerically using available optimization techniques, including interior-point methods, mirror-descent algorithms (Nemirovski and Yudin 1983, Lan et al. 2012) and randomized techniques (Nazin and Polyak 2009, Ishii and Tempo 2010, Ishii et al. 2005). In huge-dimensional cases, one can employ randomized subgradient algorithms, which are specifically designed for sparse matrices by Nesterov (2014), or first-order extragradient methods, such as Mirror-Prox and dual interpolation algorithms (Nemirovski 2004, Nesterov 2007, Juditsky and Nemirovski 2012). Moreover, one can use less accurate but extremely fast numerical methods proposed by Juditsky and Polyak (2012), Lei (2014), and Polyak and Timonina (2011). In our article, we propose some techniques that allow us to smoothen the optimization problem. We propose a projected subgradient algorithm, avoiding all nonsmooth points, and, furthermore, we develop a coordinate-wise descent method with a near-optimal step size. We also address high-dimensional algorithms in our article.

In Section 2, we propose an l₂-norm robust reformulation of the PageRank model with a growing network. Further, we state a separability lemma that plays a crucial role for the robustness: one can conclude under which conditions on the uncertainty the estimated ranking is robust to the network growth. In Section 3, we assess the approximation quality and demonstrate cases for the optimality gap reduction, including the case with zero optimality gap. In Section 4, we propose a coordinate-wise descent method with near-optimal step size. In Sections 5–7, we analyze the convergence and assess the impact of the network growth on ranking using real-world data sets on online movie repositories and on journals on computational complexity.

2. Mathematical Model

Let N be the total number of pages in the internet. The rank of some page $i,\hspace{0.17em}\forall i=1,\dots ,N$, denoted by x_i, is the probability for a random user to be on this page. Let ${\mathcal{L}}_i$ be the set of all web pages that refer to page i. The probability that a random user sitting on a page $j\in {\mathcal{L}}_i$ clicks on the link to page i is equal to P_ij. The interchange of ranks between pages is formulated in line with probability theory rules as $x_i=\sum _{j\in {\mathcal{L}}_i}{P}_{ij}{x}_j,\forall i=1,\dots ,N$, where the probability P_ij is assumed to be an inverse of the total number of links outgoing from page j denoted by n_j; that is, $P_{ij}=\frac{1}{n_j}$.

Next, consider a column-stochastic nonnegative transportation (or transition) matrix $P={\left\{P_{ij}\right\}}_{i,j=1,\dots ,N}\in {\mathbb{ℛ}}^{N\times N}$, which satisfies $P_{ij}\ge 0,\hspace{0.17em}\forall i,j=1,\dots ,N$, and ${\sum }_{i=1}^N{P}_{ij}=1,\forall j=1,\dots ,N$ by definition. The well-known Perron-Frobenius theorem (Horn and Johnson 1990) states that there exists a dominant (or principal) eigenvector $\overline{x}\in {\mathrm{\Sigma }}_N$, where ${\mathrm{\Sigma }}_N=\{v\in {\mathbb{ℛ}}^N:\hspace{0.17em}{\sum }_{i=1}^N{v}_i=1,\hspace{0.17em}{v}_i\ge 0\}$ is the standard simplex on ${\mathbb{ℛ}}^N$. In general, matrix P is huge-dimensional and very sparse, as the number of links outgoing from every node is much smaller than the total number of nodes in the network: the average number of links outgoing from a web page in the internet is equal to 20 for the whole web (Juditsky and Polyak 2012), while the total number of pages in the internet is approximately 10⁹. The dominant vector $\overline{x}$ describes the stationary distribution on the network structure: each element of the vector $\overline{x}$ denotes the rank of the corresponding node. Components of $\overline{x}$ can also be interpreted as proportions of time that an average user spends on a particular website.

The dominant vector $\overline{x}$ is not robust and may be highly vulnerable to small changes in matrix P, which happen both due to the changes in links and in the number of nodes in the World Wide Web. According to the Web Server Survey (https://news.netcraft.com/archives/category/web-server-survey/), the number of domains being registered per minute corresponds to the 1%– 3% average growth rate of the internet per month. That is why in this article we consider the following changes of matrix P:

where P is the column-stochastic transportation matrix describing the current state of the network with N pages; ξ is the matrix describing variations in links of the initial network; ζ, ψ, and χ are matrices describing links to and from M new pages, which may appear in the future (ξ is of the size N × N, ψ is of the size M × N, ζ is of the size N × M, and χ is of the size M × M). In the internet, $M\approx 0.03N$ per month.

We say that matrices P and Q are nonnegative if ξ, ζ, ψ, and χ satisfy the following properties:

In addition, the following properties must hold for column-stochastic matrices Q and P, with a right dominant eigenvector corresponding to the eigenvalue 1:

Similar to the work of Juditsky and Polyak (2012), the function ${\max }_{Q\in \mathrm{\Xi }}\Vert Qx-x\Vert$, where $\Vert ·\Vert$ is some norm, can be seen as a measure of “goodness” of a vector x as a common dominant eigenvector of the family Ξ, where Ξ stands for the set of perturbed stochastic matrices of the form (1) under Conditions (2) and (3).

Further, let us denote by $x^T=(\begin{array}{c}{x}^{{(1)}^T},\hspace{0.17em}{x}^{{(2)}^T}\end{array})$ a feasible point that is a candidate for the common dominant eigenvector of the family Ξ. Let $x^{(1)}$ be of the size $N\times 1$, and let $x^{(2)}$ be of the size $M\times 1$. Hence, x is of the size $(N+M)\times 1$. Note that the vector x must belong to the standard simplex ${\mathrm{\Sigma }}_{N+M}=\{v\in {\mathbb{ℛ}}^{N+M}:{\sum }_{i=1}^{N+M}{v}_i=1,v_i\hspace{0.17em}\ge \hspace{0.17em}0\}$, which means $x_i^{(1)}\ge 0,\hspace{0.17em}\forall i=1,\dots ,N,\hspace{0.17em}{x}_j^{(2)}\hspace{0.17em}\ge \hspace{0.17em}0,\forall j=1,\dots ,M$, and ${\sum }_{i=1}^N{x}_i^{(1)}+{\sum }_{j=1}^M{x}_j^{(2)}=1$ (i.e., $\Vert x^{(1)}{\Vert }_1+\Vert x^{(2)}{\Vert }_1=1$).

We say that the vector $\widehat{x}$ is a robust solution of the eigenvector problem on Ξ if

where $\Vert ·\Vert$ is some norm (further, we consider the l₂-norm robust formulation, but the results can be extended to the case with l₁-norm).

The reasonable choice of the uncertainty set Ξ would impose some bounds on the column-wise norms of matrices ξ, ζ, ψ, and χ, meaning that the perturbation in links of the current and future states of the network would be bounded; that is, $\Vert {[\xi ]}_j\Vert \le {\mathit{\epsilon }}_j^{(\xi )},\hspace{0.17em}\Vert {[\psi ]}_j\Vert \le {\mathit{\epsilon }}_j^{(\psi )},\Vert {[\zeta ]}_j\Vert \le {\mathit{\epsilon }}_j^{(\zeta )}$, and $\Vert {[\chi ]}_j\Vert \le {\mathit{\epsilon }}_j^{(\chi )}$, where ${[·]}_j$ denotes the jth column of a matrix. Moreover, the total uncertainty budget for matrices ξ, ζ, ψ, and χ could be fixed (Juditsky and Polyak 2012): this would imply constraints on the overall possible perturbations of the transportation matrix P.

By solving the optimization Problem (4), one would make the rank vector $\widehat{x}$ robust to link or node perturbations influencing the transportation matrix P. Currently, the Google scores vector is being updated periodically without accounting for the 1%–3% average monthly growth rate (e.g., Nazin and Polyak 2009 and Internet Live Stats since 2014). By solving the optimization problem of the type (4), one could decrease the number of updates of the score vector and, therefore, reduce the underlying personnel and machinery costs. Furthermore, the resulting ranking would be useful for cutting-edge applications in operations research, operations management, and marketing (Palmer 2002, Dhar et al. 2014, Pant and Sheng 2015, Lappas et al. 2016, De los Santos and Koulayev 2017, Negahban et al. 2017) by providing robust evaluation of the page quality.

We say that the vector ${\widehat{x}}^{(l_2)}$ is the l₂-robust solution of the eigenvector problem on the set of perturbed matrices ${\mathrm{\Xi }}^{(l_2)}$ if

where we require column stochasticity of the matrix Q and impose constraints on jth column ${[·]}_j$ of perturbation matrices ξ, ζ, ψ, and χ. Also, we impose second-order “total uncertainty budget” constraints on the perturbation matrices bounding corresponding Frobenius norms, which we denote by $\Vert ·{\Vert }_F$.

The optimization Problem (5) accounts for uncertainties in links between current and future pages: these uncertainties are incorporated via matrices ξ, ψ, ζ, and χ, whose worst-case realizations lead to the robust PageRank. The optimization problem differs from the robust formulation proposed by Juditsky and Polyak (2012), who studied uncertainties implied by the matrix $P+\xi$, describing variations in links of existing pages only. This formulation implies a lower bound on our model.

Theorem 1

(Lower Bound). Given ${\mathit{\epsilon }}_j^{(\chi )}\ge 1\hspace{0.17em}\forall j,\hspace{0.17em}{\mathit{\epsilon }}^{(\chi )}\ge 1$, the fixed-size network model imposes a lower bound on the optimization Problem (5) in the following sense:

$$\underset{\begin{array}{l} \text{Conditions}\hspace{0.17em}(2)–(3)\\ \Vert {[\xi ]}_j{\Vert }_1\le {\mathit{\epsilon }}_j^{(\xi )}\hspace{0.17em}\forall j,\hspace{0.17em}\Vert \xi {\Vert }_F\le {\mathit{\epsilon }}^{(\xi )}\end{array}}{\max }\Vert (P+\xi )x^{(1)}-x^{(1)}{|}_2\le \varphi (x)==\underset{Q\in {\mathrm{\Xi }}^{(l_2)}}{\max }\Vert Qx-x{|}_2.$$(6)

Proof.

See Online Appendix A.1 for the proof. Note that we do not require ψ = 0, and, thus, we allow columns of $P+\xi$ to sum up to values less than or equal to 1 in the left-hand side of the inequality.

Importantly, the network described by Juditsky and Polyak (2012) imposes column-wise conditions ${\mathrm{𝟙}}_N^T{[\xi ]}_j=0$ but not the nonnegativity constraints $P_{ij}+{\xi }_{ij}\ge 0,\forall i,j$. In line with Online Appendix A.1, we can relax these constraints in both sides of Inequality (6). Now, denoting by s the total importance of current N pages, that is, ${\mathrm{𝟙}}_N^T{x}^{(1)}=s$ with $s\in [0,1]$, we obtain

$$\underset{\begin{array}{l} x^{(1)}\ge 0\\ {\mathrm{𝟙}}_N^T{x}^{(1)}=s\end{array}}{\min }\underset{\begin{array}{l} {\mathrm{𝟙}}_N^T{[\xi ]}_j=0\\ \Vert {[\xi ]}_j{\Vert }_1\le {\mathit{\epsilon }}_j^{(\xi )},\forall j,\Vert \xi {\Vert }_F\le {\mathit{\epsilon }}^{(\xi )}\end{array}}{\max }\Vert (P+\xi )x^{(1)}-x^{(1)}{|}_2\le \le \underset{\begin{array}{l} \hspace{0.17em} x^{(1)},x^{(2)}\ge 0\\ {\mathrm{𝟙}}_N^T{x}^{(1)}=s,{\mathrm{𝟙}}_M^T{x}^{(2)}=1-s\end{array}}{\min } \underset{Q\in {\mathrm{\Xi }}^{(l_2)}/\{P_{ij}+{\xi }_{ij}\ge 0,\hspace{0.17em}\forall i,j\}}{\max }\Vert Qx-x{|}_2,$$

for any fixed s. $\square$

As the optimization Problem (5) is not convex-concave, the function $\varphi (x)={\max }_{Q\in {\mathrm{\Xi }}^{(l_2)}}\Vert Qx-x{\Vert }_2$ cannot be computed efficiently. For the computable robust dominant eigenvector of the family of perturbed matrices ${\mathrm{\Xi }}^{(l_2)}$, we formulate the following theorem, denoting ${\mathit{\epsilon }}_1={\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}$ and ${\mathit{\epsilon }}_2={\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1$.

Theorem 2

(Upper Bound). The following upper bound holds:

$$\varphi (x)=\underset{Q\in {\mathrm{\Xi }}^{(l_2)}}{\max }\Vert Qx-x{\Vert }_2\le \Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert x^{(1)}{\Vert }_{(a)}++{\mathit{\epsilon }}_2\Vert x^{(2)}{\Vert }_{(b)},$$(7)

where $\Vert x^{(1)}{\Vert }_{(a)}={\min }_{\lambda +\mu =x^{(1)}}\left\{\Vert \lambda {\Vert }_2+{\sum }_{j=1}^N\frac{{\mathit{\epsilon }}_j^{(\xi )}+{\mathit{\epsilon }}_j^{(\psi )}}{{\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}}|{\mu }_j|\right\}$ and $\Vert x^{(2)}{\Vert }_{(b)}={\min }_{\lambda +\mu =x^{(2)}}\left\{\Vert \lambda {\Vert }_2+{\sum }_{j=1}^M\frac{{\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1}{{\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1}|{\mu }_j|\right\}$.

Proof.

See Online Appendix A.2 for the complete proof. As a consequence, the optimal value $\varphi ({\widehat{x}}^{(l_2)})$ of the nonconvex optimization Problem (5) can be bounded from above as

$$\varphi ({\widehat{x}}^{(l_2)})\le \tilde{\varphi }(\tilde{x})==\underset{x\in {\mathrm{\Sigma }}_{N+M}}{\min }\{\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert x^{(1)}{\Vert }_{(a)}+{\mathit{\epsilon }}_2\Vert x^{(2)}{\Vert }_{(b)}\}.$$(8)

The parameter ${\mathit{\epsilon }}_1={\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}$ describes the total uncertainty in links of current N pages, while the parameter ${\mathit{\epsilon }}_2={\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1$ describes the uncertainty in links of new M pages. Importantly, both parameters can be chosen independently from the growth M, which makes the resulting vector robust to variations in the number of new nodes. We refer to the vector $\tilde{x}$ as the (computable) robust dominant eigenvector of the corresponding family of perturbed matrices ${\mathrm{\Xi }}^{(l_2)}$ if $\tilde{x}\in {\mathrm{Argmin}}_{x\in {\mathrm{\Sigma }}_{N+M}}\{\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert x^{(1)}{\Vert }_{(a)}+{\mathit{\epsilon }}_2\Vert x^{(2)}{\Vert }_{(b)}\}$.

In general, this problem is not strictly convex, as the objective function reduces to a constant plus a term proportional to the l₁-norm for a fixed x and λ. Nevertheless, functions $\Vert x^{(1)}{\Vert }_{(a)}$ and $\Vert x^{(2)}{\Vert }_{(b)}$ are equal to strictly convex l₂-norms for large-enough uncertainty parameters ${\mathit{\epsilon }}_j^{(\xi )}+{\mathit{\epsilon }}_j^{(\psi )}$ and ${\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1$, respectively. When not equal, the objective function can be bounded from above by $\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_2\Vert x^{(2)}{\Vert }_2$, which is strictly convex as the sum of a strictly convex and a convex function and, thus, has a unique minimum.

Note that the formulation in (8) coincides with the upper bound proposed by Juditsky and Polyak (2012) if M = 0 (i.e., if there is no growth in the network). In general, for M > 0, the robust reformulation in (8) differs from the l₂-bound imposed by the fixed-size network. The exception is the trivial case with ψ = 0 and M > 0, that is, the transition matrix with no links from old to new web pages. This case implies that Q is reducible, meaning that a random user who starts from some new page (among $i=N+1,\dots ,N+M$) and keeps clicking on links eventually results on one of pages $i=1,\dots ,N$ as $\zeta \ne 0$ and is not able to return to the starting page by clicking on links as ψ = 0. In this setup, the difference between the fixed-size and the growing network models is fully imposed by the links from new pages, and the ranks $x^{(2)}$ are a priori lower than $x^{(1)}$. Furthermore, as one cannot say which of pages $i=N+1,\dots ,N+M$ have higher ranks and which have lower ones before the internet structure is updated, the ranks $x^{(2)}$ could be assumed to be zero without loss of generality, and the problem could be solved numerically using algorithms proposed by Juditsky and Polyak (2012).

In this article, we focus on parameters ${\mathit{\epsilon }}^{(\psi )}>0,\hspace{0.17em}{\mathit{\epsilon }}^{(\zeta )}>0,{\mathit{\epsilon }}^{(\chi )}>0$, and we test the robust setup with ${\mathit{\epsilon }}^{(\zeta )}=0$ in Section 8, where we analyze the influence of the uncertainty budget ${\mathit{\epsilon }}^{(\chi )}=\sqrt{M}$ on ranking of current pages. Clearly, positive uncertainty budgets provide a necessary but not sufficient condition for Q to be irreducible, and, thus, the uncertainty set ${\mathrm{\Xi }}^{(l_2)}$ still includes reducible, as well as periodic transition matrices, for which the eigenvector corresponding to the maximal eigenvalue is not unique. Importantly, as we conclude from Theorem 2, the worst-case realization of Q leads to a formulation that can be bounded from above by a convex optimization problem regardless of possible reducibility or periodicity of Q. Furthermore, it can be bounded by a strictly convex function $\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert x^{(1)}{\Vert }_2+{\mathit{\epsilon }}_2\Vert x^{(2)}{\Vert }_2$ that has a unique minimum.

Lastly, recall that $\Vert {[\zeta ]}_j{\Vert }_1+\Vert {[\chi ]}_j{\Vert }_1=1$ and $\Vert {[P+\xi ]}_j{\Vert }_1+\Vert {[\psi ]}_j{\Vert }_1=1$ due to column-stochasticity of matrix Q. This leads to the requirement ${\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}\ge 1$, as the optimization Problem (5) would be infeasible otherwise. Also, matrix $P+\xi$ needs to be nonnegative with none of its columns summing up to values larger than 1. At maximum, $\Vert {[\xi ]}_j{\Vert }_1=2$ and $\Vert {[\xi ]}_j{\Vert }_1+\Vert {[\psi ]}_j{\Vert }_1=2$, which happens when matrix ξ nullifies all elements of P redistributing the values among other rows. At minimum, $\Vert {[\xi ]}_j{\Vert }_1+\Vert {[\psi ]}_j{\Vert }_1=0$, because P is column-stochastic by itself.

2.1. Robustness and Ambiguity Set

The optimal solution of the optimization Problem (8) is robust to variations in links and to the network growth. Excluding the point ${\tilde{x}}^{(1)}=0,\hspace{0.17em}{\tilde{x}}^{(2)}>0$ from the feasible region in Problem (8), one would logically imply that the uncertainty about current pages is not larger than the uncertainty about future (not yet existent) pages. Instead of cutting this point directly out of the feasible region, we find the exact cases when it may become the optimal solution.

Lemma 1.

Let ${\tilde{y}}^{(1)}$ and ${\tilde{y}}^{(2)}$ be solutions of the optimization Problems (9) and (10):

$${\tilde{y}}^{(1)}\in \underset{y^{(1)}\in {\mathrm{\Sigma }}_N}{\mathrm{Argmin}}\{\Vert P{y}^{(1)}-y^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert y^{(1)}{\Vert }_{(a)}\},$$(9)

$${\tilde{y}}^{(2)}\in \underset{y^{(2)}\in {\mathrm{\Sigma }}_M}{\mathrm{Argmin}}\{{\mathit{\epsilon }}_2\Vert y^{(2)}{\Vert }_{(b)}\},$$(10)

where ${\mathrm{\Sigma }}_N=\{\nu \in {\mathbb{ℛ}}^N:\hspace{0.17em}{\sum }_{i=1}^N{\nu }_i=1,\hspace{0.17em}{\nu }_i\ge 0\},\hspace{0.17em}{\mathrm{\Sigma }}_M=\{\nu \in {\mathbb{ℛ}}^M:{\sum }_{i=1}^M{\nu }_i=1,{\nu }_i\ge 0\}$. Then the following hold:

1. If $\Vert P{\tilde{y}}^{(1)}-{\tilde{y}}^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert {\tilde{y}}^{(1)}{\Vert }_{(a)}<{\mathit{\epsilon }}_2\Vert {\tilde{y}}^{(2)}{\Vert }_{(b)}$, then the optimal solution to Problem (8) is ${\tilde{x}}^T=(\begin{array}{c}{\tilde{x}}^{{(1)}^T}\hspace{0.17em}{\tilde{x}}^{{(2)}^T}\end{array})$ with ${\tilde{x}}^{(1)}={\tilde{y}}^{(1)}$ and ${\tilde{x}}^{(2)}=0$.

2. If $\Vert P{\tilde{y}}^{(1)}-{\tilde{y}}^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert {\tilde{y}}^{(1)}{\Vert }_{(a)}>{\mathit{\epsilon }}_2\Vert {\tilde{y}}^{(2)}{\Vert }_{(b)}$, then the optimal solution is ${\tilde{x}}^{(1)}=0$ and ${\tilde{x}}^{(2)}={\tilde{y}}^{(2)}$.

3. If $\Vert P{\tilde{y}}^{(1)}-{\tilde{y}}^{(1)}{\Vert }_2+{\mathit{\epsilon }}_1\Vert {\tilde{y}}^{(1)}{\Vert }_{(a)}={\mathit{\epsilon }}_2\Vert {\tilde{y}}^{(2)}{\Vert }_{(b)}$, then there are infinitely many solutions to the Problem (8).

Proof.

See Online Appendix A.3 for the complete proof and note that $\Vert {\tilde{y}}^{(2)}{\Vert }_{(b)}\ne 0$, as the solution of the minimization Problem (10) is being searched for on the simplex Σ_M. Nevertheless, ${\tilde{x}}^{(2)}=0$ is optimal in Problem (8) under point 1 of Lemma 1. From this follows that the robust setup assigns zero ranks to new nodes before they appear (if $\Vert x^{(2)}{\Vert }_1$ is unknown) and ranks currently existing pages only. $\square$

According to Lemma 1, the optimization Problem (8) is separable in ${\tilde{x}}^{(1)}$ and ${\tilde{x}}^{(2)}$: Problem (9) is a convex nonsmooth optimization problem, which we can only solve numerically, whereas Problem (10) is solvable theoretically. Due to this property, we are able to state conditions on ${\mathit{\epsilon }}^{(\xi )},\hspace{0.17em}{\mathit{\epsilon }}^{(\psi )},\hspace{0.17em}{\mathit{\epsilon }}^{(\zeta )}$, and ${\mathit{\epsilon }}^{(\chi )}$ implying the robustness of an estimated vector of ranks and guaranteeing that the point ${\tilde{x}}^{(1)}=0,\hspace{0.17em}{\tilde{x}}^{(2)}>0$ is never optimal, thus bounding the ambiguity set. To avoid the optimal solution ${\tilde{x}}^{(1)}=0,\hspace{0.17em}{\tilde{x}}^{(2)}>0$, one needs to make point 2 of Lemma 1 infeasible. In Statement 1, we explicitly solve Problem (10), and we state conditions on parameters ${\mathit{\epsilon }}^{(\xi )},\hspace{0.17em}{\mathit{\epsilon }}^{(\psi )},\hspace{0.17em}{\mathit{\epsilon }}^{(\zeta )}$, and ${\mathit{\epsilon }}^{(\chi )}$, which imply the feasibility of points 1 and 3 of Lemma 1.

Statement 1.

Let s be the total importance ofcurrent N pages; that is, $s={\mathrm{𝟙}}_N^T{x}^{(1)},\hspace{0.17em}s\in [0,1]$. Consider the reformulation (8) and change the variables as $x^{(1)}=s{y}^{(1)},\hspace{0.17em}{x}^{(2)}=(1-s)y^{(2)}$. In this case, $\Vert y^{(2)}{\Vert }_{(b)}={\min }_{\lambda +\mu =y^{(2)}}\left\{\Vert \lambda {\Vert }_2+{\sum }_{j=1}^M\frac{{\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1}{{\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1}|{\mu }_j|\right\}$ and the following holds:

$$\underset{y^{(2)}\in {\mathrm{\Sigma }}_M}{\min }\Vert y^{(2)}{\Vert }_{(b)}=[\begin{array}{c}\begin{array}{l}\frac{1}{\sqrt{M}},\hspace{0.17em}\text{if}\hspace{0.17em}\frac{{\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1}{{\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1}\ge \frac{1}{\sqrt{M}},\hspace{0.17em}\forall j=1,\dots ,M\hfill \\ \underset{j}{\min }\left\{\frac{{\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1}{{\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1}\right\},\hspace{0.17em}\text{if}\hspace{0.17em}\underset{j}{\min }\left\{\frac{{\mathit{\epsilon }}_j^{(\zeta )}+{\mathit{\epsilon }}_j^{(\chi )}+1}{{\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1}\right\}<\frac{1}{\sqrt{M}}.\hfill \end{array}\end{array}$$(11)

Proof.

See Online Appendix A.4 for the proof. $\square$

Thus, points 1 or 3 of Lemma 1 are satisfied if

For any ${\mathit{\epsilon }}_1={\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}$, the minimal value of Problem (9) is smaller than ϵ₁. This is due to the simplex constraints on $y^{(1)},\hspace{0.17em}\Vert y^{(1)}{\Vert }_{(a)}\le \Vert y^{(1)}{\Vert }_1$ and the Perron-Frobenius theorem, which implies that there is such $y^{(1)}$ that $\Vert P{y}^{(1)}-y^{(1)}{\Vert }_2=0$. Due to the same reasons, the minimal value of the objective function in (9) is smaller than ${\sum }_{j=1}^N({\mathit{\epsilon }}_j^{(\xi )}+{\mathit{\epsilon }}_j^{(\psi )})y_j^{(1)}$.

Now, if ${\mathit{\epsilon }}_1={\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}\le 1$, then the minimal value of Problem (9) is clearly smaller than 1. In this case, the Conditions (12) are sufficiently satisfied for ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}\ge \sqrt{M}-1$. This implies that, as the number M of new pages increases, the uncertainty about future pages ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}$ should grow faster than $\sqrt{M}-1$ to guarantee that the solution of Problem (9) is robust to variations in links and to the growth in the number of nodes. For M = 1, this is always satisfied. For higher dimensions (i.e., M > 1), we test this condition by computing Google vectors and the corresponding objective values in (9) using subnetworks of a real-world Movies data set with 4,754 nodes (see Section 5).

Next, for a given network of the size N, we denote the objective value of (9) evaluated at the Google vector x_G by ${\tilde{\varphi }}_G$: note that the Google vector is not necessarily the minimizer of (9); that is, ${\tilde{\varphi }}_G\ge \tilde{\varphi }$ and ${\tilde{\varphi }}_G=\Vert P{x}_G-x_G{\Vert }_2+{\mathit{\epsilon }}_1\Vert x_G{\Vert }_{(a)}$. In line with the first inequality in (12), the following must hold true for the vector x_G robust to variations in network size and structure:

To check the robustness of the Google vector x_G for different M, we model the growth of the Movies network. For this, we take a random sequence of increasing submatrices of its adjacency matrix. We start with low values of N and proceed to the actual size of the network (i.e., n = 4,754). We set a 3% growth rate (i.e., $M=0.03N$ nodes) in Figure 1(a) and a 15% growth rate (i.e., $M=0.15N$ nodes) in Figure 1(b). By computing the Google vector for each of the resulting networks, we verify (13) by plotting the function ${\tilde{\varphi }}_G\sqrt{M}-1$ with respect to the number of new pages M, where the value of ${\tilde{\varphi }}_G$ depends on the network size N. By this, the green area in Figure 1 identifies the values of the uncertainty parameter ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}\ge 0$ satisfying (13).

If the function ${\tilde{\varphi }}_G\sqrt{M}-1$ takes only negative values as in Figure 1(a), then the robustness criterion in (13) is validated for all feasible values of the uncertainty parameter, that is, for ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}\ge 0$. Conversely, the robustness is not validated for the area below the line (red area for feasible parameters in Figure 1(b)). Overall, the considered sequence of Google vectors satisfies the robustness criterion in (13) at the 3% growth rate of the Movies data set, being, however, not robust at the 15% growth rate for some feasible uncertainty parameters.

Further, note that the applicable growth rate obviously depends on the context of the application. In the web ranking context, if robustness is not validated at higher growth rates but is validated at lower rates, as in Figure 1, then search engines can aim to perform more frequent updates of their rankings using the Google approach in order to match the vector with the growth rate (e.g., five times more frequent updates for the case as in Figure 1). Nevertheless, this is not always possible due to the limitations in computing power, as well as high growth rates, and, more importantly, due to the influence of the chosen method on the position of currently existing pages. In particular, the Google vector may imply a very weak bound on the optimal value of Problem (9) and stop distinguishing between ranks of high-importance pages, as shown by Polyak and Timonina (2011). Further benefits of our approach are demonstrated numerically in Sections 6 and 7 after we describe the solution method in Section 4.

Before we proceed with numerical insights, we theoretically analyze the optimality gap between the optimization Problem (5) and its upper bound implied by Theorem 2 in Section 3.

3. Approximation Quality

Several corollaries of Theorems 1 and 2 provide information on the quality of approximation of the non-convex-concave Problem (5). The quality can be measured both via the distance between optimal values and, importantly, via the difference in optimal solutions, that is, in ranking vectors.

Denoting by Δ the optimality gap between Problem (5) and its approximation (8), we, first of all, state a straightforward bound $\mathrm{\Delta }\le \tilde{\varphi }$ directly implied by the optimal value $\tilde{\varphi }=\tilde{\varphi }(\tilde{x})$ of Problem (8) due to the nonnegativity of objective functions. Note that the optimal value $\tilde{\varphi }$ is, in turn, bounded from above by ${\mathit{\epsilon }}_1={\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )}$ and ${\mathit{\epsilon }}_2={\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1$. This is due to simplex constraints on the ranking vector. Furthermore, it is bounded from above by the value ${\mathit{\epsilon }}_1\Vert {\overline{x}}^{(1)}{\Vert }_{(a)}$, where ${\overline{x}}^{(1)}$ is the dominant eigenvector of the matrix P. Additionally, $\tilde{\varphi }\le {\sum }_{i=1}^N({\mathit{\epsilon }}_j^{(\xi )}+{\mathit{\epsilon }}_j^{(\psi )})$. In general, all of these gaps allow us to measure the approximation quality roughly without solving any additional optimization problems.

Next, the optimality gap can be bounded tighter for problems with larger uncertainty about new pages in case the fraction of ranks ${\mathrm{𝟙}}_N^T{x}^{(1)}$ devoted to old pages is known. This fraction is denoted by $s\in [0,1]$ and can be estimated from data similar to the parameter M: as soon as we know how many pages appear on average during the month, we can estimate the fraction of ranks that they achieve on average. Also, this fraction can be evaluated using the approach described in Section 4.2. Given s, the optimality gap can be represented as the difference in two optimization problems (see Corollary 1).

Corollary 1.

Let $s\in [0,1],\hspace{0.17em}{\mathit{\epsilon }}_j^{(\chi )}\ge 1\hspace{0.17em}\forall j$ and ${\mathit{\epsilon }}^{(\chi )}\ge \sqrt{M}$. The optimality gap $\mathrm{\Delta }(s)$ between the initial Problem (5) under constraint ${\mathrm{𝟙}}_N^T{x}^{(1)}=s$ and its convex approximation (8) under the same constraint can be bounded by $v(s)=s\tilde{\varphi }+(1-s){\mathit{\epsilon }}_2{\min }_{y^{(2)}\ge 0,\hspace{0.17em}{\mathrm{𝟙}}_M^T{y}^{(2)}=1}\Vert y^{(2)}{\Vert }_{(b)}-\underset{\sim }{\varphi }$ as $\mathrm{\Delta }(s)\le v(s)$, where $\underset{\sim }{\varphi }$ is

$$\underset{\sim }{\varphi }=\underset{x\ge 0,\hspace{0.17em}{\mathrm{𝟙}}_N^T{x}^{(1)}=s}{\min }\{c\Vert (P+\xi -I_N)x^{(1)}{\Vert }_2+(1-c)\Vert \psi x^{(1)}{\Vert }_2\},\forall c\in [0,1]\hspace{0.17em} \text{and}\hspace{0.17em}\xi ,\hspace{0.17em}\psi \in {\mathrm{\Xi }}^{(l_2)}.$$(14)

Proof.

See Online Appendix A.5 for the proof. Importantly, the same bound v(s) holds for the optimality gap Δ between Problems (5) and (8) if s underestimates the true importance of current pages $s^{*}$ optimal in Problem (5); that is, $\mathrm{\Delta }=\mathrm{\Delta }(s^{*})\le v(s),\hspace{0.17em}\forall s\le s^{*}$. Also, the bound in (A18) in Online Appendix A.5 holds $\forall s^{*}\in [0,1]$. $\square$

Note that the proposed bound v(s) can be interpreted as the difference between two ranking approaches. In the first one, corresponding to the upper bound $s\tilde{\varphi }+(1-s){\mathit{\epsilon }}_2{\min }_{y^{(2)}\ge 0,\hspace{0.17em}{\mathrm{𝟙}}_M^T{y}^{(2)}=1}\Vert y^{(2)}{\Vert }_{(b)}$, we address both structural and dimensionality changes of the network under the condition that the value of s is known. Conversely, the lower bound $\underset{\sim }{\varphi }$ taken with c = 1 generalizes the problem described by Juditsky and Polyak (2012) to the case with $s\in [0\hspace{0.17em}1]$. Therefore, it describes the network with changes in its structure and with a known fraction of ranks $1-s$ corresponding to new appearing pages. Given this, the optimality gap is equal to $\mathrm{\Delta }(1)\le \tilde{\varphi }-\underset{\sim }{\varphi }$ if s = 1 and measures the difference between robustness with and without changes in the network size. Also, if we relax the constraint on s in the lower bound $\underset{\sim }{\varphi }$, then its value will be equal to zero given the feasible solution $x^{(1)}={\mathcal{O}}_N$. If we further relax the constraint on s in the upper bound, then the optimality gap will reduce to $\mathrm{\Delta }\le \tilde{\varphi }$ given point 1 in Lemma 1, that is, to the trivial optimality gap between Problems (5) and (8). Additionally, selecting ψ = 0 and given the feasible matrix ξ, we obtain $\underset{\sim }{\varphi }=0$ due to the fact that the dominant eigenvector exists for any column-stochastic matrix $P+\xi$. Thus, the optimality gap decomposes to the trivial bound $\mathrm{\Delta }(1)\le \tilde{\varphi }$ for s = 1, ψ = 0, and a given $\xi \in {\mathrm{\Xi }}^{(l_2)}$. Further, the optimality gap can be written as

We test the performance of the bound in (15) numerically in Section 8. For this, given samples for matrices ξ and ψ, one could try to solve the nonsmooth convex optimization Problem (14) via subgradient methods. Nevertheless, as the solution would be an approximation providing an upper bound on Problem (14), this would likely result in an underestimation of the gap. As a solution to this numerical issue, we state Propositions 1 and 2, which have closed-form solutions.

Proposition 1.

Let c = 1. The optimal value of Problem (14) satisfies $\underset{\sim }{\varphi }\ge {\underset{\sim }{\varphi }}_1={\min }_{{\mathrm{𝟙}}_N^T{x}^{(1)}=s}\Vert (P+\xi -I_N)x^{(1)}{\Vert }_2$, where the lower bound has a closed-form solution $\forall \xi ,\hspace{0.17em}\psi \in {\mathrm{\Xi }}^{(l_2)}$ such that $\psi \ne 0$ and $(P+\xi -I_N)$ is invertible. The solution is $x^{(1)}=s\frac{{\tilde{P}}^{-1}{\mathrm{𝟙}}_N}{{\mathrm{𝟙}}_N^T{\tilde{P}}^{-1}{\mathrm{𝟙}}_N}$ with $\tilde{P}={(P+\xi -I_N)}^T(P+\xi -I_N)$.

Proof.

By setting c = 1 but keeping $\psi \ne 0$ in the bound in (14), we derive the statement of Proposition 1, where we relax the constraint $x\ge 0$. The optimal solution is a special case of the well-known least-squares problem under linear constraints (e.g., chapter 16.2 in Boyd and Vandenberghe 2018). Note that the optimal value is not necessarily zero, as $P+\xi$ is not column-stochastic on ${\mathrm{\Xi }}^{(l_2)}$, having columns that sum up to values less than or equal to 1 due to nonnegativity of ψ. $\square$

Proposition 2.

Let $c=\frac{1}{2}$. The optimal value of Problem (14) satisfies $\underset{\sim }{\varphi }\ge {\underset{\sim }{\varphi }}_2={\min }_{{\mathrm{𝟙}}_N^T{x}^{(1)}=s}\frac{1}{2}\sqrt{\Vert (P+\xi -I_N)x^{(1)}{\Vert }_2^2+\Vert \psi x^{(1)}{\Vert }_2^2}$, where the lower bound has a closed-form solution $\forall \xi ,\hspace{0.17em}\psi \in {\mathrm{\Xi }}^{(l_2)}$ such that $\psi \ne 0$ and $\tilde{P}={(P+\xi -I_N)}^T(P+\xi -I_N)+{\psi }^T\psi$ is invertible. The solution has the closed form $x^{(1)}=s\frac{{\tilde{P}}^{-1}{\mathrm{𝟙}}_N}{{\mathrm{𝟙}}_N^T{\tilde{P}}^{-1}{\mathrm{𝟙}}_N}$.

Proof.

This statement follows from the nonnegativity of the objective function in (14) as ${(\frac{1}{2}\Vert (P+\xi -I_N)x^{(1)}{\Vert }_2+\frac{1}{2}\Vert \psi x^{(1)}{\Vert }_2)}^2\ge \frac{1}{4}\Vert (P+\xi -I_N)x^{(1)}{\Vert }_2^2+\frac{1}{4}\Vert \psi x^{(1)}{\Vert }_2^2$. We also relax the constraint $x\ge 0$ and use the optimal solution of the constrained least-squares problem (chapter 16.2 in Boyd and Vandenberghe 2018) representing the objective function as $\frac{1}{4}\Vert {(\begin{array}{c}{(P+\xi -I_N)}^T,{\psi }^T\end{array})}^T{x}^{(1)}{\Vert }_2^2$. $\square$

In both propositions, we require invertibility of the matrix $\tilde{P}$. Note that the matrix $P-I_N$ is not invertible, as it is rank-deficient, having columns summing up to zeros. Nevertheless, as invertible matrices are dense in a set of N × N matrices for any norm, one can always find an invertible $P+\xi -I_N$ for nonzero uncertainty parameters ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}_j^{(\xi )}\hspace{0.17em}\forall j$. For this, one requires $\psi \ne 0$. Otherwise, the matrix $P+\xi$ would be column-stochastic by itself, thus, implying the rank deficiency of $P+\xi -I_N$.

Next, given the invertibility of $P+\xi -I_N$, the matrix ${(P+\xi -I_N)}^T(P+\xi -I_N)+{\psi }^T\psi$ is invertible. This is due to the fact that ${(P+\xi -I_N)}^T(P+\xi -I_N)+{\psi }^T\psi =A^{T}A$, where $A^T=({(P+\xi -I_N)}^T,\hspace{0.17em}{\psi }^T)$. Matrix A has rank N, as $P+\xi -I_N$ is invertible. This implies that matrix A^TA is symmetric and positive definite.

Using Propositions 1 and 2 to compute the lower bounds, we compare the value of v(s) imposed by Corollary 1 with the trivial bound $\tilde{\varphi }$. Thus, we derive the conditions on the total importance of old pages $s={\mathrm{𝟙}}_N^T{x}^{(1)}$ leading to the optimality gaps $\mathrm{\Delta }(s)$ lower than $\tilde{\varphi }$. Figure 2 demonstrates these conditions for different growth rates. For the analysis, we use the Computational Complexity data set with n = 789 nodes described in Section 5.1. We observe that Proposition 1 results in a broader range of the importance parameter s, providing optimality gaps lower than $\tilde{\varphi }$ given large growth rates M. Nevertheless, a small number of new pages M with high uncertainty and lower values of s suggest using Proposition 2.

The improvement is achieved in approximately 10% of cases (i.e., s > 0.9) for large M and ϵ in Figure 2(a), whereas this percentage increases as we increase the uncertainty parameters. Similarly, the smallest value of s allowing the reduction of the gap in approximately 20% of cases (i.e., s > 0.8), is obtained at an approximately 12% growth rate in Figure 2(b).

For s = 1, Figure 3 demonstrates the reduction in the ratio $\frac{v(1)}{\tilde{\varphi }}$ and, thus, the improvement of the bound v(1), which is dependent on the values of the uncertainty parameters ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}^{(\psi )}$. In Figure 3, (a) and (b), we observe that the optimality gap bound implied by Proposition 1 (respectively, Proposition 2) decreases by approximately 2% (respectively, 4%) in case of a 15% network growth given low values of ${\mathit{\epsilon }}^{(\xi )}$ and high values of ${\mathit{\epsilon }}^{(\psi )}$. Indeed, the objective value of Problem (20), which we employ for the numerical solution (see the upper bound in (19)), increases in ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}^{(\psi )}$ and is symmetrical for $\mathit{\epsilon }={\mathit{\epsilon }}^{(\xi )}$ and $\mathit{\epsilon }={\mathit{\epsilon }}^{(\psi )}$. However, due to the nonlinear dependency of the lower bound $\underset{\sim }{\varphi }$ on the values of ξ and ψ, the percent reduction of the optimality gap is asymmetrical in the perturbation parameters ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}^{(\psi )}$. This result is analogous to the one presented in Figure 18 in Section 8.

Further, Corollary 2 provides an implicit bound for the solution vector ${\widehat{x}}^{(l_2)}={({\widehat{x}}^{(l_2),(1)},\hspace{0.17em}{\widehat{x}}^{(l_2),(2)})}^T$ of the initial Problem (5) given the optimal value and the solution of the upper bound in (8). Clearly, besides the statement of the Corollary 2, simplex constraints hold for the solution ${\widehat{x}}^{(l_2)}$, implying that $\Vert {\widehat{x}}^{(l_2),(2)}{\Vert }_1=1-\Vert {\widehat{x}}^{(l_2),(1)}{\Vert }_1$.

Corollary 2.

Given ${\mathit{\epsilon }}_j^{(\chi )}\ge 1\hspace{0.17em}\forall j$ and ${\mathit{\epsilon }}^{(\chi )}\ge 1$, let $\tilde{x}={({\tilde{x}}^{(1)},{\tilde{x}}^{(2)})}^T$ be the optimal solution of the reformulation (8), and let $\tilde{\varphi }$ be its optimal value. The following implicit inequality holds for the true solution vector ${\widehat{x}}^{(l_2)}={({\widehat{x}}^{(l_2),(1)},{\widehat{x}}^{(l_2),(2)})}^T$:

$$\underset{\begin{array}{l} \text{Conditions}\hspace{0.17em}(2)-(3)\\ \Vert {[\xi ]}_j{\Vert }_1\le {\mathit{\epsilon }}_j^{(\xi )}\hspace{0.17em}\forall j,\hspace{0.17em}\Vert \xi {\Vert }_F\le {\mathit{\epsilon }}^{(\xi )}\end{array}}{\max }{({\tilde{x}}^{(1)})}^T{(P+\xi -I_N)}^T(P+\xi -I_N){\widehat{x}}^{(l_2),(1)}\le {\tilde{\varphi }}^2.$$(16)

Proof.

See Online Appendix A.6 for the proof. One can also state inequalities $\Vert (P+\xi -I_N){\widehat{x}}^{(l_2),(1)}{\Vert }_2\le \tilde{\varphi }$ and $\Vert (\chi -I_M){\widehat{x}}^{(l_2),(2)}{\Vert }_2\le \tilde{\varphi }$, which hold true for all ξ and χ in the uncertainty set. The first inequality follows directly from Theorems 1 and 2. The second one is shown in Online Appendix A.6. Additionally, if the uncertainty set is large enough (i.e., if ${\mathit{\epsilon }}^{(\chi )}\ge \sqrt{M}$), then one can deduce the bound ${({\tilde{x}}^{(1)})}^T{\psi }^T[\psi {\widehat{x}}^{(l_2),(1)}+(\chi -I_M){\widehat{x}}^{(l_2),(2)}]\le {\tilde{\varphi }}^2$, which combines components corresponding to old and new pages in one statement. $\square$

As Inequality (16) depends solely on the uncertainty matrix ξ, it provides some possibilities for numerical tests demonstrated in Section 8. Conversely, the bounds involving ranks of new pages ${\widehat{x}}^{(l_2),(2)}$ require a generation of random matrices χ and become more complex to validate as the number M of new pages grows. Note also that the worst-case perturbations in the left-hand side of the proposed inequalities cannot be computed explicitly because of the functional form and constraints of the problems.

Lastly, we demonstrate a degenerate case that implies zero optimality gap between a modified Problem (5) and its approximation. For this, we consider the set of constraints ${\mathrm{\Xi }}^{(F)}=\{{\mathrm{𝟙}}_{N+M}^T{[Q]}_j=1\hspace{0.17em}\forall j,\hspace{0.17em}\Vert \xi {\Vert }_F\le {\mathit{\epsilon }}^{(\xi )},\hspace{0.17em}\Vert \psi {\Vert }_F\le {\mathit{\epsilon }}^{(\psi )},\hspace{0.17em}\Vert \zeta {\Vert }_F\le {\mathit{\epsilon }}^{(\zeta )},\hspace{0.17em}\Vert \chi {\Vert }_F\le {\mathit{\epsilon }}^{(\chi )}\}$ and the following optimization problem,

with the optimal solution ${\widehat{x}}^{(F)}={({\widehat{x}}^{(1,F)},\hspace{0.17em}{\widehat{x}}^{(2,F)})}^T$, where ${\widehat{x}}^{(1,F)}$ is an $N\times 1$ vector and ${\widehat{x}}^{(2,F)}$ is an $M\times 1$ vector.

Corollary 3.

Given ${\mathit{\epsilon }}^{(\chi )}\ge 1$, the following bounds hold for the optimization Problem (17):

$$\begin{array}{l}\underset{Q\in {\mathrm{\Xi }}^{(F)}}{\max }\Vert Qx-x{\Vert }_2\le \Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+({\mathit{\epsilon }}^{(\xi )}+{\mathit{\epsilon }}^{(\psi )})\Vert x^{(1)}{\Vert }_2+({\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}+1)\Vert x^{(2)}{\Vert }_2,\\ \underset{Q\in {\mathrm{\Xi }}^{(F)}}{\max }\Vert Qx-x{\Vert }_2\ge \Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}^{(\xi )}\Vert x^{(1)}{\Vert }_2+{\mathit{\epsilon }}^{(\zeta )}\Vert x^{(2)}{\Vert }_2.\end{array}$$

Furthermore, if ${\mathit{\epsilon }}^{(\psi )}=0$ and ${\min }_{x^{(1)}\in {\mathrm{\Sigma }}_N}(\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}^{(\xi )}\Vert x^{(1)}{\Vert }_2)\le \frac{{\mathit{\epsilon }}^{(\zeta )}}{\sqrt{M}}$, then the optimality gap ${\mathrm{\Delta }}^{(F)}={\varphi }^{(F)}-{\tilde{\varphi }}^{(F)}=0$, where

$${\tilde{\varphi }}^{(F)}=\underset{x\in {\mathrm{\Sigma }}_N}{\min }\{\Vert P{x}^{(1)}-x^{(1)}{\Vert }_2+{\mathit{\epsilon }}^{(\xi )}\Vert x^{(1)}{\Vert }_2\}$$(18)

implies the optimality of ${\widehat{x}}^{(F)}={({\tilde{x}}^{(1,F)},\hspace{0.17em}{\mathcal{O}}_M)}^T$, where ${\tilde{x}}^{(1,F)}$ is the optimal solution of Problem (18).

Proof.

See Online Appendix A.7 for the proof and a small-dimensional example. $\square$

In Section 8, we numerically analyze Corollaries 1 and 2 for real-world data sets, considering pessimistic and optimistic scenarios with a boost growth in the number of new pages M. This allows us to assess theoretical and empirical values of optimality gaps, as well as to observe possible true solution patterns.

In Section 4, we propose solution methods for middle- and high-dimensional problems. For this, recall that Problem (8) can be subdivided into two separate smaller-size optimization problems in line with Lemma 1¸whereas Problem (10) can be solved explicitly. Thus, we focus on the convex but nonsmooth optimization Problem (9) of the type described in the works of Ben-Tal and Nemirovski (1998), Boyd and Vandenberghe (2004), and Juditsky and Polyak (2012) with $\mathit{\epsilon }={\mathit{\epsilon }}_1$:

Similar to Juditsky and Polyak (2012), we develop numerical algorithms for the right-hand side in (19), as it bounds the initial Problem (5) and is strictly convex and provides zero optimality gap with (17) under the conditions of Corollary 3.

4. Solution Method

In middle- to large-dimensional cases (i.e., $N={10}^3–{10}^6$), the convex nonsmooth optimization Problem (19) can be solved numerically using available optimization techniques, including interior-point methods, mirror-descent algorithms (Nemirovski and Yudin 1983, Lan et al. 2012) and randomized techniques (Ishii et al. 2005, Nazin and Polyak 2009, Ishii and Tempo 2010). In huge-dimensional cases (i.e., $N={10}^6–{10}^9$), one can employ randomized subgradient algorithms, which are specifically designed for sparse matrices by Nesterov (2014), or first-order extragradient methods, such as mirror-prox and dual interpolation algorithms (Nemirovski 2004, Nesterov 2007, Juditsky and Nemirovski 2012). Moreover, one can use less accurate but extremely fast numerical methods proposed by Juditsky and Polyak (2012), Lei (2014), and Polyak and Timonina (2011). In this section, we propose techniques that allow us to smoothen the optimization Problem (19), and we develop numerical algorithms for the solution of Problem (19) for middle- to high-dimensional cases. Dependent on the application area, different algorithms may be of interest.

4.1. Middle-Size Networks

In middle-dimensional cases, we focus on the upper bound of the optimization Problem (19), which is also valid for the Frobenius-norm formulation and which implies a nonsmooth nonseparable problem:

The optimization Problem (20) is strictly convex given $\mathit{\epsilon }>0$ and, thus, has a unique minimum. We solve Problem (20) using the following steps (see Online Appendix A.8).

Step 1 (Normalization).

Apply a nonstandard normalization instead of simplex constraints $x\in {\mathrm{\Sigma }}_N$. Assuming that one page (e.g., page N) has the highest rank (Polyak and Timonina 2011), one uses the constraints $\mathbb{X}=\{x\in {\mathbb{ℛ}}^N,\hspace{0.17em}{x}_N=1,\hspace{0.17em}x\ge 0\}$ with the normalization x_N = 1 instead of ${\sum }_{i=1}^N{x}_i=1$. For example, updating the ranking compared with the instance of the previous period (e.g., month), one can use the page with the highest score from this instance and leave it constant for several subsequent periods. In this case, updating the highest rank would be due to the power method employed once in several months, whereas other ranks would be due to our robust algorithm.

Step 2 (Smoothing).

Adjust Problem (20) so that avoiding nonsmooth points becomes possible. For this, we consider a perturbed Google matrix $G=\alpha P+(1-\alpha )S$, where S is the doubly stochastic matrix with entries equal to $\frac{1}{N}$ and where α is the damping factor (Haveliwala and Kamvar 2003, Langville and Meyer 2006). This matrix was initially proposed by Brin and Page (1998) and is used by Google in the power method $y^{(k+1)}=G{y}^{(k)}$.

First, by using the Google matrix G instead of P in the subgradient of Problem (20), we reduce the number of nonsmooth points to one, as for the matrix G, one can guarantee the uniqueness of the maximal (in absolute value) eigenvalue ($|\lambda |=1$) and, thus, of the eigenvector $\overline{y}$ corresponding to it (Polyak and Timonina 2011).

Second, we set a rank of a spam page (e.g., x₁) to zero. This goes in line with the practice of many search engines that manually check for obvious instances of spamdexing and remove suspect pages from their indexes in response to user complaints of false matches. Note that the numeration of pages is irrelevant, while the spam page x₁, similarly to the high-ranked page x_N, must correspond to the outline of the network and can be relabeled if necessary. Now, as the Perron-Frobenius theorem for positive matrices implies that the principal eigenvector of the matrix G has only positive entries, we guarantee that $Gx\ne x$ for each chosen x by setting the importance of some nodes to zero. By this, we enforce the uniqueness of the subgradient at each feasible point. Thus, we use the subgradient implied by the problem

which is strictly convex, similar to (20), and we note thatwhere the upper bound follows directly from the triangle inequality, and the lower bound is deduced using the duality of the norm. Bounds (22) imply that the minimal value of the function g(x) on a standard simplex approaches the minimal value of the function f(x) as $\alpha \to 1$. Furthermore, as shown by Polyak and Timonina (2011) for a smooth function, by choosing a specific sequence ${\left\{{\alpha }_k\right\}}_{k=1}^K$, which converges to 1, and updating x_k iteratively, one could aim for the convergence of $g(x_k)$ to $f(x_k)$. Clearly, choosing a constant ${\alpha }_k=\alpha$ leads to the optimal value of (21) in some distance to (20). Nevertheless, this distance is bounded in line with (22).

Now, instead of using subgradient methods that yield δ-accuracy in $O\left(\frac{1}{{\delta }^2}\right)$ iterations, one can employ accelerated gradient descent algorithms that result in δ-accuracy in $O\left(\frac{1}{\sqrt{\delta }}\right)$, a significant improvement compared with the nonsmooth case (Nesterov 2005, Beck and Teboulle 2012, Beck 2017).

Step 3 (Solution).

We reduce the objective of Problem (20), moving along the subgradient $\partial g(x)=\frac{{(G-I)}^T(G-I)x}{\Vert Gx-x{\Vert }_2}+\mathit{\epsilon }\frac{x}{\Vert x{\Vert }_2}$.

For the solution of Problem (20), we propose a coordinate-wise descent algorithm (Algorithm 1) with a near-optimal step size, adjusted at each iteration. For this, we note that the optimal choice of $x_i^{(k+1)}$ would imply the following implicit equation for the step size arising by equating the first derivative of the objective function over the ith coordinate to zero:

where ${[G]}_i,\hspace{0.17em}{[P]}_i$ and e_i are the ith columns of matrices G, P, and I_N, respectively, and where we denote ${\mathrm{\Delta }}_G^{(k)}=G{x}^{(k)}-x^{(k)}$ and ${\mathrm{\Delta }}_P^{(k)}=P{x}^{(k)}-x^{(k)}$. While $x^{(k)}$ converges to the optimality of (21), the distance between $x^{(k)}$ and $x^{(k+1)}$ decreases as $k\to \infty$. We introduce a function ${\lambda }^{(k+1)}=1+\mathit{\epsilon }\frac{\Vert {\mathrm{\Delta }}_G^{(k)}+({[G]}_i-e_i)(x_i^{(k+1)}-x_i^{(k)}){\Vert }_2}{\Vert {[G]}_i-e_i{\Vert }_2^2·\Vert x^{(k)}+e_i(x_i^{(k+1)}-x_i^{(k)}){\Vert }_2}$, which converges to $1+\mathit{\epsilon }\frac{\Vert {\mathrm{\Delta }}_G^{(k)}{\Vert }_2}{\Vert {[G]}_i-e_i{\Vert }_2^2·\Vert x^{(k)}{\Vert }_2}$ as $x^{(k+1)}-x^{(k)}\to 0$. If an oracle would provide the exact value of ${\lambda }^{(k+1)}$ already at the iteration k, then one would have the optimal step size and obtain the value of $x_i^{(k+1)}$ from (23). However, one does not have the exact value of ${\lambda }^{(k+1)}$ as, by itself, it depends on $x_i^{(k+1)}$. Thus, one can only use approximations. Computing the function ${\lambda }^{(k)}$ dependent on $x_i^{(k)}$ and $x_i^{(k-1)}$ at the iteration k, we approximate the left-hand side of Equation (23) as $x_i^{(k+1)}{\lambda }^{(k)}$, which therefore results in an explicit expression for $x_i^{(k+1)}$. We call the chosen step size near-optimal because ${\lambda }^{(k+1)}-{\lambda }^{(k)}\to 0$ as $x^{(k+1)}-x^{(k)}\to 0$, which is in line with Equation (23). Also, as ${\lambda }^{(k)}$ can still be higher or lower than ${\lambda }^{(k+1)}$ at each iteration, we keep track of the objective function value. In particular, we are interested in decreasing the objective value of the nonsmooth strictly convex optimization Problem (20) rather than of the smoothened Problem (21). Thus, we add the requirement $\Vert {\widehat{\mathrm{\Delta }}}_P^{(k+1)}{\Vert }_2+\mathit{\epsilon }\Vert {\widehat{x}}^{(k+1)}{\Vert }_2\le \Vert {\mathrm{\Delta }}_P{}^{(k)}{\Vert }_2+\mathit{\epsilon }\Vert x^{(k)}{\Vert }_2$ at each iteration, where ${\widehat{\mathrm{\Delta }}}_P^{(k+1)}$ and ${\widehat{x}}^{(k+1)}$ are the updated values of ${\mathrm{\Delta }}_P^{(k)}$ and $x^{(k+1)}$, respectively.

Algorithm 1 does not consider any nonsmooth points of the optimized function and does not iterate over pages with zero ranks. Further, Algorithm 1 does not update coordinates, for which the condition $\Vert {\widehat{\mathrm{\Delta }}}_P^{(k+1)}{\Vert }_2+\mathit{\epsilon }\Vert {\widehat{x}}^{(k+1)}{\Vert }_2\le \Vert {\mathrm{\Delta }}_P^{(k)}{\Vert }_2+\mathit{\epsilon }\Vert x^{(k)}{\Vert }_2$ is not satisfied, and, thus, we can break iterations if vector $x^{(k)}$ stops changing. Alternatively, it is possible to implement line search in Algorithm 1 to avoid stagnation: for this, one needs to decrease the step size if the stagnation occurs and continue iterations. For efficiency reasons, we do not use line search in our article.

Algorithm 1

(Coordinate-Wise Descent Algorithm with Near-Optimal Step Size)

Choose starting point $x^{(0)}$ with $x_1^{(0)}=0,\hspace{0.17em}{x}_N^{(0)}=1,x_i^{(0)}\in (0,\hspace{0.17em}1),\forall i\in \overline{2,N-1}$; set the maximal number of iterations K, damping factor $\alpha \in (0,\hspace{0.17em}1)$, and Google matrix $G=\alpha P+(1-\alpha )S$. Let $\widehat{x}=x$.

Set ${\mathrm{\Delta }}_G^{(0)}=G{x}^{(0)}-x^{(0)},\hspace{0.17em}{\mathrm{\Delta }}_P^{(0)}=P{x}^{(0)}-x^{(0)}$ and ${\lambda }_i^{(0)}=1+\mathit{\epsilon }\frac{\Vert {\mathrm{\Delta }}^{(0)}{\Vert }_2}{\Vert x^{(0)}{\Vert }_2·\Vert {[G]}_i-e_i{\Vert }_2^2},\hspace{0.17em}\forall i=\overline{2,N-1}$.

for $k=0,\dots ,K-1$ do $\forall i:$$x_i^{(k)}\ne 0,\hspace{0.17em}i\in \overline{2,N-1}$   

$$\begin{array}{l}{\widehat{x}}_i^{(k+1)}=\max \left\{\frac{1}{{\lambda }_i^{(k)}}\left(x_i^{(k)}-\frac{({[G]}_i-e_i,{\mathrm{\Delta }}_G^{(k)})}{\Vert {[G]}_i-e_i{\Vert }_2^2}\right),\hspace{0.17em}0\right\},\\ {\widehat{\mathrm{\Delta }}}_G^{(k+1)}={\mathrm{\Delta }}_G^{(k)}+({[G]}_i-e_i)({\widehat{x}}_i^{(k+1)}-x_i^{(k)})\hspace{0.17em} \text{and}\hspace{0.17em}{\widehat{\mathrm{\Delta }}}_P^{(k+1)}={\mathrm{\Delta }}_P^{(k)}+({[P]}_i-e_i)({\widehat{x}}_i^{(k+1)}-x_i^{(k)}),\\ \mathbf{i}\mathbf{f}\hspace{0.17em}\Vert {\widehat{\mathrm{\Delta }}}_P^{(k+1)}{\Vert }_2+\mathit{\epsilon }\Vert {\widehat{x}}^{(k+1)}{\Vert }_2\le \Vert {\mathrm{\Delta }}_P^{(k)}{\Vert }_2+\mathit{\epsilon }\Vert x^{(k)}{\Vert }_2\hspace{0.17em}\mathbf{s}\mathbf{e}\mathbf{t}\\ \hspace{1em}{x}_i^{(k+1)}={\widehat{x}}_i^{(k+1)},\hspace{0.17em}{\mathrm{\Delta }}_G^{(k+1)}={\widehat{\mathrm{\Delta }}}_G^{(k+1)}\hspace{0.17em} \text{and}\hspace{0.17em}{\mathrm{\Delta }}_P{}^{(k+1)}={\widehat{\mathrm{\Delta }}}_P^{(k+1)},\\ \hspace{1em}{\lambda }_i^{(k+1)}=1+\mathit{\epsilon }\frac{\Vert {\mathrm{\Delta }}_G^{(k)}+({[G]}_i-e_i)(x_i^{(k+1)}-x_i^{(k)}){\Vert }_2}{\Vert x^{(k+1)}{\Vert }_2·\Vert {[G]}_i-e_i{\Vert }_2^2}.\end{array}$$

$$\mathbf{e}\mathbf{lse}\hspace{0.17em}{\widehat{x}}_i^{(k+1)}=x_i^{(k)}\hspace{0.17em}\mathbf{end}\mathbf{if}$$

end for

Besides Algorithm 1, one can solve the optimization Problem (21) using the projected subgradient algorithm in (24). Both algorithms do not iterate over elements x₁ and x_N. The uniqueness of the subgradient at each iteration is implied by $x_1=0$, similar to Algorithm 1. The chosen step size t_k is constant for all vector coordinates and, thus, is different from the near-optimal coordinate-specific step size of Algorithm 1:

In general, subgradient methods are not necessarily descending in the objective function due to the difficulty in choosing the right step size. Therefore, one keeps track of the best function value at every iteration k; that is, $g^{\mathit{best}}={\min }_{i\in \{1,2,\dots ,k\}}g(x^{(i)})$. This, however, adds complexity and reduces efficiency: the subgradient methods achieve δ-accuracy in $O\left(\frac{1}{{\delta }^2}\right)$ iterations.

In Section 6, we test the stopping rule $g(x^{(k+1)})>g(x^{(k)})$ instead of keeping track of the best function value. Also, we do not iterate over pages with zero ranks, which allows us to enhance the method’s efficiency. The advantage of Algorithm 1 compared with (24) is the near-optimal step size, which is updated iteratively following (23), and the rule $\Vert {\widehat{\mathrm{\Delta }}}_P^{(k+1)}{\Vert }_2+\mathit{\epsilon }\Vert {\widehat{x}}^{(k+1)}{\Vert }_2\le \Vert {\mathrm{\Delta }}_P^{(k)}{\Vert }_2+\mathit{\epsilon }\Vert x^{(k)}{\Vert }_2$, which allows us to reduce the objective value of Problem (20).

4.2. High-Dimensional Networks

Algorithm 1, as a first-order method, is adapted for middle-dimensional cases. However, more efficient but less accurate approaches, which avoid differentials, are necessary in high dimensions. As stated earlier in the text, the method employed in practice uses power iterations with the Google matrix $G=\alpha P+(1-\alpha )S$ of the size N × N (Brin and Page 1998). The parameter α describes the probability of following links provided on a website j, while $1-\alpha$ is the probability of leaving to any random internet page, including staying on the site j. Matrix S is the doubly stochastic matrix with entries $\frac{1}{N}$. The probability $\alpha =0.85$ used by the Google algorithm implies that 15% of users do not follow the links provided on websites.

Nevertheless, a variety of modifications of the Google algorithm are possible in huge dimensionalities. Assuming random transitions to all websites except r pages, which are provided as links on the current website j, one would estimate the transition matrix $\tilde{G}$ as $\tilde{G}=\alpha P+(1-\alpha )\tilde{S}$, where $\tilde{S}$ is a stochastic matrix with entry ${\tilde{S}}_{ij}=\frac{1}{N-r}$ if there is no direct link from j to i and 0 otherwise $\forall i,j$ (Algorithm 2). The transition matrices G and $\tilde{G}$ have only positive entries. Therefore, their eigenvectors, corresponding to the eigenvalue λ = 1 are unique, and power methods converge according to the Perron-Frobenius theorem for positive matrices. The convergence is linear with the rate $\frac{1}{{\lambda }_2}$, where λ₂ is the second eigenvalue of the transition matrix G or $\tilde{G}$. Compared with the Google method, Algorithm 2 distorts less transition probabilities and, thus, is more suitable for less sparse transportation matrices, for example, for applications in scientometrics.

Algorithm 2

(Power Method with Modified Matrix $\tilde{G}$)

Choose a starting point $x^{(0)}$ from the standard simplex. Set maximal number of iterations K.

for $k=0,\dots ,K-1$ do $x^{(k+1)}=\alpha P{x}^{(k)}+(1-\alpha )\tilde{S}{x}^{(k)}$ end for

Another modification of the Google algorithm allows us to take network growth roughly into account. For this, let L be the number of users on the website j, and suppose that there are r direct links from this site to pages $\{1,2,\dots ,r\}$. The users follow the links independently of each other. Suppose further, that there are two additional possibilities for the users: (i) they can leave the known internet completely, or (ii) they can proceed to a site that is not yet ranked by the search engine (i.e., to a very recently appeared website). One can consider possibilities (i) and (ii) as dummy links N + 1 and N + 2 from the site (Figure 4). Note that the number of dummy nodes can be chosen arbitrarily and that the sample space is partitioned into a disjoint union of cells. An extended multinomial model accounting for nondisjoint cells is presented by Dong and Yin (2018).

The probability $\mathcal{P}(l_1,\dots ,l_r,l_{N+1},l_{N+2})$ that exactly l_i users proceed from the site j to i $\forall i$ is

where $p_i>0,\hspace{0.17em}\forall i\in \{1,\dots ,r,N+1,N+2\}$. For a fixed point $({\widehat{l}}_1,\dots ,{\widehat{l}}_r,{\widehat{l}}_{N+1},{\widehat{l}}_{N+2})$, the maximum likelihood estimator of the transition probability p_i is ${\widehat{p}}_i=\frac{{\widehat{l}}_i}{L}$, while the corresponding value of the probability mass function is denoted by $\widehat{\mathcal{P}}$. Importantly, the probabilities ${\widehat{p}}_{N+1}$ and ${\widehat{p}}_{N+2}$ are unobservable, and, thus, we would like to express them through other variables. For this, let us denote the probability ${\sum }_{i=1}^r{\widehat{p}}_i$ (i.e., the ratio of users who follow the existing links on a website j) by α_j and the ratio $\frac{1}{{\left(2\pi L\right)}^{r+1}{\prod }_{i=1}^r{\widehat{p}}_i}$ by β_j. We use the subscript j in α_j to stress the dependency on the current page, but, at the same time, we keep the analogy with α of the Google matrix, as it has a similar sense. Importantly, both α_j and β_j can be estimated from the website’s statistics: the probability ${\widehat{p}}_i$ is obtained by setting a counter of users who follow the links, and L can be evaluated by the number of different IP addresses visiting the site.

By using Stirling’s approximation, ${\widehat{l}}_i!\sim \sqrt{2\pi {\widehat{l}}_i}\hspace{0.17em}{({\widehat{l}}_i)}^{{\widehat{l}}_i}{e}^{-{\widehat{l}}_i}$ for large L and ${\widehat{l}}_i$ in (25), we arrive at the equation ${\widehat{p}}_{N+1}{\widehat{p}}_{N+2}=\frac{{\beta }_j}{{\widehat{\mathcal{P}}}^2}$. Combined with the probability addition rule ${\widehat{p}}_{N+1}+{\widehat{p}}_{N+2}=1-{\sum }_{i=1}^r{\widehat{p}}_i$ and an assumption ${\widehat{p}}_{N+1}>{\widehat{p}}_{N+2}$, it results in the expression for ${\widehat{p}}_{N+1}$, which can be bounded from above by ${\widehat{p}}_{N+1}\le \le \frac{1-{\alpha }_j}{2}\left(1+\sqrt{1-\frac{4{\beta }_j}{{(1-{\alpha }_j)}^2}}\right)$. Note that the overestimation takes the smallest value at the mean of the multinomial distribution (i.e., where $\widehat{\mathcal{P}}$ is at maximum). Now we can construct the jth column of an extended transition matrix $\overline{P}$ of the size $(N+2)\times (N+2)$; that is, ${[\overline{P}]}_j^T=\left(\begin{array}{cccccccc}\frac{{\widehat{l}}_1}{L},& \cdots & \frac{{\widehat{l}}_r}{L},& 0,& \cdots & 0,& {\widehat{p}}_{N+1},& {\widehat{p}}_{N+2}\end{array}\right)$.

Algorithm 3

(Power Method with the Modified Matrix $\overline{P}$)

Construct matrix $\overline{P}$, assuming equally probable transitions from the nodes N + 1 and N + 2 to any node.

Choose a starting point $x^{(0)}$ from the standard simplex. Set maximal number of iterations K.

for $k=0,\dots ,K-1$ do $x^{(k+1)}=\overline{P}{x}^{(k)}$ end for

After building the transportation matrix column by column with an assumption of equally probable transitions from the dummy nodes N + 1 and N + 2 to any node in the network, one computes the dominant eigenvector. Note that the constructed matrix is nonnegative and aperiodic. Furthermore, it is irreducible, as every state in the network is accessible from every other state. Therefore, power iterations described in Algorithm 3 converge. The resulting vector can be used as a starting point in middle-size cases for the estimation of the sum $s={\mathrm{𝟙}}_N^{T}x$, which is an important factor for the optimality gap reduction, as shown in Section 3. Furthermore, it can help to account for the network growth in huge-dimensional cases, where the use of subgradient methods is prohibited and the efficiency is of the utmost importance. In this case, due to the necessity store, maintain, and generate large amounts of data about each website (e.g., number of users and user ratio following existing links), some simplifications of Algorithm 3 can be beneficial. For example, an assumption about equally possible transitions among existing links would replace probabilities $\frac{{\widehat{l}}_i}{L}$ by ${\alpha }_j\frac{1}{r}$ for each website. Further, imposing ${\widehat{p}}_{N+1}={\widehat{p}}_{N+2}$ in (25) would result in ${\widehat{p}}_{N+1}=\frac{1-{\alpha }_j}{2}$. Additionally, if there would be M + 1 dummy nodes, then the transition probability from a website j to any of M new pages would be equal to $\frac{1-{\alpha }_j}{2M}$. In Figure 5, we use Models 1 and 2, as well as the real-world data sets Computational Complexity and Movies (see Section 5) to observe the influence of the ratio of users following existing links on the total importance of N pages and on the uncertainty budget $\Vert \xi {\Vert }_F$. For this, we distort matrix P, multiplying it by $\alpha \in [0,\hspace{0.17em}1]$, and we add rows with probabilities ${\widehat{p}}_{N+1}={\widehat{p}}_{N+2}=\frac{1-\alpha }{2}\hspace{0.17em}\forall j$.

We observe that values of α close to 1 result in high values of the sum s (importance of existing pages) in Figure 5(a): α beyond 0.85 leads to s > 0.9 for all of the tested models. This is especially important, as high values of s result in lower optimality gaps, as shown in Section 8. Clearly, high values of the ratio α also imply low values of the uncertainty $\Vert \xi {\Vert }_F$, which is higher for the Computational Complexity network due to its dimensionality (4,754 compared with 625 of Models 1 and 2 and 789 of the Movies data set; Figure 5(b)).

5. Data

We use both simulated and real-world data sets for the analysis of the ranking solution robust to variations in links and to the network growth. For simulated data, we use two network models introduced in the works of Timonina (2009), Polyak and Timonina (2011), and Juditsky and Polyak (2012). Using real-world data sets on link analysis ranking described in by Borodin et al. (2001), we test the proposed techniques and find implications of network growth for real-world information systems.

5.1. Simulated Data

Figure 6 demonstrates two network models introduced in the works of Timonina (2009), Polyak and Timonina (2011), and Juditsky and Polyak (2012). The models are very convenient for high-dimensional analysis due to their simplified structure allowing realistic sparsity. For both models, nodes are shown as circles, while links are denoted by arrows. In Model 1, node (n, n) represents a dangling vertex, which makes the transition matrix P reducible, leading to nonuniqueness of the eigenvector corresponding to the eigenvalue λ = 1. We assume the ability of equally probable transitions from the node (n, n) to any node in Model 1 similarly to the approach used by search engines.

The transition matrix P corresponding to Model 1 with equally probable transitions from the vertex (n, n) becomes irreducible and aperiodic, which guarantees the uniqueness of the eigenvector corresponding to the eigenvalue λ = 1, as well as the convergence of the well-known power method $x^{(k+1)}=P{x}^{(k)}$ (Timonina 2009, Polyak and Timonina 2011, Juditsky and Polyak 2012). The transition matrix P corresponding to Model 1 is clearly highly sparse.

Model 2 differs from Model 1, as it has only one link from the node (n, n). As the transition matrix corresponding to the Model 2 is periodic, the power method $x^{(k+1)}=P{x}^{(k)}$ does not converge for it (Timonina 2009, Polyak and Timonina 2011).

Both models have closed-form solutions, and, therefore, the performance of numerical algorithms can be easily tested in comparison with the true ranks (Timonina 2009, Polyak and Timonina 2011). We use a standard normalization $x_{i,j}^{*}=\frac{x_{i,j}}{{\sum }_{i,j}{x}_{i,j}}$ for both Model 1 and Model 2 (Figure 7). Note that the ranks correspond to the dominant eigenvector of the matrix P. However, they are not necessarily the optimal solution of the optimization Problem (20).

5.2. Real-World Data Sets

We use Computational Complexity (CompCom) and Movies (online film repositories) data sets, downloaded from http://www.cs.toronto.edu/t~sap/. CompCom data represent 789 journals on computational complexity with a total of 1,449 links, while Movies has 4,754 online film repositories and 17,198 links: transition matrices derived based on the data sets are both very sparse and suitable for our purposes. The optimal value of Problem (20) with the transition matrix implied by CompCom data and ϵ = 1 is equal to 0.0587, while the optimal value implied by the Movies data is 0.0288: these values and their corresponding solution vectors are used as benchmarks in Sections 6 and 7.

6. Convergence Results

We have performed extensive numerical analysis for the convergence of the proposed algorithms and their variations (Figures 8–12). For this, we have employed both simulated and real-world data sets. In this section, we provide the most representative convergence results.

We compare a high-accuracy solution of Problem (20) computed using MOSEK optimization software with (i) the scores implied by Algorithm 1, (ii) ranks estimated by the algorithm of Juditsky and Polyak (2012), as well as with (iii) the approximation obtained by the PageRank algorithm. As expected, Algorithm 1 is monotonically decreasing in the function value of Problem (20) and provides high-quality optimal value approximation results. Figure 8 demonstrates iterations of the algorithm for the Movies data set: in Figure 8(a), one can see that the optimal value is lower than both that of the Google algorithm and of Juditsky and Polyak (2012) after four iterations; Figure 8(b) demonstrates the increasing number of zero components, making the algorithm more efficient from iteration to iteration.

As Algorithm 1 is descending with near-optimal step size, it aims for the approximation of the optimal solution of Problem (20) by reducing the objective value at each iteration. Furthermore, the value of the modified objective in (21) can be decreasing in parallel due to the choice of the subgradient (Figure 9(a)). This, however, does not imply that the monotonicity condition $\Vert {\widehat{\mathrm{\Delta }}}_P^{(k+1)}{\Vert }_2+\mathit{\epsilon }\Vert {\widehat{x}}^{(k+1)}{\Vert }_2\le \Vert {\mathrm{\Delta }}_P^{(k)}{\Vert }_2+\mathit{\epsilon }\Vert x^{(k)}{\Vert }_2$ from Algorithm 1 can easily be omitted: neglecting this condition, one does not necessarily obtain the convergence (Figure 9(b)).

Further, we analyze the subgradient algorithm in (24) on Models 1 and 2 with N = 10,000. The subgradient method is not necessarily descending in the objective value, and we test different possible step sizes for its convergence. In Figures 10(a) and 11(a), we consider the square summable but not summable step sizes satisfying ${\sum }_{k\to \infty }{t}_k^2<\infty ,\hspace{0.17em}{\sum }_{k=1}^{\infty }{t}_k=\infty$. In Figures 10(b) and 11(b), we analyze the Polyak step size of the type $t_k=\frac{1}{k\Vert \partial g(x^{(k)}){\Vert }_2^2}$. We use the stopping rule $g(x^{(k+1)})>g(x^{(k)})$ in both Figures 10 and 11, and we demonstrate the performance of the algorithm without it in Figure 10(a) for the Polyak step size. Note that the stopping rule allows us to avoid keeping track of the best function value, which is less efficient.

As seen in Figure 11, the subgradient method in (24) does not distinguish between high-importance pages, similar to the Google algorithm (Polyak and Timonina 2011). Furthermore, it results in an objective value of (20) equal to approximately 0.0107, which is approximately 5% higher than the value provided by Algorithm 1, that is, 0.0102 in Figure 12(a).

Next, we test Algorithms 1–3 on Model 2 with N = 10,000 in Figure 12. As expected, Algorithm 1 results in the smallest objective value of (20). The vector minimizing (20) corresponds to the ranking that is robust to network perturbations in size and in structure. Importantly, comparing the solutions of Algorithms 1–3 with the dominant eigenvector $\overline{x}$ of matrix P in Figure 12(b), we observe that the robust vector has the largest l₂-distance to $\overline{x}$. Indeed, the larger the parameter ϵ, the further the two vectors are. Overall, Algorithm 1 is suggested for the computation of the robust vector in middle-dimensional cases due to its accuracy. In high dimensions, the robust vector can be approximated using Algorithms 2 and 3, as well as the Google algorithm. In this case, the robustness holds for small values of ϵ and regular matrices P (Juditsky and Polyak 2012): the conditions that guarantee the proximity of the true solution to $\overline{x}$.

In general, the problem of computing the dominant eigenvector of P is different from the problem of finding the robust ranking. In order to recover the dominant eigenvector of the unperturbed matrix P in high dimensions, one can implement fast iterative algorithms $x^{(k+1)}={\overline{P}}_k{x}^{(k)}$ using a sequence of transition matrices ${\overline{P}}_k$ converging to P. Such iterative regularization schemes are discussed in the works of Timonina (2009), Polyak and Timonina (2011), and Juditsky and Polyak (2012), and they remind methods for solving variational inequalities are discussed in the work of Bakushinskij and Polyak (1974).

7. Implications for Information Systems

Using the real-world data sets Movies and CompCom, we analyze the impact of the network growth on ranking. The network Movies contains 4,754 web repositories for films, while the data set CompCom contains 789 journals on computational complexity. As there is no historical data on when the nodes (i.e., movie repositories or journals) are created, we consider multiple growth paths by selecting random submatrices of increasing dimension from the transportation matrix. Afterward, we use both the Google algorithm and our Algorithm 1 to obtain ranking vectors for the networks. Figure 13 demonstrates that, already after the first iteration, Algorithm 1 provides a finer approximation for the robust vector as compared with the Google algorithm. This is true for all the growth values except very low ones, for which both methods perform similarly (Figure 13(a)). Furthermore, the objective value obtained via Algorithm 1 is lower for each surplus in the number of nodes, not only for a particular path in Figure 13(a), but also for the average among 100 random growth trajectories in Figure 13(b).

Next, comparing ranking vectors estimated via different methods with the true robust ranking obtained via MOSEK software, we see that both the Google algorithm and the algorithm of Juditsky and Polyak (2012) underestimate ranks of high-importance pages and neglect the clustering effect, which can be observed for the true ranks (Figure 14, (a), (b), and (d)). We also observe that the vector robust to the network growth does not coincide with dominant eigenvectors of G or P (Figure 14, (a) and (d)). Thus, the robust vector cannot be found simply by power iterations. Also, one cannot apply the highly efficient algorithms proposed in Polyak and Timonina (2011).

Figure 14(c) demonstrates the ranking obtained via Algorithm 1. The ranking vector clearly distinguishes between three groups of nodes: (i) high-importance nodes, (ii) nodes of medium importance, and (iii) unimportant nodes receiving zero ranks. This corresponds to the clustering observed for the true ranks in Figure 14(d). Similar effects can be observed for the CompCom data set.

Using the Movies data set, the optimal value of Problem (20) is 0.0288. After four iterations, the approximate solution supplied by Juditsky and Polyak (2012) attains the value 0.0379; Algorithm 1 provides the vector with the function value 0.0356, while the Google algorithm results in 0.0424. Thus, Algorithm 1 results in a 6% improvement compared with the solution of Juditsky and Polyak (2012) and in 16% improvement compared with the Google algorithm. In the CompCom experiment, the optimal value of Problem (20) is 0.0587. The corresponding value for the solution supplied by Juditsky and Polyak (2012) (after three iterations) is 0.0756. Algorithm 1 provides the vector with the function value 0.0703, while the Google algorithm results in 0.0814. The percentage improvements are 7% and 14%, respectively.

Lemma 1 and Statement 1 play an important role here. They allow us to claim which of the estimated vectors are robust to the network growth via computing the values for the uncertainty parameters ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}$ satisfying the robustness. Similar to Figure 1, Figure 15 demonstrates the areas where the estimated ranks are robust to variations in links and to the network growth and where the robustness is not validated. We observe that, for ${\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}\ge 0$ (i.e., for all feasible values of uncertainty), the vector obtained via Algorithm 1 satisfies the robustness criterion in (13) at a 15% growth rate (Figure 15(a)). However, the Google vector is not robust to variations in links and to the same growth in the number of nodes for a large portion of uncertainty parameters (Figure 15(b)). For example, the Google vector does not satisfy robustness criteria for the network with 260 nodes and uncertainty $0\le {\mathit{\epsilon }}^{(\zeta )}+{\mathit{\epsilon }}^{(\chi )}\le 0.5$ at 15% expected growth. Furthermore, the objective function value corresponding to the vector supplied by Algorithm 1 can be much lower than the value corresponding to the Google vector (see Figure 11). Thus, Algorithm 1 provides a finer approximation of the robust ranking vector than the Google algorithm.

8. Numerical Analysis of the Approximation Quality

In this section, we numerically assess the approximation quality of our methods. Invoking Corollary 1, we have a theoretical foundation to bound the optimality gap $\mathrm{\Delta }(s)$ given the total importance of current pages $s={\mathrm{𝟙}}_N^T{x}^{(1)}$. We denote this bound by v(s). Importantly, the optimality gap between Problems (5) and (8) can be bounded by v(s) if the value of s is underestimated with respect to the true importance of current pages $s^{*}$; that is, $\mathrm{\Delta }(s^{*})<v(s)$ if $s<s^{*}$ (see Online Appendix A.5). Conversely, if the importance s is overestimated, that is, if $s>s^{*}$, then the true optimality gap $\mathrm{\Delta }(s^{*})$ can be bounded by a value in some distance δ_v to v(s). The value of δ_v would converge to zero in the number of new pages M if the value of the uncertainty parameter ${\mathit{\epsilon }}_2={\mathit{\epsilon }}^{(\chi )}+{\mathit{\epsilon }}^{(\zeta )}+1$ would stay constant. This, however, is not the case due to the conditions of Corollary 1.

Overall, the optimality gap between Problems (5) and (8) can be bounded in line with Inequality $\text{A}18$ in Online Appendix A.5, for which the condition of Lemma 1 (point 1) should hold. Its lowest value is suggested for high growth rates M.

Further, we consider pessimistic and optimistic scenarios, which we distinguish in terms of over- or underestimation of the optimality gap given the total importance of current pages s.

Pessimistic Scenario.

We use Propositions 1 and 2 and guarantee that $\underset{\sim }{\varphi }\ge {\underset{\sim }{\varphi }}_1$ and $\underset{\sim }{\varphi }\ge {\underset{\sim }{\varphi }}_2$. Thus, the optimality gap $\mathrm{\Delta }(s)$ is overestimated. We refer to this case as the pessimistic one.

Optimistic Scenario.

We relax the lower bounds in Propositions 1 and 2, allowing perturbation matrices ξ and ψ to take values outside of the feasibility region ${\mathrm{\Xi }}^{(l_2)}$, keeping the column-stochasticity constraints only. As the optimality gap can be underestimated via this approach, we refer to this case as the optimistic one.

In both scenarios, we search for such values of s that lead to $v(s)\le \tilde{\varphi }$ (Figure 16 for the optimistic case using Propositions 1). Otherwise, the trivial bound $\tilde{\varphi }$ is stronger than the bound imposed by Corollary 1.

Further, we consider the boost growth setup, in which we require $M>>N$ (Figure 16(b)). We analyze the behavior of the optimality gap bounds for the case when ϵ₂ depends on M, as required by Corollary 1. Specifically, we impose ${\mathit{\epsilon }}^{(\chi )}=\sqrt{M},\hspace{0.17em}{\mathit{\epsilon }}^{(\psi )}=0.5$, and ${\mathit{\epsilon }}^{(\zeta )}=0$, which implies that new pages do not give a part of their possible rank to current pages. By this, we satisfy the conditions of Corollary 1, using minimal possible values for the analysis. We also take into account that only the sum ${\mathit{\epsilon }}^{(\chi )}+{\mathit{\epsilon }}^{(\zeta )}$ plays a role for the value of the optimality gap in (A18) and in Problem (8). Additionally, ${\mathit{\epsilon }}^{(\zeta )}=0$ stress-tests our ranking approach, as dangling vertices imply the reducibility of the transition matrix Q (Polyak and Timonina 2011).

Figures 2, 3, and 17 refer to the pessimistic scenario. In Figure 17(a), we consider the overestimated total importance of current pages s = 1 and observe the convergence of the ratio $\frac{v(1)}{\tilde{\varphi }}$ of two upper bounds for the optimality gap $\mathrm{\Delta }(1)$ in the number of new pages M. The convergence is not to zero due to the dependency of the uncertainty parameter ϵ₂ on M, as suggested in Online Appendix A.5.

We observe that Proposition 1 leads to an approximately 5% lower value for v(1) compared with the bound $\tilde{\varphi }$ in the boost growth case. The bound increases in relative terms if the total importance of current pages s decreases (Figure 17(b)).

Moreover, overestimating $s^{*}$ by s = 1 in Figure 17(c), we bound the optimality gap $\mathrm{\Delta }(s^{*})$ in line with Inequality (A18) and demonstrate the area where it belongs with yellow. Importantly, the resulting bound appears to be stronger than the bound $v(s^{*})$ estimated via the use of Proposition 1. This depends on the chosen set of parameters and follows from the fact that $v(1)\le v(s^{*})$, $\forall s^{*}\in [0,1]$ (Section 3 and Online Appendix A.5).

Further, Figure 18 demonstrates the dependency of the optimistic optimality gap on uncertainty parameters ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}^{(\psi )}$ in the boost growth setup. By solving the optimization Problem (20) for given ${\mathit{\epsilon }}^{(\xi )}$ and ${\mathit{\epsilon }}^{(\psi )}$, we generate matrices ξ and ψ, so that Q stays column-stochastic, and compute the value ${\underset{\sim }{\varphi }}_1$ in Proposition 1. We show the estimate in absolute and in relative terms in Figure 18, (a) and (b), respectively. In Figure 18(a), we see that larger uncertainty parameters clearly lead to higher gaps in absolute terms. Similar to Figure 3 for the pessimistic case, we observe an asymmetrical behavior suggesting the largest gap reduction for low values of the parameter ${\mathit{\epsilon }}^{(\xi )}$ and high values of ${\mathit{\epsilon }}^{(\psi )}$ (Figure 18(b)). Compared with the maximal reduction of approximately 30% in Figure 18(b), the maximal reduction achieved by Proposition 2 is approximately 15%.

Beyond the assessment of pessimistic and optimistic scenarios, we analyze the difference in optimal solutions between the initial Problem (5) and its approximation (19) for the Movies data set (Figure 19). Using Corollary 2, we generate 10,000 pairs $(\xi ,{\widehat{x}}^{(l_2),(1)})$ and verify whether they satisfy Inequality (16). We can exclude all those ${\widehat{x}}^{(l_2),(1)}$ that do not satisfy (16) at least for one generated ξ.

Figure 19 demonstrates that the ranks of current pages start equalizing on average as a result of an increase in uncertainty levels. For this, we suppose that 100%, 10%, or 1% of high ranks $\tilde{x}$ are overestimated. The score equalization would also be the result of an increase in the uncertainty parameter ${\mathit{\epsilon }}^{(\xi )}$ (Juditsky and Polyak 2012). Importantly, the only point that definitely satisfies (16) for all perturbation matrices ξ is our robust solution ${\tilde{x}}^{(1)}$.

9. Conclusion

The World Wide Web is constantly changing: the amount of information on the web is growing, new pages appear, and some websites (old or spam ones) disappear. According to the Web Server Survey, the average speed of internet growth is around 1%–3% per month. The initial idea of the PageRank algorithm is based on the interchange of ranks between web pages in line with probability theory rules and results in a ranking that does not account for the rapid growth of the network and, furthermore, is not necessarily unique. In early works on the topic, the original problem is replaced with finding the principal eigenvector of a modified matrix with only positive components to impose uniqueness of the solution in the PageRank formulation. The modified matrix assumes direct transitions between pages, which do not have direct links between one another is known as the Google matrix and is broadly used in practice. However, the principal eigenvector may be significantly influenced by this difference in matrices, and, importantly, for specific network structures, the eigenvector of the Google matrix may stop distinguishing between high-ranked pages that build the core of search engine results. In this article, we account for the uncertainty in the transition matrix, and we formulate the internet ranking problem as an optimization problem robust to variations in links and to the network growth.

We propose an l₂-norm robust reformulation of the PageRank problem for uncertain networks with fixed growth rates, and we bound this reformulation from above by a separable convex but nonsmooth optimization problem. The separability of the upper bound plays a crucial role for the robustness: one can conclude under which conditions on the uncertainty about future pages the estimated ranking is robust to variations in links and in the network structure and under which conditions it is not. We test the robustness of the Google vector, the ranking proposed by Juditsky and Polyak, and the ranking obtained via our algorithm.

We demonstrate that the formulation robust to network growth imposes an upper bound on the l₂-norm formulation robust to perturbations solely in links. We consider smoothing techniques for the subgradient method, avoiding all nonsmooth points, in order to solve the formulated optimization problem in middle-dimensional cases via a coordinate-wise descent method with near-optimal step size. We also address high-dimensional cases and propose two variations of the Google algorithm. Their performance depends on the unperturbed network structure and, thus, on the context of application. This provides a variety of opportunities for future research. In this article, we use simulated data and the real-world data sets Movies and CompCom to analyze our algorithms and assess the impact of the network growth on ranking.

By assessing the approximation quality from both theoretical and numerical perspectives, we observe that the optimality gap takes its smallest value in case of a rapid growth in the number of pages, while the size of the optimality gap reduction depends on the total importance of current pages. Also, we demonstrate a case when the optimality gap takes zero value: this happens given a specific structure of the uncertainty matrix.

The PageRank problem has a wide range of possible applications (Agarwal and Venkatesh 2002; Pant and Srinivasan 2010, 2013; Xu et al. 2011; Sundararajan et al. 2013; Dhar et al. 2014; Zhu et al. 2014; Pant and Sheng 2015; Van Vlasselaer et al. 2017; Shah et al. 2020). From theoretical and numerical perspectives, future research can focus on different types of structural perturbations of transition matrices (Juditsky and Polyak 2012) and randomized techniques for the computation of the robust stationary distribution in high-dimensional cases (Polyak 2001, Polyak and Tempo 2001, Nazin and Polyak 2009, Bomze et al. 2021).