Iterative Collaborative Filtering for Sparse Matrix Estimation

Introduction

As a prototype for such a problem, consider a noisy observation of a social network where observed interactions are signals of true underlying connections. We might want to predict the probability that two users would choose to connect if recommended by the platform, e.g. LinkedIn. As a second example, consider a recommendation system where we observe movie ratings provided by users, and we may want to predict the probability distribution over ratings for specific movie-user pairs. A popular collaborative filtering approach suggests using “similarities” between pairs of users to estimate the probability that a connection is formed or the probability a user likes a particular movie. Traditionally, the similarities between pair of users in a social network is computed by comparing the set of their friends, or in the context of movie recommendation, by comparing commonly rated movies. In the sparse setting, most pairs of users have no common friends, or most pairs of users have no commonly rated movies; thus there is insufficient data to compute the traditional similarity metrics.

In this work, the primary interest is to provide a principled way to extend the simple, intuitive approach of computing similarities between pair of users or items in order to perform sparse matrix estimation via nearest neighbor collaborative filtering. We propose to do so by incorporating information within a larger radius neighborhood of the data graph rather than restricting only to immediate neighbors. This variation of collaborative filtering and its analysis in this work can be viewed as a natural extension of the work by in the context of stochastic block model and for traditional collaborative filtering.

The primary contribution of this work is an analysis of an iterative collaborative filtering algorithm in the sparse regime. We consider the setting of a latent variable model where the matrix $F=[F(u,v)]$ can be described by a latent function $f$ evaluated over latent variables associated to the coordinates. In particular, we assume that $F(u,v)=f(\theta_{u},\theta_{v})$ where $f$ is a piece-wise Lipschitz function, and $\theta_{u},\theta_{v}\in$ are coordinate latent variables sampled uniformly at random. Details of the model are described in Section 2.

Algorithmically and methodologically, our work builds on , which estimates clusters of the stochastic block model by computing distances from local neighborhoods around vertices. We improve upon their algorithm and analysis to provide bounds on the maximum entrywise estimation error for the general latent variable model with finite spectrum. This includes a larger class of generative models such as mixed membership stochastic block models, in contrast to their work which focuses on the stochastic block model with non-overlapping communities. We note that the algorithm considered in this work, uses the knowledge of which entries are observed and which are not, in line with the literature on matrix estimation. In the setting of clustering cf. , such a knowledge is absent from the purview of the algorithm.

Withthe exception of a few recent results, by and large the literature on matrix estimation has focused on providing estimation error bounds with respect to the normalized Frobenius norm. In contrast, we provide bounds on the max entry-wise estimation error which is a lot more challenging. Our bounds are restricted to the latent variable model, while the traditional matrix estimation literature considers the underlying matrix to be an arbitrary instance from the family of (approximately) low-rank matrices with ‘incoherence’-like conditions. Indeed, understanding the relationship between these two seemingly different model classes remains an important direction for future work.

A weaker version of this result was published in the NeurIPS conference as . In contrast, this paper provides sharper bounds for both the MSE and max-norm error that improves the exponent in the convergence rates. We have also included a perturbation analysis of the algorithm that shows under “adversarial” bounded noise, the error scales gracefully with the bound on the noise. This enables analysis of our work for the approximately low-rank setting. We have also included a modified algorithm that achieves the same rates with a reduced computational complexity, and we have shown extensions of our results to relaxed modeling assumptions on the latent variable model. We have added empirical evaluation of our method compared with state-of-art methods.

2 Related Work

The related work includes that of matrix estimation or completion, collaborative filtering, and graphon estimation arising from the asymptotic theory of graphs. We provide a brief overview of prior works for each of these topics.

In the context of matrix estimation or completion, there has been much progress under the low-rank assumption and additive noise model. Most theoretically founded methods are based on spectral decompositions or minimizing a loss function with respect to spectral constraints, c.f. . In a nutshell, this collection of works establishes that if the underlying matrix has rank $r$, then it can be estimated so that the estimator has normalized Mean Squared Error (MSE) going to as $n\to\infty$ as long as $p=\Omega(rn^{-1}\log n)$. Furthermore, showed that $\omega(rn^{-1})$ samples are necessarily required for such a guarantee. These near optimal sample complexity results hold when the noise in each entry of the matrix is independent and identically distributed. For the setting of generic noise and the general latent variable model where the latent function is analytic, provide an estimator for which the MSE decays to as $n\to\infty$ as long as $p=\Omega(n^{-1}{\sf poly}(\log n))$.

The guarantee with respect to MSE does not necessarily guarantee recovery of all entries accurately. Indeed, bounding max entrywise error provides such a guarantee as established by our result. In parallel with our work, there has been recent progress on developing matrix estimation methods that provide max entrywise bounds for matrices with rank $r$. In particular, for sufficiently ‘nice’ rank $r$ matrices, establish that a simple spectral algorithm can recover the matrix with max entrywise error decaying to as long as $p=\Omega(\log n/n)$. Indeed, improving such max entrywise guarantee has been actively pursued over the past few years witnessed in the growing body of works, cf. , , , and .

The collaborative filtering method has been successfully employed across industry applications (Netflix, Amazon, Youtube) due to its simplicity and scalability, c.f. ; however the theoretical results have been relatively sparse. We call special attention to the recent works by which provide a non-parametric statistical perspective for the traditional collaborative filtering method. In particular, they suggest that the practical success of these methods across a variety of applications may be due to its ability to capture local structure like the classical nearest neighbor or kernel regression method. They establish that as long as the latent function $f$ is Lipschitz, the MSE of the resulting estimator decays to as $n\to\infty$ as long as $p=\omega(n^{-\frac{1}{2}})$. A key limitation of this approach is that it requires a dense dataset with sufficient entries in order to compute similarity metrics, requiring that each pair of rows or columns has a growing number of overlapped observed entries, which does not hold when $p=o(n^{-1/2})$.

Graphons emerged as the limiting object of a sequence of large dense graphs, c.f. , with recent work extending the theory to sparse graphs, c.f. . In the graphon estimation problem, one observes a single instance of a random graph sampled from an underlying latent variable model, and the goal is to estimate the function that governs the edge probabilities of the graph. provide minimax optimal rates for graphon estimation; however a majority of the proposed estimators are not computable in polynomial time, since they require optimizing over an exponentially large space (e.g. least squares or maximum likelihood), c.f. . provides a polynomial time method based on degree sorting in the special case when the expected degree function is monotonic. analyzes universal singular value thresholding (USVT) for graphon estimation in settings that the spectrum decays quickly, showing convergence rates which matches the minimax optimal rate for low dimensional smooth functions.

Stochastic block model (SBM) parameter estimation is an instance of graphon estimation, where the underlying function has a specific structure. Under the SBM, each vertex is associated to one of $r$ community types, and the probability of an edge is a function of the community types of both endpoints. This implies that the edge probability function is block constant. Estimating the $n\times n$ parameter matrix becomes an instance of matrix estimation with a technical distinction – all entries are fully observed, i.e. each edge is present (1) or absent (0). In SBM, the expected matrix is at most rank $r$ due to its block structure. Precise thresholds for cluster detection (better than random) and estimation have been established by . As mentioned before, our work, both algorithmically and methodically is closely related to their work. The mixed membership stochastic block model (MMSBM) allows each vertex to be associated to a length $r$ vector, which represents its weighted membership in each of the $r$ communities. The probability of an edge is a function of the weighted community memberships vectors of both endpoints, resulting in an expected matrix with rank at most $r$. Recent work by provides an algorithm for weak detection for MMSBM with sample complexity $r^{2}n$, when the community membership vectors are sparse and evenly weighted. They provide partial results to support a conjecture that $r^{2}n$ is a computational lower bound, separated by a gap of $r$ from the information theoretic lower bound of $rn$. This gap was first shown in the simpler context of the stochastic block model . proposed a spectral clustering method for inferring the edge label distribution for a network sampled from a generalized stochastic block model. When the expected function has a finite spectrum decomposition, i.e. low rank, then they provide a consistent estimator for the sparse data regime, with $\Omega(n\log n)$ samples.

In the above discussion, we have focused primarily on the sample complexity required for consistent estimation, i.e. the scaling of the number of samples required ($pn$) such that the normalized estimation error such as the MSE or max-norm goes to . When consistent estimation is feasible, we can further consider the rate of decay of the error guarantees. To that end, we provide a brief overview of the minimax scaling with respect to boudns on the MSE. identifies a minimax lower bound on the scaling of the MSE for a generic matrix estimation task characterized by the nuclear norm of the target matrix. In particular, for symmetric matrices with nuclear norm bounded by $\delta$, the minimax MSE scaling is lower bounded by \min\big{(}\frac{\delta}{\sqrt{n^{3}p}},\frac{\delta^{2}}{n^{2}},1\big{)}; furthermore argues that the universal singular value thresholding achieves this scaling. This bound holds even in the scenario where observed entries are noiseless. This characterization however is loose for the setting of low-rank matrices. Observe that for rank $r$ symmetric matrices with entries bounded in $$, the nuclear norm can scale as$n\sqrt{r}$; resulting in a bound of$\sqrt{\frac{r}{np}}$(for small enough$p$) . For the setting of rank$r$matrices with noiseless observations, provide an estimator with MSE scaling as$\frac{r}{np}$for$p=\Omega(1/n)$. This points to the fact that the class of matrices with bounded nuclear norm is more complex than the class of rank$r$matrices with bounded entries. In the setting of low rank graphon estimation (i.e. binary observations), show a minimax lower bound on the MSE scaling as$\frac{\log r}{pn}$for small enough$p=\Omega(\log r/n)$; however the existence of a computationally efficient estimator that achieves this lower bound under the more general noise setting of graphon estimation is still an open research direction.

Setup

Assume that each $u\in[n]$ is associated to a latent feature variable $\theta_{u}\sim U$, which is drawn independently across indices $[n]$ uniformly on the unit interval. We assume that the expected data matrix can be described by the latent function $f$, i.e. $F(u,v)=f(\theta_{u},\theta_{v})$, where $f:^{2}\to$ is a symmetric bounded function. The symmetry assumption can be easily relaxed but is assumed for ease of notation in the analysis. The latent function $f$ is assumed to be fixed and independent of the dimension $n$. We additionally impose local neighborhood properties that are primarily used in the nearest neighbor portion of the analysis. We will assume that $f$ is Lipschitz, but this assumption can be relaxed as discussed in Section 2.2.

Low Rank.

We assume that the latent function $f$ has finite spectrum with rank $r$ when regarded as an integral operator, i.e. for any $\theta_{u},\theta_{v}\in$,

The finite spectrum assumption also implies that the model can be represented by latent variables in the $r$ dimensional Euclidean space, where the latent variable for node $i$ would be the vector $(q_{1}(\theta_{i}),\dots q_{r}(\theta_{i}))$, and the latent function would be bilinear, having the form

This condition also implies that the expected matrix $F$ is low rank, which includes scenarios such as the mixed membership stochastic block model and finite degree polynomials. The function $f$ is fixed with respect to $n$, the rank $r$ is assumed to be finite in the low rank setting.

Approximately Low Rank.

More generally, we shall consider approximately low-rank $f$ cf. . Specifically, for a given $\varepsilon>0$, a symmetric function $f$ is said to have $\varepsilon$-approximate rank $r$ if

2 Discussion on Latent Variable Model

The latent variable model assumes a random generative model on the underlying matrix $F$, as opposed to the typical deterministic incoherence style conditions found in the literature. The generative model assuming i.i.d. sampled latent variables and boundedness of the eigenfunctions of $f$ guarantee similar properties as incoherence with high probability, as any single row or column will not dominate the signal in a way that deviates too much from the typical values of $f$. The i.i.d. sampling assumption on the latent variables is used in analyzing the local neighborhoods of the observation graph, however this assumption can likely be replaced by regularity assumptions over the empirical distribution of the latent factors for large $n$, e.g. if the latent factors are close to a typical sample set from a well-behaved underlying distribution.

The Lipschitzness assumption of $f$ together with the assumption that $\theta_{u}\sim U$, guarantees that for any given $u\in[n]$ there are sufficiently many other coordinates $v\in[n]$ such that the observed entries are similar across both rows or columns. These assumptions can be relaxed as long as the key property of “sufficiently many similarly behaving coordinates” is maintained. As examples, a piecewise Lipschitz function $f$ or a setting with finite latent types would also satisfy the needed local neighborhood properties. Similarly, the scalar assumption on the latent variables and the uniform distribution $U$ are not crucial and can be relaxed to i.i.d. sampled random latent vectors from a larger class of distributions. The critical conditions to maintain are the finite spectrum of $f$, boundedness of eigenfunctions, and local neighborhood properties. The local measure needs to be concentrated enough relative to the rate of change in the function $f$ so that when $n$ points are sampled from the space, there are sufficently many “nearby neighbors” for whom the function behaves similarly for any given point we would want to estimate. This primarily affects the nearest neighbor portion of the algorithm and analysis. also provides a formal discussion and results for extending the nearest neighbor analysis to accommodate settings beyond scalar Lipschitz functions. Our model can also be extended to asymmetric matrix settings and categorical data. Section 6 discuss how our theorem extend to some of these model variations.

3 Goal

The goal is to produce $\hat{F}$, an estimate of $F$, using observation matrix $M$ and knowledge of $\mathcal{E}$. We measure the estimation error through the maximum entry-wise error and the mean squared error. The maximum entry-wise error or $\infty$-norm of the error matrix $\hat{F}-F$ is defined as

We will provide bounds on this that hold with high probability, that is, with probability converging to $1$ as $n\to\infty$. The mean squared error (MSE) is defined as

In measuring error either with high probability or in expectation, the randomness is considered over the data generation process.

Algorithm

We propose and analyze a variation of the similarity based collaborative filtering algorithm. At its core, the collaborative filtering algorithm attempts to produce the estimate $\hat{F}(u,v)$ by averaging over observed entries $F(u^{\prime},v^{\prime})$ for a subset of tuples $(u^{\prime},v^{\prime})$ such that $u^{\prime}$ is “similar” to $u$ and $v^{\prime}$ is “similar” to $v$.

Noisy Nearest Neighbor Algorithm. We consider the following noisy nearest neighbor algorithm described below, followed by three different subroutines to compute distances depending on the sparsity regime of the dataset. (1) Compute distances $\hat{d}(u,v)$ between pairs of coordinates $u,v\in[n]^{2}$ using $M^{\prime}$ and $M^{\prime\prime}$. (2) For each $u,v\in[n]^{2}$, produce an estimate

where $\mathcal{E}^{\prime\prime\prime}_{uv}=\{(a,b)\in\mathcal{E}^{\prime\prime\prime}~{}:~{}\hat{d}(u,a)<\eta,\hat{d}(v,b)<\eta\}$ for some small enough $\eta>0$.

We will choose the threshold $\eta=\eta(n)$ depending on the local geometry of the latent feature space with respect to $\hat{d}(u,v)$, in order to guarantee that $\eta(n)$ is small enough to drive the bias to zero, yet large enough to ensure $|\mathcal{E}^{\prime\prime\prime}_{uv}|$ diverges so that the variance due to observation noise is small. The key part of the algorithm is determining how to estimate the distances $\hat{d}(u,v)$. In what follows, we describe three variations depending upon the observation density, $p$.

Dense Regime. When $p=\omega(n^{-\frac{1}{2}})$, it is feasible to compute distances by simply looking at the overlapping entries; this is popularly done in practice as well as analyzed theoretically in the recent works . For any $(u,a)\in[n]^{2}$,

where $\mathcal{O}_{ua}=\{y\in[n]:(u,y),(a,y)\in\mathcal{E}^{\prime}\}$. This is a finite sample approximation of $\int_{0}^{1}(f(\theta_{u},y)-f(\theta_{v},y))^{2}~{}dy$. When $p=\omega(n^{-\frac{1}{2}})$, it follows that $|\mathcal{O}_{ua}|=\omega(1)$ for all $u,a\in[n]^{2}$ with high probability, so that $\hat{d}(u,a)\approx\int_{0}^{1}(f(\theta_{u},y)-f(\theta_{v},y))^{2}~{}dy$. subsequently prove that for any Lipschitz latent function $f$ the MSE decays to as $n\to\infty$ as long as $p=\omega(n^{-\frac{1}{2}})$. The arguments of can be adapted to show that the maximum entry-wise error decays to with high probability as well. However, for $p=o(n^{-\frac{1}{2}})$, for most $u,a\in[n]^{2}$, $\mathcal{O}_{ua}=\emptyset$ with high probability and hence a different approach is needed – overcoming the sparse regime is the primary interest of this work.

Sparse Regime. Consider the sparse regime where $p=n^{-1+\kappa}$ for any $\kappa\in(0,\frac{1}{2})$; in this regime the overlap is small and thus new distance estimates are required. Recall that the function $f$ has finite spectrum, i.e. $f(\theta_{u},\theta_{v})=\sum_{k}\lambda_{k=1}^{r}q_{k}(\theta_{u})q_{k}(\theta_{v})$. We propose an estimator which approximates $d(u,v)=\|\Lambda^{t}Q(e_{u}-e_{v})\|_{2}^{2}$ by comparing depth $t$ neighborhoods of $u$ and $v$ in the data graph ${\cal G}=([n],\mathcal{E}^{\prime})$. Specifically, let the weight of an edge $(a,b)\in\mathcal{E}^{\prime}$ in graph $\mathcal{G}$ be the observed value $M(a,b)$ ($=M^{\prime}(a,b)$). By assumption, in expectation this weight equals $F(a,b)=f(\theta_{a},\theta_{b})$. Therefore, the product of weights along a path from $u$ to $y$, of length $t$, denoted as $(u,x_{1},\dots,x_{t-1},y)$ with $(u,x_{1}),(x_{1},x_{2}),\dots,(x_{t-1},y)\in\mathcal{E}^{\prime}$, in expectation equals

Therefore, the product of weights along the path connecting $u$ to $y$ is a good proxy of quantity $e_{u}^{T}Q^{T}\Lambda^{t}Qe_{y}$. Recall that each entry is observed independently with probability $p$ due to our assumed Bernoulli sampling model. Therefore, for any $u\in[n]$, the number of neighbors of $u$ in $\mathcal{G}$ scale as $pn=n^{\kappa}$. More generally, for $1\leq t\leq 1/\kappa$, the number of nodes at distance $t$ from $u$ scale as $n^{\kappa t}$. We choose $t$ large enough to guarantee that for any two nodes $u$ and $v$, there is a sufficient overlap between the two subset of nodes at distance $y$ from nodes $u$ and $v$ respectively. This suggests that we choose $t$ so that $n^{\kappa t}\approx n^{\frac{1}{2}}$, which in effect aggregates enough data in the sparse regime to match the expected number of observations per row in the dense regime. We formalize this intuition in the following construction of the distance estimates.

Let $\mathcal{S}_{u,s}$ denote the set of vertices which are at distance $s$ from vertex $u$ in the graph defined by edge set $\mathcal{E}^{\prime}$. Specifically, $i\in\mathcal{S}_{u,s}$ if the shortest path in $\mathcal{G}=([n],\mathcal{E}^{\prime})$ from $u$ to $i$ has a length of $s$. Let $\mathcal{T}_{u}$ denote a breadth-first tree in $\mathcal{G}$ rooted at vertex $u$. The breadth-first property ensures that the length of the path from $u$ to $i$ within $\mathcal{T}_{u}$ is equal to the length of the shortest path from $u$ to $i$ in $\mathcal{G}$. Let $\mathcal{T}_{u}^{t}\subset\mathcal{T}_{u}$ denote the sub-tree containing all nodes and edges in $\mathcal{T}_{u}$ up to and including depth $t$. If there is more than one valid breadth-first tree rooted at $u$, choose one uniformly at random. Let $N_{u,t}\in^{n}$ denote the following vector with support on the boundary of the depth-$t$ neighborhood of vertex $u$ (we also call $N_{u,t}$ the neighborhood boundary):

2 Reducing computation by subsampling vertices

The pairwise distances can only be estimated up to a limited precision depending on the sparsity of the data and amount of noise in the observations, and furthermore we tune the nearest neighbor threshold to tradeoff between bias and variance. As a result, the performance of the algorithm can be maintained with reduced computation by clustering the coordinates so that not all $n^{2}$ pairwise distances need to be computed. This would involve adding an extra step at the beginning of the algorithm that samples sufficiently many “anchor” vertices $\mathcal{K}\subset[n]$ that cover the space well. $|\mathcal{K}|$ should be chosen large enough such that for any vertex $u\in[n]$, there exists some anchor vertex $i\in\mathcal{K}$ which is “close” to $u$ in the sense that $\|\Lambda Q(e_{u}-e_{i})\|_{2}^{2}$ is small. For all $n$ vertices, we only compute the distances to each of the $|\mathcal{K}|$ anchor vertices, and we let $\pi:[n]\to\mathcal{K}$ be a mapping from each vertex to the anchor vertex that minimizes the estimated distance $\hat{d}$ as computed in the original algorithm statement, $\pi(u)=\operatorname*{arg\,min}_{i\in\mathcal{K}}\hat{d}(u,i)$. The final estimate then is given by

where $\mathcal{E}_{\pi(u)\pi(v)}$ denotes the set of undirected edges $(a,b)$ such that $(a,b)\in\mathcal{E}_{3}$ and both $\hat{d}(\pi(u),a)$ and $\hat{d}(\pi(v),b)$ are less than some threshold $\eta$. We can compute $\mathcal{E}_{\pi(u)\pi(v)}$ by the clustering assignments and distances of all vertices to the anchor vertices.

3 Computational Complexity

To analyze the computational complexity of the algorithm, we consider each step. Growing local neighborhoods around each vertex costs at most $n|\mathcal{E}|$, since there are $n$ vertices and the BFS trees visit each edge at most once. Computing the inner product for all pairs of vertices given the local neighborhood vectors costs at most $n^{2}|\mathcal{E}|$, since there are $n^{2}$ vertex pairs and $|\mathcal{E}|$ entries in the data matrix $M$. The final nearest neighbor estimator involves a (weighted) average of the datapoints, which costs at most $n^{2}|\mathcal{E}|$, as there are $n^{2}$ entries in the matrix to estimate, and at worst the estimate would involve averaging over $|\mathcal{E}|$ datapoints. This extremely crude bound leads to a computational complexity of $O(pn^{4})$. The bottleneck of the algorithm is the final nearest neighbor estimate, which may be reduced by using approximate nearest neighbor methods.

If we instead used the modified algorithm that subsamples $|\mathcal{K}|$ anchor vertices at random and treats them as “cluster centers”, there are only $(|\mathcal{K}|^{2}+n|\mathcal{K}|)$ pairwise distances computed, for a computational cost of $(|\mathcal{K}|^{2}+n|\mathcal{K}|)|\mathcal{E}|$ instead of $n^{2}|\mathcal{E}|$. Once we cluster the vertices, the final estimate is only computed for the pairwise cluster blocks, as the final estimate is a block constant matrix with only $|\mathcal{K}|^{2}$ distinct valued estimates. This results in $|\mathcal{K}|^{2}|\mathcal{E}|$ computation for the final step of the estimation. The computational complexity reduces from $O(n^{2}|\mathcal{E}|)$ to $O((|\mathcal{K}|^{2}+n|\mathcal{K}|)|\mathcal{E}|)$. The choice of $|\mathcal{K}|$ depends on the distribution of latent variables, the shape of the latent function, and the error tolerance. In a setting with finitely many latent types, then $|\mathcal{K}|$ would be roughly linear in the number of latent types.

A practical benefit of our algorithm is that it is amenable to a distributed and parallelized implementation. The key computational step of our algorithm involves comparing the expanded local neighborhoods of pairs of vertices to find the “nearest neighbors”. As the algorithm is inherently local with respect to the data graph, it can be easily implemented for large scale datasets where the data may be stored in a distributed fashion optimized for local graph computations. The local neighborhoods can be computed in parallel, as they are independent computations. Using approximate nearest neighbor techniques and subsampling vertices to cluster will additionally reduce the computation.

4 Discussion

In practice, we may not know the model parameters, and we would use cross validation to tune the BFS tree depth $t$ and nearest neighbor threshold $\eta$. If the depth $t$ is either too small or too large, then the vector $N_{u,t}$ will be too sparse, and will not optimally aggregate the datapoints. The threshold $\eta$ trades off between bias and variance of the final estimate. When the sampled observations are not uniform across entries, the algorithm may require more modifications to properly normalize for high degree hub vertices, as the optimal choice of depth $t$ may differ depending on the local sparsity.

In our algorithm, we assumed that we observed the edge set $\mathcal{E}$. Specifically, this means that we are able to distinguish between entries of the matrix that have value zero because they are not observed, i.e. $(i,j)\notin\mathcal{E}$, or if the entry was observed to be value zero, i.e. $(i,j)\in\mathcal{E}$ and $M(i,j)=Z(i,j)=0$. This fits well for applications such as recommendations, where the system does know the information of which entries are observed or not. Some social network applications contain this information (e.g. facebook would know if they have recommended a link which was then ignored) but other network information may lack this information, e.g. we do not know if link does not exist because observations are sparse, or because observations are dense but the probability of an edge is small. The absence of this knowledge would primarily affect the normalization of the neighborhood vectors as well as the normalization in the final averaging step.

The idea of comparing vertices by looking at larger radius neighborhoods was introduced in , and has connections to belief propagation and the non-backtracking operator . The non-backtracking operator was introduced to overcome the issue of sparsity. For sparse graphs, vertices with high-degree dominate the spectrum, such that the informative components of the spectrum get hidden behind the high degree vertices. The non-backtracking operator avoids paths that immediately return to the previously visited vertex in a similar manner as belief propagation, and its spectrum has been shown to be more well-behaved, perhaps adjusting for the high degree vertices, which get visited very often by paths in the graph. In our algorithm, the neighborhood paths are defined by first selecting a rooted tree at each vertex, thus enforcing that each vertex along a path in the tree is unique. This is important in our analysis, as it guarantees that the distribution of vertices at the boundary of each subsequent depth of the neighborhood is unbiased, since the sampled vertices are freshly visited.

Results

In all of the results below, we assume the latent variable model assumptions laid out in Section 2. As a reminder, we assume uniform Bernoulli sampling with density $p$, independent bounded observation noise, and a generative latent variable model where coordinates are associated to i.i.d. sampled latent variables and the underlying matrix behaves according to a bounded latent function $f$ that is Lipschitz and low rank (or approximately low rank) with bounded eigenfunctions.

We first provide theoretical bounds for the estimation error in both sparse regimes mentioned above when $f$ has finite spectrum with rank $r$.

Sparse Regime. Theorem 4.1 shows that the maximum entrywise error of the collaborative filtering algorithm using distance function (7) converges to zero in the sparse regime when $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$.

Let $f$ have rank $r$, $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$ so that $1/\kappa$ is not an integer. Consider the estimates produced by the nearest neighbor algorithm using the distance defined in (7) for $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor$ and selecting the nearest neighbor distance threshold to satisfy $\eta=\Theta(n^{-\frac{1}{2}(\kappa-\rho)})$ for any $\rho\in(0,\kappa)$. Let $C_{f}=|\lambda_{1}|/|\lambda_{r}|$ denote the condition number of the latent function $f$. With probability $1-o(1)$,

Sparser Regime. Theorem 4.2 shows that the maximum entrywise error of the collaborative filtering algorithm using distance function (8) converges to zero in the sparser regime when $p=n^{-1}\ln^{1+\kappa}n$ for some $\kappa>0$.

Let $f$ have rank $r$, $p=n^{-1}\ln^{1+\kappa}n$ for some $\kappa>0$. Consider the estimates produced by the nearest neighbor algorithm using the distance defined in (8) for $t=\lceil\frac{\ln(0.08/p)}{\ln(0.275np)}-r^{\prime}\rceil$ and selecting the nearest neighbor distance threshold to satisfy \eta=\Theta\Big{(}(\ln n)^{-\frac{1}{2}(\kappa-\rho)}\Big{)} for any $\rho\in(0,\kappa)$. With probability $1-o(1)$,

Theorems 4.1 and 4.2 show that for symmetric sparse matrix estimation, as long as the fraction of entries observed at random scale as $\frac{\log^{1+\kappa}(n)}{n}$ for any fixed $\kappa>0$, the estimation error of our proposed iterative variant of the classical collaborative filtering algorithm with respect to the $\max$-norm decays to as $n\to\infty$ assuming the underlying matrix of interest has constant rank $r$.

2 f𝑓f has ε𝜀\varepsilon-approximate rank r𝑟r

We extend the above stated result to the setting when the latent function $f$ has $\varepsilon$-approximate rank $r$; this captures settings where $f$ may have infinite but quickly decaying spectrum. We formally state the extension in the sparse regime ($p=n^{-1+\kappa}$), but we believe that a similar result is likely to hold for the sparser regime ($p=n^{-1}\log^{1+\kappa}(n)$) as well, which we omit for simplicity of presentation.

Let $f$ have $\varepsilon$-approximate rank $r$ for some $\varepsilon>0$, $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$ so that $1/\kappa$ is not an integer. Consider the estimates produced by the nearest neighbor algorithm using the distance defined in (7) for $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor<\frac{1}{\kappa}-1$ and selecting the nearest neighbor distance threshold to satisfy $\eta=\Theta(n^{-\frac{1}{2}(\kappa-\rho)})$ for any $\rho\in(0,\kappa)$. Let $C_{f,r}=|\lambda_{1}|/|\lambda_{r}|$ denote the condition number of the rank $r$ approximation to the latent function $f$. With probability $1-o(1)$,

As we assume the function values are bounded in $$, we can assume that$\varepsilon\in$, such that the dominating terms in (13) are$O\Big{(}rC_{f,r}^{1/\kappa}n^{-\frac{1}{4}(\kappa-\rho)}\Big{)}+O\Big{(}\sqrt{\varepsilon r|\lambda_{r}|^{-\frac{2}{\kappa}}\kappa^{-1}}\Big{)}$, and the dominating terms in (14) are$O\Big{(}r^{2}C_{f,r}^{2/\kappa}n^{-\frac{1}{2}(\kappa-\rho)}\Big{)}+O\Big{(}\varepsilon r|\lambda_{r}|^{-\frac{2}{\kappa}}\kappa^{-1}\Big{)}$. While we assume the function$f$is fixed with respect to$n$, when the function$f$has infinite spectrum, we can choose$\varepsilon$to decrease with$n$in order to tradeoff between the two terms in the error bound. Note that the approximate rank$r$and the approximate condition number$C_{f,r}$also depend on the choice of$\varepsilon$. In particular the relationship between$\varepsilon$,$r$, and$C_{f,r}$will depend on the spectrum of$f$and how quickly the tail decays to zero. Choosing a larger value of$r$will increase the condition number as$|\lambda_{r}|$will be smaller, and it will decrease the approximation error$\varepsilon$. Below we present a specific example as a consequence of Theorem 4.3.

Let $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$ so that $1/\kappa$ is not an integer. Consider $f$ such that for any $r\geq 1$, it has $\varepsilon_{r}$-approximate rank $r$ with $|\lambda_{r}|$ corresponding to rank $r$ approximation with $C_{f,r}=|\lambda_{1}|/|\lambda_{r}|$ being the condition number such that $|\lambda_{1}|=O(1)$ and

Then, for any $\delta>0$, for all $n$ large enough, with probability $1-o(1)$, \|\hat{F}-F\|_{\max}=O\Big{(}\sqrt{\delta}\Big{)}. Further, \text{MSE}(\hat{F})=O\Big{(}\delta\Big{)}.

For any $\delta>0$, by (15), there exists large enough $r=r(\delta)$ such that $\varepsilon_{r}|\lambda_{r}|^{-2/\kappa}r\leq\delta$. Due to $|\lambda_{1}|=O(1)$, $C_{f,r}=O(|\lambda_{r}|^{-1})$. Given choice of $r=r(\delta)$, for $n$ large enough we have $rC_{f,r}^{1/\kappa}n^{-\frac{1}{4}(\kappa-\rho)}\leq\sqrt{\delta}$. By (13) of Theorem 4.3 it follows that $\|\hat{F}-F\|_{\max}=O(\sqrt{\delta})$ with probability at least $1-o(1)$. By (14) of Theorem 4.3, it follows that $\text{MSE}(\hat{F})=O(\delta)$. $\square$

From Corollary 4.4, it follows that $\|\hat{F}-F\|_{\max}=o(1)$ with probability $1-o(1)$ and $\text{MSE}(\hat{F})=o(1)$ when the spectrum decays in such a way that $\lim_{r\to\infty}\varepsilon_{r}|\lambda_{r}|^{-2/\kappa}r=0.$

3 Discussion

In our latent variable model, the latent function $f$ is fixed with respect to $n$, so the max norm of the truth matrix is constant $\|F\|_{\max}=\Theta(1)$, and the Frobenius norm of the truth matrix scales linear with the matrix dimension so that $\frac{1}{n^{2}}\|F\|_{Fr}^{2}=\Theta(1)$. As a result the above stated results also show the convergence rates with respect to the relative errors of the max norm and normalized Frobenius norm.

The overall proof sketch can be split into two parts. First we prove that the estimated pairwise distances concentrate to a metric computed with respect to the true latent function $f$. Second we prove that given well behaved estimated distances, the nearest neighbor estimate with properly chosen thresholds to balance mean and variance will converge at the above stated rate. This second part of the proof is straightforward and follows the standard proof for any nearest neighbor style algorithm. The crux of the proof is arguing that in sparse settings the computed distances concentrate well. This relies on the uniform sampling assumption, independence of the observation noise, regularity of the latent feature variables, and the finite spectrum assumption of the latent function. The assumptions on the specific distribution of the latent variables and the Lipschitzness of the latent function are in fact primarily used for the second nearest neighbor portion of the proof, and thus can be relaxed. The key property needed is that there are sufficiently many “nearest neighbor” coordinates; the precise distribution of the latent variables and shape of the latent function will affect the tuning of the threshold parameter to tradeoff between bias and variance. We provide formal statements for a few variations of the model in Section 6.

In addition to providing bounds on the MSE, our theorem also provides bounds on the maximum entrywise error of the estimate. The rate of our maximum entrywise error is the square root of the MSE rate, which suggests that the error is uniformly spread across all entries. This is a stronger guarantee that the typical MSE bounds found in the literature, and it can be useful for downstream results that use the estimates for decision making such as ranking and recommendations.

Thus far, we have focused on finding conditions on $p$ that allow for consistent estimation with respect to both the MSE and max entrywise error. Our results also provide the rate at which the error decays. Specifically, our bound for the mean squared error (MSE) scales as $O((pn)^{-1/2+\rho})$ for any arbitrarily small constant $\rho>0$, and our bound for the max entrywise error is $O((pn)^{-1/4+\rho})$ for any small $\rho$.

Proof Sketch for Analyzing Noisy Nearest Neighbors

As the algorithm uses a fixed radius nearest neighbor estimate, the analysis boils down to arguing that the distance functions as defined in (7) and (8) have certain desired properties that enable the classical nearest neighbor algorithm to be effective. In this section we characterize the needed properties for the convergence of noisy nearest neighbors.

Our algorithm estimates $F(u,v)$, i.e. $f(\theta_{u},\theta_{v})$, according to (4), which simply averages over datapoints $M(u^{\prime},v^{\prime})$ corresponding to tuples $(u^{\prime},v^{\prime})$ for which $u^{\prime}$ is close to $u$ and $v^{\prime}$ is close to $v$ according to the estimated distance function $\hat{d}$. This simple nearest neighbor averaging estimator suggests that the last step of the analysis involves choosing the threshold $\eta$ to tradeoff between bias and variance.

Property 5.1 follows from choosing an appropriate ideal distance function $d$. In particular we will choose $d$ with respect to the spectral representation of $f$, and the desired property and the expression for $\text{{bias}}(\eta)$ will follow from the low rank assumption as well as the boundedness of the eigenfunctions.

Showing property 5.2 is the crux of the proof and follows from the design of the algorithm along with the assumptions of uniform sampling and the latent variable model. It essentially uses all the model assumptions except for Lipschitzness of $f$.

Property 5.3 is only used for the final step of the nearest neighbor analysis. In particular, as the estimate averages datapoints within an estimated nearby region of the target coordinates, there is a bias variance tradeoff that depends on how the datapoints are locally distributed. In particular, we need to guarantee that for any $(a,b)\in[n]^{2}$, there exists sufficiently many observed pairs $(u,v)\in[n]^{2}$ such that the function behaves similarly, i.e. $f(a,b)$ is close to $f(u,v)$. This property follows from our assumption that the latent variables are sampled i.i.d. from $U$, and that the function $f$ is $L$-Lipschitz. As discussed in section 2, these assumptions can be relaxed, but alternative assumptions would need to guarantee property 5.3 for some reasonable local measure function $\text{{meas}}(\eta)$.

Given the above three properties, we can then prove Lemma 5.1, which characterizes the error of the noisy nearest neighbor algorithm as a function of the bias function, meas function, and estimation error $\Delta$. Section 8 uses Lemma 5.1 to establish Theorems 4.1, 4.2, and 4.3 by simply showing the three properties for suitable choices of $\text{{bias}},\text{{meas}},$ and $\Delta$, and tuning $\eta$ accordingly to balance between different terms of the error. Proving that the distance estimates concentrate well, i.e. property 5.2, is the most involved part of the analysis, which we defer to sections 9 and 10. Property 5.1 follows from the low rank assumption and property 5.3 arises from the latent variable model assumptions, in particular the distribution of the latent variables and shape of the latent function.

Assume that properties 5.1-5.3 hold with probability $1-\alpha$ for some $\eta,\Delta,$ and $\eta^{\prime}=\eta-\Delta$; in particular $d$ is a bias-good distance function, $\hat{d}$ as estimated from $M^{\prime}$ and $M^{\prime\prime}$ is a $\Delta$-good distance estimate for $d$, and $\{\theta_{u}\}_{u\in[n]}$ is meas-represented. The noisy nearest neighbor estimate $\hat{F}$ computed according to (4) satisfies

for any $\delta\in(0,1)$. Furthermore, for any $\delta^{\prime}\in(0,1)$,

Inequality $(a)$ follows from Properties 5.1-5.2: $|d(u,a)-\hat{d}(u,a)|\leq\Delta$ and $\hat{d}(u,a)\leq\eta\implies d(u,a)\leq\eta+\Delta$. By definition $M(a,b)\in$ for all $(a,b)$, which implies $\text{Var}[M(a,b)]\leq 1$ for all $(a,b)\in\mathcal{E}^{\prime\prime\prime}$. Define ${\mathcal{V}}_{uv}=\{(a,b)\in[n]^{2}~{}:{d}(u,a)<\eta-\Delta,~{}{d}(v,b)<\eta-\Delta\}$. Assuming property 5.3,

By the Bernoulli sampling model and sample splitting process, each tuple $(a,b)\in[n]^{2}$ belongs to $\mathcal{E}^{\prime\prime\prime}$ with probability $p/2$ independently. By a straightforward application of Chernoff’s bound, it follows that for any $\delta\in(0,1)$,

Therefore, by assuming property 5.2, it follows that with probability at least $1-\exp\left(-\frac{\delta^{2}p\left(\text{{meas}}(\eta-\Delta)n\right)^{2}}{4}\right)$,

We add an additional $\alpha$ in the final MSE bound to account for the probability that properties 5.1-5.3 are violated.

This completes the proof of Lemma 5.1. $\square$

Extensions

As mentioned in Section 3.2, we can reduce the computational complexity of the algorithm by subsampling a set of anchor vertices $\mathcal{K}$ and only computing pairwise distances relative to the anchor vertices, equivalent to computing a clustering amongst vertices and using that to estimate. For pairs of anchor vertices $(a,b)\in\mathcal{K}^{2}$ which we also refer to as cluster centers, the algorithm estimates $\hat{F}(a,b)$ according to the original stated algorithm with no modifications. For $u\notin\mathcal{K}$, we denote $\pi(u)=\operatorname*{arg\,min}_{i\in\mathcal{K}}\hat{d}(u,i)$ to be a clustering that maps from $u$ to the closest anchor vertex in $\mathcal{K}$. The final estimate for $(u,v)\notin\mathcal{K}^{2}$ is then given by the estimate of the associated anchor vertices, which act as cluster centers, $\hat{F}(u,v)=\hat{F}(\pi(u),\pi(v))$.

The original argument provides high probability bounds on $|\hat{F}(u,v)-F(u,v)|$ for cluster centers $(u,b)\in\mathcal{K}^{2}$, as nothing changed in the algorithm for the cluster centers. The only additional part of the proof is to bound the additional bias for non cluster centers, as $|\hat{F}(u,v)-F(u,v)|\leq|\hat{F}(\pi(u),\pi(v))-F(\pi(u),\pi(v))|+|F(\pi(u),\pi(v))-F(u,v)|$. The first term is directly bounded by the current analysis, and the bias from the second term will depend on the size of $|\mathcal{K}|$. Recall our latent variable model assumption that each vertex $u$ is associated to a latent variable $\theta_{u}\sum U$ such that $F(u,v)=f(\theta_{u},\theta_{v})$ and $f$ is $L$-Lipschitz with respect to the latent variables. For $|\mathcal{K}|=\frac{2}{\delta}\log(\frac{1}{\delta})$, with probability at least $1-\delta$, each interval $[(i-1)\delta,i\delta]$ for $i\in[1/\delta]$ contains at least one anchor point in $\mathcal{K}$, as the latent variables of these anchor points are chosen at random. Under this good event, then $\max_{u\in[n]}\min_{i\in\mathcal{K}}|\theta_{u}-\theta_{i}|\leq\delta$.

We discuss the results and analysis for the sparse setting when $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$, however a similar argument applies for the sparser setting of $p=n^{-1}\ln^{1+\kappa}n$ as well. Equation (20) will show that $d(\theta_{u},\theta_{v})\leq|\lambda_{1}|^{2t}L^{2}|\theta_{u}-\theta_{v}|^{2}$, so that for some $u\in[n]$, the closest anchor point $a\in\mathcal{K}$ with respect to the latent representation will also satisfy $d(\theta_{u},\theta_{a})\leq|\lambda_{1}|^{2t}L^{2}\delta^{2}$. As Property 5.2 guarantees $|d(\theta_{u},\theta_{a})-\hat{d}(u,a)|\leq\Delta$ for all estimated distances, it follows that $d(\theta_{u},\theta_{\pi(u)})\leq|\lambda_{1}|^{2t}L^{2}\delta^{2}+2\Delta$ for all $u\in[n]$. By Property 5.1, $|F(\pi(u),\pi(v))-F(u,v)|\leq\text{{bias}}(|\lambda_{1}|^{2t}L^{2}\delta^{2}+2\Delta)$. We choose $|\mathcal{K}|$ so that $\delta=\frac{\sqrt{\Delta}}{L|\lambda_{1}|^{t}}$, and we plug in the choice of $\Delta$ and $t$ from Theorem 4.1, resulting in $\delta=Br|\lambda_{1}|^{(\kappa+1)/\kappa}L^{-1}n^{-\frac{1}{4}(\kappa-\rho)}=o(1)$ so that $|\mathcal{K}|=\Theta(n^{\frac{1}{4}(\kappa-\rho)})$. This choice of $|\mathcal{K}|$ will guarantee that the extra added bias does not change the existing guarantees in Theorem 4.1 by more than a constant.

2 Local Geometry

We can generalize the latent variable model beyond scalar valued latent variables and Lipschitz latent functions. These assumptions only affect the function meas in Property 5.3, and thus it only changes the last portion of the nearest neighbor proof in which we tune the threshold $\eta$ to tradeoff between the bias and variance terms. We present two examples of extending our results to a different local geometry, illustrating the modifications for the sparse setting when $p=n^{-1+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$.

Next we discuss a higher dimensional setting. Assume the latent variables are sampled uniformly over a $m$-dimensional hypercube such that $\theta_{u}\sim U(^{m})$ and the latent function $f$ is $L$-Lipschitz with respect to an underlying metric $d_{m}$, such that the measure of a ball with radius $\delta$ is $\Theta(\delta^{m})$. Property 5.3 would instead hold for $\text{{meas}}(\eta^{\prime})=\Theta((\frac{\sqrt{\eta^{\prime}}}{\lambda^{t}L})^{m})$, resulting in a different choice of threshold $\eta$ to balance between bias and variance. If $m\leq(\kappa+2)/\kappa$, then the current $\text{{bias}}(\Delta)$ term dominates such that we would choose $\eta=\Theta(\Delta)$, and the error convergence rate will be the same as that stated in Theorem 4.1. For high dimension $m>(\kappa+2)/\kappa$, we choose the threshold $\eta=\Theta((pn^{2})^{-1/(m+1)})$ such that the MSE bound will scale as $\Theta((pn^{2})^{-1/(m+1)})=\Theta(n^{-(1+\kappa)/(m+1)})$ and the max entrywise error bound will scale as $\Theta(n^{-(1+\kappa)/2(m+1)})$.

3 Asymmetric Matrix

Even though our stated results are for symmetric models, we can transform an asymmetric latent variable model to a symmetric model as long as the row and column dimensions grow proportionally to one another. Consider an $n\times m$ matrix $F$ which we would like to learn, where $F(u,v)=f(\alpha_{u},\beta_{v})\in$, and $f$ has finite spectrum. We can construct a $(n+m)\times(n+m)$ matrix where $F$ is placed on the off-diagonal blocks and the diagonal $n\times n$ and $m\times m$ blocks are set to zero. We can argue that this constructed matrix is sampled form a symmetric latent model, so that we can apply our algorithm and analysis directly.

4 Categorical Valued Data

5 Non-Uniform Sampling

We assumed a uniform sampling model, where each entry is observed independently with probability $p$. However, in reality the probability that entries are observed may not be uniform across all pairs $(i,j)$. Our results can be extended to a setting where the sampling probability is instead a function of the latent variable, i.e. entry $(i,j)$ is observed with probability $c_{n}g(\theta_{i},\theta_{j})$ where $g$ is a Lipschitz low rank function independent of $n$ and $c_{n}$ is a scaling factor governing the density. The observed data $M(i,j)$ would then be sampled according to

We can essentially then apply our algorithm twice, first using data matrix $M$ to estimate the product $g(\theta_{i},\theta_{j})f(\theta_{i},\theta_{j})$ up to a scaling factor. Second we apply our algorithm to the binary adjacency matrix representing the sparsity of the observation set $\Omega$ in order to estimate $g(\theta_{i},\theta_{j})$ up to scaling factor. The one nuance one would have to handle is that since the set of observed entries is not uniformly sampled, the constructed BFS trees will grow non-uniformly, which will affect the normalization and scaling terms. As the model is only recoverable up to scaling, this is the best we can do. If we had data from a two-step sampling process in which we first observe binary edges sampled uniformly with probability $c_{n}$, and then subsequently observed datapoints sampled with an additional probability $g(\theta_{i},\theta_{j})$, then the model would exactly fall into our assumptions and the results could directly be applied to estimating $g(\theta_{i},\theta_{j})$ and the product $g(\theta_{i},\theta_{j})f(\theta_{i},\theta_{j})$.

Experiments

We show results on synthetic data to illustrate the performance of our algorithm. We did not do sample splitting as it is primarily introduced for the purpose of the analysis. We computed distances according to equation (7) (but again without sample splitting) for fixed radius parameters of $t\in\{0,1,2,3,4\}$. Note that the depth for expanding the BFS tree is until $t+1$. We did not specifically tune the nearest neighbor threshold $\eta$, but simply chose it to be the 70th percentile amongst all estimated distances. As a result, the expected number of entries used to compute the final weighted average estimate is $0.49pn^{2}$. We compare against a naive baseline which predicts using the column-wise mean. And we compare against the softimpute implementation in python’s fancyimpute package and alternating least squares with rank 2 from parafac algorithm in the python tensorly package (higher rank performed more poorly in the sparse setting as it overfit to noise). Nuclear norm minimization was too slow for the size of instances that we show and thus was omitted.

For a $\kappa\in(0,1]$, the density is chosen to be $p=n^{-1+\kappa}$, and each entry is observed (and thus included in sample set $\Omega$) with probability $p$ independently of all other entries. For each observed entry $(u,v)\in\Omega$, there is an added independent Gaussian noise $M(u,v)=F(u,v)+\varepsilon(u,v)$, where $\varepsilon(u,v)\sim N(0,\sigma^{2})$ where $\sigma$ is chosen to be the 40th percentile of the magnitude of entries in $F$. We show results for $n=500,1000,$ and $5000$.

We compute an adjusted mean squared error (MSE), limited to the error in predicting missing entries, and we normalize by the squared error of predicting with zeros. When the adjusted MSE is larger than 1, it means the estimate is worse than predicting all zeros.

Figure 1 shows the adjusted MSE of the algorithms with respect to the sampling probability $p$. When $p$ is very small, then our algorithm with the optimal choice of the depth parameter $t$ performs better than ALS and SoftImpute, however when it is too sparse than either the simple mean estimate or predicting with all zeros is best. Note that we did not do any tuning of the nearest neighbor parameter $\eta$, and thus there may be additional gains possible for our algorithm. If we consider the minimum density for which the algorithm performs better than the simple mean, SoftImpute requires the most dense observation. The minimum density required for our algorithm depends on optimally choosing the depth parameter $t$, but for an optimal choice of $t$, our algorithm requires less data than ALS before it performance better than the simple mean.

Figure 2 shows the adjusted MSE of the algorithms with respect to the exponent of the density parameter $\kappa$ where $p=n^{-1+\kappa}$. This rescales the $x$-axis so that the small values of $p$ are more visible. We plot only up to $\kappa=0.6$ as we are focusing on the sparse regime with little overlaps in entries between pairs of rows and columns. Notice the dependence on the performance of our algorithm with respect to the radius parameter $t$ as illustrated best in Figure LABEL:fig:5000_mse_kappa. For too small values of $t$ the alg is suboptimal as it does not aggregate data sufficiently, but for too large values of $t$ the algorithm again is suboptimal as it simply estimates zeros due to the BFS trees running out of vertices.

Figure 3 shows the time each of the algorithms took to run. We can see that our proposed algorithm is faster than SoftImpute, and this gap in speed is amplified with large $n$. Alternating Least Squares (ALS) is very fast, nearly as fast as the simple mean. Nuclear norm minimization was too slow to run on the size of instances in our example and thus was not included.

Proofs for Theorems 4.1, 4.2, and 4.3

In this section, we use the noisy nearest neighbor lemma 5.1 along with to establish Theorems 4.1, 4.2, and 4.3. Proofs of the concentration of distance estimates is deffered to sections 9 and 10.

We prove that as long as $p=n^{-1+\kappa}$ for any $\kappa\in(0,\frac{1}{2})$, with high probability, properties 5.1-5.3 hold for an appropriately chosen function $d$, and for distance estimates $\hat{d}$ computed according to (7) with $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor$. We subsequently use Lemma 5.1 to conclude Theorem 4.1. The most involved part in the proof is establishing that property 5.2 holds with high probability for an appropriately chosen $\Delta$, which is delegated to Lemma 8.1.

Good distance $d$ and Property 5.1. We start by defining the ideal distance $d$ as follows. For all $(u,v)\in[n]^{2}$, let

Recall that $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor$. Assuming $p=n^{-1+\kappa}$, $\kappa\in(0,\frac{1}{2})$

where (a) follows from assuming that $|q_{k}(\theta)|\leq B$ for all $k\in[r]$ and $\theta\in$. In summary, property 5.1 is satisfied for distance function $d$ defined according to (17) and $\text{{bias}}(\eta)=2B|\lambda_{r}|^{-t}\sqrt{r\eta}$.

Good distance estimate $\hat{d}$ and Property 5.2. We state the following Lemma when $f$ has rank $r$, whose proof is delegated to Section 9.

Let $f$ has rank $r$, $p=n^{-1+\kappa}$ for $\kappa\in(0,\frac{1}{2})$ such that $1/\kappa$ is not an integer. Consider $\hat{d}$ as computed in (7) with $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor$. For any $\rho\in(0,\kappa)$

with probability at least 1-O\Big{(}n^{2}\exp\big{(}-\Theta(n^{\min(\rho,\kappa(t-\frac{1}{2}))})\big{)}\Big{)}.

Lemma 8.1 implies that property 5.2 holds with probability $1-o(1)$ for some $\Delta=\Theta(Br|\lambda_{1}|^{2/\kappa}n^{-(\kappa-\rho)/2})$ and any $\rho\in(0,\kappa)$. The distance error bound $\Delta$ is minimized by choosing $\rho$ arbitrarily close to 0 so that $\Delta$ can be arbitrarily close to $\Theta(Br|\lambda_{1}|^{2/\kappa}n^{-\kappa/2})=\Theta(Br|\lambda_{1}|^{2/\kappa}(pn)^{-1/2})$.

The corresponding statement for $f$ that has $\varepsilon$-approximate rank $r$ is stated below.

Let $f$ have $\varepsilon$-approximate rank $r$, $p=n^{-1+\kappa}$ for $\kappa\in(0,\frac{1}{2})$ such that $1/\kappa$ is not an integer. Consider $\hat{d}$ as computed in (7) with $t=\lfloor\frac{\ln(1/p)}{\ln(np)}\rfloor$. For any $\rho\in(0,\kappa)$

with probability at least 1-O\Big{(}n^{2}\exp\big{(}-\Theta(n^{\min(\rho,\kappa(t-\frac{1}{2}))})\big{)}\Big{)}.

We assumed that the latent parameters $\{\theta_{u}\}_{u\in[n]}$ are sampled i.i.d. uniformly over $$. Therefore, for any$\theta_{u}\in$, for any$v\in[n]$and$\eta^{\prime}>0$,