Optimal Link Prediction
with Matrix Logistic Regression
Nicolai Baldin∗ and Quentin Berthet†

arXiv:1803.07054v1 [math.ST] 19 Mar 2018

University of Cambridge

Abstract. We consider the problem of link prediction, based on partial observation of a large network, and on side information associated
to its vertices. The generative model is formulated as a matrix logistic regression. The performance of the model is analysed in a highdimensional regime under a structural assumption. The minimax rate
for the Frobenius-norm risk is established and a combinatorial estimator based on the penalised maximum likelihood approach is shown to
achieve it. Furthermore, it is shown that this rate cannot be attained
by any (randomised) algorithm computable in polynomial time under
a computational complexity assumption.
Primary 62J02, 62C20; secondary 68Q17.
Key words and phrases: Link prediction, Logistic regression, Computational lower bounds, Planted clique.
INTRODUCTION
In the field of network analysis, the task of link prediction consists in predicting
the presence or absence of edges in a large graph, based on the observations of
some of its edges, and on side information. Network analysis has become a growing
inspiration for statistical problems. Indeed, one of the main characteristics of
datasets in the modern scientific landscape is not only their growing size, but
also their increasing complexity. Most phenomena now studied in the natural
and social sciences concern not only isolated and independent variables, but also
their interactions and connections.
The fundamental problem of link prediction is therefore naturally linked with
statistical estimation: the objective is to understand, through a generative model,
why different vertices are connected or not, and to generalise these observations
to the rest of the graph.
Most statistical problems based on graphs are unsupervised: the graph itself
is the sole data, there is no side information, and the objective is to recover
an unknown structure in the generative model. Examples include the planted
clique problem (Alon and Sudakov, 1998; Kučera, 1995), the stochastic block
model (Holland et al., 1983)—see Abbe (2017) for a recent survey of a very active line of work (Abbe and Sandon, 2015; Banks et al., 2016; Decelle et al., 2011;
∗

Supported by the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement No 647812)
†
Supported by an Isaac Newton Trust Early Career Support Scheme and by The Alan Turing
Institute under the EPSRC grant EP/N510129/1.

1

2

BALDIN AND BERTHET

Massoulié, 2014; Mossel et al., 2013, 2015), the Ising blockmodel (Berthet et al.,
2016), random geometric graphs – see Penrose (2003) for an introduction and
(Bubeck et al., 2014; Devroye et al., 2011) for recent developments in statistics,
or metric-based learning (Bellet et al., 2014; Chen et al., 2009) and ordinal embedings (Jain et al., 2016).
In supervised regression problems on the other hand, the focus is on understanding a fundamental mechanism, formalized as the link between two variables.
The objective is to learn how an explanatory variable X allows to predict a response Y , i.e. to find the unknown function f that best approximates the relationship Y ≈ f (X). This statistical framework is often applied to the observation
of a phenomenon measured by Y (e.g. of a natural or social nature), given known
information X: the principle is to understand said phenomenon, to explain the
relationship between the variables by estimating the function f (Hoff et al., 2002;
Holland and Leinhardt, 1981).
We follow this approach here: our goal is to learn how known characteristics
of each agent (represented by a node) in the network induce a greater or smaller
chance of connection, to understand the mechanism of formation of the graph.
We propose a model for supervised link prediction, using the principle of regression for inference on graphs. For each vertex, we are given side information, a
vector of observations X ∈ Rd . Given observations Xi , Xj about nodes i and j of
a network, we aim to understand how these two explanatory variables are related
to the probability of connection between the two corresponding vertices, such
that P(Y(i,j) = 1) = f (Xi , Xj ), by estimating f within a high-dimensional class
based on logistic regression. Besides this high-dimensional parametric modelling,
various fully non-parametric statistical frameworks were exploited in the literature, see, for example, Gao et al. (2015); Wolfe and Olhede (2013), for graphon
estimation, Biau and Bleakly (2008); Papa et al. (2016) for graph reconstruction
and Bickel and Chen (2009) for modularity analysis.
Link prediction can be useful in any application where data can be gathered
about the nodes of a network. One of the most obvious motivations is in social
networks, in order to model social interactions. With access to side information
about each member of a social network, the objective is to understand the mechanisms of connection between members: shared interests, differences in artistic
tastes or political opinion (Wasserman and Faust, 1994). This can also be applied
to citation networks, or in the natural sciences to biological networks of interactions between molecules or proteins (Madeira and Oliveira, 2004; Yu et al., 2008).
The key assumption in this model is that the network is a consequence of the
information, but not necessarily based on similarity: it is possible to model more
complex interactions, e.g. where opposites attract.
The focus on a high-dimensional setting is another aspect of this work that is
also motivated by modern applications of statistics: data is often collected without
discernment and the ambient dimension d can be much larger than the sample
size. This setting is common in regression problems: the underlying model is often
actually very simple, to reflect the fact that only a small number of measured
parameters are relevant to the problem at hand, and that the intrinsic dimension
is much smaller. This is usually handled through an assumption on the rank,
sparsity, or regularity of a parameter. Here this needs to be adapted to a model
with two covariates (explanatory variables) and a structural assumption is made

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

3

in order to reflect this nature of our problem.
We therefore decide to tackle link prediction by modelling it as matrix logistic
regression. We study a generative model for which P(Y(i,j) = 1) = σ(Xi⊤ Θ⋆ Xj ),
where σ is the sigmoid function, and Θ⋆ is the unknown matrix to estimate.
It is a simple way to model how the variables interact, by a quadratic affinity
function and a sigmoid function. In order to model realistic situations with partial
observations, we assume that Y(i,j) is only observed for a subset of all the couples
(i, j), denoted by Ω.
To convey the general idea of a simple dependency on Xi and Xj , we make
structural assumptions on the rank and sparsity of Θ⋆ . This reflects that the
affinity Xi⊤ Θ⋆ Xj is a function of the projections u⊤
ℓ X for the vectors Xi and Xj ,
for a small number of orthogonal vectors, that have themselves a small number
of non-zero coefficients (sparsity assumption). In order to impose that the inverse
problem is well-posed, we also make a restricted conditioning assumption on Θ⋆ ,
inspired by the restricted isometry property (RIP). These conditions are discussed
in Section 1. We talk of link prediction as this is the legacy name but we focus
almost entirely on the problem of estimating Θ⋆ .
The classical techniques of likelihood maximization can lead to computationally intractable optimization problems. We show that in this problem as well as
others this is a fundamental difficulty, not a weakness of one particular estimation
technique; statistical and computational complexities are intertwined.
Our contribution: This work is organized in the following manner: We give
a formal description of the problem in Section 1, as well as a discussion of our
assumptions and links with related work. Section 2 collects our main statistical
results. We propose an estimator Θ̂ based on the penalised maximum likelihood
approach and analyse its performance in Section 2.1 in terms of non-asymptotic
rate of estimation. We show that it attains the minimax rate of estimation over
simultaneously block-sparse and low-rank matrices Θ⋆ , but is not computationally
tractable. In Section 2.2, we provide a convex relaxation of the problem which is
in essence the Lasso estimator applied to a vectorised version of the problem. The
link prediction task is covered in Section 2.3. A matching minimax lower bound for
the rate of estimation is given in Section 2.4. Furthermore, we show in Section 3
that the minimax rate cannot be attained by a (randomised) polynomial-time
algorithm, and we identify a corresponding computational lower bound. The proof
of this bound is based on a reduction scheme from the so-called dense subgraph
detection problem. Technical proofs are deferred to the appendix. Our findings
are depicted in Figure 1.
Notation: For any positive integer n, we denote by [n] the set {1, . . . , n} and
by [[n]] the set of couples of [n], of cardinality n2 . We denote by R the set of
real numbers and by Sd the set of real symmetric matrices of size d. For a matrix
A ∈ Sd , we denote by kAkF its Frobenius norm, defined by
X
A2ij .
kAk2F =
i,j∈[d]

We extend this definition for B ∈ Sn and any subset Ω ⊆ [[n]] to its semi-norm

4

BALDIN AND BERTHET
computationally hard
minimax
kr
N

+

k
N

Lasso
k2
N

log( de
)
k

k2
N

log(d)

Figure 1: The computational and statistical boundaries for estimation and prediction in the matrix logistic regression model. Here k denotes the sparsity of Θ⋆
and r its rank, while N is the number of observed edges in the network.

kBkF,Ω defined by

kBk2F,Ω =

X

2
Bij
.

i,j : (i,j)∈Ω

The corresponding bilinear form playing the role of inner-product of two matrices B1 , B2 ∈ Sn is denoted as hhB1 , B2 iiF,Ω . For a matrix B ∈ Sn , we also
make use of the following matrix norms and pseudo-norms for p, q ∈ [0, ∞),
with kBkp,q = (kB1∗ kp · · · kBd∗ kp ) q , where Bi∗ denotes the ith row of B, and
kBk∞ = max(i,j)∈[d] |Bij |. For asymptotic bounds, we shall write f (x) . g(x) if
f (x) is bounded by a constant multiple of g(x).
1. PROBLEM DESCRIPTION
1.1 Generative model
For a set of vertices V = [n] and explanatory variables Xi ∈ Rd associated to
each i ∈ V , a random graph G = (V, E) is generated by the following model. For
all i, j ∈ V , variables Xi , Xj ∈ Rd and an unknown matrix Θ⋆ ∈ Sd , an edge
connects the two vertices i and j independently of the others according to the
distribution
(1.1)


P (i, j) ∈ E = σ(Xi⊤ Θ⋆ Xj ) =

1
.
1 + exp(−Xi⊤ Θ⋆ Xj )

Here we denote by σ the sigmoid, or logistic function.
Definition 1. We denote by πij : Sd → [0, 1] the function mapping a matrix
Θ ∈ Sd to the probability in (1.1). Let Σ ∈ Sn with Σij = Xi⊤ ΘXj denote the
so-called affinity matrix. In particular, we then have πij (Θ) = σ(Σij ).
Our observation consists of the explanatory variables Xi and of the observation
of a subset of the graph. Formally, for a subset Ω ⊆ [[n]], we observe an adjacency
vector Y indexed by Ω that satisfies, for all (i, j) ∈ Ω, Y(i,j) = 1 if and only if
(i, j) ∈ E (and 0 otherwise). We thus have

(1.2)
Y(i,j) ∼ Bernoulli πij (Θ⋆ ) , (i, j) ∈ Ω .
The joint data distribution is denoted by PΘ⋆ and is thus completely specified
by πij (Θ⋆ ), (i, j) ∈ Ω. For ease of notation, we write N = |Ω|, representing the
effective sample size. Our objective is to estimate the parameter matrix Θ⋆ , based
on the observations Y ∈ RN and on known explanatory variables X ∈ Rd×n .

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

5

This problem can be reformulated as a classical logistic regression problem.
2
Indeed, writing vec(A) ∈ Rd for the vectorized form of a matrix A ∈ Sd , we
have that
(1.3)

Xi⊤ Θ⋆ Xj = Tr(Xj Xi⊤ Θ⋆ ) = hvec(Xj Xi⊤ ), vec(Θ⋆ )i .

The vector of observation Y ∈ RN therefore follows a logistic distribution with
2
explanatory design matrix DΩ ∈ RN ×d such that DΩ (i,j) = vec(Xj Xi⊤ ) and
2
predictor vec(Θ⋆ ) ∈ Rd . We focus on the matrix formulation of this problem,
and consider directly matrix logistic regression in order to simplify the notation
of the explanatory variables and our model assumptions on Θ⋆ , that are specific
to matrices.
1.2 Comparaison with other models
This model can be compared to other settings in the statistical and learning
literature.
Generalised linear model. As discussed above in the remark to (1.3), this is
an example of a logistic regression model. We focus in this work on the case where
the matrix Θ⋆ is block-sparse. The problem of sparse generalised linear models,
and sparse logistic regression in particular has been extensively studied (see, e.g.
Abramovich and Grinshtein, 2016; Bach, 2010; Bunea, 2008; Meier et al., 2008;
Rigollet, 2012; van de Geer, 2008, and references therein). Our work focuses on
the more restricted case of block-sparse and low-rank parameter, establishing
interesting statistical and computational phenomena in this setting.
Graphon model. The graphon model is a model of a random graph in which
the explanatory variables associated with the vertices in the graph are unknown.
It has recently become popular in the statistical community, see Gao et al. (2015);
Klopp and Tsybakov (2015); Wolfe and Olhede (2013); Zhang et al. (2015). Typically, an objective of statistical inference is a link function which belongs to either a parametric or non-parametric class of functions. Interestingly, the minimax
lower bound for the classes of Hölder-continuous functions, obtained in Gao et al.
(2015), has not been attained by any polynomial-time algorithm.
Trace regression models. The modelling assumption (1.1) of the present
paper is in fact very close to the trace regression model, as it follows from
the representation (1.3). Thus, the block-sparsity and low-rank structures are
preserved and can well be studied by the means of techniques developed for
the trace regression. We refer to Fan et al. (2016); Koltchinskii et al. (2011);
Negahban and Wainwright (2011); Rohde and Tsybakov (2011) for recent developments in the linear trace regression model, and Fan et al. (2017) for the generalised trace regression model. However, computational lower bounds have not
been studied either and many existing minimax optimal estimators cannot to be
computed in polynomial time.
Metric learning. In the task of metric learning, observations depend on an
unknown geometric representation V1 , . . . , Vn of the variables in a Euclidean space
of low dimension. The goal is to estimate this representation (up to a rigid transformation), based on noisy observations of hVi , Vj i in the form of random evaluations of similarity. Formally, our framework also recovers the task of metric
learning by taking Xi = ei and Θ⋆ an unknown semidefinite positive matrix of

6

BALDIN AND BERTHET

small rank (here V ⊤ V ), since
⊤
hVi , Vj i = hV ei , V ej i = e⊤
i V V ej .

We refer to (Bellet et al., 2014; Chen et al., 2009) and references therein for a
comprehensive survey of metric learning methods.
1.3 Parameter space
The unknown predictor matrix Θ⋆ describes the relationship between the observed features Xi and the probabilities of connection πij (Θ⋆ ) = σ(Xi⊤ Θ⋆ Xj )
following Definition 1. We focus on the high-dimensional setting where d2 ≫ N :
the number of features for each vertex in the graph, and number of free parameters, is much greater than the total number of observations. In order to
counter the curse of dimensionality, we make the assumption that the function
(Xi , Xj ) 7→ πij depends only on a small subset S of size k of all the coefficients of
the explanatory variables. This translates to a block-sparsity assumption on Θ⋆ :
the coefficients Θ⋆ ij are only nonzero for i and j in S. Furthermore, we assume
that the rank of the matrix Θ⋆ can be smaller than the size of the block. Formally,
we define the following parameter spaces
n
o
d
Pk,r (M ) = Θ ∈ S : kΘk1,1 < M , kΘk0,0 ≤ k , and rank(Θ) ≤ r ,
for the coefficient-wise ℓ1 norm k · k1,1 on Sd and integers k, r ∈ [d]. We also
denote P(M ) = Pd,d (M ) for convenience.
Remark 2. The bounds on block-sparsity and rank in our parameter space
are structural bounds: we consider the case where the matrix Θ⋆ can be concisely
described in terms of the number of parameters. This is motivated by considering
the spectral decomposition of the real symmetric matrix Θ⋆ as
Θ⋆ =

r
X

λℓ u ℓ u ⊤
ℓ .

ℓ=1

The affinity Σij = Xi⊤ Θ⋆ Xj between vertices i and j is therefore only a function
of the projections of Xi and Xj along the axes uℓ , i.e.
Σij = Xi⊤ Θ⋆ Xj =

r
X

⊤
λℓ (u⊤
ℓ Xi )(uℓ Xj ) .

ℓ=1

Assuming that there are only a few of these directions uℓ with non-zero impact
on the affinity motivates the low-rank assumption, while assuming that there are
only few relevant coefficients of Xi , Xj that influence the affinity corresponds to a
sparsity assumption on the uℓ , or block sparsity of Θ⋆ . The effect of these projections on the affinity is weighted by the λℓ . By allowing for negative eigenvalues,
we allow our model to go beyond a geometric description, where close or similar
Xs are more likely to be connected. This can be used to model interactions where
opposites attract.
The assumption of simultaneously sparse and low-rank matrices arises naturally in many applications in statistics and machine learning and has attracted

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

7

considerable recent attention. Various regularisation techniques have been developed for estimation, variable and rank selection in multivariate regression
problems (see, e.g. Bunea et al., 2012; Richard et al., 2012, and the references
therein).
1.4 Explanatory variables
As mentioned above, this problem is different from tasks such as metric learning, where the objective is to estimate the Xi with no side information. Here they
are seen as covariates, allowing us to infer from the observation on the graph the
predictor variable Θ⋆ . For this task to be even possible in a high-dimensional
setting, we settle the identifiability issue by making the following variant of a
classical assumption on X ∈ Rd×n .
Definition 3 (Block isometry property). For a matrix X ∈ Rd×n and an
integer s ∈ [d], we define ∆Ω,s (X) ∈ (0, 1) as the smallest positive real such that


N 1 − ∆Ω,s (X) kBk2F ≤ kX⊤ B Xk2F,Ω ≤ N 1 + ∆Ω,s (X) kBk2F ,

for all matrices B ∈ Sd that satisfy the block-sparsity assumption kBk0,0 ≤ s.

Definition 4 (Restriced isometry properties). For a matrix A ∈ Rn×p and
an integer s ∈ [p], δs (A) ∈ (0, 1) is the smallest positive real such that


n 1 − δs (A) kvk22 ≤ kAvk22 ≤ n 1 + δs (A) kvk22 ,

for all s-sparse vectors, i.e. satisfying kvk0 ≤ s.

When p = d2 is a square, we define δB,s (A) as the smallest positive real such
that


n 1 − δB,s (A) kvk22 ≤ kAvk22 ≤ n 1 + δB,s (A) kvk22 ,
for all vectors such that v = vec(B), where B satisfies the block-sparsity assumption kBk0,0 ≤ s.

The first definition is due to Candes and Tao (2005), with restriction to sparse
vectors. It can be extended in general, as here, to other types of restrictions (see,
e.g. Traonmilin and Gribonval, 2015). Since the restriction on the vectors in the
second definition (s-by-s block-sparsity) is more restricting than in the first one
(sparsity), δB,s is smaller than δs2 . These different measures of restricted isometry
are related, as shown in the following lemma
2

Lemma 5. For a matrix X ∈ Rd×n , let DΩ ∈ RN ×d be defined row-wise by
DΩ (i,j) = vec(Xj Xi⊤ ) for all (i, j) ∈ Ω. It holds that
∆Ω,s (X) = δB,s (DΩ ) .
Proof. This is a direct consequence of the definition of DΩ , which yields
kX⊤ B Xk2F,Ω = kDΩ vec(B)k22 , and kvec(B)k22 = kBk2F .

8

BALDIN AND BERTHET

The assumptions above guarantee that the matrix Θ⋆ can be recovered from
observations of the affinities, settling the well-posedness of this part of the inverse
problem. However, we do not directly observe these affinities, but their image
through the sigmoid function. We must therefore further impose the following
assumption on the design matrix X that yields constraints on the probabilities
πij and in essence governs the identifiability of Θ⋆ .
Assumption 6. There exists a constant M > 0 such that for all Θ in the
class P(M ) we have max(i,j)∈Ω |Xi⊤ ΘXj | < M .
In particular, under this assumption a constant
(1.4)


L(M ) := σ ′ (M ) = σ(M ) 1 − σ(M ) ,

is lower bounded away from zero, and we have
(1.5)

inf

Θ∈P(M )

σ ′ (Xi⊤ ΘXj ) ≥ L(M ) > 0 ,

for all (i, j) ∈ Ω. Assuming that L(M ) always depends on the same M , we
sometimes write simply L.
Remark 7. Assumption 6 is necessary for the identifiability of Θ⋆ : if Xi⊤ Θ⋆ Xj
can be arbitrarily large in magnitude, πij = σ(Xi⊤ Θ⋆ Xj ) can be arbitrarily close
to 0 or 1. Since our observations only depend on Θ⋆ through its image πij , this
could lead to a very large estimation error on Θ⋆ even with a small estimation
error on the πij .
Remark 8. This assumption has already appeared in the literature on highdimensional estimation, see Abramovich and Grinshtein (2016); van de Geer (2008).
Similarly to Bach (2010), Assumption 6 can be shown to be redundant for minimax optimal prediction, because the log-likelihood function in the matrix logistic
regression model satisfies the so-called self-concordant property. Our analysis to
follow can be combined with an analysis similar to Bach (2010) to get rid of the
assumption for minimax optimal prediction.
Proposition 9. The identifiability assumption max(i,j)∈Ω |Xi⊤ ΘXj | < M is
guaranteed for all Θ ∈ P(M ) and design matrices X satisfying either of the
following
• kXj Xi⊤ k∞ ≤ 1,
• kΘk2F < M1 for some M1 > 0 and the block isometry property.
1.4.1 Random designs For random designs, we require the block isometry
property to hold with high probability. Then the results in this article carry
over directly and thus we do not discuss it in full detail. It is well known that
for sparse linear models with the dimension of a target vector p̄ and the sparsity
k̄, the classical restricted isometry property holds for some classes of random
matrices with i.i.d. entries including sub-Gaussian and Bernoulli matrices, see
Mendelson et al. (2008), provided that n̄ & k̄ log(p̄/k̄), and i.i.d. subexponential
random matrices, see Adamczak et al. (2011), provided that n̄ & k̄ log2 (p̄/k̄).

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

9

In the same spirit, the design matrices with independent entries following subGaussian, subexponential or Bernoulli distributions can be shown to satisfy the
block isometry property, cf. Wang et al. (2016a), provided that the number of
observed edges in the network satisfies N & k2 log2 (d/k) for sub-Gaussian and
subexponential designs and N & k2 log(d/k) for Bernoulli designs.
2. MATRIX LOGISTIC REGRESSION
The log-likelihood for this problem is
X
ℓY (Θ) = −
ξ(s(i,j) Xi⊤ ΘXj ) ,
(i,j)∈Ω

where s(i,j) = 2Y(i,j) − 1 is a sign variable that depends on the observations Y
and ξ : x 7→ log(1 + ex ) is a softmax function, convex on R. As a consequence,
the negative log-likelihood −ℓY is a convex function of Θ. Denoting by ℓ the
expectation IEΘ⋆ [ℓY ], we recall the classical expressions for all Θ ∈ Sd
X
ℓ(Θ) = ℓ(Θ⋆ ) −
KL(πij (Θ⋆ ), πij (Θ))
(i,j)∈Ω

= ℓ(Θ⋆ ) − KL(PΘ⋆ , PΘ ) ,
where we recall πij (Θ) = σ(Xi⊤ ΘXj ), and
ℓY (Θ) = ℓ(Θ) + hh∇ζ, ΘiiF ,
where ζ is a stochastic component
of the log-likelihood with constant gradient
P
d×d
∇ζ ∈ R
given by ∇ζ = (i,j)∈Ω (Y(i,j) − πij (Θ⋆ ))Xj Xi⊤ , which is a sum of
independent centered random variables.
2.1 Penalized logistic loss
In a classical setting where d is fixed and N grows, the maximiser of ℓY - the
maximum likelihood estimator - is an accurate estimator of Θ⋆ , provided that it
is possible to identify Θ from PΘ (i.e. if the Xi are well conditioned). We are
here in a high-dimensional setting where d2 ≫ N , and this approach is not directly possible. Our parameter space indicates that the intrinsic dimension of our
problem is truly much lower in terms of rank and block-sparsity. Our assumption
on the conditioning of the Xi is tailored to this structural assumption. In the
same spirit, we also modify our estimator in order to promote the selection of
elements of low rank and block-sparsity. Following the ideas of Birgé and Massart
(2007) and Abramovich and Grinshtein (2016), we define the following penalized
maximum likelihood estimator
n
o
(2.1)
Θ̂ ∈ argmin − ℓY (Θ) + p(Θ) ,
Θ∈P(M )

with a penalty p defined as
(2.2)

p(Θ) = g(rank(Θ), kΘk0,0 ) ,

and g(R, K) = cKR + cK log

de 
,
K

where c > 0 is a universal constant and to be specified further. The proof of
the following theorem is based on Dudley’s integral argument combined with
Bousquet’s inequality and is deferred to the Appendix.

10

BALDIN AND BERTHET

Theorem 10. Assume the design matrix X satisfies max(i,j)∈Ω |Xi⊤ Θ⋆ Xj | <
M for some M > 0 and all Θ⋆ in a given class, and the penalty term p(Θ)
satisfies (2.2) with the constants c ≥ c1 /L, c1 > 1, L given in (1.4). Then for
the penalised MLE estimator Θ̂, the following non-asymptotic upper bound on the
expectation of the Kullback-Leibler divergence between the measures PΘ⋆ and PΘ̂
holds
de 
k
kr
1
IE[KL(PΘ⋆ , PΘ̂ )] ≤ C1
+ C1 log
,
N
N
k
Θ⋆ ∈Pk,r (M ) N

(2.3)

sup

where C1 > 3c is some universal constant for all k = 1, ..., d and r = 1, ..., k.
Remark 11. Random designs with i.i.d. entries following sub-Gaussian, Bernoulli
and subexponential distributions discussed in Section 1.4.1 yield the same rate as
well. It can formally be shown using standard conditioning arguments (see, e.g.
Nickl and van de Geer, 2013).
Corollary 12. Assume the design matrix X satisfies the block isometry
property from Definition 3 and max(i,j)∈Ω |Xi⊤ Θ⋆ Xj | < M for some M > 0 and
all Θ⋆ in a given class, and the penalty term p(Θ) is as in Theorem 10. Then for
the penalised MLE estimator Θ̂, the following non-asymptotic upper bound on the
rate of estimation holds
 kr


C1
de 
k

sup
IE kΘ̂ − Θ⋆ k2F ≤
+
log
,
N
k
L(M ) 1 − ∆Ω,2k (X) N
Θ⋆ ∈Pk,r (M )
where C1 > 3c is some universal constant for all k = 1, ..., d and r = 1, ..., k.

Let us define rank-constrained maximum likelihood estimators with bounded
block size as
Θ̂k,r ∈ argmin {−ℓY (Θ)} .
Θ∈Pk,r (M )

It is intuitively clear that without imposing any regularisation on the likelihood
function, the maximum likelihood approach selects the most complex model. In
fact, the following result holds.
Theorem 13. Assume the design matrix X satisfies the block isometry property from Definition 3 and max(i,j)∈Ω |Xi⊤ Θ⋆ Xj | < M for some M > 0 and all Θ⋆
in a given class. Then for the maximum likelihood estimator Θ̂k,r , the following
non-asymptotic upper bound on the rate of estimation holds

sup
Θ⋆ ∈Pk,r (M )



IE kΘ̂k,r − Θ⋆ k2F ≤

 kr
de 
k
C3

+
log
,
N
k
L(M ) 1 − ∆Ω,2k (X) N

for all k = 1, ..., d and r = 1, ..., k and some constant C3 > 0.

Remark 14. The penalty (2.2) belongs to the class of the so-called minimal
penalties, cf. Birgé and Massart (2007). In particular, a naive MLE approach
with p(Θ) = 0 in (2.1) yields a suboptimal estimator as it follows from Theorem 13.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

11

2.2 Convex relaxation
In practice, computation of the estimator (2.1) is often infeasible. In essence,
in order to compute it, we need to compare the likelihood functions over all
possible subspaces Pk,r (M ). Sophisticated step-wise model selection procedures
allow to reduce the number of analysed models. However, they are not feasible
in a high-dimensional setting either. We here consider the following estimator
Θ̂Lasso = argmin{−ℓY (Θ) + λkΘk1,1 } ,

(2.4)

Θ∈Sd

with λ > 0 to be chosen further, which is equivalent to the logistic Lasso on
vec(Θ). Using standard arguments, cf. Example 1 in van de Geer (2008), combined with the block isometry property the following result immediately follows.
Theorem 15. Assume the design matrix X satisfies the block isometry prop⊤
erty from Definition 3 and max(i,j)∈Ω |X
√i Θ⋆ Xj | < M for some M > 0 and all
Θ⋆ in a given class. Then for λ = C4 log d, where C4 > 0 is an appropriate
universal constant, the estimator (2.4) satisfies
(2.5)

sup
Θ⋆ ∈Pk,r (M )



IE kΘ̂Lasso − Θ⋆ k2F ≤

k2
C5
 log d ,
L(M ) 1 − ∆Ω,2k (X) N

for all k = 1, ..., d and r = 1, ..., k and some universal constant C5 > 0.

As one could expect the upper bound on the rate of estimation of our feasible estimator is independent of the true rank r. It is natural, when dealing with a low-rank and block-sparse objective matrix, to combine the nuclear
penalty with either the (2, 1)-norm penalty or the (1, 1)-norm penalty of a matrix,
cf. Bunea et al. (2012); Giraud (2011); Koltchinskii et al. (2011); Richard et al.
(2012). In our setting, it can be easily shown that combining the (1, 1)-norm
penalty and the nuclear penalty yields the same rate of estimation (k2 /N ) log d.
This appears to be inevitable in view of a computational lower bound, obtained
in Section 3, which is independent of the rank as well. In particular, these findings
partially answer a question posed in Section 6.4.4 in Giraud (2014).
2.3 Prediction
In applications, as new users join the network, we are interested in predicting
the probabilities of the links between them and the existing
users. It is natural to
P
2
measure the prediction error of an estimator Θ̂ by IE
(i,j)∈Ω (πij (Θ̂)−πij (Θ⋆ ))
which is controlled according to the following result using the smoothness of the
logistic function σ.
Theorem 16. Under Assumption 6, we have the following rate for estimating
the matrix of probabilities Σ⋆ = X⊤ Θ⋆ X ∈ Rn×n with the estimator Σ̂ = X⊤ Θ̂X ∈
Rn×n :

1 
de 
k
C1  kr
sup
IE kΣ̂ − Σ⋆ k2F,Ω ≤
+
log
,
L(M ) N
N
k
Θ⋆ ∈Pk,r (M ) 2N
with the constant C1 from (2.3). The rate is minimax optimal, i.e. a minimax
lower bound of the same asymptotic order holds for the prediction error of estimating the matrix of probabilities Σ⋆ = X⊤ Θ⋆ X ∈ Rn×n .

12

BALDIN AND BERTHET

2.4 Information-theoretic lower bounds
The following result demonstrates that the minimax lower bound on the rate
of estimation matches the upper bound in Theorem 10 implying that the rate of
estimation is minimax optimal.
Theorem 17. Let the design matrix X satisfy the block isometry property.
Then for estimating Θ⋆ ∈ Pk,r (M ) in the matrix logistic regression model, the
following lower bound on the rate of estimation holds
 kr


C2
de 
k
sup
IE kΘ̂ − Θ⋆ k2F ≥
inf
+
log
,
(1 + ∆Ω,2k (X)) N
N
k
Θ̂ Θ⋆ ∈Pk,r (M )
where the constant C2 > 0 is independent of d, k, r and the infimum extends over
all estimators Θ̂.
Remark 18. The lower bounds of the same order hold for the expectation
of the Kullback-Leibler divergence between the measures PΘ⋆ and PΘ̂ and the
prediction error of estimating the matrix of probabilities Σ⋆ = X⊤ Θ⋆ X ∈ Rn×n .
3. COMPUTATIONAL LOWER BOUNDS
In this section, we investigate whether the lower bound in Theorem 17 can be
achieved with an estimator computable in polynomial time. The fastest rate of
estimation attained by a (randomised) polynomial-time algorithm in the worstcase scenario is usually referred to as a computational lower bound. Recently, the
gap between computational and statistical lower bounds has attracted a lot of attention in the statistical community. We refer to Berthet and Rigollet (2013a,b);
Chen (2015); Chen and Xu (2016); Gao et al. (2017); Hajek et al. (2015); Ma and Wu
(2015); Wang et al. (2016b); Zhang and Dong (2017) for computational lower
bounds in high-dimensional statistics based on the planted clique problem (see
below), Berthet and Ellenberg (2015) using hardness of learning parity with noise
Oymak et al. (2015) for denoising of sparse and low-rank matrices, Agarwal
(2012) for computational trade-offs in statistical learning, as well as Zhang et al.
(2014) for worst-case lower bounds for sparse estimators in linear regression, as
well as Bruer et al. (2015); Chandrasekaran and Jordan (2013) for another approach on computational trade-offs in statistical problems, as well as Berthet and Chandrasekaran
(2016); Berthet and Perchet (2017) on the management of these trade-offs. In order to establish a computational lower bound for the block-sparse matrix logistic
regression, we exploit a reduction scheme from Berthet and Rigollet (2013a): we
show that detecting a subspace of Pk,r (M ) can be computationally as hard as
solving the dense subgraph detection problem.
3.1 The dense subgraph detection problem
Although our work is related to the study of graphs, we recall for absolute
clarity the following notions from graph theory. A graph G = (V, E) is a nonempty set V of vertices, together with a set E of distinct unordered pairs {i, j}
with i, j ∈ V , i 6= j. Each element {i, j} of E is an edge and joins i to j. The
vertices of an edge are called its endpoints. We consider only undirected graphs
with neither loops nor multiple edges. A graph is called complete if every pair
of distinct vertices is connected. A graph G′ = (V ′ , E ′ ) is a subgraph of a graph

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

13

G = (V, E) if V ′ ⊆ V and E ′ ⊆ E. A subgraph C is called a clique if it is complete.
The problem of detecting a maximum clique, or all cliques, the so-called Clique
problem, in a given graph is known to be NP-complete, cf. Karp (1972).
The Planted Clique problem, motivated as an average case version of the
Clique problem, can be formalised as a decision problem over random graphs,
parametrised by the number of vertices n and the size of the subgraph k. Let
Gn denote the collection of all graphs with n vertices and G(n, 1/2) denote distribution of Erdös-Rényi random graphs, uniform on Gn , where each edge is
drawn independently at random with probability 1/2. For any k ∈ {1, ..., n}
and q ∈ (1/2, 1], let G(n, 1/2, k, q) be a distribution on Gn constructed by first
picking k vertices independently at random and connecting all edges in-between
with probability q, and then joining each remaining pair of distinct vertices by
an edge independently at random with probability 1/2. Formally, the Planted
Clique problem refers to the hypothesis testing problem of
(3.1)

H0 : A ∼ G(n, 1/2)

vs. H1 : A ∼ G(n, 1/2, k, 1) ,

based on observing an adjacency matrix A ∈ Rn×n of a random graph drawn
from either G(n, 1/2) or G(n, 1/2, k, 1).
One of the main properties of the Erdös-Rényi random graph were studied in
Erdős and Rényi (1959), as well as in Grimmett and McDiarmid (1975), who in
particular proved that the size of the largest clique in G(n, 1/2) is asymptotically
close to 2 log2 n almost surely. On the other hand, Alon et al. (1998) proposed a
√
spectral method that for k > c n detects a planted clique with high probability
in polynomial time. Hence the most intriguing regime for k is
√
(3.2)
2 log2 n ≤ k ≤ c n .
The conjecture that no polynomial-time algorithm exists for distinguishing between hypotheses in (3.1) in the regime (3.2) with probability tending to 1 as
n → ∞ is the famous Planted Clique conjecture in complexity theory. Its variations have been used extensively as computational hardness assumptions in statistical problems, see Berthet and Rigollet (2013b); Cai and Wu (2018); Gao et al.
(2017); Wang et al. (2016b).
The Planted Clique problem can be reduced to the so-called dense subgraph
detection problem of testing the null hypothesis in (3.1) against the alternative
H1 : A ∼ G(n, 1/2, k, q), where q ∈ (1/2, 1]. This is clearly a computationally
harder problem. In this paper, we assume the following variation of the Planted
Clique conjecture which is used to establish a computational lower bound in the
matrix logistic regression model.
Conjecture 19 (The dense subgraph detection conjecture). For any sequence k = kn such that k ≤ nβ for some 0 < β < 1/2, and any q ∈ (1/2, 1],
there is no (randomised) polynomial-time algorithm that can correctly identify the
dense subgraph with probability tending to 1 as n → ∞, i.e. for any sequence of
(randomised) polynomial-time tests (ψn : Gn → {0, 1})n , we have

lim inf P0 (ψn (A) = 1) + P1 (ψn (A) = 0) ≥ 1/3 .
n→∞

14

BALDIN AND BERTHET

3.2 Reduction to the dense subgraph detection problem and a
computational lower bound
Consider the vectors of explanatory variables Xi = N 1/4 ei , i = 1, ..., n and
assume without loss of generality that the observed set of edges Ω in the matrix
logistic regression model
consists of the interactions of the n nodes Xi , i.e. it

holds N = |Ω| = n2 . It follows from the matrix logistic regression modelling
assumption (1.1) that the Erdös-Rényi graph G(n, 1/2) corresponds to a random
graph associated with the matrix Θ0 = 0 ∈ Rd×d . Let Gl (k) be a subset of
Pk,1 (M ) with a fixed support l of the block. In addition, let GkαN ⊂ Pk,1 (M ) be
a subset consisting of the matrices Θl ∈ Gl (k), l = 1, ..., K, K = nk such that all
√
elements in the block of a matrix Θl equal some αN = α/ N > 0, see Figure 2.
Then we have
P (i, j) ∈ E|Xi , Xj ) =

1
−Xi⊤ ΘXj

1+e

=

1
,
1 + e−α

for all Θ ∈ GkαN . Therefore, the testing problem
(3.3)

H0 : Y ∼ PΘ0

vs. H1 : Y ∼ PΘ , Θ ∈ GkαN ,

where Y ∈ {0, 1}N is the adjacency vector of binary responses in the matrix
logistic regression model, is reduced to the dense subgraph detection problem with
q = 1/(1 + e−α ). This reduction scheme suggests that the computational lower
bound for separating the hypotheses in the dense subgraph detection problem
mimics the computational lower bound for separating the hypotheses in (3.3) in
the matrix logistic regression model. The following theorem exploits this fact in
order to establish a computational lower bound of order k2 /N for estimating the
matrix Θ⋆ ∈ Pk,r (M ).
Theorem 20. Let Fk be any class of matrices containing GkαN ∪ Θ0 from the
reduction scheme. Let c > 0 be a positive constant and f (k, d, N ) be a real-valued
function satisfying f (k, d, N ) ≤ ck 2 /N for k = kn < nβ , 0 < β < 1/2 and a
sequence d = dn , for all n > m0 ∈ IN. If Conjecture 19 holds, for some the design
X that fulfils the block isometry property from Definition 3, there is no estimator of
Θ⋆ ∈ Fk , that attains the rate f (k, d, N ) for the Frobenius norm risk, and can be
evaluated using a (randomised) polynomial-time algorithm, i.e. for any estimator
Θ̂, computable in polynomial time, there exists a sequence (k, d, N ) = (kn , dn , N ),
such that
(3.4)



1
sup IE kΘ̂ − Θ⋆ k2F → ∞ ,
f (k, d, N ) Θ⋆ ∈Fk

as n → ∞. Similarly, for any estimator Θ̂, computable in polynomial time, there
exists a sequence (k, d, N ) = (kn , dn , N ), such that
(3.5)


1
1 
sup
IE kΣ̂ − Σ⋆ k2F,Ω → ∞ ,
f (k, d, N ) Θ⋆ ∈Fk N

for the prediction error of estimating Σ⋆ = X⊤ Θ⋆ X.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

αN
GkαN

15

k
n

∋Θ=

d

Figure 2: The construction of matrices GkαN used in the reduction scheme.
Remark 21. Thus the computational lower bound for estimating the matrix
Θ⋆ in the matrix logistic regression model is of order k2 /N compared to the minimax rate of estimation of order kr/N +(k/N ) log(de/k) and the rate of estimation
(k2 /N ) log(d) for the Lasso estimator Θ̂Lasso , cf. Figure 1. Hence the computational gap is most noticeable for the matrices of rank 1. Furthermore, as a simple
consequence of this result, the corresponding computational lower bound for the
prediction risk of estimating Σ⋆ = X⊤ Θ⋆ X is k2 /N as well.
Proof. We here provide a proof of the computational lower bound on the
prediction error (3.5) for convenience. The bound on the estimation error (3.4)
is straightforward to show by utilizing the block isometry property. Assume that
there exists a hypothetical estimator Θ̂ computable in polynomial time that attains the rate f (k, d, N ) for the prediction error, i.e. such that it holds that

1  ⊤
1
sup
IE kX (Θ̂ − Θ⋆ )Xk2F,Ω ≤ b < ∞ ,
lim sup
n→∞ f (k, d, N ) Θ⋆ ∈Fk N
for all sequences (k, d, N ) = (kn , dn , N ) and a constant b. Then by Markov’s
inequality, we have
1
(3.6)
kX⊤ (Θ̂ − Θ⋆ )Xk2F,Ω ≤ uf (k, d, N ) ,
N
for some numeric constant u > 0 with probability 1 − b/u for all Θ⋆ ∈ Fk .
Following the reduction scheme, we consider the design vectors Xi = N 1/4 ei ,
i = 1, ..., n and the subset of edges Ω, such that
X
1
kX⊤ (Θ̂ − Θ⋆ )Xk2F,Ω =
(Θ̂ij − Θ⋆ij )2 = kΘ̂ − Θ⋆ k2F,Ω ,
(3.7)
N
(i,j)∈Ω

for any Θ⋆ ∈ GkαN . Note, that the design vectors Xi = N 1/4 ei , i = 1, ..., n clearly
satisfy Assumption 6. Thus, in order to separate the hypotheses
(3.8)

H0 : Y ∼ P0

vs. H1 : Y ∼ PΘ , Θ ∈ GkαN ,

it is natural to employ the following test

ψ(Y ) = 1 kΘ̂kF,Ω ≥ τd,k (u) ,

(3.9)

2 (u) = uf (k, d, N ). The type I error of this test is controlled automatiwhere τd,k
cally due to (3.6) and (3.7), P0 (ψ = 1) ≤ b/u. For the type II error, we obtain

sup PΘ (ψ = 0) = sup PΘ kΘ̂kF,Ω < τd,k (u)
α

Θ∈Gk N

α

Θ∈Gk N


2
≤ sup PΘ kΘ̂ − Θk2F,Ω > kΘk2F,Ω − τd,k
(u) ≤ b/u ,
α

Θ∈Gk N

16

BALDIN AND BERTHET

provided that
2
(u) = 2uf (k, d, N ) ,
k(k − 1)α2N /2 ≥ 2τd,k

which is true in the regime k ≤ nβ , β < 1/2, and α2 ≥ 4u/c, (hence α2N ≥
4u/(cN ) ) by the definition of the function f (k, d, N ). Putting the pieces together,
we obtain

lim sup P0 (ψ(Y ) = 1) + sup PΘ (ψ(Y ) = 0) ≤ 2b/u < 1/3 ,
n→∞

α

Θ∈Gk N

for u > 6b. Hence, the test (3.9) separates the hypotheses (3.8). This contradicts
Conjecture 19 and implies (3.5).

4. CONCLUDING REMARKS
Our results shed further light on the emerging topic of statistical and computational trade-offs in high-dimensional estimation. The matrix logistic regression
model is very natural to study the connection between statistical accuracy and
computational efficiency as the model is based on the study of a generative model
for random graphs. It is also an extension of lower bound for all statistical procedures to a model with covariates, the first of its kind.
Our fundings suggest that the block-sparsity is a limiting model selection criterion for polynomial-time estimation in the logistic regression model. That is,
imposing further structure, like an additional low-rank constraint, and thus reducing the number of studied models yields an expected gain in the minimax
rate, but that gain can never be achieved by a polynomial-time algorithm. In this
setting, this implies that with a larger parameter space, while the statistical rates
might be worse, they might be closer to those that are computationally achievable. As an illustration, both efficient and minimax optimal estimation is possible
for estimating sparse vectors in the high-dimensional linear regression model, see,
for example, SLOPE for achieving the exact minimax rate in Bellec et al. (2018);
Bogdan et al. (2015), extending upon previous results on the Danzig selector and
Lasso in Bickel et al. (2009); Candes and Tao (2007).
The logistic regression is also a representative of a large class of generalised
linear models. Furthermore, the proof of the minimax lower bound on the rate of
estimation in Theorem 17 can be extended to all generalised linear models. The
combinatorial estimator (2.1) can well be used to achieve the minimax rate. The
computational lower bound then becomes a delicate issue. A more sophisticated
reduction scheme is needed to relate the dense subgraph detection problem to
an appropriate testing problem for a generalised linear model. Approaching this
question might require notions of noise discretisation and Le Cam equivalence
studied in Ma and Wu (2015).
An interesting question is whether it is possible to adopt polynomial-time
algorithms available for detecting a dense subgraph for estimating the target
matrix in the logistic regression model in all sparsity regimes. A common idea
behind those algorithms is to search a dense subgraph over the vertices of a
high degree and thus substantially reduce the number of compared models of
subgraphs. The network we observe in the logistic regression model is generated
by a sparse matrix. We may still observe a fully connected network which is

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

17

generated by a small block in the target matrix. Therefore, it is not yet clear
how to adopt algorithms for dense subgraph detection to submatrix detection. It
remains an open question to establish whether these results can be extended to
any design matrix, and all parameter regimes.
REFERENCES
Abbe, E. (2017). Community detection and stochastic block models: recent
developments. arXiv:1703.10146 .
Abbe, E. and Sandon, C. (2015). Detection in the stochastic block model with
multiple clusters: proof of the achievability conjectures, acyclic BP, and the
information-computation gap. arXiv:1512.09080 .
Abramovich, F. and Grinshtein, V. (2016). Model selection and minimax
estimation in generalized linear models. IEEE Transactions on Information
Theory 62 3721–3730.
Adamczak, R., Litvak, A., Pajor, A. and Tomczak-Jaegermann, N.
(2011). Restricted isometry property of matrices with independent columns
and neighborly polytopes by random sampling. Constructive Approximation
34 61–88.
Agarwal,
A.
(2012).
Computational
Trade-offs in Statistical Learning.
Theses, University of California, Berkeley.
http://research.microsoft.com/en-us/um/people/alekha/thesismain.pdf.
Alon, N., Krivelevich, M. and Sudakov, B. (1998). Finding a large hidden
clique in a random graph. Random Structures and Algorithms 13 457–466.
Alon, N. K. M. and Sudakov, B. (1998). Finding a large hidden clique in a
random graph. In Proceedings of the Eighth International Conference “Random
Structures and Algorithms” (Poznan, 1997), vol. 13.
Bach, F. (2010). Self-concordant analysis for logistic regression. Electronic
Journal of Statistics 4 384–414.
Banks, J., Moore, C., Neeman, J. and Netrapalli, P. (2016).
Information-theoretic thresholds for community detection in sparse networks.
arXiv:1601.02658 .
Bellec, P., Lecué, G. and Tsybakov, A. (2018). Slope meets lasso:
improved oracle bounds and optimality. Annals of Statistics (to appear)
ArXiv:1605.08651.
Bellet, A., Habrard, A. and Sebban, M. (2014). A survey on metric learning
for feature vectors and structured data.
Berthet, Q. and Chandrasekaran, V. (2016). Resource allocation for statistical estimation. Proceedings of the IEEE 104 115–125.
Berthet, Q. and Ellenberg, J. (2015). Detection of planted solutions for flat
satisfiability problems .

18

BALDIN AND BERTHET

Berthet, Q. and Perchet, V. (2017). Fast rates for bandit optimization with
upper-confidence frank-wolfe. NIPS 2017, to appear .
Berthet, Q. and Rigollet, P. (2013a). Complexity theoretic lower bounds
for sparse principal component detection. In Conference on Learning Theory.
Berthet, Q. and Rigollet, P. (2013b). Optimal detection of sparse principal
components in high dimension. The Annals of Statistics 41 1780–1815.
Berthet, Q., Rigollet, P. and Srivastava, P. (2016). Exact recovery in the
ising blockmodel .
Biau, G. and Bleakly, K. (2008). Statistical inference on graphs. Statistics
and decisions 24 209–232.
Bickel, P. J. and Chen, A. (2009). A nonparametric view of network models and newman–girvan and other modularities. Proceedings of the National
Academy of Sciences 106 21068–21073.
Bickel, P. J., Ritov, Y. and Tsybakov, A. (2009). Simultaneous analysis of
lasso and dantzig selector. Ann. Statist. 37 1705–1732.
URL https://doi.org/10.1214/08-AOS620
Birgé, L. and Massart, P. (2007). Minimal penalties for gaussian model selection. Probability theory and related fields 138 33–73.
Bogdan, M., van den Berg, E., Sabatti, C., Su, W. and Candès, E. (2015).
Slope—adaptive variable selection via convex optimization. The Annals of
Applied Atatistics 9 1103.
Bousquet, O. (2002). A bennett concentration inequality and its application to
suprema of empirical processes. Comptes Rendus Mathematique 334 495–500.
Bruer, J. J., Tropp, J. A., Cevher, V. and Becker, S. R. (2015). Designing
statistical estimators that balance sample size, risk, and computational cost.
IEEE Journal of Selected Topics in Signal Processing 9 612–624.
Bubeck, S., Ding, J., Eldan, R. and Rácz, M. (2014). Testing for highdimensional geometry in random graphs .
URL https://arxiv.org/abs/1411.5713
Bunea, F. (2008). Honest variable selection in linear and logistic regression
models via ℓ1 and ℓ1 + ℓ2 penalization. Electronic Journal of Statistics 2 1153–
1194.
Bunea, F., She, Y. and Wegkamp, M. (2012). Joint variable and rank selection for parsimonious estimation of high-dimensional matrices. The Annals of
Statistics 40 2359–2388.
Cai, T. and Wu, Y. (2018). Statistical and computational limits for sparse
matrix detection. arXiv preprint arXiv:1801.00518 .

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

19

Candes, E. and Plan, Y. (2011). Tight oracle inequalities for low-rank matrix recovery from a minimal number of noisy random measurements. IEEE
Transactions on Information Theory 57 2342–2359.
Candes, E. and Tao, T. (2007). The dantzig selector: Statistical estimation
when p is much larger than n. Ann. Statist. 35 2313–2351.
URL https://doi.org/10.1214/009053606000001523
Candes, E. J. and Tao, T. (2005). Decoding by Linear Programming. IEEE
Trans. Information Theory 51 4203–4215.
Chandrasekaran, V. and Jordan, M. I. (2013). Computational and statistical tradeoffs via convex relaxation. Proceedings of the National Academy of
Sciences .
Chen, Y. (2015). Incoherence-optimal matrix completion. IEEE Transactions
on Information Theory 61 2909–2923.
Chen, Y., Garcia, E. K., Gupta, M. R., Rahimi, A. and Cazzanti, L.
(2009). Similarity-based classification: Concepts and algorithms. J. Mach.
Learn. Res. 10 747–776.
Chen, Y. and Xu, J. (2016). Statistical-computational tradeoffs in planted
problems and submatrix localization with a growing number of clusters and
submatrices. Journal of Machine Learning Research 17 1–57.
URL http://jmlr.org/papers/v17/14-330.html
Decelle, A., Krzakala, F., Moore, C. and Zdeborová, L. (2011). Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications. Physical Review E 84 066106.
Devroye, L., György, A., Lugosi, G. and Udina, F. (2011). Highdimensional random geometric graphs and their clique number. Electronic
Communications in Probability 16 2481–2508.
Dudley, R. (1967). The sizes of compact subsets of hilbert space and continuity
of gaussian processes. Journal of Functional Analysis 1 290 – 330.
Erdős, P. and Rényi, A. (1959). On random graphs. Publicationes Mathematicae 6 290–297.
Fan, J., Gong, W. and Zhu, Z. (2017). Generalized high-dimensional trace
regression via nuclear norm regularization. arXiv preprint arXiv:1710.08083 .
Fan, J., Wang, W. and Zhu, Z. (2016). Robust low-rank matrix recovery. arXiv
preprint arXiv:1603.08315 .
Gao, C., Lu, Y. and Zhou, H. (2015). Rate-optimal graphon estimation. The
Annals of Statistics 43 2624–2652.
Gao, C., Ma, Z. and Zhou, H. (2017). Sparse cca: Adaptive estimation and
computational barriers. arXiv preprint arXiv:1409.8565 .

20

BALDIN AND BERTHET

Giné, E. and Nickl, R. (2016). Mathematical Foundations of InfiniteDimensional Statistical Models. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press.
Giraud, C. (2011). Low rank multivariate regression. Electronic Journal of
Statistics 5 775–799.
Giraud, C. (2014). Introduction to high-dimensional statistics. Chapman &
Hall/CRC Monographs on Statistics & Applied Probability, Taylor & Francis.
Grimmett, G. and McDiarmid, C. (1975). On colouring random graphs.
In Mathematical Proceedings of the Cambridge Philosophical Society, vol. 77.
Cambridge Univ Press.
Hajek, B., Wu, Y. and Xu, J. (2015). Computational lower bounds for community detection on random graphs. In Proceedings of The 28th Conference
on Learning Theory, vol. 40 of Proceedings of Machine Learning Research.
Hoff, P. D., Raftery, A. E. and Handcock, M. S. (2002). Latent space
approaches to social network analysis. Journal of the American Statistical
Association 97 1090–1098.
Holland, P. W., Laskey, K. B. and Leinhardt, S. (1983). Stochastic blockmodels: First steps. Social Networks 5 109 – 137.
Holland, P. W. and Leinhardt, S. (1981). An exponential family of probability distributions for directed graphs. Journal of the American Statistical
Association 76 33–50.
Jain, L., Jamieson, K. and Nowak, R. (2016). Finite sample prediction and
recovery bounds for ordinal embedding .
Karp, R. (1972). Reducibility among combinatorial problems. Springer.
Klopp, O. and Tsybakov, a. V. N., A. (2015). Oracle inequalities for network models and sparse graphon estimation. Annals of Statistics (to appear)
ArXiv:1507.04118.
Koltchinskii, V., Lounici, K. and Tsybakov, A. (2011). Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion. The Annals
of Statistics 39 2302–2329.
Kučera, L. (1995). Expected complexity of graph partitioning problems. Discrete Appl. Math. 57 193–212. Combinatorial optimization 1992 (CO92) (Oxford).
Ma, Z. and Wu, Y. (2015). Computational barriers in minimax submatrix
detection. The Annals of Statistics 43 1089–1116.
Madeira, S. C. and Oliveira, A. L. (2004). Biclustering algorithms for biological data analysis: A survey. IEEE/ACM Trans. Comput. Biol. Bioinformatics
1 24–45.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

21

Massart, P. (2007). Concentration inequalities and model selection, vol. 6.
Springer.
Massoulié, L. (2014). Community detection thresholds and the weak Ramanujan property. In Proceedings of the 46th Annual ACM Symposium on Theory
of Computing. ACM.
Meier, L., van de Geer, S. and Bühlmann, P. (2008). The group lasso for
logistic regression. Journal of the Royal Statistical Society: Series B (Statistical
Methodology) 70 53–71.
Mendelson, S., Pajor, A. and Tomczak-Jaegermann, N. (2008). Uniform
uncertainty principle for bernoulli and subgaussian ensembles. Constructive
Approximation 28 277–289.
Mossel, E., Neeman, J. and Sly, A. (2013). A proof of the block model
threshold conjecture. arXiv:1311.4115 .
Mossel, E., Neeman, J. and Sly, A. (2015). Reconstruction and estimation
in the planted partition model. Probability Theory and Related Fields 162
431–461.
Negahban, S. and Wainwright, M. (2011). Estimation of (near) low-rank
matrices with noise and high-dimensional scaling. The Annals of Statistics 39
1069–1097.
Nickl, R. and van de Geer, S. (2013). Confidence sets in sparse regression.
Ann. Statist. 41 2852–2876.
Oymak, S., Jalali, A., Fazel, M., Eldar, Y. and Hassibi, B. (2015). Simultaneously structured models with application to sparse and low-rank matrices.
IEEE Transactions on Information Theory 61 2886–2908.
Papa, G., Bellet, A. and Clémençon, S. (2016). On graph reconstruction via
empirical risk minimization: Fast learning rates and scalability. In Advances
in Neural Information Processing Systems 29 (D. D. Lee, M. Sugiyama, U. V.
Luxburg, I. Guyon and R. Garnett, eds.). Curran Associates, Inc., 694–702.
Penrose, M. (2003). Random Geometric Graphs. Oxford University Press.
Richard, E., Savalle, P. and Vayatis, N. (2012). Estimation of simultaneously sparse and low-rank matrices. In Proceedings of the 29th International
Conference on Machine Learning (ICML-12).
Rigollet, P. (2012). Kullback-leibler aggregation and misspecified generalized
linear models. The Annals of Statistics 40 639–665.
Rohde, A. and Tsybakov, A. (2011). Estimation of high-dimensional low-rank
matrices. The Annals of Statistics 39 887–930.
Traonmilin, Y. and Gribonval, R. (2015). Stable recovery of low-dimensional
cones in hilbert spaces: One rip to rule them all .
URL https://arxiv.org/abs/1510.00504

22

Tsybakov, A. (2008).
Springer.

BALDIN AND BERTHET

Introduction to Nonparametric Estimation.

1st ed.

van de Geer, S. (2008). High-dimensional generalized linear models and the
lasso. The Annals of Statistics 36 614–645.
Wang, T., Berthet, Q. and Plan, Y. (2016a). Average-case hardness of rip
certification. In Proceedings of the 30th International Conference on Neural
Information Processing Systems. NIPS’16, Curran Associates Inc., USA.
URL http://dl.acm.org/citation.cfm?id=3157382.3157525
Wang, T., Berthet, Q. and Samworth, R. (2016b). Statistical and computational trade-offs in estimation of sparse principal components. The Annals
of Statistics 44 1896–1930.
Wasserman, S. and Faust, K. (1994). Social Network Analysis: Methods and
Applications. Structural Analysis in the Social Sciences, Cambridge University
Press.
Wolfe, P. and Olhede, S. (2013). Nonparametric graphon estimation. arXiv
preprint arXiv:1309.5936 .
Yu, H., Braun, P. and Yildirim, M. (2008). High quality binary protein
interaction map of the yeast interactome network. Science (New York, NY)
322 104–110.
Zhang, A. and Dong, X. (2017). Tensor svd: Statistical and computational
limits. arXiv preprint arXiv:1703.02724 .
Zhang, Y., Levina, E. and Zhu, J. (2015). Estimating network edge probabilities by neighborhood smoothing. arXiv preprint arXiv:1509.08588 .
Zhang, Y., Wainwright, M. and Jordan, M. (2014). Lower bounds on the
performance of polynomial-time algorithms for sparse linear regression 35.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

23

APPENDIX A: PROOFS
A.1 Some geometric properties of the likelihood
Let us recall the stochastic component of the likelihood function
X

Y(i,j) − πij (Θ⋆ ) Xi⊤ ΘXj ,
ζ(Θ) = ℓY (Θ) − ℓ(Θ) =
(i,j)∈Ω

which is a linear function in Θ. The deviation of the gradient ∇ζ of the stochastic
component is governed by the deviation of the independent Bernoulli random
variables εi,j = Y(i,j) − IE[Y(i,j) ] = Y(i,j) − πij (Θ⋆ ), (i, j) ∈ Ω. Let us introduce an
upper triangular matrix EΩ = (εi,j )(i,j)∈Ω with zeros on the complement set Ωc .
In this notation, we have ζ(Θ) = hhζ, ΘiiF , with
X
∇ζ =
εi,j Xj Xi⊤ = X EΩ X⊤ ∈ Rd×d .
(i,j)∈Ω

P
In particular, ∇ζ is sub-Gaussian with parameter (i,j)∈Ω kXj Xi⊤ k2F /4 = kX⊤ Xk2F,Ω /4,
i.e. it holds for the moment generating function of hhζ, BiiF for any B ∈ Rd×d
and σ 2 = 1/4,
Y




IE exp(tεi,j hhXj Xi⊤ , BiiF )
ϕhhζ,B ii (t) := IE exp(thhζ, BiiF ) =
F
(i,j)∈Ω

(A.1)

≤

Y

2 2

exp t σ

(i,j)∈Ω

2
hhXj Xi⊤ , BiiF /2




= exp tσ 2 kX⊤ BXk2F,Ω /2 .

We shall be frequently using versions of the following inequality, which is based
on the fact that ∇ℓ(Θ⋆ ) = 0 ∈ Rd×d , the Taylor expansion and (1.5), and holds
for any Θ ∈ P(M ),
ℓ(Θ⋆ ) − ℓ(Θ) =
≥

(A.2)

1 X
2

(i,j)∈Ω

σ ′ (Xi⊤ Θ0 Xj )hhXj Xi⊤ , Θ⋆ − Θii2



L X
L
hhXj Xi⊤ , Θ⋆ − Θii2F = kX⊤ (Θ⋆ − Θ)Xk2F,Ω ,
2
2
(i,j)∈Ω

where Θ0 ∈ [Θ, Θ⋆ ] element-wise. Furthermore, using that supt∈R σ ′ (t) ≤ 1/4, we
obtain for all Θ ∈ P(M )
ℓ(Θ⋆ ) − ℓ(Θ) ≤

1 ⊤
kX (Θ⋆ − Θ)Xk2F,Ω .
8

We shall also be using the bounds
(A.3)
(A.4)


max εi,j XiT (Θ − Θ⋆ )Xj ≤ kX⊤ (Θ − Θ⋆ )XkF,Ω ,

(i,j)∈Ω

 1
Var hhEΩ , X⊤ (Θ − Θ⋆ )XiiF ≤ kX⊤ (Θ − Θ⋆ )Xk2F,Ω .
4

a.s. ,

24

BALDIN AND BERTHET

A.2 Entropy bounds for some classes of matrices
Recall that an ε-net of a bounded subset K of some metric space with a
metric ρ is a collection {K1 , ..., KNε } ∈ K such that for each K ∈ K, there
exists i ∈ {1, ..., Nε } such that ρ(K, Ki ) ≤ ε. The ε-covering number N (ε, K, ρ)
is the cardinality of the smallest ε-net. The ε-entropy of the class K is defined
by H(ε, K, ρ) = log2 N (ε, K, ρ). The following statement is adapted from Lemma
3.1 in Candes and Plan (2011).
Let T0 := {Θ ∈ Rk×k : rank(Θ) ≤ r, kΘkF ≤ 1}. Then it holds

Lemma 22.
for any ε > 0


9
H(ε, T0 , k · kF ) ≤ (2k + 1)r + 1 log
.
ε

A.3 Proof of Theorem 10 and Theorem 16
It suffices to show the following uniform deviation inequality
(A.5)

sup
Θ⋆ ∈Pk,r (M )


PΘ⋆ ℓ(Θ⋆ ) − ℓ(Θ̂) + p(Θ̂) > 2p(Θ⋆ ) + Rt2 ≤ e−cRt ,

for any Rt > 0 and some numeric constant c > 0. Indeed, then taking Rt2 = p(Θ⋆ ),
it follows that ℓ(Θ⋆ ) − ℓ(Θ̂) ≤ 3p(Θ
√ ⋆ ) uniformly for all Θ⋆ in the considered class
−c p(Θ⋆ )
. The upper bound (2.3) of Theorem 10
with probability at least 1 − e
the follows directly integrating the deviation inequality (A.5), while the upper
bound on the prediction error in Theorem 16 further follows using (A.2) and the
smoothness of the logistic function, supt∈R σ ′ (t) ≤ 1/4. Define
τ 2 (Θ; Θ⋆ ) := ℓ(Θ⋆ ) − ℓ(Θ) + p(Θ) , GR (Θ⋆ ) := {Θ : τ (Θ; Θ⋆ ) ≤ R} .

(A.6)

The inequality (A.5) clearly holds on the event {τ 2 (Θ̂; Θ⋆ ) ≤ 2p(Θ⋆ )}. In view of
ℓY (Θ̂) − p(Θ̂) ≥ ℓY (Θ⋆ ) − p(Θ⋆ ), we have on the complement:
1
hhEΩ , X⊤ (Θ̂ − Θ⋆ )Xii ≥ ℓ(Θ⋆ ) − ℓ(Θ̂) + p(Θ̂) − p(Θ⋆ ) ≥ τ 2 (Θ̂; Θ⋆ ) .
2
Therefore, for any Θ⋆ ∈ Pk,r (M ), we have


PΘ⋆ τ 2 (Θ̂; Θ⋆ ) > 2p(Θ⋆ ) + Rt2 ≤ PΘ⋆

sup
τ (Θ;Θ⋆ )≥Rt

hhEΩ , X⊤ (Θ − Θ⋆ )Xii
1
.
≥
τ 2 (Θ; Θ⋆ )
2

We now apply the so-called “peeling device” (or “slicing” as it sometimes called
in the literature). The idea is to “slice” the set τ (Θ; Θ⋆ ) ≥ Rt into pieces on
which the penalty term p(Θ) is fixed and the term ℓ(Θ⋆ ) − ℓ(Θ) is bounded. It
follows,
PΘ⋆



sup
τ (Θ;Θ⋆ )≥Rt

1
hhEΩ , X⊤ (Θ − Θ⋆ )Xii
≥
τ 2 (Θ; Θ⋆ )
2

(A.7)
≤

∞
K X
d X
X

K=1 R=1 s=1

PΘ⋆



sup
Θ∈G2s Rt (Θ⋆ )
k(Θ)=K, rank(Θ)=R


1
hhEΩ , X⊤ (Θ − Θ⋆ )Xii ≥ 22s Rt2 .
8

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

25

On the set {Θ ∈ G2s Rt (Θ⋆ ), k(Θ) = K, rank(Θ) = R}, it holds by the definitions
(A.6)
ℓ(Θ⋆ ) − ℓ(Θ) ≤ 22s Rt2 − p(K, R) ,
and therefore using (A.2), this implies
(A.8) kX⊤ (Θ⋆ − Θ)XkF,Ω ≤ Z(K, R, s) ,

Z 2 (K, R, s) =


2 2s 2
2 Rt − p(K, R) .
L

Let us fix the location of the block, that is the support of a matrix Θ′ ∈ G1 :=
{Θ ∈ Rd×d : k(Θ) = K, rank(Θ) = R} belongs to the upper-left block of size
K × K. Then following the lines of the proof of Lemma 22 and using the singular
value decomposition, we derive

9B 
.
H(ε, {X⊤ Θ′ X : Θ′ ∈ G1 , kX⊤ Θ′ XkF,Ω ≤ B}, k · kF,Ω ) ≤ (2K + 1)R + 1 log
ε
Consequently, for the set T := {X⊤ (Θ − Θ⋆ )X : Θ ∈ Rd×d , rank(Θ) = R, k(Θ) =
K, kX⊤ (Θ⋆ − Θ)XkF,Ω ≤ Z(K, R, s)}, we obtain
 9Z(K, R, s) 

de 
H(ε, T, k · kF,Ω ) ≤ (2K + 1)R + 1 log
+ K log
.
ε
K
q
√

Denote t(K, R) := KR + K log de
K . By Dudley’s entropy integral bound, see
Dudley (1967) and Giné and Nickl (2016) for a more recent reference, we then
have
Z Z(K,R,s) p


⊤
′
H(ε, T, k · kΩ )dε
IE
sup
hhEΩ , X (Θ − Θ⋆ )Xii ≤ C
Θ∈G2s Rt (Θ⋆ )
k(Θ)=K, rank(Θ)=R

√ Z
≤ C kr

9Z(K,R,s)

′′

0

0

r

≤ CZ(K, R, s)t(K, R) ,

r
 9Z(K, R, s) 
de 
′′
log
dε + 9C Z(K, R, s) K log
ε
K

for some universal constant C > 0. Furthermore, by Bousquet’s version of Talagrand’s inequality, see Theorem A.10, in view of the bounds (A.3) and (A.4), we
have for all u > 0

hhEΩ , X⊤ (Θ − Θ⋆ )Xii ≥ CZ(K, R, s)t(K, R)
sup
PΘ⋆
Θ∈G2s Rt (Θ⋆ )
k(Θ)=K, rank(Θ)=R

+

r


Z(K, R, s)u 
1 2
≤ e−u .
Z (K, R, s) + 4CZ 2 (K, R, s)t(K, R) u +
2
3

Taking u(K, R, s) := L1/2 Z(K, R, s)+L−1/2 t(K, R)+2 log d and using inequalities
√
√
√
√
c1 + c2 ≤ c1 + c2 and c1 c2 ≤ 21 (c1 ε + cε2 ), which hold for any c1 , c2 , ε > 0,
we obtain

hhEΩ , X⊤ (Θ − Θ⋆ )Xii
sup
PΘ⋆
Θ∈G2s Rt (Θ⋆ )
k(Θ)=K, rank(Θ)=R

≥


1
LZ 2 (K, R, s) + C12 t2 (K, R)/L ≤ e−u(K,R,s) ,
16

26

BALDIN AND BERTHET

for some numeric constant C1 > 0. Plugging this back into (A.7) and using (A.8),
we obtain


1
hhEΩ , X⊤ (Θ − Θ⋆ )Xii ≥ 22s Rt2 ≤ e−u(K,R,s) ,
sup
PΘ⋆
8
Θ∈G2s R (Θ⋆ )
t

k(Θ)=K, rank(Θ)=R

for some numeric constant C2 > 0, provided that
1
C2
1
1
LZ 2 (K, R, s) + 1 t2 (K, R) ≤ 22s Rt2 = LZ 2 (K, R, s) + 8p(K, R) ,
16
L
8
16

(A.9)

which is satisfied for p(K, R) ≥ (C12 /L)t2 (K, R). Therefore,
PΘ⋆



sup
τ (Θ;Θ⋆ )≥Rt

K ∞
d
1  X X X −u(K,R,s)
hhEΩ , X⊤ (Θ − Θ⋆ )Xii
e
≤ e−cRt ,
≤
≥
τ 2 (Θ; Θ⋆ )
2
K=1 R=1 s=1

for some numeric constants c > 0 using (A.9), which concludes the proof.
The following prominent result is due to Bousquet (2002).
Theorem 23 (Bousquet’s version of Talagrand’s inequality). Let (B, B) be a
measurable space and let ε1 , ..., εn be independent B-valued random variables.
Let F be a countable set of measurable real-valued functions on B such that
f (εi ) ≤ b < ∞ a.s. and IEf (εi ) = 0 for all i = 1, ..., n, f ∈ F . Let
S := sup

n
X

f ∈F i=1

f (εi ) ,

v := sup

n
X

f ∈F i=1

IE[f 2 (εi )].

Then for all u > 0, it holds that

p
bu 
≤ e−u .
(A.10)
P S − IE[S] ≥ 2(v + 2bIE[S])u +
3
A.4 Proof of Theorem 13
For the MLE Θ̂k,r , it clearly holds ℓY (Θ̂k,r ) ≥ ℓY (Θ⋆ ) implying
ℓ(Θ⋆ ) − ℓ(Θ̂k,r ) ≤ hhEΩ , X⊤ (Θ̂k,r − Θ⋆ )Xii.
Furthermore, in view of (A.2), we derive
(A.11)

hhEΩ , X⊤ (Θ̂k,r − Θ⋆ )Xii
L ⊤
kX (Θ̂k,r − Θ⋆ )XkF,Ω ≤
2
kX⊤ (Θ̂k,r − Θ⋆ )XkF,Ω

(A.12)

≤

hhEΩ , X⊤ (Θ − Θ⋆ )Xii
.
⊤
Θ∈Pk,r (M ) kX (Θ − Θ⋆ )XkF,Ω
sup

Following the lines of Section A.3, by Dudley’s integral we next obtain
r

√
hhEΩ , X⊤ (Θ − Θ⋆ )Xii 
de 
≤ c kr + c k log
,
IE
sup
⊤
k
Θ∈Pk,r (M ) kX (Θ − Θ⋆ )XkF,Ω
for some universal constant c > 0. Plugging this bound back into (A.12) and
using the block isometry property yields the desired assertion.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

27

A.5 Proof of Theorem 17
Proof. The proof is split into two parts. First, we show a lower bound of the
order kr and then a lower bound of the order k log(de/k). A simple inequality
(a + b)/2 ≤ max{a, b} for all a, b > 0 then completes the proof. Both parts of the
proof exploit a version of remarkable Fano’s inequality given in Proposition 24 to
follow, cf. Section 2.7.1 in Tsybakov (2008).
1. A bound kr. The proof of this bound is similar to the proof of a minimax
lower bound for estimating a low-rank matrix in the trace-norm regression model
given in Theorem 5 in Koltchinskii et al. (2011). For the sake of completeness,
we provide the details here. Consider a subclass of matrices

C = A ∈ Rk×r : ai,j = {0, αN }, 1 ≤ i ≤ k, 1 ≤ j ≤ r ,
γ log 2
r
α2N =
,
1 + ∆Ω,2k (X) 2kN
where γ > 0 is a positive constant, ∆Ω,2k (X) > 0 is the the block isometry
constant from Definition 3 and ⌊x⌋ denotes the integer part of x. Further define
B(C) =

1
(A + A⊤ ) : A = (Ã| · · · |Ã|O) ∈ Rk×k , Ã ∈ C ,
2

where O denotes the k × (k − r⌊k/r⌋) zero matrix. By construction, any matrix
Θ ∈ B(C) is symmetric, has rank at most r with entries bounded by αN . Applying
a standard version of the Varshamov-Gilbert lemma, see Lemma 2.9 in Tsybakov
(2008), there exists a subset B ◦ ⊂ B(C) of cardinality card(B ◦ ) ≥ 2kr/16 + 1 such
that
kr αN 2  k 
≤ kΘu − Θv k2F ≤ k2 α2N ,
16 2
r
for all Θu , Θv ∈ B ◦ . Thus B ◦ is a 2δ-separated set in the Frobenius metric with
αN 2 k
⌊ r ⌋. The Kullback-Leibler divergence between the measures PΘu
δ2 = kr
64 2
and PΘv , Θu , Θv ∈ B ◦ , u 6= v, is upper bounded as
1 X
KL(PΘu , PΘv ) = IEPΘu [ℓY (Θu )] − IEPΘu [ℓY (Θv )] ≤
hhXj Xi⊤ , Θu − Θv ii2F
8
(i,j)∈Ω

1 + ∆Ω,2k (X) 2 2
≤
k αN N .
8
Taking γ > 0 small enough, we obtain
kr
1 + ∆Ω,2k (X) 2 2
kr
k αN N + log 2 =
γ log 2 + log 2 = log(2 16 γ+1 ) < log(2kr/16 + 1) ,
8
16
which, in view of Proposition 24, yields
 the desired lower bound.
2. A bound k log(de/k). Let K = kd and consider the set GkαN ⊂ Pk,1 (M ) from
the reduction scheme in Section 3.2 with
de 
4γ log 2
log
,
α2N =
kN (1 + ∆Ω,2k (X))
k

where γ > 0 is a positive constant. Using simple calculations, we then have
(2k − 1)α2N ≤ kΘu − Θv k2F ≤ 2k2 α2N for all Θu , Θv ∈ GkαN , u 6= v. Furthermore,
according to Lemma 25 to follow, there exists a subset GkαN ,0 ⊂ GkαN such that
c0 k2 α2N ≤ kΘu − Θv k2F ≤ 2k2 α2N ,

28

BALDIN AND BERTHET

and of cardinality card(GkαN ,0 ) ≥ 2ρk log(de/k) + 1 for some ρ > 0 depending on a
constant c0 > 0 and independent of k and d. Thus GkαN ,0 is a 2δ-separated set
in the Frobenius metric with δ2 = c0 k2 α2N /4. The Kullback-Leibler divergence
between the measures PΘu and PΘv , Θu , Θv ∈ GkαN ,0 , u 6= v, is upper bounded
as
1 X
hhXj Xi⊤ , Θu − Θv ii2F
KL(PΘu , PΘv ) = IEPΘu [ℓY (Θu )] − IEPΘu [ℓY (Θv )] ≤
8
(i,j)∈Ω

≤

1 + ∆Ω,2k (X) 2 2
k αN N ,
4

for all u 6= v and ∆Ω,2k (X) > 0 from Definition 3. As in the first part of the proof,
taking γ > 0 small enough, we obtain
1 + ∆Ω,2k (X) 2 2
de 
k αN N + log 2 = kγ log(2) log
+ log 2 = log(2kγ log(de/k)+1 )
4
k
< log(2ρk log(de/k) + 1) .
The desired lower bound then follows from Proposition 24.
Proposition 24 (Fano’s method). Let {Θ1 , ..., ΘJ } be a 2δ-separated set in
in the Frobenius metric, meaning that kΘk − Θl kF ≥ 2δ for all elements
Θk , Θl , l 6= k in the set. Then for any increasing and measurable function F :
[0, ∞) → [0, ∞), the minimax risk is lower bounded as
P
2


u,v KL(PΘu , PΘv )/J + log 2 
inf sup IEPΘ F (kΘ̂ − ΘkF ) ≥ F (δ) 1 −
.
log J
Θ̂ Θ

Rd×d

Lemma 25 (Variant of the Varshamov-Gilbert lemma). Let G ⊂ Pk,1 (M ) be
a set of {0, 1}d×d symmetric block-sparse matrices with the
size of the block k,
d
where k ≤ αβd for
some
α,
β
∈
(0,
1).
Denote
K
=
the
cardinality of G
k
P
′
′
and ρH (E, E ) = i,j 1(Ei,j 6= Ei,j ) the Hamming distance between two matrices
E, E ′ ∈ G. Then there exists a subset G 0 = {E (0) , ..., E (J) } ⊂ G of cardinality
log J := log(card(G 0 )) ≥ ρk log(
where ρ =

α
− log(αβ) (− log β

de
),
k

+ β − 1) such that
ρH (E (k) , E (l) ) ≥ ck 2 ,

for all k 6= l where c = 2(1 − α2 ) ∈ (0, 2).
Proof. Let E (0) = {0}k×k , D = ck 2 , and construct the set E1 = {E ∈ G :
ρH (E (0) , E) > D}. Next, pick any E (1) ∈ E1 and proceed iteratively so that for a
matrix E (j) ∈ Ej we construct the set
Ej+1 = {E ∈ Ej : ρ(E (j) , E) > D} .
Let J denote the last index j for which Ej 6= ∅. It remains to bound the cardinality
J of the constructed set G 0 = {E (0) , ..., E (J) }. For this, we consider the cardinality
nj of the subset {Ej \ Ej+1 }:
nj := #{Ej \ Ej+1 } ≤ #{E ∈ G : ρH (E (j) , E) ≤ D}.

OPTIMAL LINK PREDICTION WITH MATRIX LOGISTIC REGRESSION

29

For all E, E ′ ∈ G, we have
ρH (E, E ′ ) = 2(k2 − (k − m)2 ) ,
where m ∈ [0, k] corresponds to the number of distinct columns of E (or E ′ ).
Solving the quadratic equation (A.5) for ρH (E, E ′ ) = D = ck 2 we obtain
p
mD = k(1 − 1 − c/2) ,
for the maximum number of distinct columns of a block-sparse matrix E (and
E ′ ) such that ρH (E, E ′ ) ≤ D = ck 2 for c ∈ [0, 2]. For instance, in order to get the
distance between matrices 2k2 , i.e. c = 2 we need to shift all the k columns (and
consequently rows) and so the number of distinct columns of a matrix is m = k,
and in order to get the minimal possible distance 4k − 2, i.e. c = (4k − 2)/k2 we
need to shift only one column and a corresponding row, i.e. m = 1. Therefore,
for nj in (A.5), we have
nj ≤ #{E ∈ G : ρH (E

(j)

, E) ≤ D} =

Together with an evident equality
k
X

i=k−mD


mD  
X
k
d−k
i=0

PJ

j=0 nj

i

i

=K=

d
k ,

=

k
X

i=k−mD

 

k
d−k
.
i
k−i

this implies

 
  
k
d−k
d
1
.
/
≥
J +1
i
k−i
k

Note that taking mD = k, which as we have seen corresponds to c = 2, we
have a trivial bound J + 1 ≥ 1 using Vandermonde’s convolution. Furthermore, the expression on the left-hand side
p in (A.5) is exactly the probability
P(X ≥ k − mD ) = P(X ≥ kα) for α = 1 − c/2, where the variable X follows
the hypergeometric distribution H(d, k, k/d). The rest of the proof is based on
applying Chernoff’s inequality and follows the scheme of the proof of Lemma 4.10
in Massart (2007).
Nicolai Baldin
Department of Pure Mathematics
and Mathematical Statistics,
University of Cambridge
Wilbeforce Road
Cambridge, CB3 0WB, UK
(n.baldin@statslab.cam.ac.uk)

Quentin Berthet
Department of Pure Mathematics
and Mathematical Statistics,
University of Cambridge
Wilbeforce Road
Cambridge, CB3 0WB, UK
(q.berthet@statslab.cam.ac.uk)

