Nearest Neighbors for Matrix Estimation Interpreted as Blind Regression for Latent Variable Model

TL;DR

This paper introduces a nearest-neighbor-based algorithm for blind matrix regression that effectively estimates missing entries using latent features, with proven consistency and competitive performance on real datasets and tensor completion tasks.

Contribution

The paper presents a novel nearest-neighbor approach for blind matrix estimation, providing theoretical guarantees and demonstrating practical effectiveness on real-world data and tensor completion.

Findings

Algorithm is consistent under Lipschitz conditions.

Performs better than basic collaborative filtering.

Competitive with state-of-the-art tensor completion methods.

Abstract

We consider the setup of nonparametric {\em blind regression} for estimating the entries of a large $m \times n$ matrix, when provided with a small, random fraction of noisy measurements. We assume that all rows $u \in [m]$ and columns $i \in [n]$ of the matrix are associated to latent features $x_{\text{row}}(u)$ and $x_{\text{col}}(i)$ respectively, and the $(u,i)$-th entry of the matrix, $A(u, i)$ is equal to $f(x_{\text{row}}(u), x_{\text{col}}(i))$ for a latent function $f$. Given noisy observations of a small, random subset of the matrix entries, our goal is to estimate the unobserved entries of the matrix as well as to "de-noise" the observed entries. As the main result of this work, we introduce a nearest-neighbor-based estimation algorithm, and establish its consistency when the underlying latent function $f$ is Lipschitz, the underlying latent space is a bounded diameter Polish space, and the random fraction of observed entries in the matrix is at least $\max \left( m^{-1 + \delta}, n^{-1/2 + \delta} \right)$, for any $\delta > 0$. As an important byproduct, our analysis sheds light into the performance of the classical collaborative filtering algorithm for matrix completion, which has been widely utilized in practice. Experiments with the MovieLens and Netflix datasets suggest that our algorithm provides a principled improvement over basic collaborative filtering and is competitive with matrix factorization methods. Our algorithm has a natural extension to the setting of tensor completion via flattening the tensor to matrix. When applied to the setting of image in-painting, which is a $3$-order tensor, we find that our approach is competitive with respect to state-of-art tensor completion algorithms across benchmark images.

Equations

$$N_{\text{eff}}(\mathcal{X},r,\epsilon )\triangleq \underset{I,S\subset \mathcal{X}}{\inf }\left\{|I|\text{ s.t. }S\subset {\cup }_{x\in I}B(x,r)\text{ and }\mu (\mathcal{X}\setminus S)\le \epsilon \right\}.$$

$$Z(u,i)=\{\begin{array}{ll}A(u,i)+N(u,i)& \text{when }M(u,i)=1,\\ \text{unknown}& \text{when }M(u,i)=0.\end{array}$$

$$A(u,i)=f(x_{\text{row}}(u),x_{\text{col}}(i))$$

$$\displaystyle\big{|}f(\alpha_{1},\alpha_{1})-f(\beta_{1},\beta_{2})\big{|}$$

$$\displaystyle\qquad\leq L\max\big{(}d_{\mathcal{X}_{\textrm{row}}}(\alpha_{1},\beta_{1}),d_{\mathcal{X}_{\textrm{col}}}(\alpha_{2},\beta_{2})\big{)}.$$

$$M(u,i)=\{\begin{array}{ll}1& \text{with probability }p,\\ 0& \text{with probability }1-p.\end{array}$$

$$\text{MSE}(\hat{A})=\mathbb{E}\left[\frac{1}{mn}\sum _{u=1}^m\sum _{i=1}^n{\left(\hat{A}(u,i)-A(u,i)\right)}^2\right].$$

$$\mathbf{Z}(u,i)\stackrel{d}{=}\mathbf{Z}\left(\sigma (u),\tau (i)\right)\text{ for all }(u,i),$$

$$\displaystyle{\bf Z}\left(u,i\right)\overset{d}{=}f_{\theta}\big{(}\theta_{\textrm{row}}(u),\theta_{\textrm{col}}(i),\theta_{\text{entry}}(u,i)\big{)}~{}\text{ for all }(u,i)$$

$$\displaystyle Z(u,i)=f_{\theta}\big{(}\theta_{\textrm{row}}(u),\theta_{\textrm{col}}(i),\theta_{\text{entry}}(u,i)\big{)}.$$

$$Z(u,i)=f(x_{\text{row}}(u),x_{\text{col}}(i))+\eta (u,i)$$

$$Z(u,i)$$

$$=f(x_{\text{row}}(u),x_{\text{col}}(i))+\eta (u,i).$$

$${\text{dissim}}_{\text{row}}(u,v)=\frac{1}{|{\mathcal{B}}_{\text{adj}}(u)\cap {\mathcal{B}}_{\text{adj}}(v)|}$$

$$\displaystyle\qquad\qquad\times\sum_{j\in\mathcal{B}_{\text{adj}}(u)\cap\mathcal{B}_{\text{adj}}(v)}\big{(}Z(u,j)-Z(v,j)\big{)}^{2},$$

$${\mathcal{B}}_{\text{adj}}(v)\triangleq \{j\in [n]:(v,j)\in {\Omega }_1\},$$

$${\mathcal{B}}_{\text{est}}(u,i)$$

$$\displaystyle\quad\textrm{dissim}_{\text{row}}(u,v)\leq\eta~{}\text{and}~{}(v,j)\in\Omega_{2}\big{\}}.$$

$$\hat{A}(u,i)=\frac{1}{|{\mathcal{B}}_{\text{est}}(u,i)|}\sum _{(v,j)\in {\mathcal{B}}_{\text{est}}(u,i)}Z(v,j)$$

$$\displaystyle\mathbb{E}_{\textrm{cond}}\left[\big{(}\widehat{A}(u,i)-A(u,i)\big{)}^{2}\right]$$

$$\displaystyle=\mathbb{E}_{\textrm{cond}}\left[\bigg{(}\frac{1}{N_{\text{base}}}\sum_{(v,i)\in\mathcal{B}_{\text{est}}(u,i)}\big{[}A(v,i)-A(u,i)+N(v,i)\big{]}\bigg{)}^{2}\right]$$

$$\displaystyle\leq\max_{(v,i)\in\mathcal{B}_{\text{est}}(u,i)}\big{\|}f(x_{\text{row}}(u),\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}+\frac{2\sigma^{2}}{N_{\text{base}}(\eta)}.$$

$$\displaystyle\mathbb{E}\left[\big{(}Z(u,j)-Z(v,j)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(u),x_{\text{row}}(v)\right]$$

$$\displaystyle\qquad=\mathbb{E}\left[\big{(}A(u,j)-A(v,j)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(u),x_{\text{row}}(v)\right]$$

$$\displaystyle\qquad\qquad+\mathbb{E}\left[\big{(}N(u,j)-N(v,j)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(u),x_{\text{row}}(v)\right]$$

$$\displaystyle\qquad=\big{\|}f(x_{\text{row}}(u),\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}+2\sigma^{2}.$$

$$w_{ui}(v,j)=\{\begin{array}{ll}{e}^{-\lambda \text{dissim}((u,i),(v,j))}& \text{if }(v,j)\in {\mathcal{B}}_{\text{est}}(u,i),\\ 0& \text{else},\end{array}$$

$$\hat{A}(u,i)=\frac{{\sum }_{(v,j)}{w}_{ui}(v,j)Z(v,j)}{{\sum }_{(v,j)}{w}_{ui}(v,j)}.$$

$$\textrm{MSE}(\hat{A})=O\Bigg{(}(D_{f}+\sigma)^{2}\bigg{(}\sqrt{\frac{\log(m-1)}{(n-1)p^{2}}}\vee\frac{|\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})|}{(m-1)p}\bigg{)}\Bigg{)}.$$

$$\displaystyle\phi(x,r)=\mathbb{P}_{x_{\text{row}}(v)\sim\mu_{\mathcal{X}_{\textrm{row}}}}\left(\big{\|}f(x,\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}\leq r\right).$$

$$\text{MSE}(\hat{A})$$

$$\displaystyle=\mathbb{E}\left[\big{(}\widehat{A}(1,1)-A(1,1)\big{)}^{2}\right].$$

$$\displaystyle\mathbb{E}\left[\big{(}\widehat{A}(1,1)-A(1,1)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(1)\right]$$

$$\displaystyle\leq\big{(}\eta-2\sigma^{2}\big{)}+CK\bigg{[}\frac{1}{\sqrt{(n-1)p^{2}}}+\exp\big{(}-c(n-1)p^{2}\big{)}\bigg{]}$$

$$\displaystyle+2\mathbb{E}_{x_{\text{row}}(1)}\Bigg{[}\bigg{[}(m-1)p\cdot\phi\Big{(}x_{\text{row}}(1),\eta^{\prime}-2\sigma^{2}\Big{)}\bigg{]}^{-1}\Bigg{]}$$

$$\displaystyle+\big{(}D_{f}^{2}+1\big{)}$$

$$\displaystyle\qquad\times\mathbb{E}_{x_{\text{row}}(1)}\Bigg{[}\exp\bigg{(}-\frac{(m-1)p}{8}\cdot\phi\Big{(}x_{\text{row}}(1),\eta^{\prime}-2\sigma^{2}\Big{)}\bigg{)}\Bigg{]}$$

Full text

I Introduction

The problem of matrix completion has received enormous attention in the past decade: consider an $m\times n$ matrix $A$ of interest. Suppose we observe a subset of the entries of an $m\times n$ matrix $Z$, which is a noisy version of $A$, such that each $(u,i)$-th entry $Z(u,i)$ is a random variable with $\mathbb{E}[Z(u,i)]=A(u,i)$ for $u\in[m],~{}i\in[n]$111We shall utilize notation $[m]=\{1,\dots m\}$.. The goal of matrix completion is to recover matrix $A$ given partial observations from $Z$.

I-A Our Contributions

We provide a similarity-based nearest neighbor algorithm akin to popular collaborative filtering with theoretical performance guarantees under the latent variable model. In addition, we extend the analysis of our algorithm to the setting of tensor completion through flattering tensor to matrix. To our knowledge, this is the first theoretical analysis for similarity-based collaborative filtering algorithms, shedding insight into the widespread success of this popular heuristic for the past two decades. The algorithm we introduce is a simple variant of classical collaborative filtering, in which we compute similarities between pairs of rows and pairs of columns by comparing their common overlapped entries. Our model assumes that each row and column is associated to a latent variable, i.e. a vector of hidden features, and that the data entry is in expectation equal to some unknown function of those latent variables. We assume the latent space is a complete, separable metric space aka Polish space equipped with a Borel probability measure. The key regularity condition that we require is that the image of the latent space by the latent function has a small effective covering number with respect to the push-forward measure222Let $(\mathcal{X},d,\mu)$ denote the latent space. Let $B(x,r)=\{z\in\mathcal{X}:d(x,z)\leq r\}$. Given $\varepsilon,r>0$, the effective covering number of $\mathcal{X}$ with respect to $\mu$

$N_{\textrm{eff}}(\mathcal{X},r,\varepsilon)\triangleq\inf_{I,S\subset\mathcal{X}}\left\{|I|\text{ s.t. }S\subset\cup_{x\in I}B(x,r)\text{ and }\mu(\mathcal{X}\setminus S)\leq\varepsilon\right\}.$

which Borel measure over Polish space naturally satisfies.

Given this latent variable model, we prove that the estimate produced by this algorithm is consistent as long as the fraction of entries that are observed is at least $\max(m^{-1+\delta},n^{-1/2+\delta})$ for some $\delta>0$ for an $m\times n$ matrix (for precise statement, See Corollary 3 and its implication). We provide experiments using our method to predict ratings in the MovieLens and Netflix datasets. The results suggest that our algorithm improves over basic collaborative filtering and is competitive with factorization-based methods.

We also discuss that the algorithm and analysis can be extended to tensor completion by flattening the tensor to a matrix. We implemented our method for predicting missing pixels in image in-painting, which showed that our method is competitive with existing spectral methods used for tensor completion.

The algorithm that we propose has similarities to classical non-parametric nearest neighbor method, cf. [1] and kernel regression, which also relies on approximations by local smoothing, cf. [2, 3]. However, since kernel regression and other similar methods use explicit knowledge of the input features, their analysis and proof techniques do not extend to our context. Instead of using distance in the unknown latent space, the algorithm weighs datapoints according to similarities that are computed from the data itself. Our analysis shows that although the similarities between the data points may not reflect the distance between the latent features, they essentially reflect the functional distances (in the $L^{2}$ sense) between the latent function restricted to the pair of rows (or columns) associated with the data points, which is sufficient to guarantee that the datapoints with high similarities are indeed similar in value.

I-B Related Literature

The primary methods used to solve the problem in the literature include neighbor-based approaches, such as collaborative filtering, and spectral approaches, which include low-rank matrix factorization or minimization of a loss function with respect to spectral constraints.

Spectral Methods

In the recent years, there have been exciting intellectual developments in the context of spectral approaches such as matrix factorization. All matrices admit a singular-value decomposition, such that they can be uniquely factorized. The goal of the factorization-based method is to recover row and column singular vectors accurately from the partially observed, noisy matrix $Z$ and subsequently estimate the matrix $A$. [4] was one of the earliest works to suggest the use of low-rank matrix approximation in this context. Subsequently, statistically efficient approaches were suggested using optimization-based estimators, proving that matrix factorization can fill in the missing entries with sample complexity as low as $rm\log m$ for an $m\times m$ matrix, where $r$ is the rank of the matrix [5, 6, 7, 8, 9]. There has been an exciting line of ongoing work to make the resulting algorithms faster and scalable [10, 11, 12, 13, 14, 15].

[16] proposed a spectral clustering method for inferring the edge label distribution for a network sampled from a generalized stochastic block model. The model is similar to the proposed latent variable model introduced in Section II, except that the edges are labeled by one of finitely many labels in a symmetric setup with $m=n$, and the goal is to estimate the label distribution in addition to the expected label. When the expected function has a finite spectrum decomposition, i.e. low rank, then they provide a consistent estimator for the sparse data regime, with $\Omega(m\log m)$ samples. When the function is only approximately low rank (e.g. the class of general Lipschitz functions), for a fixed rank $r$ approximation, the mean squared error bounds converge to a positive constant which captures the low rank approximation gap. That is, $\Omega(m\log m)$ samples are not sufficient to guarantee consistent estimation for the entire class of Lipschitz functions.

Many of these approaches are based on the structural assumption that the underlying matrix is low-rank and the matrix entries are reasonably “incoherent”. Unfortunately, the low-rank assumption may not hold in practice. The recent work [17] makes precisely this observation, showing that a simple non-linear, monotonic transformation of a low-rank matrix could easily produce an effectively high-rank matrix, despite few free model parameters. They provide an algorithm and analysis specific to the form of their model, which achieves sample complexity of $O((mn)^{2/3})$ for an $m\times n$ matrix. However, their algorithm only applies to functions $f$ which are a nonlinear monotonic transformation of the inner product of the latent features. [18] propose an algorithm for estimating locally low rank matrices, however their algorithm assumes prior knowledge of the “correct” kernel function between pairs of rows and columns which is not known a priori.

[19] proposes the universal singular value thresholding estimator (USVT) inspired by low-rank matrix approximation. Somewhat interestingly, it argues that under the latent variable model considered in this work (see Section II), the USVT algorithm provides an accurate estimate for any Lipschitz function. However, to guarantee consistency of the USVT estimator for an $m\times m$ (i.e. $m=n$) matrix, it requires observing $\Omega\left(m^{\frac{2(d+1)}{(d+2)}}\right)$ many entries out of the $m^{2}$ total entries, where $d$ is the dimension of the latent space in which the row and column latent features belong. In recent work, [20] extends the analysis of USVT for graphon estimation, assuming a generative latent variable model for binary observation matrices representing networks. If the latent function is $\alpha$-Holder smooth, author establishes that the spectrum decays polynomially, and thus the MSE of the USVT estimator is bounded above by O\big{(}(mp)^{-\frac{2\alpha}{2\alpha+d}}\big{)}, which converges to zero as long as $p=\omega(m^{-1})$ (for an $m\times m$ symmetric matrix).

Collaborative Filtering

The term collaborative filtering was coined by [21], and this technique is widely used in practice due to its simplicity and ability to scale. There are two main paradigms in neighborhood-based collaborative filtering: the user-user paradigm and the item-item paradigm. To recommend items to a user in the user-user paradigm, one first looks for similar users, and then recommends items liked by those similar users. In the item-item paradigm, in contrast, items similar to those liked by the user are found and subsequently recommended. Much empirical evidence exists that the item-item paradigm performs well in many cases [22, 23, 24]. There have also been many heuristic improvements upon the basic algorithm, such as normalizing the data, combining neighbor methods with spectral methods, combining both user and item neighbors, and additionally optimizing over interpolation weights given to each datapoint within the neighborhood when computing the final prediction [25, 26, 27].

Despite the widespread success of similarity-based collaborative filtering heuristics, the theoretical understanding of these method is very limited. In recent works, latent mixture models have been introduced to explain the collaborative filtering algorithm as well as the empirically observed superior performance of item-item paradigms, c.f. [28, 29]. However, these results assume binary ratings and a specific parametric model, such as a mixture distribution model for preferences across users and movies. We hope that by providing an analysis for collaborative filtering within a nonparametric model, we can provide a better understanding of collaborative filtering.

Within the context of dense graphon estimation, when the entries in the data matrix are binary and the sample probability $p$ is constant $\Theta(1)$, there have been a few theoretical results that prove convergence of the mean squared error for similarity-based methods [30, 31]. They hinge upon computing similarities between rows or columns by comparing commonly observed entries, similar to collaborative filtering. Similar to our result, they are able to prove convergence for the class of Lipschitz functions. However, [30] assumes that the algorithm is given multiple instances of the sampled dataset, which is not available in our formulation. [31] is weaker than our result in that it assumes $p=\Theta(1)$, however they are able to handle a more general noise model, when the entries are binary. The similarity between a pair of vertices is computed from the maximum difference between entries in the associated rows of the second power of the data matrix, which is computationally more expensive than directly comparing rows in the original data matrix.

Tensor Completion

Recently there have been efforts to extend decomposition methods or neighborhood-based approaches to the context of tensor completion, however this has proven to be significantly more challenging than matrix completion due to the complication that tensors do not have a canonical decomposition such as the singular value decomposition (SVD) for a matrix. This property makes obtaining a decomposition for a tensor challenging. The survey [32] elaborates on these challenges. There have been recent developments in obtaining efficient tensor decompositions in form of rank-1 tensors (tensors obtained from one vector), presented in [33]. This has been especially effective in learning latent variable models and estimating missing data as shown in, for example [34, 35].

Many results in tensor estimation take the approach of flattening the tensor to a matrix and subsequently apply matrix estimation algorithms [36, 37, 38, 39]. A $d$-order tensor where each dimension is length $n$ would be flattened to a $n^{\lfloor d/2\rfloor}\times n^{\lceil d/2\rceil}$ matrix, resulting in a sample complexity of $\Omega(n^{\lceil d/2\rceil}\text{polylog}(n))$. Subsequently there has been a line of work extending spectral methods, local iterative methods, or the sum of squares method to the specific tensor structure to obtain improved sample complexities of $\Omega(n^{d/2}\text{polylog}(n))$ [34, 40, 41, 42, 43, 44, 45].

Beyond tensor decomposition, there have been recent developments in the context of learning latent variable models or mixture distributions also called non-negative matrix factorization, c.f. [46, 47].

I-C Organization of the Paper

In Section II we setup the formal model and problem statement and discuss the assumptions needed for our analysis. In Section III we introduce the basic form of our algorithm, which is similar to the user-user variant of collaborative filtering. We present heuristic variants of the algorithm that perform well in practice. In Section IV we present the main theoretical results of our paper as they pertain to matrix completion, showing provable convergence of the user-user (and by symmetry item-item) variant of our algorithm. In Section V we provide a discussion of our results and comparison with other models. In Section VI we discuss how to extend the algorithm and analysis to tensor completion. In Section VII we present experimental results from applying our methods to both matrix completion in the context of predicting movie ratings, and tensor completion in the context of image inpainting. The detailed proofs are presented in the Appendix.

II Setup

II-A Our Model

Suppose that there is an unknown $m\times n$ matrix $A$ which we would like to estimate. We observe only a fraction of the total $mn$ entries of $A$ with some noise added. Let $\Omega\subset[m]\times[n]$ denote the index set of observed entries. Specifically, we observe entries of data matrix $Z$ that is generated as follows. Let $M\in\{0,1\}^{m\times n}$ be a binary matrix, which we call the masking matrix. We let $A\in\mathbb{R}^{m\times n}$ denote the signal matrix and $N\in\mathbb{R}^{m\times n}$ denote the noise matrix. For each $(u,i)\in[m]\times[n]$,

[TABLE]

For later use, we let $\Omega\triangleq\{(u,i)\in[m]\times[n]:M(u,i)=1\}$ denote the set of index pairs of the observed entries.

II-A1 Latent Variable Model

We assume $A(u,i)$ is generated by the following model.

• •

Nonparametric model: there exists a latent function $f$ such that

[TABLE]

for all $(u,i)\in[m]\times[n]$. Here, $x_{\text{row}}(u),x_{\text{col}}(i)$ denote latent variables associated with row $u$ and column $i$, respectively.

• •

Regularity Assumptions:

• –

For all $u\in[m]$, $x_{\text{row}}(u)\in\mathcal{X}_{\textrm{row}}$, where $(\mathcal{X}_{\textrm{row}},d_{\mathcal{X}_{\textrm{row}}})$ is a complete, separable metric space a.k.a. Polish space equipped with a Borel probability measure $\mu_{\mathcal{X}_{\textrm{row}}}$ and $x_{\text{row}}(u)$ is drawn i.i.d. according to $\mu_{\mathcal{X}_{\textrm{row}}}$.

• –

For all $i\in[n]$, $x_{\text{col}}(i)\in\mathcal{X}_{\textrm{col}}$, where $(\mathcal{X}_{\textrm{col}},d_{\mathcal{X}_{\textrm{col}}})$ is a complete, separable metric space aka Polish space equipped with a Borel probability measure $\mu_{\mathcal{X}_{\textrm{col}}}$ and $x_{\text{col}}(i)$ is drawn i.i.d. according to $\mu_{\mathcal{X}_{\textrm{col}}}$.

• –

Latent function $f$ is bounded, i.e. for all $\alpha\in\mathcal{X}_{\textrm{row}}$ and $\beta\in\mathcal{X}_{\textrm{col}}$, \big{|}f(\alpha,\beta)\big{|}\leq D_{f}.

• –

Without loss of generality, we shall assume that there exists $\alpha\in\mathcal{X}_{\textrm{row}}$ and $\beta\in\mathcal{X}_{\textrm{col}}$ such that $f(\alpha,\beta)=0$.

Our results will additionally require that the local measure is well-behaved such that there is a sufficient mass of nearest neighbors. In particular we introduce two concrete models that have good local neighborhood properties:

1. (a)

Finite Types: Let $\mathcal{X}_{\textrm{row}}$ be equipped with the discrete metric333$d_{\mathcal{X}_{\textrm{row}}}(x_{1},x_{2})=1$ if and only if $x_{1}\neq x_{2}$. topology and suppose $\mu_{\mathcal{X}_{\textrm{row}}}$ has finite support in $\mathcal{X}_{\textrm{row}}$ with $\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})$ denoting the support of $\mu_{\mathcal{X}_{\textrm{row}}}$. 2. (b)

Lipschitz Latent Function: Assume that the function $f$ is $L$-Lipschitz in the sense that

[TABLE]

II-A2 Noise

We assume $N(u,i)$ is a centered, sub-gaussian random variable for all $(u,i)\in[m]\times[n]$ such that

• •

Centered: $\mathbb{E}N(u,i)=0$ for all $(u,i)\in[m]\times[n]$.

• •

Sub-gaussian444The Orlicz $\psi_{2}$-norm of a random variable $X$ is defined as \|X\|_{\psi_{2}}\triangleq\inf_{t>0}\{\mathbb{E}\left[\exp\big{(}\frac{|X|^{2}}{t^{2}}\big{)}\right]\leq 2\}.: there exists $\sigma>0$ such that $\|N(u,i)\|_{\psi_{2}}\leq\sigma$ for all $(u,i)\in[m]\times[n]$.

• •

Independent and identically distributed (i.i.d.): $N(u,i)$’s are independent and identically distributed.

II-A3 Masking

We assume $M$ is a random matrix with each entry drawn as per $\textrm{Bernoulli}(p)$ for some $p\in(0,1]$, i.i.d. That is, for each $(u,i)\in[m]\times[n]$,

[TABLE]

Our algorithm will involve one threshold parameter which will need to be tuned. As a result of our analysis we can specify the optimal choice of the threshold as a function of $n,m,p,\sigma,L,D_{f}$ which trades off between the bias and variance. In practice, all of these parameters are assumed to be unknown and the threshold parameter can be chosen via cross-validation.

II-B Problem Statement: Blind Regression

Note that the latent variable model representation in (2) is not the unique representation, as there exists multiple equivalent representations. Suppose that one applies a measure-preserving transformation $T$ on the latent feature space $\mathcal{X}_{\textrm{row}}\times\mathcal{X}_{\textrm{col}}$, and take the push-forward of $f$ with respect to $T$ ($f\circ T^{-1}$) as the new latent function. This new representation – the pair of the latent space and the latent function – yields the same data generation process. Therefore, the question of estimating the function $f$ itself is not well posed, and we focus our energy on predicting the values $A(u,i)$.

Problem 1** (Blind Regression).**

Given $Z\in\mathbb{R}^{m\times n}$ that is partially observed on $\Omega\subset[m]\times[n]$ as described in Section II-A, we want to estimate the underlying matrix $A$.

We call the problem of interest Blind Regression for the following reason. In the setting of Regression, one observes data containing features and associated labels; the goal is to learn the functional relationship (or model) between features and labels assuming that labels are noisy observations. In our setting, tuples $(x_{\text{row}}(u),x_{\text{col}}(i))$ are the relevant (but unobserved) features and $Z(u,i)$ are noisy observations of associated labels $A(u,i)=f(x_{\text{row}}(u),x_{\text{col}}(i))$. We want to predict the value of $f(x_{\text{row}}(u),x_{\text{col}}(i))$ for all pairs $(u,i)$ for $u\in[m],i\in[n]$. In the sense that we want to predict the function $f$ evaluated on new points given previous data, the task has the feel of Regression. However, the features $(x_{\text{row}}(u),x_{\text{col}}(i))$ are latent; and thus we use the term Blind Regression.

Given an estimator $\hat{A}$ for the unknown matrix $A\in\mathbb{R}^{m\times n}$ of interest, we use the mean-squared error (MSE) to evaluate the performance of the estimator, defined as

[TABLE]

The expectation here is taken over all sources of randomness in the data generation process: (i) realization of the latent variables; (ii) realization of the noise variables; and (iii) masking. An estimator $\hat{A}$ is called consistent if $\lim_{m,n\to\infty}\textrm{MSE}(\hat{A})=0$. For a consistent estimator, we also want to establish an upper bound on the rate of convergence rate for the mean-squared error. Now we pose a follow-up question.

Problem 2**.**

Can we achieve a consistent estimator $\hat{A}$ for $A$ in the setup of Problem 1? If so, how fast does $\textrm{MSE}(\hat{A})$ decay to [math] for given $p=p(m,n)$ as the problem size $m,n\to\infty$?

II-C Exchangeability and Latent Variable Model

The latent variable model is well motivated and arises as a canonical representation for row and column exchangeable data, cf. [48] and [49]. Suppose that our data matrix $Z$ is a particular realization of the first $m\times n$ entries of a random array ${\bf Z}=\left\{{\bf Z}(u,i)\right\}_{(u,i)\in\mathbb{N}\times\mathbb{N}}$, which satisfies

[TABLE]

for every pair of permutations555The permutations over $\mathbb{N}$ are defined in the usual manner where only finitely many indices are permuted. $\sigma,\tau$ of $\mathbb{N}$. We use $\overset{d}{=}$ to denote that the joint distribution over $\{Z(u,i)\}_{u,i}$ is equivalently distributed as the joint distribution over $\{{\bf Z}\left(\sigma(u),\tau(i)\right)\}_{u,i}$, i.e. the random variables on both sides have the same distribution. Random array ${\bf Z}$ satisfying (4) is called exchangeable666To be precise, separately row and column exchangeable.. For an interested reader, [50] and [51] present overviews of exchangeable arrays.

In practice, the use of exchangeable arrays as a model is appropriate for variety of reasons. For example, in the setting of a recommendation system with anonymized data, this property may be reasonable if the order of the users in the system does not intrinsically carry information about the type of user; or in other words, if a user in the system could equally likely have been located in any row of the dataset.

In addition to exchangeability being quite a reasonable property for a wide variety of applications, it also leads to a convenient latent variable representation. The Aldous-Hoover representation theorem provides a succinct characterization for such exchangeable arrays. According to the theorem (see Corollary 3.3 in [51] for example), a random data array ${\bf Z}$ is exchangeable if and only if it can also be represented as

[TABLE]

where $\theta$, \big{\{}\theta_{\textrm{row}}(u)\big{\}}_{u\in\mathbb{N}}, \big{\{}\theta_{\textrm{col}}(i)\big{\}}_{i\in\mathbb{N}}, \big{\{}\theta_{\text{entry}}(u,i)\big{\}}_{(u,i)\in\mathbb{N}\times\mathbb{N}} are independent random variables drawn uniformly at random from the unit interval $[0,1]$, and $f_{\theta}$ is a measurable function indexed by the realization of $\theta$. As described in [51], this suggests the following generative model:

Sample an instance of $\theta\sim U[0,1]$ determining the governing function $f_{\theta}$. 2. 2.

For every row $u\in[m]$ and every column $i\in[n]$, independently sample uniform random variables $\theta_{\textrm{row}}(u)\sim U[0,1]$, $\theta_{\textrm{col}}(i)\sim U[0,1]$, $\theta_{\text{entry}}(u,i)\sim U[0,1]$. 3. 3.

Compute the realized data matrix $Z$ according to

[TABLE]

By comparing the model from (5) with the latent variable model described in Section II-A1 (cf. (1), (2)), we can see that the latent variable model considered in this work is a restricted subclass of exchangeable models that additionally impose an additive noise model and regularity conditions on the function $f_{\theta}$ in exchange for a more general latent space and associated probability measure.

Equivalence of Models

In our model we have conditioned on the universal index777Equivalently, we may think that $f_{\theta}=f$ for all $\theta\in[0,1]$. Note that we do not know what $f$ is a priori. $\theta$, such that given partial observations from matrix $Z$ for a particular $f_{\theta}$, our goal is to learn predicted outcomes of the realized $f_{\theta}$. Our model takes the form of

[TABLE]

where $\{x_{\text{row}}(u)\}_{u\in[m]},\{x_{\text{col}}(i)\}_{i\in[n]},\{\eta(u,i)\}_{(u,i)\in[m]\times[n]}$ are sampled independently. We can transform this to the form of (5) by considering $f$ to be equal to the realized function $f_{\theta}$, considering the latent variables $x_{\text{row}}(u)\sim P_{\mathcal{X}_{1}}$ and $x_{\text{col}}(i)\sim P_{\mathcal{X}_{2}}$ to be higher dimensional representations of $\theta_{\textrm{row}}(u)$ and $\theta_{\textrm{col}}(i)$ in spaces $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$, and considering the noise term $\eta_{\text{entry}}(u,i)$ to be generated by applying some transformation to the variable $\theta(u,i)$. Given these transformations, it becomes equivalent that

[TABLE]

Instead of allowing $f$ to be any arbitrary measurable function over $[0,1]\times[0,1]$, our model additionally imposes regularity conditions on $f$ by requiring it to be bounded and either finite types or Lipschitz continuous with respect to a higher dimensional representation $\mathcal{X}_{1}\times\mathcal{X}_{2}$. From a modeling perspective, we are effectively transferring the model complexity from a potentially complex measurable latent function over $[0,1]\times[0,1]$ to a simpler (e.g. Lipschitz) latent function over a potentially more complex latent variable space $\mathcal{X}_{1}\times\mathcal{X}_{2}$. The simple functional form provides analytic tractability for establishing theoretical results.

II-D Comparison with Other Models

In this section we discuss how our model relates to other models considered in the literature. In particular we would like to clarify which models are captured in the latent variable model and which are not. The key assumptions of our model are (1) bounded or Lipschitz latent function, (2) additive i.i.d. sub-gaussian noise, and (3) latent space being finite or more generally bounded Polish space.

II-D1 Low-rank Assumption

Suppose that the latent spaces $\mathcal{X}_{\textrm{row}}$ and $\mathcal{X}_{\textrm{col}}$ are equal dimensional vector spaces and the latent function $f$ is a bilinear form. This is a low-rank model with its rank being equal to $\dim\mathcal{X}_{\textrm{row}}=\dim\mathcal{X}_{\textrm{col}}$. If we assume $\mathcal{X}_{\textrm{row}}$ and $\mathcal{X}_{\textrm{col}}$ have bounded diameters and are equipped with some metric topology, then this model is a specific case of our model with Lipschitz constant being \big{(}\textrm{diam}~{}\mathcal{X}_{\textrm{row}}+\textrm{diam}~{}\mathcal{X}_{\textrm{col}}\big{)} (up to a multiplicative constant). Our latent variable model additionally assumes that the latent feature vectors are randomly sampled i.i.d. according to some probability measure, whereas typical low rank results allow for arbitrarily chosen latent feature vectors that satisfy incoherence.

Note that we utilize a different measure of ‘model complexity’ than the rank of the parameter matrix. When the underlying function is truly bilinear and \min\big{\{}\dim\mathcal{X}_{\textrm{row}},\dim\mathcal{X}_{\textrm{row}}\big{\}}\ll\textrm{diam}~{}\mathcal{X}_{\textrm{row}}+\textrm{diam}~{}\mathcal{X}_{\textrm{row}}, then rank of the data matrix could be a better measure to capture the essential complexity of the model; note that the rank of the data matrix does not scale as the diameter of the latent space increases. However, if the underlying model is nonlinear, the data matrix is likely to have full rank, e.g., when $f(x_{\text{row}}(u),x_{\text{col}}(i))=g(\left<x_{\text{row}}(u),x_{\text{col}}(i)\right>)$ for some monotone increasing nonlinear function $g$. In particular, when the size of the matrix far exceeds the ‘intrinsic model complexity’ (and hence so does the rank of the data matrix), it would make more sense to use our model.

Suppose that there is an underlying low rank matrix, yet the observation is an entrywise nonlinear monotone transformation of the low rank matrix; this is also known as single index models. [17] pointed out that adding nonlinearity could easily result in the underlying matrix no longer being low rank, requiring more complex approaches to estimate. In contrast, our model very easily handles nonlinearities as long as they satisfy local smoothness conditions.

II-D2 Biclustering, Submatrix Detection, Stochastic Block Model

The models used in biclustering, submatrix detection, and planted clique assume that there is some submatrix for which the data is significantly shifted in expectation compared to the rest of the matrix which is assumed to be uniform plus noise [52, 53]. This would correspond to “finite types” in our model, which can be modeled as a piecewise constant function. The desired goal is to detect or identify the deviant submatrix corresponding to a subset of the rows and columns. In contrast our task focuses on estimating the expected matrix. Given an estimation algorithm that could guarantee max error bounds entrywise or for each row/column, we could simply threshold to obtain an estimate for the deviant submatrix. Another distinction is that the literature in submatrix detection focuses on understanding thresholds of the minimal size submatrix that is detectable. However, in our setting as the latent feature vectors are sampled from a fixed distribution, the size of the deviant submatrix that our model studies would always be proportional to $O(mn)$ as long as the probability of sampling a deviant row or column is bounded below by a fixed constant. The stochastic block model also assumes finite types modeled through a piecewise constant function, which would fit within our latent variable model. It is commonly used to study the task of clustering, where the goal is to recover a true underlying partitioning of the rows and columns rather than only estimate the expected matrix [54]. Their setting specifically assumes a binary observation model, while our model assumes identically distributed additive sub-gaussian noise.

II-D3 Strong Stochastic Transitivity, Statistical Seriation, Low Permutation-Rank Matrices

There have been many models introduced to study matrices that arise from pairwise comparisons for the purposes of ranking or estimation. Many parametric models such as Bradley-Terry-Luce or Thurstone, can be described by a smooth function of scalar latent row and column variables. Our model adds the additional assumption that the latent variables must be sampled i.i.d. from a bounded diameter subspace. There have also been a series of nonparametric models introduced, in particular the strong stochastic transitivity assumption, which only requires that there exists some permutation of the rows and columns such that the entries are monotonically increasing along the permuted rows and columns [55, 56]. Statistical seriation assumes that there is some permutation of the rows and columns such that after the permutation, all rows and columns have the same shape, which could include monotonicity or other shape constraints [57]. Low permutation-rank matrices assume that a mixture model over strong stochastic transitive matrices, i.e. a mixture over permutation matrices [58].

Our latent variable model can encompass a subset of such matrices by imposing monotonicity assumptions on our latent function. However our model additionally assumes latent space being Polish space endowed with Borel measure, which are not required for permutation matrices. Similarly, a subset of mixtures of permutation matrices can be modeled by latent functions which are a mixture of different monotonic functions; however this may not be able to fully encompass all low permutation-rank models. If one were to impose boundedness on the entries, then it is possible that as the size of the matrix grows, strong stochastic transitivity might imply the existence of good nearest neighbor rows, which would then allow our algorithm to also perform well, but that remains an important direction to explore.

III Algorithm

Our algorithm builds on intuition from local approximation methods such as kernel regression. Therefore it takes the form of a similarity-based method, which first defines a kernel, i.e. similarity between pairs of rows or columns, and then computes the estimate for each matrix entry by averging over datapoints that are determined to be close to the entry of interest according to the ‘similarity.’

We present the general form of our algorithm in Section III-A. In Section III-B, we describe the basic user-user fixed radius nearest neighbor algorithm888This is equivalent to a variant of the classical similarity based collaborative filtering methods. as an example of the algorithm with a concrete choice of similarity and the averaging scheme to define the estimate. In Section III-C, we describe other variations of our algorithm, e.g., a variation of the algorithm that combines both row and column similarities to compute the kernel between datapoints.

The user-user nearest neighbor algorithm serves as our prime example and we provide theoretical guarantees for a vanishing upper bound on its mean squared error later in Section IV. Also, we show experimental results that suggest combining row and column similarities improve the quality of estimates; see Section VII.

III-A General Form

Our algorithm takes (i) the data matrix $Z$, (ii) a rule for sample splitting, and (iii) a thresholding parameter $\eta\geq 0$ as its inputs. The algorithm outputs $\widehat{A}\in\mathbb{R}^{m\times n}$.

Specifically, our algorithm determines the set of reliable neighbors by considering the ‘behavioral’ similarity in the function values. We need to define a similarity/dissimilarity statistic that can capture such similarities as well as can be computed from the data matrix. In Section III-B, we describe a version of our algorithm that utilizes the squared $\ell_{2}$ distance as a measure of dissimilarity, cf. Algorithm 2. This version of algorithm is the prime example we consider in this work and we provide a theoretical guarantee on its performance in Section IV.

We remark that some similarity functions can be preferred over others, depending on the model assumptions. For example, the squared $\ell_{2}$ distance is a natural choice when the latent space is assumed to be Euclidean, while the cosine similarity (= the angle between the latent variables) would be a more faithful measure of similarity when the latent space is the projective space (= spherical). We briefly discuss this matter in Section III-C with additional examples of variations of our algorithm.

Remark 1*.*

There are many variations in (1) how we determine the set of reliable neighbors (e.g., by defining similarity/dissimilarity functions), and (2) how we compute the final estimate based on the observations at those neighbors.

Remark 2*.*

We note that sample splitting is done for the ease of analysis, and is not essential in executing the algorithm.

III-B Prime Example: User-user Fixed Radius Nearest Neighbor (of the [math]-th order)

In this section, we describe our prime example of the generic algorithm described in Algorithm 1, namely, the user-user fixed radius neighbor algorithm. In this version of algorithm,

• •

we measure the dissimilarity between two rows with the squared $\ell_{2}$ distance between the rows at the overlapping column indices; specifically, we define the dissimilarity between two rows $u,v\in[m]$ as

[TABLE]

with the convention $0/0=\infty$, and

• •

for each $(u,i)\in[m]\times[n]$, we estimate $\widehat{A}(u,i)$ by averaging $Z(v,i)$ for $v\in[m]\setminus\{u\}$ such that $M(v,1)=1$ and $v$ is similar to the row $u$. We define that $v$ is similar to $u$ if the squared $\ell_{2}$ distance between their rows is no greater than some threshold $\eta\geq 0$, which is a tunable parameter input to the algorithm.

See Algorithm 2 for the full description. We state the algorithm for estimating a single entry $(u,i)$, which affects the sample splitting rule. The sample splitting is used for a cleaner analysis, but we believe that the results should extend without sample splitting as well.

Remark 3*.*

Note that Algorithm 2 is equivalent to the classical user-user fixed radius nearest neighbor collaborative filtering algorithm. The algorithm analyzed in the preliminary version of this work, cf. [59], is similar but not identical to Algorithm 2. The algorithm in [59] is motivated by the kernel regression of the first order and it is asymptotically equivalent to a mean-adjusted variant of the user-user $k$-nearest neighbor collaborative filtering algorithm. On the other hand, Algorithm 2 does not adjust the estimate the empirical means and it is kernel regression of the zeroth order in that sense.

Intuition behind the Algorithm

The following equations provide intuition for the reason why Algorithm 2 works. Assuming the row latent variables are instantiated as $x_{\text{row}}(1),\ldots,x_{\text{row}}(m)$ and $|\mathcal{B}_{\text{est}}(u,i)|=N_{\text{base}}(\eta)$, the mean squared error of the resulting estimate is given as999In fact, $\mathcal{B}_{\text{est}}(u,i)$ is not deterministic but random. We provide a complete analysis with formal proofs in Section IV and Appendix A.

[TABLE]

Here, \mathbb{E}_{\textrm{cond}}\left[~{}\cdot~{}\right]=\mathbb{E}\big{[}~{}\cdot~{}\big{|}~{}x_{\text{row}}(1),\ldots,x_{\text{row}}(m),|\mathcal{B}_{\text{est}}(u,i)|$$=N_{\text{base}}(\eta)\big{]} denotes the conditional expectation. This expression shows that the expected squared error conditioned on row latent variables is directly related to the squared $L^{2}$ distance between the slices of the latent function $f$ associated to rows $u$ and $v$. The good news is that the $L^{2}$ distance can in fact be estimated from the data itself. For any pair of row indices $u,v\in[m]$,

[TABLE]

Because of the concentration of measure, we expect the dissimilarity between rows $u$ and $v$ defined in (6) to be close to \big{\|}f(x_{\text{row}}(u),\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}+2\sigma^{2} as the size of the overlap $|\mathcal{B}_{\text{adj}}(u)\cap\mathcal{B}_{\text{adj}}(v)|$ increases.

Remark 4*.*

For fixed $m$ rows and $n$ columns, choosing a large $\eta$ leads to the increase in the size of $\mathcal{B}_{\text{est}}(u,i)$, which results in the increase in the first term in (7) and the decrease in the second term in (7). This demonstrates a bias-variance tradeoff associated to the choice of algorithmic parameter $\eta\geq 0$.

Remark 5*.*

Our algorithm does not require any prior knowledge on the model parameters such as the regularity parameters $L$, diameter of ${\mathcal{X}_{\textrm{row}}}$ and the noise variance $\sigma$. The algorithmic parameter $\eta\geq 0$ can be chosen arbitrarily. However, we need $\eta$ to be in a certain range to achieve good theoretical guarantee for the upper bound on its MSE. See Theorem 1.

III-C Additional Examples of Our Algorithm: Other Variations

In this section, we exhibit other variations of our algorithm and provide intuition for them. However, we do not formally prove error bounds for these variations in this paper.

III-C1 Item-item and User-item Variants

In Section III-B, we described a user-user fixed radius nearest neighbor algorithm that utilizes the (dis-)similarity between rows, $\textrm{dissim}_{\text{row}}(u,v)$ for $u,v\in[m]$. We can apply the same idea to use the similarity between columns, $\textrm{dissim}_{\text{col}}(i,j)$ for $i,j\in[n]$. Moreover, once the similarity between rows and columns are identified, we can define similarity between any pair of index tuples $(u,i),(v,j)\in[m]\times[n]$ based on the $\textrm{dissim}_{\text{row}}(u,v)$ and $\textrm{dissim}_{\text{col}}(i,j)$. For example, we may define $\textrm{dissim}((u,i),(v,j))=\max\{\textrm{dissim}_{\text{row}}(u,v),\textrm{dissim}_{\text{col}}(i,j)\}$.

III-C2 Kernels for Weighting

Once the (dis-)similarities between rows and columns are identified, we need a weighting scheme to define the estimates. For example, in Algorithm 2, we use the the hard-thresholding kernel in which the weight is $1$ if the row dissimilarity is less than $\eta$ and the weight is [math] otherwise. However, this choice is not necessary and our algorithm can be implemented with any choice of kernels (i.e. weighting scheme).

Example 1** (Gaussian Kernel Weights).**

Given $(u,i)\in[m]\times[n]$, suppose that we have computed dissimilarity between $(u,i)$ and $(v,j)$ for all $(v,j)\in\mathcal{B}_{\text{est}}(u,i)$, e.g., by defining $\textrm{dissim}((u,i),(v,j))=\max\{\textrm{dissim}_{\text{row}}(u,v),\textrm{dissim}_{\text{col}}(i,j)\}$. Then, given a bandwith parameter $\lambda\in\mathbb{R}_{+}$, we define the weights according to a Gaussian kernel101010One may doubt why we call this Gaussian kernel instead of Laplacian kernel. The reason is that we treat \textrm{dissim}((u,i),(v,j))\big{)} as a proxy of the squared $L^{2}$ distance. as: for each $(v,j)\in[m]\times[n]$,

[TABLE]

and then define

[TABLE]

When $\lambda\to\infty$, the estimate $\widehat{A}(u,i)$ only depends on the ‘nearest’ neighbor $(v,j)$. On the other hand, when $\lambda=0$, the algorithm equally averages all the $Z(v,j)$ for $(v,j)\in\mathcal{B}_{\text{est}}(u,i)$. Empirically, this variant of the algorithm seems to perform very well with an appropriate selection of the bandwidth parameter $\lambda$, which can be tuned using cross validation.

III-C3 Other Ways to Define Dissimilarity

We suggest other alternatives to define similarities.

Higher Order Information

Algorithm 2 detects reliable neighbors $(v,j)\in\mathcal{B}_{\text{est}}(u,i)$ and estimates $\widehat{A}(u,i)$ by averaging $Z(v,j)$. This procedure is equivalent to traditional local smoothing, or kernel regression of the [math]-th order, except that we needed to estimate the ‘proximity’ of the latent features from the data matrix. If we consider $Z(v,j)$ as a constant function, then $\widehat{A}(u,i)$ can be viewed as the evaluation of the average of the constant function $Z(v,j)$ at $(u,i)$.

The next question naturally arises: “Can we come up with an algorithm that is equivalent to higher-order kernel regression algorithm?” For example, one can build linear regression estimators centered at $(v,j)$ for each $(v,j)\in\mathcal{B}_{\text{est}}(u,i)$, instead of considering the constant function at the level of $Z(v,j)$. Evaluating the average of such linear estimators at $(u,i)$ will yield a different, probably more refined, estimate $\widehat{A}(u,i)$. This turns out to be equivalent to kernel regression of the $1$-st order with latent features. This version of algorithm is analyzed in the preliminary version of this paper, cf. [59].

Beyond the Euclidean Latent Space

In our proposed algorithm, we defined the dissimilarities between the rows and the columns using empirical $L^{2}$ distance of the slices of the latent function. Depending on the geometry of the latent space, other ways of defining dissimilarities can be favorable. For example, suppose that only the ‘direction’ of the rows in the data matrix matters, while the amplitude does not carry much information. Then it would make sense to measure the dissimilarity between rows by estimating their ‘angles’ instead of their Euclidean distance111111To be precise, we want to define the dissimilarity between the equivalent classes of the rows. In this example, the equivalence relation is defined by positive scaling. In fact, this is what cosine similarity measures, which is commonly used in classical collaborative filtering.

More importantly, the canonical latent space in the data generation process can be non-Euclidean. For example, the latent variables could be drawn from a sphere121212more generally, from a Riemannian manifold with non-zero curvature. Cosine similarity makes use of this knowledge, while the $L^{2}$ dissimilarity does not – it essentially isometrically embeds the sphere to a higher-dimensional Euclidean space and then uses the metric of the ambient space. This could be problematic because one may need a much larger dimension for Euclidean embedding than the intrinsic dimension. Additionally in some cases it may be impossible to find an isometric embedding to the Euclidean space of any finite dimension; in particular, when the latent space is negatively curved. One prominent example of such latent space is the tree equipped with the geodesic distance; it is impossible to find an isometric embedding of an infinitely growing $3$-regular tree to $\mathbb{R}^{d}$ for any $d<\infty$.

Both the question of finding a good similarity/dissimilarity function that utilizes the structure of the latent space and the question of systematically exploiting higher-order information of the latent function remain important direction of future research.

III-D Computational Complexity

In order to compute $\textrm{dissim}_{\text{row}}(u,v)$, we need to average differences over the set $|\mathcal{B}_{\text{adj}}(u)\cap\mathcal{B}_{\text{adj}}(v)|$. As the sparsity of each row is in expectation $np$, the cost of the similarity computation is in expectation bounded by $O(np)$. We need to compute $\textrm{dissim}_{\text{row}}(u,v)$ for all $\binom{m}{2}$ pairs. For each row $u$, we need to compute the set $\{v:\textrm{dissim}_{\text{row}}(u,v)\leq\eta\}$, which takes $O(m)$. To compute the final estimate for each entry $(u,i)$ we need to average over the set $\mathcal{B}_{\text{est}}(u,i)$, which in expectation is bounded above by $mp$. Therefore, the computational complexity is bounded by $O(m^{2}(np)+m^{2}+mn(mp))=O(m^{2}np)$.

The algorithm can be parallelized into two phases. In phase 1, we can in parallel compute $\textrm{dissim}_{\text{row}}(u,v)$ for all $\binom{m}{2}$ pairs, each taking computation time $O(np)$. In phase 2, we can compute the estimates $\hat{A}(u,i)$ in parallel for all $mn$ indices of the matrix, each taking computation time $O(m+mp)$. Therefore assuming we have $O(m^{2}+mn)$ processors that can compute simultaneously, the parallelized computational complexity is $O(np+m)$.

IV Main Results

IV-A Theorem Statement

We present a theorem which upper bounds the MSE of the estimate produced by the user-user fixed radius nearest neighbor algorithm presented in Section III-B (cf. Algorithm 2). Recall from our problem statement in Section II that $p$ is the probability each entry is observed, $D_{f}$ is an upper bound on magnitude of the latent function, and $\sigma^{2}$ is an upper bound on the noise variance131313In fact, $\|N(u,i)\|_{\psi_{2}}\leq\sigma$. We may identify $\sigma^{2}$ with an upper bound on the variance, up to an absolute constant..

IV-A1 Informally Stated Asymptotic Upper Bounds

We present an informal version of our main theorem under a simplified setting for ease of exposition. The general form of main result is stated in Theorem 1 which is rather involved. To that end, consider a simplified setting where there are only finite types of row latent variables. In that case, we obtain the following result.

Theorem** (Informal Statement of Corollary 1).**

Suppose that the measure for the latent row space has finite support. If p\gg\max\Big{(}\frac{1}{m},\sqrt{\frac{\log m}{n}}\Big{)}, then

[TABLE]

Note that the prefix constant $D_{f}+\sigma$ is the measure of variability (complexity) for the latent function (up to noise).

There are two main sources of the estimation error: (i) the error in estimation of similarities (step 2 of Algorithm 2) and (ii) the error in approximation by smoothing (step 3 of Algorithm 2). Intuitively, empirical estimates of similarity/dissimilarity between two rows becomes more accurate as $(n-1)p^{2}$, which is equal to the expected size of the overlap, increases. This is captured by the term $\sqrt{\frac{\log(m-1)}{(n-1)p^{2}}}$.

Given the similarities are sufficiently accurate, the estimation error is dominated by the approximation error. There are $|\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})|$ types of rows and $(m-1)p$ number of available samples, hence, $|\mathcal{B}_{\text{est}}(u,i)|\approx\frac{(m-1)p}{|\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})|}$ in expectation. This captures the second term in the error. This also reflects the impact of local geometry of the measure $\mu_{\mathcal{X}_{\textrm{row}}}$ on the estimation error.

IV-A2 Formal Statement of the Main Theorem

We define the function $\phi$ for $x\in\mathcal{X}_{\textrm{row}}$ and $r\geq 0$ as

[TABLE]

Theorem 1 presents an upper bound on the mean-squared error (cf. (3)) of the user-user Fixed radius nearest neighbor algorithm described in Algorithm 2. Before presenting the theorem statement, we remark that

[TABLE]

due to the exchangeability of the model and the linearity of the expectation. In Theorem 1, we provide an upper bound on \mathbb{E}\left[\big{(}\widehat{A}(1,1)-A(1,1)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(1)\right], conditioned on $x_{\text{row}}(1)$, which is the latent feature of the first row.

Theorem 1** (Main Theorem).**

Let $\widehat{A}$ be the estimator returned by the Algorithm 2. Let constant K\triangleq\Big{(}\frac{2D_{f}}{\sqrt{\ln 2}}+2\sigma\Big{)}^{2}. Suppose that our algorithm uses threshold parameter \eta\geq\eta^{\prime}+K\max\bigg{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\bigg{)} for some $\eta^{\prime}\geq 2\sigma^{2}$. Then

[TABLE]

In the above expression, $C,c>0$ are absolute constants, and \phi(x,r)=\mathbb{P}_{x_{\text{row}}(v)\sim\mu_{\mathcal{X}_{\textrm{row}}}}\left(\big{\|}f(x,\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}\leq r\right) for any $x\in\mathcal{X}_{\textrm{row}}$ and $r>0$.

The threshold $\eta\geq 0$ is a tunable parameter that our algorithm uses. Our analysis captures the bias-variance tradeoff associated with the parameter $\eta$. When $\eta$ is too large, our upper bound becomes loose because the algorithm utilizes information from less reliable neighbors; on the other hand, when $\eta$ is chosen too small, the resulting estimate count only on a few neighbors and suffers from large variance.

The first two terms of the MSE bound come from bounding the bias of the estimator. The first term \big{(}\eta-2\sigma^{2}\big{)} reflects the bias that is unavoidable due to the selection of the threshold $\eta$ and the error in estimating $\textrm{dissim}_{\text{row}}(u,v)$, and the second term bounds the tail of the bias along with the bad event that $|\mathcal{B}_{\text{adj}}(u)\cap\mathcal{B}_{\text{adj}}(v)|$ is too small to produce a good estimate of $\textrm{dissim}_{\text{row}}(u,v)$. The remaining three terms bound the variance of the estimator; the error bound for the trivial estimates (when $\mathcal{B}_{\text{est}}(1,1)=\emptyset$), given in the form of $D_{f}^{2}\mathbb{P}\left(\mathcal{B}_{\text{est}}(1,1)=\emptyset\right)$ is also subsumed in these three terms.

Taking expectation of the bound in Theorem 1 with respect to $x_{\text{row}}(1)\sim\mu_{\mathcal{X}_{\textrm{row}}}$, we obtain the desired bound on MSE, i.e.

[TABLE]

IV-B Implications

Theorem 1 provides an upper bound on MSE for Algorithm 2, which depends on the underlying probability measure $\mu_{\mathcal{X}_{\textrm{row}}}$ and the latent function $f$ through $\phi$. Here we evaluate this implicit bound for three special examples:

• •

The latent space has finitely many elements, or equivalently $\mu_{\mathcal{X}_{\textrm{row}}}$ has finite support.

• •

The latent space is unit hypercube in a finite dimensional space, latent function is Lipschitz.

• •

The latent space is a complete, separate metric space aka Polish space with bounded diameter, latent function is Lipschitz.

IV-B1 Finite Types

We state the following Corollary of Theorem 1 when $\mu_{\mathcal{X}_{\textrm{row}}}$ has finite support.

Corollary 1** (Finite support).**

Suppose that $\mathcal{X}_{\textrm{row}}$ is equipped with the discrete metric141414$d_{\mathcal{X}_{\textrm{row}}}(x_{1},x_{2})=1$ if and only if $x_{1}\neq x_{2}$. topology and $\mu_{\mathcal{X}_{\textrm{row}}}$ has finite support in $\mathcal{X}_{\textrm{row}}$ with $\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})$ denoting the support of $\mu_{\mathcal{X}_{\textrm{row}}}$. Let K=\Big{(}\frac{2D_{f}}{\sqrt{\ln 2}}+2\sigma\Big{)}^{2}, $\eta^{\prime}=2\sigma^{2}$, and

[TABLE]

If p\geq\max\bigg{(}\frac{8}{m-1},~{}\Big{(}4\vee\sqrt{\frac{2}{c}}\Big{)}\cdot\sqrt{\frac{\log(m-1)}{n-1}}\bigg{)}, then

[TABLE]

where $C,c>0$ are absolute constants.

That is, as long as p\gg\max\Big{(}\frac{1}{m},\sqrt{\frac{\log m}{n}}\Big{)}, under discrete measure with finite support

[TABLE]

Proof of Corollary 1.

Our interest is in bounding

[TABLE]

and

[TABLE]

from the statement of Theorem 1 to obtain the desired result. Since we are considering discrete metric with measure having finite support, for any $r\geq 0$ and $x\in\textrm{supp}(\mu_{\mathcal{X}_{\textrm{row}}})$, $\phi(x,r)\geq\mu_{\mathcal{X}_{\textrm{row}}}(x)$. Given choice of $\eta^{\prime}$, $\eta^{\prime}-2\sigma^{2}\geq 0$. Therefore, (9) can be written as

[TABLE]

Similarly, (10) reduces to

[TABLE]

In above, we have used an easy to verify fact that \sup_{\theta\in[0,1]}\theta\exp\bigg{(}-\frac{(m-1)p}{8}\theta\bigg{)} is achieved for $\theta=\frac{8}{(m-1)p}\in[0,1]$. Replacing (IV-B1) and (IV-B1) in the statement of Theorem 1, and realizing that other terms in the bound of Theorem 1 are asymptotically smaller order, we obtain (1). ∎

IV-B2 Uniform Measure over $[0,1]^{d}$, Lipschitz Latent Function

Let $L$ denote the Lipschitz constant of the latent function $f$. We consider the setting where $\mu_{\mathcal{X}_{\textrm{row}}}$ is uniform Lebesgue measure on $[0,1]^{d}$, the unit cube in $d$ dimension. We define $B(x,r)=\{x^{\prime}\in\mathcal{X}_{\textrm{row}}:d_{\mathcal{X}_{\textrm{row}}}(x,x^{\prime})\leq r\}$ for any $r>0$. By Lipschitzness, $d_{\mathcal{X}_{\textrm{row}}}(x_{\text{row}}(1),x_{\text{row}}(v))\leq\frac{\sqrt{\eta^{\prime}-2\sigma^{2}}}{L}$ implies \big{\|}f(x_{\text{row}}(1),\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}\leq\eta^{\prime}-2\sigma^{2}. Therefore,

[TABLE]

There exists universal constants $\alpha,\beta>0$ such that for any $d\geq 1$, and $x\in[0,1]^{d}$ and $r>0$,

[TABLE]

We shall assume that $\mathcal{X}_{\textrm{col}}=[0,1]^{d}$ as well. Therefore,

[TABLE]

and hence $D_{f}\leq L\sqrt{d}$ since there exists $\alpha\in\mathcal{X}_{\textrm{row}},\beta\in\mathcal{X}_{\textrm{col}}$ such that $f(\alpha,\beta)=0$. Now, we state the following implication of Theorem 1 in this setting.

Corollary 2** (Unit cube and the Lebesgue measure).**

Let $\mu_{\mathcal{X}_{\textrm{row}}}$ be the uniform measure on $\mathcal{X}_{\textrm{row}}=[0,1]^{d}$. Let K=\big{(}\frac{2L\sqrt{d}}{\sqrt{\ln 2}}+2\sigma\big{)}^{2}, \eta^{\prime}=2\sigma^{2}+\alpha^{2/d}\beta^{2}L^{2}\big{(}(m-1)p\big{)}^{-\frac{2}{d+2}}, and

[TABLE]

If p\geq\max\bigg{(}\frac{8}{m-1},~{}\Big{(}4\vee\sqrt{\frac{2}{c}}\Big{)}\cdot\sqrt{\frac{\log(m-1)}{n-1}}\bigg{)}, then

[TABLE]

where $C,C^{\prime},c>0$ are absolute constants.

That is, as long as p\gg\max\Big{(}\frac{1}{m},\sqrt{\frac{\log m}{n}}\Big{)}, under the uniform measure,

[TABLE]

Proof of Corollary 2.

Similar to proof of Corollary 1, our interest is in bounding (9) and (10). To that end, under uniform Lebesgue measure $\mu_{\mathcal{X}_{\textrm{row}}}$, on $d$ dimensional unit cube $\mathcal{X}_{\textrm{row}}$, it follows that for any $x\in\mathcal{X}_{\textrm{row}}$ and $r>0$

[TABLE]

Therefore, by choice of $\eta^{\prime}$

[TABLE]

Similarly,

[TABLE]

From (IV-B2) and (IV-B2), and statement of Theorem 1, and realizing that other terms in the bound of Theorem 1 are asymptotically smaller order, we conclude (2). ∎

IV-B3 Bounded Polish Space, Lipschitz Latent Function

Let $L$ denote the Lipschitz constant of the latent function $f$. We define $B(x,r)=\{x^{\prime}\in\mathcal{X}_{\textrm{row}}:d_{\mathcal{X}_{\textrm{row}}}(x,x^{\prime})\leq r\}$ for any $r>0$. By Lipschitzness, $d_{\mathcal{X}_{\textrm{row}}}(x_{\text{row}}(1),x_{\text{row}}(v))\leq\frac{\sqrt{\eta^{\prime}-2\sigma^{2}}}{L}$ implies \big{\|}f(x_{\text{row}}(1),\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}\leq\eta^{\prime}-2\sigma^{2}. Therefore,

[TABLE]

We consider the row latent space $\mathcal{X}_{\textrm{row}}$ to be a complete, separable metric space, i.e. a Polish space151515Recall that a metric space is separable if it has a countable dense subset, and complete if every Cauchy sequence converges within the space.. Let the diameter of the space be bounded, i.e. there exists finite $D>0$ such $\sup_{x,x^{\prime}\in\mathcal{X}_{\textrm{row}}}d_{\mathcal{X}_{\textrm{row}}}(x,x^{\prime})\leq D$. We shall also assume that the diameter of $\mathcal{X}_{\textrm{col}}$ is bounded by $D$ as well. This implies

[TABLE]

since there exists $\alpha\in\mathcal{X}_{\textrm{row}},\beta\in\mathcal{X}_{\textrm{col}}$ such that $f(\alpha,\beta)=0$.

An important consequence of $\mathcal{X}_{\textrm{row}}$ being Polish is that $\mu_{\mathcal{X}_{\textrm{row}}}$ is tight, i.e. for any $\delta>0$, there exists a compact set $S_{\delta}\subseteq\mathcal{X}_{\textrm{row}}$ such that $\mu_{\mathcal{X}_{\textrm{row}}}(S_{\delta})\geq 1-\delta$. Due to compactness of $S_{\delta}$ and the space being Polish, there exists a finite number of balls of any given radius $\varepsilon>0$ of choice such that they cover $S_{\delta}$. That is, for any $\varepsilon,\delta>0$, the effective covering number $N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)$ is always finite. Let $B(x_{i},\varepsilon)$ for $i\in[N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)]$ denote the collection of balls so that $S_{\delta}\subseteq\cup_{i\in[N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)]}B(x_{i},\varepsilon)$. By construction, \mu_{\mathcal{X}_{\textrm{row}}}\Big{(}\cup_{i\in[N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)]}B(x_{i},\varepsilon)\Big{)}\geq\mu_{\mathcal{X}_{\textrm{row}}}(S_{\delta})\geq 1-\delta. Define $\text{BAD}=\{i\in[N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)]~{}:~{}\mu_{\mathcal{X}_{\textrm{row}}}(B(x_{i},\varepsilon))<\delta/N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)\}$. Let $\text{GOOD}=\{i\in[N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)]~{}:~{}i\notin\text{BAD}\}$. Therefore,

[TABLE]

since $|\text{BAD}|\leq N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)$. Therefore,

[TABLE]

For any $x\in\cup_{i\in\text{GOOD}}B(x_{i},\varepsilon)$, there exists $i\in\text{GOOD}$ such that $x\in B(x_{i},\varepsilon)$. Hence, $B(x_{i},\varepsilon)\subseteq B(x,2\varepsilon)$. Therefore, for any $x\in S^{\text{GOOD}}(\delta,\varepsilon)\equiv\cup_{i\in\text{GOOD}}B(x_{i},\varepsilon)$,

[TABLE]

where \mu_{\mathcal{X}_{\textrm{row}}}\Big{(}S^{\text{GOOD}}(\delta,\varepsilon)\Big{)}~{}\geq~{}1-2\delta. By (16) it follows that for $x_{\text{row}}(1)\in S^{\text{GOOD}}(\delta,\varepsilon)$ with \mu_{\mathcal{X}_{\textrm{row}}}\Big{(}S^{\text{GOOD}}(\delta,\varepsilon)\Big{)}\geq 1-2\delta,

[TABLE]

for any choice of $\varepsilon>0,~{}\delta\in(0,1/2)$. Let us choose $\varepsilon=\varepsilon(\delta)$ where $4\varepsilon^{2}(\delta)L^{2}=\delta$. Replacing this choice of $\varepsilon$ in (17), we obtain that for any $\delta>0$, there exists a set $S^{\prime}(\delta)\equiv S^{\text{GOOD}}(\delta,\varepsilon(\delta))$ such that (i) \mu_{\mathcal{X}_{\textrm{row}}}\big{(}S^{\prime}(\delta)\big{)}\geq 1-2\delta and (ii) for all $x_{\text{row}}(1)\in S^{\prime}(\delta)$,

[TABLE]

For all $t\geq N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{1}{2L},1)$, we define function $\delta^{\star}(t)$ to be

[TABLE]

For any $\delta>0$, $N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{\sqrt{\delta}}{2L},\delta)$ is finite, and thus for $t\geq\delta^{-2}N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{\sqrt{\delta}}{2L},\delta)$, $\delta^{\star}(t)\leq\delta$. We can verify that $\delta^{\star}(t)$ is also monotonically non-increasing. Therefore, it follows that

[TABLE]

Using the above developed machinery, now we are ready to bound $\textrm{MSE}(\hat{A})$ as summarized in Corollary 3.

Corollary 3** (Bounded Polish Space).**

Let $\mu_{\mathcal{X}_{\textrm{row}}}$ be any measure on $\mathcal{X}_{\textrm{row}}$, which is assumed to be a bounded diameter Polish space. Let the diameter of $\mathcal{X}_{\textrm{row}}$ and $\mathcal{X}_{\textrm{col}}$ be bounded above by $D$. Let K=\big{(}\frac{2L\sqrt{d}}{\sqrt{\ln 2}}+2\sigma\big{)}^{2}, $\eta^{\prime}=2\sigma^{2}+\delta^{\star}((m-1)p)$, and

[TABLE]

If p\geq\max\bigg{(}\frac{8N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{1}{2L},1)}{m-1},~{}\Big{(}4\vee\sqrt{\frac{2}{c}}\Big{)}\cdot\sqrt{\frac{\log(m-1)}{n-1}}\bigg{)}, then

[TABLE]

where $C,C^{\prime},c>0$ are absolute constants, and $\delta^{\star}(t)\to 0$ as $t\to\infty$.

As an immediate consequence of Corollary 3, it follows that as long as p\gg\max\Big{(}\frac{1}{m},\sqrt{\frac{\log m}{n}}\Big{)}, for any measure on a bounded Polish space,

[TABLE]

This implies that our estimator is consistent for any latent variable model with bounded Polish space as long as $p\geq\max(m^{-1+\delta},n^{-\frac{1}{2}+\delta})$ for any $\delta>0$ and $\log m=o(n^{\delta})$.

If $\mathcal{X}_{\textrm{row}}=[0,1]^{d}$, for any $\delta\in(0,1)$, $N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\varepsilon,\delta)=O(\varepsilon^{-d})$. Therefore, $\delta^{*}(t)=O(t^{-\frac{2}{d+4}}L^{\frac{2d}{d+4}})$. We can conclude that as long as p\gg\max\Big{(}\frac{1}{m},\sqrt{\frac{\log m}{n}}\Big{)}, and $\mathcal{X}_{\textrm{row}}=[0,1]^{d}$,

[TABLE]

with universal constant $C>0$. It is worth noting the similarity between (2) and (IV-B3) – the only difference is the $d+4$ instead of $d+2$ in the denominator of the exponent for term $(m-1)p$. This is precisely the minimal cost of generalizing from the specific uniform distribution to any arbitrary distribution.

Proof.

Since $p\geq\frac{8N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{1}{2L},1)}{m-1}$, as argued before, $\delta^{\star}((m-1)p)$ is well defined and we shall choose $\delta=\delta^{\star}((m-1)p)$. From the choice of $\eta^{\prime}$, for any $x_{\text{row}}(1)\in S^{\prime}(\delta)$, it follows from (18) that $\phi(x_{\text{row}}(1),\delta)\geq\frac{\delta}{N_{\text{eff}}(\mathcal{X}_{\textrm{row}},\frac{\sqrt{\delta}}{2L},\delta)}$. Therefore, for any $x_{\text{row}}(1)\in S^{\prime}(\delta)$,

[TABLE]

where the last inequality follows from the fact that $\delta=\delta^{\star}((m-1)p)$. Similarly,

[TABLE]

where again the last inequality follows from the fact that $\delta=\delta^{\star}((m-1)p)$. Replacing (IV-B3) and (IV-B3) in the statement of Theorem 1, and realizing that the other terms in the bound of Theorem 1 are asymptotically smaller order and \exp\big{(}-\frac{1}{8\delta}\big{)}\leq\delta for any $\delta\in(0,1)$, we obtain

[TABLE]

By Cauchy-Schwarz inequality and the fact that $\mu_{\mathcal{X}_{\textrm{row}}}(S^{\prime}(\delta^{\star}((m-1)p)))\geq 1-2\delta^{\star}((m-1)p))$, we have

[TABLE]

In what follows, we shall argue that

[TABLE]

By (IV-B3), (IV-B3) and (25), the main claim (3) follows. With that in mind, we shall now establish (25) to conclude the proof. Recall from Algorithm 2 that

[TABLE]

in our user-user fixed radius nearest neighbor algorithm when $|\mathcal{B}_{\text{est}}(1,1)|\geq 1$ and $\widehat{A}(1,1)=0$ when $|\mathcal{B}_{\text{est}}(1,1)|=0$. Therefore,

[TABLE]

By the assumption that the magnitude of the latent function $f$ is bounded by $D_{f}$,

[TABLE]

By introducing the notations

[TABLE]

we can rewrite the expected squared error as

[TABLE]

where we have used the independence of terms in $X$ and $Y$. It immediately follows that $X$ is a bounded random variable with $|X|\leq 2D_{f}$. Therefore,

[TABLE]

The randomness influencing the selection of $\mathcal{B}_{\text{est}}(1,1)$ is independent of the random variables in the summation of term $Y$. Conditioned on $|\mathcal{B}_{\text{est}}(1,1)|=\ell$ for any $\ell\geq 1$, $Y$ is simply a summation of $\ell$ i.i.d. random variables, each with ${\psi_{2}}$ norm bounded by $\sigma$ and zero mean. Therefore, it follows that conditioned on $|\mathcal{B}_{\text{est}}(1,1)|=\ell$ for any $\ell\geq 1$,

[TABLE]

From (IV-B3)-(28), the inequality (25) follows. This completes the proof of Corollary 3. ∎

V Discussion

V-A Intuition

V-A1 Structural Assumptions and Sample Complexity

Without any structure assumed on the latent function $f$ and the latent spaces $\mathcal{X}_{\textrm{row}},\mathcal{X}_{\textrm{col}}$, one would need $mn$ number of samples to recover the matrix. Our framework relies on two key assumptions to reduce the sample complexity: (i) (most of) the latent spaces can be covered by balls centered at a few representative points; and (ii) the latent function $f$ is regular (Lipschitz) and hence the proximity of two points in the latent space results in the similarity of the function values, which is observable from data.

Given a latent space $(\mathcal{X}_{\textrm{row}},d_{\mathcal{X}_{\textrm{row}}},\mu_{\mathcal{X}_{\textrm{row}}})$, we want to cover the entire space $\mathcal{X}_{\textrm{row}}$ minus a small fraction that has negligible measure $\varepsilon>0$, with a collection of balls of radius $\delta>0$. In other words, our model considers the metric entropy161616The smallest number of bits that suffices to specify every point $x$ in the set with accuracy $\delta$ in the metric $d_{\mathcal{X}_{\textrm{row}}}$. of $\mathcal{X}_{\textrm{row}}$ minus $\varepsilon$-fraction as a measure of complexity. However, we don’t have direct access to the latent features, but the distance between two points must be estimated from the pattern of the associated function values. If the latent function were an isometry, the distance in the image would faithfully reflect the distance in the domain. Lipschitzness assumption on $f$ is a robust analogue of the isometry assumption; if the latent function $f$ is $L$-Lipschitz, then the distance in the domain (= latent space) cannot be inflated more than $L$ times in the image; therefore, the image of the latent function for close neighbors in the latent space will behave in a similar fashion. With the Lipschitz assumption, the cost of indirect measurement is no more than $L$, and it suffices to consider a $\frac{\delta}{L}$-covering of the latent space instead of a $\delta$-covering.

V-A2 MSE Upper Bounds and the Rate of Convergence

The minimax optimal rate for nonparametric regression is $O(N^{-2/(2+d)})$ where $N$ is the number of observations uniformly distributed in a $d$-dimensional hypercube. In our algorithm we use the data across columns to learn similarities between the rows, yet when we actually compute the final estimates, we only use the datapoints in each column separately. The number of observations per column is $mp$ in expectation. In the case when row latent features are sampled from a uniform measure over a $d$-dimensional cube, the second term of our MSE bound from Corollary 2 is $O((mp)^{-2/(2+d)})$. This indicates that the second term of our MSE bound is optimal for any estimator that uses entries in each column separately. In particular, even if the algorithm were given oracle knowledge of the row latent features it could not do better as long as it constrains itself to estimating each column separately. An algorithm which would average both similar rows and columns could be able to improve beyond this bound. The first term of our MSE bound comes from the step of estimating $\textrm{dissim}_{\text{row}}(u,v)$ for all pairs of rows; however this first term is quickly dominated by the second term as $d$ increases, suggesting that the second term dominates the MSE. This suggests that our analysis is tight for our estimator.

Our result can be compared with the upper bound on the MSE for the UVST estimator as presented in Theorem 2.7 of [19]. For simplicity, consider the setting of a square matrix, i.e. $m=n$. For a matrix sampled from the latent variable model with latent variable dimension $d$, their theorem guarantees that

[TABLE]

for some constant $C$ as long as $p\geq m^{-1+\delta}$. This upper bound is meaningful only when $p>m^{-\frac{2}{d+2}}$, because the MSE bound in (29) is bounded below by $C$ when $p\leq m^{-\frac{2}{d+2}}$. However, requiring $p>m^{-\frac{2}{d+2}}$ can be too restrictive when the latent dimension $d$ is large since it means that we need to sample almost every entry to achieve a nontrivial bound. Chatterjee’s result stems from showing that a Lipschitz function can be approximated by a piecewise constant function, which upper bounds the rank of the (approximate) target matrix. This global discretization results in a large penalty with regards to the dimension of the latent space.

In recent work, [20] extends the analysis of USVT for graphon estimation, when the observation matrix is binary. He shows if the latent function is $\alpha$-Hölder smooth, the spectrum decays polynomially, and thus the MSE of the USVT estimator is bounded above by

[TABLE]

This refined bound shows that at least for the binary observation case, the USVT estimator is consistent as long as $p=\omega(n^{-1})$, independent from the latent dimension $d$, relying on the quick decay of the spectrum.

Our algorithm and analysis provides a vanishing upper bound on the MSE whenever $p\geq\max\left\{m^{-1+\delta},n^{-1/2+\delta}\right\}$, also independent of the latent dimension. In fact even as $d$ grows with $m$, as long as $d=o(\log m)$, our analysis guarantees that our algorithm achieves a vanishing MSE. Our analysis relies on the “local” structure of the latent space; even if the latent dimension increases, we only need to ensure that there are sufficiently many close neighbor rows so that nearest neighbor averaging produces a good estimate.

V-B Future Work

Our analysis seems tight due to the comparison with minimax rates for noparametric regression, however it is likely that a model which could average entries across both rows and columns would be able to improve the MSE bounds. In particular, the minimax optimal rates for nonparametric regression would imply a lower bound of

[TABLE]

in the setting where both row and column latent features are sampled uniformly from a $d$-dimensional hypercube. This is achievable with an oracle estimator that were given knowledge of the latent features, but we do not know of information theoretic lower bounds for the MSE that are specific to the latent variable model where features are unobserved. For specific settings such as when the function $f$ when considered as an integral operator has finite spectrum, it is equivalent to low-rank models, for which lower bounds have been characterized. For specific noise models such as the binary observation model which corresponds to the graphon generative model for random graphs, [60, 61] show that variants of the least squares estimator achieve optimal rates, but unfortunately they are not polynomial time computable.

From an implementation perspective, the similarity based algorithm proposed in this work, similar to classical collaborative filtering methods, is easy to implement and scales extremely well to large datasets, as it naturally enjoys a parallelizable implementation. Furthermore, the operation of finding $k$ nearest neighbors can benefit from computational advances in building scalable approximate nearest neighbor indices, cf. [62, 63].

Next we discuss some natural extensions and directions for future work. In our model, the latent function $f$ is assumed to be Lipschitz. However, the proof only truly utilizes the fact that “locally” the function value does not oscillate too wildly. Intuitively, this suggests that the result may extend to a broader class of functions, beyond Lipschitz functions. For example, a function with bounded Fourier coefficients does not oscillate too wildly, and thus it may behave well for the purposes of analyzing our algorithm.

Another possible direction for extension is related to the measurement of similarity and the sample complexity. Our current algorithm measures the similarity of rows $u$ and $v$ from their overlapping observed entries, which critically determines the sample complexity requirement of $np^{2}\gg 1$. However, for sparser datasets without overlaps, we may be able to reveal the similarity by instead comparing distribution signatures such as moments or comparing them through their “extended” neighborhoods [64].

As a concluding remark, we would like to mention that the latent variable model is a fairly general model and there is a large body of related applications. Some of the popular recent examples, which are special cases of latent variable model, include Stochastic blockmodels for community detection, the Bradley-Terry model for ranking from pair-wise comparison data and the Dawid-Skene model for low-cost crowd sourcing. Another prominent example of latent variable model is the generative model for random graphs referred to as a Graphon, which has been shown to be the limit of a sequence of graphs. We refer interested readers to [19, Section 2.4] for an excellent overview on the broad applicability of the latent variable model.

VI Extending Beyond Matrices to Tensors

A natural extension beyond matrix completion is the problem of completing a tensor of higher ($>2$) order. Given an unknown tensor $T$ of order $\tau$ with dimensions $n_{1}\times\dots\times n_{\tau}$, suppose that we observe a fraction of its entries corrupted by noise. Similar to matrix completion, the goal in tensor completion is to estimate the missing entries in the tensor from the noisy partial observations, as well as to “de-noise” the observed entries.

VI-A Short Background

The tensor completion problem is important within a wide variety of applications, including recommendation systems, multi-aspect data mining [65, 66], and machine vision [67, 68, 69].

Although tensor completion has been widely studied, there is still a wide gap in understanding, unlike matrix completion. This gap partially stems from the hardness of tensor decomposition, as most recovery methods rely on retrieving hidden algebraic structure through the framework of low-rank factorization. Tensors do not have a canonical decomposition such as the singular value decomposition (SVD) for a matrix.

There is a factorization scheme, namely the CANDECOMP/PARAFAC (CP) decomposition, which factorizes the tensor as a sum of rank-1 tensors (outer product of vectors). However, it is known that finding the rank of a tensor is NP-Complete, which makes it computationally intractable. Also, there are known ill-posedness [70] issues with CP-based low-rank approximation.

There are other kinds of decompositions such as the Tucker decomposition. Approaches based on Tucker decomposition essentially unfold (matricize or flatten) the tensor, and make use of matrix completion theory and methods [71, 72, 73, 67, 74].

VI-B Setup for Tensor Completion

VI-B1 Revisiting Our Model

The nonparametric model presented in Section II-A naturally extends beyond the bivariate case that corresponds to matrices, to multivariate setups encompassing higher-order tensors. We recap the latent variable model for a tensor following similar assumptions.

Suppose that there is an unknown $n_{1}\times\dots\times n_{\tau}$ tensor $T_{A}$ that we would like to estimate. We observe only a fraction of the total $\prod_{i=1}^{\tau}n_{i}$ entries of $T_{A}$ with some noise added. Let $\Omega\subset[n_{1}]\times\cdots\times[n_{\tau}]$ denote the index set of observed entries.

Precisely, we consider the following data generation model. Our data tensor $T_{Z}$ is an $n_{1}\times\dots\times n_{\tau}$ tensor that is generated as follows. Let $T_{M}\in\{0,1\}^{n_{1}\times\dots\times n_{\tau}}$ be a binary matrix, which we call the masking tensor. We let $T_{A}\in\mathbb{R}^{n_{1}\times\dots\times n_{\tau}}$ denote the signal tensor and $T_{N}\in\mathbb{R}^{n_{1}\times\dots\times n_{\tau}}$ denote the noise tensor. For each $\vec{\alpha}=(\alpha_{1},\ldots,\alpha_{\tau})\in[n_{1}]\times\cdots\times[n_{\tau}]$,

[TABLE]

Similar to the model assumptions for the matrix case in Section II-A1, we assume $T_{A}(\vec{\alpha})$ is generated by the following latent variable model equipped with certain regularity assumptions.

• •

Nonparametric model: there exists a latent function $f$ such that

[TABLE]

for all $\vec{\alpha}\in[n_{1}]\times\dots\times[n_{\tau}]$. Here, $x_{1}(\alpha_{1}),\dots,x_{\tau}(\alpha_{\tau})$ denote latent variables associated with index $\alpha_{i}$ in the $i$-th coordinate of the tensor $T_{A}$, respectively.

• •

Regularity Assumptions

• –

For each $i\in[\tau]$, the latent variables $x_{i}(\alpha_{i})\in\mathcal{X}_{i}$ for all $\alpha_{i}\in[n_{i}]$, where $(\mathcal{X}_{i},d_{i})$ is a metric space such that $\textrm{diam}~{}\mathcal{X}_{i}=\sup_{\alpha,\beta\in\mathcal{X}_{i}}d_{i}(\alpha,\beta)\leq D_{i}$. Moreover, $\mathcal{X}_{i}$ is equipped with a Borel probability measure $\mu_{i}$ and $x_{i}(\alpha_{i})$ is drawn i.i.d. according to $\mu_{i}$.

• –

Latent function $f$ is bounded, specifically that there exists a constant $D_{f}$ such that for all $\vec{\alpha}$, $|f(\vec{\alpha})|\leq D_{f}$.

• –

Latent function $f$ is $L$-Lipschitz in the sense that

[TABLE]

Note that our model assumptions for the noise matrix and the masking matrix are stated in an entrywise fashion and readily extends to their tensor analogues, cf. Section II-A.

VI-B2 Flattening a Tensor to a Matrix

A $\tau$-order tensor $T\in\mathbb{R}^{n_{1}\times\cdots\times n_{\tau}}$ can be viewed as a $\tau$-dimensional array of numbers. It is possible to ‘flatten’ the tensor $T$ to a matrix (i.e. 2-dimensional array) by rearranging the numbers in the $\tau$-dimensional array. Formally, given a set $U\subset[\tau]$, we define $\Phi_{U}(T)$ to be a $\prod_{i\in U}n_{i}$ by $\prod_{j\in[\tau]\setminus U}n_{j}$ matrix obtained by flattening the original tensor. Without loss of generality171717by taking transpose of $T$, we may assume $U=[\upsilon]$ for some $0\leq\upsilon\leq\tau$. We index the rows and the columns of $\Phi_{U}(T)$ using a $\upsilon$-tuple and a $(\tau-\upsilon)$-tuple, respectively. That is to say, for any $(\alpha_{1},\ldots,\alpha_{\upsilon})\in[n_{1}]\times\cdots\times[n_{\upsilon}]$ and $(\bar{\alpha}_{1},\ldots,\bar{\alpha}_{\tau-\upsilon})\in[n_{\upsilon+1}]\times\cdots\times[n_{\tau}]$, we have

[TABLE]

VI-C Tensor Completion Algorithm

We describe our generic recipe for tensor completion in Algorithm 3. We remark here that any matrix estimation algorithm can be used as the matrix estimation subroutine in Step 2.

VI-D Analysis

Suppose that Algorithm 2 is used as the matrix estimation subroutine in Algorithm 3. We can obtain an MSE upper bound for the tensor completion algorithm, which is similar to that stated in Theorem 1. We do not include a formal theorem statement and its proof here, but we briefly point out what remains unchanged and what needs to be modified in the proof of Theorem 1 (cf. Appendix A) to obtain its counterpart for tensor completion.

Given a $\tau$-order tensor $T_{Z}\in\mathbb{R}^{n_{1}\times\dots\times n_{\tau}}$ and an index set $U\subset[\tau]$, we consider the flattened matrix, $\Phi_{U}(T_{Z})\in\mathbb{R}^{\prod_{i\in U}n_{i}\times\prod_{j\in[\tau]\setminus U}n_{j}}$. Since this matrix is obtained by flattening $T_{Z}$, the ‘rows’ and the ‘columns’ of $\Phi_{U}(T_{Z})$ are not fully exchangeable – they satisfy only ‘partial’ exchangeability induced by the exchangeability in $T_{Z}$. As a result, we cannot assume the latent variables associated with $\prod_{i\in U}n_{i}$ rows of $\Phi_{U}(T_{Z})$ are drawn i.i.d. from a latent space in the current setup. In fact, there are only $\sum_{i\in U}n_{i}$ number of independent latent random variables ($n_{i}$ from $\mathcal{X}_{i}$ for each $i\in U$) associated to the $\prod{i\in U}n_{i}$ rows of the flattened matrix. Similarly, the latent variables associated to the columns of $\Phi_{U}(T_{Z})$ are also not sampled i.i.d. as they involve shared coordinates in the original tensor. This difference adds a complication to the analysis of Algorithm 3, but its effect is limited to the Step 1 in the proof of Lemma 2.

Following the discussion in Appendix A-A, we observe that there are $\sum_{i=1}^{\tau}n_{i}+2\prod_{i=1}^{\tau}n_{i}$ independent sources of randomness in our model for $T_{Z}$: \big{\{}x_{i}(\alpha_{i})\big{\}}_{i\in[\tau]\atop\alpha_{i}\in[n_{i}]}, \big{\{}N(\vec{\alpha})\big{\}}_{\vec{\alpha}\in[n_{1}]\times\dots\times[n_{\tau}]}, \big{\{}M(\vec{\alpha})\big{\}}_{\vec{\alpha}\in[n_{1}]\times\dots\times[n_{\tau}]}. For the sake of simplicity, we may assume $U=[\upsilon]$ for some $0\leq\upsilon\leq\tau$. For $\vec{\alpha}^{(1)}\in[n_{1}]\times\dots\times[n_{\upsilon}]$ and $\vec{\alpha}^{(2)}\in[n_{\upsilon+1}]\times\dots\times[n_{\tau}]$, we let \vec{x}_{\text{row}}(\vec{\alpha}^{(1)})=\big{(}x_{1}(\alpha^{(1)}_{1}),\dots,x_{\upsilon}(\alpha^{(1)}_{\upsilon})\big{)} and \vec{x}_{\text{col}}(\vec{\alpha}^{(2)})=\big{(}x_{\upsilon+1}(\alpha^{(2)}_{1}),\dots,x_{\tau}(\alpha^{(2)}_{\tau-\upsilon})\big{)}, respectively. Also, we let $\vec{1}$ denote a sequence of an appropriate length with all coordinates being $1$.

By the same exchangeability argument, we can upper bound the MSE by the sum of signal component and the noise component as in the proof of Theorem 1 – see Eq. (31), Lemmas 1 and 3 in Appendix A):

[TABLE]

It remains to bound each of the two terms in (30), which can be accomplished by concentration inequalities.

The second term in (30) can be bounded by a very similar argument as in Lemma 4, but we need to adjust for the fact that $\vec{x}_{\text{row}}(\cdot)$ is associated to a tuple of $U$ latent variables, that are shared across different rows. In Step 1 of the proof or Lemma 4, we would use Chernoff bound separately for each dimension of the tensor $i\in U$ to argue that for each entry $x_{1}(\alpha^{(1)}_{i})$, there are sufficiently many “nearest neighbor” coordinates with similar latent variables in the $i$-th dimension of the tensor. The number of nearest neighbor rows to $\vec{x}_{\text{row}}(\vec{\alpha}^{(1)})$ would then be lower bounded by the product of the number of nearest neighbors for each coordinate $i\in U$.

We also need a small modification in bounding the first term in (30). In Step 1 of the proof of Lemma 2, we use Bernstein’s inequality to prove the concentration of

[TABLE]

We observe that we cannot use Bernstein’s inequality in the current setup of tensor completion, because the summands, \big{(}T_{Z}(\vec{1},\vec{j})-T_{Z}(\vec{v},\vec{j})\big{)}^{2}, are no longer independent, unlike the setup of matrix completion.

Nevertheless, the dependence between the summands is still reasonably weak. We may consider $\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})$ as a function of the independent random variables, \big{\{}x_{i}(\alpha_{i})\big{\}}_{i\in[\tau]\atop\alpha_{i}\in[n_{i}]}, \big{\{}N(\vec{\alpha})\big{\}}_{\vec{\alpha}\in[n_{1}]\times\dots\times[n_{\tau}]}, \big{\{}M(\vec{\alpha})\big{\}}_{\vec{\alpha}\in[n_{1}]\times\dots\times[n_{\tau}]}. Then we are able to prove concentration of $\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})$ to its expectation, using a different tool, e.g., by modified log-Sobolev inequality. Once we show \textrm{dissim}_{\text{row}}(\vec{1},\vec{v})\approx\big{\|}f(x_{\text{row}}(\vec{v}),\cdot)-f(x_{\text{row}}(\vec{1}),\cdot)\big{\|}_{L^{2}}^{2}+2\sigma^{2} in a similar form as in (33), the rest of the proof of Lemma 2 can be reused.

Remark 6*.*

In order to quickly get a sense of the appropriate concentration result, we consider the case where $n_{i}=n$ ($n\geq 2$) for all $i\in[\tau]$, $p=1$, $|U|=\upsilon$ and the noise is bounded, i.e., $|N(\vec{\alpha})|\leq\gamma$ for all $\vec{\alpha}\in[n]^{\tau}$. We also assume $\textrm{diam}~{}\mathcal{X}_{i}\leq D$ for all $i=\upsilon+1,\dots,\tau$. Given $x_{\text{row}}(\vec{1})$ and $x_{\text{row}}(\vec{v})$, we can observe that $\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})$ is a function of \big{\{}x_{i}(\alpha_{i})\big{\}}_{i=\upsilon+1,\ldots,\tau\atop\alpha_{i}\in[n]}, \big{\{}N(\vec{1},\vec{\beta})\big{\}}_{\vec{\beta}\in[n]^{\tau-\upsilon}\atop\vec{\beta}\neq\vec{1}}, and \big{\{}N(\vec{v},\vec{\beta})\big{\}}_{\vec{\beta}\in[n]^{\tau-\upsilon}\atop\vec{\beta}\neq\vec{1}}. Note that $|\mathcal{B}_{\text{adj}}(\vec{1})\cap\mathcal{B}_{\text{adj}}(\vec{v})|=n^{\tau-\upsilon}-1$.

We consider the influence that each independent random variable exerts on the function $\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})$. It is easy to verify that $\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})$ satisfies the bounded difference property. To verify the property, let us introduce shorthand notations g\big{[}x_{i}(\alpha_{i})\big{]}\triangleq\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})\big{[}x_{i}(\alpha_{i});~{}\textrm{the rest}\big{]} and h\big{[}N(\vec{1},\vec{\beta})\big{]}\triangleq\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})\big{[}N(\vec{1},\vec{\beta});~{}\textrm{the rest}\big{]}. Then we observe that

[TABLE]

Applying the bounded difference inequality (e.g., McDiarmid’s), it follows that for any $t>0$,

[TABLE]

Note that we may assume $\tau-\upsilon\geq 1$ and therefore $\frac{2n}{n^{\tau-\upsilon}-1}\leq 4$ (because $n\geq 2$). As a result, we can expect \big{|}\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})-\mathbb{E}\textrm{dissim}_{\text{row}}(\vec{1},\vec{v})\big{|}\lesssim\frac{(L^{2}D^{2}+\gamma^{2})\sqrt{\tau-\upsilon}}{\sqrt{n}}. This concentration inequality is $\sqrt{\tau-\upsilon}$ times weaker than the result (Bernstein’s inequality) from Step 1 of the proof of Lemma 2, but its dependence on $n$ remains the same.

Remark 7*.*

The sample complexity of our matrix estimation algorithm requires that $p\geq\max(m^{-1+\delta},n^{-1/2+\delta})$, where $m$ is the number of rows and $n$ is the number of columns. Therefore, the optimal flattening of the tensor that would minimize sample complexity is the flattening that such that

[TABLE]

If $n_{i}=n$ for all $i\in[\tau]$, then the optimal flattening would result in matrix dimensions of $n^{\lfloor\tau/3\rfloor}\times n^{\lceil 2\tau/3\rceil}$. In this case the natural extension of our analysis should result in MSE convergence with sample complexity of $p\geq n^{-\lfloor\tau/3\rfloor+\delta}$.

VII Experiments

In this section we present experimental results from applying the User-Item Gaussian Kernel variant of our algorithm to real datasets. We state in Algorithm 4 the specific algorithm variant that is used in the experiments. We did not sample split as it was primarily used for the analysis but is not necessary in practice. The implemented algorithm adds a few modifications:

• •

$\textrm{dissim}_{\text{row}}(u,v)$ is computed according to the sample variance of the differences between two rows

• •

we use Gausian kernel weights, and we combine both row and column dissimilarities by taking the maximum

• •

we use an estimator $A(u,i)\approx(Z(v,i)+Z(u,j)-Z(v,j))$ motivated by first order Taylor series approximation

The resulting algorithm is similar to the mean-adjusted variant of collaborative filtering, except for the addition of combining both row and column dissimilarities.

VII-A Matrix completion

We evaluated the performance of our algorithm on predicting user-movie ratings for the MovieLens 1M and Netflix datasets. We chose the overlap parameter $\beta=2$ to ensure the algorithm is able to compute an estimate for all missing entries. When $\beta$ is larger, the algorithm enforces rows (or columns) to have more commonly rated movies (or users). Although this increases the reliability of the estimates, it also reduces the fraction of entries for which the estimate is defined.

We compared our method with user-user collaborative filtering, item-item collaborative filtering, and SoftImpute from [75]. We chose the classic mean-adjusted collaborative filtering method, in which the weights are proportional to the cosine similarity of pairs of users or items (i.e. movies). SoftImpute is a matrix-factorization-based method which iteratively replaces missing elements in the matrix with those obtained from a soft-thresholded SVD.

The MovieLens 1M data set contains about 1 million ratings by 6000 users of 4000 movies. The Netflix data set consists of about 100 million movie ratings by 480,189 users of about 17,770 movies. For both MovieLens and Netflix data sets, the ratings are integers from 1 to 5. From each dataset, we generated 100 smaller user-movie rating matrices, in which we randomly subsampled 2000 users and 2000 movies. For each rating matrix, we randomly select and withhold a percentage of the known ratings for the test set, while the remaining portion of the data set is revealed to the algorithm for computing the estimates (or training). After the algorithm computes its predictions for all the missing user-movie pairs, we evaluate the Root Mean Squared Error (RMSE) of the predictions compared to the ratings from the withheld test set. Figure 1 plots the RMSE of our method along with classic collaborative filtering and SoftImpute evaluated against $10\%$, $30\%$, $50\%$, and $70\%$ withheld test sets. The RMSE is averaged over 100 subsampled rating matrices, and $95\%$ confidence intervals are provided.

Figure 1 suggests that our algorithm achieves a systematic improvement over classical user-user and item-item collaborative filtering. SoftImpute performs worse than all methods on the MovieLens dataset, but it performs better than all methods on the Netflix dataset.

VII-B Tensor completion

We consider the problem of image inpainting for evaluating the performance of tensor completion algorithm. Inpainting is the process of reconstructing lost or deteriorated parts of image or videos. Such methods, in particular, have revitalized the process of recovery old artifacts in museum world which was historically done by conservators or art restorers. An interested reader is referred to a recent survey [69] for summary of the state of art on methods and techniques. We compare performance of our algorithm against existing methods in the literature on the image inpainting problem. Figure 2 shows a sample of the image inpainting results for the facade and pepper images when $70\%$ of the pixels are removed.

An image can be represented as a $3^{\text{rd}}$-order tensor where the dimensions are rows $\times$ columns $\times$ RGB. In particular we used three images (Lenna, Pepper, and Facade) of dimensions 256 $\times$ 256 $\times$ 3. For each image, a percentage of the pixels are randomly removed, and the missing entries are filled in by various tensor completion algorithms.

For the implementation of our tensor completion method, we collapsed the last two dimensions of the tensor (columns and RGB) to reduce the image to a matrix, and applied our method. We set the overlap parameter $\beta=2$. We compared our method against fast low rank tensor completion (FaLRTC) [76], alternating minimization for tensor completion (TenAlt) [34], and fully Bayesian CP factorization (FBCP) [77], which extends the CANDECOMP/PARAFAC(CP) tensor factorization with automatic tensor rank determination.

To evaluate the outputs produced by each method, we computed the relative squared error (RSE), defined as

[TABLE]

where $\bar{Z}$ is the average value of the true entries. Figure 3 plots the RSE achieved by each tensor completion method on the three images, as a function of the percentage of pixels removed. The results demonstrate that our tensor completion method is competitive with existing tensor factorization based approaches, while maintaining a naive simplicity. In short, a simple algorithm works nearly as good as the best algorithm for this problem!

Acknowledgment

This work is supported in parts by ARO under MURI award 133668-5079809, by NSF under grants CMMI-1462158, CMMI-1634259 and TRIPODS, by DARPA under grant W911NF-16-1-0551, and additionally by a Samsung Scholarship, Siebel Scholarship, NSF Graduate Fellowship, and Claude E. Shannon Research Assistantship.

Appendix A Proof of Theorem 1

A-A Outline of the Proof

There are $m+n+2mn$ independent sources of randomness in our model, \big{\{}x_{\text{row}}(u)\big{\}}_{u\in[m]}, \big{\{}x_{\text{col}}(i)\big{\}}_{i\in[n]},\big{\{}N(u,i)\big{\}}_{(u,i)\in[m]\times[n]}, \big{\{}M(u,i)\big{\}}_{(u,i)\in[m]\times[n]}. We let $\Theta$ denote the collection of these random variables.

Due to the exchangeability of the model and the linearity of the expectation,

[TABLE]

In Theorem 1, we provide an upper bound on \mathbb{E}\left[\big{(}\widehat{A}(1,1)-A(1,1)\big{)}^{2}~{}\big{|}~{}x_{\text{row}}(1)\right], conditioned on the latent feature of the first row.

Recall from Algorithm 2 that

[TABLE]

in our user-user fixed radius nearest neighbor algorithm when $|\mathcal{B}_{\text{est}}(1,1)|\geq 1$ and $\widehat{A}(1,1)=0$ when $|\mathcal{B}_{\text{est}}(1,1)|=0$. Therefore,

[TABLE]

By the assumption that the magnitude of the latent function $f$ is bounded by $D_{f}$,

[TABLE]

and that

[TABLE]

where $C>0$ is an absolute constant (coming from $\|A\|_{\psi_{2}}\leftrightarrow\text{Var}(A)$). The inequality in (a) follows from Lemmas 1 and 3.

We need local properties of the probability measure on the latent space to upper bound the two terms in the last line; see Lemmas 2 and 4. We define the function $\phi$ for $r\geq 0$ and $x\in\mathcal{X}_{\textrm{row}}$

[TABLE]

We prove in Lemma 5 that

[TABLE]

A-B Upper Bounding the Contribution of Signal on MSE

Lemma 1**.**

[TABLE]

Proof.

Recall from Section A-A that we let $\Theta$ denote the collection of $m+n+2mn$ independent sources of randomness in our model, \big{\{}x_{\text{row}}(u)\big{\}}_{u\in[m]}, \big{\{}x_{\text{col}}(i)\big{\}}_{i\in[n]},\big{\{}N(u,i)\big{\}}_{(u,i)\in[m]\times[n]}, \big{\{}M(u,i)\big{\}}_{(u,i)\in[m]\times[n]}.

By the tower property of expectation,

[TABLE]

We investigate the conditional expectation, again using the tower property and the fact that $|\mathcal{B}_{\text{est}}(1,1)|$ is fully determined when we condition on $\Theta\setminus x_{\text{col}}(1)$. For the sake of readability, \mathbb{I}\big{(}{|\mathcal{B}_{\text{est}}(1,1)|\geq 1}\big{)} is omitted in the lines below.

[TABLE]

Equation (a) follows from the conditional independence; (b) follows from Fubini’s theorem and the independence between $\theta_{\textrm{col}}{1}$ and $|\mathcal{B}_{\text{est}}(1,1)|$; and (c) follows from Cauchy-Schwarz inequality.

Consequently, it follows that

[TABLE]

∎

Lemma 2**.**

Let $\eta\geq 2\sigma^{2}$. The following inequality holds for user-user fixed radius neighbor algorithm:

[TABLE]

In the above expression, K\triangleq\Big{(}\frac{D_{f}}{\sqrt{\ln 2}}+2\sigma\Big{)}^{2} and $C,c,c^{\prime}>0$ are absolute constants.

Proof.

Step 1: Choose $v\in[m]\setminus\{1\}$ such that $(v,1)\in\mathcal{B}_{\text{est}}(1,1)$. Observe that

[TABLE]

Note that Z(1,j)-Z(v,j)=\big{(}A(1,j)-A(v,j)\big{)}+\big{(}N(1,j)-N(v,j)\big{)}. By the model assumptions, we have

[TABLE]

Therefore (cf. [78, Lemma 2.7.6]),

[TABLE]

For notational conciseness, we let $f_{1v}^{2}$ denote the conditional expectation of $\textrm{dissim}_{\text{row}}(1,v)$:

[TABLE]

Now, suppose that $|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|=N_{v}$. We observe that $\mathbb{E}\left[\textrm{dissim}_{\text{row}}(1,v)|x_{\text{row}}(1),x_{\text{row}}(v),|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|=N_{v}\right]=f_{1v}^{2}$ regardless of $N_{v}$. Then it follows from Bernstein’s inequality (cf. [78, Theorem 2.8.2]) that for any $t\geq 0$,

[TABLE]

Note that we provided an upper bound on the conditional probability in (33) that holds for any realization of $x_{\text{row}}(1)$ and $x_{\text{row}}(v)$, and thus we could remove the conditioning on $x_{\text{row}}(v)$. By plugging in the value of $f_{1v}^{2}$, this leads to the following inequality: for any $t\geq 0$,

[TABLE]

Step 2: Let $v^{*}\in[m]\setminus\{1\}$ denote the maximizer such that $(v^{*},1)\in\mathcal{B}_{\text{est}}(1,1)$ and

[TABLE]

According to the description of our user-user fixed radius nearest neighbor algorithm, $\textrm{dissim}_{\text{row}}(1,v^{*})\leq\eta$.

[TABLE]

Now we observe that

[TABLE]

Recall that for $a>0$ and $b\in\mathbb{R}$,

[TABLE]

Step 3: Observe that $|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|\sim\textrm{Binomial}(n-1,p^{2})$. By the binomial Chernoff theorem,

[TABLE]

Therefore, with the shorthand notation $N_{v}=|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|$, we can see that

[TABLE]

The inequality in (a) follows from that $N_{v}\geq 1$ and $\frac{Ke^{-cN_{v}}}{cN_{v}}$ is maximized when $N_{v}=1$. ∎

A-C Upper Bounding the Contribution of Noise on MSE

Lemma 3**.**

[TABLE]

Proof.

By the tower property of expectation,

[TABLE]

(a) follows from \mathbb{E}\Big{[}|N(v,1)|^{2}~{}\big{|}~{}x_{\text{row}}(1),|\mathcal{B}_{\text{est}}(1,1)|\Big{]}=\mathbb{E}\big{[}|N(v,1)|^{2}\big{]}\leq C\sigma^{2}, as a result of $\|N(v,1)\|_{\psi_{2}}\leq\sigma$. ∎

Lemma 4**.**

Let $\eta^{\prime}\geq 2\sigma^{2}$ and \eta\geq\eta^{\prime}+K\max\Big{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\Big{)}. The following inequality holds for user-user fixed radius nearest neighbor algorithm:

[TABLE]

where \phi(x,r)=\mathbb{P}_{x_{\text{row}}(v)\sim\mu_{\mathcal{X}_{\textrm{row}}}}\left(\big{\|}f(x,\cdot)-f(x_{\text{row}}(v),\cdot)\big{\|}_{L^{2}}^{2}\leq r\right) for any $x\in\mathcal{X}_{\textrm{row}}$ and $r>0$.

Proof.

Step 1: Our interest is in bounding conditional expectation of \frac{1}{|\mathcal{B}_{\text{est}}(1,1)|}~{}\mathbb{I}\big{(}|\mathcal{B}_{\text{est}}(1,1)|\geq 1\big{)} given $x_{\text{row}}(1)\in\mathcal{X}_{\textrm{row}}$. To that end, let us fix $x_{\text{row}}(1)$. Let \eta^{\prime}=\eta-K\max\Big{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\Big{)}. For each $v\in[m]\setminus\{1\}$, the following two conditions

[TABLE]

are satisfied with success probability181818Note that the first condition solely depends on $M(v,1)$, while the second condition depends only on $x_{\text{row}}(1),x_{\text{row}}(v)$, hence, they are independent events. p\cdot\phi\big{(}x_{\text{row}}(1),\eta^{\prime}-2\sigma^{2}\big{)}.

Let $\Omega_{\text{good row}}\subset[m]\setminus\{1\}$ denote the set of row indices that satisfy the two conditions described in (37). By the binomial Chernoff bound,

[TABLE]

Step 2: Next, we want to show that $\{(v,1):v\in\Omega_{\text{good row}}\}\subset\mathcal{B}_{\text{est}}(1,1)$ with high probability. For that purpose, we first require $|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|$ to be sufficiently large for all $v\in\Omega_{\text{good row}}$. Observe that for each $v\in\Omega_{\text{good row}}$, $|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v)|$ is distributed following $\textrm{Binomial}(n-1,p^{2})$. Again by the binomial Chernoff bound (as in (36)), we have

[TABLE]

Since $|\Omega_{\text{good row}}|\leq m-1$, it follows from the union bound that

[TABLE]

By construction, for all $v\in\Omega_{\text{good row}}$,

[TABLE]

Note that we can obtain the following concentration inequality by similar arguments191919This provides a probabilistic tail bound on the opposite side of that in (33). The proof remains the same. as in (33):

[TABLE]

With another application of the union bound, this implies that

[TABLE]

Therefore202020Recall that $\eta^{\prime}=\eta-t$ with t=K\max\Big{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\Big{)}.,

[TABLE]

Observe that

[TABLE]

which implies that $|\mathcal{B}_{\text{est}}(1,1)|\geq|\Omega_{\text{good row}}|$.

Step 3: Let

[TABLE]

By the law of total probability,

[TABLE]

where (a) follows from \frac{1}{|\mathcal{B}_{\text{est}}(1,1)|}~{}\mathbb{I}\big{(}|\mathcal{B}_{\text{est}}(1,1)|\geq 1\big{)}\leq 1; and (b) follows from that (42) and

[TABLE]

∎

A-D Upper Bounding the Probability of $|\mathcal{B}_{\text{est}}(1,1)|=0$

Lemma 5**.**

Let $\eta^{\prime}\geq 2\sigma^{2}$ and \eta\geq\eta^{\prime}+K\max\Big{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\Big{)}. The following inequality holds for user-user fixed radius nearest neighbor algorithm:

[TABLE]

Proof.

Recall the definition of $\Omega_{\text{good row}}$ from (37):

[TABLE]

Note that we obtain in (38) that

[TABLE]

by applying the binomial Chernoff bound.

Now we observe that if there exists at least one $v\in\Omega_{\text{good row}}$, for which

[TABLE]

then $|\mathcal{B}_{\text{est}}(1,1)|\geq 1$ due to the Lipschitzness assumption on $f$. Therefore, we have

[TABLE]

where $v_{0}=\min\{v\in\Omega_{\text{good row}}\}$.

We bound the two terms in (43) separately. First, as long as \phi\Big{(}x_{\text{row}}(1),\eta^{\prime}-2\sigma^{2}\Big{)}>0, it follows from (38) that

[TABLE]

Second, by the usual trick of total probability (recall from (32) that $f_{1v_{0}}^{2}$ denotes $\mathbb{E}\left[\textrm{dissim}_{\text{row}}(1,v_{0})|x_{\text{row}}(1),x_{\text{row}}(v_{0})\right]$),

[TABLE]

where (a) follows from (33) and the independence between $|\mathcal{B}_{\text{adj}}(1)\cap\mathcal{B}_{\text{adj}}(v_{0})|$ and $x_{\text{row}}(1)$; and (b) follows from (36) and the assumption that \eta-\eta^{\prime}\geq K\max\Big{(}\sqrt{\frac{4\log(m-1)}{c(n-1)p^{2}}},\frac{4\log(m-1)}{c(n-1)p^{2}}\Big{)}.

All in all,

[TABLE]

∎