ESTIMATION OF LOW RANK DENSITY MATRICES: BOUNDS IN
SCHATTEN NORMS AND OTHER DISTANCES

arXiv:1604.04600v1 [stat.ML] 15 Apr 2016

By Dong Xia

∗

and Vladimir Koltchinskii

†

School of Mathematics
Georgia Institute of Technology
Let Sm be the set of all m × m density matrices (Hermitian positively semi-definite matrices of unit trace). Consider a problem of
estimation of an unknown density matrix ρ ∈ Sm based on outcomes
of n measurements of observables X1 , . . . , Xn ∈ Hm (Hm being the
space of m × m Hermitian matrices) for a quantum system identically prepared n times in state ρ. Outcomes Y1 , . . . , Yn of such measurements could be described by a trace regression model in which
Eρ (Yj |Xj ) = tr(ρXj ), j = 1, . . . , n. The design variables X1 , . . . , Xn
are often sampled at random from the uniform distribution in an
orthonormal basis {E1 , . . . , Em2 } of Hm (such as Pauli basis). The
goal is to estimate the unknown density matrix ρ based on the data
(X1 , Y1 ), . . . , (Xn , Yn ). Let
Ẑ :=

n
m2 X
Yj Xj
n
j=1

and let ρ̌ be the projection of Ẑ onto the convex set Sm of density matrices. It is shown that for estimator ρ̌ the minimax lower
bounds in classes of low rank density matrices (established earlier)
are attained up logarithmic factors for all Schatten p-norm distances,
p ∈ [1, ∞] and for Bures version of quantum Hellinger distance. Moreover, for a slightly modified version of estimator ρ̌ the same property
holds also for quantum relative entropy (Kullback-Leibler) distance
between density matrices.

1. Introduction. Let Mm be the set of all m × m matrices with entries in C.
For A ∈ Mm , let A∗ denote its conjugate transpose and let tr(A) denote the trace of
A. The complex linear space Mm of dimension m2 will be equipped with the HilbertSchmidt inner product hA, Bi = tr(AB ∗ ), A, B ∈ Mm . In what follows, the sign ⊗
denotes the tensor product of vectors or matrices (linear transformations). For instance,
for u, v ∈ Cm , u ⊗ v is a linear transformation from Cm into itself defined as follows:
(u ⊗ v)x = uhx, vi, x ∈ Cm .
Let
Hm := {A ∈ Mm : A = A∗ }
∗
†

Supported in part by NSF Grants DMS-1207808 and CCF-1523768
Supported in part by NSF Grants DMS-1509739 and CCF-1523768

1

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

2

be the set of all Hermitian matrices. Clearly, Hm is a linear space of dimension m2 over
the field of real numbers. For A ∈ Hm , the notation A < 0 means that A is positively
semi-definite. A density matrix is a positively semi-definite Hermitian matrix of unit
trace. The set of all m × m density matrices will be denoted by
Sm := {S ∈ Hm : S < 0, tr(S) = 1}.
Density matrices are used in quantum mechanics to characterize the states of quantum
systems. More generally, the states are represented by self-adjoint positively semidefinite
operators of unit trace acting in an infinite-dimensional Hilbert space. In this case,
density matrices of a large dimension m could be used to approximate the states of the
system.
The goal of quantum state tomography is to estimate the density matrix for a system
prepared in an unknown state based on specially designed measurements. Let X ∈ Hm be
P ′
a Hermitian matrix (an observable) with spectral representation X = m
j=1 λj Pj , where
m′ ≤ m, λj ∈ R, j = 1, . . . , m′ being the distinct eigenvalues of X and Pj , j = 1, . . . , m′
being the corresponding eigenprojections. For a system prepared in state ρ ∈ Sm , possible
outcomes of a measurement of observable X are the eigenvalues λj , j = 1, . . . , m′ and
they occur with probabilities pj := tr(ρPj ), j = 1, . . . , m′ . If Y is a random variable
representing such an outcome, then
Eρ Y = tr(ρX) = hρ, Xi.
In a simple model of quantum state tomography considered in this paper, an observable
X is sampled at random from some probability distribution Π in Hm , Eρ (Y |X) = hρ, Xi
and Y = hρ, Xi + ξ with noise ξ such that Eρ (ξ|X) = 0. Given a sample X1 , . . . , Xn of n
i.i.d. copies of X, n measurements of observables X1 , . . . , Xn are performed for a system
identically prepared n times in the same unknown state ρ ∈ Sm resulting in outcomes
Y1 , . . . , Yn . This leads to the following trace regression model
(1.1)

Yj = hρ, Xj i + ξj , j = 1, . . . , n

with design variables Xj , j = 1, . . . , n, response variables Yj , j = 1, . . . , n and noise ξj , j =
1, . . . , n satisfying the assumption Eρ (ξj |Xj ) = 0, j = 1, . . . , n and Eρ (Yj |Xj ) = hρ, Xj i.
The goal is to estimate the target density matrix ρ based on the data (X1 , Y1 ), . . . , (Xn , Yn ),
with the estimation error being measured by one of the statistically meaningful distances
between density matrices such as the Schatten p-norm distances for p ∈ [1, ∞] or quantum versions of Hellinger and Kullback-Leibler distances.
This version of the problem of quantum state tomography has been intensively studied in the recent years. The noiseless case (quantum compressed sensing) was considered
in [12] and [11]. In these papers, sharp bounds on the number n of measurements needed
to recover a density matrix of rank r were obtained based on a subtle argument (so called
“golfing scheme”) utilizing matrix Bernstein type inequalities. These developments were
related to an earlier work on low rank matrix completion [7]. In the noisy case, trace

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

3

regression problems have been studied by many authors (see, e.g., [15] and references
therein). The main focus was on nuclear norm penalized least squares estimator (matrix
LASSO) and related methods such as matrix Dantzig selector (see [6], [19], [25], [14]). In
[21], sharp bounds for matrix LASSO and matrix Dantzig selector, in particular, for Pauli
measurements in quantum state tomography were obtained. Most of the results in these
papers included upper bounds on the estimation error in Hilbert–Schmidt (Frobenius)
norm as well as low rank oracle inequalities ([19], [15], [18]). In [19], an upper bound on
the operator norm error of a nuclear norm penalized modified least squares estimator
was also proved. This result was further developed in [22]. In [16], upper bounds and
low rank oracle inequalities for von Neumann entropy penalized least squares estimators
were studied (including the bounds on the error in Bures distance and quantum relative
entropy distance). A rank penalized estimator of density matrix was studied in [1]. The
minimax lower bounds on the Frobenius norm error for matrix completion problems in
classes of matrices of rank r were obtained in [19] (the operator norm version could be
found in [22]). In [23], a method of deriving lower bounds for unitary invariant matrix
norms (including Schatten p-norms) was developed and, among other matrix estimation
problems, such bounds were obtained for matrix completion. Minimax lower bounds on
the nuclear norm error in density matrix estimation were obtained in [9], where it was
also shown that these bounds are attained (up to logarithmic factors) for the matrix
versions of LASSO and Dantzig selector. In our recent paper [20], we derived minimax
lower bounds in classes of low rank density matrices for the whole range of Schatten
p-norm distances as well as for Bures (quantum Hellinger) and quantum relative entropy
distance. We also showed that these minimax bounds are attained (up to logarithmic
factors) for von Neumann entropy penalized least squares estimators introduced in [16]
simultaneously for Bures, relative entropy and Schatten p-norm distances for p ∈ [1, 2].
The current paper could be viewed as a continuation of [20]. Our main goal is to
study a minimal distance estimator ρ̌ of ρ (initially proposed in [17]) defined as the
projection of a simple unbiased estimator
Ẑ =

n
m2 X
Yj Xj
n j=1

onto the convex set of density matrices Sm . We show that the minimax error rates
established in [20] for the classes of low rank density matrices are attained for this
estimator up to logarithmic factors in the whole range of Schatten p-norm distances for
p ∈ [1, ∞] as well as for Bures and relative entropy distance. The proof of these results
relies on simple properties of projections of Hermitian matrices onto the convex set Sm
of density matrices (see theorems 7 and 8) that might be of independent interest.
Throughout the paper, h·, ·i denotes either Hilbert–Schmidt inner product (defined
above), or (with a little abuse of notation) the canonical inner product of Cm . The
corresponding norm in Cm is denoted by | · |. For A, B ≥ 0, the notation A . B means
that A ≤ CB for a numerical constant C > 0, A & B means that B . A and A ≍ B
means that B . A . B. If needed, these signs might be provided with subscripts

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

4

indicating that the constant is allowed to depend on parameters. Say, A .γ B would
mean that A ≤ CB with C depending on γ.
2. Preliminaries.
2.1. Distances between density matrices. The Schatten p-norm of a matrix A ∈ Hm
is defined as


kAkp :=

m
X

j=1

1/p

|λj (A)|p

, p ∈ [1, +∞],

where λ1 (A) ≥ · · · ≥ λm (A) are the eigenvalues of A arranged in a non-increasing order.
For p = 1, the norm kAk1 is called the nuclear or the trace norm; for p = 2, kAk2 is the
Hilbert–Schmidt (generated by the Hilbert–Schmidt inner product) or Frobenius norm;
for p = +∞, kAk∞ = max1≤j≤m |λj (A)| is called the operator or the spectral norm.
Note that, for all A ∈ Hm , [1, ∞] ∋ p 7→ kAkp is a non-increasing function. The following
interpolation inequality is well known and can be easily deduced from a similar result for
1
ℓp -norms. Let 1 ≤ p < q < r ≤ ∞ and let µ ∈ [0, 1] be such that µp + 1−µ
r = q , then
kAkq ≤ kAkµp kAk1−µ
, A ∈ Hm .
r

(2.1)

In addition to the distances generated by the Schatten p-norms, the following two distances (extending well known distances between probability distributions used in the
classical statistics) are of importance in quantum statistics: Bures distance and KullbackLeibler divergence. The Bures distance is a quantum version of Hellinger distance and it
is defined as follows:
q

1/2

1/2

H 2 (S1 , S2 ) := 2 − 2tr S1 S2 S1 , S1 , S2 ∈ Sm .
q

1/2

1/2

The quantity tr S1 S2 S1 is called the fidelity of states S1 , S2 (a quantum version of
Hellinger affinity). Note that 0 ≤ H 2 (S1 , S2 ) ≤ 2 and that H(S1 , S2 ) defines a metric
in the space Sm . The non-commutative Kullback-Leibler divergence, or relative entropy
distance is defined as
K(S1 kS2 ) := tr(S1 log S1 − S1 log S2 ), S1 , S2 ∈ Sm .
If S2 is a density matrix of rank strictly smaller than m, log S2 is not well defined and
K(S1 kS2 ) := +∞. Clearly, K(S1 kS2 ) is not a metric (in particular, it is not symmetric).
It is well known that K(S1 kS2 ) is the supremum of classical Kullback-Leibler divergences
between the distributions of outcomes of all possible measurements (represented by positive operator valued measures (POVM)) for the system prepared in states S1 and S2 .
Similar property holds also for the Bures (Hellinger) distance and for the nuclear norm
distance kS1 − S2 k1 which is the supremum of classical total variation distances between

5

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

the distributions of outcomes of all measurements (see [26], [13]). These observations
easily imply the following inequalities:
1
kS1 − S2 k21 ≤ H 2 (S1 , S2 ) ≤ K(S1 kS2 ) ∧ kS1 − S2 k1
4

(2.2)
(see also [16]).

2.2. Sampling from an orthonormal basis. Uniform sampling from an orthonormal basis is a model of design distribution in trace regression (1.1) that has been
frequently used in the literature on quantum compressed sensing (see, [12], [11]). Let
E := {E1 , . . . , Em2 } be an orthonormal basis of the space Hm of Hermitian matrices. Let
U := max kEj k∞ .
1≤j≤m2

Clearly, U ≤ 1 and

1 = max 2 kEj k2 ≤ m1/2 U,
1≤j≤m

implying that U ≥ m−1/2 . In what follows, it will be assumed that Π is a uniform
distribution on the basis E. As a result, the response variables Yj , j = 1, . . . , n of trace
regression model (1.1) could be viewed as noisy measurements of n randomly picked
Fourier coefficients of the target density matrix ρ in basis E. This model includes, in
particular, the so called Pauli measurements, an important approach to quantum state
tomography (see, e.g., [12], [11]).
Example: Pauli bases and Pauli measurements. The space of observables for
a single qubit system is the space H2 of 2 × 2 Hermitian matrices. Let
σ0 :=

1 0
0 1

!

,

σ1 :=

0 1
1 0

!

,

σ2 :=

0 i
−i 0

!

,

σ3 :=

1 0
0 −1

!

.

The matrices σ1 , σ2 , σ3 (often denoted σx , σy , σz ) are called Pauli matrices. The matrices
Wi = √12 σi , i = 0, 1, 2, 3 form an orthonormal basis of the space H2 (the Pauli basis). For

a system consisting of k qubits, the space of observables is Hm , where m = 2k . The Pauli
basis of this space is defined by tensorizing the Pauli basis of H2 : it consists of m2 = 4k
tensor products Wi1 ⊗. . .⊗Wik , (i1 , . . . , ik ) ∈ {0, 1, 2, 3} k . Let E1 = W0 ⊗. . .⊗W0 and let
E2 , . . . , Em2 be the rest of the matrices of the Pauli basis of Hm . It is straightforward to
check that E1 = √1m Im , where Im denotes m × m identity matrix (thus, √1m is the only
eigenvalue of E1 ). Matrices E2 , . . . , Em2 have eigenvalues ± √1m . Therefore, kEj k∞ =

m−1/2 , implying that, for the Pauli basis, U = m−1/2 . The fact that the matrices of this
basis have the smallest possible operator norms has been used in quantum compressed
sensing (see [12], [11], [21]). Matrices Ej have the following spectral representations:
Ej = √1m Pj+ − √1m Pj− with eigenprojections Pj+ , Pj− , j = 1, . . . , m2 (for E1 , P1− = 0). A
measurement of Ej for a k qubit system prepared in state ρ results in a random outcome

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

6

τj with two possible values ± √1m taken with probabilities ρ, Pj± . For random variable
τj , Eρ τj = hρ, Ej i. The density matrix ρ admits the following representation in the Pauli
basis:
m2
X
α
√ j Ej
ρ=
m
j=1
with α1 = 1 and with some αj ∈ R, j = 2, . . . , m2 . This implies that Eρ τj =

α
√j ,
m

n
1 o 1 ± αj
=
P ρ τj = ± √
m
2
1−α2j
m . Note
2, . . . , m2 , |αj | < 1

and Varρ (τj ) =

n

that, for j = 1, α1 = 1, Pρ τ1 =

√1
m

o

= 1 and Varρ (τ1 ) = 0.

For j =
and Varρ (τj ) > 0.
Let ν be picked at random from the set {1, . . . , m2 } (with the uniform distribution)
and let X = Eν , Y = τν (which corresponds to random sampling from the Pauli basis
with a subsequent measurement of observable X resulting in the outcome Y ). Then
2
ν
Eρ (Y |X) = hρ, Xi and Var ρ (Y |X) = 1−α
m . Moreover, we have
2

m
n
α2j
2kρk22
1o
2 X
1 o
=
.
= P α2ν ≥
≤ 2Eα2ν =
P Varρ (Y |X) ≤
2m
2
m j=1 m
m

n

2
Since, for ρ ∈ Sm , kρk2 ≤ 1, this means that, for m > 2 with probability at least 1 − m
,
1
1
2
Var ρ (Y |X) > 2m . In other words, the number of j = 1, . . . , m such that Var ρ (τj ) > 2m
1
.
is at least m2 − 2m implying that, for the most of the values of j, Varρ (τj ) ≍ m
The variance could be further reduced by repeating the measurement of the observable X K times (for a system identically prepared in state ρ) and averaging the
outcomes of the resulting K measurements. In this case, the response variable becomes
2
ν
Y = hρ, Xi + ξ, where Eρ (ξ|X) = 0 and Eρ (ξ 2 |X) = Varρ (Y |X) = 1−α
Km .

2.3. Minimax lower bounds. In [20], the problem of density matrix estimation
was studied in the case of trace regression model (1.1) with i.i.d. random design variables X1 , . . . , Xn sampled from the uniform distribution in an orthonormal basis E =
{E1 , . . . , Em2 } in two different settings: trace regression with Gaussian noise and trace
regression with a bounded response. In both cases, minimax lower bounds on the estimation error of the unknown target density matrix ρ of rank at most r were obtained
for the Schatten p-norm distances (p ∈ [1, +∞]) as well as for the Bures version of quantum Hellinger distance and for the quantum Kullback-Leibler (relative entropy) distance.
These results of [20] are stated below.
Denote by Sr,m the set of all density matrices of rank at most r (1 ≤ r ≤ m).
Assumption 1 (Trace regression with Gaussian noise). Let (X, Y ) be a random
couple with X being a random matrix sampled from the uniform distribution Π in an

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

7

orthonormal basis E = {E1 , . . . , Em2 } ⊂ Hm . Suppose that, for some density matrix
ρ ∈ Sm , Y = hρ, Xi + ξ, where ξ is a mean zero normal random variable with variance
σξ2 independent of X. Let (X1 , Y1 ), . . . , (Xn , Yn ) be n i.i.d. copies of (X, Y ).
In this model, the level of the noise ξ is characterized by its variance which should
be involved in the error bound (this could be viewed as a normal approximation of the
noise in the case when repeated measurements are performed for each observable Xj
with averaging of the outcomes).
Theorem 1. Suppose Assumption 1 holds. For all p ∈ [1, +∞], there exist constants c, c′ > 0 such that, the following bounds hold:1
(2.3)

(2.4)





inf sup Pρ kρ̂ − ρkp ≥ c r
ρ̂ ρ∈Sr,m

1/p σξ m

√

3
2

n

1
^ σξ m3/2 1− p ^ 

√
n

1

≥ c′ ,

3



σξ m 2 ^
inf sup Pρ H 2 (ρ̂, ρ) ≥ c r √
1
≥ c′ ,
ρ̂ ρ∈Sr,m
n

and
(2.5)

inf sup Pρ
ρ̂ ρ∈Sr,m

3

σξ m 2 ^
1
≥ c′ ,
K(ρkρ̂) ≥ c r √
n





where inf ρ̂ denotes the infimum over all estimators ρ̂ in Sm based on the data (X1 , Y1 ), . . . , (Xn , Yn )
satisfying the Gaussian trace regression model with noise variance σξ2 .
The trace regression model with a bounded response is characterized by the size
U of the range of response variable Y, which usually coincides with the bound on the
operator norms of the basis matrices Ej . It includes, in particular, Pauli measurements
discussed above (for which U = m−1/2 ).
Assumption 2 (Trace regression with a bounded response). Let (X, Y ) be a random couple with X being a random matrix sampled from the uniform distribution Π in
an orthonormal basis E = {E1 , . . . , Em2 } ⊂ Hm with U := max1≤j≤m2 kEj k∞ and Y
being a random variable with values in the interval [−U, U ]. Suppose that, for some density matrix ρ ∈ Sm , E(Y |X) = hρ, Xi a.s. Let (X1 , Y1 ), . . . , (Xn , Yn ) be n i.i.d. copies of
(X, Y ).
Let Pr,m (U ) denote the class of all distributions P of (X, Y ) such that Assumption
2 holds for some U > 0 and E(Y |X) = hρP , Xi for some ρP ∈ Sr,m . For a given P ∈
Pr,m (U ), PP denotes the corresponding probability measure such that (X1 , Y1 ), . . . , (Xn , Yn )
are i.i.d. copies of (X, Y ).
1

Here Pρ denotes a probability measure such that Assumption 1 is satisfied with density matrix ρ.

8

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

Theorem 2.

Suppose Assumption 2 is satisfied and, for some constant γ ∈ (0, 1),
tr(Ek ) ≤ (1 − γ)U m, k = 1, . . . , m2 .

(2.6)

Then, for all p ∈ [1, +∞], there exist constants cγ , c′γ > 0 such that the following bounds
hold:
(2.7)

(2.8)

inf

sup





PP kρ̂ − ρP kp ≥ cγ r

ρ̂ P ∈Pr,m (U )

inf

sup

ρ̂ P ∈Pr,m (U )



1/p U m

√

2

PP H (ρ̂, ρP ) ≥ cγ

3
2

n

1
^ U m3/2 1− p ^ 

√
n

1

≥ c′γ ,

3

Um2 ^
r √
1
≥ c′γ ,
n



and
(2.9)

inf

sup

ρ̂ P ∈Pr,m (U )



PP K(ρP kρ̂) ≥ cγ

3

Um2 ^
1
≥ c′γ ,
r √
n



where inf ρ̂ denotes the infimum over all estimators ρ̂ in Sm based on the i.i.d. data
(X1 , Y1 ), . . . , (Xn , Yn ) sampled from P.
As it was pointed out in [20] (see Remark 12), if γ in condition (2.6) is small enough
(say, γ < 1 − √12 ), then, in a given orthonormal basis E, there exists at most one matrix
Ej such that tr(Ej ) > (1 − γ)U m. In the case of Pauli basis, such a matrix indeed exists
and it is E1 = W0 ⊗ · · · ⊗ W0 . Thus, Theorem 2 does not apply directly to the Pauli
measurement model. However, the following result does hold (see [20], Theorem 10).
Theorem 3. Let {E1 , . . . , Em2 } be the Pauli basis in the space Hm of m × m Hermitian matrices and let X1 , . . . , Xn be i.i.d. random variables sampled from the uniform
distribution in {E1 , . . . , Em2 }. Let Y1 , . . . , Yn be outcomes of measurements of observables
X1 , . . . , Xn for the system being identically prepared n times in state ρ. The corresponding probability measure will be denoted by Pρ . Then, for all p ∈ [1, +∞], there exist
constants c, c′ > 0 such that the following bounds hold:
(2.10)

(2.11)





inf sup Pρ kρ̂ − ρkp ≥ c r
ρ̂ ρ∈Sr,m

inf sup Pρ
ρ̂ ρ∈Sr,m

1/p

m ^ m
√
√
n
n


1− p1 ^ 

1

m ^
H (ρ̂, ρ) ≥ c r √
1
n



≥ c′ ,

m ^
1
K(ρkρ̂) ≥ c r √
n



≥ c′ ,



2



≥ c′ ,

and
(2.12)

inf sup Pρ
ρ̂ ρ∈Sr,m





where inf ρ̂ denotes the infimum over all estimators ρ̂ in Sm based on the data (X1 , Y1 ), . . . , (Xn , Yn ).

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

9

It was also shown in [20] that, in the case of Schatten p-norm distances for p ∈ [1, 2],
Bures distance and Kullback-Leibler distance, the minimax lower bounds of theorems
1, 2 and 3 are attained up to logarithmic factors in m and n for a penalized least
squares estimator with von Neumann entropy penalty introduced in [16]. In the current
paper, our main goal is to show that the minimax optimal rates are attained up to
logarithmic factors for a very simple minimal distance estimator (that does not require
any penalization) in the whole range of Schatten p-norms, p ∈ [1, ∞], as well as for Bures
and Kullback-Leibler distances.
3. Main Results. For the model of uniform sampling from an orthonormal basis
E = {E1 , . . . , Em2 }, the following simple estimator of unknown state ρ ∈ Sm is unbiased:
Ẑ :=

n
m2 X
Yj Xj .
n j=1

Indeed,
Eρ Ẑ = m2 Eρ (Y X) = m2 E(Eρ (Y |X)X) = m2 Etr(ρX)X
2

m
1 X
hρ, Ej iEj = ρ.
= m Ehρ, XiX = m 2
m j=1
2

2

Clearly, Ẑ is not necessarily a density matrix.
We will now define the minimal distance estimator ρ̌ as the projection of Ẑ onto the
convex set Sm of all density matrices. More precisely, for an arbitrary Z ∈ Hm , define
πSm (Z) := argminS∈Sm kZ − Sk22 .

(3.1)

Clearly, πSm (Z) is the closest density matrix to Z with respect to the Hilbert–Schmidt
norm distance (that is, the projection of Z onto Sm ; such a closest density matrix exists in
view of compactness of Sm and it is unique in view of strict convexity of S 7→ kZ − Sk22 ).
Let
ρ̌ := πSm (Ẑ).
Remark 1.

This definition is equivalent to the following


ρ̌ := argminS∈Sm − n2
(3.2)

argminS∈Sm



1
n

Pn

Pn

2
j=1 Yj



−2
2
j=1 Yj hS, Xj i + m kSk2 =

−

2
n

Pn

j=1 Yj hS, Xj i +

m−2 kSk22



that was considered in [17] (in [19], similar estimators involving nuclear norm penalty
were studied). Note that replacing the term m−2 kSk22 in the right hand side of (3.2) by
P
its unbiased estimator n−1 nj=1 hS, Xj i2 yields the usual least squares estimator

(3.3)



ρ̂ := argminS∈Sm n

−1

n
X

2



(Yj − hS, Xj i) .

j=1

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

10

Note that we also have


ρ̂ := argminS∈Sm n

(3.4)

−1

n
X

2

(Yj − hS, Xj i) + εkSk1

j=1



since, for S ∈ Sm , kSk1 = tr(S) = 1. Thus, ρ̂ coincides with the nuclear norm penalized
least squares estimator (also called the matrix LASSO estimator) for any value of the
regularization parameter ε.
We will show that the upper bounds on the error rates in Schatten p-norm distances
for p ∈ [1, ∞] and in Bures distance that match the minimax lower bounds of theorems
1, 2 and 3 up to logarithmic factors hold for the estimator ρ̌. We will then introduce a
simple modification of this estimator for which a matching upper bound holds also for
Kullback-Leibler distance.
First, we consider the case of Gaussian trace regression model (Assumption 1). We
need an additional assumption that σξ ≥ mU1/2 (the variance of the noise is not too small).
Theorem 4. Suppose Assumption 1 holds and σξ ≥ mU1/2 . For all p ∈ [1, +∞],
there exists a constant C > 0 such that, for all A ≥ 1 the following bounds hold:
(3.5)
p
3p



 1

2
A log(2m) ^ σξ m3/2 A log(2m) 1− p ^
1/p σξ m
√
√
sup Pρ kρ̌−ρkp ≥ C r
1
≤ (2m)−A
n
n
ρ∈Sr,m
and


2

sup Pρ H (ρ̌, ρ) ≥ c r

(3.6)

ρ∈Sr,m

If σξ <



U
,
m1/2

3

σξ m 2

p

A log(2m) ^
√
1
n



the bounds still hold with σξ replaced by

≤ (2m)−A .

U
.
m1/2

Similarly, in the case of trace regression with a bounded response, the following
result holds.
Theorem 5. Suppose Assumption 2 is satisfied. Then, for all p ∈ [1, +∞], there
exists a constant C > 0 such that, for all A ≥ 1, the following bounds hold:
(3.7)
p
1
3p


^ U m3/2 A log(2m) 1− p ^ 
1/p U m 2 A log(2m)
√
√
sup PP kρ̌−ρP kp ≥ C r
1
≤ (2m)−A
n
n
P ∈Pr,m (U )
and
(3.8)

sup
P ∈Pr,m (U )



2



PP H (ρ̌, ρP ) ≥ C r

3

Um2

p

A log(2m) ^
√
1
n



≤ (2m)−A .

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

11

For completeness, we state also the upper bounds in the case of Pauli measurements
(that immediately follow from Theorem 5).
Theorem 6. Suppose the assumptions of Theorem 3 hold. Then, for all p ∈
[1, +∞], there exists a constant C such that, for all A ≥ 1, the following bounds hold:
(3.9)
p
p
1


^ m A log(2m) 1− p ^ 
1/p m A log(2m)
√
√
sup Pρ kρ̌ − ρkp ≥ c r
1
≤ (2m)−A
n
n
ρ∈Sr,m
and
(3.10)

sup Pρ

ρ∈Sr,m

m ^
1
H (ρ̌, ρ) ≥ c r √
n







2

≤ (2m)−A .

The proof of these results relies on the following fact that might be of independent
interest and that essentially shows that πSm (Z) is the closest density matrix to Z not
only in the Hilbert–Schmidt norm distance, but also in the operator norm distance.
Theorem 7.

For all Z ∈ Hm ,
kZ − πSm (Z)k∞ = inf kZ − Sk∞ .
S∈Sm

The proof of this theorem will be given in Section 4. Here we use it to establish the
next result that is the main ingredient of the proofs of theorems 4, 5 and 6.
Theorem 8.

Let p ∈ [1, +∞]. For all Z ∈ Hm and all S ∈ Sr,m ,




kπSm (Z) − Skp ≤ min 23/p+1 r 1/p kZ − Sk∞ , 2kZ − Sk1−1/p
.
∞
The proof relies on Theorem 7 and on a simple lemma stated below.
Lemma 1.

Let S, S ′ ∈ Sm and rank(S) = r. Then, for all p ∈ [1, ∞],




kS ′ − Skp ≤ min (8r)1/p kS ′ − Sk∞ , 21/p kS ′ − Sk1−1/p
.
∞
Proof. Let S = rj=1 λj (φj ⊗φj ) be the spectral decomposition of S with eigenvalues λj and eigenvectors φj . Let L := supp(S) be the linear span of vectors φ1 , . . . , φr ∈
Cm . Denote by PL , PL⊥ the orthogonal projection operators onto subspace L and its orthogonal complement L⊥ , respectively. We will need the following projection operators
PL , PL⊥ : Hm 7→ Hm :
P

PL⊥ (A) = PL⊥ APL⊥ , PL (A) = A − PL⊥ APL⊥ , A ∈ Hm .

12

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

The following bounds are obvious:
kSk1 = 1 = kS ′ k1 = kS ′ − S + Sk1 = kPL (S ′ − S) + PL⊥ (S ′ − S) + Sk1
≥ kPL⊥ (S ′ − S) + Sk1 − kPL (S ′ − S)k1 .
Since S = PL SPL , we can use the pinching inequality for unitary invariant norm k · k1
(see [4], p. 97) to get:
kPL⊥ (S ′ − S) + Sk1 = kPL SPL + PL⊥ (S ′ − S)PL⊥ k1
= kPL SPL k1 + kPL⊥ (S ′ − S)PL⊥ k1 = kSk1 + kPL⊥ (S ′ − S)k1 .
Therefore,
kSk1 ≥ kSk1 + kPL⊥ (S ′ − S)k1 − kPL (S ′ − S)k1 ,
implying that
kPL⊥ (S ′ − S)k1 ≤ kPL (S ′ − S)k1 .
It follows from the last bound that
kS ′ − Sk1 = kPL (S ′ − S) + PL⊥ (S ′ − S)k1 ≤ 2kPL (S ′ − S)k1 .
Since dim(L) = r, the matrix PL (S ′ − S) is of rank at most 2r. This implies that
kPL (S ′ − S)k1 ≤ 2rkPL (S ′ − S)k∞
≤ 2r(k(S ′ − S)PL k∞ + kPL (S ′ − S)PL⊥ k∞ ) ≤ 4rkS ′ − Sk∞ .

Therefore, kS ′ − Sk1 ≤ 8rkS ′ − Sk∞ , and since also kS ′ − Sk1 ≤ 2, S, S ′ ∈ Sm , we
conclude that
kS ′ − Sk1 ≤ min(8rkS ′ − Sk∞ , 2).
Together with interpolation inequality this yields that for all p ∈ [1, ∞]


1/p



1−1/p
≤ min (8r)1/p kS ′ − Sk∞ , 21/p kS ′ − Sk1−1/p
kS ′ − Skp ≤ kS ′ − Sk1 kS ′ − Sk∞
.
∞

Proof. We now prove Theorem 8. It immediately follows from Theorem 7 that,
for all S ∈ Sm ,
kπSm (Z) − Sk∞ ≤ kπSm (Z) − Zk∞ + kZ − Sk∞ ≤ 2kZ − Sk∞ .
If S ∈ Sm is a density matrix of rank r, the last bound could be combined with the
bound of Lemma 1 to get that for all p ∈ [1, +∞]




1−1/p
kπSm (Z) − Skp ≤ min 23/p+1 r 1/p kZ − Sk∞ , 2kZ − Sk∞
.

13

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

Proof. We now turn to the proof of theorems 4, 5 and 6. To this end, we use the
bound of Theorem 8 with Z = Ẑ and S = ρ ∈ Sr,m that yields:




1−1/p
kρ̌ − ρkp ≤ min 23/p+1 r 1/p kẐ − ρk∞ , 2kẐ − ρk∞
.

(3.11)
The control of

n
m2 X
Yj Xj − ρ
n j=1

kẐ − ρk∞ =

∞

is based on a standard application of matrix Bernstein type inequalities. We give a
detailed argument for completeness. Note that kρ̌ − ρkp in the left-hand side of bound
(3.11) is upper bounded by 2, so, if Bernstein bound on kẐ − ρk∞ is larger than 1 (or
even & 1), it could be replaced by the trivial bound equal to 1. In the case of Theorem 5,
we use the following version of Bernstein inequality for i.i.d. bounded random matrices
(see, e.g., [28]).
Lemma 2. Let V, V1 , . . . , Vn be i.i.d. random matrices in Hm with EV = 0. Suppose
that, for some constant U > 0, kV k∞ ≤ U a.s. Let σ 2 := kEV 2 k∞ . Then, for all t > 0
with probability at least 1 − e−t ,
V1 + · · · + Vn
n

∞



≤2 σ

s

t + log(2m) _ t + log(2m)
U
.
n
n


For V = Y X − E(Y X), we get, under Assumption 2, that
σ 2 = kEV 2 k∞ ≤ kE(Y 2 X 2 )k∞ ≤ U 2 kEX 2 k∞ .
It is also well known that, under the same assumption, kEX 2 k∞ = m−1 . [Indeed, if
{ej , j = 1, . . . , m} is an orthonormal basis of Cm , then
kEX 2 k∞ =

=

sup
v∈Cm ,|v|≤1

E

sup

EhX 2 v, vi =

v∈Cm ,|v|≤1

m
X

j=1

E|Xv|2 =

sup
v∈Cm ,|v|≤1

2

|hX, v⊗ej i| =

sup
v∈Cm ,|v|≤1

sup
v∈Cm ,|v|≤1

m−2

m
X

j=1

m
X

j=1

sup

E

v∈Cm ,|v|≤1

m
X

j=1

|hXv, ej i|2

2

m

−2

m
X

k=1

2

|hEk , v⊗ej i| =

sup
v∈Cm ,|v|≤1

m

−2

m
X

j=1

|v|2 |ej |2 = m−1 ].

We use the bound of Lemma 2 with t = A log(2m), A ≥ 1 to get that with probability
at least 1 − (2m)−A ,
n
m2 X
Yj Xj − ρ
n j=1

∞



≤ C Um

3/2

s

A log(2m) _ U 2 m2 A log(2m)
n
n



kv⊗ej k22 =

14

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

with some absolute constant C ≥ 1. If
U 2 m2 A log(2m)
≥ U m3/2
n
q

s

A log(2m)
,
n

q

then U m1/2 A log(2m)
≥ 1 implying that U m3/2 A log(2m)
≥ 1. Thus, when the bound
n
n
on kẐ q
− ρk∞ is substituted in bound (3.11), it is enough to keep only the first term
, the second term could be dropped. This implies that with some conU m3/2 A log(2m)
n
′
stant C > 0 (that does not depend on ρ ∈ Sr,m ) the inequality
kρ̌ − ρkp ≤ C

′



r

3

1/p U m 2

A log(2m) ^ U m3/2 A log(2m)
√
√
n
n

p

p



1− p1 ^ 

1

holds with probability at least 1 − (2m)−A , implying the first bound of Theorem 5. The
second bound immediately follows from the inequality H 2 (ρ̌, ρ) ≤ kρ̌ − ρk1 (see (2.2)).
Theorem 6 is an immediate consequence of Theorem 5.
The proof of Theorem 4 is very similar. In this case, Assumption 1 holds and it is
natural to split Ẑ − ρ into two parts
(3.12)

Ẑ − ρ =

n
n
m2 X
m2 X
hρ, Xj iXj − ρ +
ξj Xj .
n j=1
n j=1

and to bound kẐ − ρk∞ by triangle inequality. For the first part, an application of matrix
Bernstein inequality of Lemma 2 yields the bound
(3.13)

n
m2 X
hρ, Xj iXj − ρ
n j=1

∞



≤ C Um

s

A log(2m) _ U 2 m2 A log(2m)
n
n



that holds for some absolute constant C ≥ 1 with probability at least 1−(2m)−A . Indeed,
in this case V = hρ, XiX − Ehρ, XiX and
σ 2 ≤ kEhρ, Xi2 X 2 k∞ ≤ U 2 Ehρ, Xi2 =

U 2 kρk22
U2
≤
,
m2
m2

khρ, XiXk∞ ≤ kρk1 kXk2∞ ≤ kXk2∞ ≤ U 2 ,
and Lemma 2 implies (3.13). As before, if
U 2 m2 A log(2m)

U 2 m2 A log(2m)
n

≥ Um

q

A log(2m)
,
n

then U m

q

A log(2m)
n

could be dropped when the bound on kẐ − ρk∞
1. Thus, the second term
n
(for which the right hand side of (3.13) is a part) is substituted in (3.11).
As to the second part of representation (3.12) that involves normal random variables
ξj , it is bounded using another version of matrix Bernstein inequality for not necessarily
bounded random matrices (see [16], [15], [18]).

≥

15

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

Lemma 3. Let V, V1 , . . . , Vn be i.i.d. random matrices in Hm with EV = 0. Suppose
that, for some α ≥ 1, U (α) := 2 kV k∞ ψα < +∞. 2 Let σ 2 := kEV 2 k∞ . Then, for all
t > 0 with probability at least 1 − e−t ,
V1 + · · · + Vn
n

∞



≤C σ

s

U (α) t + log(2m)
t + log(2m) _ (α)
U log1/α
.
n
σ
n






We apply the bound of Lemma 3 in the case when V := ξX, α = 2 and t = A log(2m)
for A ≥ 1. By an easy computation,
2

σ =

σξ2
=
m

σξ2 kEX 2 k∞

and
U (2) = 2 ξkXk∞

ψ2

≤ 2U kξkψ2 ≤ 4σξ U.

This yields the following bound
(3.14)

n
m2 X
ξj Xj
n j=1

∞



≤ C σξ m3/2

s

√ 
m2 A log(2m) log1/2 (4U m)
A log(2m) _
σξ U
n
n

that holds with probability at least 1 − (2m)−A and with some absolute constant C ≥ 1.
If the second
term in the maximum in the right hand side of (3.14) is dominant, then
q
√
A log(2m)
1/2
log1/2 (4U m) ≥ 1. Under the condition that σξ ≥ U m−1/2 , this implies
Um
n
q

& 1. Thus, when the bound in the right hand side of (3.14)
that also σξ m3/2 A log(2m)
n
(used to control kẐ − ρk∞ ) is substituted in (3.11), it is enough to keep only the first
term in the maximum. Finally, under the assumption σξ ≥ U m−1/2 , the first term of
bound (3.14) dominates the first term of (3.13), so, only this term is needed to control
kẐ − ρk∞ in bound (3.11). These considerations imply the bound
kρ̌ − ρkp ≤ C

′



r

1/p σξ m

3
2

A log(2m) ^ σξ m3/2 A log(2m)
√
√
n
n

p

p



1− p1 ^ 

1

that holds with some constant C ′ > 0 (that does not depend on ρ ∈ Sr,m ) and with
probability at least 1 − (2m)−A . The first bound of Theorem 4 now follows for all p ∈
[1, ∞] (which also implies the second bound in view of (2.2)).
2

Here k · kψα denotes the ψα Orlicz norm in the space of random variables defined as follows:



kηkψα := inf c > 0 : E exp

n |η|α o
cα



≤2 .

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

16

It turns out that for a slightly modified version of estimator ρ̌, minimax lower
bounds are also attained (up to logarithmic factors) in the case of Kullback-Leibler
m
. Clearly, Sδ ∈ Sm . Let
distance. For S ∈ Sm and δ ∈ [0, 1], define Sδ = (1 − δ)S + δ Im
Sm,δ := {Sδ : S ∈ Sm }. Define πSm,δ (Z) the projection of Z ∈ Hm onto the convex set
Sm,δ :
πSm,δ (Z) := argminS∈Sm,δ kZ − Sk22 .
Let
ρ̌δ := πSm,δ (Ẑ)
with ρ̌0 = ρ̌. We will prove the following versions of theorems 4, 5 and 6 for the estimator
ρ̌δ .
Suppose Assumption 1 holds, σξ ≥

Theorem 9.

3

δ≤

σξ m 2

U
m1/2

and

p

log(2m) ^
√
1.
n

Then bounds (3.5) and (3.6) hold for estimator ρ̌δ . Moreover, for A ≥ 1, define
λ :=

rσξ m5/2

q

A log(2m)
n

δ

V

m

.

Then, for some constant c > 0,
(3.15)

3



sup Pρ K(ρkρ̌δ ) ≥ c r

ρ∈Sr,m

If σξ <



U
,
m1/2

p

A log(2m) ^
√
1 log(1 + cλ) ≤ (2m)−A .
n

σξ m 2



the bounds still hold with σξ replaced by

Theorem 10.



U
.
m1/2

Suppose Assumption 2 is satisfied and
3

δ≤

Um2

p

log(2m) ^
√
1.
n

Then (3.7) and (3.8) hold for estimator ρ̌δ . Moreover, for A ≥ 1, define
λ :=

rU m5/2

q

A log(2m)
n

δ

V

m

.

Then, for some constant c > 0,
(3.16)





sup Pρ K(ρkρ̌δ ) ≥ c r

ρ∈Sr,m

3

Um2

p

A log(2m) ^
√
1 log(1 + cλ) ≤ (2m)−A .
n




D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

Theorem 11.

17

Suppose the assumptions of Theorem 3 hold and
p

m log(2m) ^
√
1.
δ≤
n
Then (3.9) and (3.10) hold for estimator ρ̌δ . Moreover, for A ≥ 1, define
λ :=

rm2

q

A log(2m)
n

δ

V

m

.

Then, for some constant c > 0,
(3.17)

sup Pρ

ρ∈Sr,m

Remark 2.

p

m A log(2m) ^
√
K(ρkρ̌δ ) ≥ c r
1 log(1 + cλ) ≤ (2m)−A .
n









If, under the assumptions of Theorem 10, we choose
3

δ=

Um2

p

log(2m) ^
√
1,
n

then the logarithmic factor in bound (3.16) satisfies the inequality
√
log(1 + cλ) ≤ log(1 + crm A),
so it is of the order log m. Under the assumptions of Theorem 9, this would require the
choice of δ
3p
σξ m 2 log(2m) ^
√
δ=
1,
n

so δ would depend on an unknown parameter σξ . Replacing σξ in the definition of δ by


√
σξ
the lower bound U m−1/2 would result in a logarithmic factor . log 1 + crm A U m−1/2
.

Proof. We start with the following modification of Theorem 8.
Lemma 4.
holds:

Let p ∈ [1, ∞]. For all Z ∈ Hm and all S ∈ Sr,m , the following bound








kπSm,δ (Z)−Skp ≤ min 23/p+1 r 1/p kZ−Sk∞ +2δ , 2(1−δ)1/p kZ−Sk∞ +2δ

1−1/p 

+2δ.

18

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

Proof. The following formula is straightforward: for δ ∈ [0, 1),
πSm,δ (Z) = (1 − δ)πSm



δ Im
Z
−
1−δ 1−δ m



+δ

Im
.
m

′

Indeed, πSm,δ (Z) coincides with (1 − δ)S + δ Imm , where
S ′ := argminS∈Sm Z − (1 − δ)S − δ
= argminS∈Sm

Z
δ Im
−
−S
1−δ 1−δ m

2
2

= πS m



implying the claim.
Let S ∈ Sr,m . Then, for p ∈ [1, ∞],
(3.18)

Im
m

2
2

Z
δ Im
−
,
1−δ 1−δ m


kπSm,δ (Z) − Skp ≤ kπSm,δ (Z) − Sδ kp + kSδ − Skp
≤ (1 − δ) πSm



Z
1−δ

−

δ Im
1−δ m



−S

+ 2δ.
p

To control the first term in the right hand side, we use the bound of Theorem 8, which
Z
δ Im
requires bounding 1−δ
− 1−δ
. We have
m −S
∞

Z
1−δ

(3.19)
≤

1
1−δ kZ

−

δ Im
1−δ m

− Sk∞ +

−S

1
1−δ kS

=
∞

1
1−δ kZ

− Sδ k∞ ≤

− Sδ k∞

1
1−δ kZ

− Sk∞ +

2δ
1−δ .

Using bounds (3.18), (3.19) along with the bound of Theorem 8, we get the bound of
the lemma.

We will use the bound of Lemma 4 to control kρ̌δ − ρkp for ρ ∈ Sr,m . To this end, we
need to bound kẐ −ρk∞ using matrix Bernstein inequalities exactly as it was done in the
¯ such
proof of theorems 4, 5 and 6 (under assumptions of these theorems). Denote by ∆
−A
an upper bound
that
q
q on kẐ − ρk∞ that holds with probability a least 1 − (2m) . Recall
A
log(2m)
A
log(2m)
¯ ≍ U m3/2
¯ ≍ σξ m3/2
under the conditions of Theorem 4 and ∆
∆
n

n

under the conditions of Theorem 5 (it is the same under the conditions of Theorem 6
¯ ∧ 1, we get from the bound of Lemma 4 that
with U = m−1/2 ). Setting ∆ = ∆


3/p+1 1/p

kρ̌δ − ρkp ≤ min 2

r





1/p

∆ + 2δ , 2(1 − δ)



∆ + 2δ

1−1/p 

+ 2δ

¯ by
that holds with the same probability at least 1 − (2m)−A . Recall that we replace ∆
∆ since the left hand side kρ̌δ − ρkp ≤ 2; for the same reason, we can and do drop the

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

19

“exponential parts” of matrix Bernstein bounds leaving in the definition of ∆ only the
“Gaussian parts”. For δ . ∆, we get
kρ̌δ − ρkp . min(r 1/p ∆, ∆1−1/p ).
Exactly as in the proof of theorems 4, 5 and 6, this implies that bounds (3.5), (3.6),
(3.7), (3.8), (3.9) and (3.10) hold for estimator ρ̌δ .
The bound on the Kullback-Leibler divergence K(ρkρ̌δ ) is an immediate consequence
of the bound on kρ̌δ − ρk1 and the next lemma that follows from Corollary 1 in [3].
Lemma 5. Let S1 , S2 ∈ Sm be density matrices and let β := λmin (S2 ) be the
smallest eigenvalue of S2 . Suppose that β > 0. Then


K(S1 kS2 ) ≤ kS1 − S2 k1 log 1 +

kS1 − S2 k1
.
2β


We apply Lemma 5 to S1 = ρ, S2 = ρ̌δ , observing that ρ̌δ ∈ Sm,δ and λmin (ρ̌δ ) ≥
δ/m. We then use the bound on kρ̌δ −ρk1 to complete the proof of the bound on K(ρkρ̌δ ).
We conclude this section with a simple result concerning the least squares estimator
ρ̂ defined by (3.3). It shows that the estimators ρ̂ and ρ̌ are close in the Hilbert-Schmidt
norm. As a result, the bounds of the previous theorems could be applied to estimator ρ̂
as well (at least, under some additional assumptions).
Theorem 12. Under the assumption that i.i.d. design variables X1 , . . . , Xn are
sampled from the uniform distribution Π in an orthonormal basis E = {E1 , . . . , Em2 },
the following bound holds with some constant C > 0 for all A ≥ 1 with probability at
least 1 − (2m2 )−A :
s
kρ̌ − ρ̂k2 ≤ Cm

A log(2m)
.
n

Proof. Note that the gradient (and subgradient) of convex function S 7→ kS − Zk22
is equal to 2(S − Z). By a necessary condition of minimum in convex minimization
problem (3.1), for ρ̌ = πSm (Ẑ), Ẑ − ρ̌ should belong to the normal cone NSm (ρ̌) of
the convex set Sm at point ρ̌ (see [2], Proposition 5, Chapter 4, Section 1). Since both
ρ̌, ρ̂ ∈ Sm , this implies that
(3.20)

hρ̌ − Ẑ, ρ̌ − ρ̂i ≤ 0.

Similar analysis of convex optimization problem (3.3) shows that



n
m2 X
(hρ̂, Xj i − Yj )Xj , ρ̌ − ρ̂ ≥ 0,
n j=1

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

20

which could be rewritten as follows:


(3.21)


n
m2 X
hρ̂, Xj iXj − Ẑ, ρ̌ − ρ̂ ≥ 0.
n j=1

Subtracting (3.21) from (3.20) yields


ρ̌ −


n
m2 X
hρ̂, Xj iXj , ρ̌ − ρ̂ ≤ 0,
n j=1

implying that
(3.22)

kρ̌ −

We will now write

ρ̂k22

= hρ̌ − ρ̂, ρ̌ − ρ̂i ≤




n
m2 X
hρ̂, Xj iXj − ρ̂, ρ̌ − ρ̂ .
n j=1

3

n
n 

m2 X
m2 X
hρ̂, Xj iXj − ρ̂ =
hρ̂, Xj iXj − Ehρ̂, XiX
n j=1
n j=1

=m

2



n
1X
(Xj ⊗ Xj − E(X ⊗ X)i ρ̂.
n j=1



It follows from (3.22) that
kρ̌ − ρ̂k22 ≤ m2

n
1X
Xj ⊗ Xj − E(X ⊗ X)
n j=1

op

kρ̂k2 kρ̌ − ρ̂k2 .

Since kρ̂k2 ≤ 1, we get
kρ̌ − ρ̂k2 ≤ m2

(3.23)

n
1X
Xj ⊗ Xj − E(X ⊗ X)
n j=1

.
op

It remains to control the operator norm in the right hand side for which we can again
use matrix Bernstein inequality of Lemma 2 applying it to V = X ⊗ X − E(X ⊗ X). In
this case,
σ 2 = kEV 2 kop ≤ kE(X⊗X)2 kop = sup Eh(X⊗X)2 U, U i = sup Eh(X⊗X)U, (X⊗X)U i
kU k2 ≤1

kU k2 ≤1

kU k22
1
= 2
2
m
kU k2 ≤1 m

= sup E|hU, Xi|2 kXk22 ≤ sup E|hU, Xi|2 = sup
kU k2 ≤1

kU k2 ≤1

3
Here we view the tensor product A ⊗ B of operators A, B ∈ Mm as an operator acting from the
space Mm of m × m matrices equipped with Hilbert-Schmidt inner product h·, ·i into itself as follows:
(A ⊗ B)C = AhC, Bi. Let k · kop denote the operator norm of linear operators from Mm into itself, which
corresponds to the k · k∞ in the case of m × m matrices.

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

21

and
kV kop ≤ kX ⊗ Xkop + EkX ⊗ Xkop = kXk22 + EkXk22 ≤ 2.

Bound (3.23) along with the bound of Lemma 2 with t = A log(2m2 ), A ≥ 1 yield the
following inequality
kρ̌ − ρ̂k2 . m

s

A log(2m) _ 2 A log(2m)
m
n
n

that holds with probability at least 1 − (2m2 )−A . Since kρ̌ − ρ̂k2 ≤ 2, the second term
in the right hand side could be dropped (if this term is dominant, the bound
m2 A log(2m)
n
is & 1). This completes the proof of the theorem.
Since kρ̌ − ρ̂k∞ ≤ kρ̌ − ρ̂k2 , the bound of Theorem 12 also holds for kρ̌ − ρ̂k∞ .
Combining this with the bound of Theorem 5 for p = ∞, it is easy to conclude that
under conditions of this theorem
kρ̂ − ρk∞ . U m3/2

s

A log(2m)
n

and that the last bound holds (with a proper choice of constant in relationship .) with
probability at least 1 − (2m)−A . In view of Lemma 1, this immediately implies that all
the bounds of Theorem 5 also hold for the least squares estimator ρ̂. In a special case of
Pauli measurements, this means that Theorem 6 holds for the estimator ρ̂. Concerning
Theorem 9, the same conclusion is true under the additional assumption that σξ ≥ m−1/2 .
Moreover, if ρ̂δ is the following modification of estimator ρ̂
(3.24)



ρ̂δ := argminS∈Sm,δ n−1

n
X



(Yj − hS, Xj i)2 ,

j=1

then the statements of theorems 9, 10 and 11 hold for the estimator ρ̂δ (in the case of
Theorem 9, under the additional assumption that σξ ≥ m−1/2 ).
4. Proof of Theorem 7. Recall that
πSm (Z) := argminS∈Sm kZ − Sk22 , Z ∈ Hm
defines the projection of Z onto Sm . The mapping Hm ∋ Z →
7 πSm (Z) ∈ Sm possesses a
d the set of all
couple of simple properties stated in the next proposition. Denote by Sm
diagonal density matrices.
Proposition 1.

1. For all m × m unitary matrices U,
πSm (U −1 ZU ) = U −1 πSm (Z)U, Z ∈ Hm .

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

22

d.
2. If D ∈ Hm is a diagonal matrix, then πSm (D) ∈ Sm

Proof. To prove the first claim, note that, by the unitary invariance of the Hilbert–
Schmidt norm,
kU −1 ZU − Sk22 = kU −1 (Z − U SU −1 )U k22 = kZ − U SU −1 k22 .
In addition, the mapping S 7→ U SU −1 is a bijection from the set Sm onto itself. This
immediately implies that
πSm (U −1 ZU ) = argminS∈Sm kZ − U SU −1 k22 = U −1 πSm (Z)U.
d
For an m × m matrix A = (aij )m
i,j=1 ∈ Hm , let A be the diagonal matrix with
diagonal entries aii , i = 1, . . . , m. It is easy to see that if A is a density matrix, then Ad
is also a density matrix. Moreover, it is also obvious that, for a diagonal matrix D,

kD − Ad k22 ≤ kD − Ak22 , A ∈ Sm ,
with a strict inequality if A is not diagonal. These observations immediately imply the
second claim.

We will now state and prove a vector version of Theorem 7 in which the role of the
set of density matrices Sm is played by the simplex
n

∆m := u = (u1 , . . . , um ) ∈ Rm : uj ≥ 0,

m
X

j=1

o

uj = 1

in Rm (this is equivalent to considering the set of diagonal density matrices). We will
then show that the matrix version of the problem reduces to the vector case.
Define
, z ∈ Rm .
π∆m (z) := argminu∈∆m kz − uk2ℓm
2
Since the function u 7→ kz − uk2ℓm is strictly convex and ∆m is a compact convex set,
2
such a minimizer exists and is unique. In other words, π∆m (z) is the projection of the
point z ∈ Rm onto simplex ∆m (the closest point to z in the set ∆m with respect to the
Euclidean ℓm
2 -distance). The next lemma shows that the same point also minimizes the
ℓm
-distance
from z to the simplex ∆m .
∞
Lemma 6.

For all z ∈ Rm ,
.
= min kz − vkℓm
kz − π∆m (z)kℓm
∞
∞
v∈∆m

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

23

Proof. Without loss of generality, assume that z = (z1 , . . . , zm ) ∈ Rm is a point
with z1 ≥ · · · ≥ zm . Denote
z̄j :=

z1 + · · · + zj
, j = 1, . . . , m.
j

Clearly, z̄1 = z1 and z̄j ≥ zj , j = 1, . . . , m. Let


k := max j ≤ m : z̄j ≤ zj +

1
.
j


Note that if k > 1, then, for all j < k, z̄j ≤ zj + 1j . Indeed,
z̄j =

kz̄k −

Pk

i=j+1 zi

j

≤

kzk + 1 − (k − j)zk
jzk + 1
1
1
=
= zk + ≤ zj + .
j
j
j
j

On the other hand, if k < m, then z̄k > zk+1 + k1 . Indeed, if z̄k ≤ zk+1 + k1 , then
z̄k+1 =

kz̄k + zk+1
kzk+1 + 1 + zk+1
1
≤
= zk+1 +
,
k+1
k+1
k+1

which would contradict the definition of k.
Let λ = (λ1 , . . . , λm ), where λj = zj − z̄k + k1 for j = 1, . . . , k and λj = 0 for
j = k + 1, . . . , m. Since z̄k ≤ zk + k1 ≤ zj + k1 for all j ≤ k, we have λj ≥ 0, j = 1, . . . , m
and
 X
k
m
k 
X
X
1
=
zj − kz̄k + 1 = 1.
λj =
zj − z̄k +
k
j=1
j=1
j=1
Thus, λ ∈ ∆m . It turns out that π∆m (z) = λ.
z − λ ∈ N∆m (λ), where

4

To prove this it is enough to show that

N∆m (λ) := {u ∈ Rm : hu, v − λi ≤ 0, v ∈ ∆m }
is the normal cone of the convex set ∆m at point λ (see, e.g., [2], Proposition 5, Chapter
4, Section 1). Let t := z̄k − k1 . Clearly, we have zk+1 < t ≤ zk if k < m and t ≤ zm if
k = m. For k = m, z − λ = (t, . . . , t) and
hz − λ, v − λi =

m
X
i=1

t(vi − λi ) = t

X
m
i=1

vi −

m
X
i=1



λi = 0

since v, λ ∈ ∆m . For k < m, note that
z − λ = (t, . . . t, zk+1 , . . . , zm )
4
The computation of the projection onto a simplex occurs in many applications and has been studied
before: see, e.g., [24] and [27]. See also [8], where an explicit expression for the projection was derived.
For completeness, we provide our version of the proof below.

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

and, for v ∈ ∆m ,
hz − λ, v − λi =
Using the facts that

Pm

i=1 vi

= 1 and

= −t

vi +

i=k+1

i=1

t(vi − λi ) +

Pk

i=1 λi

hz − λ, v − λi = t
m
X

k
X

X
k
i=1

m
X

zi vi .

i=k+1

= 1, we get

vi −

zi vi =

i=k+1

m
X

24

k
X



λi +

i=1

m
X

m
X

zi vi

i=k+1

(zi − t)vi ≤ 0,

i=k+1

where we also used that, for all i = k + 1, . . . , m, zi − t ≤ zk+1 − t ≤ 0 and vi ≥ 0. Thus,
z − λ ∈ N∆m (λ) and, by the uniqueness of the minimum, λ = π∆m (z).
Note that
kz − λkℓm
= max(|t|, |zk+1 |, . . . , |zm |).
∞
For any v ∈ ∆m ,
t = z̄k −

k
m
k
k
k
1X
1X
1X
1X
1X
1
zi −
vi ≤
zi −
vi =
(zi − vi ) ≤ kz − vkℓm
.
=
∞
k
k i=1
k i=1
k i=1
k i=1
k i=1

On the other hand,
zm ≥ zm − vm ≥ −kz − vkℓm
.
∞
Since
t = z̄k −

1
≥ zk+1 ≥ · · · ≥ zm ,
k

we conclude that, for all v ∈ ∆m ,
kz − λkℓm
≤ kz − vkℓm
.
∞
∞

We now turn to the proof of Theorem 7.
Proof. Any matrix Z ∈ Hm admits spectral representation Z = U −1 DU, where
D is the diagonal matrix with real entries d1 , . . . , dm on the diagonal and U is a unitary
m × m matrix. Let d = (d1 , . . . , dm ) ∈ Rm . Given v = (v1 , . . . , vm ) ∈ ∆m , the diagonal
matrix V with entries v1 , . . . , vm is a density matrix. This defines a bijection ∆m ∋ v 7→
d of all diagonal m×m density matrices.
V = J(v) between the simplex ∆m and the set Sm
d : kJ(v) − J(u)k2 = ku − vk2 , u, v ∈ ∆ .
Moreover, J is an isometry of ∆m and Sm
m
2
ℓm
2
We will now prove the following lemma.

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

25

Lemma 7. Let Z = U −1 DU with a unitary m × m matrix U and diagonal matrix
D with d = (d1 , . . . , dm ) ∈ Rm being the vector of its diagonal entries. Then
πSm (Z) = U −1 J(π∆m (d))U.
Proof. This is an immediate consequence of Proposition 1 and the following simple
fact:



2
2
argminA∈Sm
d kD − Ak2 = J argminv∈∆ kJ(d) − J(v)k2
m



J argminv∈∆m kd − vk2ℓm
2





= J(π∆m (d)).

To complete the proof of Theorem 7, observe that, In view of lemmas 6, 7,
kZ − πSm (Z)k∞ = kU −1 (J(d) − J(π∆m (d)))U k∞
.
= inf kd − vkℓm
= kJ(d) − J(π∆m (d))k∞ = kd − π∆m (d)kℓm
∞
∞
v∈∆m

Without loss of generality, assume that d1 ≥ · · · ≥ dm . Let S ∈ Sm be a density matrix
with eigenvalues v1 ≥ · · · ≥ vm . Clearly, v = (v1 , . . . , vm ) ∈ ∆m . Therefore,
kZ − πSm (Z)k∞ ≤ kd − vk∞ ≤ kZ − Sk∞ ,
where to get the last bound we used Weyl’s perturbation inequality (see [4], Corollary
III.2.6).

5. Comments on computational aspects of the problem. An advantage of
minimal distance estimator ρ̌ = πSm (Ẑ) is the simplicity of its computational implemen2 Pn
2
tation. The computation of the matrix Ẑ = mn
i=1 Yi Xi requires O(nm ) operations.
It is followed by an eigen-decomposition of Z that requires O(m3 ) operations(see [10]);
there exist efficient software packages designed for this kind of tasks, for instance, LINPACK and PROPACK, etc.). As it is shown in the previous section, the problem of
computing πSm (Ẑ) then reduces to projecting of the vector of eigenvalues of Z arranged
in a non-increasing order onto the simplex ∆m . The last problem has been studied in the
literature (see [24], [27], [8]) and it has an explicit solution of computational complexity
proportional to m (see the proof of Lemma 6). Thus, the computational implementation
of the minimal distance estimator ρ̌ requires O((n + m)m2 ) operations.
The matrix version of LASSO estimator for density matrices is equivalent to solving
the following optimization problem
(5.1)

ρ̂ := arg min
S∈Sm

n 
2
1X
Yi − S, Xi
n i=1

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

26

that results in the least squares estimator. Clearly, there is no explicit solution for this
optimization problem and it is usually solved by iterative algorithms. For example, a
well know iterative singular value thresholding (SVT) algorithm was proposed in [5],
and also implemented in quantum compressed sensing in [9]. The main idea is that (5.1)
is equivalent to the following optimization problem: for any τ > 0,
ρ̂ :=

n 
2
m2 X
Yi − Z, Xi
+ τ kS − Zk22 .
S∈Sm ,Z∈Hm ,S=Z n i=1

arg min

The proposed algorithm updates Z and S alternatively, with the only constraint for S
being that S ∈ Sm . Therefore, the main ingredient of SVT is the following iterative
updating rule (with initial Z0 = 0): for k = 1, 2, . . .,
(5.2)

(

Sk = πSm (Zk−1 )
Zk = Sk + δk Ẑ −

m2
n

Pn

i=1

Sk , Xi Xi



with certain pre-determined step sizes δk > 0. The algorithm terminates at some step
k = N and outputs SN ∈ Sm when kSN − SN −1 k2 ≤ ǫ for some numerical threshold
ǫ > 0. It is clear that the minimal distance estimator ρ̌ can be produced by the above
algorithm with one iteration and the initialization Z0 = Ẑ, δ1 = 0. When the number
of qubits k is not small (for instance, about 20) and the dimension m is very large, the
iterative algorithm (5.2) is much more computationally expensive than the algorithm for
the minimal distance estimator (since every iteration requires the eigen-decomposition
of a high dimensional matrix).
References.
[1] P. Alquier, C. Butucea, M. Hebiri, K. Meziani, and T. Morimae. Rank penalized estimation of a
quantum system. Physical Reviews A, 88:032113, 2013.
[2] J.-P. Aubin and I. Ekeland. Applied Nonlinear Analysis. Courier Corporation, 2006.
[3] K. Audenaert and J. Eisert. Continuity bounds on the quantum relative entropy - ii. Journal of
Mathematical Physics, 52(112201), 2011.
[4] R. Bhatia. Matrix Analisis. Springer, 1997.
[5] J.-F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding algorithm for matrix completion.
SIAM Journal on Optimization, 20(4):1956–1982, 2010.
[6] E. J. Candès and Y. Plan. Tight oracle inequalities for low-rank matrix recovery from a minimal
number of noisy random measurements. IEEE Transactions on Information Theory, 57(4):2342–
2359, 2011.
[7] E. J. Candès and T. Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE
Transactions on Information Theory, 56(5):2053–2080, 2010.
[8] Y. Chen and X. Ye. Projection onto a simplex. arXiv preprint arXiv:1101.6081, 2011.
[9] S. T. Flammia, D. Gross, Y.-K. Liu, and J. Eisert. Quantum tomography via compressed sensing:
error bounds, sample complexity and efficient estimators. New Journal of Physics, 14(9):095022,
2012.
[10] G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
[11] D. Gross. Recovering low-rank matrices from few coefficients in any basis. IEEE Transactions on
Information Theory, 57(3):1548–1566, 2011.

D. XIA AND V. KOLTCHINSKII/ESTIMATION OF DENSITY MATRICES

27

[12] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert. Quantum state tomography via
compressed sensing. Physical Review Letters, 105(15):150401, 2010.
[13] H. Klauck, A. Nayak, A. Ta-Shma, and D. Zuckerman. Interaction in quantum communication.
IEEE Transactions on Information Theory, 53(6):1970–1982, 2007.
[14] O. Klopp. Noisy low-rank matrix completion with general sampling distribution. Bernoulli,
20(1):282–303, 2014.
[15] V. Koltchinskii. Oracle Inequalities in Empirical Risk Minimization and Sparse Recovery Problems:
École d’Été de Probabilités de Saint-Flour XXXVIII-2008. Springer, 2011.
[16] V. Koltchinskii. von Neumann entropy penalization and low-rank matrix estimation. The Annals
of Statistics, 39(6):2936–2973, 2011.
[17] V. Koltchinskii. A remark on low rank matrix recovery and noncommutative Bernstein type inequalities. In From Probability to Statistics and Back: High-Dimensional Models and Processes–A
Festschrift in Honor of Jon A. Wellner, pages 213–226. Institute of Mathematical Statistics, 2013.
[18] V. Koltchinskii. Sharp oracle inequalities in low rank estimation. In Empirical Inference, pages
217–230. Springer, 2013.
[19] V. Koltchinskii, K. Lounici, and A. B. Tsybakov. Nuclear-norm penalization and optimal rates for
noisy low-rank matrix completion. The Annals of Statistics, 39(5):2302–2329, 2011.
[20] V. Koltchinskii and D. Xia. Optimal estimation of low rank density matrices. Journal of Machine
Learning Research, 16(Sep):1757–1792, 2015.
[21] Y.-K. Liu. Universal low-rank matrix recovery from Pauli measurements. In Advances in Neural
Information Processing Systems, pages 1638–1646, 2011.
[22] K. Lounici. Optimal spectral norm rates for noisy low-rank matrix completion. arxiv:1110.5346,
2011.
[23] Z. Ma and Y. Wu. Volume ratio, sparsity, and minimaxity under unitarily invariant norms. IEEE
Transactions on Information Theory, 61(12):6939–6956, 2015.
[24] C. Michelot. A finite algorithm for finding the projection of a point onto the canonical simplex of
Rn . Journal of Optimization Theory and Applications, 50(1):195–200, 1986.
[25] S. Negahban and M. J. Wainwright. Estimation of (near) low-rank matrices with noise and highdimensional scaling. The Annals of Statistics, 39(2):1069–1097, 2011.
[26] M. Nielsen and I. Chuang. Quantum Computation and Quantum Information. Cambridge University
Press, 2000.
[27] S. Shalev-Shwartz and Y. Singer. Efficient learning of label ranking by soft projections onto polyhedra. Journal of Machine Learning Research, 7:1567–1599, 2006.
[28] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational
Mathematics, 12(4):389–434, 2012.
School of Mathematics, Georgia Institute of Technology, Atlanta, GA, 30332, USA.
E-mail: vlad@math.gatech.edu; dxia7@math.gatech.edu

