Robust Max Entrywise Error Bounds for Tensor Estimation from Sparse Observations via Similarity Based Collaborative Filtering

TL;DR

This paper introduces a similarity-based collaborative filtering method for tensor estimation from sparse noisy observations, achieving near-optimal sample complexity and robustness to noise, with theoretical guarantees on error decay.

Contribution

The paper presents a novel tensor estimation algorithm that nearly matches the conjectured lower bound on sample complexity and demonstrates robustness to arbitrary bounded noise.

Findings

Achieves tensor recovery with error decaying as sample probability increases

Robust to arbitrary bounded noise, with error degrading polynomially

Sample complexity matches the connectivity threshold of the similarity graph

Abstract

Consider the task of estimating a 3-order $n \times n \times n$ tensor from noisy observations of randomly chosen entries in the sparse regime. We introduce a similarity based collaborative filtering algorithm for estimating a tensor from sparse observations and argue that it achieves sample complexity that nearly matches the conjectured computationally efficient lower bound on the sample complexity for the setting of low-rank tensors. Our algorithm uses the matrix obtained from the flattened tensor to compute similarity, and estimates the tensor entries using a nearest neighbor estimator. We prove that the algorithm recovers a finite rank tensor with maximum entry-wise error (MEE) and mean-squared-error (MSE) decaying to $0$ as long as each entry is observed independently with probability $p = \Omega(n^{-3/2 + \kappa})$ for any arbitrarily small $\kappa > 0$. More generally, we establish robustness of the estimator, showing that when arbitrary noise bounded by $\varepsilon \geq 0$ is added to each observation, the estimation error with respect to MEE and MSE degrades by $\text{poly}(\varepsilon)$. Consequently, even if the tensor may not have finite rank but can be approximated within $\varepsilon \geq 0$ by a finite rank tensor, then the estimation error converges to $\text{poly}(\varepsilon)$. Our analysis sheds insight into the conjectured sample complexity lower bound, showing that it matches the connectivity threshold of the graph used by our algorithm for estimating similarity between coordinates.

Equations

$$F=\sum _{k=1}^r{u}_k\otimes v_k\otimes w_k,$$

$$F=(Q_1\otimes \cdots Q_d)\cdot (\Lambda ):=\sum _{\mathbf{k}\in [r_1]\times [r_2]\cdots \times [r_d]}\Lambda (\mathbf{k})Q_1(\cdot ,k_1)\otimes Q_2(\cdot ,k_2)\cdots \otimes Q_d(\cdot ,k_d),$$

$$M(u,v,w)$$

$$\text{MSE}(\hat{F}):=\mathbb{E}\left[\frac{1}{n^3}{\sum }_{(u,v,w)\in [n{]}^3}(\hat{F}(u,v,w)-F(u,v,w){)}^2\right].$$

$$\Vert F-\hat{F}{\Vert }_{\max }:=\underset{(u,v,w)\in [n{]}^3}{\max }|\hat{F}(u,v,w)-F(u,v,w)|.$$

$$f(u,v,w)={\sum }_{k=1}^r{\lambda }_k{q}_k({\theta }_u)q_k({\theta }_v)q_k({\theta }_w),$$

$$F={\sum }_{k=1}^r{\lambda }_k(Q^T{e}_k)\otimes (Q^T{e}_k)\otimes (Q^T{e}_k).$$

$$\underset{{\theta }_u,{\theta }_v,{\theta }_w\in [0,1]}{\sup }|f({\theta }_u,{\theta }_v,{\theta }_w)-f_r({\theta }_u,{\theta }_v,{\theta }_w)|\le \mathbf{\epsilon }$$

$$F_r(u,v,w)=f_r({\theta }_u,{\theta }_v,{\theta }_w)={\sum }_{k=1}^r{\lambda }_k{q}_k({\theta }_u)q_k({\theta }_v)q_k({\theta }_w),$$

$$\Vert F-F_r{\Vert }_{\max }$$

$$\mathbb{E}[(M(a,b,c)-M(u,b,c){)}^2|{\theta }_a,{\theta }_u]$$

$$=\mathbb{E}[(F(a,b,c)-F(u,b,c){)}^2|{\theta }_a,{\theta }_u]+\mathbb{E}[{ϵ}_{abc}^2+{ϵ}_{ubc}^2]$$

$$=\mathbb{E}[({\sum }_k{\lambda }_k(q_k({\theta }_a)-q_k({\theta }_u))q_k({\theta }_b)q_k({\theta }_c){)}^2|{\theta }_a,{\theta }_u]+2{\sigma }^2$$

$$=\mathbb{E}[{\sum }_k{\lambda }_k^2(q_k({\theta }_a)-q_k({\theta }_u){)}^2{q}_k({\theta }_b{)}^2{q}_k({\theta }_c{)}^2|{\theta }_a,{\theta }_u]+2{\sigma }^2$$

$$={\sum }_k{\lambda }_k^2(q_k({\theta }_a)-q_k({\theta }_u){)}^2+{\sigma }^2$$

$$=\Vert \Lambda Q(e_a-e_u){\Vert }_2^2+2{\sigma }^2,$$

$$d({\theta }_u,{\theta }_a)$$

$$\mathbb{E}\left[M(u,e_1^1,e_2^1)\left(\prod _{i=1}^{t-1}M(e_1^i,e_2^i,x^i)M(x^i,e_1^{i+1},e_2^{i+1})\right)M(e_1^t,e_2^t,a)|{\theta }_u,{\theta }_a\right]$$

$$=\mathbb{E}\left[f({\theta }_u,{\theta }_{e_1^1},{\theta }_{e_2^1})\left(\prod _{i=1}^{t-1}f({\theta }_{e_1^i},{\theta }_{e_2^i},{\theta }_{x^i})f({\theta }_{x^i},{\theta }_{e_1^{i+1}},{\theta }_{e_2^{i+1}})\right)f({\theta }_{e_1^t},{\theta }_{e_2^t},{\theta }_a)|{\theta }_u,{\theta }_a\right]$$

$$=\sum _{k=1}^r{\lambda }_k^{2t}{q}_k({\theta }_u)q_k({\theta }_a)=e_u^T{Q}^T{\Lambda }^{2t}Q{e}_a.$$

$${\mathcal{B}}_{u,s}={\cup }_{h\in \lfloor s/2\rfloor }{\mathcal{S}}_{u,h}{\cup }_{l\in \lceil s/2\rceil }{\mathcal{U}}_{u,l}.$$

$${\mathcal{B}}_{u,s}={\cup }_{h\in \lfloor s/2\rfloor }{\mathcal{S}}_{u,h}\cup \left\{x\in [n]\text{s.t.}\exists (y,z)\in {\cup }_{l\in \lceil s/2\rceil }{\mathcal{U}}_{u,l}\text{ satisfying }x\in \{y,z\}\right\}.$$

$$N_{u,s}(a)$$

$$W_{u,s}(b,c)$$

$$t$$

$$\text{dist}(u,v)$$

$$=\frac{1}{|{\mathcal{V}}_B|p^2}\sum _{(\alpha ,\beta )\in {\mathcal{V}}_B}\sum _{a,b\in [n{]}^2}({\tilde{N}}_{u,t}(a)-{\tilde{N}}_{v,t}(a))M_B(a,\alpha ,\beta )M_B(b,\alpha ,\beta )({\tilde{N}}_{u,t}(b)-{\tilde{N}}_{v,t}(b))$$

$$\text{dist}(u,v)$$

$$Z_{uv}$$

$${\mathcal{V}}_B(u,v,t)$$

$$T_{uv}(\alpha ,\beta )$$

$$\hat{F}(u,v,w)$$

$$\eta$$

$$=O(n^{-(\kappa -\psi )})=O\left(\frac{n^{\psi }}{(p^2{n}^3{)}^{1/2}}\right),$$

$$\Vert F-\hat{F}{\Vert }_{\max }$$

$$\eta$$

$$=O(n^{-(\kappa -\psi )}+t\mathbf{\epsilon }(1+\mathbf{\epsilon }{)}^{2t-1}+t^2{\mathbf{\epsilon }}^2(1+\mathbf{\epsilon }{)}^{4t-2})$$

Full text

1 Introduction

Tensor estimation involves the task of predicting underlying structure in a high-dimensional tensor structured dataset given only a sparse subset of observations. We call this “tensor estimation” rather than the conventional “tensor completion” as the goal is not only to fill missing entries but also to estimate entries whose noisy observations are available. Whereas matrices represent data associated to two modes, rows and columns, tensors represent data associated to general $d$ modes. For example, a datapoint collected from a user-product interaction an e-commerce platform may be associated to a user, product, and date/time, which could be represented in a 3-order tensor where the three modes would correspond to users, products, and date/time. Image data is also naturally represented in a 3-order tensor format, with two modes representing the location of the pixel, and the third mode representing the RGB color components. Video data furthermore introduces a fourth mode indexing the time. Dynamic network data can also be represented in a tensor with one mode indexing the time and the other two modes indexing the nodes in the network.

There are many applications in which the dataset inherently has a lot of noise or is very sparsely observed. For example, e-commerce data is typically very sparse as the typical number of products a user interacts with is very small relative to the entire product catalog; furthermore the timepoints at which the user interacts with the platform may be sparse. When the dataset can be represented as a matrix, equivalent to a 2-order tensor, there has been a significant amount of research in designing practical algorithms and studying statistical limits for matrix estimation, a critical step in data pre-processing. Under conditions on uniform sampling and incoherence, the minimum sample complexity for estimation has been tightly characterized and achieved by simple algorithms. It is a natural and relevant question then to consider whether the techniques developed can extend to higher order tensors as well.

The previous literature has primarily focused on attaining consistency with respect to the mean squared error (MSE). Unfortunately as this is aggregated over the error in the full tensor, it does not translate to consistent bounds on entrywise error, as the error on a single entry could be very large despite the MSE being small due to averaging over many entries. However, entrywise bounds are important in practice as the results of tensor estimation are often used subsequently for decisions that involve comparisons between the estimates of individual entries.

In this work we focus on attaining consistent max entrywise error bounds by extending similarity based collaborative filtering algorithms to tensor estimation. Similarity based collaborative filtering is widely used in industry due to its simplicity, interpretability, and amenability to distributed and parallelized implementations. In the analysis of our proposed algorithm we show that it achieves a sample complexity that nearly matches a conjectured lower bound for computationally efficient algorithms. Perhaps most notably, our theoretical guarantees provide high probability bounds on the maximum entrywise error of the estimate, which is significantly stronger than the typical mean squared error style bounds found in the literature for other algorithms. We also provide error bounds under arbitrary bounded noise, which has implications towards approximately low rank settings.

1.1 Related Literature

Algorithms for analyzing sparse low rank matrices (equivalent to 2-order tensors) where the observations are sampled uniformly randomly have been well-studied. The algorithms consist of spectral decomposition or matrix factorization [1, 2, 3], nuclear norm minimization [4, 5, 6, 7, 8, 9], gradient descent [1, 2, 10, 11, 12], alternating minimization [13, 14], and nearest neighbor style collaborative filtering [15, 16, 17, 18, 19]. These algorithms have been shown to be provably consistent as long as the number of observations is $\Omega(rn\,{\sf poly}(\log n))$ for the noiseless setting where $r$ is the rank and $n$ is the number of rows and columns [1, 4]; similar results have been attained under additive Gaussian noise [2, 5] and generic bounded noise [3, 18]. Lower bounds show that $\Omega(rn)$ samples are necessary for consistent estimation, and $\Omega(rn\log(n))$ samples are necessary for exact recovery [5, 6], implying that the proposed algorithms are nearly sample efficient order-wise up to the information theoretic lower bounds.

There are results extending matrix estimation algorithms to higher order tensor estimation, assuming the tensor is low rank and that observations are sampled uniformly at random. The earliest approaches simply flatten or unfold the tensor to a matrix and subsequently apply matrix estimation algorithms [20, 21, 22, 23]. A $d$-order tensor where each dimension is length $n$ would be unfolded to a $n^{\lfloor d/2\rfloor}\times n^{\lceil d/2\rceil}$ matrix, resulting in a sample complexity of $O(n^{\lceil d/2\rceil}\,{\sf poly}(\log n))$, significantly larger than the natural statistical lower bound that is linear with $n$ due to the model being parameterized by linear in $n$ latent variables. When $d$ is odd, for example $d=3$ the sample complexity for this naive approach scales as $O(n^{2}\,{\sf poly}(\log n))$.

Subsequent works have improved upon this sample complexity, requiring only $\Omega(n^{3/2}\,{\sf poly}(\log n))$ observed entries for a 3-order tensor [24, 25, 26, 27, 28, 29, 30, 31]. [24, 25] analyzes the alternating minimization algorithm for exact recovery of the tensor given noiseless observations and finite rank $r=\Theta(1)$. [27, 28] use the sum of squares (SOS) method, and [29] introduces a spectral method. Both of these latter algorithms can handle noisy observations and overcomplete tensors where the rank is larger than the dimension. [30, 31] furthermore characterize the minimax optimal rate for the MSE and achieve it using spectral initialization followed by power iteration. For a general $d$-order tensor these results translate into a sample complexity scaling as $O(n^{d/2})$, improving upon $O(n^{\lceil d/2\rceil}$. [32] prove that tensor nuclear norm minimization can recover the underlying low-rank $d$-order tensor with $O(n^{3/2}\,{\sf poly}(\log n))$ samples in the noiseless setting; however, the algorithm is not efficiently computable as computing tensor nuclear norm is NP-hard [33].

[27] conjecture that any polynomial time estimator for a 3-order tensor must require $\Omega(n^{3/2})$ samples, based on a reduction between tensor estimation for a rank-1 tensor to the random 3-XOR distinguishability problem. They argue that if using the sum of squares hierarchy to construct relaxations for tensor rank, any result that achieves a consistent estimator with fewer than $n^{3/2}$ samples will violate a conjectured hardness of random 3-XOR distinguishability. Information theoretic bounds imply that one needs at least $\Omega(drn)$ observations to recover a $d$-order rank $r$ tensor, consistent with the degrees of freedom or number of parameters in the model. Interestingly, this implies a conjectured gap between the computational and statistically achievable sample complexities, highlighting how tensor estimation is distinctly more difficult than matrix estimation.

The majority of the results in tensor estimation provide bounds on the mean squared error, which aggregates errors across entries. In contrast our results will also provide bounds on the maximum entrywise error. There has been recent interest on developing matrix estimation methods that provide max entrywise bounds using a leave one out analysis, cf. [34], [35], [36], [37], [38] and [39]. Subseqeuntly [40, 41] extended the leave one out analysis to the tensor setting to obtain entrywise error bounds for a gradient descent algorithm with spectral initialization. The analysis and algorithm in our paper is significantly different than their work, as it results from showing high probability guarantees on the similarity computation, which is akin to a spectral algorithm, followed by a nearest neighbor analysis. Our results suggests that a combination of spectral analysis and nearest neighbor smoothing can achieve entrywise consistent estimates without further gradient descent refinements. As nearest neighbor methods are still widely used in industry, understanding their theoretical performance is of interest.

1.2 Contribution

Our results answer the following unresolved questions in the literature.

Is there a computationally efficient estimator that can provide a consistent estimation of low-rank tensor with respect to maximum entry-wise error (MEE) with minimal sample complexity of $\Omega(n^{\frac{3}{2}})$ in the presence of noise?

• 2.

Is there an extension of matrix estimation collaborative filtering algorithm for the setting of tensor estimation that can provide consistent estimation with such minimal sample complexity?

• 3.

Can the estimator be robust to adversarial bounded noise in the observations?

To begin with, we propose an algorithm for a symmetric $3$-order tensor estimation which generalizes a nearest neighbor collaborative filtering algorithm for sparse matrix estimation introduced in [18]. Naïvely applying the matrix estimation algorithm in [18] to the $n\times n^{2}$ matrix obtained by unfolding the $3$-order tensor would require $\Omega(n^{2})$ samples, far more than the desired sample complexity of $\Omega(n^{\frac{3}{2}})$. However, we argue that such a matrix obtained from the unfolded tensor can be used, after non-trivial modification, to compute the similarities between rows accurately using $\Omega(n^{\frac{3}{2}+\kappa})$ samples for any positive $\kappa>0$. After computing these similarities we can achieve consistent estimation via a nearest neighbor estimator by additionally using the tensor structure.

Specifically, we establish that the mean squared error (MSE) in the estimation converges to [math] as long as $\Omega(n^{3/2+\kappa})$ random samples are observed for any $\kappa>0$ for tensor with rank $r=\Theta(1)$. We further establish a stronger guarantee that the maximum entry-wise error (MEE) converge to [math] with high probability with similar sample complexity of $\Omega(n^{3/2+\kappa})$. Thus, this simple iterative collaborative filtering algorithm nearly achieves the conjectured computational sample complexity lower bound of $\Omega(n^{3/2})$ for tensor estimation. While we present the results for symmetric tensors, our method and analysis can extend to asymmetric tensors, which we discuss in Section 5.4.

Beyond low-rank tensors, our results hold for tensors with potentially countably infinite rank as long as they can be well approximated by a low-rank tensor. Specifically, if the tensor can be approximated with $\boldsymbol{\varepsilon}\geq 0$ with respect to max-norm by a rank $r=\Theta(1)$ tensor, then the MSE converges to ${\sf poly}(\boldsymbol{\varepsilon})$ and MEE converges to ${\sf poly}(\boldsymbol{\varepsilon})$ with high probability as long as $\Omega(n^{3/2+\kappa})$ random samples are observed for any $\kappa>0$. This follows as a consequence of the robustness property of the algorithm that we establish: if arbitrary noise bounded by $\boldsymbol{\varepsilon}\geq 0$ is added to each observation, then the estimation error with respect to MEE and MSE degrades by ${\sf poly}(\boldsymbol{\varepsilon})$.

To establish our results, the key analytic tool is utilizing certain concentration properties of a bilinear form arising from the local neighborhood expansion of any given coordinate for an asymmetric matrix with dimensions $n\times n^{2}$. This generalizes the analysis of a similar property for symmetric matrices in the prior work of [18]. Specifically, establishing the desired concentration requires handling dependencies arising in the local neighborhood expansion of the $3$-order tensor that was absent in the matrix setting considered in [18]. Subsequently, we require a novel analytic method compared to the prior work. In particular we believe that the proof techniques in Lemma 7.7 may be useful to other settings in which one may desire a tighter concentration on sums of sparse random variables. As a consequence, we also establish performance guarantees for matrix estimation for asymmetric matrices having dimensions of different order, generalizing beyond of [18].

The algorithm and analysis also sheds insight on the conjectured lower bound for $3$-order tensor. In particular, the threshold of $n^{3/2}$ is precisely the density of observations needed for the connectivity in the associated graph that is utilized to calculate similarities. If the graph is disconnected, the similarities can not be computed, while if the graph is connected, we are able to show that similarity calculations yield an excellent estimator. Understanding this relationship further remains an interesting open research direction.

A benefit of our algorithm is that it can be implemented in a parallelized manner where the similarities between pair of indices are computed in parallel. This lends itself to a distributed, scalable implementation. A naive bound on the computational complexity of our algorithm for 3-order tensor is at most $pn^{6}$. As discussed in Section 5.4, with use of approximate nearest neighbors, these can be further improved and made truly implementable. 333A weaker abbreviated version of this result appeared in [42] without any proofs or discussion. Since the preliminary results, the convergence rates of the error have improved, and we have new results showing a perturbation analysis under arbitrary bounded noise, which extends our results to the approximately low rank setting. We also present a modification of the algorithm that significantly improves the overall computational complexity of the algorithm.

2 Preliminaries

Tensor estimation from sparse observations hinges on an assumption that the true model exhibits low dimensional structure despite the high dimensional representation. However, there is not a unique definition of rank in the tensor setting, as natural generalizations of matrix rank lead to different quantities when extended to higher order tensors. We will focus on two commonly used definitions of tensor rank, the CP rank and the Tucker or multilinear rank.

For a $d$-order tensor $F\in\mathbb{R}^{n^{d}}$, we can decompose $F$ into a sum of rank-1 tensors. For example if $d=3$, then

[TABLE]

where $\{u_{k},v_{k},w_{k}\}_{k\in[r]}$ is a collection of length $n$ vectors. The CP-rank is the minimum number $r$ such that $F$ can be written as a sum of $r$ rank-1 tensors, which we refer to as a CP-decomposition. The CP-rank may in fact be larger than the dimension $n$, and furthermore the latent vectors need not be orthogonal as is the case in the matrix setting.

An alternate notion of tensor rank is defined according to the dimension of subspaces corresponding to each mode. Let $F_{(y)}$ denote the unfolded tensor along the $y$-th mode, which is a matrix of dimension $n\times n^{d-1}$. Let columns of $F_{(y)}$ be referred to as mode $y$ fibers of tensor $F$ as depicted in Figure 1.

The Tucker rank, or multilinear rank, is a vector $(r_{1},r_{2},\dots r_{d})$ such that for each mode $\ell\in[d]$, $r_{\ell}$ is the dimension of the column space of $F_{(y)}$. The Tucker rank is also the minimal values of $(r_{1},r_{2},\dots r_{d})$ such that the tensor $F$ can be decomposed according to a multilinear multiplication of a core tensor $\Lambda\in\mathbb{R}^{r_{1}\times r_{2}\times\dots r_{d}}$ with latent factor matrices $Q_{1}\dots Q_{d}$ for $Q_{\ell}\in\mathbb{R}^{n_{\ell}\times r_{\ell}}$, denoted as

[TABLE]

and depicted in Figure 2. The higher order SVD (HOSVD) specifies a unique Tucker decomposition in which the factor matrices $Q_{1}\dots Q_{d}$ are orthonormal and correspond to the left singular vectors of the unfolded tensor along each mode [43].

If the CP-rank is $r$, the Tucker-rank is bounded above by $(r,r,\dots r)$ by constructing a superdiagonal core tensor. If the Tucker rank is $(r_{1},r_{2},\dots r_{d})$, the CP-rank is bounded by the number of nonzero entries in the core tensor, which is at most $r_{1}r_{2}\cdots r_{d}/(\max_{\ell}r_{\ell})$ [43]. While the latent factors of the HOSVD are orthogonal, the latent factors corresponding to the minimal CP-decomposition may not be orthogonal. For simplicity of presentation, we will consider a limited setting where there exists a decomposition of the tensor into the sum of orthogonal rank-1 tensors. This is equivalent to enforcing that the core tensor $\Lambda$ associated to the Tucker decomposition is superdiagonal, or equivalently enforcing that the latent factors in the minimal CP-decomposition are orthogonal. There does not always exist such an orthogonal CP-decomposition, however this class still includes all rank 1 tensors which encompasses the class of instances used to construct the hardness conjecture in [27]. Our results also extend beyond to general tensors as well, though the presentation is simpler in the orthogonal setting.

3 Problem Statement and Model

Consider an $n\times n\times n$ symmetric tensor $F$ generated as follows: For each $u\in[n]$, sample $\theta_{u}\sim U[0,1]$ independently. Let the true underlying tensor $F$ be described by a Lipschitz function $f$ evaluated over the latent variables, $F(u,v,w)=f(\theta_{u},\theta_{v},\theta_{w})$ for $u,v,w\in[n]$. Without loss of generality, we shall assume that $\sup_{u,v,w\in[0,1]}|f(\theta_{u},\theta_{v},\theta_{w})|\leq 1$. For example, if the coordinates of one mode of the tensor represent users or products in an e-commerce platform, one can view the latent variables associated to a coordinate as the unknown “type” of the user/product, which can be thought of as sampled i.i.d. from an underlying population distribution. The latent function $f$ would then describe the expected observed interaction between units of type $\theta_{u},\theta_{v}$, and $\theta_{w}$.

Let $M$ denote the observed symmetric data tensor, and let $\Omega\subseteq[n]^{3}$ denote the set of observed indices. Due to the symmetry, it is sufficient to restrict the index set to triplets $(u,v,w)$ such that $u\leq v\leq w$, as the datapoint is identical for all other permutations of the same triplet. The datapoint at each of these distinct triplets $\{(u,v,w):u\leq v\leq w\}$ is observed independently with probability $p\in(0,1]$, where we assume the observation is corrupted by mean zero independent additive noise terms. For $(u,v,w)\in\Omega$,

[TABLE]

and for $(u,v,w)\notin\Omega$, $M(u,v,w)=\star$444The notation of $\star$ is used to denote the missing observation. When convenient, we shall replace $\star$ by [math] for the purpose of computation. . We shall assume that $|M(u,v,w)|\leq 1$ with probability $1$. We allow $\epsilon_{uvw}$ to have different distributions for different distinct triplets $(u,v,w)$ as long as it is uniformly bounded so as to satisfy the boundedness constraints on $|M(u,v,w)|$ and $|F(u,v,w)|$. When the observations $M(u,v,w)$ are binary, this model is equivalent to the $3$-uniform simple lipschitz hypergraphon in [44], which states a generative model for hypergraphs where the hyperedges consist of size 3 vertex sets. In this setting, $F(u,v,w)\in[0,1]$ would represent the probability of observing the hyperedge $(u,v,w)$, and $M(u,v,w)\in\{0,1\}$ would indicate the presence of the hyperedge $(u,v,w)$. The noise term $\epsilon_{uvw}$ is clearly bounded since the observations are binary. For any application in which the observations $M(u,v,w)$ are bounded, then $F(u,v,w)=\mathbb{E}[M(u,v,w)]$ would also be bounded, such that the boundedness on the noise $\epsilon_{uvw}$ would be reasonable. However, for applications in which $M(u,v,w)$ may not be bounded, or the almost sure bound is very large, we can extend our analysis to allow for sub-Gaussian noise $\epsilon_{uvw}$ rather than uniformly bounded noise, which is further discussed in section 6.

The goal is to recover the underlying tensor $F$ from the incomplete noisy observation $M$ so that the mean squared error (MSE) is small, where MSE for an estimate $\hat{F}$ is defined as

[TABLE]

We will also be interested in the maximum entry-wise error (MEE) defined as

[TABLE]

3.1 Finite spectrum

Consider the setting where the function $f$ has finite spectrum. That is,

[TABLE]

where $r=\Theta(1)$ and $q_{k}(.)$ denotes the orthonormal $\ell_{2}$ eigenfunctions, satisfying $\int_{0}^{1}q_{k}(\theta)^{2}d\theta=1$ and $\int_{0}^{1}q_{k}(\theta)q_{h}(\theta)d\theta=0\text{ for }k\neq h$. Assume that the eigenfunctions are bounded, i.e. $|q_{k}(\theta)|\leq B$ for all $k\in[r]$.

Let $\Lambda$ denote the diagonal $r\times r$ matrix where $\Lambda_{kk}=\lambda_{k}$. Let $Q$ denote the $r\times n$ matrix where $Q_{ka}=q_{k}(\theta_{a})$. Let $\mathcal{Q}$ denote the $r\times\binom{n}{2}$ matrix where $\mathcal{Q}_{kb}=q_{k}(\theta_{b_{1}})q_{k}(\theta_{b_{2}})$ for some $b\in\binom{n}{2}$ that represents the pair of vertices $(b_{1},b_{2})$ for $b_{1}<b_{2}$. The finite spectrum assumption for $f$ implies that the sampled tensor $F$ is such that,

[TABLE]

That is, $F$ has CP-rank at most $r$. In above and in the remainder of the paper, $e_{k}$ denotes a vector with all [math]s but $k$th entry being $1$ of appropriate dimension (here it is $r$).

3.2 Approximately finite spectrum

In general, $f$ may not have finite spectrum, e.g. a generic analytic function $f$. For such a setting, we shall consider $f$ with approximately finite spectrum. Specifically, a function $f:[0,1]^{3}\to\mathbb{R}$, it is said to have $\boldsymbol{\varepsilon}$-approximate finite spectrum with rank $r$ for $\boldsymbol{\varepsilon}\geq 0$ if there exists a symmetric function $f_{r}:[0,1]^{3}\to\mathbb{R}$ such that

[TABLE]

where $r=\Theta(1)$ and $q_{k}(.)$ denotes the orthonormal $\ell_{2}$ eigenfunctions as before. That is, they satisfy $\int_{0}^{1}q_{k}(\theta)^{2}d\theta=1$, $\int_{0}^{1}q_{k}(\theta)q_{h}(\theta)d\theta=0\text{ for }k\neq h$ and $|q_{k}(\theta)|\leq B$ for all $k\in[r]$.

The above describe property of $f$ implies that the sampled tensor $F$ is has $\boldsymbol{\varepsilon}$-approximate rank $r$ such that $F_{r}=\textstyle\sum_{k=1}^{r}\lambda_{k}(Q^{T}e_{k})\otimes(Q^{T}e_{k})\otimes(Q^{T}e_{k})$ and

[TABLE]

3.3 Extensions Beyond Orthogonal CP-rank

The orthogonality conditions on our latent variable decomposition imply that the tensor $F$ can be written as a sum of $r$ rank-1 tensors, where the latent factors are approximately orthogonal. Alternately, this would suggest a Tucker decomposition of the tensor where the core tensor is superdiagonal. Not all tensors admit an orthogonal CP-rank decomposition, but this assumption has also been used in the literature as in [24]. This assumption of the existence of an orthogonal CP-rank decomposition can be relaxed as the main property that our algorithm and analysis use is the orthogonal decomposition of the unfolded tensors along each mode. Our algorithm and analysis will still extend to tensors with Tucker rank bounded by $(r,r,r)$. For a general $(r,r,r)$ Tucker rank tensor, we would instead carry out the analysis with respect to the latent orthogonal factors corresponding to the SVD of the unfolded tensor into a matrix, and the algorithm would use the same procedure to estimate similarities along each of the three modes separately. The presentation is stated for the orthogonal symmetric setting for simplicity.

3.4 Comparision of Assumptions with Literature

In the decomposition of the model $f$ when it has finite spectrum, we assume that the functions $q_{k}$ are orthonormal. This induces a decomposition of tensor $F$ in terms of $Q\in\mathbb{R}^{r\times n}$ with respect to the sampled latent features $\theta\sim U[0,1]$. The rank $r$ of the underlying decomposition is assumed to be $\Theta(1)$. We compare and contrast these with those assumed in the tensor estimation literature.

Most literature on tensor estimation do not impose a distribution on the underlying latent variables, but instead assume deterministic ‘incoherence’ style conditions on the latent singular vectors associated to the underlying tensor decomposition. This plays a similar role to our combined assumption of $q_{k}$ being orthonormal and the latent variables sampled from a uniform distribution so that the mass in the singular vector matrix is roughly uniformly spread. For example, the notion of incoherence used in [27] imposes that the entries of the latent factors are bounded by a constant when the norm of the latent factor vectors scales as $\Theta(\sqrt{n})$. As $\theta_{u}\sim U[0,1]$, it holds that the latent factor vector $(q_{k}(\theta_{1}),q_{k}(\theta_{2}),\dots q_{k}(\theta_{n}))$ will have norm scaling as $\sqrt{n}$. Due to the the boundedness assumption that $|q_{k}(\theta)|\leq B$, our model will satisfy incoherence as defined in [27]. Some of the literature on tensor estimation allows for overcomplete tensors, i.e. $r>n$. While our finite spectrum setup requires $r=\Theta(1)$, the approximately finite spectrum can allow for potentially countably infinite spectrum but with sharply decaying spectrum so that it has $\boldsymbol{\varepsilon}$-approximate rank being $r=\Theta(1)$.

In order to establish our result for the approximately finite spectrum setting, we perform a perturbation analysis wherein each observed entry is perturbed arbitrarily bounded by $\boldsymbol{\varepsilon}$ in magnitude: we shall establish that the resulting estimation error is changed by ${\sf poly}(\boldsymbol{\varepsilon})$, both with respect to the MSE and max-norm. That is, with respect to arbitrary bounded noise in the observations, we are able to characterize the error induced by our method, which is of interest in its own right.

We remark on the Lipschitz property of $f$: the Lipschitz assumption implies that the tensor is “smooth”, and thus there are sets of rows and columns that are similar to one another. As our algorithm is based on a nearest neighbor style approach we need that for any coordinate $u$, there is a significant mass of other coordinates $a$ that are similar to $u$ with respect to the function behavior. Other regularity conditions beyond Lipschitz that would also guarantee sufficiently many “nearest neighbors” would lead to similar results for our algorithm. Lipschitzness also implies approximate low rankness as a Lipschitz function can be approximated by a piecewise constant function, where the number of pieces would then upper bound the rank.

4 Algorithm

The algorithm is a nearest neighbor style in which the first phase is to estimate a distance function between coordinates, denoted $\text{{\sf dist}}(u,a)$ for all $(u,a)\in[n]^{2}$. Given the similarities, for some threshold $\eta$, the algorithm estimates by averaging datapoints from coordinates $(a,b,c)$ for which $\text{{\sf dist}}(u,a)\leq\eta$, $\text{{\sf dist}}(v,b)\leq\eta$, and $\text{{\sf dist}}(w,c)\leq\eta$.

The entry $F(a,b,c)$ depends on a coordinate $a$ through its representation in the eigenspace, given by $Qe_{a}$. Therefore $f(a,b,c)\approx f(u,v,w)$ as long as $Qe_{u}\approx Qe_{a}$, $Qe_{v}\approx Qe_{b}$, and $Qe_{w}\approx Qe_{c}$. Ideally we would like our distance function $\text{{\sf dist}}(u,a)$ to approximate $\|Qe_{u}-Qe_{a}\|_{2}$, but these are hidden latent features that we do not have direct access to.

Let’s start with a thought experiment supposing that the density of observations were $p=\omega(n^{-1})$ and the noise variance is $\sigma^{2}$ for all entries. For a pair of coordinates $u$ and $a$, the expected number of pairs $(b,c)$ such that both $(u,b,c)$ and $(a,b,c)$ are observed is on the order of $p^{2}n^{2}=\omega(1)$. For fixed $\theta_{a},\theta_{u}$, and for randomly sampled $\theta_{b},\theta_{c}$, the expected squared difference between the two corresponding datapoints reflects the distance between $Qe_{a}$ and $Qe_{u}$ along with the overall level of noise,

[TABLE]

where we use the fact that $q_{k}(\cdot)$ are orthonormal. This suggests that approximating $\text{{\sf dist}}(u,a)$ with the average squared difference between datapoints corresponding to pairs $(b,c)$ for which both $(u,b,c)$ and $(a,b,c)$ are observed.

This method does not attain the $p=n^{-3/2}$ sample complexity, as the expected number of pairs $(b,c)$ for which $(a,b,c)$ and $(u,b,c)$ are both observed will go to zero for $p=o(n^{-1})$. This limitation arises due to the fact that when $p=o(n^{-1})$, the observations are extremely sparse. Consider the $n\times\binom{n}{2}$ “flattened” matrix of the tensor where row $u$ correspond to coordinates $u\in[n]$, and columns correspond to pairs of indices, e.g. $(b,c)\in[n]\times[n]$ with $b\leq c$. For any given row $u$, there are very few other rows that share observations along any column with the given row $u$, i.e. the number of ‘neighbors’ of any row index is few. If we wanted to exploit the intuition of the above simple calculations, we have to somehow enrich the neighborhood. We do so by constructing a graph using the non-zero pattern of the matrix as an adjacency matrix. This mirrors the idea from [18] for matrix estimation, which approximates distances by comparing expanded depth $2t+1$ local neighborhoods in the graph representing the sparsity pattern of the unfolded or flattened tensor. In particular, we will construct a statistic $\text{{\sf dist}}(u,a)$ such that with high probability it concentrates around $d(u,a)$ for

[TABLE]

As $F(u,v,w)=f(\theta_{u},\theta_{v},\theta_{w})=\sum_{k=1}^{r}\lambda_{k}q_{k}(\theta_{u})q_{k}(\theta_{v})q_{k}(\theta_{w})$, we can show that if $d(\theta_{u},\theta_{a}),d(\theta_{v},\theta_{b}),$ and $d(\theta_{w},\theta_{c})$ are small, then $F(u,v,w)$ will be close to $F(a,b,c)$. The remaining challenge thus how to approximate $d(u,a)$. Consider a length $2t$ path in the bipartite graph from $u$ to $a$, denoted by $(u,e^{1},x^{1},e^{2},x^{2},\dots x^{t-1},e^{t},a)$, where $x^{1},\dots x^{t-1}$ are distinct coordinates in $[n]\setminus\{u,a\}$, and $e^{1},\dots e^{t}$ consist of pairs in $[n]^{2}$ such that the coordinates represented in these pairs are distinct from each other as well as $\{u,a,x^{1},\dots x^{t-1}\}$. Let us denote the pair $e^{i}=(e^{i}_{1},e^{i}_{2})$. Then the product of weights along this path in expectation conditioned on $\theta_{u},\theta_{a}$ is equal to

[TABLE]

Therefore, the product of weights along the path connecting $u$ to $a$ is a good proxy of quantity $e_{u}^{T}Q^{T}\Lambda^{2t}Qe_{a}$, for paths which do not revisit coordinates. The algorithm first constructs the local neighborhood of depth $2t$ centered at each coordinate $u$, and then connects these neighborhoods to form paths of length $4t+2$ in total, where $t$ is chosen such that for every $u,a$ there are sufficiently many paths used to construct the statistic $\text{{\sf dist}}(u,a)$ to guarantee concentration. The tensor setting requires an important modification of how one constructs the local breadth-first-search (BFS) trees due to the shared latent variables across different modes, as described in step 3 below.

4.1 Formal Description

We provide a formal description of the algorithm below. The crux of the algorithm is to compute similarity between any pair of indices using the matrix obtained by flattening the tensor, and then using a nearest neighbor estimator using these similarities between indices over the tensor structure. Details are as follows.

Step 1: Sample Splitting. Let us assume for simplicity of the analysis that we obtain $2$ independent fresh observation sets of the data, ${\Omega_{1}}$ and ${\Omega_{2}}$. Tensors $M_{1}$ and ${M_{2}}$ contain information from the subset of the data in $M$ associated to ${\Omega_{1}}$ and ${\Omega_{2}}$ respectively. $M_{1}$ is used to compute pairwise similarities between coordinates, and ${M_{2}}$ is used to average over datapoints for the final estimate. Furthermore, we take the coordinates $[n]$ and split it into two sets, $[n]=\{1,2,\dots,n/2\}\cup\{n/2+1,n/2+2,\dots n\}$. Without loss of generality, let’s assume that $n$ is even. Let $\mathcal{V}_{A}$ denote the set of coordinate pairs within set 1 consisting of distinct coordinates, i.e. $\mathcal{V}_{A}=\{(b,c)\in[n/2]^{2}~{}s.t.~{}b<c\}$. Let $\mathcal{V}_{B}$ denote the set of coordinate pairs within set 2 consisting of distinct coordinates, i.e. $\mathcal{V}_{B}=\{(b,c)\in([n]\setminus[n/2])^{2}~{}s.t.~{}b<c\}$. The sizes of $|\mathcal{V}_{A}|$ and $|\mathcal{V}_{B}|$ are both equal to $\binom{n/2}{2}$. We define ${M_{A}}$ to be the $n$-by-$\binom{n/2}{2}$ matrix taking values ${M_{A}}(a,(b,c))=M_{1}(a,b,c)$, where each row corresponds to an original coordinate of the tensor, and each column corresponds to a pair of coordinates $(b,c)\in\mathcal{V}_{A}$ from the original tensor. We define ${M_{B}}$ to be the $n$-by-$\binom{n/2}{2}$ matrix taking values ${M_{B}}(a,(b,c))=M_{1}(a,b,c)$, where each row corresponds to an original coordinate of the tensor, and each column corresponds to a pair of coordinates $(b,c)\in\mathcal{V}_{B}$ from the original tensor. A row-column pair in the matrix corresponds to a triplet of coordinates in the original tensor. We will use matrices ${M_{A}}$ and ${M_{B}}$ to compute similarities or distances between coordinates, and we use tensor ${M_{2}}$ to compute the final estimates via nearest neighbor averaging. The data in ${M_{A}}$ is used to construct depth $2t$ local neighborhoods rooted at each coordinate $u$, and ${M_{B}}$ is used to connect the neighborhoods to form $2t+2$ length paths between any two coordinates $u$ and $a$, which are then used to estimate the similarity between $u$ and $a$. This approach is akin to the technique of “sprinkling” used in random graph analysis, in which we first analyze local neighborhoods formed with the edges in ${M_{A}}$, and then “sprinkle” the edges in ${M_{B}}$ to connect these neighborhoods and argue that there are sufficiently many paths then that connect any two coordinates $u$ and $v$.

Step 2: Construct Bipartite Graph from ${\Omega_{1}},{M_{A}}$. We define a bipartite graph corresponding to the flattened matrix ${M_{A}}$. Construct a graph with vertex set $[n]\cup\mathcal{V}_{A}$. There is an edge between vertex $a\in[n]$ and vertex $(b,c)\in\mathcal{V}_{A}$ if $(a,b,c)\in{\Omega_{1}}$, and the corresponding weight of the edge is $M_{1}(a,b,c)$. Recall that we assumed a symmetric model such that triplets that are permutations of one another will have the same data entry and thus the same edge weight in the associated graph. Figure 3(a) provides a concrete example of a bipartite graph constructed from tensor observations.

Step 3: Expanding the Neighborhood. Consider the graph constructed from ${\Omega_{1}},{M_{A}}$. For each vertex $u\in[n]$, we construct a breadth first search (BFS) tree rooted at vertex $u$ such that the vertices for each depth of the BFS tree consists only of new and previously unvisited coordinates, i.e. if vertex $a\in[n]$ is first visited at depth 4 of the BFS tree, then no vertex corresponding to $(a,b)$ for any $b\in[n]$ can be visited in any subsequent depths greater than 4. Similarly, if $(a^{\prime},b^{\prime})$ is visited in the BFS tree at depth 3, then vertices that include either of these coordinates, i.e. $a^{\prime}$, $b^{\prime}$, $(a^{\prime},c)$, or $(b^{\prime},c)$ for any $c\in[n]$, can not be visited in subsequent depths greater than 3. This restriction is only across different depths; we allow $(a,b)$ and $(a,c)$ to be visited at the same depth of the BFS tree.

There may be multiple valid BFS trees due to different ordering of visiting edges at the same depth. For example, if a vertex at depth $s$ has edges to two different vertices at depth $s-1$ (i.e. two potential parents), only one of the edges can be chosen to maintain the tree property, but either choice is equally valid. Let us assume that when there is more than one option, one of the valid edges are chosen uniformly at random. Figure 3(c) shows valid BFS trees for a bipartite graph constructed from an example tensor.

The graph is bipartite so that each subsequent layer of the BFS tree alternates between the vertex sets $[n]$ and $\mathcal{V}_{A}$. Consider a valid BFS tree rooted at vertex $u\in[n]$ which respects the constraint that no coordinate is visited more than once. We will use $\mathcal{U}_{u,s}\subseteq\mathcal{V}_{A}$ to denote the set of vertices at depth $(2s-1)$ of the BFS tree, and we use $\mathcal{S}_{u,s}\subseteq[n]$ to denote the set of vertices at depth $2s$ of the BFS tree. Let $\mathcal{B}_{u,s}\subset[n]\cup\mathcal{V}_{A}$ denote the set of vertices which are visited in the first $s$ layers of the BFS tree,

[TABLE]

We will overload notation and sometimes use $\mathcal{B}_{u,s}$ to denote the subset of coordinates in $[n]$ visited in the first $s$ layers of the BFS tree, including both visited single coordinate vertices or coordinates in vertices $\mathcal{V}_{A}$, i.e.

[TABLE]

Let $\mathcal{G}(\mathcal{B}_{u,s})$ denote all the information corresponding to the subgraph restricted to the first $s$ layers of the BFS tree rooted at $u$. This includes the vertex set $\mathcal{B}_{u,s}$, the latent variables $\{\theta_{a}\}_{a\in\mathcal{B}_{u,s}}$ and the edge weights $\{M_{1}(a,b,c)\}_{a,(b,c)\in\mathcal{B}_{u,s}}$.

We define neighborhood vectors which represent the different layers of the BFS tree. Let $N_{u,s}\in[0,1]^{n}$ be associated to set $\mathcal{S}_{u,s}$, where the $a$-th coordinate is equal to the product of weights along the path from $u$ to $a$ in the BFS tree for $a\in\mathcal{S}_{u,s}$. Similarly, let $W_{u,s}\in[0,1]^{\mathcal{V}_{A}}$ be associated to set $\mathcal{U}_{u,s}$, where the $(b,c)$-th coordinate is equal to the product of weights along the path from $u$ to $(b,c)$ in the BFS tree for $(b,c)\in\mathcal{U}_{u,s}$. For $a\in[n]$, let $\pi_{u}(a)$ denote the parent of $a$ in the valid BFS tree rooted at vertex $u$. For $(b,c)\in\mathcal{V}_{A}$, let $\pi_{u}(b,c)$ denote the parent of $(b,c)$ in the BFS tree rooted at vertex $u$. We can define the neighborhood vectors recursively,

[TABLE]

and $N_{u,0}=e_{u}$. Let $\tilde{N}_{u,s}$ denote the normalized vector $\tilde{N}_{u,s}=N_{u,s}/|\mathcal{S}_{u,s}|$ and let $\tilde{W}_{u,s}$ denote the normalized vector $\tilde{W}_{u,s}=W_{u,s}/|\mathcal{U}_{u,s}|$. Figure 3(d) illustrates the neighborhood sets and vectors for a valid BFS tree. Recall from Eq. (4.2) that conditioned on $\theta_{u}$ and $\theta_{a}$, $\mathbb{E}[N_{u,s}(a)~{}|~{}\theta_{u},\theta_{a}]=\mathbb{P}\left(a\in\mathcal{S}_{us}\right)e_{u}^{T}Q^{T}\Lambda^{2s}Qe_{a}$. Furthermore the event that $a\in\mathcal{S}_{us}$ only depends on the presence of the edges as determined by the Bernoulli uniform sampling, and is thus independent from the latent variable $\theta_{a}$. We will show in a subsequent Lemma 7.2 that $e_{k}^{T}Q\tilde{N}_{u,t}\approx e_{k}^{T}\Lambda^{2t}Qe_{u}$, implying that the neighborhood vector $\tilde{N}_{u,t}$, which is constructed from products of weights over length $2s$ paths originating at $u$, is a statistic that is approximates $\Lambda^{2t}Qe_{u}$. Similar calculations show that $\mathbb{E}[W_{u,s}(b,c)~{}|~{}\theta_{u},\theta_{b},\theta_{c}]=\mathbb{P}\left((b,c)\in\mathcal{U}_{us}\right)\sum_{k=1}^{r}\lambda_{k}^{2t-1}q_{k}(\theta_{u})q_{k}(\theta_{b})q_{k}(\theta_{c})$.

Step 4: Computing the distances using ${M_{B}}$. Let

[TABLE]

A heuristic for the distance would be

[TABLE]

For technical reasons that facilitate cleaner analysis, we use the following distance calculations. There are two deviations from the equation in (4.4). First we exclude $a=b$ from the summation. Second we exclude coordinates for $\alpha$ or $\beta$ that have been visited previously in $\mathcal{B}_{u,2t}$ or $\mathcal{B}_{v,2t}$. Define distance as

[TABLE]

We will show in Lemma 6.2 that $\text{{\sf dist}}(u,v)\approx d(u,v):=\|\Lambda^{2t+1}Q(e_{u}-e_{v})\|_{2}^{2}$. The estimate is constructed by averaging over the product of weights over paths from $u$ to $v$, where the term $N_{u,t}(a)N_{v,t}(b){M_{B}}(a,(\alpha,\beta)){M_{B}}(b,(\alpha,\beta))$ in Eq. (4.6) is the product of weights over the path that goes from $u$ to $a$ to $(\alpha,\beta)$ to $b$ to $v$, as $N_{u,t}(a)$ represents the products of weights on the path from $u$ to $a$ and $N_{v,t}(b)$ represents the products of weights on the path from $b$ to $v$. The parameter $t$ is chosen such that there are sufficiently many paths that we are averaging over in order to reduce the noise. In particular, the choice of $t\geq\ln(n)/2\ln(p^{2}n^{3})$ implies that $|\mathcal{S}_{u,t}|\geq(p^{2}n^{3})^{t}=\Omega(n^{1/2})$, such that the number of paths the estimator averages over is approximately $n^{2}p^{2}|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|=\Omega(p^{2}n^{3})$. This rough calculation highlights that $p$ must be $\omega(n^{-3/2})$ to guarantee that the number of paths used to compute $\text{{\sf dist}}(u,v)$ is increasing with $n$.

Step 5: Averaging datapoints to produce final estimate. Let ${\Omega_{2}}_{uvw}$ denote the set of indices $(a,b,c)$ such that $a\leq b\leq c$, $(a,b,c)\in{\Omega_{2}}$, and the estimated distances $\text{{\sf dist}}(u,a)$, $\text{{\sf dist}}(v,b)$, $\text{{\sf dist}}(w,c)$ are all less than some chosen threshold parameter $\eta$. The final estimate averages the datapoints corresponding to indices in ${\Omega_{2}}_{uvw}$,

[TABLE]

4.2 Difference between tensor and matrix setting

The modifications in the construction of the breadth-first-search (BFS) tree for the tensor setting relative to the matrix setting are critical to the analysis. If we simply considered the classical construction of a BFS tree in the associated bipartite graph (as the matrix setting uses), this would lead to higher variance and bias due to the correlations of vertices sharing common latent variables associated to the same underlying coordinates of the tensor. Alternatively, if one constructed a BFS tree by not allowing any coordinate of the tensor to be visited more than once, this would also lead to suboptimal results as it would throw away too many entries, limiting the computed statistic to only order $n$ data points. Our final algorithm, which allows for vertices with shared coordinates in the same depth of the BFS but not across different depths, is carefully chosen in order to break dependencies across different depths of the BFS tree, while still allowing for sufficient expansion in each depth.

To extend the algorithm to $d$-order tensors for $d>3$, we will compute similarities between $u,v$ via a similar computation as described in Steps 2-4 above, except it would be applied to the $n\times n^{d-1}$ matrix and associated bipartite graph corresponding to an unfolding of the tensor. Given the similarity estimates for pairwise coordinates, the final estimate would result from a standard nearest neighbor estimator over the high dimensional tensor. The primary part of the proof that would need to be modified is the analysis of the neighborhood vectors $N_{u,s}$ and $W_{us}$ in step 3, which may involve constraining the growth of the BFS trees such that they extend deep enough before exhausting the visited coordinates.

5 Main Result

We provide an upper bound on the mean squared error (MSE) as well as the max entry-wise error (MEE) for the algorithm, showing that both the MSE and the MEE converge to zero as long as $p=n^{-3/2+\kappa}$ for some $\kappa>0$. Our result implies that the simple variant of collaborative filtering algorithm based on estimating similarities produces a consistent estimator when the tensor latent function has finite spectrum or low rank. Further we show that it is robust to arbitrary, additive perturbation in that the estimation error increases gracefully in the amount of perturbation. To the best of our knowledge, such robustness to arbitrary bounded additive noise with respect to max-norm estimation is first of its kind in the literature on tensor estimation.

5.1 Finite spectrum

We establish consistency of our estimator with respect to MSE and max-norm error of the algorithm when the underlying $f$ has finite spectrum, i.e. rank $r$ model with $r=\Theta(1)$.

Theorem 5.1**.**

We assume that the function $f$ is rank $r$, $L$-Lipschitz and that $\theta\sim U[0,1]$. Assume that $p=n^{-3/2+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$. Let $t$ be defined according to (4.3). For any arbitrarily small $\psi\in(0,\min(\kappa,\frac{3}{8}))$, choose the threshold

[TABLE]

The algorithm produces estimates so that,

[TABLE]

and

[TABLE]

with probability $1-O\left(n^{4}\exp(-\Theta(n^{2\psi}))\right)$.

The error bounds in Theorem 5.1 imply that our estimator is consistent as long as $p=n^{-3/2+\kappa}$ for some $\kappa>0$, with a MSE that scales as $O(1/pn^{3/2})$. The threshold of $p=\Omega(n^{-3/2})$ is optimal, and furthermore this requirement is precisely the threshold at which the constructed bipartite graph in Step 2 of the algorithm is fully connected. Below the connectivity threshold, the graph will be disconnected into small components with insufficient information to recover the expected value of edges across disconnected components.

In comparison to the literature, [31] prove that the minimax optimal MSE is $O(1/pn^{2})$, which is achieved via spectral initialization followed by power iteration [31] or gradient descent [40] as long as $p=\Omega(n^{-3/2})$. While our result achieves a similar sample complexity threshold, our MSE rate is suboptimal by a factor of $\sqrt{n}$. A limitation of neighborhood smoothing is that we do not achieve exact recovery under the noiseless setting, and we do not achieve the minimax optimal rates. It is unclear whether the gap is due to a limitation in the analysis or the algorithm. A benefit of our analysis in contrast to the literature is that the neighborhood smoothing approach can more easily deal with approximate low rank settings as arise under smoothness, as presented in Theorem 5.2. While low rank is often a useful modeling concept for real world datasets, in reality most real-world applications are likely only approximately low rank rather than exactly low rank.

5.2 Approximately finite spectrum

For approximate rank $r$ model, we establish a natural perturbation result for the algorithm. Specifically, if the underlying model has $\boldsymbol{\varepsilon}$-approximate rank $r$, then we argue that the result of Theorem 5.1 remain true, both with respect to MSE and max-norm error, with perturbation amount of ${\sf poly}(\boldsymbol{\varepsilon})$.

Theorem 5.2**.**

We assume that the function $f$ has $\boldsymbol{\varepsilon}$-approximate rank $r$, $L$-Lipschitz and that $\theta\sim U[0,1]$. Assume that $p=n^{-3/2+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$. Choosing $t$ according to (4.3), it follows that $t=\lceil\frac{1}{4\kappa}\rceil$. For any arbitrarily small $\psi\in(0,\min(\kappa,\frac{3}{8}))$, choose the threshold

[TABLE]

The algorithm produces estimates so that,

[TABLE]

and

[TABLE]

with probability $1-O\left(n^{4}\exp(-\Theta(n^{2\psi}))\right)-O(n^{-2})$.

As the entries of $F$ are normalized such that $\|F\|_{\max}\leq 1$, the bound is meaningful when $\boldsymbol{\varepsilon}<1$, in which case the dominating term of the additional error due to the perturbation is linear in $\boldsymbol{\varepsilon}$, as $t$ is a constant. The proof of Theorem 5.2 relies on the following observation: the distribution of the data under the setting where the latent function $f$ has $\boldsymbol{\varepsilon}$-rank $r$ is equivalent to the distribution of data generated according to the rank $r$ approximation of $f$ and then adding a deterministic perturbation to each observation accounting for the difference between $f$ and its rank $r$ approximation $f_{r}$, which is entrywise bounded by $\boldsymbol{\varepsilon}$. In particular, the proof of Theorem 5.2 shows that under arbitary deterministic perturbation of a rank $r$ model where the perturbation is bounded by $\boldsymbol{\varepsilon}$, the estimation error is perturbed by at most ${\sf poly}(\boldsymbol{\varepsilon})$. As a byproduct, our result proves that our estimator is robust to arbitrary deterministic bounded noise in the observations.

The approximation guarantee depends on the spectral decay. Since the analysis allows any arbitrary adversarial model for the $\boldsymbol{\varepsilon}$ deviation from a low rank model, the resulting guarantee depending polynomially in $\boldsymbol{\varepsilon}$ is qualitatively best one can hope for. By definition, the minimal error must be lower bounded by $\boldsymbol{\varepsilon}$, but determining the best achievable error as a function of $\boldsymbol{\varepsilon}$ is an important open question for future investigation. It is worth noting that no other prior work addresses such a robust error model.

5.3 Reducing Computational Complexity

The computational complexity can be estimated by analyzing steps 3-5 of the algorithm. Step 3 costs $O(pn^{4})$, as there are $n$ BFS trees to construct, which each take at most $pn^{3}$ edge traversals as there are at most order $pn^{3}$ edges in the constructed graph. Step 4 costs $O(p^{2}n^{6})$ as there are order $n^{2}$ pairwise distances to compute, and each computed distance involves sums over terms indexed by $a,b,\alpha,\beta\in[n]^{4}$ where $(a,\alpha,\beta)$ and $(b,\alpha,\beta)$ are in the observation set. As the sparsity of the dataset is $p$, this results in order $p^{2}n^{4}$ nonzero terms in the summation, each of which is the product of 4 quantities, taking $O(1)$ to compute. Step 5 costs $O(pn^{6})$ as there are $\Theta(n^{3})$ triplets we need to estimate, and each involves averaging at most $O(pn^{3})$ datapoints. In summary, the computation cost of the entire method, for $p=n^{-\frac{3}{2}+\kappa}$ is $O(pn^{4}+p^{2}n^{6}+pn^{6})$ where the cost in Step 5 dominates.

This computation cost can be improved drastically. For example, as explained in [18], by use of ‘representative’ or ’anchor’ vertices chosen as random, the algorithm can instead cluster the vertices with respect to these anchor vertices and learn a block constant estimate, significantly reducing the involved computation. If there are $y$ anchor vertices, then Step 4 reduces to only computing pairwise distances between $\binom{y}{2}+ny$ pairs of vertices, as non-anchor vertices are only compared to the small set of $y$ anchor vertices. Step 5 reduces to only estimating $\binom{y}{3}$ entries of the tensor corresponding to combinations of the anchor vertices, and then extrapolating the estimate to other vertices assigned to the same cluster. This would result in a computational cost of $O(pn^{4}+(y^{2}+ny)p^{2}n^{4}+y^{3}pn^{3})$. When $p=n^{-\frac{3}{2}+\kappa}$, our proof indicates that by choosing $y=\Theta((p^{2}n^{3})^{1/4})=\Theta(n^{\kappa/2})$, the corresponding block constant estimator would achieve the same rates on the MSE and MEE as presented in Theorem 5.1, while requiring a reduced computational complexity of $O(n^{5/2+\kappa}+n^{2+5\kappa/2})$.

Corollary 5.3**.**

We assume that the function $f$ is rank $r$, $L$-Lipschitz and that $\theta\sim U[0,1]$. Assume that $p=n^{-3/2+\kappa}$ for some $\kappa\in(0,\frac{1}{2})$. Let $t$ be defined as per (4.3). For any arbitrarily small $\psi\in(0,\min(\kappa,\frac{3}{8}))$, choose the threshold

[TABLE]

The modified algorithm which subsamples $y=\Omega((p^{2}n^{3})^{1/4})=\Omega(n^{\kappa/2})$ anchor vertices at random and uses them to cluster the vertices to learn a block constant estimate will achieve

[TABLE]

and

[TABLE]

with probability $1-O\left(n^{4}\exp(-\Theta(n^{2\psi}))\right)$.

5.4 Discussion of Assumptions

We assumed in our algorithm and analysis that we had two fresh samples of the dataset, $M_{1}$ and $M_{2}$. The dataset $M_{1}$ is used to estimate distances between coordinates, and the dataset $M_{2}$ is used to compute the final nearest neighbor estimates. Given only a single dataset, the same theoretical results can also be shown by simply splitting the samples uniformly into two sets, one used to estimate distances and one used to compute the nearest neighbor estimates. as we are considering the sparse regime with $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$, the two subsets after sample splitting will be nearly independent, such that the analysis only needs to be slightly modified. This is formally handled in the paper on collaborative filtering for matrix estimation in [18].

Our model and analysis assumes that the latent variables $\{\theta_{u}\}_{u\in[n]}$ are sampled uniformly on the unit interval, and that the function $f$ is Lipschitz with respect to $\theta$. This assumption can in fact be relaxed significantly, as it is only used in the final step of the proof in analyzing the nearest neighbor estimator. Proving that the distance estimates concentrate well does not require these assumptions, in particular it primarily uses the low rank assumption. Given that the distance estimate concentrates well, the analysis of the nearest neighbor estimator depends on the local measure, i.e. what fraction of other coordinates have similar function values so that the estimated distance is small. We used Lipschitzness and uniform distribution on the unit interval in order to lower bound the fraction of nearby coordinates, however many other properties would also lead to such a bound. The dependence of the noisy nearest neighbor estimator on the local measure is discussed in detail in [16]. Similar extensions as presented in [16] would apply for our analysis here, leading to consistency and convergence rate bounds for examples including when

• •

the latent space has only finitely many elements, or equivalently the distribution of $\theta$ has finite support;

• •

the latent space is the unit hypercube in a finite dimensional space and the latent function is Lipschitz;

• •

the latent space is a complete, separable metric space, i.e. Polish space, with bounded diameter and the latent function is Lipschitz.

Although our stated results assume a symmetric tensor, the results naturally extend to asymmetric $(n_{1}\times n_{2}\times n_{3})$ tensors as long as $n_{1},n_{2},$ and $n_{3}$ are proportional to one another. Our analysis can be modified for the asymmetric setting, or one can reduce the asymmetric tensor to a $(n\times n\times n)$ symmetric tensor where $n=n_{1}+n_{2}+n_{3}$, and the coordinates of the new tensor consists of the union of the coordinates in all three dimensions of the asymmetric tensor. The results applied to this larger tensor would still hold with adjustments of the model allowing for piecewise Lipschitz functions.

In the proof sketch that follows below, we show that for the 3-order tensor, the sample complexity threshold of $p=\omega(n^{-3/2})$ directly equals the density of observations needed to guarantee the bipartite graph is connected with high probability. Although our stated results assume a 3-order tensor, our algorithm and analysis can be likely extended to general $d$-order tensors. The proof would instead require analysis of the $n\times n^{d-1}$ matrix and associated bipartite graph corresponding to the unfolding of the tensor. The bipartite graph would consist of vertex sets $[n]$ and $[n]^{d-1}$.

Remark 5.1*.*

In order for the vertices in $[n]$ to be fully connected to each other, $p$ needs to be $\Omega(n^{d/2})$, which can be proved using the standard branching process analysis as is used to prove the connectivity threshold of an Erdos Renyi graph [45]. Let $X_{u}$ denote the set of vertices $v\in[n]$ such that the distance between $u$ and $v$ in the bipartite graph is 2. It follows that $\mathbb{P}\left(v\in X_{u}\right)=1-(1-p^{2})^{n^{d-1}}$, such that $\mathbb{E}[|X_{u}|]=(n-1)(1-(1-p^{2})^{n^{d-1}})=\Theta(p^{2}n^{d})$. As a result, if $p=o(n^{d/2})$, it follows that $\mathbb{P}\left(X_{u}=\emptyset\right)=1-\mathbb{P}\left(|X_{v}|\geq 1\right)\geq 1-\mathbb{E}[|X_{u}|]=1-o(1)$, such that with probability tending to 1 the vertex $u$ will be isolated, i.e. not connected to any other vertex $v\in[n]$. To prove that the graph is connected for $p=\tilde{\Omega}(n^{d/2})$, one would relate the growth of a local neighborhood in the graph to an appropriate branching process, where the expected number of descendents alternates between $pn$ and $pn^{d-1}$. To formalize the argument, one would then argue that the branching process survives to infinity if $p^{2}n^{d}$ is sufficiently large, e.g $\text{polylog}(n)$, and also argue that the branching process is a reasonable approximation for the neighborhood growth of the graph until a linear number of vertices are visited. The requirement for graph connectivity in our algorithm and analysis arises from the similarity computation $\text{{\sf dist}}(u,v)$, which involves products of weights over paths in the graph that connect $u$ and $v$. The fact that $n^{d/2}$ also corresponds to the connectivity threshold in the corresponding bipartite graph sheds light on the computational lower bound for tensor completion in [27], giving another way to explain the $\Omega(n^{d/2})$ lower bound.

6 Proof

In this section, we present the proof for Theorem 5.1. The proof outline is similar to the matrix setting in [18], in that the core of the analysis is proving that the distance function as defined in (4.5) concentrates appropriately and captures an appropriate notion of distance that enables the classical “nearest neighbor” algorithm to be effective. However, due to high-dependencies across latent factors associated with columns that share tensor coordinates, the concentration of the BFS neighborhood expansion in section 7.2 requires a new argument beyond the simple martingale argument in the matrix setting. This involves a careful application of the concentration of U-statistics. Furthermore, the concentration of the distance calculation in Eq (4.6) as analyzed in section 7.4 requires a new argument relating the computed statistic to a thresholded variant more amenable to analysis. This is due to both the dependencies in the latent factors along with the lopsidedness in the dimensions so that straightforward applications of standard concentration results are too weak and insufficient to drive the error to zero.

While the proof is stated for bounded observations, i.e. bounded noise, the result can be extended to sub-Gaussian noise rather than uniformly bounded noise. This would involve showing that the norms of the neighborhood vectors $N_{us}$ and $W_{us}$ are well-controlled such that the application of Hoeffding’s inequality used in Lemmas 7.4 and 7.5 still hold. Additionally the proof of 7.7 naively bounds the products of weights over a path in absolute value by 1; if the noise were not bounded but sub-Gaussian, one would have to additionally argue that the product of the weights would be sufficienty controlled with high probability.

The critical lemma that the proof hinges on shows that the computed similarities, i.e. $\text{{\sf dist}}(u,a)$, concentrates around the function $d(u,a)=\|\Lambda^{2t+1}Q(e_{u}-e_{a})\|_{2}^{2}$. This then implies that if $d(u,a)$, $d(v,b)$ and $d(w,c)$ are small, the function value $F(u,v,w)$ would be close to $F(a,b,c)$. Additionally, we use Lipschitzness of the latent function $f$ along with the assumption that $\theta_{u}$ are sampled independently from $U[0,1]$ to argue that for any $u$, there is a sufficiently large set of other coordinates $a$ such that $d(u,a)$ is small. If these properties hold, then a simple analysis of nearest neighbor averaging using $\text{{\sf dist}}(u,a)$ to determine the neighbors will result in a bias variance tradeoff that can be tuned to show our final results. As such, the complexity of the proof revolves around showing that the computed $\text{{\sf dist}}(u,a)$ concentrates around $d(u,a)$. This involves a delicate analysis which involves first arguing that the normalized neighborhood vectors $\tilde{N}_{u,t}$ satisfy $e_{k}^{T}Q\tilde{N}_{u,t}\approx e_{k}^{T}\Lambda^{2t}Qe_{u}$, which involves martingale concentration as well as concentration of appropriately defined U-statistics. Subsequently we need to argue that the statistic

[TABLE]

A challenge in the analysis is that each term in the sum has very small probability of being nonzero such that the sum is sparse enough that the standard concentration inequalities are not tight enough. Thus we relate the sum to a thresholded variant and use a tighter approximation of the binomial cdf to obtain the desired bound.

6.1 Analyzing Noisy Nearest Neighbors

We start by stating an important Lemma 6.1, adapted from [18] that characterizes the error of the noisy nearest neighbor algorithm. Recall that our algorithm estimates $F(u,v,w)$, i.e. $f(\theta_{u},\theta_{v},\theta_{w})$, according to (4.7), which simply averages over data-points ${M_{2}}(a,b,c)$ corresponding to tuples $(a,b,c)$ for which $a$ is close to $u$, $b$ is close to $v$ and $c$ is close to $w$ according to the estimated distance function. The choice of parameter $\eta$ allows for tradeoff between bias and variance of the algorithm.

We first argue that the data-driven distance estimates dist will concentrate around an ideal data-independent distance $d(\theta_{u},\theta_{v})$ for $d:[0,1]^{2}\to\mathbb{R}_{+}$. We subsequently argue that the nearest neighbor estimate produced by (4.7) using $d(\theta_{u},\theta_{v})$ in place of $\text{{\sf dist}}(u,v)$ will yield a good estimate by properly choosing the threshold $\eta$ to tradeoff between bias and variance. The bias will depend on the local geometry of the function $f$ relative to the distances defined by $d$. The variance depends on the measure of the latent variables $\{\theta_{u}\}_{u\in[n]}$ relative to the distances defined by $d$, i.e. the number of observed tuples $(a,b,c)\in{\Omega_{2}}$ such that $d(\theta_{u},\theta_{a})\leq\eta$, $d(\theta_{v},\theta_{b})\leq\eta$ and $d(\theta_{w},\theta_{c})\leq\eta$ needs to be sufficiently large. We formalize the above stated desired properties.

Property 6.1* (Good Distance).*

We call an ideal distance function $d:[0,1]^{2}\to\mathbb{R}_{+}$ to be a bias-good distance function for some $\text{{\sf bias}}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ if for any given $\eta>0$ it follows that $|f(\theta_{a},\theta_{b},\theta_{c})-f(\theta_{u},\theta_{v},\theta_{w})|\leq\text{{\sf bias}}(\eta)$ for all $(\theta_{a},\theta_{b},\theta_{c},\theta_{u},\theta_{v},\theta_{w})\in[0,1]^{4}$ such that $d(\theta_{u},\theta_{a})\leq\eta$, $d(\theta_{v},\theta_{b})\leq\eta$ and $d(\theta_{w},\theta_{c})\leq\eta$.

Property 6.2* (Good Distance Estimation).*

For some $\Delta>0$, we call distance $\hat{d}:[n]^{2}\to\mathbb{R}_{+}$ a $\Delta$-good estimate for ideal distance $d:[0,1]^{2}\to\mathbb{R}_{+}$, if $|d(\theta_{u},\theta_{a})-\hat{d}(u,a)|\leq\Delta$ for all $(u,a)\in[n]^{2}$.

Property 6.3* (Sufficient Representation).*

The collection of coordinate latent variables $\{\theta_{u}\}_{u\in[n]}$ is called meas-represented for some $\text{{\sf meas}}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ if for any $u\in[n]$ and $\eta^{\prime}>0$, $\frac{1}{n}\sum_{a\in[n]}\mathbb{I}_{(d(u,a)\leq\eta^{\prime})}\geq\text{{\sf meas}}(\eta^{\prime})$.

Lemma 6.1**.**

Assume that property 6.1 holds with probability $1$, property 6.2 holds for any given pair $u,a\in[n]$ with probability $1-\alpha_{1}$, and property 6.3 holds with probability $1-\alpha_{2}$ for some $\eta,\Delta,$ and $\eta^{\prime}=\eta-\Delta$; in particular $d$ is a bias-good distance function, $\hat{d}=\text{{\sf dist}}$ as estimated from ${M_{A}}$ is a $\Delta$-good distance estimate for $d$, and $\{\theta_{u}\}_{u\in[n]}$ is meas-represented. Then noisy nearest neighbor estimate $\hat{F}$ computed according to (4.7) satisfies

[TABLE]

for any $\delta\in(0,1)$. Furthermore, for any $\delta^{\prime}\in(0,1)$ and $(u,v,w)\in[n]^{3}$,

[TABLE]

with probability at least

[TABLE]

The proof of Lemma 6.1 is a modification from [18] and is included in the Appendix.

6.2 Proofs of Theorems 5.1 and 5.2

Proof.

We prove that as long as $p=n^{-3/2+\kappa}$ for any $\kappa\in(0,\frac{1}{2})$, with high probability, properties 6.1-6.3 hold for an appropriately chosen function $d$, and for distance estimates $\hat{d}=\text{{\sf dist}}$ computed according to (4.5) with $t$ defined in (4.3). We subsequently use Lemma 6.1 to conclude Theorem 5.1 and Theorem 5.2. The proofs of Properties 6.1 and 6.3 are identical in Theorem 5.1 and Theorem 5.2, while that of property 6.2 differ. For Theorem 5.1, we utilize Lemma 6.2 while for Theorem 5.2, we utilize Lemma 6.3. The proof of Theorem 5.2 follows nearly the same argument, where $f$ will be replaced by the rank $r$ approximation $f_{r}$, c.f. (3.2).

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

[TABLE]

Recall that $t$ is defined in (4.3). Since $p=n^{-3/2+\kappa}$ and $\kappa\in(0,\frac{1}{2})$, we have that

[TABLE]

We want to show that there exists $\text{{\sf bias}}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ so that $|(f(\theta_{a},\theta_{b},\theta_{c})-f(\theta_{u},\theta_{v},\theta_{w}))|\leq\text{{\sf bias}}(\eta)$ for any $\eta>0$ and $(u,a,v,b,w,c)\in[n]^{3}$ such that $d(\theta_{u},\theta_{a})\leq\eta$, $d(\theta_{v},\theta_{b})\leq\eta$ and $d(\theta_{w},\theta_{c})\leq\eta$. Consider

[TABLE]

Now

[TABLE]

In above, (a) follows from the $\|q_{k}(\cdot)\|_{\infty}\leq B$ for all $k$. Repeating this argument to bound the other terms in (6.3), we obtain that

[TABLE]

In summary, property 6.1 is satisfied for distance function $d$ defined according to (6.1) and $\text{{\sf bias}}(\eta)=3B^{2}|\lambda_{r}|^{-2t}\sqrt{r\eta}$.

Good distance estimate $\hat{d}$ and Property 6.2. We state the following Lemma whose proof is delegated to Section 7.

Lemma 6.2**.**

Given $f$ with rank $r$, assume that $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$. Let $\hat{d}=\text{{\sf dist}}$ as defined in (4.5). Then for any $(u,a)\in[n]^{2}$, for any $\psi\in(0,\kappa)$,

[TABLE]

with probability at least 1-O\Big{(}\exp(-n^{2\psi}(1-o(1)))\Big{)}.

Lemma 6.2 implies that property 6.2 holds with probability $1-o(1)$ for $\Delta=\Theta\left(r\lambda_{\max}^{4t}n^{-(\kappa-\psi)}\right)$ when $f$ has rank $r$.

Lemma 6.3**.**

Given $f$ with $\boldsymbol{\varepsilon}$-approximate rank $r$ for $\boldsymbol{\varepsilon}\geq 0$, assume that $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$. Let $\hat{d}=\text{{\sf dist}}$ as defined in (4.5). Then for any $(u,a)\in[n]^{2}$, for any $\psi\in(0,\kappa)$,

[TABLE]

with probability at least 1-O\Big{(}\exp(-n^{2\psi}(1-o(1)))\Big{)}-O\Big{(}n^{-6}\Big{)}.

Lemma 6.3 implies that property 6.2 holds with probability $1-o(1)$ for

[TABLE]

when $f$ has $\boldsymbol{\varepsilon}$-approximate rank $r$.

Sufficient representation and Property 6.3. Since $f$ is $L$-Lipschitz, the distance $d$ as defined in (6.1) is bounded above by the squared $\ell_{2}$ distance:

[TABLE]

where in (a) we have used the fact that $q_{k}(\cdot),k\in[r]$ are orthonormal with respect to uniform distribution over $[0,1]$. We assumed that the latent parameters $\{\theta_{u}\}_{u\in[n]}$ are sampled i.i.d. uniformly over $[0,1]$. Therefore, for any $\theta_{u}\in[0,1]$, for any $v\in[n]$ and $\eta^{\prime}>0$,

[TABLE]

Let us define

[TABLE]

for all $\eta^{\prime}\in(0,|\lambda_{1}|^{4t}L^{2})$. By an application of Chernoff’s bound and a simple majorization argument, it follows that for all $\eta^{\prime}\in(0,|\lambda_{1}|^{4t}L^{2})$ and $\delta\in(0,1)$,

[TABLE]

By using union bound over all $n$ indices, it follows that for any $\eta^{\prime}\in(0,|\lambda_{1}|^{4t}L^{2})$, with probability at least $1-n\exp\left(-\frac{\delta^{2}(n-1)\sqrt{\eta^{\prime}}}{2|\lambda_{1}|^{2t}L}\right)$, property 6.3 is satisfied with meas as defined in (6.8).

Concluding Proof of Theorem 5.1. In summary, property 6.1 holds with probability $1$, by Lemma 6.2 property 6.2 holds for a given tuple $(u,a)\in[n]^{2}$ with probability $1-\alpha_{1}$ where \alpha_{1}=O\Big{(}\exp(-n^{2\psi}(1-o(1)))\Big{)} for $\psi\in(0,\min(\kappa,\frac{3}{8}))$ and $\kappa\in(0,\frac{1}{2})$, property 6.3 holds with probability $1-\alpha_{2}$ where $\alpha_{2}=n\exp\left(-\frac{\delta^{2}(n-1)\sqrt{\eta-\Delta}}{2|\lambda_{1}|^{2t}L}\right)$ with distance estimate $\hat{d}=\text{{\sf dist}}$ defined in (4.5) with

[TABLE]

for any $\eta>0$, $\delta\in(0,1)$ and $\eta^{\prime}=\eta-\Delta\in(0,|\lambda_{1}|^{4t}L^{2})$. By substituting the expressions for bias, meas, and $\alpha$ into Lemma 6.1, it follows that

[TABLE]

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

[TABLE]

with probability at least

[TABLE]

By selecting \eta=\Theta\big{(}\Delta\big{)}=\Theta(r\lambda_{\max}^{4t}n^{-(\kappa-\psi)}) with a large enough constant so that $\eta-\Delta=\Theta(\eta)$, it follows that by the conditions that $\psi>0$ and $\kappa<\frac{1}{2}$,

[TABLE]

By substituting this choice of $\eta$ and $\delta=\frac{1}{2}$, it follows that

[TABLE]

By choosing $\delta^{\prime}=n^{-(\kappa-\psi)/2}$ such that $\delta^{\prime}=\Theta(\sqrt{\eta})$ and $\delta^{\prime 2}p(\sqrt{\eta-\Delta}n)^{3}=\Omega(n^{\frac{3}{4}})=\Omega(n^{2\psi})$ because $\psi<\frac{3}{8}$. Therefore, by substituting into (6.10), it follows that for any given $(u,v,w)\in[n]^{3}$, with probability 1-O(n\exp(-\Theta(n^{2\psi}))\Big{)},

[TABLE]

Using union bound over choices of $(u,v,w)\in[n]^{3}$, it follows that the maximum entry-wise error is bounded above by O\Big{(}n^{-(\kappa-\psi)/2}\Big{)} with probability 1-O(n^{4}\exp(-\Theta(n^{2\psi}))\Big{)}. This completes the proof of Theorem 5.1. ∎

Concluding Proof of Theorem 5.2. We follow similar line of argument as for proof of Theorem 5.1. As noted earlier, property 6.1 holds with probability $1$, by Lemma 6.3 property 6.2 holds for a given tuple $(u,a)\in[n]^{2}$ with probability $1-\alpha_{1}$ where \alpha_{1}=O\Big{(}\exp(-n^{2\psi}(1-o(1)))+n^{-6}\Big{)} for $\psi\in(0,\min(\kappa,\frac{3}{8}))$ and $\kappa\in(0,\frac{1}{2})$, property 6.3 holds with probability $1-\alpha_{2}$ where $\alpha_{2}=n\exp\left(-\frac{\delta^{2}(n-1)\sqrt{\eta-\Delta}}{2|\lambda_{1}|^{2t}L}\right)$ with distance estimate $\hat{d}=\text{{\sf dist}}$ defined in (4.5) with

[TABLE]

for any $\eta>0$, $\delta\in(0,1)$ and $\eta^{\prime}=\eta-\Delta\in(0,|\lambda_{1}|^{4t}L^{2})$. By substituting the expressions for bias, meas, and $\alpha$ into Lemma 6.1, it follows that for any $\delta^{\prime}\in(0,1)$,

[TABLE]

with probability at least

[TABLE]

By selecting

[TABLE]

it follows by the conditions $\psi>0$ and $\kappa<\frac{1}{2}$ that

[TABLE]

By choosing $\delta^{\prime}=n^{-(\kappa-\psi)/2}$ such that $\delta^{\prime}=O(\sqrt{\eta})$ and $\delta^{\prime 2}p(\sqrt{\eta-\Delta}n)^{3}=\Omega(n^{\frac{3}{4}})=\Omega(n^{2\psi})$ because $\psi<\frac{3}{8}$. Therefore, by substituting into (6.14), it follows that for any given $(u,v,w)\in[n]^{3}$, with probability 1-O(n\exp(-\Theta(n^{2\psi}))\Big{)}-O(n^{-5}),

[TABLE]

where the bias between $f_{r}$ and $f$ is bounded by $\boldsymbol{\varepsilon}$, and dominated by the bound between $\hat{F}$ and $f_{r}$. The final result follows from a union bound over $(u,v,w)\in[n]^{3}$.

The bound on MSE also follows by substituting $\delta=\frac{1}{2}$ and the same choice of $\eta$ from (6.15) into Lemma 6.1, and again noting that the bias between $F$ and $F_{r}$ is dominated by the error between $\hat{F}$ and $F_{r}$ such that

[TABLE]

This completes the proof of Theorem 5.2.

6.3 Proof of Corollary 5.3

Proof.

The proof follows the same format as the proof of Theorem 5.1. Let us denote the set of anchor vertices as $\mathcal{Y}$ such that $|\mathcal{Y}|=y$, and they are assumed to be chosen uniformly at random amongst all vertices. For a pair of vertices $(a,b)\in\mathcal{Y}^{2}$, the estimate $\hat{F}(a,b)$ follows the same exact computation as described in Section 4.1. As a result it follows from Theorem 5.1 that with high probability,

[TABLE]

Next we need to show the error is not degraded for non-anchor vertices $(u,v,w)\in([n]\setminus\mathcal{Y})^{3}$. Let $\zeta:[n]\to\mathcal{Y}$ denote the function that maps from each vertex to the closest anchor vertex as determined by the true distances $d$,

[TABLE]

and let $\hat{\zeta}:[n]\to\mathcal{Y}$ denote the data-dependent function that maps from each vertex to the closest anchor vertex as determined by the computed distances $\hat{d}$,

[TABLE]

The estimate for non-anchor vertices is then taken to be the estimate computed for the corresponding closest anchor vertices,

[TABLE]

such that

[TABLE]

By Theorem 5.1, as $(\hat{\zeta}(u),\hat{\zeta}(v),\hat{\zeta}(w))\in\mathcal{Y}^{3}$, the first term is bounded by $O(n^{-(\kappa-\psi)/2})$ with high probability. By property 6.1,

[TABLE]

The modified algorithm computes distances using Step 4 of the described algorithm between all pairs of anchor vertices, as well as all pairs $(u,a)$ such that $u\in[n]$ and $a\in\mathcal{Y}$. For each computed distance between a pair $(u,a)$, by Lemma 6.2, property 6.2 holds for $\Delta=\Theta(n^{-\kappa+\psi})$ with probability $1-\alpha_{1}$ where \alpha_{1}=O\Big{(}\exp(-n^{2\psi}(1-o(1)))\Big{)} for $\psi\in(0,\min(\kappa,\frac{3}{8}))$ and $\kappa\in(0,\frac{1}{2})$.

In order to bound $\max_{u\in[n]}d(\theta_{u},\theta_{\hat{\zeta}(u)})$, we argue that for every $u\in[n]$, with high probability

[TABLE]

where (a) and (c) hold with high probability for $\Delta=\Theta(n^{-\kappa+\psi})$ as a result of property 6.2, and (b) and (d) follow from the definition of the functions $\hat{\zeta}$ and $\zeta$.

To bound $\min_{a\in\mathcal{A}}d(\theta_{u},\theta_{a})$, we use (6.7) from property 6.3 to show that for any $u\in[n]$, $\eta=\Theta(n^{-\kappa+\psi})$, and $y=\Omega((p^{2}n^{3})^{1/4})=\Omega(n^{\kappa/2})$

[TABLE]

As a result, the max entrywise error is bounded by $O(n^{-(\kappa-\psi)/2})$ with high probability, which can be used to show the MSE bound of $O(n^{-(\kappa-\psi)})$. ∎

7 Proving distance estimate is close

In this section we argue that the distance estimate as defined in (4.5) is close to an ideal distance as claimed in the Lemma 6.2.

7.1 Regular enough growth of breadth-first-search (BFS) tree

The distance estimation algorithm of interest constructs a specific BFS tree for each vertex $u\in[n]$ with respect to the bipartite graph between vertices $[n]$ and $\mathcal{V}_{A}$ where recall that $\mathcal{V}_{A}=\{(b,c)\in[n/2]^{2}~{}\text{s.t.}~{}b<c\}$. The BFS tree construction is done so that vertices at different levels do not share coordinates, i.e. if vertex $a\in[n]$ is visited in an earlier layer of the BFS tree, then no vertex corresponding to $(a,b)$ for any $b\in[n]$ can be visited subsequently. Similarly, if $(a,b)$ is visited in the BFS tree, then no subsequent vertices including either coordinates $a$ or $b$ can be visited. The restriction is placed across different depths, whereas pairs of vertices $(a,b)$ and $(a,c)$ can be visited in the same depth. Amongst various valid BFS trees, the algorithm chooses one arbitrarily (for example, see Figure 3(c)).

We recall some notations. Consider a valid BFS tree rooted at vertex $u\in[n]$ which respects the constraint that no coordinate is visited more than once. Recall that for any $s\geq 1$, $\mathcal{U}_{u,s}\subseteq\mathcal{V}_{A}$ denotes the set of vertices at depth $(2s-1)$ and $\mathcal{S}_{u,s}\subseteq[n]$ denotes the set of vertices at depth $2s$ of the BFS tree, $\mathcal{B}_{u,s}=\cup_{l\in\lceil s/2\rceil}\mathcal{U}_{u,l}\cup_{h\in\lfloor s/2\rfloor}\mathcal{S}_{u,h}$, $\mathcal{G}(\mathcal{B}_{u,s})$ denotes all the information corresponding to the subgraph restricted to the first $s$ layers of the BFS tree which includes $\mathcal{B}_{u,s}$, the latent variables $\{\theta_{a}\}_{a\in\mathcal{B}_{u,s}}$ and the edge weights $\{M_{1}(a,b,c)\}_{a,(b,c)\in\mathcal{B}_{u,s}}$. The vector $N_{u,s}\in[0,1]^{n}$ is such that the $a$-th coordinate is equal to the product of weights along the path from $u$ to $a$ in the BFS tree for $a\in\mathcal{S}_{u,s}$, and the vector $W_{u,s}\in[0,1]^{\mathcal{V}_{A}}$ is such that the $(b,c)$-th coordinate is equal to the product of weights along the path from $u$ to $(b,c)$ in the BFS tree for $(b,c)\in\mathcal{U}_{u,s}$. The normalized vectors are $\tilde{N}_{u,s}=N_{u,s}/|\mathcal{S}_{u,s}|$ and $\tilde{W}_{u,s}=W_{u,s}/|\mathcal{U}_{u,s}|$ for $u\in[n],~{}s\geq 1$.

In a valid BFS tree rooted at vertex $u$, $\pi_{u}(a)$ denotes the parent of $a\in[n]$, and $\pi_{u}(b,c)$ denotes the parent of $(b,c)\in\mathcal{V}_{A}$. The neighborhood vectors satisfy recursive relationship,

[TABLE]

with $N_{u,0}=e_{u}$. We state the following result regarding regularity in the growth of the BFS tree.

Lemma 7.1**.**

Let $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$. Let $t$ be as defined in (4.3). For a given $\delta\in(0,\frac{1}{2})$ and for any $u\in[n]$, with probability 1-O\Big{(}n\exp\big{(}-\Theta(n^{2\kappa})\big{)}\Big{)}, or all $s\in[t-1]$,

[TABLE]

for $s=t$,

[TABLE]

and for $s\in[t]$,

[TABLE]

The set of single coordinate vertices visited within depth $2t$ is $o(n)$,

[TABLE]

Proof.

First observe that if $t$ is as defined in (4.3) with $\kappa\in(0,\frac{1}{2})$, then

[TABLE]

such that

[TABLE]

Note that $t$ is constant with respect to $n$.

For any $s\in[t]$, we study the growth of $|\mathcal{S}_{u,s}|$ and $|\mathcal{U}_{u,s}|$ conditioned on $\mathcal{B}_{u,2s-1}\cup\mathcal{U}_{u,s}$ and $\mathcal{B}_{u,2(s-1)}\cup\mathcal{S}_{u,s-1}$ respectively. To that end, conditioned on the set $\mathcal{B}_{u,2s-1}$ and the set $\mathcal{U}_{u,s}$, any vertex $i\in[n]\setminus\mathcal{B}_{u,2s-1}$ is in $\mathcal{S}_{u,s}$ independently with probability $(1-(1-p)^{|\mathcal{U}_{u,s}|})$. Thus the number of vertices in $\mathcal{S}_{u,s}$ is distributed as a binomial random variable. By Chernoff’s bound,

[TABLE]

Similarly, conditioned on the sets $\mathcal{B}_{u,2(s-1)}$ and $\mathcal{S}_{u,s-1}$, the set of vertices in $\mathcal{U}_{u,s}$ is equivalent to the number of edges in a graph with vertices $[n/2]\setminus\mathcal{B}_{u,2(s-1)}$ and an edge between $(i,j)$ if there is some $h\in\mathcal{S}_{u,s-1}$ such that $(i,j,h)\in{\Omega_{1}}$. This is an Erdos-Renyi graph, as each edge is independent with probability $(1-(1-p)^{|\mathcal{S}_{u,s-1}|})$. By Chernoff’s bound,

[TABLE]

Let us define the events

[TABLE]

Since $p\in(0,1)$ and hence $1-(1-p)^{x}\leq px$ for all $x\geq 1$, we have that under events $\mathcal{A}^{1}_{u,s}(\delta)\cap\mathcal{A}^{2}_{u,s}(\delta)$,

[TABLE]

which together implies that conditioned on event \cap_{h=1}^{s}\big{(}\mathcal{A}^{1}_{u,h}(\delta)\cap\mathcal{A}^{2}_{u,h}(\delta)\big{)}, for all $s\in[t]$

[TABLE]

and

[TABLE]

Therefore, for any $s\in[t-1]$ such that $s\leq\frac{1}{4\kappa}$ by the definition of $t$,

[TABLE]

With a similar argument we can show that

[TABLE]

The last step follows from checking that when $\kappa\in[\frac{1}{4},\frac{1}{2})$, $t=1$ such that $n^{2\kappa t}=o(n)$, and when $\kappa\in(0,\frac{1}{4})$, from $t\leq\frac{1}{4\kappa}+1$, it follows such that $n^{2\kappa t}=O(n^{\frac{1}{2}+2\kappa})=o(n)$ as $\kappa<\frac{1}{4}$. Recall that we split the coordinates such that $\cup_{\ell=1}^{t}\mathcal{U}_{u,\ell}\subset\mathcal{V}_{A}$, and the coordinates represented in $(a,b)\in\mathcal{V}_{A}$ are such that $a\in[n/2]$ and $b\in[n/2]$. Therefore by (7.13),

[TABLE]

Using (7.12), we establish lower bounds on $|\mathcal{S}_{u,s}|$ and $|\mathcal{U}_{u,s}|$ next. Note that, for $p\in(0,1)$, $1-p\leq e^{-p}$ and for any $x\in(0,1)$, $e^{-x}\leq 1-x+x^{2}$. It follows that $1-(1-p)^{x}\geq px(1-px)$. For $s\in[t]$ we can show that

[TABLE]

For $s\in[t-1]$ we can show that

[TABLE]

and for $s=t$,

[TABLE]

Then for $s\in[t]$,

[TABLE]

for $s\in[t-1]$,

[TABLE]

and for $s=t$, $|\mathcal{S}_{u,t}|\geq(1-\delta)^{2t}2^{-3t-1}n^{2\kappa t}(1-o(1))$.

To conclude the proof of Lemma 7.1, we need to argue that \cap_{s=1}^{t}\big{(}\mathcal{A}^{1}_{u,s}(\delta)\cap\mathcal{A}^{2}_{u,s}(\delta)\big{)} holds with high probability. To that end,

[TABLE]

We bound the each of the two summation terms on the right hand side in the last inequality next. Using (7.6) and (7.10), we have

[TABLE]

Similarly, using (7.7) and (7.11), we have

[TABLE]

Putting it all together, we have that

[TABLE]

since $\kappa\in(0,\frac{1}{2})$. By union bound over all $u\in[n]$, we obtain the desired bound on the probability of error. This concludes the proof of Lemma 7.1. ∎

7.2 Concentration of Quadratic Form One

Let $\mathcal{A}^{3}_{u,t}(\delta)$ denote the event that (7.1) holds for all $s\in[t-1]$, (7.2) holds, (7.3) holds for all $s\in[t]$, and (7.4) holds. Lemma 7.1 established that this event holds with high probability. Conditioned on the event $\mathcal{A}^{3}_{u,t}(\delta)$, we prove that a specific quadratic form concentrates around its mean. This will be used as the key property to eventually establish that the distance estimates are a good approximation to the ideal distances.

Lemma 7.2**.**

Let $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$, $t$ as defined in (4.3), $\delta\in(0,\frac{1}{2})$, and $\psi\in(0,\kappa)$. For any $u\in[n]$, with probability $1-2\exp(-n^{2\psi}(1-o(1)))$,

[TABLE]

Proof.

Recall that conditioning on event $\mathcal{A}^{3}_{u,t}(\delta)$ simply imposes the restriction that the neighborhood of $u\in[n]$ grows at a specific rate. This event is independent from latent parameters $\{\theta_{a}\}_{a\in[n]}$, the precise entries in $\Omega_{1}$ as well as associated values, i.e. $M_{1}$.

Conditioned on $\mathcal{A}^{3}_{u,t}(\delta)$, let $\mathcal{F}_{u,s}$ for $0\leq s\leq 2t$ denote the sigma-algebra containing information about the latent parameters, edges and the values associated with nodes in the bipartite graph up to distance $s$ from $u$, i.e. nodes $\mathcal{S}_{u,h^{\prime}}$ for $h^{\prime}\leq\lfloor s/2\rfloor$, $\mathcal{U}_{u,h^{{}^{\prime\prime}}}$ for $h^{{}^{\prime\prime}}\leq\lceil s/2\rceil$, associated latent parameters as well as edges of $\Omega_{1}$. Specifically, $\mathcal{F}_{u,0}$ contains information about latent parameter $\theta_{u}$ associated with $u\in[n]$; $\mathcal{F}_{u,s}$ contains information about latent parameters $\cup_{h=1}^{\lfloor s/2\rfloor}\{\theta_{a}\}_{a\in\mathcal{S}_{u,h}}\cup_{h=1}^{\lceil s/2\rceil}\{\theta_{b},\theta_{c}\}_{(b,c)\in\mathcal{U}_{u,h}}$ and all the associated edges and observations. This implies that $\mathcal{F}_{u,0}\subset\mathcal{F}_{u,1}\subset\mathcal{F}_{u,2}$, etc.

Recall that $Q$ denotes the $r\times n$ matrix where $Q_{ka}=q_{k}(\theta_{a}),k\in[r],a\in[n]$. We modify the notation due to the sample splitting, and we let $\mathcal{Q}$ denote the $r\times\binom{n/2}{2}$ matrix where $\mathcal{Q}_{kb}=q_{k}(\theta_{b_{1}})q_{k}(\theta_{b_{2}})$ for some $b\in\mathcal{V}_{A}$ that represents the pair of coordinates $(b_{1},b_{2})$ for $b_{1}<b_{2}\in[n/2]$.

We shall consider a specific martingale sequence with respect to the filtration $\mathcal{F}_{u,s}$ that will help establish the desired concentration of $e_{k}^{T}Q\tilde{N}_{u,t}-e_{k}^{T}\Lambda^{2t}Qe_{u}$. For $1\leq s\leq 2t$, define

[TABLE]

Note that $\tilde{N}_{u,0}=e_{u}$, and $Y_{u,s}$ is measurable with respect to $\mathcal{F}_{u,s}$ because $e_{k}^{T}\Lambda^{2t-s}Q\tilde{N}_{u,s/2}$ and $e_{k}^{T}\Lambda^{2t-s}Q\tilde{N}_{u,(s+1)/2}$ only depend on observations in the BFS tree within depth $s$.

By Lemmas 7.4 and 7.5, it follows that $Y_{u,s}$ is martingale with respect to $\mathcal{F}_{u,s}$ for $1\leq s\leq t$, i.e.

[TABLE]

Furthermore, for properly chosen $\nu_{s}$ as specified in Lemmas 7.4 and 7.5,

[TABLE]

almost surely for any $\lambda\in\mathbb{R}$.

We can then apply Proposition 7.3 with any arbitrarily small $\alpha_{*}$ such that for any $x\geq 0$,

[TABLE]

where for $n$ large enough,

[TABLE]

and $1+4\pi\leq 16$.

For $\psi\in(0,\kappa)$, we choose $x=\frac{16\lambda_{k}^{2t-2}n^{\psi}}{(1-\delta)n^{\kappa}}=o(1)$, such that with probability $1-2\exp(-n^{2\psi}(1-o(1)))$,

[TABLE]

∎

We recall the following concentration inequality for Martingale difference sequence, cf. [46, Theorem 2.19]:

Proposition 7.3**.**

Let $\{D_{k},\mathcal{F}_{k}\}_{k\geq 1}$ be a martingale difference sequence such that $\mathbb{E}[e^{\lambda D_{k}}|\mathcal{F}_{k-1}]\leq e^{\lambda^{2}\nu_{k}^{2}/2}$ almost surely for all $\lambda\in\mathbb{R}$. Then for all $x>0$,

[TABLE]

Lemma 7.4**.**

For any $s\in[t]$,

[TABLE]

Let $\nu=\sqrt{\frac{Q(1+4\pi)}{2}}$, and

[TABLE]

For any $\lambda\in\mathbb{R}$,

[TABLE]

Proof.

Recall $\mathcal{F}_{2s-1}$ contains all information in the depth $2s-1$ neighborhood of vertex $u$. In particular this includes the vertex set

[TABLE]

the vertex latent variables $\{\theta_{i}\}_{i\in\mathcal{B}_{u,2s-1}}$ and the edges and corresponding weights. Let us additionally condition on the set $\mathcal{S}_{u,s}$. As $2s$ is even,

[TABLE]

Let us define

[TABLE]

The randomness in $X_{i}$ only depends on $\theta_{i},\epsilon_{a_{1}a_{2}i},\mathbb{I}_{(a=\pi(i))}$. Note that we already conditioned on $\theta_{a_{1}},\theta_{a_{2}}$ for $a\in\mathcal{U}_{u,s}\subset\mathcal{B}_{2s-1}$. $X_{i}$ is independent from $X_{j}$ because the vertices and edges are disjoint, and $\pi(i)$ is independent from $\pi(j)$ as different vertices are allowed to have the same (or different) parents. First we compute the mean of $X_{i}$ (conditioned on $i\in\mathcal{S}_{u,s}$). For any vertex $i\in\mathcal{S}_{u,s}$, it must have exactly one parent in $\mathcal{U}_{u,s}$ due to the BFS tree constraints. The parent is equally likely to be any vertex in $\mathcal{U}_{u,s}$ due to the symmetry in the randomly sampled observations. Because the additive noise terms are mean zero, the eigenfunctions are orthonormal, and $\pi(i)$ is equally likely to be any $a\in\mathcal{U}_{u,s}$,

[TABLE]

Furthermore, $|X_{i}|\leq B$ almost surely as we assumed $|q_{k}(\theta)|\leq B$. By Hoeffding’s inequality, it follows that

[TABLE]

If we condition on the event $\mathcal{A}^{3}_{u,t}(\delta)$,

[TABLE]

for $s\in[t]$. Therefore,

[TABLE]

We finish the proof by using Lemma A.2 with $c=2$ and $Q=\frac{B^{2}\lambda_{k}^{4t-4s}2^{3s}(1+o(1))}{(1-\delta)^{2s}n^{2\kappa s}}$. ∎

Lemma 7.5**.**

For any $s\in[t]$,

[TABLE]

Let $\nu=\sqrt{\frac{Q(1+16\pi)}{2}}$, and

[TABLE]

For any $\lambda\in\mathbb{R}$,

[TABLE]

Proof.

As $2s-1$ is odd,

[TABLE]

Recall $\mathcal{F}_{2s-2}$ contains all information in the depth $2s-2$ neighborhood of vertex $u$. In particular this includes the vertex set

[TABLE]

the vertex latent variables $\{\theta_{i}\}_{i\in\mathcal{B}_{u,2(s-1)}}$ and the edges and corresponding weights $\{M_{1}(i,a)\}_{i,a\in\mathcal{B}_{u,2(s-1)}}$. Consider

[TABLE]

where we define for $i=(i_{1},i_{2})\in\mathcal{U}_{u,s}$,

[TABLE]

Conditioned on $\mathcal{F}_{u,2(s-1)}$ the randomness in $X_{i}$ only depends on $\theta_{i_{1}},\theta_{i_{2}}$, $\epsilon_{\pi(i)i_{1}i_{2}}$, and $\mathbb{I}_{(v=\pi(i))}$. Conditioned on $\mathcal{U}_{u,s}$ and $\{\theta_{i_{1}},\theta_{i_{2}}\}_{i\in\mathcal{U}_{u,s}}$, the random variables $X_{i}$ are independent as $\epsilon_{\pi(i)i_{1}i_{2}}$ and $\mathbb{I}_{(v=\pi(i))}$ are independent. The parent of $i=(i_{1},i_{2})\in\mathcal{U}_{u,s}$ is equally likely to be any vertex in $\mathcal{S}_{u,s-1}$, and the parent of $i=(i_{1},i_{2})\in\mathcal{U}_{u,s}$ is independent from the parent of $j=(j_{1},j_{2})\in\mathcal{U}_{u,s}$ with $j\neq i$ as different vertices are allowed to have the same (or different) parent. First we compute the mean of $X_{i}$ conditioned on $i\in\mathcal{U}_{u,s}$ and $\theta_{i_{1}},\theta_{i_{2}}$. Because the additive noise terms are mean zero and the parent of $i$ is equally likely to be any $v\in\mathcal{S}_{u,s-1}$,

[TABLE]

Furthermore, $|X_{i}|\leq B^{2}$ almost surely as we assumed $|q_{k}(\theta)|\leq B$. By Hoeffding’s inequality, it follows that

[TABLE]

Next we consider concentration with respect to the random subset $\mathcal{U}_{u,s}$ out of the $\mathcal{V}_{A}\setminus\mathcal{B}_{u,2(s-1)}$ possible vertices. In particular we would like to argue that with high probability,

[TABLE]

Next, we formalize it. To that end, conditioned on the size $|\mathcal{U}_{u,s}|$, the set $\mathcal{U}_{u,s}$ is a uniform random sample of the possible set of vertices $\mathcal{V}_{A}\setminus\mathcal{B}_{u,2(s-1)}$. The above expression on the left is thus the mean of a random sample $\mathcal{U}_{u,s}$ without replacement from $\mathcal{V}_{A}\setminus\mathcal{B}_{u,2(s-1)}$. Due to negative dependence, it concentrates around its means no slower than assuming that they were a sample of the same size from the same population with replacement, cf. [47, Theorem 4]. Therefore, using

[TABLE]

we can apply Hoeffding’s inequality to argue that

[TABLE]

Finally, we need to account for the randomness in $\{\theta_{i_{1}},\theta_{i_{2}}\}_{i\in\mathcal{V}_{A}}$, arguing that with high probability

[TABLE]

To formalize this, we start by recalling that

[TABLE]

Let $n_{u,s}=|[n/2]\setminus\mathcal{B}_{u,2(s-1)}|$, then $|\mathcal{V}_{A}\setminus\mathcal{B}_{u,2(s-1)}|=\binom{|n_{u,s}|}{2}$. Then the above summation can be written as a pairwise U-statistic,

[TABLE]

where $g$ is a symmetric function and each term $g(\theta_{i_{1}},\theta_{i_{2}})$ is bounded in absolute value by $B^{2}$. Furthermore,

[TABLE]

by the orthogonality model assumption. Therefore, by Lemma A.3

[TABLE]

By putting together all calculations, it also follows that

[TABLE]

and for $z_{1},z_{2},z_{3}>0$, with probability at least

[TABLE]

it holds that

[TABLE]

since $\|N_{u,s-1}\|_{\infty}\leq 1$. Conditioned on $\mathcal{A}^{3}_{u,t}(\delta)$, $n_{u,s}=n/2(1-o(1))$, and

[TABLE]

As a result, for $z_{1}=z_{2}=z_{3}$, the expression in (7.2) and (7.2) asymptotically dominate the expression in (7.2). It follows that, with appropriate choice of $z_{1}=z_{2}=z_{3}$ in the above,

[TABLE]

We finish the proof by using Lemma A.2 with $c=4$ and $Q=\frac{72\lambda_{k}^{4t-4s+2}B^{4}2^{3s}(1+o(1))}{(1-\delta)^{2s-1}n^{\min\{1,\frac{1}{2}+\kappa(2s-1)\}}}$.

∎

7.3 Concentration of Quadratic Form Two

Lemma 7.2 suggests the following high probability events: for any $u\in[n],k\in[r]$, $t$ as defined in (4.3), i.e. t=\Big{\lceil}\frac{1}{4\kappa}\Big{\rceil}, $\delta\in(0,1)$, and

[TABLE]

define

[TABLE]

Now, we state a useful concentration that builds on the above condition holding. It will be useful step towards establishing Lemma 6.2.

Lemma 7.6**.**

Let $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$, $t$ as defined in (4.3), and $\delta\in(0,\frac{1}{2})$. For any $u,v\in[n]$, conditioned on \cap_{k=1}^{r}\big{(}\mathcal{A}^{4}_{u,k,t}(x,\delta)\cap\mathcal{A}^{4}_{v,k,t}(x,\delta)\big{)},

[TABLE]

Proof.

Proof of Lemma 7.6. Assuming event \cap_{k=1}^{r}\big{(}\mathcal{A}^{4}_{u,k,t}(x,\delta)\cap\mathcal{A}^{4}_{v,k,t}(x,\delta)\big{)} holds,

[TABLE]

In above, we have simply used the fact that for two vectors $a,b\in\mathbb{R}^{r}$, $a^{T}b=\sum_{k}a_{k}b_{k}=\sum_{k}(e_{k}^{T}a)(e_{k}^{T}b)$. Now, consider the first term on the right hand side of the last inequality. If $\cap_{k=1}^{r}\mathcal{A}^{4}_{u,k,t}(x,\delta)$ holds, then $|(e_{k}^{T}Q\tilde{N}_{u,t}-e_{k}^{T}\Lambda^{2t}Qe_{u})|\leq x$. And if $\cap_{k=1}^{r}\mathcal{A}^{4}_{v,k,t}(x,\delta)$ holds, then $|(e_{k}^{T}\Lambda^{2}Q\tilde{N}_{v,t}-e_{k}^{T}\Lambda^{2(t+1)}Qe_{v})|\leq\lambda_{k}^{2}x$. Similar application to other terms and the fact that $|e_{k}^{T}Qe_{u}|,|e_{k}^{T}Qe_{v}|\leq\|q_{k}(\cdot)\|_{\infty}\leq B$, we conclude that

[TABLE]

∎

7.4 Concentration of Quadratic Form Three

We establish a final concentration that will lead us to the proof of good distance function property. For any $u\in[n]$, define event

[TABLE]

Lemma 7.7**.**

Let $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$, $t$ as defined in (4.3), $\delta\in(0,\frac{1}{2})$, and

[TABLE]

Let $S\equiv S_{u,v,t}=[n]\backslash(\mathcal{B}_{u,2t}\cup\mathcal{B}_{v,2t}\cup[n/2])$. Then, under event $A^{\prime}_{u,v,t}(x,\delta)$,

[TABLE]

with probability at least $1-4\exp(-n^{2\psi}(1-o(1)))-O(n^{-6})$ with $\psi\in(0,\kappa)$.

Proof.

First, note that $A^{\prime}_{u,v,t}(x,\delta)$ includes events $\mathcal{A}^{3}_{u,t}(\delta)$ and $\mathcal{A}^{3}_{v,t}(\delta)$. This implies that $|S|=\frac{n}{2}-o(n)=\frac{n(1-o(1))}{2}$. Furthermore, it implies that $|\mathcal{S}_{u,t}|$ and $|\mathcal{S}_{v,t}|$ are both greater than or equal to $(1-\delta)^{2t}2^{-3t-1}n^{2\kappa t}(1-o(1))$. As a result,

[TABLE]

The asymptotic relationships follow from the choice of $t\geq\frac{1}{4\kappa}$, and the fact that $\delta$ and $t$ are both constants.

Recall that ${M_{B}}(a,(\alpha,\beta))=\mathbb{I}_{((a,\alpha,\beta)\in{\Omega_{1}})}(F(a,\alpha,\beta)+\epsilon_{a\alpha\beta})$ for

[TABLE]

There are 3 sources of randomness: the sampling of entries in ${\Omega_{1}}$, the observation noise terms $\epsilon_{a\alpha\beta}$, and the latent variables $\theta_{a},\theta_{\alpha},\theta_{\beta}$. Since we enforce that $\alpha$ and $\beta$ are in the complement of $\mathcal{B}_{u,2t}\cup\mathcal{B}_{v,2t}$, the sampling, observations, and latent variables involved in ${{M_{B}}}$ are independent from $N_{u,t}$ and $N_{v,t}$.

Let us define the quantity

[TABLE]

for $\phi=\lceil 16/(1-2\kappa)\rceil$ where recall that

[TABLE]

Trivially, due to this thresholding, $|\tilde{T}(\alpha,\beta)|\leq\phi^{2}$ such that $|\tilde{T}(\alpha,\beta)-\mathbb{E}[\tilde{T}(\alpha,\beta)]|\leq 2\phi^{2}$.

To begin with, $N_{u,t}(a)=0$ if $a\notin\mathcal{S}_{u,t}\subset\mathcal{B}_{u,2t}$ and $N_{v,t}(b)=0$ if $b\notin\mathcal{S}_{v,t}\subset\mathcal{B}_{v,2t}$. Further, conditioned on event $A^{\prime}_{u,v,t}(x,\delta)$, all the information associated with $\mathcal{B}_{u,2t}$ and $\mathcal{B}_{v,2t}$ is revealed; however, information about $[n]\backslash(\mathcal{B}_{u,2t}\cup\mathcal{B}_{v,2t})$ is not. Let $\mathcal{F}(u,v,t,x,\delta)$ denote all the information revealed such that event $A^{\prime}_{u,v,t}(x,\delta)$ holds.

Let’s prove concentration in two steps. In step one, we condition on $\mathcal{F}(u,v,t,x,\delta)$ and the latent variables $\{\theta_{i}\}_{i\in[n]}$. The sampling process (edges in ${\Omega_{1}}$) and the observation noise are independent for distinct pairs $(\alpha,\beta)$ and $(\alpha^{\prime},\beta^{\prime})$. As a result, $T(\alpha,\beta)$ and $T(\alpha^{\prime},\beta^{\prime})$ are conditionally independent as long as $\{\alpha,\beta\}\cap\{\alpha^{\prime},\beta^{\prime}\}\neq 2$, i.e. they are not the exact same pair. The correlations across $T(\alpha,\beta)$ and $T(\alpha^{\prime},\beta^{\prime})$ are due only to the latent variables if $\alpha,\beta,\alpha^{\prime},\beta^{\prime}$ share any values. We will bound the variance of $T(\alpha,\beta)$ in Lemma 7.8, and by combining it with the conditional independence property across $T(\alpha,\beta)$, it follows that (using notation $\mathcal{F}=\mathcal{F}(u,v,t,x,\delta)$)

[TABLE]

The variables $\tilde{T}(\alpha,\beta)$ are also independent across $(\alpha,\beta)$ conditioned on the latent variables $\{\theta_{i}\}_{i\in[n]}$, and their variance is bounded above by the corresponding variances of $T(\alpha,\beta)$. Using the boundedness of $\tilde{T}(\alpha,\beta)$, by applying Bernstein’s inequality with the choice of $z=2n^{\psi}\left(\binom{|S|}{2}p^{2}|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|\right)^{1/2}$ for $\psi\in(0,\kappa)$, it follows that

[TABLE]

The last equality arises from the observation that $t$ is chosen such that conditioned on $\mathcal{F}$, we can plug in (7.27) to show that for our choice of $z$, it holds that

[TABLE]

In Lemma 7.9, we will show a bound on $\mathbb{P}\left(|T(\alpha,\beta)|\geq\phi^{2}\right)$, which translates to a bound on $\mathbb{E}[\mathbb{I}_{|T(\alpha,\beta)|\geq\phi^{2}}(|T(\alpha,\beta)|-\phi^{2})~{}|~{}\mathcal{F}]$, which then upper bound the difference between the conditional expectations of $T$ and $\tilde{T}$ according to

[TABLE]

Using this bound from Lemma 7.9 along with the conditions from $\mathcal{F}$ that guarantee $|S|=\Theta(n)$ and naively $|\mathcal{S}_{u,t}\cup\mathcal{S}_{v,t}|=O(n)$, it follows that

[TABLE]

We choose $\phi=\lceil\frac{16}{1-2\kappa}\rceil\geq\frac{16}{1-2\kappa}$ so that this difference between the expectations of $T$ and $\tilde{T}$ is $O(n^{-6})$.

By plugging in our choice of $\phi$ into Lemma 7.9, it also follows that

[TABLE]

By combining (7.28), (7.29), and (7.30), with probability at least $1-2\exp(-n^{2\psi}(1-o(1)))-O(n^{-6})$,

[TABLE]

where the first term will dominate the second term.

Finally we want to show concentration of the following expression with respect to the latent variables,

[TABLE]

The expression can be written as a pairwise U-statistic,

[TABLE]

where $g$ is a symmetric function, and

[TABLE]

It follows by boundedness of entries in $F$ and the fact that $\|N_{u,t}\|_{\infty}\leq 1$ and $\|N_{u,t}\|_{0}=|\mathcal{S}_{u,t}|$, that $|g(\theta_{\alpha},\theta_{\beta})|\leq 1$ almost surely. Therefore, by Lemma A.3 and choosing $z=\sqrt{8}n^{\psi}|S|^{-1/2}$,

[TABLE]

The expected value with respect to the randomness in the latent variables is

[TABLE]

Furthermore,

[TABLE]

By combining (7.31), (7.32), and (7.33), it follows that conditioned on $\mathcal{F}(u,v,s,\ell,x,\delta)$, with probability

[TABLE]

it holds that

[TABLE]

Note that the third term is dominated by the first term as $|S|p=o(1)$. This completes the proof of Lemma 7.7.

∎

Lemma 7.8**.**

Let $\mathcal{F}=\mathcal{F}(u,v,t,x,\delta)$ denote all the information revealed such that event $A^{\prime}_{u,v,t}(x,\delta)$ holds.

[TABLE]

Proof.

To compute the variance of $T(\alpha,\beta)$ conditioned on $\mathcal{F},\{\theta_{i}\}_{i\in[n]}$, note that there is correlation in the terms within the sum of $T(\alpha,\beta)$ as there may be pairs $(a,b)$ and $(a^{\prime},b^{\prime})$ that share coordinates. In particular because the observation noise and sampling randomness for ${M_{B}}(a,b,c)$ is independent across different entries $(a,b,c)$, then conditioned on $\{\theta_{i}\}_{i\in[n]}$, for $a\neq b$ and $a^{\prime}\neq b^{\prime}$, if all four coordinates $\{a,b,a^{\prime},b^{\prime}\}$ are distinct,

[TABLE]

if $|\{a,b\}\cap\{a^{\prime},b^{\prime}\}|=2$, i.e. $(a^{\prime},b^{\prime})=(a,b)$ or $(a^{\prime},b^{\prime})=(b,a)$,

[TABLE]

and if $\{a,b\}\cup\{a^{\prime},b^{\prime}\}=\{x,y,z\}$ such that $\{a,b\}\cap\{a^{\prime},b^{\prime}\}=\{x\}$, then

[TABLE]

The inequalities follow from the property that every entry of ${M_{B}}$ has absolute value bounded by 1, and takes value 0 with probability $(1-p)$ in the event it is not observed.

We use this to expand the variance calculation, and use the properties that for every entry $a$, $|N_{u,t}(a)|\leq\mathbb{I}_{(a\in\mathcal{S}_{u,t})}$. We have dropped the conditioning notation due to the length of the expressions.

[TABLE]

The first term dominates because $p|\mathcal{S}_{u,t}|\leq pn=o(1)$ and $p|\mathcal{S}_{v,t}|=o(1)$. ∎

Lemma 7.9**.**

[TABLE]

As a result,

[TABLE]

and

[TABLE]

Proof.

Let us define

[TABLE]

Furthermore, because $|{M_{B}}(\cdot,\cdot,\cdot)|\leq 1$ and $\|N_{u,s}\|_{\infty}\leq 1$, for any $a,b\in[n]$, it follows that

[TABLE]

which implies

[TABLE]

Note that $|Z_{u}(\alpha,\beta)\cup Z_{v}(\alpha,\beta)|\sim\text{Binomial}(|\mathcal{S}_{u,t}\cup\mathcal{S}_{v,t}|,p)$. It follows then that

[TABLE]

where we used the fact that $|\mathcal{S}_{u,t}\cup\mathcal{S}_{v,t}|p=o(1)$.

We use the bound on the tail probabilities to show that

[TABLE]

∎

7.5 Proof of Lemma 6.2

Proof.

Now we are ready to bound the difference between $d(u,v)$ and $\hat{d}(u,v)$ for any $u,v\in[n]$. Recall,

[TABLE]

and according to (4.5),

[TABLE]

for $S\equiv S_{u,s,t}=n\setminus(\mathcal{B}_{u,t}\cup\mathcal{B}_{v,t}\cup[n/2])$ and

[TABLE]

By Lemma 7.1, event $\mathcal{A}^{3}_{u,t}(\delta)$ holds with probability at least 1-O\Big{(}n\exp\big{(}-\Theta(n^{2\kappa})\big{)}\Big{)}. By Lemmas 7.2 and 7.6, conditioned on $\mathcal{A}^{3}_{u,t}(\delta)$, for

[TABLE]

event $A^{\prime}_{u,v,t}(x,\delta)$ holds with probability at least $1-4r\exp(-n^{2\psi}(1-o(1)))$, implying

[TABLE]

By Lemma 7.7, conditioned on $A^{\prime}_{u,v,t}(x,\delta)$, with probability $1-4\exp(-n^{2\psi}(1-o(1)))$,

[TABLE]

where $|S_{u,s,t}|=\Theta(n)=\Omega(n^{2\kappa})$, by event $A^{\prime}_{u,v,t}(x,\delta)$ and $t\geq\frac{1}{4\kappa}$,

[TABLE]

To put it all together, for $\psi\in(0,\kappa)$, with probability at least

[TABLE]

it holds that

[TABLE]

This completes the proof of Lemma 6.2. ∎

8 Proof of Lemma 6.3: perturbation analysis of distance

We establish the proof of Lemma 6.3 here. To do so, we establish a perturbation property of dist here, which combined with Lemma 6.2 will result into the proof of Lemma 6.3.

We study the perturbation in the dist estimate when each noisy observed entry is arbitrarily perturbed. Specifically, for any $(u,v,w)\in[n]^{3}$, $M_{1}(u,v,w)$ is observed with probability $p$. If observed, according to (3.1), $M_{1}(u,v,w)=F(u,v,w)+\epsilon_{uvw}=F_{r}(u,v,w)+\boldsymbol{\varepsilon}_{uvw}+\epsilon_{uvw}$, where $F_{r}$ is the best rank $r$ approximation to $F$. This expression shows that we can interpret the deviation from a rank $r$ model as a deterministic perturbation of $\boldsymbol{\varepsilon}_{uvw}$, bounded in absolute value by $\boldsymbol{\varepsilon}$. Note that $\boldsymbol{\varepsilon}_{uvw}$ can be any arbitrary (or adversarial), unknown deterministic quantity satisfying $|\boldsymbol{\varepsilon}_{uvw}|\leq\boldsymbol{\varepsilon}$.

Lemma 8.1 provides a bound on the perturbation in the distance estimate, dist, that results from these entrywise perturbabtions of the observations.

Lemma 8.1**.**

Let $p=n^{-3/2+\kappa}$ for $\kappa\in(0,\frac{1}{2})$, $t$ as defined in (4.3), $\delta\in(0,\frac{1}{2})$, and

[TABLE]

For any $u\in[n]$, recall the event

[TABLE]

Let event $A^{\prime}_{u,v,t}(x,\delta)$ hold. Let each observed entry of $M_{1}$ be perturbed by adding arbitrary, deterministic quantity bounded by $\boldsymbol{\varepsilon}\geq 0$. Then for any $u,v\in[n]^{2}$, the distance estimate $\text{{\sf dist}}(u,v)$ is perturbed by at most $O(t\boldsymbol{\varepsilon}(1+\boldsymbol{\varepsilon})^{2t-1}+t^{2}\boldsymbol{\varepsilon}^{2}(1+\boldsymbol{\varepsilon})^{4t-2})$ with probability at least 1-\exp\big{(}-\Omega(n^{2\kappa})\big{)}-O(n^{-8}).

Proof.

Recall definition of dist in (4.5):

[TABLE]

We shall bound the perturbation on $Z_{uv}$. Similar bounds will follow for the other three terms which will conclude the main results. Our interest is in understanding how does $Z_{uv}$ change if each observed entry is changed by arbitrary quantity bounded by $\boldsymbol{\varepsilon}\geq 0$. This will induce a bound on the changes in $T_{uv}(\cdot,\cdot)$ which will help bound the change in $Z_{uv}$. By assumption (8.1), $A^{\prime}_{u,v,t}(x,\delta)$ holds. Conditioned on event $A^{\prime}_{u,v,t}(x,\delta)$, all the information associated with $\mathcal{B}_{u,2t}$ and $\mathcal{B}_{v,2t}$ is revealed; however, information about $[n]\backslash(\mathcal{B}_{u,2t}\cup\mathcal{B}_{v,2t})$ is not. Let $\mathcal{F}(u,v,t,x,\delta)$ denote all the information revealed such that event $A^{\prime}_{u,v,t}(x,\delta)$ holds.

Under $A^{\prime}_{u,v,t}(x,\delta)$, by definition $\mathcal{A}^{3}_{u,t}(\delta)$ and $\mathcal{A}^{3}_{v,t}(\delta)$ holds. This implies that for $S\equiv\mathcal{V}_{B}(u,v,t)$, $|S|=\frac{n}{2}-o(n)=\frac{n(1-o(1))}{2}$. Furthermore, it implies that $|\mathcal{S}_{u,t}|$ and $|\mathcal{S}_{v,t}|$ are both greater than or equal to $(1-\delta)^{2t}2^{-3t-1}n^{2\kappa t}(1-o(1))$. As shown in (7.27),

[TABLE]

results from the choice of $t\geq\frac{1}{4\kappa}$, and the fact that $\delta$ and $t$ are both constants.

For given $\alpha\neq\beta\in\mathcal{V}_{B}(u,v,t)$, $T_{uv}(\alpha,\beta)$ is summation over terms, indexed by $a\neq b\in[n]$, containing product $N_{u,t}(a)N_{v,t}(b){M_{B}}(a,(\alpha,\beta)){M_{B}}(b,(\alpha,\beta))$. Now $N_{u,t}(a)=0$ if $a\notin\mathcal{S}_{u,t}$, $N_{v,t}(b)=0$ if $b\notin\mathcal{S}_{v,t}$. For $a\in\mathcal{S}_{u,t}$, $N_{u,t}(a)$ is product of $2t$ terms, each bounded in absolute value by $1$: let $N_{u,t}(a)=\prod_{i=1}^{2t}w_{i}$ with $|w_{i}|\leq 1$ for all $i\leq 2t$. Let $\varepsilon_{i}$ be arbitrary, deterministic quantity added to $w_{i}$ with $|\varepsilon_{i}|\leq\boldsymbol{\varepsilon}$ for $i\leq 2t$. Then change in $N_{u,t}(a)$ is bounded as

[TABLE]

That is, $N_{u,t}(a)$ changes by at most $\Delta(t,\boldsymbol{\varepsilon})$. Similarly $N_{v,t}(b)$ changes by at most $\Delta(t,\boldsymbol{\varepsilon})$. Therefore, $N_{u,t}(a)N_{v,t}(b)$ can change at most by $O(\Delta(t,\boldsymbol{\varepsilon})+\Delta(t,\boldsymbol{\varepsilon})^{2})$.

By definition $|{M_{B}}(a,(\alpha,\beta))|,|{M_{B}}(b,(\alpha,\beta))|\leq 1$. Further, ${M_{B}}(a,(\alpha,\beta)){M_{B}}(b,(\alpha,\beta))\neq 0$ only if $\mathbb{I}_{((a,\alpha,\beta)\in{\Omega_{1}})}\mathbb{I}_{((b,\alpha,\beta)\in{\Omega_{1}})}=1$. Therefore, we can bound change in the term $N_{u,t}(a)N_{v,t}(b){M_{B}}(a,(\alpha,\beta)){M_{B}}(b,(\alpha,\beta))$ as $\mathbb{I}_{((a,\alpha,\beta)\in{\Omega_{1}})}\mathbb{I}_{((b,\alpha,\beta)\in{\Omega_{1}})}O(\Delta(t,\boldsymbol{\varepsilon})+\Delta(t,\boldsymbol{\varepsilon})^{2})$. Therefore, we can bound the change in $Z_{uv}$ by

[TABLE]

where

[TABLE]

To conclude the Lemma, it will be sufficient to argue that $\sum_{\alpha,\beta\in S}X_{\alpha\beta}=O(|S|^{2}p^{2}|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|)$ with high probability given $\mathcal{F}$. We use a similar argument as the proof of Lemma 7.7. Given $\mathcal{F}\equiv\mathcal{F}(u,v,t,x,\delta)$, $\{X_{\alpha\beta}\}_{\alpha,\beta\in S^{2}}$ are conditionally independent random variables. By the same argument as that in Lemma 7.9, it follows that

[TABLE]

We define $\tilde{X}_{\alpha\beta}=\min(X_{\alpha\beta},\phi^{2})$ for $\phi=\lceil 16/(1-2\kappa)\rceil$ so that $|\tilde{X}_{\alpha\beta}-\mathbb{E}[\tilde{X}_{\alpha\beta}]|\leq\phi^{2}$. By (8.6), and the choice of $\phi$ along with the conditions from $\mathcal{F}$ that guarantee $|S|=\Theta(n)$ and naively $|\mathcal{S}_{u,t}\cup\mathcal{S}_{v,t}|=O(n)$,

[TABLE]

By (8.7),

[TABLE]

By the same argument as that in Lemma 7.8, it follows that

[TABLE]

By Bernstein’s inequality, for $z=2n^{\psi}|S|p(|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|)^{1/2}$ for $\psi\in(0,\kappa)$,

[TABLE]

By (8.5), (8.8), (8.9), and (8.10), given $A^{\prime}_{u,v,t}(x,\delta)$ holds, with probability 1-\exp\big{(}-n^{2\kappa}(1-o(1))\big{)}-O(n^{-6}), the change in $Z_{uv}$ is bounded above by

[TABLE]

By (8.2), this choice of $2n^{\psi}|S|p(|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|)^{1/2})=o(|S|^{2}p^{2}|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|)$ for $\psi\in(0,\kappa)$. Finally we use the bound that $\mathbb{E}[\sum_{\alpha,\beta\in S}X_{\alpha\beta}~{}|~{}\mathcal{F}]=\binom{|S|}{2}p^{2}|\mathcal{S}_{u,t}||\mathcal{S}_{v,t}|$ to argue that with high probability the change in $Z_{uv}$ is bounded above by

[TABLE]

This completes the proof of Lemma 8.1.

∎

8.1 Completing proof of Lemma 6.3

Under the setup of Lemma 6.3, as argued in the proof of Lemma 6.2, $A^{\prime}_{u,v,t}(x,\delta)$, with appropriate choice of $x,\delta$ as considered in statement of Lemma 8.1, holds with probability at least $1-4r\exp(-n^{2\psi}(1-o(1)))$. And dist (without perturbation), is within $O\left(n^{-(\kappa-\psi)}\right)$ for any pair of $u,v\in[n]$. By Lemma 8.1, under event $A^{\prime}_{u,v,t}(x,\delta)$, the dist is further perturbed by $O\left(t\boldsymbol{\varepsilon}(1+\boldsymbol{\varepsilon})^{2t-1}+t^{2}\boldsymbol{\varepsilon}^{2}(1+\boldsymbol{\varepsilon})^{4t-2}\right)$ with probability at least 1-O\Big{(}\exp(-n^{2\psi}(1-o(1)))\Big{)}-O\Big{(}n^{-6}\Big{)}. Putting these together, we conclude the claim of Lemma 6.3.

Acknowledgements

We gratefully acknowledge funding from the NSF under grants CCF-1948256 and CNS-1955997. Christina Lee Yu is also supported by an Intel Rising Stars Award.

Appendix A Useful Lemmas and Omitted Proofs

We present the Proof of Lemma 6.1 below.

Proof.

[Lemma 6.1] We assumed the algorithm has access to two fresh samples, where $M_{1}$ is used to compute $\hat{d}$, and ${M_{2}}$ is used to compute the final estimate $\hat{F}$. Alternatively one could effectively obtain two sample sets by sample splitting. For some $(a,b,c)\in{\Omega_{2}}$, the observation ${M_{2}}(a,b,c)$ is independent of $\hat{d}$, and $\mathbb{E}[{M_{2}}(a,b,c)]=f(\theta_{a},\theta_{b},\theta_{c})$. Conditioned on ${\Omega_{2}}$, by definition of $\hat{F}$ and by assuming properties 6.1 and 6.2, it follows that

[TABLE]

Inequality $(a)$ follows from property 6.1 and property 6.2 for all $3n$ tuples $\{(u,a):a\in[n]\}\cup\{(v,b):b\in[n]\}\cup\{(w,c):c\in[n]\}$: $|d(u,a)-\hat{d}(u,a)|\leq\Delta$ and $\hat{d}(u,a)\leq\eta\implies d(u,a)\leq\eta+\Delta$, similarly $d(v,b),d(w,c)\leq\eta+\Delta$. As per (3.1), we have that $\text{Var}[{M_{2}}(a,b,c)]\leq\sigma^{2}$ for all $(a,b,c)\in{\Omega_{2}}$. Define ${\mathcal{V}}_{uvw}=\{(a,b,c)\in[n]^{3}~{}:{d}(u,a)<\eta-\Delta,~{}{d}(v,b)<\eta-\Delta,~{}{d}(w,c)<\eta-\Delta\}$. Assuming property 6.3,

[TABLE]

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

[TABLE]

Therefore, by assuming property 6.2 for $3n$ tuples $\{(u,a):a\in[n]\}\cup\{(v,b):b\in[n]\}\cup\{(w,c):c\in[n]\}$, it follows that with probability at least $1-\exp\left(-\frac{\delta^{2}p\left(\text{{\sf meas}}(\eta-\Delta)n\right)^{3}}{2}\right)$,

[TABLE]

Define the event $\mathcal{H}=\{|{\Omega_{2}}_{uvw}|\geq(1-\delta)p\left(\text{{\sf meas}}(\eta-\Delta)n\right)^{3}|\}$. It follows that

[TABLE]

By definition, $F(u,v,w)=f(\theta_{u},\theta_{v},\theta_{w})\in[0,1]$ for all $u,v,w\in[n]$. Therefore,

[TABLE]

We add an additional $3n\alpha_{1}+\alpha_{2}$ in the final MSE bound: $3n\alpha_{1}$ for violation of property 6.2 for any of the $3n$ tuples $\{(u,a):a\in[n]\}\cup\{(v,b):b\in[n]\}\cup\{(w,c):c\in[n]\}$, and $\alpha_{2}$ for violation of property 6.3.

To obtain the high-probability bound, note that ${M_{2}}(a,b,c)$ are independent across indices $(a,b,c)\in{\Omega_{2}}$ as well as independent of observations in $\Omega_{1}$. Additionally, the model assumes that $F(a,b,c),{M_{2}}(a,b,c)\in[0,1]$, and $\mathbb{E}[{M_{2}}(a,b,c)]=F(a,b,c)$ for observed tuples $(a,b,c)$. By an application of Hoeffding’s inequality for bounded, zero-mean independent variables, for any $\delta^{\prime}\in(0,1)$ it follows that assuming property 6.1, property 6.2 for $3n$ tuples $\{(u,a):a\in[n]\}\cup\{(v,b):b\in[n]\}\cup\{(w,c):c\in[n]\}$, and property 6.3 hold, we have

[TABLE]

Therefore,

[TABLE]

with probability at least

[TABLE]

This completes the proof of Lemma 6.1. ∎

Lemma A.1**.**

The following inequalities hold:

• (a)

For any $\rho\geq 2$ and integer $r\geq 1$,

[TABLE]

• (b)

For any $\rho\geq 2$ and non-negative integer $s$,

[TABLE]

• (c)

Further, if $\exp(-a\rho)\leq\frac{1}{2}$ for some $a>0$, then

[TABLE]

Proof.

To prove (a), note that for any $\rho\geq 2$,

[TABLE]

To prove (b), first check that it trivially holds for $s=0$ and $s=1$. The inequality holds for $s=2$ iff $\rho\geq 2$. The inequality hold for $s$ iff $\rho\geq s^{1/(s-1)}$. We can verify that $s^{1/(s-1)}$ is a decreasing function in $s$, such that if the inequality holds for $s=2$, it will also hold for $s\geq 2$. To prove (c), further consider $\exp(-a\rho)\leq\frac{1}{2}$,

[TABLE]

∎

Lemma A.2**.**

If $\mathbb{P}\left(|X|\geq z\right)\leq c\exp(-\frac{z^{2}}{Q})$, then for all $\lambda\in\mathbb{R}$,

[TABLE]

with $\nu=\sqrt{\frac{Q(1+c^{2}\pi)}{2}}$.

Proof.

[TABLE]

Using the fact that $\sqrt{x}<e^{x/5}$ for all $x\geq 0$, it follows that

[TABLE]

Therefore, for all $\lambda\in\mathbb{R}$,

[TABLE]

∎

Lemma A.3**.**

Let $X_{1},\dots,X_{n}$ be i.i.d. random variables taking values in $\mathcal{X}$. Let $g:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ be a symmetric function. Consider U-statistics with respect to $g$ of $X_{1},\dots,X_{n}$ defined as

[TABLE]

Let $\|g\|_{\infty}\leq b$ for some $b>0$. Then,

[TABLE]

The proof follows directly from an implication of Azuma-Hoeffding’s inequality. For example, see [46, Example 2.23] for a proof.