Poisson Matrix Completion

0.88 MB·English

by  Yang Cao

Poisson Matrix Completion
Yang Cao Yao Xie
H. Milton Stewart School of H. Milton Stewart School of
Industrial and Systems Engineering Industrial and Systems Engineering
Georgia Institute of Technology Georgia Institute of Technology
[email protected] [email protected]
Abstract—We extend the theory of matrix completion to the establish a lower bound. Another related work [22] (not in the
5 casewherewemakePoissonobservationsforasubsetofentriesof matrix completion setting) considers the case where all entries
1 a low-rank matrix. We consider the (now) usual matrix recovery of the low-rank matrix are observed and the observations are
0 formulationthroughmaximumlikelihoodwithproperconstraints Poisson counts of the entries of the underlying matrix. In the
2 on the matrix M, and establish theoretical upper and lower
compressed sensing literature, there is a line of research for
bounds on the recovery error. Our bounds are nearly optimal
r sparse signal recovery in the presence of Poisson noise [23]–
a up to a factor on the order of O(log(d d )). These bounds are
1 2 [25] and the corresponding performance bounds. The recently
M obtained by adapting the arguments used for one-bit matrix
developed SCOPT [26], [27] algorithm can also be used to
completion [1] (although these two problems are different in
5 nature) and the adaptation requires new techniques exploiting solve the Poisson compressed sensing problems.
2 properties of the Poisson likelihood function and tackling the
difficulties posed by the locally sub-Gaussian characteristic of In this paper, we extend the theory of matrix completion

] the Poisson distribution. Our results highlight a few important to the case of Poisson observations. We study recovery based
L distinctions of Poisson matrix completion compared to the prior
on maximum likelihood with proper constraints on a matrix
M workinmatrixcompletionincludinghavingtoimposeaminimum
signal-to-noise requirement on each observed entry. We also M with ra√nk less than or equal to r (nuclear norm bound
t. develop an efficient iterative algorithm and demonstrate its good (cid:107)M(cid:107)∗ ≤ α rd1d2 for some constant α and bounded entries
a performance in recovering solar flare images. β ≤ Mij ≤ α). Note that the formulation differs from the
t one-bit matrix completion case in that we also require a lower
s Keywords—matrix completion, Poisson noise, high-dimensional
[ bound on each entry of the matrix. This is consistent with
statistics, information theory
an intuition that the value of each entry can be viewed as the
6
signal-to-noiseratio(SNR)foraPoissonobservation,andhence
v
I. INTRODUCTION this essentially poses a requirement for the minimum SNR. We
3
also establish upper and lower bounds on the recovery error,
4 Matrix completion, with a goal of recovering a low-rank
2 matrixM ∈Rd1×d2 fromobservationsofasubsetofitsentries, by adapting the arguments used for one-bit matrix completion
6 [1]. The upper and lower bounds nearly match up to a factor
attracts much interests recently due to its important real world
0 applications including the famous Netflix problem [2]. Poisson on the order of O(log(d1d2)), which shows that the convex
. relaxation formulation for Poisson matrix completion is nearly
1 matrix completion, where the observations are Poisson counts
optimal. (We conjecture that such a gap is inherent to the
0 of a subset of the entries, is an important instance in its own
5 as it occurs from a myriads of applications including optical Poisson problem). Moreover, we also highlight a few important
1 distinctions of Poisson matrix completion compared to the
imaging, nuclear medicine, low-dose x-ray imaging [3], and
: prior work on matrix completion in the absence of noise and
v network traffic analysis [4].
with Gaussian noise: (1) Although our arguments are adapted
i
X Recently, much success has been achieved in solving the from one-bit matrix completion (where the upper and lower
matrix completion problem using nuclear norm minimization, bounds nearly match), in the Poisson case there will be a
r
a partly inspired by the theory of compressed sensing [5]. It has gap between the upper and lower bounds, possibly due to the
been shown that when M is low rank, it can be recovered fact that Poisson distribution is only locally sub-Gaussian. In
from only a few observations on its entries (see, e.g. [6]–[14]). our proof, we notice that the arguments based on bounding
Earlier work on matrix completion typically assume that the all moments of the observations, which usually generate tight
observationsarenoiseless,i.e.,wemaydirectlyobserveasubset bounds for prior results with sub-Gaussian observations, do not
of entries of M. In the real world, however, the observations generate tight bounds here; (2) We will need a lower bound
are noisy, which is the focus of the subsequent work [15]– on each matrix entry in the maximum likelihood formulation,
[20], most of which consider a scenario where M is the sum which can be viewed as a requirement for the lowest signal-
of a low-rank matrix with a Gaussian random matrix, i.e., to-noise ratio (since the signal-to-no√ise ratio (SNR) of a
the observations are a subset of entries of M contaminated PoissonobservationwithintensityI is I).Comparedwiththe
with Gaussian noise. Recently there has also been work which more general framework for M-estimator [28], our results are
consider the more general noise models, including noisy 1-bit specific to the Poisson case, which may possible be stronger
observations [1], which may be viewed as a case where the but do not apply generally. We also develop several simple
observations are Bernoulli random variables whose parameters yet efficient algorithms, including proximal and accelerated
depend on a underlying low-rank matrix. The other method proximal gradient descent algorithms, and an algorithm which
[21] is developed for Poisson matrix completion but it does not is based on singular value thresholding that we examine in
details. This algorithm can be viewed as a consequence of parameters p and q, where p,q ∈ R is given by D(p(cid:107)q) (cid:44)
+
approximating the log likelihood function by its second order plog(p/q)−(p−q);theHellingerdistancebetweentwoPoisson
Taylor expansion and invoking a theorem for exact solution distributions with parameters p and q with p,q ∈R is given
of a nuclear norm regularized problem [12]. Our algorithm byd2 (p,q)(cid:44)2−2exp(cid:16)−1(cid:0)√p−√q(cid:1)2(cid:17),Wefur+therdefine
is related to [29]–[31] and can be viewed as a special case H 2
the average KL divergence and Hellinger distance for entries
where a simple closed form solution for the algorithm exists.
of two matrices P, Q∈Rd1×d2, where each entry corresponds
We further demonstrate the good performance of the algorithm +
to the parameter of a Poisson random variable:
in recovering solar flare images.
1 (cid:88)
Ourformulationandresultsareinspiredbytheseminalwork D(P(cid:107)Q)= d d D(Pij(cid:107)Qij),
of one-bit matrix completion [1], yet with several important 1 2 i,j
distinctions. In one-bit matrix completion, the value of each
1 (cid:88)
observation Yij is binary-valued and hence bounded, whereas d2H(P,Q)= d d d2H(Pij,Qij).
in our problem, each observation is a Poisson random variable 1 2 i,j
which is unbounded; hence, the arguments involve bounding
measurements have to be changed. In particular, we need to
II. FORMULATION
bound max Y when Y is a Poisson random variable with
ij ij ij
intensity Mij. Moreover, the Poisson likelihood function is non Suppose we observe a subset of entries of a matrix M ∈
Lipschitz (due to a bad point when Mij tends to zero), and Rd+1×d2 on the index set Ω ⊂ [d1] × [d2]. The indices are
hence we need to introduce a lower bound on each entry of the randomly selected with E|Ω|=m. In other words, I
{(i,j)∈Ω}
matrix Mij, which can be interpreted as the lowest required are i.i.d. Bernoulli random variables with parameter m/(d1d2).
SNR.Otherdistinctionsalsoincludeanalysistakingintoaccount The observations are Poisson counts of the observed matrix
of the property of the Poisson likelihood function, and using entries
Kullback-Leibler (KL) divergence as well as Hellinger distance Y ∼Poisson(M ), ∀(i,j)∈Ω. (1)
ij ij
that are different from those for the Bernoulli random variable
as used in [1]. Our goal is to recover the matrix M from the Poisson
observations {Y } .
ij (i,j)∈Ω
While working on this paper we realize a parallel work [32]
which also studies performance bounds for low rank matrix We make the following assumptions. First, we set an upper
completion with exponential family noise under more general bound α > 0 for the entries of M to entail the recovery
assumptions and using a different approach for proof (Poisson problem is well-posed [18]. This assumption is also reasonable
noiseisaspecialcaseoftheirs).TheirupperboundfortheMSE in practice; for instance, M may represent an image which
per entry is on the order of O(log(d +d )rmax{d ,d }/m) is usually not too spiky. Second, assume the rank of M is
1 2 1 2
(ourupperboundisO(cid:16)log(d1d2)(cid:112)r(d1+d2)/m(cid:17)),andtheir less than or equal to a positive integer r ≤min{d1,d2} (this
assumption is not restrictive in that we only assume an upper
lower bound is on the order of O(rmax{d ,d }/m) (versus
(cid:16)(cid:112) (cid:17) 1 2 bound on the rank). The third assumption is characteristic to
our lower bound is O r(d1+d2)/m ). Poisson matrix completion: we set a lower bound β > 0 for
eachentryM .Thisentry-wiselowerboundisrequiredforour
ij
The rest of the paper is organized as follows. Section II sets
later analysis, and it also has an interpretation of a minimum
up the formalism for Poisson matrix completion. Section III
required signal-to-noise ratio (S√NR), as the SNR of a Poisson
presents the matrix recovery based on constrained maximum
observation with intensity I is I.
likelihood and establishes the upper and lower bounds for the
recovery accuracy. Section IV presents an efficient iterative We recover the matrix M using a regularized maximum
algorithm that solves the maximum likelihood approximately likelihood formulation. Note that the log-likelihood function
and demonstrates its performance on recovering solar flare for the Poisson observation model (1) is proportional to
images. All proofs are delegated to Appendix.
(cid:88)
F (X)= Y logX −X , (2)
The notation in this paper is standard. In particular, R Ω,Y ij ij ij
+
(i,j)∈Ω
denotes the set of positive real numbers; [d] = {1,2,...,d};
I is the indicator function for an event ε; |A| denotes the where the subscript Ω and Y indicate the random quantities
[ε]
number of elements in a set A; diag{λ } denotes a diagonal involved in the maximum likelihood function F. Based on our
i
matrixwithasetofnumbers{λ }onitsdiagonal;1 denotes assumptions, we may define a set of candidate estimators
i n×m
an n-by-m matrix of all ones. Let entries of a matrix M be (cid:110) (cid:112)
denoted by Mij. Let (cid:107)M(cid:107) be the spectral norm which is the S (cid:44) X ∈R+d1×d2 :(cid:107)X(cid:107)∗ ≤α rd1d2, (3)
(cid:113)
largest absolute singular value, (cid:107)M(cid:107) = (cid:80) M2 be the β ≤X ≤α,∀(i,j)∈[d ]×[d ]}.
F i,j ij ij 1 2
Frobenius norm, (cid:107)M(cid:107)∗ be the nuclear norm which is the sum Here the upper bound on the nuclear norm (cid:107)M(cid:107) comes from
∗
of the singular values, and finally (cid:107)M(cid:107)∞ = maxij|Mij| be combining the assumptions (cid:107)M(cid:107) ≤ α and rank(M) ≤ r,
∞ √
the infinity norm. Let rank(M) denote the rank of a matrix M. (cid:112)
since (cid:107)M(cid:107) ≤ rank(M)(cid:107)M(cid:107) and (cid:107)M(cid:107) ≤ d d (cid:107)M(cid:107)
We say that a random variable X follows Poisson distribution ∗ √ F F 1 2 ∞
with parameter λ (or X ∼ Poisson(λ) if its probability mass lead to (cid:107)M(cid:107)∗ ≤ α rd1d2. An estimator M(cid:99) for M can be
function P(X = k) = e−λλk/(k!)). We also define the KL obtained by solving the following convex optimization problem:
divergence and Hellinger distance for Poisson distribution as M(cid:99)=argmaxFΩ,Y(X). (4)
follows: the KL divergence of two Poisson distributions with X∈S
2
III. PERFORMANCEBOUNDS any algorithm which, for any M ∈S, returns an estimator M(cid:99).
Then there exists M ∈ S such that with probability at least
In the following, we establish an upper bound and an 3/4,
information theoretic lower bound on the mean square error
(cid:40) (cid:114) (cid:41)
(MSE) per entry (cid:107)M(cid:99)−M(cid:107)2F/(d1d2) for the estimator in (4). d1d (cid:107)M −M(cid:99)(cid:107)2F ≥min C1,C2α3/2 rmaxm{d1,d2}
Theorem 1 (Upper bound). Assume M ∈ S, rank(M) = r, 1 2
Ω is chosen at random following our sampling model with (7)
E|Ω|=m,andM(cid:99)isthesolutionto(4).Thenwithaprobability as long as the right-hand side of (7) exceeds rα2/min{d1,d2},
where C ,C ,C are absolute constants.
exceeding (1−C/(d d )), we have 0 1 2
1 2
√
(cid:18) (cid:19)
1 (cid:107)M −M(cid:99)(cid:107)2 ≤C(cid:48) 8αT (α r)· Similarto[1],[33],proofofTheorem2reliesoninformation
d d F 1−e−T β theoretic arguments outlined as follows. First we find a set of
1 2
(cid:114) (cid:114) matricesχ∈S sothatthedistancebetweenanyX(i),X(j) ∈χ,
(cid:0)α(e2−2)+3log(d d )(cid:1) d1+d2 1+ (d1+d2)log(d1d2). identified as (cid:107)X(i)−X(j)(cid:107) , is sufficiently large. Then, for
1 2 m m F
(5) any X ∈S and the recovered X(cid:98), if we assume that they are
sufficiently close to each other with high probability, then we
If m≥(d1+d2)log(d1d2) then (5) simplifies to canclaimthatX istheelementinthesetS thatisclosesttoX(cid:98).
√
1 √ (cid:18) 8αT (cid:19)(cid:18)α r(cid:19) Finally, by applying a generalized Fano’s inequality involving
(cid:107)M −M(cid:99)(cid:107)2 ≤ 2C(cid:48) · KL divergence, we claim that the probability for the event that
d d F 1−e−T β
1 2 (6) X is the matrix in set S closest to X(cid:98) must be small, which
(cid:114)
(cid:0)α(e2−2)+3log(d d )(cid:1) d1+d2. leads to a contraction and hence proves our lower bound.
1 2 m
Remark3. TheassumptionsinTheorem2canbeachieved,for
Above, T,C(cid:48),C are absolute constants. where T,C,C(cid:48) are example, by the following construction. First, choose an α such
absolute constants. that α≥max{1,2β}, and then an r ≥4. Then, for d (or d )
1 2
sufficiently large, the conditions that α2rmax{d ,d } ≥ C
1 2 0
and the right-hand side of (7) exceeds rα2/min{d ,d } are
The proof of Theorem 1 is an extension of the ingenious 1 2
met. Since r ≤ O(min{d ,d }/α2), M ∈ S, what has been
arguments for one-bit matrix completion [1]. The extension 1 2
chosen is approximately low-rank. In other words, no matter
for Poisson case here is nontrivial for various aforementioned
how large r is, we can always find d (or d ) large enough so
reasons (notably the non sub-Gaussian and only locally sub- 1 2
that the assumptions in Theorem 2 are satisfied and thus there
Gaussian nature of the Poisson observations). An outline of
exist an M which can not be recovered with arbitrarily small
our proof is as follows. First, we establish an upper bound for
error by any method.
the KL divergence D(M(cid:107)X) for any M,X ∈S by applying
Lemma 2 given in the appendix. Second, we find an upper Remark 4. When m ≥ (d + d )log(d d ) and m =
1 2 1 2
bound for the Hellinger distance d2 (M,M(cid:99)) using the fact that O(r(d +d )logδ(d d )) with δ > 2, the ratio between the
H 1 2 1 2
theKLdivergencecanbeboundedfrombelowbytheHellinger upper bound in (6) and the lower bound in (7) is on the order
distance. Finally, we bound the mean squared error in Lemma of O(log(d d )). Hence, the lower bound matches the upper
1 2
3 via the Hellinger distance. bound up to a logarithmic factor.
Remark 1. Fixing d ,d ,m, α and β, the upper bound in
1 2
Theorem 1 increases as r increases. This is consistent with the IV. ALGORITHMS
intuitionthatourmethodisbetteratdealingwithapproximately
low-rank matrices (than with nearly full rank matrices). On the The matrix completion problem formulated in (4) is a
otherhand,fixingd ,d ,α,β andr,theupperbounddecreases Semidefinite program (SDP), since it is a nuclear norm
1 2
as m increases, which is also consistent with our intuition that minimization problem with a convex feasible domain. Hence,
M is supposed to be recovered more accurately with more wemaysolvedit,forexample,viatheinterior-pointmethod[34].
observations. Although the interior-point method returns an exact solution to
(4), it does not scale well with the dimensions of the matrix d
1
Remark 2. In the upper bound (5), the mean-square-error per
and d .
2
entrycanbearbitrarilysmall,inthesensethattheupperbound
goes to zero as d and d go to infinity when the number of Inthefollowing,wewilldevelopasetofiterativealgorithms
1 2
the measurements m = O((d +d )logδ(d d )) (m ≤ d d ) that solves the problem approximately and are more efficient
1 2 1 2 1 2
for δ >2 when r is fixed, or for δ >3 when r is sublinear on than solving the problem as SDP. In doing so, we use the
the order of o(log(d d )). framework of proximal algorithms to solve (4). At first, we
1 2
rewrite search space S as the intersection of two closed and
convex set in Rd1×d2:
The following theorem establishes an information theoretic
lower bound and demonstrates that there exists an M ∈S such Γ (cid:44){M ∈Rd1×d2 :(cid:107)M(cid:107) ≤r(cid:112)d d } and
1 ∗ 1 2
that any recovery method cannot achieve a mean square error
per entry less than the order of O(rmax{d1,d2}/m). Γ2 (cid:44){M ∈Rd1×d2 :β ≤Mij ≤α,∀(i,j)∈[d1]×[d2]},
Theorem 2 (Lower bound). Fix α, r, d , and d to be such where the first set is a nuclear norm ball and the second set
1 2
thatα,d ,d ≥1,r ≥4,α≥2β,andα2rmax{d ,d }≥C . is a high-dimensional box. Let f(M) (cid:44) −F (M) be the
1 2 1 2 0 Ω,Y
Let Ω be any subset of [d ]×[d ] with cardinality m. Consider negative log-likelihood function, then optimization problem (4)
1 2
3
is equivalent to
M(cid:99)=arg min f(M). (8)
M∈Γ1(cid:84)Γ2 The remaining of the problem is then to deal with the
projection onto the search space S. Since S is an intersection
(cid:84)
Noticing that the search space S = Γ1 Γ2 is closed and oftwoconvexsets,wemayusealternatingprojectionalgorithm
convex and f(M) is a convex function, we can use proximal to compute a sequence that converges to this intersection of
gradient methods to solve (8). Let IΓ(M) be an indicator Γ1 and Γ2, which is stated in Algorithm 3. Algorithm 3 is
function that takes value zero if M ∈Γ and is ∞ if M ∈Γc.
Then problem (8) is also equivalent to
Algorithm 3 Alternating Projection Algorithm
M(cid:99)=argM∈mRdi1n×d2f(M)+IΓ1(cid:84)Γ2(M). (9) 12:: Ifonritijali=ze1:,U20,.i.s.tdheomatrix needed to be projected onto S.
To guarantee the convergence of proximal gradient method, 3: Vj =ΠΓ1(Uj−1)
we need the Lipschitz constant L>0. In our case, Lipschitz 4: Uj =ΠΓ2(Vj)
constant L is a positive number satisfying 5: If (cid:107)Vj −Uj(cid:107)F ≤10−6 then return Uj.
6: end for
(cid:107)∇f(X)−∇f(Y)(cid:107) ≤L(cid:107)X−Y(cid:107) ,∀X,Y ∈S, (10)
F F
and hence L=α/β2 by the definition of our problem. Define efficient if some closed forms of projection onto the convex
sets can be achieved. Fortunately, computation of the projection
the projection of Y onto Γ as
onto Γ in our case is quite simple. Based on the definition
2
Π (Y)=argmin(cid:107)X−Y(cid:107)2. of Frobenius norm, [Π (Y)] is: β if Y <β; α if Y >α;
Γ X∈Γ F Y if β ≤Y ≤α. EvΓe2n if thijere is no cliojsed form expirjession
ij ij
for projection onto Γ , we can use TFOCS, a matlab package,
1
Algorithm 1 ProximalGradientforPoissonMatrixCompletion to implement this step.
1: Initialize: [M0]ij = Yij for (i,j) ∈ Ω and [M0]ij = (α+ Similar to the construction in [30], we may rewrite (4) as
β)/2 otherwise; the maximum number of iterations K.
2: for k =1,2,...K do M(cid:99)=arg min f(M)+λ(cid:107)M(cid:107)∗, (11)
3: Mk =ΠS(Mk−1−(1/L)∇f(Mk−1)) M∈Γ2
4: end for where λ is a regularizing parameter that balances the goodness
of data fit versus regularization.
Algorithm1haslinearconvergencerate,whichisestablished
The PMLSV algorithm can be derived as follows (in the
in the following theorem:
same spirit as [29], [19]). Let f(M) (cid:44) −F (M) be the
Ω,Y
Theorem3. Let{Mk}bethesequencegeneratedbyAlgorithm negative log-likelihood function. In the kth iteration, we may
1. Then for any k >1, we have form a Taylor expansion of f(M) around M , keep up to
k−1
second term and then solve
f(Mk)−f(M(cid:99))≤ L(cid:107)M02−kM(cid:99)(cid:107)2F. Mk =argMm∈iΓn2[Qtk(M,Mk−1)+λ(cid:107)M(cid:107)∗], (12)
with
Although Algorithm 1 can be implemented easily, its linear
Q (M,M )(cid:44)f(M )+(cid:104)M −M ,∇f(M )(cid:105)
convergence rate is not sufficiently if the Lipschitz constant L tk k−1 k−1 k−1 k−1
t
is large. In such scenarios, we prefer Nesterov’s accelerated + k(cid:107)M −M (cid:107)2, (13)
method for solving this problem which is our Algorithm 2. 2 k−1 F
Algorithm 2 has faster convergence, as stated in the following where ∇f is the gradient of f, t is the reciprocal of the
k
step size in the kth iteration, which we will specify later. By
Algorithm2AcceleratedProximalGradientforPoissonMatrix dropping and introducing terms independent of M whenever
Completion needed (more details can be found in [35]), (12) is equivalent
1: Initialize: [M0]ij = Yij for (i,j) ∈ Ω and [M0]ij = to
(α+β)/2 otherwise; Z = M ; the maximum number
0 0 M =
of iterations K. k
2: for k =1,2,...K do (cid:34)1(cid:13)(cid:13) (cid:18) 1 (cid:19)(cid:13)(cid:13)2 λ (cid:35)
3: Mk =ΠS(Zk−1−(1/L)∇f(Zk−1)) argmMin 2(cid:13)(cid:13)M − Mk−1− tk∇f(Mk−1) (cid:13)(cid:13)F + tk(cid:107)M(cid:107)∗ .
4: Zk =Mk+((k−1)/(k+2))(Mk−Mk−1) (14)
5: end for
Using a theorem proved in [12], we may show (in Appendix A)
Theorem 4. that the exact solution to (14) is given by a form of Singular
Value Thresholding (SVT):
Theorem4. Let{M }bethesequencegeneratedbyAlgorithm
k (cid:18) (cid:19)
2. Then for any k >1, we have 1
M =D M − ∇f(M ) , (15)
k λ/tk k−1 t k−1
f(Mk)−f(M(cid:99))≤ 2L(cid:107)(Mk0+−1)M(cid:99)2 (cid:107)2F. where Dτ(Σ)(cid:44)diag{(σi−τ)+} ankd (x)+ =max{x,0}.
4
Algorithm 4 PMLSV for Poisson Matrix Completion 50% and 30% of the image are observed. The results show that
our algorithm can recover the original image accurately when
1: Initialize: [M0]ij = Yij for (i,j) ∈ Ω and is (α+β)/2
50% or above of the image entries are observed. In the case of
otherwise,themaximumnumberofiterationsK,parameters
only 30% of the image entries are observed, our algorithm still
η, and L.
captures the main features in the image. The PMLSV algorithm
2: for k =1,2,...K do
is very efficient: the running time on a laptop with 2.40Hz two
3: C =Mk−1−(1/L)∇f(Mk−1)
4: C =UDVT {singular value decomposition} core CPU and 8GB RAM for all three examples are less than
1.2 seconds (much faster than solving SDP).
5: Dnew =diag((diag(D)−λ/L)+)
6: Mk =ΠΓ2(cid:0)UDnewVT(cid:1) Rank = 10
7: If f(Mk)>QL(Mk,Mk−1) then L=ηL, go to 4.
8: If |f(Mk)−QL(Mk,Mk−1)|<0.5/K then exit;
9: end for
Rank = 10
(a) p=0.8. (b) λ=0.1,K =2000.
Fig. 2: (a) Observed image with 80% of entries known (dark spots
represent missing entries). (b) Recovered image with λ = 0.1 and
Fig. 1: Solar flare image of size 48-by-48 with rank 10. no more than 2000 iterations, where the elapsed time is 1.176595
seconds.
The PMLSV algorithm is summarized in Algorithm 4. In
the algorithm description, L is the reciprocal of the step size,
η >1 is a scale parameter to change the step size, and K is
the maximum number of iterations, which is user specified:
a larger K leads to more accurate solution, and a small K
obtains the coarse solution quickly. If the cost function value
does not decrease, the step size is shortened to change the
singular values more conservatively. The algorithm terminates
whentheabsolutedifferenceinthecostfunctionvaluesbetween
two consecutive iterations is less than 0.5/K.
(a) p=0.5. (b) λ=0.1,K =2000.
Under the assumption that the box constraint is not binding,
Fig. 3: (a) Observed image with 50% of entries known (dark spots
Algorithm4isthesameasthatin[29]andconvergenceanalysis
represent missing entries). (b) Recovered image with λ = 0.1 and
can be found there.
no more than 2000 iterations, where the elapsed time is 1.110226
Despite of its simplicity, the PMLSV algorithm has a seconds.
surprisingly good performance. With the simple initialization
for M , the magnitude of the gradient is typically small at each
0
iteration. Hence, we can ensure M to be belong to or be close
k
to Γ by choosing an appropriate step size in the kth iteration.
The complexity of PMLSV is on the order of O(d2d +d3)
1 2 2
(which comes from the most expensive step of performing
singular value decomposition). This is much lower than the
complexity of solving an SDP, which is on the order of
O(d3+d d3+d2d2).Inparticular,forad-by-dmatrix,PMLSV
1 1 2 1 2
algorithm has complexity O(d2) versus solving the SDP has
(a) p=0.3. (b) λ=0.1,K =2000.
complexity O(d3).
Fig. 4: (a) Observed image with 30% of entries known (dark spots
V. NUMERICALEXAMPLE represent missing entries). (b) Recovered image with λ = 0.1 and
no more than 2000 iterations, where the elapsed time is 1.097281
We demonstrate the good performance of our estimator in seconds.
recovering a solar flare image. The solar flare image is of size
48-by-48. We break the image into 8-by-8 patches, then collect
the vectorized patches into a 64-by-36 matrix: such a matrix is
well approximated by a low-rank matrix, as demonstrated in
Fig. 1. ACKNOWLEDGEMENT
Suppose entries are observed using our sampling model TheauthorswouldliketothankProf.YuejieChi,Prof.Mark
with E|Ω|=m. Let p(cid:44)m/(d d ), then we observe (100p)% Davenport, and Prof. Yaniv Plan for stimulating discussions
1 2
of entries. We use L = 10−4 and η = 1.1 in the PMLSV and inspiring comments. This work is partially supported by
algorithm. Fig. 2 to Fig. 4 show the recovery result when 80%, NSF grant CCF-1442635.
5
REFERENCES poisson inverse problems with physical constraints,” arXiv preprint
arXiv:1403.6532,2014.
[1] M.A.Davenport,Y.Plan,E.v.d.Berg,andM.Wootters,“1-bitmatrix
[26] A.K.Q.Tran-DinhandV.Cevher,“Aproximalnewtonframeworkfor
completion,”InformationandInference,2014.
compositeminimization:Graphlearningwithoutcholeskydecomposition
[2] A.SIGKDD,“Netflix,”inProceedingsofkddcupandworkshop,2007. andmatrixinversions,”Proc.30thInt.Conf.MachineLearning(ICML),
[3] D. J. Brady, Optical imaging and spectroscopy. John Wiley & Sons, 2013.
2009. [27] A.K.Q.Tran-DinhandV.Cevher,“Compositeself-concordantmini-
[4] J.A.Bazerque,G.Mateos,andG.B.Giannakis,“Inferenceofpoisson mization,”J.MachineLearningResearch(LMLR),2014.
countprocessesusinglow-ranktensordata,”inAcoustics,Speechand
[28] M.J.Wainwright,“Structuredregularizersforhigh-dimensionalproblems:
SignalProcessing(ICASSP),IEEEInternationalConferenceon,pp.5989
Statisticalandcomputationalissues,”AnnualReviewofStatisticsandits
–5993,2013.
Applications,pp.233–253,2014.
[5] D. L. Donoho, “Compressed sensing,” Information Theory, IEEE
[29] S.JiandJ.Ye,“Anacceleratedgradientmethodfortracenormmini-
Transactionson,vol.52,no.4,pp.1289–1306,2006.
mization,”inProceedingsofthe26thAnnualInternationalConference
[6] E. J. Cande`s and B. Recht, “Exact matrix completion via convex onMachineLearning,pp.457–464,ACM,2009.
optimization,”FoundationsofComputationalmathematics,vol.9,no.6,
[30] M.J.Wainwright,“Structuredregularizersforhigh-dimensionalproblems:
pp.717–772,2009.
Statisticalandcomputationalissues,”AnnualReviewofStatisticsand
[7] R.H.Keshavan,A.Montanari,andS.Oh,“Matrixcompletionfroma ItsApplication,vol.1,pp.233–253,2014.
fewentries,”InformationTheory,IEEETransactionson,vol.56,no.6,
pp.2980–2998,2010. [31] A.Agarwal,S.Negahban,andM.J.Wainwright,“Fastglobalconver-
genceratesofgradientmethodsforhigh-dimensionalstatisticalrecovery,”
[8] E.J.Cande`sandT.Tao,“Thepowerofconvexrelaxation:Near-optimal inAdvancesinNeuralInformationProcessingSystems,pp.37–45,2010.
matrixcompletion,”InformationTheory,IEEETransactionson,vol.56,
no.5,pp.2053–2080,2010. [32] J.Lafond,“Lowrankmatrixcompletionwithexponentialfamilynoise,”
arXivpreprintarXiv:1502.06919,2015.
[9] W. Dai and O. Milenkovic, “Set: an algorithm for consistent matrix
completion,”inAcousticsSpeechandSignalProcessing(ICASSP),2010 [33] E.J.CandesandM.A.Davenport,“Howwellcanweestimateasparse
IEEEInternationalConferenceon,pp.3646–3649,IEEE,2010. vector?,”AppliedandComputationalHarmonicAnalysis,vol.34,no.2,
pp.317–323,2013.
[10] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank
solutions of linear matrix equations via nuclear norm minimization,” [34] Z.LiuandL.Vandenberghe,“Interior-pointmethodfornuclearnorm
SIAMreview,vol.52,no.3,pp.471–501,2010. approximationwithapplicationtosystemidentification,”SIAMJournal
[11] B.Recht,“Asimplerapproachtomatrixcompletion,”TheJournalof on Matrix Analysis and Applications, vol. 31, no. 3, pp. 1235–1256,
MachineLearningResearch,vol.12,pp.3413–3430,2011. 2009.
[12] J.-F. Cai, E. J. Cande`s, and Z. Shen, “A singular value thresholding [35] Y.CaoandY.Xie,“Low-rankmatrixrecoveryinpoissonnoise,”arXiv
algorithmformatrixcompletion,”SIAMJournalonOptimization,vol.20, preprintarXiv:1407.0726,2014.
no.4,pp.1956–1982,2010. [36] M.LedouxandM.Talagrand,ProbabilityinBanachSpaces:isoperimetry
[13] Z.Lin,A.Ganesh,J.Wright,L.Wu,M.Chen,andY.Ma,“Fastconvex andprocesses,vol.23. Springer,1991.
optimization algorithms for exact recovery of a corrupted low-rank [37] D. Pollard, A user’s guide to measure theoretic probability, vol. 8.
matrix,”ComputationalAdvancesinMulti-SensorAdaptiveProcessing CambridgeUniversityPress,2002.
(CAMSAP),vol.61,2009.
[38] B.Yu,“Assouad,fano,andlecam,”inFestschriftforLucienLeCam,
[14] R. Mazumder, T. Hastie, and R. Tibshirani, “Spectral regularization
pp.423–435,Springer,1997.
algorithms for learning large incomplete matrices,” The Journal of
MachineLearningResearch,vol.11,pp.2287–2322,2010.
[15] R.Keshavan,A.Montanari,andS.Oh,“Matrixcompletionfromnoisy
entries,”inAdvancesinNeuralInformationProcessingSystems,pp.952–
960,2009.
APPENDIXA
[16] E.J.CandesandY.Plan,“Matrixcompletionwithnoise,”Proceedings SINGULARVALUETHRESHOLDING
oftheIEEE,vol.98,no.6,pp.925–936,2010.
[17] S. Negahban, M. J. Wainwright, et al., “Estimation of (near) low-
Consider the following problem
rank matrices with noise and high-dimensional scaling,” The Annals
ofStatistics,vol.39,no.2,pp.1069–1097,2011. (cid:26)1 (cid:27)
[18] wS.eiNghegteadhbmaantraixndcoMm.pJl.etWioani:nOwpritgimhta,l“bRoeusntdrisctwedithstnroonisge,c”oTnhveexJiotyurannadl Y∈mRdi1n×d2 2(cid:107)Y −X(cid:107)2F +τ(cid:107)Y(cid:107)∗ , (16)
ofMachineLearningResearch,vol.13,no.1,pp.1665–1697,2012. where X ∈ Rd1×d2 is given and τ is the regularization
[19] A.Rohde,A.B.Tsybakov,etal.,“Estimationofhigh-dimensionallow- parameter.ForamatrixX ∈Rd1×d2 withrankr,letitssingular
r2a0n1k1.matrices,” The Annals of Statistics, vol. 39, no. 2, pp. 887–930, value decomposition be X = UΣVT, where U ∈ Rd1×r,
[20] A.Soni,S.Jain,J.Haupt,andS.Gonella,“Errorboundsformaximum V ∈Rd2×r,Σ=diag({σi},i=1,2,...,r),andσi isasingular
value of the matrix X. For each τ ≥ 0, define the singular
likelihoodmatrixcompletionundersparsefactormodels,”
value thresholding operator as:
[21] A.Soni,S.Jain,J.Haupt,andS.Gonella,“Noisymatrixcompletion
undersparsefactormodels,”arXivpreprintarXiv:1411.0282,2014. D (X)(cid:44)UD (Σ)VT. (17)
τ τ
[22] A.SoniandJ.Haupt,“Estimationerrorguaranteesforpoissondenoising
withsparseandstructureddictionarymodels,” The solution to (16) is given by singular value thresholding
[23] M. Raginsky, R. M. Willett, Z. T. Harmany, and R. F. Marcia, according to the following theorem
“Compressedsensingperformanceboundsunderpoissonnoise,”Signal
Processing,IEEETransactionson,vol.58,no.8,pp.3990–4002,2010. Theorem 5 (Theorem 2.1 in [12]). For each τ ≥ 0, and
[24] M.Raginsky,S.Jafarpour,Z.T.Harmany,R.F.Marcia,R.M.Willett,and X ∈Rd1×d2:
R. Calderbank, “Performance bounds for expander-based compressed (cid:26) (cid:27)
1
sensing in poisson noise,” Signal Processing, IEEE Transactions on, D (X)=arg min (cid:107)Y −X(cid:107)2 +τ(cid:107)Y(cid:107) . (18)
vol.59,no.9,pp.4139–4153,2011. τ Y∈Rd1×d2 2 F ∗
[25] X. Jiang, G. Raskutti, and R. Willett, “Minimax optimal rates for
6
 (cid:12)
(cid:12)
APPPREONODFIXS B =2hEsup (cid:12)(cid:12)(cid:12)(cid:88)(cid:15)ijI[(i,j)∈Ω](Yij(−logXij))
X∈S(cid:12)i,j
(cid:12)h
In the following, Lemma 1 is used in proving Lemma 2, (cid:12)
and Lemma 4 corresponds to Lemma 3 in [1]. +(cid:88)(cid:15)ijI[(i,j)∈Ω]Xij(cid:12)(cid:12)(cid:12) 
Lemma 1. Assuming Y ∼ Poisson (λ) is a Poisson random i,j (cid:12)
variable with λ ≤ α. Then P(Y −λ ≥ t) ≤ e−t,∀t ≥ t0 for  (cid:12) (cid:12)h
t (cid:44)α(e2−3). (cid:12) (cid:12)
0 ≤22h−1Esup (cid:12)(cid:12)(cid:12)(cid:88)(cid:15)ijI[(i,j)∈Ω](Yij(−logXij))(cid:12)(cid:12)(cid:12) 
X∈S(cid:12)i,j (cid:12)
Proof: (20)
 (cid:12) (cid:12)h
(cid:12) (cid:12)
We introduce θ ≥0, +22h−1Esup (cid:12)(cid:12)(cid:12)(cid:88)(cid:15)ijI[(i,j)∈Ω]Xij(cid:12)(cid:12)(cid:12) ,
X∈S(cid:12)i,j (cid:12)
P(Y −λ≥t)=P(Y ≥t+λ)
=P(θY ≥θ(t+λ))=P(exp(θY)≥exp(θ(t+λ))). where the expectation are over both Ω and Y.
Using Markov inequality, we can have
P(Y −λ≥t)≤exp(−θ(t+λ))E(cid:0)eθY(cid:1).
=exp(−θ(λ+t))·exp(cid:0)λ(eθ−1)(cid:1)
Letting θ =2,
P(Y −λ≥t) =exp(cid:0)−t+λ(e2−3)(cid:1).
exp(−t)
Inthefollowing,wewill usecontractionprincipletofurther
Define that bound the first term of (20). We let φ(t)=−βlog(t+1). We
t (cid:44)α(e2−3), know φ(0) = 0 and |φ(cid:48)(t)| = |β/(t+1)|, so |φ(cid:48)(t)| ≤ 1 if
0
t ≥ β −1. Setting Z = X −1 , then we have Z ≥
to make P(Y−λ≥t) ≤ 1, we derive that P(Y −λ≥t) ≤ e−t d1×d2 √ √ ij
exp(−t) β−1,∀(i,j)∈[d1]×[d2] and (cid:107)Z(cid:107)∗ ≤α rd1d2+ d1d2 by
when triangle inequality. Therefore, φ(Z ) is a contraction and it
ij
vanishes at 0. By Theorem 4.12 in [36] and using the fact that
t≥t ≥λ(e2−3)
0 |(cid:104)A,B(cid:105)|≤(cid:107)A(cid:107)(cid:107)B(cid:107) , we have
∗
 (cid:12) (cid:12)h
(cid:12) (cid:12)
(L2e)mamndaS2.beLettheFΩse,Yt d(Xefi)nebdeitnhe(3li)k,etlhihenood function defined in 22h−1EXsu∈pS(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)i,j (cid:15)ijI[(i,j)∈Ω](Yij(−logXij))(cid:12)(cid:12)(cid:12)(cid:12) 
(cid:26)
P sup |FΩ,Y(X)−EFΩ,Y(X)| (cid:20) (cid:21) (cid:12)(cid:12) (cid:12)(cid:12)h
≥X∈CS(cid:48)(cid:0)α√r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)· ≤22h−1E mi,ajxYihj Xsu∈pS(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)i,j (cid:15)ijI[(i,j)∈Ω]((−logXij))(cid:12)(cid:12)(cid:12)(cid:12) 
(cid:16)(cid:112) (cid:17)(cid:111) (19)
m(d1+d2)+d1d2log(d1d2)  (cid:12) (cid:12)h
(cid:20) (cid:21) (cid:12) (cid:18) (cid:19)(cid:12)
≤ d1Cd2, =22h−1E mi,ajxYihj EXsu∈pS(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)i,j (cid:15)ijI[(i,j)∈Ω] β1φ(Zij) (cid:12)(cid:12)(cid:12)(cid:12) 
where C(cid:48) and C are absolute positive constants and the  (cid:12) (cid:12)h
probability and the expectation are both over Ω and Y. ≤22h−1(cid:18)β2(cid:19)hE(cid:20)mi,ajxYihj(cid:21)EXsu∈pS(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:88)i,j (cid:15)ijI[(i,j)∈Ω]Zij)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 
Proof: In order to prove the lemma, we let (cid:15)ij are i.i.d. (cid:18)2(cid:19)h (cid:20) (cid:21) (cid:20) (cid:21)
Rademacher random variables. In the following derivation, the =22h−1 E maxYh E sup |(cid:104)∆ ◦E,Z(cid:105)|h
firstinequalityisduetheRadamachersymmetrizationargument β i,j ij X∈S Ω
(Lemma 6.3 in [36]) and the second inequality is due to the (cid:18)2(cid:19)h (cid:20) (cid:21) (cid:20) (cid:21)
panodwehr≥me1a.nTihneenquwaelithya:v(ea+b)h ≤2h−1(ah+bh) if a,b>0 ≤22h−1 β E mi,ajxYihj E Xsu∈pS(cid:107)E◦∆Ω(cid:107)h(cid:107)Z(cid:107)h∗
E(cid:20)sup |FΩ,Y(X)−EFΩ,Y(X)|h(cid:21) ≤22h−1(cid:18)β2(cid:19)h(cid:0)α√r+1(cid:1)h(cid:16)(cid:112)d1d2(cid:17)hE(cid:20)mi,ajxYihj(cid:21)E(cid:2)(cid:107)E◦∆Ω(cid:107)h(cid:3),
X∈S (21)
 (cid:12) (cid:12)h
≤2hEsup (cid:12)(cid:12)(cid:12)(cid:12)(cid:88)(cid:15)ijI[(i,j)∈Ω](YijlogXij −Xij)(cid:12)(cid:12)(cid:12)(cid:12)  wdehneorteesEthedeinndoitceastotrhemamtraitxrifxorwΩithanedn◦trideesnogtievsenthebyHa(cid:15)dija,m∆arΩd
X∈S(cid:12)i,j (cid:12) product.
7
Similarly,thesecondtermof(20)canbeboundedasfollows: Above, firstly we use triangle inequality and power mean
 (cid:12) (cid:12)h inequality, then along with independence, we use (24) in the
(cid:12) (cid:12) third inequality. By standard computations for exponential
22h−1Esup (cid:12)(cid:12)(cid:12)(cid:88)(cid:15)ijI[(i,j)∈Ω]Xij(cid:12)(cid:12)(cid:12)  random variables,
X∈S(cid:12)i,j (cid:12) (cid:20) (cid:21)
(cid:20) (cid:21) (22) E maxWihj ≤2h!+logh(d1d2). (26)
≤22h−1E sup (cid:107)E◦∆ (cid:107)h(cid:107)X(cid:107)h i,j
Ω ∗
X∈S Thus, we have
≤22h−1(cid:0)α√r(cid:1)h(cid:16)(cid:112)d d (cid:17)hE(cid:2)(cid:107)E◦∆ (cid:107)h(cid:3). (cid:20) (cid:21) (cid:16) (cid:17)
1 2 Ω E maxYh ≤22h−1 αh+(t )h+2h!+logh(d d ) .
ij 0 1 2
i,j
(27)
Plugging (21) and (22) into (20), we have
(cid:20) (cid:21)
E sup |F (X)−EF (X)|h
Ω,Y Ω,Y
X∈S
≤22h−1(cid:0)α√r+1(cid:1)h(cid:16)(cid:112)d1d2(cid:17)hE(cid:2)(cid:107)E◦∆Ω(cid:107)h(cid:3)· (23)
(cid:32)(cid:18)2(cid:19)h (cid:20) (cid:21) (cid:33)
E maxYh +1 .
β i,j ij
Therefore, combining (27) and (23), we have
To bound E(cid:2)(cid:107)E◦∆ (cid:107)h(cid:3), we can use the result from [1] (cid:20) (cid:21)
Ω E sup |F (X)−EF (X)|h
if we take h=log(d d )≥1: Ω,Y Ω,Y
1 2
X∈S
E(cid:2)(cid:107)E◦∆Ω(cid:107)h(cid:3) ≤24h−1(cid:0)α√r+1(cid:1)h(cid:16)(cid:112)d d (cid:17)hE(cid:2)(cid:107)E◦∆ (cid:107)h(cid:3)· (28)
1 2 Ω
(cid:115) h
≤C0(cid:16)2(1+√6)(cid:17)h m(d1+d2)+d1dd12d2log(d1d2) (cid:18)β2(cid:19)h(cid:16)αh+(t0)h+2h!+logh(d1d2)(cid:17).
Then,
for some constant C . Therefore, the only term we need to (cid:18) (cid:20) (cid:21)(cid:19)1
bound is E(cid:2)max Y0h(cid:3). E sup |F (X)−EF (X)|h h
i,j ij Ω,Y Ω,Y
X∈S
From Lemma 1, if t≥t0, then for any (i,j)∈[d1]×[d2], ≤16(cid:0)α√r+1(cid:1)(cid:16)(cid:112)d d (cid:17)E(cid:2)(cid:107)E◦∆ (cid:107)h(cid:3)h1 ·
the following inequality holds since t >α: 1 2 Ω
0
(cid:18) (cid:19)
2
P(|Yij −Mij|≥t)=P(Yij ≥Mij +t)+P(Yij ≤Mij −t) β (α+t0+2h+log(d1d2))
≤exp(−t)+0
=P(W ≥t), ≤16(cid:18)2(cid:19)(cid:0)α√r+1(cid:1)(cid:16)(cid:112)d d (cid:17)E(cid:2)(cid:107)E◦∆ (cid:107)h(cid:3)h1 ·
ij β 1 2 Ω
(24)
(cid:0)α(e2−2)+3log(d d )(cid:1)
wvahreiarbeleWs.ij are independent standard exponential random ≤128(cid:16)1+√6(cid:17)Ch1 (cid:18)α1√2r(cid:19)(cid:0)α(e2−2)+3log(d d )(cid:1)·
0 β 1 2
Below we use the fact that for any positive random variable (cid:16)(cid:112) (cid:17)
X, we can write EX =(cid:82)∞P(X ≥t)dt, allowing us to bound m(d1+d2)+d1d2log(d1d2) .
0
(cid:20) (cid:21) (29)
E maxYh
i,j ij whereweusethefactthat(ah+bh+ch+dh)1/h ≤a+b+c+dif
(cid:18) (cid:20) (cid:21)(cid:19) a,b,c,d>0inthefirstinequalityandwetakeh=log(d d )≥
1 2
≤22h−1 αh+E max|Yij −Mij|h 1 in the second inequality.
i,j
(cid:18) (cid:90) ∞ (cid:18) (cid:19) (cid:19)
=22h−1 αh+ P max|Y −M |h ≥t dt
ij ij
0 i,j
(cid:32) (cid:90) ∞ (cid:18) (cid:19) (cid:33)
≤22h−1 αh+(t )h+ P max|Y −M |h ≥t dt
0 ij ij
(t0)h i,j
(cid:32) (cid:90) ∞ (cid:18) (cid:19) (cid:33)
≤22h−1 αh+(t )h+ P maxWh ≥t dt √
0 (t0)h i,j ij Moreover when C(cid:48) ≥128(cid:0)1+ 6(cid:1)e,
≤22h−1(cid:18)αh+(t0)h+E(cid:20)mi,ajxWihj(cid:21)(cid:19) C (cid:32)128(1+√6)(cid:33)log(d1d2) ≤ C0
(25) 0 C(cid:48) d1d2
8
Therefore we can use Markov inequality to see that 1. For all X ∈χ, each entry has |X |=αγ.
ij
(cid:26)
P sup |F (X)−EF (X)|
Ω,Y Ω,Y
≥X∈CS(cid:48)(cid:0)α√r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)· α2γ22.d1Fdo2r/2a.ll X(i),X(j) ∈ χ, i (cid:54)= j, (cid:107)X(i) − X(j)(cid:107)2F >
(cid:16)(cid:112)m(d +d )+d d log(d d )(cid:17)(cid:111) Lemma 5. For x,y >0, D(x(cid:107)y)≤(y−x)2/y.
1 2 1 2 1 2
(cid:26)
=P sup |F (X)−EF (X)|h
Ω,Y Ω,Y
X∈S √ Proof: First assume x ≤ y. Let z = y−x. Then z ≥ 0
≥C(cid:48)(cid:0)α r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)· and D(x(cid:107)x+z) = xlogx+xz +z. Taking the first derivative
(cid:16)(cid:112) (cid:17)(cid:111) of this with respect to z, we have ∂ D(x(cid:107)x + z) = z .
m(d1+d2)+d1d2log(d1d2) ∂z x+z
Thus, by Taylor’s theorem, there is some ξ ∈ [0,z] so that
(cid:20) (cid:21)
≤E sup |FΩ,Y(X)−EFΩ,Y(X)|h / D(x(cid:107)y)=D(x(cid:107)x)+z·x+ξξ.Sincetheright-hand-sideincreases
X∈S√ in ξ, we may replace ξ with z and obtain D(x(cid:107)y)≤ (y−x)2.
{(cid:0)C(cid:48)(cid:0)α r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)· For x > y, with the similar argument we may concludeythat
(cid:16)(cid:112) (cid:17)(cid:17)h for z =y−x<0 there is some ξ ∈[z,0] so that D(x(cid:107)y)=
m(d1+d2)+d1d2log(d1d2) } D(x(cid:107)x)+z· ξ . Since z <0 and ξ/(x+ξ) increases in ξ,
x+ξ
C thentheright-hand-sideisdecreasinginξ.Wemayalsoreplace
≤ ,
d d ξ with z and this proves the lemma.
1 2
√
where C(cid:48) ≥128(1+ 6)e and C are absolute constants.
Proof of Theorem 1: Lemma 1, Lemma 2, and Lemma
3 are used in the proof. In the following, the expectation are
taken with respect to both Ω and {Y }. First, note that
Lemma 3. Let β ≤Mij,M(cid:99)ij ≤α,∀(i,j)∈[d1]×[d2], then ij
(cid:20) (cid:18) (cid:19) (cid:21)
d2H(M,M(cid:99))≥ 1−4αeT−T (cid:107)Md−dM(cid:99)(cid:107)2F, FΩ,Y(X)−FΩ,Y(M)=(i(cid:88),j)∈Ω Yijlog MXiijj −(Xij −Mij) .
1 2
where T = 1 (α−β)2. Then for any X ∈S,
8β
E[F (X)−F (M)]
Ω,Y Ω,Y
(cid:20) (cid:18) (cid:19) (cid:21)
Proof: Assuming x is any entry in M and y is any entry = m (cid:88) M log Xij −(X −M )
in M(cid:99), then β ≤ x,y ≤ α and 0 ≤ |x−y| ≤ α−β. By the d d ij M ij ij
1 2 ij
mean value theorem there exists an ξ(x,y)∈[β,α] such that i,j
(cid:20) (cid:18) (cid:19) (cid:21)
1 √ √ 1(cid:32) 1 (cid:33)2 1 =−dmd (cid:88) Mijlog MXij −(Mij −Xij) (32)
( x− y)2 = (x−y) = (x−y)2 ≤T. 1 2 i,j ij
(cid:112)
2 2 2 ξ(x,y) 8ξ(x,y) m (cid:88)
=− D(M (cid:107)X )=−mD(M(cid:107)X).
d d ij ij
The function f(z) = 1−e−z is concave in [0,+∞], so if 1 2 i,j
z ∈[0,T], we may bound it from below with a linear function
For M ∈ S, we know M(cid:99) ∈ S and FΩ,Y(M(cid:99)) ≥ FΩ,Y(M).
1−e−T Thus we write
1−e−z ≥ z. (30)
√ √ T 0≤FΩ,Y(M(cid:99))−FΩ,Y(M)
Pluggingz = 12( x− y)2 = 8ξ(1x,y)(x−y)2 in(30),wehave =FΩ,Y(M(cid:99))+EFΩ,Y(M(cid:99))−EFΩ,Y(M(cid:99))
(cid:18) 1 √ √ (cid:19) 1−e−T 1 +EFΩ,Y(M)−EFΩ,Y(M)−FΩ,Y(M)
2−2exp − ( x− y)2 ≥ (x−y)2 (cid:104) (cid:105)
2 T 4ξ(x,y) ≤E FΩ,Y(M(cid:99))−FΩ,Y(M) +
1−e−T 1 (cid:12) (cid:12)
≥ T 4α(x−y)2. (cid:12)(cid:12)FΩ,Y(M(cid:99))−EFΩ,Y(M(cid:99))(cid:12)(cid:12)+|FΩ,Y(M)−EFΩ,Y(M)|
(31)
≤−mD(M(cid:107)M(cid:99))+2 sup |FΩ,Y(X)−EFΩ,Y(X)|.
Note that (31) holds for any x and y. This concludes the proof. X∈S
Applying Lemma 2, we obtain that with probability at least
√
Lemma 4. Let H (cid:44) (cid:8)M :(cid:107)M(cid:107) ≤α rd d ,(cid:107)M(cid:107) ≤α(cid:9) (1−C/(d1d2)),
∗ 1 2 ∞
and γ ≤1 be such that r is an integer. Suppose r/γ2 ≤d ,
γ2 1 0≤−mD(M(cid:107)M(cid:99))
then we may construct a set χ∈H of size +2C(cid:48)(cid:0)α√r/β(cid:1)(cid:0)α(e2−2)+3log(d d )(cid:1)·
(cid:18) (cid:19) 1 2
|χ|≥exp rd2 (cid:16)(cid:112)m(d +d )+d d log(d d )(cid:17).
16γ2 1 2 1 2 1 2
√
with the following properties: After rearranging terms and applying the fact that d d ≤
1 2
9
d1+d2, we obtain Now consider an algorithm that for any X ∈S returns X(cid:98)
√ such that
D(M(cid:107)M(cid:99))≤2(cid:16)C(cid:112)(cid:48)(cid:0)α r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)(cid:17)· (33) d1d (cid:107)X−X(cid:98)(cid:107)2F <(cid:15)2 (36)
m(d +d )+(d +d )2log(d d ) . 1 2
1 2 1 2 1 2
with probability at least 1/4. Next, we will show this leas to
Note that the KL divergence can be bounded below by the an contradiction. Let
Hellinger distance (Chapter 3 in [37]):
X∗ =arg min (cid:107)X(i)−X(cid:98)(cid:107)2,
d2 (x,y)≤D(x(cid:107)y). X(i)∈χ F
H
by the same argument as that in [1], we have X∗ = X as
Thus from (33), we obtain
long as (36) holds. Using the assumption that (36) holds with
√
d2H(M,M(cid:99))≤2C(cid:48)(cid:0)α r/β(cid:1)(cid:0)α(e2−2)+3log(d1d2)(cid:1)· probability at least 1/4, we have
(cid:16)(cid:112)m(d1+d2)+(d1+d2)2log(d1d2)(cid:17). P(X∗ (cid:54)=X)≤ 3. (37)
(34) 4
Using a generalized Fano’s inequality for the KL divergence in
Finally, Theorem 1 is proved by applying Lemma 3.
[38], we have
max D(X(k)(cid:107)X(l))+1
P(X∗ (cid:54)=X)≥1− X(k)(cid:54)=X(l) . (38)
Proof of Theorem 2:
log|χ|
We will prove by contradiction. Lemma 4 and Lemma 5 are Define D (cid:44) D(X(k)(cid:107)X(l)) = (cid:80) D(X(k)(cid:107)X(l)). We
used in the proof. Without loss of generality, assume d ≥d . (i,j)∈Ω ij ij
2 1 know that each term in the sum is either 0, D(α(cid:107)α(cid:48)), or
Choose (cid:15)>0 such that
D(α(cid:48)(cid:107)α). From Lemma 5, since α(cid:48) <α, we have
(cid:40) (cid:114) (cid:41)
1 rd
(cid:15)2 =min ,C α3/2 2 , m(γα)2 64m(cid:15)2
256 2 m D ≤ ≤ .
α(cid:48) α(cid:48)
where C2 is an absolute constant that will be be specified later. Combining (37) and (38), we have that
First, choose γ such that r is an integer and
γ2 1 D+1
√ ≤1−P(X (cid:54)=X∗)≤
4 2(cid:15) 8(cid:15) 1 4 log|χ|
α ≤γ ≤ α ≤ 2. (cid:32)64m(cid:15)2 +1(cid:33) (cid:32)64m(cid:15)2 +1(cid:33) (39)
≤16γ2 α(cid:48) ≤1024(cid:15)2 α(cid:48) .
We may make such a choice because rd α2rd
2 2
α2r r α2r
≤ ≤ Suppose 64m(cid:15)2 ≤α(cid:48), then with (39), we have
64(cid:15)2 γ2 32(cid:15)2
1 2
and ≤1024(cid:15)2 ,
α2r α2r α2r 4 α2rd2
− = >4α2r >1.
32(cid:15)2 64(cid:15)2 64(cid:15)2 which implies that α2rd ≤32. Then if we set C >32, this
2 0
leads to a contradiction. Next, suppose 64m(cid:15)2 >α(cid:48), then with
Furthermore, since we have assumed that (cid:15)2 is larger than
(39), we have
Crα2/d , r/γ2 ≤ d for an appropriate choice of C. Let
1 1
χ(cid:48) bethesetdefinedinLemma4,byreplacingαwithα/2 1 (cid:18) 128m(cid:15)2 (cid:19)
α/2,γ <1024(cid:15)2 .
and with this choice of γ. Then we can construct a packing set 4 (1−γ)α3rd
2
χ of the same size as χ(cid:48) by defining
α/2,γ Since 1−γ >1/2, we have
(cid:110) (cid:16) γ(cid:17) (cid:111)
χ(cid:44) X(cid:48)+α 1− 1 :X(cid:48) ∈χ(cid:48) . (cid:114)
2 d1×d2 α/2,γ (cid:15)2 > α3/2 rd2.
The distance between pairs of elements in χ is bounded since 1024 m
α2γ2d d Setting C2 ≤1/4096, this leads to a contradiction. Therefore,
(cid:107)X(i)−X(j)(cid:107)2 ≥ 1 2 ≥4d d (cid:15)2. (35) (36) must be incorrect with probability at least 3/4. This
F 4 2 1 2
concludes our proof.
Define α(cid:48) (cid:44)(1−γ)α, then every entry of X ∈χ has X ∈
ij
{α,α(cid:48)}. Since we have assumed r ≥4, for every X ∈χ, we
have
(cid:16) γ(cid:17) γ (cid:112)
(cid:107)X(cid:107) =(cid:107)X(cid:48)+α 1− 1 (cid:107) ≤(cid:107)X(cid:48)(cid:107) +α(1− ) d d
∗ 2 d1×d2 ∗ ∗ 2 1 2Lemma 6. If f is a closed convex function satisfying Lipschitz
α(cid:112) (cid:112) (cid:112) condition (10), then for any X,Y ∈S, the following inequality
≤ rd d +α d d ≤α rd d ,
2 1 2 1 2 1 2 holds:
for some X(cid:48) ∈χ(cid:48) . Since the γ we choose is less than 1/2, L
α/2,γ f(Y)≤f(X)+(cid:104)∇f(X),Y −X(cid:105)+ (cid:107)Y −X(cid:107)2.
α(cid:48) is greater than α/2. Therefore, from the assumption that 2 F
β ≤α/2, we conclude that χ⊂S.
10