Matrix estimation by Universal Singular Value Thresholding

Introduction

Consider a statistical estimation problem where the unknown parameter is not a single value or vector, but an $m\times n$ matrix $M$. Given an estimator $\hat{M}$, one choice for a measure of the error in estimation is the mean-squared error, defined as

Here, $\hat{m}_{ij}$ and $m_{ij}$ denote the $(i,j)$th elements of $\hat{M}$ and $M$, respectively. If we have a sequence of such problems, and $M_{n}$ and $\hat{M}_{n}$ denote the parameter and the estimator in the $n$th problem, then by usual statistical terminology we may say that the sequence of estimators $\hat{M}_{n}$ is consistent if

The problem of estimating the entries of a large matrix from incomplete and/or noisy entries has received widespread attention ever since the proliferation of large data sets. Early work using spectral analysis was done by a number of authors in the engineering literature, for example, by Azar et al. azaretal01 and Achlioptas and McSherry achlioptas01 . This was followed by a sizable body of work on spectral methods, the main pointers to which may be found in the important recent papers of Keshavan, Montanari and Oh kmo10a , kmo10b . Nonspectral methods also appeared, for example, in renniesrebro05 .

In a different direction, statisticians have worked on matrix completion problems under a variety of modeling assumptions. Possibly the earliest works are due to Fazel fazel02 and Rudelson and Vershynin rv07 . The emergence of compressed sensing donoho06 , candesrombergtao06 has led to an explosion in activity in the field of matrix estimation and completion, beginning with the work of Candès and Recht candesrecht09 . The pioneering works of Emmanuel Candès and his collaborators candesrecht09 , candestao10 , candesplan10 , ccs introduced the technique of matrix completion by minimizing the nuclear norm under convex constraints, which is a convex optimization problem tractable by standard algorithms. This method has the advantage of exactly, rather than approximately, recovering the entries of the matrix when a suitable low rank assumption is satisfied, together with a certain other assumption called “incoherence.”

Since the publication of candesrecht09 , a number of statistics papers have attacked the matrix completion problem from various angles. Some notable examples are negahban , mht , klt , rohdetsybakov11 , kol2012 , davenport . In a different direction, a paper that seems to have a close bearing on the analytical aspects of this paper is a manuscript of Oliveira oliveira09 .

In addition to the theoretical advances, a large number of algorithms for matrix completion and estimation have emerged. The main ones are nicely summarized and compared in mht .

The purpose of this paper is to introduce a new estimator that is capable of solving a variety of matrix estimation problems that are not tractable by existing tools (at least in a mathematically provable sense). The estimator and its properties are described in this introductory section. Section 2 focuses on applications, which include applications to low rank matrices, stochastic blockmodels, distance matrices, latent space models, positive definite matrices, graphons and generalized Bradley–Terry models. All proofs are in Section 3. An expanded version (version 5) of the paper containing more theorems, examples and simulation results is available on arXiv at the URL: http://arxiv.org/pdf/1212.1247v5.pdf.

For interesting new developments that appeared after the first draft of this paper was posted on arXiv, see choi , dg , nadakuditi . Further references and citations are given in subsequent sections.

Similarly, one can define the “skew-symmetric model,” where the difference $X-M$ is skew-symmetric, with independence on and above the diagonal as in the symmetric model. This model is used for analyzing the nonparametric Bradley–Terry model in Section 2.7.

2 The USVT estimator

In the above models, we construct an estimator $\hat{M}$ of $M$ based on the observed entries of $X$ along the following steps. Tentatively, I call this the Universal Singular Value Thresholding (USVT) algorithm. {longlist}[5.]

For each $i,j$, let $y_{ij}=x_{ij}$ if $x_{ij}$ is observed, and let $y_{ij}=0$ if $x_{ij}$ is unobserved. Let $Y$ be the matrix whose $(i,j)$th entry is $y_{ij}$.

Let $Y=\sum_{i=1}^{m}s_{i}u_{i}v_{i}^{T}$ be the singular value decomposition of $Y$. (In the symmetric and skew-symmetric models, $m=n$.)

Let $\hat{p}$ be the proportion of observed values of $X$. In the symmetric and skew-symmetric models, let $\hat{p}$ be the proportion of observed values on and above the diagonal.

Choose a small positive number $\eta\in(0,1)$ and let $S$ be the set of “thresholded singular values,” defined as

[Note: (a) In simulations, the method described below seemed to work even if $\eta$ was taken to be exactly equal to zero; but the mathematical proof that I have requires $\eta$ to be positive. In practice, one may choose $\eta$ a priori to be some arbitrary small positive number, say, $0.01$; but a data-dependent choice is not allowed. (b) If it is known that $\operatorname{Var}(x_{ij})\leq\sigma^{2}$ for all $i,j$, where $\sigma$ is a known constant${}\leq 1$, then the threshold $(2+\eta)\sqrt{n\hat{p}}$ may be improved to $(2+\eta)\sqrt{n\hat{q}}$, where $\hat{q}:=\hat{p}\sigma^{2}+\hat{p}(1-\hat{p})(1-\sigma^{2})$.]

Let $w_{ij}$ denote the $(i,j)$th element of $W$. Define

Let $\hat{M}$ be the matrix whose $(i,j)$th entry is $\hat{m}_{ij}$.

If the entries of $M$ and $X$ are known to belong to an interval $[a,b]$ instead of $$, then subtract$(a+b)/2$from each entry of$X$and divide by$(b-a)/2$, so that the entries are forced to lie in$$, then apply the above procedure, and finally multiply the end-result by$(b-a)/2$and add$(a+b)/2$to get the estimate of$M$.

If $m>n$, then one should work with $M^{T}$ and $X^{T}$ instead of $M$ and $X$, so that the number of rows is forced to be $\leq$ the number of columns.

3 Main result

Recall that the nuclear norm of $M$, written $\|M\|_{*}$, is defined as the sum of the singular values of $M$. Recall also the definition (1) of the mean squared error of a matrix estimator. The following theorem gives an error bound for the estimator $\hat{M}$ in terms of the nuclear norm of $M$. This is the main result of this paper.

Let $\hat{M}$ and $M$ be as above. Let $\operatorname{MSE}(\hat{M})$ be defined as in (1). Suppose that $p\geq n^{-1+\varepsilon}$ for some $\varepsilon>0$. Then

where $C$ and $c$ are positive constants that depend only on the choice of $\eta$ and $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$. The same result holds for the symmetric and skew-symmetric models, after putting $m=n$.

Moreover, if in the same setting as above, we know that $\operatorname{Var}(x_{ij})\leq\sigma^{2}$ for all $i,j$ for some known $\sigma^{2}\leq 1$, and the threshold is set at $(2+\eta)\sqrt{n\hat{q}}$ (see step $4$ of the algorithm), the same result holds under the condition that $q\geq n^{-1+\varepsilon}$, where $q:=p\sigma^{2}+p(1-p)(1-\sigma^{2})$. In this case the exponential term in the error changes to $C(\varepsilon)e^{-cnq}$ and the term $\|M\|_{*}/(m\sqrt{np})$ improves to $\|M\|_{*}\sqrt{q}/(m\sqrt{n}p)$.

Incidentally, the proof shows that the condition $p>n^{-1+\varepsilon}$ may be improved to $p>n^{-1}(\log n)^{6+\varepsilon}$ (see Theorem 18), but I prefer to retain the present version for aesthetic reasons, especially considering that it is not a real improvement from any practical point of view.

It should be emphasized that although singular value thresholding has been used in a number of papers on matrix completion and estimation (see, e.g., azaretal01 , achlioptas01 , ccs , kmo10a , kmo10b and references therein), the above algorithm has the unique feature that the threshold is universal. In the literature, it is usually assumed that the matrix $M$ has a rank $r$ that is known, and uses the value of $r$ while thresholding. The USVT algorithm manages to cut off the singular values at the “correct” level, depending on the structure of the unknown parameter matrix. The adaptiveness of the USVT threshold is somewhat similar in spirit to that of the SureShrink algorithm of Donoho and Johnstone dj95 . SureShrink performs function estimation by estimating Fourier coefficients in some suitable basis, and then thresholds the coefficients at a threshold that automatically adapts to the smoothness of the unknown function. Analogously, the USVT algorithm computes the eigenvalues of the observed matrix, and then thresholds the eigenvalues at a universal threshold that is automatically adaptive in nature, because it picks out as much “structure” as is available and throws out all the randomness. This point will become more clear from the examples discussed in Section 2.

One limitation of USVT is the requirement that the entries should lie in a bounded interval. One may relax this requirement by assuming, for example, that the errors $x_{ij}-m_{ij}$ are distributed as normal random variables with mean zero and variance $\sigma^{2}$. If $\sigma^{2}$ is known, then I believe that one can modify the USVT algorithm by thresholding at $(2+\eta)\sigma\sqrt{n}$ and obtain the same theorems. The rationale behind this belief is as follows: if $A$ is a large symmetric random matrix whose entries on and above the diagonal are independent, have zero mean, and are bounded by $1$ in absolute value, then the spectral norm of $A$ is less than $2+\eta$ with high probability. This is the key ingredient in the proof of Theorem 1. But such a result continues to be true, after replacing $2+\eta$ with $(2+\eta)\sigma$, if the entries are normally distributed with mean zero and variance bounded by $\sigma^{2}$. Therefore, it is conceivable that the proof of Theorem 1 may be modified to accommodate this altered situation. However, if $\sigma^{2}$ is unknown, I do not know how to proceed. In reality, $\sigma^{2}$ will not be known; this is why I have not worked with the normality assumption. Also, the situation of normally distributed entries but with a large proportion missing, seems to be trickier.

4 Minimax lower bound

where $c$ is a positive universal constant. Moreover, if $p<1/2$ then $X$ and $M$ may be chosen such that $X=M$. The same lower bound holds in the symmetric case and in the skew-symmetric case.

It is worth noting that the exponentially small discrepancy is necessary. For example, if $\delta=0$, then the minimax error is obviously zero. However, there is still an exponentially small chance that $\hat{M}$ may be nonzero. It is also worth noting that if $\delta$ is not too small (e.g., if $\delta>\sqrt{m/p}$), then the exponential discrepancy does not matter, and the combination of Theorems 1 and 2 gives the correct minimax error up to a universal multiplicative constant.

An examination of the proof of Theorem 1 indicates that with slight modifications, one may obtain bounds on tail probabilities instead of an upper bound on the mean squared error. I have retained the present version for aesthetic reasons.

Incidentally, two notable recent papers, namely, Koltchinskii et al. klt and Davenport et al. davenport , have suggested matrix estimation by nuclear norm penalization and proved minimax optimality results that match up to logarithmic factors. Davenport et al. davenport , Theorem 3, show (in the notation of our Theorem 2) that if the entries of $X$ belong to $\{-1,1\}$ and if $\delta\geq 4\sqrt{mn}$, then the minimax error is bounded below by a universal constant times $\min\{\delta/(m\sqrt{np}),1\}$, provided that this quantity is bigger than $\delta^{2}/(m^{2}n)$. This is almost the same as the conclusion of Theorem 2, except that it does not cover the case of $\delta$ smaller than $4\sqrt{mn}$. Section 3.1 of davenport gives a matrix estimation algorithm based on nuclear norm penalization that achieves this minimax rate up to a logarithmic factor. However, the implementation of this algorithm requires that the user has a reasonable estimate for the nuclear norm of the unknown matrix $M$, since that is used as the regularization parameter. USVT has no such requirement. Another advantage that USVT has over the algorithm of davenport is that it may be easier to implement, especially for very large matrices, because it does not involve convex optimization.

The estimator of Koltchinskii et al. klt is also based on nuclear norm penalization: translating to our notation, they estimate $M$ by minimizing $\|X-\hat{M}\|_{F}^{2}+\lambda\|\hat{M}\|_{*}$ over all $\hat{M}$, where $\|\cdot\|_{F}$ is Frobenius norm, $\|\cdot\|_{*}$ is nuclear norm, and $\lambda$ is a regularization parameter. It is shown in klt that this problem is actually equivalent to soft singular value thresholding, where the threshold depends on the parameter $\lambda$. A conservative choice of $\lambda$ (albeit with an unspecified constant) and a minimax lower bound that matches the upper bound up to a logarithmic factor are given in klt . The minimax bound is computed over the set of all matrices with rank less than a given number and, therefore, is not directly comparable to the minimax bound in Theorem 2. With a suitable choice of $\lambda$—but again with unspecified constants—the upper bound in klt , Theorem 3, becomes (up to a logarithmic factor) essentially equal to $\|M\|_{*}/(m\sqrt{np})$. Note that this is the same as the main term in Theorem 1. However, if we additionally know that $M$ has low rank, then the upper bound in klt , Theorem 3, becomes substantially better (see Section 2.1).

5 Practical issues and warnings

I do not consider the USVT algorithm as presented above to be in a form that may implemented “as is.” This is mainly for the following reasons: {longlist}[(a)]

USVT is minimax optimal only up to a constant factor. In fact, it is very likely that one may be able to build a better estimator by taking into account the ratio $m/n$, and getting improved bounds when this ratio is small. Although Theorem 2 shows that the improvement will be limited to multiplication by a constant factor, such an improvement may be important for practical purposes. The recent paper dg has explored the issue of attaining the minimax error all the way up to the correct constant.

The number $\eta$ is a “tuning parameter” for this algorithm, that may be chosen by the implementer. The theorem is valid with any choice of $\eta$ in the interval $(0,1)$, although the constants in the error bounds blow up as $\eta$ tends to zero. I have noticed in simulations that taking $\eta=0$ works quite well, but I do not know how to prove that. Choosing $\eta$ to be a small but fixed positive number such as $0.01$ is consistent with the requirements of Theorem 1 and seemed to give good results in simulations. Choosing $\eta$ in a data-dependent manner is, however, not covered by Theorem 1.

Note that in practice, any data matrix may be centered and scaled so that the entries are forced to lie in the interval $$. However, if the centering and scaling are done in a data-dependent manner, then the assertion of Theorem 1 is no longer guaranteed to be true.

6 An impossibility theorem for error estimates

Theorem 1 gives an upper bound on the mean squared error of $\hat{M}$. The estimate involves the nuclear norm of parameter matrix $M$. A natural question is: Is it possible to estimate the true MSE of $\hat{M}$ from the data?

A straightforward approach is to use parametric bootstrap. Having estimated $M$ using $\hat{M}$, one may choose a large number $K$, generate $K$ copies of the data using $\hat{M}$ as the parameter matrix, compute the estimates $\hat{M}^{(i)}$, $i=1,\ldots,K$ for the $K$ simulations, and estimate the MSE of $\hat{M}$ using the bootstrap estimator

For the validity of the bootstrap estimate of the MSE, it is essential that the original $\hat{M}$ is an accurate estimate of $M$. In other words, we need to know a priori that $\operatorname{MSE}(\hat{M})$ is small to be able to claim that the bootstrap estimator of $\operatorname{MSE}(\hat{M})$ is accurate. Theorem 1 implies that if we know that $\|M\|_{*}$ is small enough from assumptions, this is true.

In the above setting, the following theorem establishes the impossibility of the existence of a good estimator for the MSE.

There cannot exist a good procedure for estimating the mean squared error of a nontrivial estimator.

Applications

Throughout this section, $m$, $n$, $M$, $X$, $p$ and $\hat{M}$ will be as in Section 1. Just to remind the reader, $M$ is an $m\times n$ matrix where $1\leq m\leq n$. The entries of $M$ are assumed to be bounded by $1$ in absolute value. The matrix $X$ is a random matrix whose entries are independent, and the $(i,j)$th element $x_{ij}$ has expected value equal to $m_{ij}$, the $(i,j)$th entry of $M$. Moreover, they satisfy $|x_{ij}|\leq 1$ with probability one. In particular, $X$ may be exactly equal to $M$, with no randomness. Each entry of $X$ is observed with probability $p$ and unobserved with probability $1-p$, independently of other entries. Occasionally, we will assume the symmetric model, where $m=n$, and the matrices $M$ and $X$ are symmetric. In the special case of the Bradley–Terry model in Section 2.7, we will assume the skew-symmetric model, where $X-M$ is skew-symmetric.

We will now work out various specific cases where Theorem 1 gives useful results.

Estimating low rank matrices has been the focus of the vast majority of prior work azaretal01 , achlioptas01 , fazel02 , renniesrebro05 , rv07 , candesrecht09 , candesplan10 , candestao10 , ccs , negahban , kmo10a , kmo10b , klt , kol2012 , mht . Theorem 1 works for low rank matrices. The following theorem, which is a simple corollary of Theorem 1, shows that $\hat{M}$ is a good estimate whenever the rank of $M$ is small compared to $mp$ (after assuming, as in Theorem 1, that $p\geq n^{-1+\varepsilon}$).

Suppose that $M$ has rank $r$. Suppose that $p\geq n^{-1+\varepsilon}$ for some $\varepsilon>0$. Then

where $C$ and $c$ depend only on $\eta$ and $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$. Moreover, the same result holds when $M$ and $X$ are symmetric.

The term $1/np$ in the error bound is necessary to take care of the case $r=0$. Even if $M$ is identically zero, the estimator $\hat{M}$ will incur some error due to the (possible) randomness in $X$.

Let us now inspect how the condition $r\ll mp$ compares with available bounds. In a notable sequence of papers, Keshavan, Montanari and Oh kmo10a , kmo10b obtain the same condition but only if $m$ and $n$ are comparable and the rank is known. Theorem 4, on the other hand, works even for “very rectangular” matrices where $m\ll n$ and the rank is unknown.

Candès and Tao candestao10 obtain the condition $r\ll mp$ with an extra poly-logarithmic term in the error. Moreover, they too require that $m$ and $n$ be comparable, and additionally they need the so-called “incoherence condition”. However, as noted before, the incoherence condition allows exact recovery, while our approach only gives approximate recovery.

The recent important work of Davenport et al. davenport gives an estimator with an error bound that is almost the same as that given by Theorem 4, but with a complicated optimization algorithm.

Theorem 4, however, is probably not an optimal result. It has been shown by Koltchinskii et al. klt , Theorems 3 and 5, that the true minimax error rate for a closely related problem is actually $r/mp$, up to a logarithmic factor.

The following theorem shows that the condition $r\ll mp$ is necessary for estimating $M$.

where $C$ is a positive universal constant and $[m/r]$ is the integer part of $m/r$.

2 The stochastic blockmodel

Consider an undirected graph on $n$ vertices. A stochastic blockmodel assumes that the vertices $1,\ldots,n$ are partitioned into $k$ blocks, and the probability that vertex $i$ is connected to vertex $j$ by an edge depends only on the blocks to which $i$ and $j$ belong. As usual, edges are independent of each other. Let $M$ be the matrix whose $(i,j)$th element is the probability of an edge existing between vertices $i$ and $j$. The matrix $X$ here is the adjacency matrix of the observed graph. Here, all elements of $X$ are observed, so $p=1$.

This is commonly known as the stochastic blockmodel. It was introduced by Holland, Laskey and Leinhardt hll83 as a simple stochastic model of social networks. It has become one of the most successful and widely used models for community structure in networks, especially after the advent of large data sets.

Early analysis of the stochastic blockmodel was carried out by Snijders and Nowicki sn97 , sn01 , who provided consistent parameter estimates when there are exactly two blocks. This was extended to a finite but fixed number of blocks of equal size by Condon and Karp condonkarp01 . Bickel and Chen bickelchen09 were the first to give consistent estimates for finite number of blocks of unequal size. It was observed by Leskovec et al. leskovecetal08 that in real data, the number of blocks often seem to grow with the number of nodes. This situation was rigorously analyzed for the first time in Rohe et al. rcy , and was followed up shortly thereafter by bcl11 , cwa , mossel12 , chaudhurietal12 with more advanced results.

However, all in all, I am not aware of any estimator for the stochastic blockmodel that works whenever the number of blocks is small compared to the number of nodes. The best result till date is in the very recent manuscript of Rohe et al. roheetal12 , who prove that a penalized likelihood estimator works whenever $k$ is comparable to $n$ “up to log factors.” The following theorem says that the USVT estimator $\hat{M}$ gives a complete solution to the estimation problem in the stochastic blockmodel if $k\ll n$, with no further conditions required. (The method will not work very well for sparse graphs, however; for recent advances on estimation in sparse graphs, see chenetal12 .)

For a stochastic blockmodel with $k$ blocks,

where $C$ is a constant that depends only on our choice of $\eta$.

Note that estimating the stochastic blockmodel is a special case of low rank matrix estimation with noise. It is not difficult to prove that the estimation problem is impossible when $k$ is of the same order as $n$. We will not bother to write down a formal proof.

3 Distance matrices

Suppose that $K$ is a compact metric space with metric $d$. Let $x_{1},\ldots,x_{n}$ be arbitrary points from $K$, and let $M$ be the $n\times n$ matrix whose $(i,j)$th entry is $d(x_{i},x_{j})$. Such matrices are called “distance matrices”. Since $K$ is a compact metric space, the diameter of $K$ with respect to the metric $d$ must be finite. Scaling $d$ by a constant factor, we may assume without loss of generality that the diameter is bounded by $1$, so that the entries of $M$ are bounded by $1$ as required by Theorem 1.

Completing a distance matrix with missing entries has been a popular problem in the engineering and social sciences for a long time; see, for example, sd74 , bj95 , alfakihetal99 , biswasetal06 , singer08 , singer10 . It has become particularly relevant in engineering problems related to sensor networks. It is also an important issue in multidimensional scaling borggroenen10 . For some recent theoretical advances, see ohetal10 , jm11 .

In general, distance matrices need not be of low rank. Therefore, much of the literature on matrix estimation and completion does not apply to distance matrices. Surprisingly, Theorem 1 gives a complete solution of the distance matrix completion and estimation problem.

Suppose that $p\geq n^{-1+\varepsilon}$ for some $\varepsilon>0$. If $M$ is a distance matrix as above, then

where $c$ depends only on $\eta$, $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$, and $C(K,d,n)$ is a number depending only on $K$, $d$, $n$ and $\eta$ such that

The above theorem is not wholly satisfactory, since it does not indicate how fast $p$ can go to zero as $n\rightarrow\infty$ so that $\hat{M}$ is still consistent. To understand that, we need to know more about the structure of the space $K$. The following theorem gives a quantitative estimate.

Suppose that for each $\delta>0$, $N(\delta)$ is a number such that $K$ may be covered by $N(\delta)$ open $d$-balls of radius $\delta$. Then

where $C$ and $c$ depend only on $\eta$ and $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$.

(Note that the exponential term need not appear because the main term is bounded below by a positive constant if $p<n^{-2/3}$.) Thus, $\hat{M}$ is a consistent estimate as long as $p$ goes to zero slower than $n^{-2/3}$ as $n\rightarrow\infty$.

4 Latent space models

where $\varepsilon_{ij}$ are independent errors with zero mean, satisfying the restriction that $|x_{ij}|\leq 1$ almost surely. For example, $X$ may be the adjacency matrix of a random graph where the probability of an edge existing between vertices $i$ and $j$ is $f(\beta_{i},\beta_{j})$. This is one context where latent space models are widely used, starting with the work of Hoff, Raftery and Handcock hrh02 . A large body of work applying the latent space approach to real data has grown in the last decade. On the theoretical side, it was observed in bickelchen09 , bcl11 that the latent space model arises naturally from an exchangeability assumption due to the Aldous–Hoover theorem aldous81 , hoover82 . Note that distance matrices and stochastic blockmodels are both special cases of latent space models.

There have been various attempts to estimate parameters in the latent space models (e.g., hrh02 , hrt07 , airoldietal08 ). Almost all of these approaches rely on heuristic arguments and justification through simulations. The problem is that in addition to the vectors $\beta_{1},\ldots,\beta_{n}$, the function $f$ itself is an unknown parameter. If either $\beta_{i}$’s are known, or $f$ is known, the estimation problem is tractable. For example, when $f(x,y)$ is of the form $e^{x+y}/(1+e^{x+y})$, the problem was solved in cds . However, when both $f$ and $\beta_{i}$’s are unknown, the problem becomes seemingly intractable. In particular, there is an identifiability issue because $f(x,y)$ may be replaced by $h(x,y):=f(g(x),g(y))$ and $\beta_{i}$ by $g^{-1}(\beta_{i})$ for any invertible function $g$ without altering the model.

In view of the above discussion, it is a rather surprising consequence of Theorem 1 that it is possible to estimate the numbers $f(\beta_{i},\beta_{j})$, $i,j=1,\ldots,n$ from a single realization of the data matrix, under no additional assumptions than the stated ones.

Suppose that $p\geq n^{-1+\varepsilon}$. If $M$ is as above, then

where $c$ depends only on $\eta$, $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$, and $C(K,k,f,n)$ depends only on $K$, $k$, $f$, $n$ and $\eta$ such that

The problem with Theorem 9, just like Theorem 7 in Section 2.3, is that it does not give an explicit error bound, which makes it impossible to determine how fast $p$ can go to zero with $n$ so that consistency holds. Again, this is easy to fix by assuming smoothness properties of $f$ and applying Lemma 20. As a particular example, suppose that $f$ is Lipschitz with Lipschitz constant $L$, in the sense that

for all $x,y,x^{\prime},y^{\prime}\in K$.

where $C(K,k,L)$ is a constant depending only on $K$, $k$, $L$ and $\eta$.

5 Positive definite matrices

Assume that $m=n$ and $M$ is positive semi-definite. (In the statistical context, this is the same as saying that $M$ is a covariance matrix. When the diagonal entries are all $1$, $M$ is a correlation matrix.)

Completing positive definite matrices with missing entries has received a lot of attention in the linear algebra literature groneetal84 , johnson90 , bhatia08 , although most of the techniques are applicable only for relatively small matrices or when a sizable fraction of the entries are observed. In the engineering sciences, estimation of covariance matrices from a small subset of observed entries arises in the field of remote sensing (see candesplan10 , candesrecht09 , candestao10 for brief discussions).

The statistical matrix completion literature cited in Section 1 applies only to low rank positive definite matrices. It is therefore quite a surprise that the completion problem may be solved for any positive definite matrix whenever we get to observe a large number of entries from each row.

Suppose that $m=n$ and $M$ is positive semi-definite. Suppose that $p\geq n^{-1+\varepsilon}$. Then

where $C$ and $c$ depend only on $\eta$ and $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$.

What if $p$ is of order $1/n$ or less? The following theorem shows that it is impossible to estimate $M$ in this situation.

where $C$ is a positive universal constant.

6 Graphon estimation

A graphon is a measurable function $f$ from $^{2}$ into $$that satisfies$f(x,y)\equiv f(y,x)$. The term “graphon” was coined by Lovász and coauthors in the growing literature on limits of dense graphs borgsetal06 , borgsetal08 , borgsetal07 , lovaszszegedy06 , lovaszbook . Such functions also arise in the related study of weakly exchangeable random arrays diaconisjanson08 , austin08 , aldous81 , hoover82 . They have also appeared recently in large deviations cv , cv2 , lubetzky12 and mathematical statistics cd3 , radinyin .

In the graph limits literature, graphons arise as limits of graphs with increasing number of nodes. Conversely, graphons are often used to generate random graphs in a natural way. Take any $n$ and let $U_{1},\ldots,U_{n}$ be i.i.d. $\operatorname{Uniform}$ random variables. Construct a random undirected graph on $n$ vertices by putting an edge between vertices $i$ and $j$ with probability $f(U_{i},U_{j})$, doing this independently for all $1\leq i<j\leq n$. This procedure is sometimes called “sampling from a graphon” (see borgsetal08 , Section 4.4).

The statistical question is the following: Suppose that we have a random graph on $n$ vertices that is sampled from a graphon. Is it possible to estimate the graphon from a single realization of the graph? More precisely, is it possible to accurately estimate the numbers $f(U_{i},U_{j})$, $1\leq i<j\leq n$, from a single realization of the random graph? The question is similar to the one investigated in Section 2.3, but the difference is that here we are not allowed to assume any regularity on $f$ except measurability.

Taking things back to our usual setting, let $M$ be the matrix whose $(i,j)$th element is $f(U_{i},U_{j})$. Note that unlike our previous examples, $M$ is now random. So the definition of MSE should be modified to take expectation over $M$ as well.

where $C(f,n)$ is a constant depending only on $f$, $n$ and $\eta$, such that

Incidentally, after the first version of this paper was put up on arXiv, several papers (e.g., wolfe1 , yang ) on graphon estimation, advocating a number of different techniques and demonstrating applications in the statistical study of networks, have appeared in the literature.

7 Nonparametric Bradley–Terry model

Suppose there are $n$ teams playing against each other in a tournament. Every team plays against every other team at least once (often, exactly once). Suppose that $p_{ij}$ is the probability that team $i$ wins against team $j$ in a match between $i$ and $j$. Then $p_{ji}=1-p_{ij}$.

The Bradley–Terry model bt52 , originally proposed by Zermelo zermelo29 , assumes that $p_{ij}$ is of the form $a_{i}/(a_{i}+a_{j})$ for some unknown nonnegative numbers $a_{1},\ldots,a_{n}$. It is known how to estimate the parameters $a_{1},\ldots,a_{n}$ if we assume that the outcomes of all games are independent—which, in this case, is a reasonable assumption.

The Bradley–Terry model has found great success among practitioners. For an old survey of the literature on the model dating back to 1976, see df76 . Numerous extensions and applications have been proposed, for example, ht98 , agresti90 , raokupper67 , plackett75 , luce59 , luce77 , huangetal06 . The monographs of David david88 and Diaconis diaconis88 , Chapter 9, explain the statistical foundations of these models. More recently, several authors have proposed to perform Bayesian inference for (generalized) Bradley–Terry models adams05 , gormleymurphy08 , gormleymurphy09 , goruretal06 , guiversnelson09 , carondoucet12 .

For the basic Bradley–Terry model, it is possible to find the maximum likelihood estimate of the $a_{i}$’s using a simple iterative procedure zermelo29 , hunter04 , langeetal00 . The maximum likelihood estimate was shown to be jointly consistent for all $n$ parameters by Simons and Yao simonsyao99 .

We now generalize the Bradley–Terry model as follows. Suppose, as before, that $p_{ij}$ is the probability that team $i$ beats team $j$. Suppose that the teams have a particular ordering in terms of strength that is unknown to the observer. Assume that if team $i$ is stronger than team $j$, then $p_{ik}\geq p_{jk}$ for all $k\neq i,j$. Do not assume anything else about the $p_{ij}$’s; in particular, do not assume any formula for the $p_{ij}$’s in terms of hidden parameters. This is what we may call a “nonparametric Bradley–Terry model.” Note that the usual Bradley–Terry model is a special case of the nonparametric version.

In the nonparametric Bradley–Terry model, is it possible to estimate all the $p_{ij}$’s from a tournament where every team plays against every other exactly once? Is it possible to estimate the $p_{ij}$’s if only a randomly chosen fraction of the games are played? How small can this fraction be, so that accurate estimation is still possible? The following theorem provides some answers.

Consider the nonparametric Bradley–Terry model defined above. Let $M$ be the matrix whose $(i,j)$th entry is $p_{ij}$ if $i\neq j$ and if $i=j$. Let $X$ be the data matrix whose $(i,j)$th entry is $1$ if team $i$ won over team $j$, if team $j$ won over team $i$ and recorded as missing if team $i$ did not play versus team $j$. If team $i$ has played against team $j$ multiple times, let the $(i,j)$th entry of $X$ be the proportion of times that $i$ won over $j$. (Draws are not allowed.) Let all diagonal entries of $X$ be zero. Given $p\in$, suppose that for each $i$ and $j$, the game between $i$ and $j$ takes place with probability $p$ and does not take place with probability $1-p$, independent of other games. Let $\hat{M}$ be the estimate of $M$ based on the data matrix $X$. Then

where $C$ depends only on our choice of $\eta$. In particular, the estimation problem is solvable whenever $p\gg n^{-1/2}$.

A natural question is whether the threshold $p\gg n^{-1/2}$ is sharp. I do not know the answer to this question.

Proofs

We need to recall some background material before embarking on the proof of Theorem 1.

Let $A=(a_{ij})_{1\leq i\leq m,1\leq j\leq n}$ be an $m\times n$ real matrix with singular values $\sigma_{1},\ldots,\sigma_{k}$, where $k=\min\{m,n\}$. The following matrix norms are widely used in this proof.

The nuclear norm or the trace norm of $A$ is defined as

The Frobenius norm, also called the Hilbert–Schmidt norm, is defined as

The spectral norm or the operator norm of $A$ is defined as

The spectral norm may be alternatively expressed as

In particular, the spectral norm is a Lipschitz function of the matrix entries (with Lipschitz constant $1$), if the entries are collectively considered as a vector of length $mn$.

The triangle inequality for the spectral norm also implies that the map $A\mapsto\|A\|$ is convex. Indeed, for any $0\leq t\leq 1$,

Perturbation of singular values

The following perturbative result from matrix analysis is used several times in this manuscript. Let $A$ and $B$ be two $m\times n$ matrices. Let $k=\min\{m,n\}$. Let $\sigma_{1},\ldots,\sigma_{k}$ be the singular values of $A$ in decreasing order and repeated by multiplicities, and let $\tau_{1},\ldots,\tau_{k}$ be the singular values of $B$ in decreasing order and repeated by multiplicities. Let $\delta_{1},\ldots,\delta_{k}$ be the singular values of $A-B$, in any order but still repeated by multiplicities.

The above result follows, for example, from a combination ofTheorem III.4.4 and Exercise II.1.15 in bhatia97 . It may also be derived as a consequence of Wielandt’s minimax principle bhatia97 , Section III.3, or Lidskii’s theorem bhatia97 , Exercise III.4.3. The case $p=2$ is sometimes called the Hoffman–Wielandt theorem agz , Lemma 2.1.19 and Remark 2.1.20, and the inequality involving the maximum is sometimes called Weyl’s perturbation theorem bhatia97 , Corollary III.2.6.

Bernstein’s inequality

The following inequality is known as “Bernstein’s inequality.”

Suppose that $X_{1},\ldots,X_{n}$ are independent random variables with zero mean, and $M$ is a constant such that $|X_{i}|\leq M$ with probability one for each $i$. Let $S:=\sum_{i=1}^{n}X_{i}$ and $v:=\operatorname{Var}(S)$. Then for any $t\geq 0$,

This inequality was proved by Bernstein bernstein . For a discussion of Bernstein’s inequality and improvements, see Bennett bennett62 .

Talagrand’s concentration inequality

The following concentration inequality is one of the several striking inequalities that are collectively known as “Talagrand’s concentration inequalities.”

For a proof of Theorem 17, see talagrand96 , Theorem 6.6.

The above inequality has a number of uses in the proof of Theorem 1.

Spectral norms of random matrices

The following bound on spectral norms of random matrices is a crucial ingredient for this paper. The proof follows from a combinatorial argument of Vu vu07 (which is itself a refinement of a classical argument of Füredi and Komlós furedikomlos81 ), together with Talagrand’s inequality (4).

Take any two numbers $m$ and $n$ such that $1\leq m\leq n$. Suppose that $A=(a_{ij})_{1\leq i\leq m,1\leq j\leq n}$ is a matrix whose entries are independent random variables that satisfy, for some $\sigma^{2}\in$,

Suppose that $\sigma^{2}\geq n^{-1+\varepsilon}$ for some $\varepsilon>0$. Then for any $\eta\in(0,1)$,

where $C_{1}(\varepsilon)$ depends only on $\varepsilon$ and $\eta$ and $C_{2}$ depends only on $\eta$. The same result is true when $m=n$ and $A$ is symmetric or skew-symmetric, with independent entries on and above the diagonal, all other assumptions remaining the same. Lastly, all results remain true if the assumption $\sigma^{2}\geq n^{-1+\varepsilon}$ is changed to $\sigma^{2}\geq n^{-1}(\log n)^{6+\varepsilon}$.

First assume that $m=n$ and $A$ is symmetric. Note that for any even number $k$,

Thus, if $W(n,k,p)$ is the number of tours of length $k$ that visit exactly $p$ vertices and traverse each of its edges at least twice, then

Vu vu07 , equation (5), proves that if $p>k/2$ then $W(n,k,p)=0$ and if $p\leq k/2$ then

Using this bound, one can proceed as in vu07 , Section 2, to arrive at the conclusion that if $k$ is largest even number $\leq\sigma^{1/3}n^{1/6}$, then

This shows that if $\sigma^{2}\geq n^{-1+\varepsilon}$ [or if $\sigma^{2}\geq n^{-1}(\log n)^{6+\varepsilon}$], then there is a constant $C(\varepsilon)$ depending only on $\varepsilon$ and $\eta$ such that if $n\geq C(\varepsilon)$ then

Since $a_{ij}$ are independent and $|a_{ij}|\leq 1$ almost surely for all $i,j$, and the spectral norm is a convex Lipschitz function of matrix entries with Lipschitz constant $1$ (by the discussion about matrix norms at the beginning of this section), therefore one can apply Talagrand’s inequality [Theorem 17 and inequality (4)] together with (8) and the assumption that $\sigma^{2}\geq n^{-1+\varepsilon}$ to conclude that there is a constant $C(\varepsilon)$ such that if $n\geq C(\varepsilon)$ then

where $C_{1}$ and $C_{2}$ depend only on $\eta$. Replacing $C_{1}$ by a large enough constant $C_{1}(\varepsilon)$, the condition $n\geq C(\varepsilon)$ may be dropped. It is clear from the argument that it goes through in the skew-symmetric case as well.

Let us now drop the assumption of symmetry, but retain the assumption that $m=n$. Let $a_{ij}^{\prime}:=a_{ji}$. Then inequality (5) must be modified to say that for any even $k$,

As before, the term inside the sum is zero for any tour that traverses an edge exactly once. (In fact, there are more terms that are zero now; a term may be zero even if a tour traverses all of its edges at least twice.) Similarly, inequalities (6) and (7) continue to hold and, therefore, so does the rest of the argument.

Lastly, consider the case $m<n$. Augment the matrix $A$ by adding an extra $n-m$ rows of zeros to make it an $n\times n$ matrix that satisfies all the conditions of the theorem. Clearly, the new matrix has the same spectral norm as the old one. This completes the proof.

The key lemma

Suppose that $A$ and $B$ are two $m\times n$ matrices, where $m\leq n$. Let $a_{ij}$ be the $(i,j)$th entry of $A$ and $b_{ij}$ be the $(i,j)$th entry of $B$. It is easy to see from definition that

Thus, if $\|A-B\|$ is small enough, then the entries of $A$ are approximately equal to the entries of $B$, on average. In other words, the matrix $A$ is an estimate of the matrix $B$.

The goal of this section is to show that if in addition to the smallness of $\|A-B\|$, we also know that the nuclear norm $\|B\|_{*}$ is not too large, it is possible to get a better estimate of $B$ based on $A$.

Let $A=\sum_{i=1}^{m}\sigma_{i}x_{i}y_{i}^{T}$ be the singular value decomposition of $A$. Fix any $\delta>0$ and define

where $K(\delta)=(4+2\delta)\sqrt{2/\delta}+\sqrt{2+\delta}$.

Let $B=\sum_{i=1}^{m}\tau_{i}u_{i}v_{i}^{T}$ be the singular value decomposition of $B$. Without loss of generality, assume that $\sigma_{i}$’s and $\tau_{i}$’s are arranged in decreasing order. Let $S$ be the set of $i$ such that $\sigma_{i}>(1+\delta)\|A-B\|$. Define

Note that by the definition of $\hat{B}$, the largest singular value of $A-\hat{B}$ is bounded above by $(1+\delta)\|A-B\|$. In other words,

Since $\hat{B}$ and $G$ both have rank $\leq|S|$, the difference $\hat{B}-G$ has rank at most $2|S|$. Using this and (14), we have

Combining (LABEL:abmain) and (18), the proof is complete.

Finishing the proof of Theorem 1

We will prove the theorem only for the asymmetric model. The only difference in the proofs for the symmetric model and the skew-symmetric model is that we need to use the symmetric and skew-symmetric parts of Theorem 18 instead of the asymmetric part.

Throughout this proof, $C(\varepsilon)$ will denote any constant that depends only on $\varepsilon$ and $\eta$, and $C$ and $c$ will denote constants that depend only on $\eta$. The values of $C(\varepsilon)$, $C$ and $c$ may change from line to line or even within a line. We will use the fact that $\eta\in(0,1)$ without mention on many occasions.

Let $\hat{p}$ be the proportion of observed entries. Define two events $E_{1}$ and $E_{2}$ as

By Bernstein’s inequality (Theorem 16), for any $t\geq 0$,

Let $K(\delta)$ be the constant in the statement of Lemma 19. It is easy to see that there is a constant $C$ depending only on $\eta$ such that if $\delta\geq\eta/5$, then $K(\delta)\leq C\sqrt{1+\delta}$. Therefore, by Lemma 19, if $E_{1}$ and $E_{2}$ both happen, then

By the definition of $\hat{M}$, it is obvious that $|\hat{m}_{ij}-m_{ij}|\leq|w_{ij}-m_{ij}|$ for all $i$ and $j$. Together with (3.1), this shows that under $E_{1}\cap E_{2}$,

Dividing throughout by $mn$, we arrive at the inequality

First, suppose that $\|M\|_{*}>\eta\sqrt{n/p}/20$. Then

and so (24) follows from (23). Therefore, assume that $\|M\|_{*}\leq\eta\sqrt{n/p}/20$. Then in particular, $\|M\|\leq\eta\sqrt{n/p}/20$. Therefore, if $E_{1}\cap E_{2}$ happens, then

This implies that there is no singular value of $Y$ that exceeds $(2+\eta)\sqrt{n\hat{p}}$, and therefore $\hat{M}=0$. Consequently,

Thus, if $\|M\|_{*}\leq\eta\sqrt{n/p}/20$, then by (20) and (21),

Combining the above steps and observing that $\operatorname{MSE}(\hat{M})\leq 1$ due to the boundedness of the entries of $M$ and $\hat{M}$, we get

To remove the $1/np$ term, note that if that term indeed matters, then we are in a situation where

But this inequality, on the other hand, implies that

Therefore, the $1/np$ term can be removed from the above bound. This completes the proof of Theorem 1 if no nontrivial bound on $\operatorname{Var}(x_{ij})$ is known.

If $\sigma^{2}\leq 1$ is a known constant such that $\operatorname{Var}(x_{ij})\leq\sigma^{2}$ for all $i,j$, then the estimate (19) may be improved to

The maximum must be attained at one of the four vertices of $R$. An easy verification shows that the maximum is always attained at the vertex $(1-\sigma^{2},\sigma^{2})$, which gives the upper bound

If $E_{1}\cap E_{2}\cap E_{3}$ happens, then the subsequent steps remain the same, but with some suitable modifications that replace the term $\|M\|_{*}/(m\sqrt{np})$ by the improved term $\|M\|_{*}\sqrt{q}/(m\sqrt{n}p)$.

2 Proof of Theorem 2 (Minimax optimality)

Throughout this proof, $C$ will denote any positive universal constant, whose value may change from line to line.

Take any $\delta\in[0,m\sqrt{n}]$ and let $\theta:=\delta/(m\sqrt{n})$. We will first work out the proof under the assumption that $p<1/2$. Under this assumption, three situations are considered. First, suppose that

Let $k:=[m\theta\sqrt{p}]$. Clearly, $k\leq m$. Let $M$ be an $m\times n$ random matrix whose first $k$ rows consist of i.i.d. $\operatorname{Uniform}$ random variables, and copy this block $[1/p]$ times. This takes care of $k[1/p]$ rows. [This is okay, since $k/p\leq m\theta/\sqrt{p}\leq m$ by (25).] Declare the remaining rows, if any, to be zero. Then note that $M$ has rank $\leq k\leq m\theta\sqrt{p}$. Therefore, by inequality (2),

Let $X=M$. Let $D$ be our data, that is, the observed values of $X$. One can imagine $D$ as a matrix whose $(i,j)$th entry is $x_{ij}$ if $x_{ij}$ is observed, and a question mark if $x_{ij}$ is unobserved. For any $(i,j)$ belonging to the nonzero portion of the matrix $M$, $M$ contains $[1/p]$ copies of $m_{ij}$. Since the $X$-value at the location of each copy is observed with probability $p$, independent of the other copies, and $p<1/2$, therefore, the chance that none of these copies are observed is bounded below by a positive universal constant. If none of the copies are observed, then the data contains no information about $m_{ij}$. Using this, it is not difficult to write down a formal argument that shows

Combining the last two displays, we see that

The argument that led to the above lower bound is a typical example of the classical Bayesian argument for obtaining minimax lower bounds, and will henceforth be referred to as the “standard minimax argument” to avoid repetition of details.

Let $M$ be an $m\times n$ matrix whose first row consists of i.i.d. random variables uniformly distributed over the interval $[-m\theta\sqrt{p},m\theta\sqrt{p}]$, and this row is copied $[1/p]$ times, and all other rows are zero. Then $M$ has rank $\leq 1$, and therefore by inequality (2),

In particular, under (27), there exists $M$ with $\|M\|_{*}\leq\delta$ such that

Let $M$ be an $m\times n$ matrix whose first $[mp]$ rows consist of i.i.d. random variables uniformly distributed over $$, and this block is copied$[1/p]$times. Then the rank of$M$is$\leq[mp]$, and so by (28) and (2),

This complete the proof of Theorem 2 for the asymmetric model. For the symmetric model, simply observe that the singular values of any square matrix $M$ are the same as those of the symmetric matrix

with multiplicity doubled. It is now clear how the minimax arguments for the asymmetric model may be carried over to the symmetric case by considering the same Bayesian models for $M$ and working with the corresponding symmetrized matrices. For the skew-symmetric case, replace the $M^{T}$ by $-M^{T}$ in the above matrix.

3 Proof of Theorem 3 (Impossibility of error estimation)

Suppose, without loss of generality, that all the data matrices are defined on the same probability space. Then taking a subsequence if necessary, we may assume that in addition to (29) and (30), we also have

where the last step follows from the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$ and the triangle inequality for the Frobenius norm. Taking expectation on both sides gives

In particular, since mean squared errors are uniformly bounded by $1$,

4 Proof of Theorem 4 (Upper bound for low rank matrix estimation)

5 Proof of Theorem 5 (Lower bound for low rank matrix estimation)

Let $M$ be an $m\times n$ random matrix whose first $r$ rows consist of i.i.d. $\operatorname{Uniform}$ random variables, and copy this block $[m/r]$ times. Declare the remaining rows, if any, to be zero. Then note that $M$ has rank${}\leq r$.

Let $D$ be our data, that is, the observed values of $M$. One can imagine $D$ as a matrix whose $(i,j)$th entry is $m_{ij}$ if $m_{ij}$ is observed, and a question mark if $m_{ij}$ is unobserved. For any $(i,j)$ belonging to the nonzero portion of the matrix $M$, $M$ contains $[m/r]$ copies of $m_{ij}$. Since the $M$-value at the location of each copy is observed with probability $p$, independent of the other copies, the chance that none of these copies are observed is equal to $(1-p)^{[m/r]}$. If none of the copies are observed, then the data contains no information about $m_{ij}$. Using this, it is not difficult to write down a formal argument that shows

Combining the last two displays, we see that

6 Proof of Theorem 6 (Block model estimation)

If two vertices $i$ and $j$ are in the same block, then the $i$th and $j$th rows of $M$ are identical. Therefore, $M$ has at most $k$ distinct rows and so the rank of $M$ is $\leq k$. An application of Theorem 4 completes the proof.

7 Proofs of Theorems 7 and 8 (Distance matrix estimation)

The proofs of Theorems 7 and 8 follow from a more general lemma that will also be useful later for other purposes. Suppose that $S=\{x_{1},\ldots,x_{n}\}$ is a finite set and $f\dvtx S\times S\rightarrow$ is an arbitrary function. Suppose that for each $\delta>0$, there exists a partition $\mathcal{P}(\delta)$ of $S$ such that whenever $x,y,x^{\prime},y^{\prime}$ are four points in $S$ such that $x,x^{\prime}\in P$ for some $P\in\mathcal{P}(\delta)$ and $y,y^{\prime}\in Q$ for some $Q\in\mathcal{P}(\delta)$, then $|f(x,y)-f(x^{\prime},y^{\prime})|\leq\delta$. Let $M$ be the $n\times n$ matrix whose $(i,j)$th element is $f(x_{i},x_{j})$.

where $C$ and $c$ depend only on $\eta$, and $C(\varepsilon)$ depends only on $\varepsilon$ and $\eta$.

Fix some $\delta>0$. Let $T$ be a subset of $S$ consisting of exactly one point from each member of $\mathcal{P}(\delta)$. For each $x\in S$, let $p(x)$ be the unique element of $T$ such that $x$ and $p(x)$ belong to the same element of $\mathcal{P}(\delta)$. Let $N$ be the matrix whose $(i,j)$th element is $f(p(x_{i}),p(x_{j}))$. Then

By the triangle inequality for the nuclear norm, the inequality (2) and the above inequality,

Now, if $x_{i}$ and $x_{j}$ belong to the same element of $\mathcal{P}(\delta)$, then $p(x_{i})=p(x_{j})$, and hence the $i$th and $j$th rows of $N$ are identical. This shows that $N$ has at most $|\mathcal{P}(\delta)|$ distinct rows and, therefore, has rank${}\leq|\mathcal{P}(\delta)|$. Therefore, by the inequality (2),

The proof is completed by applying Theorem 1.

Using Lemma 20, it is easy to prove Theorems 7 and 8. {pf*}Proof of Theorem 8 Let all notation be as in Theorem 8. To apply Lemma 20, let $S$ be the set $\{x_{1},\ldots,x_{n}\}$. From the definition of $N(\delta)$, it is easy to see that there is a partition $\mathcal{P}(\delta)$ of $S$ of size $\leq N(\delta/4)$, such that any two points belonging to the same element of the partition are at distance $\leq\delta/2$ from each other. Consequently, if $x,x^{\prime}\in P$ and $y,y^{\prime}\in Q$ for some $P,Q\in\mathcal{P}(\delta)$, then by the triangle inequality for the metric $d$,

Putting $f=d$ in Lemma 20, the proof is complete.

Proof of Theorem 7 Since $K$ is compact, there exists a finite number $N(\delta)$ for each $\delta>0$ such that $K$ may be covered by $N(\delta)$ open $d$-balls of radius $\delta$. By Theorem 8, this shows that for any sequence $\delta_{n}$ decreasing to zero,

To complete the proof, choose $\delta_{n}$ going to zero so slowly that $N(\delta_{n}/4)=o(n)$ as $n\rightarrow\infty$.

8 Proof of Theorem 9 (Latent space models: General case)

We will apply Lemma 20. Let $S$ be the set $\{\beta_{1},\ldots,\beta_{n}\}$. Since $f$ is continuous on $K$ and $K$ is compact, $f$ must be uniformly continuous. This shows that for each $\delta>0$ we can find a partition $\mathcal{P}(\delta)$ of $S$ satisfying the condition required for Lemma 20, such that the size of $\mathcal{P}(\delta)$ may be bounded by a constant $N(\delta)$ depending only on $K$, $k$, $f$ and $\delta$. Choosing $\delta_{n}\rightarrow 0$ slowly enough so that $N(\delta_{n}/4)=o(n)$ and applying Lemma 20 completes the proof.