arXiv:1802.06292v2 [stat.ML] 6 Mar 2018

Nonparametric Estimation of Low Rank
Matrix Valued Function
Fan Zhou∗
School of Mathematics
Georgia Institute of Technology
Atlanta, GA 30332-0160
e-mail: fzhou40@math.gatech.edu
Abstract: Let A : [0, 1] → Hm (the space of Hermitian matrices) be
a matrix valued function which is low rank with entries in Hölder class
Σ(β, L). The goal of this paper is to study statistical estimation of A based
on the regression model E(Yj |τj , Xj ) = hA(τj ), Xj i, where τj are i.i.d.
uniformly distributed in [0, 1], Xj are i.i.d. matrix completion sampling
matrices, Yj are independent bounded responses. We propose an innovative
nuclear norm penalized local polynomial estimator and establish an upper
bound on its point-wise risk measured by Frobenius norm. Then we extend
this estimator globally and prove an upper bound on its integrated risk
measured by L2 -norm. We also propose another new estimator based on
bias-reducing kernels to study the case when A is not necessarily low rank
and establish an upper bound on its risk measured by L∞ -norm. We show
that the obtained rates are all optimal up to some logarithmic factor in
minimax sense. Finally, we propose an adaptive estimation procedure based
on Lepski’s method and the penalized data splitting technique which is
computationally efficient and can be easily implemented and parallelized.
MSC 2010 subject classifications: Primary 62G05, 62G08; secondary
62H12.
Keywords and phrases: Nonparametric estimation, low rank, matrix
completion, nuclear norm penalization, local polynomial estimator, model
selection.

1. Introduction
Let A : [0, 1] → Hm (the space of Hermitian matrices) be a matrix valued
function. The goal of this paper is to study the problem of statistical estimation
of A based on the regression model
E(Yj |τj , Xj ) = hA(τj ), Xj i, j = 1, . . . , n,

(1.1)

where τj are i.i.d. time design variables uniformly distributed in [0, 1], Xj are
i.i.d. matrix completion sampling matrices, Yj are independent bounded random
responses. Sometimes, it will be convenient to write model (1.1) in the form
Yj = hA(τj ), Xj i + ξj , j = 1, . . . , n,

(1.2)

where the noise variables ξj = Yj − E(Yj |τj , Xj ) are independent and have zero
means. In particular, we are interested in the case when A is low rank and
∗ Supported

in part by NSF Grants DMS-1509739 and CCF-1523768.
1

imsart-generic ver. 2014/10/16 file: matrix_func.tex date: March 7, 2018

satisfies certain smoothness condition. When A(t) = A0 with some A0 ∈ Hm
and for any t ∈ [0, 1], such problem coincides with the well known matrix completion/recovery problem which has drawn a lot of attention in the statistics
community during the past few years, see [7, 5, 6, 8, 16, 13, 18, 27, 25, 9]. The
low rank assumption in matrix completion/estimation problems has profound
practical background. For instance, when [20] introduced their famous work on
matrix factorization techniques for recommender systems, they considered temporal dynamics, see [19]. Another very common example is Euclidean distance
matrix (EDM) which contains the distance information of a large set of points
like molecules which are in low dimensional spaces such as R2 or R3 . To be more
specific, given m points p1 , ..., pm in Rd , the EDM D ∈ Rm×m formed by them
has entries Dij = kpi − pj k22 . Clearly, this matrix has rank at most d + 1 regardless of its size m. If m ≫ d, then the recovery problem falls into the low rank
realm. Similar topics in cases when points are fixed (see [34]) or in rigid motion
(see [28]) have been studied. While points are moving in smooth trajectories,
the EDMs are naturally high dimensional low rank matrix valued functions.
An appealing way to address the low rank issue in matrix completion/estimation
is through nuclear norm minimization, see [26]. In section 3, we inherit this idea
and propose a local polynomial estimator (see [11]) with nuclear norm penalization:
τ − t 
E2
DX
1 X  τj − t0 
j
0
Si p i
+ εkSk1 .
K
Yj −
, Xj
Sbh = arg min
S∈D nh
h
h
i=0
j=1
ℓ

n

(1.3)
where D ⊂ H(ℓ+1)m is a closed subset of block diagonal matrices with Sj ∈ Hm
on its diagonal, and {pi } is a sequence of orthogonal polynomials with nonnegative weight function K. The solution to the convex optimization problem (1.3)
P
induces a pointwise estimator of A(t0 ): Sbh (t0 ) := ℓi=0 Sbih pi (0). We prove that
2β/(2β+1) 
on the
under mild conditions, Sbh (t0 ) achieves a rate of O mr log n
n

2
pointwise risk measured by m12 Sbh (t0 ) − A(t0 ) 2 over Hölder class Σ(β, L) with
low rank parameter r, where k·k2 denotes the Frobenius norm of a matrix. In secb based on the local results and prove
tion 4, we propose a new global estimator A
2β/(2β+1) 
mr
log
n
b achieves a rate of O
on the integrated risk measured
that A
n
R1
2
1
b − A(t) dt. Then we study another naive kernel
by L2 -norm, i.e. 2
A(t)
m

0

2

estimator Ã which can be used to estimate matrix valued functions which are
not necessarily low rank. This estimator is associated with another popular approach to deal with low rank recovery which is called singular value thresholding,


n 2β/(2β+1)
see [5, 18, 9]. We prove that Ã achieves a rate of O m log
measured
n
2

by supt∈[h,1−h] m12 Ã(t) − A(t) , where k · k denotes the matrix operator norm.
Note that those rates coincide with that of classical matrix recovery/estimation
setting when the smoothness parameter β → ∞.
An immediate question is whether the above rates are optimal. In section 5,
we prove that all the rates we established are optimal up to some logarithmic
factor in the minimax sense, which essentially verified the effectiveness of our
2

methodology.
As one may have noticed, there is an adaptation issue involved in (1.3).
Namely, one needs to choose a proper bandwidth h and a proper order of degree
ℓ of polynomials. Both parameters are closely related to the smoothness of A
which is unknown to us in advance. In section 6, we propose a model selection
procedure based on Lepskii’s method ([23]) and the work of [3] and [33]. We
prove that this procedure adaptively selects an estimator that achieves a rate
on the integrated risk measured by L2 -norm which is the smallest among all
candidates plus a negligible term. What is more important, such a procedure is
computationally efficient, feasible in high dimensional setting, and can be easily
parallelized.
The major contribution of our paper is that we generalized the recent developments of matrix completion/estimation theory to low rank matrix valued
function setting by proposing a new optimal estimation procedure. To our best
knowledge, no one has ever thoroughly studied such problems from a theoretical
point of view.
2. Preliminaries
In this section, we introduce some important definitions, basic facts, and notations for the convenience of presentation.
2.1. Notations
For any Hermitian matrices A, B ∈ Hm , denote hA, Bi = tr(AB) which is known
as the Hilbert-Schmidt inner product. Denote hA, BiL2 (Π) = EhA, XihB, Xi,
where Π denotes the distribution of X. The corresponding norm kAk2L2 (Π) is
given by kAk2L2 (Π) = EhA, Xi2 .
We use k·k2 to denote the Hilbert-Schimidt norm (Frobenus norm or Schatten
2-norm) generated by the inner product h·, ·i; k · k to denote the operator norm
(spectral norm) of a matrix: the largest singular value; k · k1 to denote the trace
norm (Schatten 1-norm or nuclear norm), i.e. the sum of singular values; |A| to
denote the nonnegative matrix with entries |Aij | corresponding to A.
We denote
2
σ 2 := Eξ 2 , σX
:= n−1

n
X
j=1

EXi2 , UX := kXk

L∞

.

2.2. Matrix Completion and statistical learning setting
The matrix completion setting refers to that the random sampling matrices Xj
are i.i.d. uniformly distributed on the following orthonormal basis X of Hm :
X := {Ekj : k, j = 1, ..., m},
3

where Ekk := ek ⊗ ek , k = 1, ..., m; Ejk := √12 (ek ⊗ ej + ej ⊗ ek ), 1 ≤ k < j ≤ m;
Ekj := √i2 (ek ⊗ ej − ej ⊗ ek ), 1 ≤ k < j ≤ m with {ej }m
j=1 being the canonical
m
basis of C . The following identities are easy to check when the design matrices
are under matrix completion setting:
A

2
L2 (Π)

=

1
2
A 2,
m2

2
σX
≤

2
,
m

UX = 1.

(2.1)

The statistical learning setting refers to the bounded response case: there exists
a constant a such that
max |Yj | ≤ a, a.s.
(2.2)
j=1,...n

In this paper, we will consider model (1.1) under both matrix completion and
statistical learning setting.
2.3. Matrix Valued Functions
Let A : [0, 1] → Hm be a matrix valued function. One should notice that we
consider the image space to be Hermitian matrix space for the convenience
of presentation. Our methods and results can be readily extended to general
rectangular matrix space. Now we define the rank of a matrix valued function.
Let rankA (t) := rank(A(t)), ∀t ∈ [0, 1].
Definition 1. Let β and L be two positive real numbers. The Hölder class
Σ(β, L) on [0, 1] is defined as the set of ℓ = ⌊β⌋ times differentiable functions
f : [0, 1] → R with derivative f (ℓ) satisfying
|f (ℓ) (x) − f (ℓ) (x′ )| ≤ L|x − x′ |β−ℓ , ∀x, x′ ∈ [0, 1].

(2.3)

In particular, we are interested in the following assumptions on matrix valued
functions.
A1 Given a measurement matrix X and for some constant a,
sup hA(t), Xi ≤ a.

t∈[0,1]

A2 Given a measurement matrix X and for some constant a, the derivative
matrices A(k) of A satisfy
sup hA(k) (t), Xi ≤ a, k = 1, ..., ℓ.

t∈[0,1]

A3 The rank of A, A′ , ...,A(ℓ) are uniformly bounded by a constant r,
sup rankA(k) (t) ≤ r, k = 0, 1, ..., ℓ.

t∈[0,1]

A4 Assume that for ∀i, j, aij is in the Hölder class Σ(β, L).
4

3. A local polynomial Lasso estimator
In this section, we study the estimation of matrix valued functions that are
low rank. The construction of our estimator is inspired by localization of nonparametric least squares and nuclear norm penalization. The intuition of the
localization technique originates from classical local polynomial estimators, see
[11]. The intuition behind nuclear norm penalization is that whereas rank function counts the number of non-vanishing singular values, the nuclear norm sums
their amplitude. The theoretical foundations behind nuclear norm heuristic for
the rank minimization was proved by [26]. Instead of using the trivial basis
{1, t, t2 , ..., tℓ } to generate an estimator, we use orthogonal polynomials which
fits our problem better. Let {pi }∞
i=0 be a sequence of orthogonal polynomials
with nonnegative weight function K compactly supported on [−1, 1], then
Z 1
K(u)pi (u)pj (u)du = δij .
−1

There exist an invertible linear transformation T ∈ R(ℓ+1)×(ℓ+1) such that
(1, t, t2 /2!, ..., tℓ /ℓ!)T = T (p0 , p1 , ..., pℓ )T .

Apparently, T is lower triangular. We denote R(T ) = max1≤j≤ℓ+1

ℓ+1
P
i=1

|Tij |. Note

that in some literature, R(T ) is denoted as kT k1 as the matrix ”column norm”.
Since we already used k · k1 to denote the nuclear norm, R(T ) is used to avoid
any ambiguity. Denote
)
(


⊂ Hm(ℓ+1)
D := Diag S0 S1 . . . Sℓ−1 Sℓ
the set of block diagonal matrices with Sk ∈ Hm satisfying |Sij | ≤ R(T )a. With
observations (τj , Xj , Yj ), j = 1, ..., n from model (1.1), define Sbh as
τ − t 
E2
DX
1 X  τj − t0 
j
0
Si p i
+ εkSk1 .
K
Yj −
, Xj
Sbh = arg min
S∈D nh
h
h
i=0
j=1
ℓ

n

(3.1)

Sbh naturally induces a local polynomial estimator of order ℓ around t0 :
Sbh (τ ) :=

τ − t  n τ − t
o
0
0
Sbih pi
I
≤1 .
h
h
i=0

ℓ
X

(3.2)

The point estimate at t0 is given by

Sbh (t0 ) :=

ℓ
X
i=0

Sbih pi (0).

(3.3)

In the following theorem, we establish an upper bound on the point-wise risk
of Sbh (t0 ) when A(t) is in the Hölder class Σ(β, L) with ℓ = ⌊β⌋.
5

Theorem 3.1. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a, a.s. for some constant a > 0. Let A be a matrix
valued function satisfying A1, A2, A3, and A4.
√
Denote Φ = maxi=0,...,ℓ k Kpi k∞ , and ℓ = ⌊β⌋. Take
s
1
 (ℓ3 (ℓ!)2 Φ2 R(T )2 a2 mr log n  2β+1
log 2m
b
,
hn = C1
, ε = DℓaΦ
L2 n
nmb
hn

for some numerical constants C1 and D. Then for any b
h n ≤ t0 ≤ 1 − b
hn , the
1
estimator defined in (3.3) satisfies with probability at least 1 − nmr
,
1 bh
S (t0 ) − A(t0 )
m2

2
2

2β
 mr log n  2β+1
,
≤ C1 (a, Φ, ℓ, L)
n

(3.4)

where C1 (a, Φ, ℓ, L) is a constant depending on a, Φ, ℓ and L.

The proof of Theorem 3.1 can be found in section 7.1. One should notice
that when β → ∞, bound (3.4) coincides with similar result in classical matrix
completion. In section 5, we prove that bound (3.4) is minimax optimal up to a
logarithmic factor.
4. Global estimators and upper bounds on integrated risk
In this section, we propose two global estimators and study their integrated risk
measured by L2 -norm and L∞ -norm.
4.1. From localization to globalization
Firstly, we construct a global estimator based on (3.2). Take
1
 ℓ3 (ℓ!)2 Φ2 R(T )2 a2 mr log n  2β+1
b
hn = C1
, M = ⌈1/b
hn ⌉.
L2 n

Without loss of generality, assume that M is even. Denote Sbkh (t) the local polynomial estimator around t2k−1 as in (3.2) by using orthogonal polynomials with
K = I[−1,1] , where t2k−1 = 2k−1
M , k = 1, 2, ..., M/2 and I is the indicator function. Denote
M/2

b =
A(t)

X

k=1

Sbkh (t)I(t2k−1 −bhn ,t2k−1 +bhn ] , t ∈ (0, 1).

(4.1)

Note that the weight function K is not necessary to be I[−1,1] . It can be replaced
by any K that satisfies K ≥ K0 > 0 on [−1, 1]. The following result characterizes the global performance of estimator (4.1) under matrix completion setting
measured by L2 -norm.
6

b be
Theorem 4.1. Assume that the conditions of Theorem 3.1 hold, and let A
1
,
an estimator defined as in (4.1). Then with probability at least 1 − nmr−1
1
m2

Z

0

1

b − A(t)
A(t)

2
dt
2

2β
 mr log n  2β+1
,
≤ C2 (a, Φ, ℓ, L)
n

(4.2)

where C2 (a, Φ, ℓ, L) is a constant depending on a, Φ, ℓ, L.
The proof of Theorem 4.1 can be found in section 7.2. Compared with the
integrated risk measured by L2 -norm of real valued functions in Hölder class,
the result in (4.2) has an excess log n term, which is introduced by the matrix
Bernstein inequality, see [30]. In section 5, we show that bound (4.2) is minimax
optimal up to a logarithmic factor.
4.2. Bias reduction through higher order kernels
If A(t) is not necessarily low rank, we propose an estimator which is easy to
implement and prove an upper bound on its risk measured by L∞ -norm. Such
estimators are related to another popular approach parallel to local polynomial
estimators for bias reduction, namely, using high order kernels to reduce bias.
They can also be applied to another important technique of low rank estimation
or approximation via singular value thresholding, see [5] and [9]. The estimator
proposed by [18] is shown to be equivalent to soft singular value thresholding of
such type of estimators.
The kernels we are interested in satisfy the following conditions:
K1
K2
K3
K4
K5

K(·) is symmetric, i.e. K(u) = K(−u).
K(·) isRcompactly supported on [−1, 1].
∞
RK = −∞ K 2 (u)du < ∞.
K(·) is of order ℓ for some ℓ ∈ N∗ .
K(·) is Lipschitz continuous with Lipschitz constant 0 < LK < ∞.

Consider
Ã(t) =

n
m2 X  τj − t 
K
Yj Xj .
nh j=1
h

(4.3)

Note that when K ≥ 0, (4.3) is the solution to the following optimization
problem
n
1 X  τj − t 
Ã(t) = arg min
K
(Yj − hS, Xj i)2 .
(4.4)
S∈D nh
h
j=1
In the following theorem we prove an upper bound on its global performance
measured by L∞ -norm over Σ(β, L) which is much harder to obtain for matrix
lasso problems.
Theorem 4.2. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a a.s. for some constant
7

a > 0; Kernel K satisfies K1-K5; let A be any matrix valued function satisfying
A1 and A4. Denote ℓ = ⌊β⌋. Take
1
 a2 (ℓ!)2 m log n  2β+1
,
(4.5)
h̃n := c∗ (K)
2βL2 n
Then with probability at least 1 − n−2 , the estimator defined in (4.3) satisfies
2β
 a2 (ℓ!)2 m log n  2β+1
1
2
∗
,
(4.6)
Ã(t)
−
A(t)
≤
C
(K)
sup
2
2βL2 n
t∈[h̃n ,1−h̃n ] m

where C ∗ (K) and c∗ (K) are constants depending on K.
The proof of Theorem 4.2 can be found in section 7.3. When the smoothness
parameter β tends to infinity, bound (4.6) coincides with similar bounds in
3
classical matrix completion, which is O( m nlog n ). When m degenerates to 1, the
bound coincides with that of real valued case, which is O(( logn n )2β/(2β+1) ). In
section 5, we show that this bound is minimax optimal up to a logarithmic
factor.
5. Lower bounds under matrix completion setting
In this section, we prove the minimax lower bound of estimators (3.3), (4.1)
and (4.3). In the realm of classical low rank matrix estimation, [25] studied
the optimality issue measured by the Frobenius norm on the classes defined
in terms of a ”spikeness index” of the true matrix; [27] derived optimal rates
in noisy matrix completion on different classes of matrices for the empirical
prediction error; [18] established that the rates of the estimator they propose
under matrix completion setting are optimal up to a logarithmic factor measured
by the Frobenius norm. Based on the ideas of [18], standard methods to prove
minimax lower bounds in real valued case in [31], and some fundamental results
in coding theory, we establish the corresponding minimax lower bounds of (3.4),
(4.2) and (4.6) which essentially shows that the upper bounds we get are all
optimal up to some logarithmic factor.
For the convenience of representation, we denote by inf Ab the infimum over
b of A. We denote by A(r, a) the set of matrix valued functions
all estimators A
satisfying A1, A2, A3, and A4. We denote by P(r, a) the class of distributions
of random triplet (τ, X, Y ) that satisfies model (1.1) with any A ∈ A(r, a).
Theorem 5.1. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a, a.s. for some constant a > 0; let A be any matrix valued function in A(r, a). Then there is an
absolute constant η ∈ (0, 1) such that for all t0 ∈ [0, 1]
n 1
 2β o

b 0 ) − A(t0 ) 2 > C(β, L, a) mr 2β+1 ≥ η.
sup
PPτ,X,Y
A
inf
A(t
2
b PA
m2
n
A
τ,X,Y ∈P(r,a)

(5.1)

8

where C(β, L, a) is a constant depending on β, L and a.
The proof of Theorem 5.1 can be found in section 7.4. Note that compared
with the upper bound (3.4), the lower bound (5.1) matches it that up to a
logarithmic factor. As a consequence, it shows that the estimator (3.3) achieves
a near optimal minimax rate of pointwise estimation. Although, the result of
Theorem 5.1 is under bounded response condition, it can be readily extended
to the case when the noise in (1.2) is Gaussian.
Theorem 5.2. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a, a.s. for some constant a > 0; let A be any matrix valued function in A(r, a). Then there is an
absolute constant η ∈ (0, 1) such that
n 1 Z 1
b
sup
PPτ,X,Y
A
inf
A(t)−
A(t)
2
b PA
m
A
0
∈P(r,a)
τ,X,Y

2
dt
2

2β o
 mr  2β+1
≥ η,
> C̃(β, L, a)
n

(5.2)

where C̃(β, L, a) is a constant depending on L, β and a.
The proof of Theorem 5.2 can be found in section 7.5. The lower bound
in (5.2) matches the upper bound we get in (4.2) up to a logarithmic factor.
Therefore, our estimator (4.1) achieves a near optimal minimax rate on the
integrated risk measured by L2 -norm up to a logarithmic factor. The result of
Theorem 5.2 can be readily extended to the case when the noise in (1.2) is
Gaussian.
Now we consider the minimax lower bound on integrated risk measured by
L∞ -norm for general matrix valued functions without any rank information.
Denote

A(a) = A(t) ∈ Hm , ∀t ∈ [0, 1] : |aij (t)| ≤ a, aij ∈ Σ(β, L) .

We denote by P(a) the class of distributions of random triplet (τ, X, Y ) that
satisfies model (1.1) with any A ∈ A(a).
Theorem 5.3. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a, a.s. for some constant a > 0; let A be any matrix valued function in A(a). Then there exist an
absolute constant η ∈ (0, 1) such that
inf

sup

b PA
A
τ,X,Y ∈P(a)

A
PPτ,X,Y

n

1 b
A(t)−A(t)
2
t∈(0,1) m
sup

2

2β o
 m ∨ log n  2β+1
> C̄(β, L, a)
≥ η.
n

(5.3)

where C̄(β, L, a) is a constant depending on β, L and a.
The proof of Theorem 5.3 can be found in section 7.6. Recall that in the real
valued case, the minimax lower bound measured by sup norm of Hölder class is
9

O(( logn n )2β/(2β+1) ). In our result (5.3), if the dimension m degenerates to 1, we
get the same result as in real valued case and it is optimal. While the dimension
m is large enough such that m ≫ log n, the lower bound (5.3) shows that the
estimator (4.3) achieves a near minimax optimal rate up to a logarithmic factor.
6. Model selection
Despite the fact that estimators (3.3) and (4.1) achieve near optimal minimax
rates in theory with properly chosen bandwidth h and order of degree ℓ, such
parameters depend on quantities like β and L which are unknown to us in advance. In this section, we propose an adaptive estimation procedure to choose h
and ℓ adaptively. Two popular methods to address such problems are proposed
in the past few decades. One is Lepskii’s method, and the other is aggregation
method. In the 1990s, many data-driven procedures for selecting the ”best” estimator emerged. Among them, a series of papers stood out and shaped a method
what is now called ”Lepskii’s method”. This method has been described in its
general form and in great detail in [23]. Later, [21] proposed a bandwidth selection procedure based on pointwise adaptation of a kernel estimator that achieves
optimal minimax rate of point estimation over Hölder class, and [22] proposed
a new bandwidth selector that achieves optimal rates of convergence over Besov
classes with spatially imhomogeneous smoothness. The basic idea of Lepskii’s
method is to choose a bandwidth from a geometric grid to get an estimator
not ”very different” from those indexed by smaller bandwidths on the grid. Although Lepskii’s method is shown to give optimal rates in pointwise estimation
over Hölder class in [21], it has a major defect when applied to our problem:
the procedure already requires a huge amount of computational cost when real
valued functions are replaced by matrix valued functions. Indeed, with Lepskii’s
method, in order to get a good bandwidth, one needs to compare all smaller
bandwidth with the target one, which leads to dramatically growing computational cost. Still, we have an extra parameter ℓ that needs to fit with h. As a
result, we turn to aggregation method to choose a bandwidth from the geometric grid introduced by Lepskii’s method, which is more computationally efficient
for our problem. The idea of aggregation method can be briefly summarized as
follows: one splits the data set into two parts; the first is used to build all candidate estimators and the second is used to aggregate the estimates to build
a new one (aggregation) or select one (model selection) which is as good as
the best candidate among all constructed. The model selection procedure we
use was initially introduced by [3] in classical nonparametric estimation with
bounded response. [33] generalized this method to the case where the noise can
be unbounded but with a finite p-th moment for some p > 2. One can find a
more detailed review on such penalization methods in [14].
Firstly, we introduce the geometric grid created by [21] where to conduct our
model selection procedure. Assume that the bandwidth being considered falls
into the range [hmin , hmax ]. Recall that the ”ideal” bandwidth b
hn which is given
10

as

1
 ℓ3 ℓ!ΦR(T )a2 mr log n  2β+1
b
,
(6.1)
hn = C1
L2 n
hmax , hmin can be chosen as
 ℓ3 ℓ !ΦR(T )a2 mr log n  2β 1+1
 ℓ∗3 ℓ∗ !ΦR(T )a2 mr log n  2β∗1+1
∗
∗
, hmin = C1 ∗
hmax = C1
,
2
∗2
L∗ n
L n

where [β∗ , β ∗ ] and [L∗ , L∗ ] are the possible ranges of β, L respectively. Obviously,
β is the most important parameter among all. Note that when those ranges are
not so clear, a natural upper bound of hmax is 1, and a typical choice of hmin
can be set to n−1/2 . Denote
r
r
h
h


1
max
max
.
d(h) = 1 ∨ 2 log
, dn = 2 log
, α(h) = p
h
hmin
d(h)
√
Apparently, dn = O( log n). Define grid H inductively by
n
o
hk
H = hk ∈ [hmin , hmax ] : h0 = hmax , hk+1 =
, k = 0, 1, 2, ... . (6.2)
1 + α(hk )

Note that the grid H is a decreasing sequence and the sequence becomes denser
as k grows.
We consider possible choices of ℓk for each hk . A trivial candidate set is
ℓk ∈ L := {⌊β∗ ⌋, ⌊β∗ ⌋ + 1, ..., ⌊β ∗ ⌋} ⊂ N∗ . If the size of this set is large,
one can shrink it through the correspondence (6.1) between h and β, ℓk ≤
1
n−1
mr log 2m
 (1− dlog) log
 log n−1 +log
−1 
−1 
log hk
hk
d
.
If
n
≥
m
for
some
d
>
1,
≤ ℓk ≤
2
2
n−1
 log

−1
log hk
, which indicates the more the data, the narrower the range. We
2
denote the candidate set for ℓ as L. Then the set
H̃ = H × L := {(h, ℓ) : h ∈ H, ℓ ∈ L}

indexed a countable set of candidate
estimators. Once (hk , ℓk ) is fixed, one can
q
log 2m
take εk = D(ℓk + 1)R(T )Φa nmhk .

Now we introduce our model selection procedure based on H̃. We split the
data (τj , Xj , Yj ), j = 1, ..., 2n, into two parts with equal size. The first part of the
observations {(τj , Xj , Yj ) : j ∈ ~n } contains n data points, which are randomly
drawn without replacement from the original data set. We construct a sequence
bk , k = 1, 2, ... based on the training data set ~n through (4.1) for
of estimators A
b among {A
bk }, which
each pair in H̃. Our main goal is to select an estimator A
is as good as the one that has the smallest mean square error. We introduce an
bk which serves as a penalty term.
quantity πk associated with each estimator A
We use the remaining part of the data set {(τj , Xj , Yj ) : j ∈ ℓn } to perform the
selection procedure:
1 X
bk (τj ), Xj i)2 + πk .
(Yj − hA
(6.3)
k ∗ = arg min
k n
n
j∈ℓn

11

∗

b∗ = A
bk as the adaptive estimator. In practice, we suggest one to rank
Denote A
bk according to the following rule: 1. pairs with bigger h always
all estimators A
have smaller index; 2. if two pairs have the same h, the one with smaller ℓ has
smaller index. Our selection procedure can be summarized in Algorithm 1:
1. Construct the geometric grid H defined as in (6.2) and compute the
candidate set H̃;
2. Equally split the data set (τj , Xj , Yj ), j = 1, ..., N into two parts (~n and
ℓn ) by randomly drawing without replacement;
bk defined as in (4.1) using
3. For each pair in H̃, construct an estimator A
data set ~n ;
4. Using the second data set ℓn to perform the selection rule defined as in
(6.3).
Algorithm 1: Model Selection Procedure
The selection procedure described in Algorithm 1 have several advantages:
firstly, it chooses a global bandwidth instead of a local one; secondly, since our
bk , no matrix
selection procedure is only based on computations of entries of A
computation is involved in the last step, which saves a lot of computational
cost and can be easily applied to high dimensional problems; finally, step 3 and
4 can be easily parallelized. The following theorem shows that the integrated
b∗ measured by L2 -norm can be bounded by the smallest one among all
risk of A
candidates plus an extra term of order n−1 which is negligible.
Theorem 6.1. Under model (1.1), let (τj , Xj , Yj ), j = 1, ..., 2n be i.i.d. copies
of the random triplet (τ, X, Y ) with X uniformly distributed in X , τ uniformly
distributed in [0, 1], X and τ are independent, and |Y | ≤ a, a.s. for some constant a > 0; let A be a matrix valued function satisfying A1, A2, A3, and A4;
bk } be a sequence of estimators constructed from H̃; let A
b∗ be the adaptive
let {A
1
estimator selected through Algorithm 1. Then with probability at least 1 − nmr
Z 1
Z 1
o
n
1
b∗ (t) − A(t) 2 dt ≤ 3 min 1
bk (t) − A(t) 2 dt + πk + C(a) ,
A
A
2
2
k
m2 0
m2 0
n
n
(6.4)
where C(a) is a constant depending on a.

The proof of Theorem 6.1 can be found in section 7.7. Recall that Card(H) =
O(log n), we can take πk = kmr. Then πk ≤ c1 mr log n uniformly for all k
with some numerical constant c1 . According to Lepskii’s method that at least
one candidate in H gives the optimal bandwidth associated with the unknown
smoothness parameter β, together with the result of Theorem 4.1, the following
b∗ is
corollary is a direct consequence of Theorem 4.1 and 6.1, which shows A
adaptive.
Corollary 6.1. Assume that the conditions of Theorem 6.1 hold with πk =
1
kmr, and n > mr log n. Then with probability at least 1 − nmr−1
Z 1
 2β

1
b∗ (t) − A(t) 2 dt ≤ C(a, Φ, ℓ, L) mr log n 2β+1
(6.5)
A
2
m2 0
n
12

where C(a, Φ, ℓ, L) is a constant depending on a, Φ, ℓ, and L.
7. Proofs
7.1. Proof of Theorem 3.1
Proof. Firstly, we introduce a sharp oracle inequality of ”locally integrated risk”
of estimator (3.2) in the following lemma.
Lemma 1. Assume that the condition of Theorem 3.1 holds. Then there exist
a numerical constants D > 0 such that for all
r log 2m _ (log 2m)Φ 
,
ε ≥ D(ℓ + 1)R(T )Φa
nmh
nh

and for arbitrary η > 0, the estimator (3.3) satisfies with probability at least
1 − e−η
E2
 τ − t D
1
0
A(τ ) − Sbh (τ ), X
EK
h
h
E2 D2 (ℓ + 1)2 Φ2 R(T )2 a2 (rank(S)m log 2m + η) o
 τ − t D
n1
0
A(τ ) − S(τ ), X +
.
EK
≤ inf
S∈D h
h
nh
(7.1)
ℓ

P
τ −t0
Si p i h .
where S(τ ) :=
i=0

The proof of Lemma 1 can be derived from Theorem 19.1 in [17]. To be more
specific, one just needs to rewrite (3.1) as
E2
D
1 X
+ εkSk1.
Sbh = arg min
Ỹj − S, X̃j
S∈D n
j=1
n

hr

1
hK



τj −t0
h

(7.2)

r
r
 

 

 
 i


τj −t0
τj −t0
τj −t0
τj −t0
τj −t0
1
1
p0
X
,
p
X
,
...,
p
Xj ,
K
K
j
1
j
ℓ
h
h
h
h
h
h
h

where X̃j = Diag
r


τ −t
and Ỹj = h1 K j h 0 Yj . Then Lemma 1 is simply an application of Theorem

19.1 in [17]. Indeed, when the loss function is set to squared loss, then all arguments in the proof of Lemma 1 can be reproduced from the original proof. Since
it is mostly tedious repeated arguments, we omit it here. One should notice that
in the original proof of Theorem 19.1 in [17], a matrix isometry condition needs
to be satisfied. That is kAk2L (Π̃) = µ0 kAk22 for some constant µ0 > 0 and any
2

A ∈ Hm with Π̃ being the distribution of X̃. One can easily check that it is true
with (7.2). This is also the primary reason why we abandoned the classical local polynomial estimator with the trivial polynomial basis and used orthogonal
polynomials instead.
13

Consider
 τ − t  E2
 τ − t D
X
1
0
0
Sbih pi
A(τ ) −
,X
EK
h
h
h
i=0
ℓ

ℓ
ℓ
τ − t  X
 τ − t  E2
 τ − t D
X
1
0
0
0
Si p i
(Si − Sbih )pi
A(τ ) −
+
,X
= EK
h
h
h
h
i=0
i=0

 τ − t  E2
 τ − t  E2 1
 τ − t D X
 τ − t D
X
1
0
0
0
0
Si p i
(Si − Sbih )pi
, X + EK
A(τ ) −
,X
EK
h
h
h
h
h
h
i=0
i=0
ℓ

=

ℓ

ℓ
ℓ
 τ − t  ED X
τ − t  E
 τ − t D
X
2
0
0
0
Si p i
(Si − Sbih )pi
A(τ ) −
,X
,X
EK
h
h
h
h
i=0
i=0
(7.3)
Therefore, from (7.1) and (7.3), we have for any S ∈ D

+

 τ − t  E2
 τ − t D X
1
0
0
(Si − Sbih )pi
,X
EK
h
h
h
i=0
ℓ

≤

ℓ
ℓ
 τ − t  ED X
τ − t  E
τ − t  D
X
2
0
0
0
Si p i
(Si − Sbih )pi
,X
,X
A(τ ) −
EK
h
h
h
h
i=0
i=0

D2 (ℓ + 1)2 Φ2 R(T )2 a2 (rank(S)m log 2m + η)
.
nh
ℓ
 τ − t  E2
 c4  1
 τ − t D
X
0
0
≤ 2
Si p i
A(τ ) −
,X
EK
c −1 h
h
h
i=0
 c2 n D2 (ℓ + 1)2 Φ2 R(T )2 a2 (rank(S)m log 2m + η) o
,
+ 2
c −1
nh

+

where we used the fact that for any positive constants a and b, 2ab ≤
for some c > 1. Take S such that
ℓ
X
i=0

Si p i

(7.4)
1 2
2 2
c2 a +c b

τ − t 
τ − t 
A(ℓ) (t0 )hℓ  τ − t0 ℓ
0
0
= A(t0 )+ A′ (t0 )h
+ ...+
. (7.5)
h
h
ℓ!
h

Note that this is possible since the right hand side is a matrix valued polynomial
0
of τ −t
up to order ℓ, and span{p0 , p1 , ..., pℓ } = span{1, x, ..., xℓ }. Under the
h
condition that all entries of A(k) (t) are bounded by a, then entries of Sk are
bounded by R(T )a. Thus, the corresponding S ∈ D. Obviously, rank(Si ) ≤
(ℓ + 1 − i)r. Since A ∈ Σ(β, L), we consider ℓ-th order Taylor expansion of A at
t0 to get
Ã(τ − t0 )ℓ
A(τ ) = A(t0 ) + A′ (t0 )(τ − t0 ) + ... +
,
(7.6)
ℓ!
(ℓ)

where Ã is the matrix with Ãij = aij (t0 + αij (τ − t0 )) for some αij ∈ [0, 1].
14

Then we apply the Taylor expansion (7.6) and identity (7.5) to get
 τ − t  E2
 τ − t D
 τ − t  1 LU (τ − t)β
X
1
1
0
0
0
Si p i
A(τ ) −
, X ≤ EK
EK
h
h
h
h
h
m2
ℓ!
i=0
(7.7)
where U denotes the matrix with all entries being 1. The first inequality is due
to aij ∈ Σ(β, L), and the second is due to |τ − t0 | ≤ h. Under the condition that
X is uniformly distributed in X , and the orthogonality of {pi }ℓi=0 , it is easy to
check that
ℓ

ℓ
ℓ
 τ − t  E2
 τ − t D X
1
1 X bh
0
0
h
b
(Si − Si )pi
,X = 2
kSi − Si k22
EK
h
h
h
m
i=0
i=0

2
2

≤

L2 h2β
.
(ℓ!)2

(7.8)

Note that

Sbh (t0 ) − S(t0 )

2
2

=

ℓ
X
i=0

(Sbih − Si )pi (0)

2
2

≤ (ℓ + 1)Φ2

ℓ
X
i=0

Sbih − Si

2
,
2

(7.9)

where the second inequality is due to Cauchy-Schwarz inequality and pi are
uniformly bounded on [-1,1]. Combining (7.4), (7.7), (7.8), and (7.9), we get
with probability at least 1 − e−η
 c4  2L2 h2β  c2 n D2 (ℓ + 1)2 Φ2 R(T )2 a2 (rank(S)m log 2m + η) o
1 bh
,
kS (t0 )−A(t0 )k22 ≤ 2
+ 2
2
m
c − 1 (ℓ!)2
c −1
nh

By optimizing the right hand side with respect to h and take η = mr log n, we
take
1
 ℓ3 (ℓ!)2 Φ2 R(T )2 a2 mr log n  2β+1
b
.
hn = C
L2 n
where C is a numerical constant. This completes the proof of the theorem.
7.2. Proof of Theorem 4.1
Proof. It is easy to see that
Z

1

0

b − A(t)k2 dt ≤
kA(t)
2

For each k,
Z t2k−1 +bhn
1
Sbkh (t) − A(t)
m2 t2k−1 −bhn

M/2 Z t2k−1 +b
hn
X
k=1

2
dt
2

t2k−1 −b
hn

kSbkh (t) − A(t)k22 dt.

(7.10)

D
E2
= Eτ,X I(t2k−1 −bhn ,t2k−1 +bhn ] A(τ ) − Sbh (τ ), X

By (7.1), (7.7) and arguments used to prove Theorem 3.1, we have with proba1
bility at least 1 − nmr
,
1
2
hn
m b

Z

t2k−1 +b
hn
t2k−1 −b
hn

Sbkh (t) − A(t)

2
dt
2

15

2β
 mr log n  2β+1
.
≤ C1 (a, Φ, ℓ, L)
n

We take the union bound over k, from (7.10) we get with probability at least
1
,
1 − nmr−1
1
m2

Z

0

1

 2β

b − A(t)k2 dt ≤ C2 (a, Φ, ℓ, L) mr log n 2β+1 .
kA(t)
2
n

where C2 (a, Φ, ℓ, L) is a constant depending on a, Φ, ℓ, L.
7.3. Proof of Theorem 4.2

Proof. In this proof, we use C(K) to denote any constant depending on K which
may vary from place to place. This simplifies the representation while does no
harm to the soundness of our proof. Consider
sup
t∈[h̃n ,1−h̃n ]

kÃ(t)−A(t)k ≤

sup
t∈[h̃n ,1−h̃n ]

kÃ(t)−EÃ(t)k+

sup
t∈[h̃n ,1−h̃n ]

kEÃ(t)−A(t)k.

(7.11)
The first term on the right hand side is recognized as the variance and the second
is the bias. Firstly, we deal with the bias term. Denote B(t0 ) := EÃ(t0 ) − A(t0 ),
t0 ∈ [h̃n , 1 − h̃n ]. Recall from (1.2), E(ξj |τj , Xj ) = 0 for any t0 ∈ [h̃n , 1 − h̃n ] we
have
n
τ − t 
m2
m2 X  τj − t0 
0
K
(hA(τj ), Xj i+ξj )Xj =
hA(τ ), XiX.
EK
EÃ(t0 ) = E
nh j=1
h
h
h

By applying the Taylor expansion of A(τ ) as in (7.6) and the fact that K is a
kernel of order ℓ, we get
EÃ(t0 ) = E

m2  τ − t0  (τ − t0 )ℓ
m 2  τ − t0 
hA(t0 ), XiX + E
K
K
hÃ, XiX,
h
h
h
h
ℓ!

where Ã is the same as in (7.6). It is easy to check that the first term on the
right hand side is A(t0 ). Therefore we rewrite B(t0 ) as
B(t0 ) = E

m2  τ − t0  (τ − t0 )ℓ
m2  τ − t0  (τ − t0 )ℓ
K
hÃ, XiX = E
K
hÃ−A(ℓ) (t0 ), XiX,
h
h
ℓ!
h
h
ℓ!

where the second equity is due to the fact that each element of A(t) is in Σ(β, L)
and K is a kernel of order ℓ. Then we can bound each element of matrix B(t0 )
as
Z 1
1  τ − t0  |τ − t0 |ℓ (ℓ)
(ℓ)
|Bij (t0 )| ≤
K
|aij (t0 + α(τ − t0 )) − aij (t0 )|dτ
h
h
ℓ!
0
Z 1
|uh|β
du
|K(u)|
≤L
ℓ!
0
Lhβ
≤ C(K)
.
ℓ!
16

Thus
sup
t∈[h̃n ,1−h̃n ]

kB(t)k ≤ C(K)

Lmhβ
.
ℓ!

(7.12)

On the other hand, for the variance term supt∈[h̃n ,1−h̃n ] kÃ(t) − EÃ(t)k2 , we
construct a δ − net on the interval [0, 1] with δ = 1/M , and
M = n2 ,

tj =

2j − 1
, j = 1, ..., M.
2M

Denote Sn (t) := Ã(t) − EÃ(t), then we have
sup
t∈[h̃n ,1−h̃n ]

kSn (t)k ≤ sup kSn (t)k ≤ max kSn (ti )k + sup kSn (t) − Sn (t′ )k.
i

t∈[0,1]

|t−t′ |≤δ

(7.13)
Next, we bound both terms on the right hand side respectively. For each ti ,
Sn (ti ) =

n
τ − t 

m2 X   τj − ti 
j
i
K
Yj Xj − EK
Yj Xj .
nh j=1
h
h

The right hand side is a sum of zero mean random matrices, we apply the matrix
Bernstein inequality, see [30]. Under the assumption of Theorem 4.2, one can
easily check that with probability at least 1 − e−η ,
r
 a2 (η + log 2m) _ a(η + log 2m) 
2
.
kSn (ti )k ≤ C(K)m
mnh
nh


2
2
Indeed, by setting X̄ = mh K τ h−t Y X − E mh K τ −t
h Y X, it is easy to check
2
that UX̄ . kKk∞am2 /h and σX̄
. RK a2 m3 /h. By taking the union bound
over all i and setting η = 4 log n, we get with probability at least 1 − n−2 ,
2

max Sn (ti )
i

≤ C(K)

a2 m3 log n
,
nh

As for the second term on the right hand side of (7.13), by the assumption
that K is a Lipschitz function with Lipschitz constant LK , we have
sup kSn (t) − Sn (t′ )k ≤

|t−t′ |≤δ

≤

sup k(Ã(t) − Ã(t′ ))k + sup kE(Ã(t) − Ã(t′ ))k

|t−t′ |≤δ

|t−t′ |≤δ

3

LK am
LK am
+ 2 2 .
2
2
n h
n h

Thus with probability at least 1 − n−2 ,
sup
t∈[h̃n ,1−h̃n ]

kSn (t)k2 ≤ C(K)

a2 m3 log n
nh

Together with the upper bound we get on the bias in (7.12), we have with
probability at least 1 − n−2 ,
 a2 m log n L2 h2β 
1
2
.
k
Ã(t)
−
A(t)k
≤
C(K)
+
sup
2
nh
ℓ!2
t∈[h̃n ,1−h̃n ] m
17

Choose

1
 a2 (ℓ!)2 m log n  2β+1
,
h̃n = C(K)
2
2βL n

we get

2β
 a2 (ℓ!)2 m log n  2β+1
1
2
.
k
Ã(t)
−
A(t)k
≤
C(K)
2
2βL2 n
t∈[h̃n ,1−h̃n ] m

sup

7.4. Proof of Theorem 5.1
Proof. Without loss of generality, we assume that both m and r are even numbers. We introduce several notations which are key to construct the hypothesis
set. For some constant γ > 0, denote

m
r
C = Ã = (aij ) ∈ C 2 × 2 : aij ∈ {0, γ}, ∀1 ≤ i ≤ m/2, 1 ≤ j ≤ r/2 ,
and consider the set of block matrices
(
)


m
m
×
B(C) =
Ã Ã . . . Ã O ∈ C 2 2 : Ã ∈ C ,

(7.14)

where O denotes the m/2 × (m/2 − r⌊m/r⌋/2) zero matrix. Then we consider
a subset of Hermitian matrices Sm ⊂ Hm ,
("
#
)
b
Õ A
m×m
b
Sm =
: A ∈ B(C) .
(7.15)
b∗ Õ ∈ C
A
An immediate observation is that for any matrix A ∈ Sm , rank(A) ≤ r.
Due to the Varshamov-Gilbert bound (see Lemma 2.9 in [31]), there exists
a subset A0 ⊂ Sm with cardinality Card(A0 ) ≥ 2mr/32 + 1 containing the zero
m × m matrix 0 such that for any two distinct elements A1 and A2 of A0 ,

m2
mr m 2
⌊ ⌋γ ≥ γ 2
.
(7.16)
16 r
32

0
Let fn (t) denote the function fn (t) := Lhβn f t−t
hn , t ∈ [0, 1], where hn =
1
  2β+1
c0 mr
, with some constant c0 > 0, and f ∈ Σ(β, 1/2)∩C ∞ and Supp(f ) =
n
[−1/2, 1/2]. Note that there exist functions f satisfying this condition. For instance, one can take
− 1
f (t) = αe 1−4u2 I(|u| < 1/2),
(7.17)
kA1 − A2 k22 ≥

for some sufficient small α > 0. It is easy to check that fn (t) ∈ Σ(β, L) on [0, 1].
We consider the following hypotheses of A at t0 :

b = Afn (t), t ∈ [0, 1] : A ∈ A0 .
Aβ0 := A(t)
18

The following claims are easy to check: firstly, any element in Aβ0 together with
its derivative have rank uniformly bounded by r, and the difference of any two
elements of Aβ0 satisfies the same property for fixed t0 ; secondly, the entries
of any element of Aβ0 together with its derivative are uniformly bounded by
some constant for sufficiently small chosen γ; finally, each element of A(t) ∈ Aβ0
belongs to Σ(β, L). Therefore, Aβ0 ⊂ A(r, a) with some chosen γ.
b1 (t) and A
b2 (t) of Aβ , the
According to (7.16), for any two distinct elements A
0
b
b
difference between A1 (t) and A2 (t) at point t0 is given by


2 2 2β 2
b1 (t0 ) − A
b2 (t0 )k2 ≥ γ L c0 f (0) m2 mr 2β+1 .
kA
2
32
n
2β

(7.18)

A
On the other hand, we consider the joint distributions Pτ,X,Y
such that τ ∼
U [0, 1], X ∼ Π0 where Π0 denotes the uniform distribution on X , τ and X are
independent, and
(
1
+ hA(τ ),Xi , Y = a,
PA (Y |τ, X) = 21 hA(τ4a),Xi
, Y = −a.
2 −
4a
A
One can easily check that as long as A(τ ) ∈ Aβ0 , such Pτ,X,Y
belongs to the
distribution class P(r, a). We denote the corresponding n−product probability
measure by PA . Then for any A(τ ) ∈ Aβ0 , the Kullback-Leibler Divergence
between P0 and PA is


1 − p0 (τ, X) 
p0 (τ, X)
+ (1 − p0 (τ, X)) log
,
K(P0 , PA ) = nE p0 (τ, X) log
pA (τ, X)
1 − pA (τ, X)
where pA (τ, X) = 1/2 + hA(τ ), Xi/4a. Note that PA (Y = a|τ, X) ∈ [1/4, 3/4] is
guaranteed provided that |hA(t), Xi| ≤ a. Thus by the inequality − log(1 + u) ≤
−u + u2 /2, ∀u > −1, and the fact that PA (Y = a|τ, X) ∈ [1/4, 3/4], we have
K(P0 , PA ) ≤ nE2(p0 (τ, X) − pA (τ, X))2 ≤

n
EhA(τ ), Xi2 .
8a2

Recall that A(τ ) = Afn (τ ) ∈ Aβ0 , by τ ∼ U [0, 1] and X ∼ Π0 , we have
K(P0 , PA ) ≤

L2 kf k22 c2β+1
γ2
n 1 2
2 2β+1 2 2
0
L
kf
k
h
m
γ
≤
mr.
2 n
2
2
2
8a m
8a

(7.19)

Therefore, provided the fact that Card(A0 ) ≥ 2mr/32 + 1, together with (7.19),
we have
X
1
K(P0 , PA ) ≤ α log(Card(Aβ0 ) − 1)
(7.20)
β
Card(A0 ) − 1
β
A∈A0

is satisfied for any α > 0 if γ is chosen as a sufficiently small constant. In view
of (7.18) and (7.20), the lower bound (5.1) follows from Theorem 2.5 in [31].
19

7.5. Proof of Theorem 5.2
Proof. Without loss of generality, we assume that both m and r are even numbers. Take a real number c1 > 0, define

and

1
l  n  2β+1
m
M = c1
,
mr

φj (t) = Lhβn f

t − t 
j

hn

hn =

,

1
,
2M

tj =

j = 1, ...M,

2j − 1
,
2M

t ∈ [0, 1],

where f is defined the same as in (7.17).
Meanwhile, we consider o
the set of all
n
binary sequences of length M : Ω = ω = (ω1 , ..., ωM ), ωi ∈ {0, 1} = {0, 1}M .

By Varshamov-Gilbert bound, there exists a subset Ω0 = {ω 0 , ..., ω N } of Ω such
M
8
that ω 0 = (0, ..., 0) ∈ Ω0 , and d(ω j , ω k ) ≥ M
8 , ∀ 0 ≤ j < k ≤ N, and N ≥ 2 ,
where d(·, ·) denotes the Hamming distance of two binary sequences. Then we
o
n
M
P
ωj φj (t) : ω ∈ Ω0 .
define a collection of functions based on Ω0 : E = fω (t) =
j=1

From the result of Varshamov-Gilbert bound, we know that S := Card(E) =
M
Card(Ω0 ) ≥ 2 8 + 1. It is also easy to check that for all fω , fω′ ∈ E,
Z

0

1

(fω (t) − fω′ (t))2 dt =

M
X
j=1

(ωj − ωj′ )2

= L2 hn2β+1 kf k22
≥

Z

∆j

M
X
j=1

φ2j (t)dt

(ωj − ωj′ )2

(7.21)

2
L2 h2β
n kf k2 /16,

where ∆j = [(j − 1)/M, j/M ].
In what follows, we combine two fundamental results in coding theory: one
is Varshamov-Gilbert bound ([12, 32]) in its general form of a q-ary code, the
other is the volume estimate of Hamming balls. Let Aq (n, d) denote the largest
size of a q-ary code of block length n with minimal Hamming distance d.
Proposition 7.1. The maximal size of a q − ary code of block length n with
minimal Hamming distance d = pn, satisfies
Aq (n, d + 1) ≥ q n(1−hq (p)) ,

(7.22)

where p ∈ [0, 1 − 1/q], hq (p) = p logq (q − 1) − p logq p − (1 − p) logq (1 − p) is the
q − ary entropy function.
We now have all the elements needed in hand to construct our hypotheses set.
Denote Ω1 = {ω 1 , ..., ω N }, which is a subset of Ω0 without ω 0 . We then consider
o
n
M
P
ωj φj (t) : ω ∈ Ω1 .
a subset E1 of E which is given by E1 := fω (t) =
j=1

20

Clearly, S1 := Card(E1 ) ≥ 2M/8 . Then we define a new collection of matrix
valued functions as

r
m
C = Ã = (aij ) ∈ C 2 × 2 : aij ∈ {δfω : ω ∈ Ω1 , δ ∈ C, ∀1 ≤ i ≤ m/2, 1 ≤ j ≤ r/2}.

Obviously, the collection C is a S1 -ary code of block length mr/4. Thus, we can
apply the result of Proposition 7.1. It is easy to check that for p = 1/4, and
q≥4
1 − hq (p) = 1 − p logq

1
q−1
+ (1 − p) logq (1 − p) ≥ .
p
4

(7.23)

In our case, q = S1 ≥ 2M/8 and n = mr/4. If we take p = 1/4, we know that
mr/16

AS1 (mr/4, mr/16) ≥ AS1 (mr/4, mr/16 + 1) ≥ S1

.

(7.24)

In other words, (7.24) guarantees that there exists a subset H0 ⊂ C with
Card(H0 ) ≥ 2Mmr/128 such that for any A1 , A2 ∈ H0 , the Hamming distance
between A1 and A2 is at least mr/16. Now we define the building blocks of our
hypotheses set
o
n
H := H0 ∪ O m2 × r2 ,
r
where O m2 × r2 is the m
2 × 2 zero matrix. Obviously, H has size Card(H) ≥
2Mmr/64 + 1, and for any A1 (t), A2 (t) ∈ H, the minimum Hamming distance is
still greater than mr/16. We consider the set of matrix valued functions
(
)


B(H) =
Ã Ã . . . Ã O : Ã ∈ H ,

where O denotes the m/2 × (m/2 − r⌊m/r⌋/2) zero matrix. Finally, our hypotheses set of matrix valued functions Hm is defined as
("
#
)
b
Õ A
m×m
b ∈ B(H) .
Hm =
:A
b∗ Õ ∈ C
A

By the definition of Hm and similar to the arguments in proof of Theorem 5.1,
it is easy to check that Hm ⊂ A(r, a), and also
Card(Hm ) ≥ 2Mmr/64 + 1.
Now we consider any two different hypotheses Aj (t), Ak (t) ∈ Hm .
Z 1
Z
mr j m k 1
kAj (t) − Ak (t)k22 dt ≥ γ 2
2
(fω (t) − fω′ (t))2 dt,
16
r
0
0

where ω 6= ω ′ . Based on (7.21), we have
Z 1
2β
 mr  2β+1
2
1
γ 2 L2 h2β
n kf k2
2
.
≥
c
kA
(t)
−
A
(t)k
dt
≥
∗
j
k
2
m2 0
256
n
21

(7.25)

(7.26)

(7.27)

where c∗ is a constant depending on kf k2 , L, c1 and γ.
On the other hand, we repeat the same analysis on the Kullback-Leibler
divergence K(P0 , PA ) as in the proof of Theorem 5.1. One can get
K(P0 , PA ) ≤

M Z
n
n 2X 1 2
γ 2 c2β+1
L2 M mrkf k22
2
1
EhA(τ
),
Xi
≤
γ
,
φ
(τ
)dτ
≤
j
2
2
8a
8a
8a2
j=1 0

(7.28)

where A(τ ) ∈ Hm . Combine (7.25) and (7.28) we know that
1
Card(Hm ) − 1

X

A(t)∈Hm

K(P0 , PA ) ≤ α log(Card(Hm ) − 1)

(7.29)

is satisfied for any α > 0 if γ is chosen as a sufficiently small constant. In view
of (7.27) and (7.29), the lower bound follows from Theorem 2.5 in [31].
7.6. Proof of Theorem 5.3
Proof. Without loss of generality,
assume that m is an even number.
o For some
n
m
constant γ > 0, denote V = v ∈ C 2 : ai ∈ {0, γ}, ∀ 1 ≤ i ≤ m/2 . Due to the

Varshamov-Gilbert bound (see Lemma 2.9 in [31]), there exists a subset V 0 ⊂ V
m
with cardinality Card(V 0 ) ≥ 2m/16 + 1 containing the zero vector 0 ∈ C 2 , and
such that for any two distinct elements v1 and v2 of V 0 ,
kv1 − v2 k22 ≥
Consider the set of matrices
n
v v
B(V) =

...

v



m 2
γ .
16

(7.30)

o
m
m
∈ C 2 × 2 : v ∈ V0 .

Clearly, B(V) is a collection of rank one matrices. Then we construct another
matrix set Vm ,
)
(

Õ V
m×m
∈C
: V ∈ B(V)
Vm =
V ∗ Õ
where Õ is the m/2 × m/2 zero matrix. Apparently, Vm ⊂ Hm .
On the other hand, we define the grid on [0, 1]
l 
M = c2

1
m
 2β+1
n
,
m + log n

and
φj (t) = Lhβn f

t − t 
j
,
hn
22

hn =

1
,
2M

j = 1, ...M,

tj =

2j − 1
,
2M

t ∈ [0, 1]

where
f is defined
n
o the same as in (7.17), and c2 is some constant. Denote Φ :=
b =
φj : j = 1, ...M . We consider the following set of hypotheses: AβB := {A(t)

V φj (t) : V ∈ Vm , φj ∈ Φ}. One can immediately get that the size of AβB satisfies
Card(AβB ) ≥ (2m/16 + 1)M.

(7.31)

b of Aβ has
By construction, the following claims are obvious: any element A(t)
B
β
b
rank at most 2; the entries of A(t)
∈ AB are uniformly bounded for some
bij (t) ∈ Σ(β, L). Thus Aβ ⊂ A(a).
sufficiently small γ, and A
B
Now we consider the distance between two distinct elements A(t) and A′ (t)
of AβB . An immediate observation is that
sup kA(t) − A′ (t)k2 ≥

t∈[0,1]

1
sup kA(t) − A′ (t)k22 ,
4 t∈[0,1]

due to the fact that ∀t ∈ (0, 1), rank(A(t)−A′ (t)) ≤ 4. Then we turn to get lower
bound on sup kA(t) − A′ (t)k22 . Recall that by construction of AβB , we have for
t∈(0,1)

any A 6= A′ , A(t) = A1 φj (t), A′ (t) = A2 φk (t), where A1 , A2 ∈ Vm . There are
three cases need to be considered: 1). A1 6= A2 and j = k; 2). A1 = A2 6= 0 and
j 6= k; 3). A1 6= A2 and j 6= k.
For case 1,
 m + log n  2β+1
m2 2 2 2β
γ L hn kf k2∞ ≥ c∗1 m2
,
16
n
2β

sup kA(t) − A′ (t)k22 = kA1 − A2 k22 kφj k2∞ ≥

t∈[0,1]

where c∗1 is a constant depending on kf k2∞ , β, L and γ.
For case 2,
 m + log n  2β+1
m2 2 2 2β
,
γ L hn kf k2∞ ≥ c∗2 m2
16
n
2β

sup kA(t) − A′ (t)k22 = kA1 k22 kφj − φk k2∞ ≥

t∈[0,1]

where c∗2 is a constant depending on kf k2∞ , β, L and γ.
For case 3,
 m + log n  2β+1
m2 2 2 2β
γ L hn kf k2∞ ≥ c∗3 m2
,
16
n
2β

sup kA(t) − A′ (t)k22 ≥ (kA1 k22 kφj k2∞ ∨ kA2 k22 kφk k2∞ ) ≥

t∈[0,1]

where c∗3 is a constant depending on kf k2∞ , β, L and γ.
Therefore, by the analysis above we conclude that for any two distinct elements A(t) and A′ (t) of AβB ,
 m + log n  2β+1
1
sup kA(t) − A′ (t)k22 ≥ c∗ m2
,
4 t∈[0,1]
n
t∈[0,1]
(7.32)
where c∗ is a constant depending on kf k2∞ , L, γ and β.
23
2β

sup kA(t) − A′ (t)k2 ≥

Meanwhile, we repeat the same analysis on the Kullback-Leibler divergence
K(P0 , PA ) as in the proof of Theorem 5.1. One can get that for any A ∈ AβB ,
the Kullback-Leibler divergence K(P0 , PA ) between P0 and PA satisfies

K(P0 , PA ) ≤

n
n
E|hA(τ ), Xi|2 ≤ 2 γ 2
2
8a
8a

Z

1

0

φ2j (τ )dτ ≤

γ 2 c2β+1
L2 (m + log n)kf k22
2
.
8a2
(7.33)

Combine (7.31) and (7.33) we know that
1
Card(AβB ) − 1

X

A∈Aβ
B

K(P0 , PA ) ≤ α log(Card(AβB ) − 1)

(7.34)

is satisfied for any α > 0 if γ is chosen as a sufficiently small constant. In view
of (7.32) and (7.34), the lower bound follows from Theorem 2.5 in [31].
7.7. Proof of Theorem 6.1
bk , denote the difference in empirical loss between A
bk and A by
Proof. For any A
n

n

n

X
X
1X
bk (τj ), Xj i)2 − 1
bk , A) : = 1
(Yj − hA
(Yj − hA(τj ), Xj i)2 = −
Uj ,
rn (A
n j=1
n j=1
n j=1
bk (τj ), Xj i)2 . It is easy to check that
where Uj = (Yj − hA(τj ), Xj i)2 − (Yj − hA

bk (τj ) − A(τj ), Xj i − hA
bk (τj ) − A(τj ), Xj i2 . (7.35)
Uj = 2(Yj − hA(τj ), Xj i)hA

bk (τ ) − A(τ ), Xi2 . The following concentration inWe denote r(Abk , A) := EhA
equality developed by [10] to prove Bernstein’s inequality is key to our proof.

Lemma 2. Let Uj , j = 1, ..., n be independent bounded
Pn random variables
satisfying |Uj − EUj | ≤ M with h = M/3. Set Ū = n−1 j=1 Uj . Then for all
t>0
n
nεvar(Ū ) o
t
≤ e−t ,
+
P Ū − EŪ ≥
nε
2(1 − c)
with 0 < εh ≤ c < 1.

Firstly, we bound the variance of Uj . Under the assumption that |Y | and
|hA(τ ), Xi| are bounded by a constant a, one can easily check that h = 8a2 /3.
Given E(Yj |τj , Xj ) = hA(τj ), Xj i, we know that the covariance between the
two terms on the right hand side of (7.35) is zero. Conditionally on (τ, X), the
second order moment of the first term satisfies
bk (τj ) − A(τj ), Xj i2 ≤ 4a2 r(Abk , A).
4EσY2 |τ,X hA

To see why, one can consider the random variable Ỹ with the distribution P{Ỹ =
a} = P{Ỹ = −a} = 1/2. The variance of Y is always bounded by the variance
24

bk (τj ), Xj i| are bounded
of Ỹ which is a2 under the assumption that |Yj | and |hA
by a constant a > 0. Similarly, we can get that the variance of the second term
bk (τj )− A(τj ), Xj i2 . As a result,
conditioned on (τ, X) is also bounded by 4a2 EhA
2
k
bk with
b
nvar(Ū ) ≤ 8a r(A , A). By the result of Lemma 2, we have for any A
−t
probability at least 1 − e
2
bk
bk , A) < t + 4a εr(A , A) .
r(Abk , A) − rn (A
nε
1−c

Set t = επk + log 1/δ, we get with probability at least 1 − δ/eεπk

 log 1/δ 
2
bk , A) + πk + 4a
.
(1 − α)r(Abk , A) < rn (A
n
(1 − c)α
n

where α = 4a2 ε/(1 − c) < 1. Denote

n
πk o
.
k̃ ∗ = arg min r(Abk , A) +
k
n

b∗ , we have with probability at least 1 − δ/eεbπ∗
By the definition of A
 log 1/δ 
2
bk̃∗ , A) + πk̃∗ + 4a
(1 − α)r(Ab∗ , A) < rn (A
.
n
(1 − c)α
n

(7.36)

b∗ .
where π
b∗ is the penalty terms associated with A
Now we apply the result of Lemma 2 one more time and set t = log 1/δ, we
get with probability at least 1 − δ
bk̃∗ , A) ≤ (1 + α)r(Abk̃∗ , A) +
rn (A

4a2 log 1/δ
.
(1 − c)α n

(7.37)

Apply the union bound of (7.36) and (7.37), we get with probability at least
∗
1 − δ(1 + e−εbπ )
b∗ , A) ≤
r(A

π ∗
4a2
log 1/δ
(1 + α)  bk̃∗
r(A , A) + k̃ +
.
(1 − α)
n
(1 − c)α(1 − α) n

By taking ε = 3/32a2 and c = εh,



2
bk̃∗ , A) + πk̃∗ + 64a log 1/δ .
r(Ab∗ , A) ≤ 3 r(A
n
3
n

By taking δ = 1/nmr and adjusting the constant, we have with probability at
least 1 − 1/nmr
1
m2

Z

0

1

b∗ (t)−A(t)
A

2
dt
2

≤ 3 min
k

n 1 Z 1
bk (t)−A(t)
A
m2 0

where C(a) is a constant depending on a.
25

πk
2
dt+
2
n

o
mr log n
+C(a)
n

Acknowledgements
The author would like to thank Dr. Vladimir Koltchinskii for all the guidance,
discussion and inspiration along the course of writing this paper. The author also
would like to thank NSF for its generous support through Grants DMS-1509739
and CCF-1523768.
References
[1] Rudolf Ahlswede and Andreas Winter. Strong converse for identification
via quantum channels. IEEE Transactions on Information Theory, 48(3):
569–579, 2002.
[2] Jean-Pierre Aubin and Ivar Ekeland. Applied nonlinear analysis. Courier
Corporation, 2006.
[3] Andrew R Barron. Complexity regularization with application to artificial
neural networks. Nonparametric Functional Estimation and Related Topics,
335:561–576, 1991.
[4] Olivier Bousquet. A bennett concentration inequality and its application
to suprema of empirical processes. Comptes Rendus Mathematique, 334(6):
495–500, 2002.
[5] Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen. A singular value
thresholding algorithm for matrix completion. SIAM Journal on Optimization, 20(4):1956–1982, 2010.
[6] Emmanuel J Candes and Yaniv Plan. Matrix completion with noise. Proceedings of the IEEE, 98(6):925–936, 2010.
[7] Emmanuel J Candès and Benjamin Recht. Exact matrix completion via
convex optimization. Foundations of Computational mathematics, 9(6):717,
2009.
[8] Emmanuel J Candès and Terence Tao. The power of convex relaxation:
Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5):2053–2080, 2010.
[9] Sourav Chatterjee. Matrix estimation by universal singular value thresholding. The Annals of Statistics, 43(1):177–214, 2015.
[10] Cecil C Craig. On the Tchebychef inequality of Bernstein. The Annals of
Mathematical Statistics, 4(2):94–102, 1933.
[11] Jianqing Fan and Irene Gijbels. Local polynomial modelling and its applications: monographs on statistics and applied probability 66, volume 66. CRC
Press, 1996.
[12] Edgar N Gilbert. A comparison of signalling alphabets. Bell Labs Technical
Journal, 31(3):504–522, 1952.
[13] David Gross. Recovering low-rank matrices from few coefficients in any
basis. IEEE Transactions on Information Theory, 57(3):1548–1566, 2011.
[14] Vladimir Koltchinskii. Local rademacher complexities and oracle inequalities in risk minimization. The Annals of Statistics, 34(6):2593–2656, 2006.
[15] Vladimir Koltchinskii. Oracle inequalities in empirical risk minimization
and sparse recovery problems. 2011.
26

[16] Vladimir Koltchinskii. Von Neumann entropy penalization and low-rank
matrix estimation. The Annals of Statistics, pages 2936–2973, 2011.
[17] Vladimir Koltchinskii. Sharp oracle inequalities in low rank estimation. In
Empirical Inference, pages 217–230. Springer, 2013.
[18] Vladimir Koltchinskii, Karim Lounici, and Alexandre B Tsybakov. Nuclearnorm penalization and optimal rates for noisy low-rank matrix completion.
The Annals of Statistics, 39(5):2302–2329, 2011.
[19] Yehuda Koren. Collaborative filtering with temporal dynamics. Communications of the ACM, 53(4):89–97, 2010.
[20] Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8), 2009.
[21] Oleg V Lepski and Vladimir G Spokoiny. Optimal pointwise adaptive methods in nonparametric estimation. The Annals of Statistics, pages 2512–
2546, 1997.
[22] Oleg V Lepski, Enno Mammen, and Vladimir G Spokoiny. Optimal spatial
adaptation to inhomogeneous smoothness: an approach based on kernel
estimates with variable bandwidth selectors. The Annals of Statistics, pages
929–947, 1997.
[23] Oleg V Lepskii. On a problem of adaptive estimation in gaussian white
noise. Theory of Probability & Its Applications, 35(3):454–466, 1991.
[24] Elliott H Lieb. Convex trace functions and the wigner-yanase-dyson conjecture. Advances in Mathematics, 11(3):267–288, 1973.
[25] Sahand Negahban and Martin J Wainwright. Restricted strong convexity
and weighted matrix completion: Optimal bounds with noise. Journal of
Machine Learning Research, 13(May):1665–1697, 2012.
[26] Benjamin Recht, Maryam Fazel, and Pablo A Parrilo. Guaranteed
minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 52(3):471–501, 2010.
[27] Angelika Rohde and Alexandre B Tsybakov.
Estimation of highdimensional low-rank matrices. The Annals of Statistics, 39(2):887–930,
2011.
[28] Amit Singer and Mihai Cucuringu. Uniqueness of low-rank matrix completion by rigidity theory. SIAM Journal on Matrix Analysis and Applications,
31(4):1621–1641, 2010.
[29] Michel Talagrand. New concentration inequalities in product spaces. Inventiones mathematicae, 126(3):505–563, 1996.
[30] Joel A Tropp. User-friendly tail bounds for sums of random matrices.
Foundations of computational mathematics, 12(4):389–434, 2012.
[31] Alexandre B Tsybakov. Introduction to nonparametric estimation. revised
and extended from the 2004 french original. translated by vladimir zaiats,
2009.
[32] Rom R Varshamov. Estimate of the number of signals in error correcting
codes. In Dokl. Akad. Nauk SSSR, volume 117, pages 739–741, 1957.
[33] Marten Wegkamp. Model selection in nonparametric regression. The Annals of Statistics, 31(1):252–273, 2003.
[34] Luwan Zhang, Grace Wahba, and Ming Yuan. Distance shrinkage and
27

euclidean embedding via regularized kernel estimation. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 78(4):849–867, 2016.

8. Appendix: Proof of Lemma 1
The proof of Lemma 1 follows from a similar approach introduced by [17].
Pr
Proof. For any S ∈ Hm of rank r, S = j=1 λi (ej ⊗ ej ), where λj are nonzero eigenvalues of S (repeated with their multiplicities) and ej ∈ Cm are the
r
P
sign(λi )(ej ⊗ej ).
corresponding orthonormal eigenvectors. Denote sign(S) :=
j=1

Let PL , PL⊥ be the following orthogonal projectors in the space (Hm , h·, ·i):
PL (A) := A − PL⊥ APL⊥ , PL⊥ (A) := PL⊥ APL⊥ , ∀A ∈ Hm

where PL denotes the orthogonal projector on the linear span of {e1 , ..., er },
and PL⊥ is its orthogonal complement. Clearly, this formulation provides a
decomposition of a matrix A into a ”low rank part” PL (A) and a ”high rank
part” PL⊥ (A) if rank(S) = r is small. Given b > 0, define the following cone in
the space Hm :
K(D; L; b) := {A ∈ D : kPL⊥ Ak1 ≤ bkPL (A)k1 }
which consists of matrices with a ”dominant” low rank part if S is low rank.
Firstly, we can rewrite (3.1) as
E2
D
1 X
Sbh = arg min
+ εkSk1.
Ỹj − S, X̃j
S∈D n
j=1
n

(8.1)

r
r
 

 

 
 i



hr
τj −t0
τj −t0
τ −t
τj −t0
τj −t0
τj −t0
1
1
1
p
X
,
p
X
,
...,
pℓ j h 0 Xj ,
K
K
where X̃j = Diag
0
j
1
j
h
h
h
h
h
h
hK
h
r


τ −t
and Ỹj = h1 K j h 0 Yj .
Denote the loss function as

E2

D

L Ỹ ; hS(τ ), X̃i := Ỹj − S, X̃j
,

and the risk



1  τ − t0 
P L Ỹ ; hS(τ ), X̃i := EL Ỹ ; hS(τ ), X̃i = σ 2 + E K
Y − hS(τ ), Xi)2
h
h

Since Sbh is a solution of the convex optimization problem (8.1), there exists
a Vb ∈ ∂ Sbh 1 , such that for ∀S ∈ D (see [2] Chap. 2)

E
2 X D bh
S , X̃j − Ỹj hSbh − S, X̃j i + εhVb , Sbh − Si ≤ 0.
n j=1
n

28

This implies that, for all S ∈ D,

b X̃i)hSbh − S, X̃i + εhVb , Sbh − Si
EL′ (Ỹ ; hS,
n
2 X bh
(hS , X̃j i − Ỹj )hSbh − S, X̃j i.
≤ EL′ (Ỹ ; hSbh , X̃i)hSbh − S, X̃i −
n j=1

(8.2)

where L′ denotes the partial derivative of L(y; u) with respect to u. One can
easily check that for ∀S ∈ D,

EL′ (Ỹ ; hSbh , X̃i)hSbh − S, X̃i ≥ E(L(Ỹ ; hSbh , X̃i) − L(Ỹ ; hS, X̃i)) + kSbh − Sk2L2 (Π̃) .
(8.3)
where Π̃ denotes the distribution of X̃. If EL(Ỹ ; hSbh , X̃i) ≤ EL(Ỹ ; hS, X̃i) for
∀S ∈ D, then the oracle inequality in Lemma 1 holds trivially. So we assume
that EL(Ỹ ; hSbh , X̃i) > EL(Ỹ ; hS, X̃i) for some S ∈ D. Thus, inequalities (8.2)
and (8.3) imply that
EL(Ỹ ; hSbh , X̃i) + kSbh − Sk2L2 (Π̃) + εhVb , Sbh − Si

≤ EL(Ỹ ; hS, X̃i) + EL′ (Ỹ ; hSbh , X̃i)hSbh − S, X̃i −

n

2 X bh
(hS , X̃j i − Ỹj )hSbh − S, X̃j i.
n j=1

(8.4)
According to the well known representation of subdifferential of nuclear norm,
see [15] Sec. A.4, for any V ∈ ∂kSk1 , we have
V := sign(S) + PL⊥ (W ), W ∈ Hm , kW k ≤ 1.

By the duality between nuclear norm and operator norm
hPL⊥ (W ), Sbh − Si = hPL⊥ (W ), Sbh i = hW, PL⊥ (Sbh )i = kPL⊥ (Sbh )k1 .

Therefore, by the monotonicity of subdifferentials of convex function k · k1 , for
any V := sign(S) + PL⊥ (W ) ∈ ∂kSk1 , we have
hV, Sbh − Si = hsign(S), Sbh − Si + kPL⊥ (Sbh )k1 ≤ hVb , Sbh − Si,

(8.5)

we can use (8.5) to change the bound in (8.4) to get
EL(Ỹ ; hSbh , X̃i) + kS − Sbh k2L2 (Π̃) + εkPL⊥ (Sbh )k1

≤ EL(Ỹ ; hS, X̃i) + εhsign(S), S − Sbh i + EL′ (Ỹ ; hSbh , X̃i)hSbh − S, X̃i
n
2 X bh
(hS , X̃j i − Ỹj )hSbh − S, X̃j i.
−
n j=1

(8.6)

For the simplicity of representation, we use the following notation to denote the
empirical process:
(P − Pn )(L′ (Ỹ ; hSbh , X̃i))hSbh − S, X̃i :=
n
2 X bh
EL′ (Ỹ ; hSbh , X̃i)hSbh − S, X̃i −
(hS , X̃j i − Ỹj )hSbh − S, X̃j i.
n j=1
29

(8.7)

The following part of the proof is to derive an upper bound on the empirical
process (8.7). Before we start with the derivation, let us present several vital
ingredients that will be used in the following literature. For a given S ∈ D and
for δ1 , δ2 , δ3 , δ4 ≥ 0, denote
A(δ1 , δ2 ) := {A ∈ D : A − S ∈ K(D; L; b), kA − SkL2 (Π̃) ≤ δ1 , kPL⊥ Ak1 ≤ δ2 },

Ã(δ1 , δ2 , δ3 ) := {A ∈ D : kA − SkL2 (Π̃) ≤ δ1 , kPL⊥ Ak1 ≤ δ2 , kPL (A − S)k1 ≤ δ3 },
and

Ǎ(δ1 , δ4 ) := {A ∈ D : kA − SkL2 (Π̃) ≤ δ1 , kA − Sk1 ≤ δ4 },
αn (δ1 , δ2 ) := sup{|(P − Pn )(L′ (Ỹ ; hA, X̃i))hA − S, X̃i| : A ∈ A(δ1 , δ2 )},

α̃n (δ1 , δ2 , δ3 ) := sup{|(P − Pn )(L′ (Ỹ ; hA, X̃i))hA − S, X̃i| : A ∈ Ã(δ1 , δ2 , δ3 )},
α̌n (δ1 , δ4 ) := sup{|(P − Pn )(L′ (Ỹ ; hA, X̃i))hA − S, X̃i| : A ∈ Ǎ(δ1 , δ4 )}.

Given the definitions above, Lemma 3 below shows upper bounds on the three
quantities αn (δ1 , δ2 ), α̃n (δ1 , δ2 , δ3 ), α̌n (δ1 , δ4 ). The proof of Lemma 3 can be
found in section 8.1. Denote
n
X
(8.8)
εj X̃j
Ξ := n−1
j=1

where εj are i.i.d. Rademacher random variables.
Lemma 3. Suppose 0 < δk− < δk+ , k = 1, 2, 3, 4. Let η > 0 and
η̄ := η +

2
X

log([log2 (

δk+
)] + 2) + log 3,
δk−

3
X

log([log2 (

δk+
)] + 2) + log 3,
δk−

k=1

η̃ := η +

k=1

η̌ := η +

X

k=1,k=4

log([log2 (

δk+
)] + 2) + log 3.
δk−

Then with probability at least 1 − e−η , for all δk ∈ [δk− , δk+ ], k=1,2,3

p
2(ℓ + 1)R(T )Φaη̄
C1 (ℓ + 1)R(T )Φa n
√
√
EkΞk( rank(S)mδ1 +δ2 )+
+δ1
h
n h
(8.9)
r
2(ℓ + 1)R(T )Φaη̃
η̃ o
C2 (ℓ + 1)R(T )Φa n
√
√
EkΞk(δ2 + δ3 ) +
α̃n (δ1 , δ2 , δ3 ) ≤
+ δ1
n
h
n h
(8.10)
r o
η̌
C3 (ℓ + 1)R(T )Φa n
2(ℓ + 1)R(T )Φaη̌
√
√
(8.11)
+ δ1
α̌n (δ1 , δ4 ) ≤
EkΞkδ4 +
n
h
n h
where C1 , C2 , and C3 are numerical constants.
30
αn (δ1 , δ2 ) ≤

r o
η̄
n

Since both Sbh and S are in D, by the definition of α̃ and α̌, we have

(P − Pn )(L′ (Ỹ ; hSbh , X̃i))hSbh − S, X̃i ≤ α̃(kSbh − SkL2 (Π̃) ; kPL⊥ Sbh k1 ; kPL (Sbh − S)k1 ),
(8.12)
and
(P − Pn )(L′ (Ỹ ; hSbh , X̃i))hSbh − S, X̃i ≤ α̌(kSbh − SkL2 (Π̃) ; kSbh − Sk1 ), (8.13)

If Sbh − S ∈ K(D; L; b), by the definition of α, we have

(P − Pn )(L′ (Ỹ ; hSbh , X̃i))hSbh − S, X̃i ≤ α(kSbh − SkL2 (Π̃) ; kPL⊥ Sbh k1 ), (8.14)

Assume for a while that

kSbh − SkL2 (Π̃) ∈ [δ1− , δ1+ ], kPL⊥ Sbh k1 ∈ [δ2− , δ2+ ], kPL⊥ (Sbh − S)k1 ∈ [δ3− , δ3+ ].
(8.15)
By the definition of subdifferential, for any Vb ∈ ∂kSbhk1 ,
hVb , S − Sbh i ≤ kSk1 − kSbh k1 .

Then we apply (8.13) in bound (8.4) and use the upper bound on α̌n (δ1 , δ4 ) of
Lemma 3, and get with probability at least 1 − e−η ,
P (L(Ỹ ; hSbh , X̃i)) + kSbh − Sk2L2 (Π̃)

≤ P (L(Ỹ ; hS, X̃i)) + ε(kSk1 − kSbh k1 ) + α̌n (kSbh − SkL2 (Π̃) , kSbh − Sk1 )

≤ P (L(Ỹ ; hS, X̃i)) + ε(kSk1 − kSbh k1 )
r o
η̌
C3 (ℓ + 1)R(T )Φa n
2(ℓ
+
1)R(T
)Φaη̌
h
h
√
√
.
+
EkΞkkSb − Sk1 +
+ kSb − SkL2 (Π̃)
n
h
n h
(8.16)
Assuming that
C(ℓ + 1)R(T )Φa
√
ε>
EkΞk,
(8.17)
h
where C = C1 ∨ 4C2 ∨ C3 . From (8.16)
C3 (ℓ + 1)2 R(T )2 Φ2 a2 η̃
.
P (L(Ỹ ; hSbh , X̃i)) ≤ P (L(Ỹ ; hS, X̃i)) + 2εkSk1 +
nh
(8.18)
We now apply the upper bound on α̃n (kSbh − SkL2 (Π̃) , kPL⊥ Sbh )k1 , kPL (Sbh −
S)k1 ) to (8.6) and get with probability at least 1 − e−η ,
P (L(Ỹ ; hSbh , X̃i)) + kSbh − Sk2L2 (Π̃) + εkPL⊥ (Sbh )k1

≤ P (L(Ỹ ; hS, X̃i)) + εkPL (Sbh − S)k1 + α̃n (kSbh − SkL2 (Π̃) , kPL⊥ Sbh )k1 , kPL (Sbh − S)k1 )

≤ P (L(Ỹ ; hS, X̃i)) + εkPL (Sbh − S)k1
o C (ℓ + 1)2 R(T )2 Φ2 a2 η̃
C2 (ℓ + 1)R(T )Φa n
2
√
,
EkΞk(kPL⊥ Sbh )k1 + kPL (Sbh − S)k1 ) +
+
nh
h
(8.19)
31

where the first inequality is due to the fact that
|hsign(S), S−Sbh i| = |hsign(S), PL (S−Sbh )i| ≤ ksign(S)kkPL (S−Sbh )k1 ≤ kPL (S−Sbh )k1 .

With assumption (8.17) holds, we get from (8.19)
P L(Ỹ ; hSbh , X̃i) + εkPL⊥ (Sbh )k1

5ε
ε
C2 (ℓ + 1)2 R(T )2 Φ2 a2 η̃
kPL (Sbh − S)k1 + kPL⊥ (Sbh )k1 +
.
4
4
nh
(8.20)
If the following is satisfied:
≤ P L(Ỹ ; hS, X̃i) +

5ε
ε
C2 (ℓ + 1)2 R(T )2 Φ2 a2 η̃
≥ kPL (Sbh − S)k1 + kPL⊥ (Sbh )k1 ,
nh
4
4

(8.21)

we can just conclude that

C2 (ℓ + 1)2 R(T )2 Φ2 a2 η̃
P (L(Ỹ ; hSbh , X̃i)) ≤ P (L(Ỹ ; hS, X̃i)) +
,
nh

(8.22)

which is sufficient to meet the bound of Lemma 1. Otherwise, by the assumption
that P (L(Ỹ ; hSbh , X̃i)) > P (L(Ỹ ; hS, X̃i)), one can easily check that
kPL⊥ (Sbh − S)k1 ≤ 5kPL (Sbh − S)k1 ,

which implies that Sbh − S ∈ K(D; L; 5). This fact allows us to use the bound on
αn (δ1 , δ2 ) of Lemma 3. We get from (8.6)
P (L(Ỹ ; hSbh , X̃i)) + kSbh − Sk2L2 (Π̃) + εkPL⊥ (Sbh )k1

≤ P (L(Ỹ ; hS, X̃i)) + εhsign(S), S − Sbh i
p
C1 (ℓ + 1)R(T )Φa
C1 (ℓ + 1)2 R(T )2 Φ2 a2 η̄
√
+
EkΞk( rank(S)mkSbh − SkL2 (Π̃) + kPL⊥ (Sbh )k1 ) +
.
nh
h
(8.23)
By applying the inequality
p
hsign(S), Sbh − Si ≤ m rank(S)kSbh − SkL2 (Π̃) ,

and the assumption (8.17), we have with probability at least 1 − e−η ,

C1 (ℓ + 1)2 R(T )2 Φ2 a2 η̄
.
nh
(8.24)
To sum up, the bound of Lemma 1 follows from (8.18), (8.22) and (8.24)
provided that condition (8.17) and condition (8.15) hold.
We still need to specify δk− , δk+ , k = 1, 2, 3, 4 to establish the bound of the
theorem. By the definition of Sbh , we have
P (L(Ỹ ; hSbh , X̃i)) ≤ P (L(Ỹ ; hS, X̃i)) + ε2 m2 rank(S) +

Pn (L(Ỹ ; hX, Sbh i)) + εkSbh k1 ≤ Pn (L(Ỹ ; hX, 0i)) ≤ Q,
32

Q
⊥ bh
bh
bh
implying that kSbh k1 ≤ Q
ε . Next, kPL S k1 ≤ kS k1 ≤ ε and kPL (S − S)k1 ≤
2Q
h
h
b
b
2kS − Sk1 ≤ ε + 2kSk1. Finally, we have kS − SkL2 (Π̃) ≤ 2a. Thus, we can
2Q
Q
+
+
take δ1+ := 2a, δ2+ := Q
ε , δ3 := ε + 2kSk1 , δ4 := ε + kSk1 . With these choices,
+
δk , k = 1, 2, 3, 4 are upper bounds on the corresponding norms in condition
2

δ+

2

δ+

2

δ+

a
a
a
(8.15). We choose δ1− := √an , δ2− := nε
∧ 22 , δ3− := nε
∧ 23 , δ4− := nε
∧ 24 .
Let η ∗ := η + 3 log(B log2 (kSk1 ∨ n ∨ ε ∨ a−1 ∨ Q)). It is easy to verify that
η̄ ∨ η̃ ∨ η̃ ≤ η ∗ . for a proper choice of numerical constant B in the definition
of η ∗ . When condition (8.15) does not hold, which means at least one of the
numbers δk− , k = 1, 2, 3, 4 we chose is not a lower bound on the corresponding
norm, we can still use the bounds

(P − Pn )(L′ (Y ; hSbh , X̃i))hSbh − S, X̃i
≤ α̃(kSbh − Sk
∨ δ − ; kPL⊥ Sbh k1 ∨ δ − ; kPL (Sbh − S)k1 ∨ δ − ),
L2 (Π̃)

and

1

2

(8.25)

3

−
−
b ε − Sk
bh
(P − Pn )(L′ (Y ; hSbh , X̃i))hSbh − S, X̃i ≤ α̌(kA(t)
L2 (Π̃) ∨ δ1 ; kS − Sk1 ∨ δ4 ),
(8.26)
instead of (8.12), (8.13). In the case when Sbh − S ∈ K(D; L; 5), we can use the
bound

(P − Pn )(L′ (Y ; hSbh , X̃i))hSbh − S, X̃i ≤ α(kSbh − SkL2 (Π̃) ∨ δ1− ; kPL⊥ Sbh k1 ∨ δ2− ),
(8.27)
instead of bound (8.14). Then one can repeat the arguments above with only
minor modifications. By the adjusting the constants, the result of Lemma 1
holds.
The last thing we need to specify is the size of ε which controls the nuclear
norm penalty. Recall that from condition (8.17), the essence is to control EkΞk.
Here we use a simple but powerful noncommutative matrix Bernstein inequalities. The original approach was introduced by [1]. Later, the result was improved
by [30] based on the classical result of [24]. We give the following lemma which
is a direct consequence of the result proved by [30], and we omit the proof here.
Lemma 4. Under the model (1.1), Ξ is defined as in (8.8) with τj are i.i.d.
uniformly distributed in [0,1], and εj are i.i.d. Rademacher random variables,
and Xj are i.i.d uniformly distributed in X . Then for any η > 0, with probability
at least 1 − e−η
r
 (η + log 2m) _ (η + log 2m)Φ 
√
,
kΞk ≤ 4
nm
n h
and by integrating its exponential tail bounds
r log 2m _ (log 2m)Φ 
√
EkΞk ≤ C
nm
n h
where C is a numerical constant.

33

Together with (8.17), we know for some numerical constant D > 0,
r
Φa(ℓ + 1)R(T )  log 2m _ (log 2m)Φ 
√
√
ε≥D
.
nm
h
n h
which completes the proof of Lemma 1.
8.1. Proof of Lemma 3
Proof. We only prove the first bound in detail, and proofs of the rest two bounds
are similar with only minor modifications. By Talagrand’s concentration inequality [29], and its Bousquet’s form [4], with probability at least 1 − e−η ,
√
24(ℓ + 1)2 R(T )2 Φ2 a2 η 12(ℓ + 1)R(T )Φaδ1 η
√
+
αn (δ1 , δ2 ) ≤ 2Eαn (δ1 , δ2 ) +
.
nh
nh
(8.28)
By standard Rademacher symmetrization inequalities, see [15], Sec. 2.1, we can
get
n
n 1X
o
Eαn (δ1 , δ2 ) ≤ 4E sup
εj (hA, X̃j i − Ỹj )hA − S, X̃j i : A ∈ A(δ1 , δ2 ) ,
n j=1

(8.29)
where {εj } are i.i.d. Rademacher random variables independent of {(τj , Xj , Ỹj )}.
√ and |v|, |u| ≤
Then we consider the function f (u) = (u−y+v)u, where |y| ≤ 2Φa
h
)Φa
√
Clearly, this function has a Lipschitz constant 6(ℓ+1)R(T
. Thus
h
by comparison inequality, see [15], Sec. 2.2, we can get
n
o
n
X
εj (hA, X̃j i − Ỹj )hA − S, X̃j i : A ∈ A(δ1 , δ2 )
E sup n−1
2(ℓ+1)R(T )Φa
√
.
h

j=1

n
o
n
X
6(ℓ + 1)R(T )Φa
√
≤
εi hA − S, X̃j i : A ∈ A(δ1 , δ2 ) .
E sup n−1
h
j=1

(8.30)

As a consequence of (8.29) and (8.30), we have
n
o
n
X
12(ℓ + 1)R(T )Φa
−1
√
εi hA − S, X̃j i : A ∈ A(δ1 , δ2 ) .
E sup n
Eαn (δ1 , δ2 ) ≤
h
j=1
(8.31)
n
P
εi hA − S, X̃j i . Recall that
The next step is to get an upper bound on n−1
j=1

Ξ := n

−1

n
P

j=1

check that

εj X̃j , then we have n

−1

n
P

j=1

εi hA − S, X̃j i = hA − S, Ξi. One can

|hA − S, Ξi| ≤ |hPL (A − S), PL Ξi| + |hPL⊥ (A − S), Ξi|

≤ kPL Ξk2 kPL (A − S)k2 + kΞkkPL⊥ Ak1
p
≤ m 2rank(S)kΞkkA − SkL2 (Π̃) + kΞkkPL⊥ Ak1 .
34

The second line of this inequality is due to Hölder’s inequality and the third
line is due to the facts that (A − S) ∈ K(D; L; 5), rank(PL (Ξ)) ≤ 2rank(S),
p
2
kPL Ξk2 ≤ 2 rank(PL (Ξ))kΞk, and A − S L2 (Π̃) = m12 A − S 2 . Therefore,
n
o
n
X
12(ℓ + 1)R(T )Φa
√
εi hA − S, X̃j i : A ∈ A(δ1 , δ2 )
E sup n−1
h
j=1

p
12(ℓ + 1)R(T )Φa
√
≤
EkΞk(2 2rank(S)mδ1 + δ2 ).
h

(8.32)

It follows from (8.28), (8.31) and (8.32) that with probability at least 1 − e−η ,
 12(ℓ + 1)R(T )Φa

p
√
EkΞk( rank(S)mδ1 + δ2 )
h
√
2
12(ℓ + 1)R(T )Φaδ1 η
24(ℓ + 1) R(T )2 Φ2 a2 η
√
+
.
+
nh
nh

αn (δ1 , δ2 ) ≤

Now similar to the approach in [17], we make this bound uniform in δk ∈ [δk− , δk+ ].
Let δkjk = δk+ 2−jk , jk = 0, ..., [log2 (δk+ /δk− )] + 1, k = 1, 2. By the union bound,
with probability at least 1 − e−η /3, for all jk = 0, ..., [log2 (δk+ /δk− )] + 1, k = 1, 2,
we have
 12(ℓ + 1)R(T )Φa

p
√
EkΞk( rank(S)mδ1j1 + δ2j2 )
αn (δ1 , δ2 ) ≤
h
√
2
12(ℓ + 1)R(T )Φaδ1j1 η
24(ℓ + 1) R(T )2 Φ2 a2 η
√
+
.
+
nh
nh
which implies that for all δk ∈ [δk− , δk+ ], k = 1, 2,
 12(ℓ + 1)R(T )Φa

p
√
EkΞk( rank(S)mδ1 + δ2 )
h
√
24(ℓ + 1)2 R(T )2 Φ2 a2 η̄
12(ℓ + 1)R(T )Φaδ1 η̄
√
.
+
+
nh
nh

αn (δ1 , δ2 ) ≤

The proofs of the second and the third bounds are similar to this one, we omit
the repeated arguments.

35

