# Local Stochastic Factored Gradient Descent
for Distributed Quantum State Tomography

Junhyung Lyle Kim⋆, Mohammad Taha Toghani‡, César A. Uribe‡, Anastasios
Kyrillidis⋆ ⋆Department of Computer Science, Rice University, Houston, TX,
USA. {jlylekim<EMAIL_ADDRESS>of Electrical and Computer
Engineering, Rice University, Houston, TX, USA. {mttoghani<EMAIL_ADDRESS>
###### Abstract

We propose a distributed Quantum State Tomography (QST) protocol, named Local
Stochastic Factored Gradient Descent (Local SFGD), to learn the low-rank
factor of a density matrix over a set of local machines. QST is the canonical
procedure to characterize the state of a quantum system, which we formulate as
a stochastic nonconvex smooth optimization problem. Physically, the estimation
of a low-rank density matrix helps characterizing the amount of noise
introduced by quantum computation. Theoretically, we prove the local
convergence of Local SFGD for a general class of restricted strongly
convex/smooth loss functions. Local SFGD converges locally to a small
neighborhood of the global optimum at a linear rate with a constant step size,
while it locally converges exactly at a sub-linear rate with diminishing step
sizes. With a proper initialization, local convergence results imply global
convergence. We validate our theoretical findings with numerical simulations
of QST on the Greenberger-Horne-Zeilinger (GHZ) state.

##  I Introduction

A fully-functional fault-tolerant quantum computer faces many technical
hurdles. For instance, using superconducting materials technology, quantum
computers must remain cooled at a very low temperature—almost absolute zero—to
preserve coherence [1]. Moreover, environmental noise from the electronics
controlling the quantum system can disrupt the coherence of its qubits. Thus,
the behavior of current quantum computer implementations needs to be
characterized, verified, and rigorously certified, before their widespread
commercial use [2].

Quantum State Tomography (QST) is the canonical procedure to characterize the
state of a quantum system at various steps of a given computation [3]. In
particular, experimental quantum physicists design quantum circuits that in
theory lead to a specific target pure state; then, they compare the prepared
(input) state and the reconstructed (output) state. To do so, measurements are
taken on independently prepared copies of the state of a quantum system, and
then used to estimate the unknown state by post-processing the data [3].
However, the description complexity of a quantum state grows _exponentially
with the number of qubits_ , leading to challenging data acquisition,
processing, and storage. Therefore, as the number of qubits and quantum gates
increases, so does the need for efficient, robust, and experimentally-
accessible protocols to benchmark quantum information processors.

A quantum state can be represented by a density matrix $\rho$ which is a
complex, positive semi-definite (PSD) matrix with unit trace. The goal of QST
is to design protocols that estimate $\rho$. For an $n$-qubit mixed
state111Mixed states are the most general way to express a quantum state.
$\Psi\in\mathbb{C}^{2^{n}}$, its density matrix can be written as a mixture of
$r$ pure states:
$\rho=\sum_{k}^{r}p_{k}\Psi_{k}\Psi_{k}^{\dagger}\in\mathbb{C}^{2^{n}\times
2^{n}},$ where $(\cdot)^{\dagger}$ denotes the complex conjugate operator.
Here, $p_{k}$ is the probability of finding $\rho$ in the pure state
$\Psi_{k}$. Given these definitions, QST can be formulated as the estimation
of a low-rank density matrix $\rho^{\star}\in\mathbb{C}^{d\times d}$ on an
$n$-qubit Hilbert space with dimension $d=2^{n}$, through the following
$\ell_{2}$-norm optimization problem:

$\displaystyle\min_{\rho\in\mathbb{C}^{d\times d}}$ $\displaystyle
F(\rho):=\tfrac{1}{2m}\|\mathcal{A}(\rho)-y\|_{2}^{2}$ (1) subject to
$\displaystyle\rho\succeq 0,~{}\text{rank}(\rho)\leqslant r,$

where $\mathcal{A}:\mathbb{C}^{2^{n}\times 2^{n}}\rightarrow\mathbb{R}^{m}$ is
the linear sensing map such that $\mathcal{A}(\rho)_{k}=\text{Tr}(A_{k}\rho)$,
for $k=1,\dots,m$. The sensing map used in QST has a particular structure: it
is the Kronecker product of Pauli matrices $A_{k}$, and is closely related to
how quantum computers take measurements in practice [4].

The exponential dependency on the number of qubits implies that $\rho$ has
more than a trillion entries for a $20$-qubit system. Storing this matrix
demands tens of terabytes of memory, which is only available as distributed
memory in sizable clusters. Analogously, a quantum system with more than 30
qubits would demand $100\times$ more memory than the one present in the
world’s fastest supercomputer.

To alleviate these challenges, we study the following distributed optimization
problem to be jointly solved over a set of $M$ machines:

$\displaystyle\min_{X\in\mathbb{R}^{d\times d}}$
$\displaystyle\bigg{\\{}f(X)=\frac{1}{M}\sum_{i=1}^{M}f_{i}(X)\bigg{\\}}$ (2)
subject to $\displaystyle X\succeq 0,{~{}\text{rank}(X)\leqslant r},$

In (2),
$f_{i}(X):=\mathbb{E}_{j\sim\mathcal{D}_{i}}\big{[}f_{i}^{j}(X)\big{]}$, and
$f_{i}^{j}(X)$ is the loss function evaluated at the $j$-th observation of the
locally stored dataset of machine $i$, which follows the distribution
$\mathcal{D}_{i}.$222We assume the homogeneous data case where
$\mathcal{D}_{i}=\mathcal{D}$ for all $i$. The function
$f_{i}:\mathbb{R}^{d\times d}\rightarrow\mathbb{R}$ is a (restricted) strongly
convex/smooth differentiable function, and $X\succeq 0$ is the set of positive
semi-definite matrices with $\text{rank}(X)\leqslant r$.333We provide theory
for the real case; extensions to complex domains can be obtained with complex
conversions and Wirtinger derivatives [5].

To solve (2), we introduce the Local Stochastic Factored Gradient Descent
(Local SFGD) algorithm, and prove its convergence. To the best of our
knowledge, this is the first work that studies Local SGD in the non-convex
factorized objective, and provides convergence in terms of the distance to the
optimal model parameter. Our contributions can be summarized as:

  * •

We introduce a distributed problem setup for QST as an instance of (2).

  * •

We propose Local SFGD, a distributed algorithm that uses matrix factorization
and utilizes local stochastic gradient steps for the minimization of a non-
convex function.

  * •

We provide local convergence guarantees for Local SFGD for restricted strongly
convex/smooth losses, which is of independent interest, and subsumes the QST
problem as a special case.

  * •

We corroborate our theoretical findings with numerical simulations of QST for
the Greenberger-Horne-Zeilinger (GHZ) state.

This paper is organized as follows. Section II reviews the QST protocol, and
sets up the non-convex distributed objective. Section III introduces the Local
SFGD algorithm, followed by Section IV where we provide the main theoretical
results along with the proofs. Lastly, in Section V, we use Local SFGD for the
reconstruction of the GHZ state.

##  II Preliminaries

Classically, the sample complexity $m$ for reconstructing
$\rho^{\star}\in\mathbb{C}^{d\times d}$ is $O(d^{2}),$ where $d$ itself grows
exponentially with $n$. To address such large sample complexity requirements,
we use low-rankness as prior, as many lab-constructed density matrices have
low-rank structure, including the maximally-entangled Greenberger-Horne-
Zeilinger (GHZ) state [6]. While the low-rank constraint is non-convex, it
provides a significant reduction in the sample complexity. Under appropriate
assumptions, a rank-$r$ density matrix can be reconstructed with $m=O(r\cdot
d\cdot\text{poly}\log(d))$ measurements, instead of $m=O(d^{2})$ [7].

We propose to solve a factorized version of (1) to efficiently handle its low-
rank constraint, following [8, 9]:

$\min_{U\in\mathbb{C}^{d\times
r}}~{}G(U):=F(UU^{\dagger})=\tfrac{1}{2m}\|\mathcal{A}(UU^{\dagger})-y\|_{2}^{2}.$
(3)

In (3), we parametrize the low-rank density matrix $\rho$ by its factor
$U\in\mathbb{C}^{d\times r}$. By rewriting $\rho=UU^{\dagger}$, both the PSD
and the low-rank constraints are automatically satisfied, leading to the
unconstrained non-convex formulation in (3). Moreover, working in the factored
space improves time and space complexities [8, 9, 10]. However, even with the
reduced sample complexity $m=O(r\cdot d\cdot\text{poly}\log(d))$, linear
dependency on $d=2^{n}$ makes computation infeasible, e.g., for $n=20$ and
rank $r=100$, the reduced sample complexity still reaches $2.02\times
10^{10}$.

To handle this explosion of data, we consider the setting where the
measurements $y\in\mathbb{R}^{m}$ and the sensing matrices
$\mathcal{A}:\mathbb{C}^{d\times d}\rightarrow\mathbb{R}^{m}$ from a central
quantum computer _are locally stored across $M$ different classical machines_.
These classical machines perform some local operations based on their local
data, and communicate back and forth with the central quantum server to
reconstruct a density matrix.

The _distributed_ QST problem can be written as:

$\begin{split}\min_{U\in\mathbb{C}^{d\times
r}}\Big{\\{}g(U)&=\frac{1}{M}\sum_{i=1}^{M}g_{i}(U)\Big{\\}},\\\
\text{where}\quad
g_{i}(U)&:=\mathbb{E}_{j\sim\mathcal{D}_{i}}\|\mathcal{A}_{i}^{j}(UU^{\dagger})-y_{i}^{j}\|_{2}^{2},\end{split}$
(4)

with $j$ being a random variable that follows a distribution $\mathcal{D}_{i}$
for machine $i$. In the next section, we introduce our approach to solve (4),
which can be more generally applied to (2).

##  III Algorithms

We now introduce the _Local Stochastic Factored Gradient Descent_ (Local SFGD)
algorithm. We review the _Factored Gradient Descent_ (FGD) algorithm [8, 9,
10] and its stochastic variant [11], on which the Local SFGD is based.

$\diamond$ Factored Gradient Descent (FGD). A common approach to solve the
factorized non-convex objective in (3) in _centralized_ settings is to use
gradient descent on the factor $U$:

$\displaystyle U_{t+1}$ $\displaystyle=U_{t}-\eta_{t}\nabla
G(U_{t})=U_{t}-\eta_{t}\nabla F(U_{t}U_{t}^{\dagger})\cdot U_{t}$
$\displaystyle={U_{t}-\frac{\eta_{t}}{m}\Big{(}\sum_{k=1}^{m}\left\\{\text{Tr}(A_{k}U_{t}U_{t}^{\dagger})-y_{k}\right\\}A_{k}\Big{)}\cdot
U_{t},}$ (5)

where $\eta_{t}>0$ is the step size. From (5), we can see that a pass over
full data is required to compute the gradient on every iteration. This can be
computationally challenging or even infeasible when $m$ is large, which is
almost always the case for QST, even for moderate number of qubits $n$.

$\diamond$ Stochastic Factored Gradient Descent (SFGD). A simple and effective
way to mitigate this burden is to use _Stochastic Factored Gradient Descent_
(SFGD), which replaces the true gradient $\nabla G$ with an unbiased estimator
$H$. For instance, one can use the following SFGD update:

$\displaystyle U_{t+1}$ $\displaystyle=U_{t}-\eta_{t}\cdot H(U_{t})$
$\displaystyle={U_{t}-\frac{\eta_{t}}{b}\Big{(}\sum_{k=1}^{b}\left\\{\text{Tr}(A_{k}U_{t}U_{t}^{\dagger})-y_{k}\right\\}A_{k}\Big{)}\cdot
U_{t}},$ (6)

which simply uses $b$ measurements instead of $m\gg b$ to approximate $\nabla
G$, where the hyperparameter $b$ is the batch size. In [11], the convergence
of SFGD was shown for (restricted) strongly convex/smooth functions.

From (6), one can see that SFGD is amenable to parallelization, simply by
replacing $H(U_{t})$ with an average of stochastic gradients that are computed
independently from the local machines. This scheme is often the state-of-the-
art in distributed learning problems [12, 13]. However, it exhibits a major
drawback: on every iteration, each machine has to send the local (stochastic)
gradient to the server, and receive back the aggregated model parameter. Such
communication is much more expensive than—typically about 3 orders of
magnitude—the local computations that each machine has to perform [14].

$\diamond$ Our approach: Local SFGD. There are two main approaches to resolve
the aforementioned communication overhead. One is to reduce the number of
transmitted bits via gradient compression schemes, such as quantization [15]
or sparsification [16]. The other is to increase the amount of local
iterations performed on each machine, in order to reduce the total
communication rounds. The latter approach is called Local SGD, and was shown
to outperform (parallel) SGD in some settings [17, 18, 19, 20].

In this work, we introduce the Local SFGD to estimate the low-rank factor of a
density matrix over a set of local machines. Although our main application is
to solve the distributed QST objective in (4), Local SFGD is more generally
applicable to the distributed objective in (2); see Section IV for details.
Local SFGD is summarized in Algorithm 1. While there are non-convex results on
Local SGD [21, 22], they consider a different problem setting, and only
provide convergence in terms of the norm of the gradient. To the best of our
knowledge, this is the first work that studies Local SGD in the non-convex
factorized objective, and provides convergence in terms of the distance to the
optimal model parameters.

Local SFGD produces $M$ sequences in parallel, where $M$ is the number of
machines. If a synchronization step happens at time $t$, i.e., $t=t_{p}$ for
some $p\in\mathbb{N},$ then the local parameters at each machine $U_{t}^{i}$
are sent to the central server, and their average is computed (line 6).
Otherwise, each machine performs (possibly many iterations of) SFGD without
communicating with the central server (line 8). An important metric to
consider for Local SFGD is the maximum time interval between two
synchronization time steps: $\max_{p}|t_{p}-t_{p+1}|,$ which we assume is
bounded by $h\geqslant 1$; see also Theorems 2 and 4. If communication happens
on every iteration, i.e., $h=1,$ then Algorithm 1 reduces to the (parallel)
SFGD in (6).

Algorithm 1 Local SFGD

1:Set number of iterations $T>0$, synchronization time steps
$t_{1},t_{2},\dots$, and initialize $U_{0}^{i}=U_{0}$ as below:
$U_{0}^{i}=\texttt{SVD}\Big{(}-\sum_{i=1}^{M}\tfrac{m_{i}}{m}\nabla
f_{i}(0)\Big{)}\quad\forall i\in[M],\vspace{-2mm}$ (7) where SVD denotes the
singular value decomposition.

2:for each round $t=0,\ldots T$ do

3: for in parallel for $i\in[M]$ do

4: Sample $j_{t}$ uniformly at random from $[m_{i}]$.

5: if $t=t_{p}$ for some $p\in\mathbb{N}$ then

6: $U_{t+1}^{i}=\frac{1}{M}\sum_{i=1}^{M}\big{(}U_{t}^{i}-\eta_{t}\nabla
g_{i}^{j_{t}}(U_{t}^{i})\big{)}$

7: else

8: $U_{t+1}^{i}=U_{t}^{i}-\eta_{t}\nabla g_{i}^{j_{t}}(U_{t}^{i})$

9: end if

10: end for

11:end for

12:return $\hat{U}_{T+1}:=\frac{1}{M}\sum_{i=1}^{M}U_{T+1}^{i}.$

##  IV The Convergence of Local SFGD

We provide local convergence guarantees of Local SFGD in Algorithm 1 for
_restricted $\mu$-strongly convex/$L$-smooth objectives_. Similarly to (3), as
we parametrize $X=UU^{\top},$ Problem (2) becomes non-convex, where:

$\min_{U\in\mathbb{R}^{d\times
r}}\Big{\\{}g(U)=\frac{1}{M}\sum_{i=1}^{M}g_{i}(U)\Big{\\}},$ (8)

which now is unconstrained, as both the PSD and the low-rank constraints are
automatically satisfied.

We assume $f_{i}$ is a symmetric function: $f_{i}(X)=f_{i}(X^{\top}).$ Then,
the gradient of $g_{i}(U)=f_{i}(UU^{\top})$ simplifies to:444Without loss of
generality, we absorb 2 into $\eta_{t}$ to use $\nabla g_{i}(U)$ and $\nabla
f_{i}(UU^{\top})U$ interchangeably.

$\displaystyle\nabla g_{i}(U)=\left(\nabla f_{i}(UU^{\top})+\nabla
f_{i}(UU^{\top})^{\top}\right)U=2\nabla f_{i}(UU^{\top})U.$

We now state key assumptions used in our main results.

######  Assumption 1.

The function $f_{i}$ is $\mu$-restricted strongly convex and $L$-restricted
smooth. That is, $\forall X,Y\succeq 0$ and $\forall i\in[M]$, it holds that

$\displaystyle f_{i}(Y)\geqslant f_{i}(X)$ $\displaystyle+\langle\nabla
f_{i}(X),Y{-}X\rangle+\tfrac{\mu}{2}\|X-Y\|_{F}^{2},$ (I-a)
$\displaystyle\text{and}\quad\|\nabla f_{i}(X)$ $\displaystyle-\nabla
f_{i}(Y)\|_{F}\leqslant L\|X-Y\|_{F}.$ (I-b)

######  Assumption 2.

The stochastic gradient $\nabla g_{i}^{j}$ is unbiased, has a bounded
variance, and is bounded in expectation, $\forall i\in[M].$ That is,

$\displaystyle\mathbb{E}_{j}\big{[}\nabla g_{i}^{j}(U)\big{]}=\nabla
g_{i}(U),$ (II-a) $\displaystyle\mathbb{E}_{j}\big{[}\|\nabla
g_{i}^{j}(U)-\nabla
g_{i}(U)\|_{F}^{2}\big{]}\leqslant\sigma^{2},\quad\text{and}$ (II-b)
$\displaystyle{\mathbb{E}_{j}\big{[}\|\nabla
g_{i}^{j}(U)\|_{F}^{2}\big{]}\leqslant G^{2}},$ (II-c)

where $j$ follows a uniform distribution.

Assumptions (I-a) and (I-b) respectively state that $\mu$-strong convexity and
$L$-smoothness hold when we restrict the space of $d\times d$ matrices to the
set of PSD matrices. Such assumptions have become standard in optimization
analysis, and are significantly weaker than assuming global strong convexity.
Importantly, note that we only assume $f_{i}(X)$ to have such structures—the
transformed function $g_{i}(U)$ in (8) typically does not satisfy restricted
strong convexity/smoothness [23].

Assumptions (II-a) and (II-b) respectively imply that the stochastic gradient
is unbiased and has a bounded variance, and both are standard assumptions in
stochastic optimization [24]. Assumption (II-c) states that the stochastic
gradient has a bounded norm, in expectation. This assumption may seem strong
when the objective is (unconstrained) strongly convex [25]; however, note that
Assumption (I-a) is restricted to PSD matrices, and the original Problem (2)
is constrained.

Apart from (8) being non-convex, another difficulty that arises by the
parametrization $X=UU^{\top}$ is that the solution can become non-
unique.555Consider reconstructing $X^{\star}=\begin{bmatrix}1&1\\\
1&1\end{bmatrix}$. It can be seen that both $U^{\star}=[1~{}1]^{\top}$ and
$\tilde{U}^{\star}=-[1~{}1]^{\top}$ satisfy
$U^{\star}U^{\star\top}=\tilde{U}^{\star}\tilde{U}^{\star\top}=X^{\star}.$ We
can remove this ambiguity by defining the following rotation invariant
distance metric.

######  Definition 1 (Eq. ($3.1$) in [10]).

For any $U,V\in\mathbb{R}^{d\times r},$ let
$D(U,V):=\min_{R\in\mathcal{O}}\|U-VR\|_{F},$ where
$\mathcal{O}\subseteq\mathbb{R}^{r\times r}$ is the set of orthonormal
matrices such that $R^{\top}R=\mathbb{I}_{r\times r}.$

######  Remark 1.

Definition 1 regards all $U,V\in\mathbb{R}^{d\times r}$ to be in the same
distance such that $D(U,VR)=D(UR,V)=D(U,V).$ Hence, it defines the equivalence
classes $\\{UR:R^{\top}R=\mathbb{I}_{r\times r}\\}$ and
$\\{VR:R^{\top}R=\mathbb{I}_{r\times r}\\}$ [23].

A crucial component for our convergence analysis is the following lemma, which
replaces the role of (strong) convexity in classical convergence analysis of
gradient descent:

######  Lemma 1 (Lemma 14 in [26]).

Let Assumption 1 hold. Assume that
$D^{2}(U_{0}^{i},U^{\star})\leqslant\frac{\sigma_{r}(X^{\star})}{100\cdot\kappa\cdot\sigma_{1}(X^{\star})}$,
where $\sigma_{k}(X^{\star})$ is the $k$-th singular value of
$X^{\star}$,666Without loss of generality, singular values are sorted in
descending order. and $\kappa:=\frac{L}{\mu}.$ Then, the following inequality
holds:

$\begin{split}&\big{\langle}U_{t}^{i}-U^{\star}R^{\star},\nabla
g_{i}(U_{t}^{i})\big{\rangle}\\\
&\quad\quad\geqslant\tfrac{2\eta_{t}}{3}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+\tfrac{3\mu}{20}\sigma_{r}(X^{\star})\cdot
D^{2}(U_{t}^{i},U^{\star}).\end{split}$ (11)

######  Remark 2.

The initialization scheme (7) in Algorithm 1 is modified from [26, Theorem 11]
to distributed version, and satisfies the initialization condition in Lemma 1
for small enough $\kappa$; for the QST problem, the Pauli sensing matrices
$A_{k}$ satisfy the Restriced Isometry Property (RIP) [4, 27], implying
$\kappa\approx\tfrac{1+\delta}{1-\delta}$, where $\delta\in(0,1)$ is the RIP
constant. Hence, by using right prior information (e.g., low-rankness), we can
apply the compressed sensing results, implying that $X^{\star}$ is unique and
can be recovered exactly [28].

We are now ready to present the main theoretical results. We first show in
Theorem 2 that the Local SFGD converges locally at a linear rate to a small
neighborhood of the global optimum with a constant step size. Then, in Theorem
4, we show the exact local convergence by using an appropriately diminishing
step size, at the expense of reducing the convergence rate to a sub-linear
rate.

######  Theorem 2 (Local linear convergence with constant step size).

Let Assumptions 1, 2, and the initialization condition of Lemma 1 hold.
Moreover, let $\eta_{t}=\eta<\tfrac{1}{\alpha}$ for $t\in[0:T]$ and
$\max_{p}|t_{p}-t_{p+1}|\leqslant h$. Then, the output of Algorithm 1 has the
following property:

$\begin{split}\mathbb{E}\big{[}D^{2}(\hat{U}_{T+1},U^{\star})\big{]}&\leqslant\left(1-\eta\alpha\right)^{T+1}D^{2}(\hat{U}_{0},U^{\star})\\\
&\quad+\eta\left(\tfrac{(h-1)^{2}G^{2}}{\alpha}+\tfrac{\sigma^{2}}{M\alpha}\right),\end{split}$
(12)

where $X^{\star}$ is the optimum of $f$ over the set of PSD matrices such that
$\text{rank}(X^{\star})=r$, $U^{\star}$ is such that
$X^{\star}=U^{\star}U^{\star\top}$, and
$\alpha=\tfrac{3\mu}{10}\sigma_{r}(X^{\star})$ is a global constant.

######  Remark 3.

In (12), the expectation is with respect to the previous iterates,
$\\{\hat{U}_{t}\\}_{t=0}^{T}$. We make a few remarks about Theorem 2. First,
notice the last variance term $\tfrac{\sigma^{2}}{M\alpha}$, which disappears
in the noiseless case, is reduced by the number of machines $M$. Second, we
assume a single-batch is used in the proof; by using batch size $b>1$, this
term can be further divided by $b$. Lastly, by plugging in $h=1$ (i.e.,
synchronization happens on every iteration), the first variance term
disappears, exhibiting similar local linear convergence to SFGD [11].

###### Proof.

We start with the following auxiliary lemma.

######  Lemma 3.

Let Assumptions 1 and (II-c) hold. Then, the output of Algorithm 1 with
$\max_{p}|t_{p+1}-t_{p}|\leqslant h$ satisfies:

$\displaystyle{\frac{1}{M}\sum_{i=1}^{M}\mathbb{E}\big{[}\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}\big{]}\leqslant(h-1)^{2}\eta_{t}^{2}G^{2}.}$
(13)

The proof is almost identical to [17, Lemma 3.3], and hence is omitted.
Throughout the proof, we use the notations:

$\displaystyle
U_{t+1}^{i}=U_{t}^{i}-\eta_{t}g_{t}^{i}\quad\text{and}\quad\hat{U}_{t+1}=\hat{U}_{t}-\eta_{t}g_{t},$

where $\hat{U_{t}}=\frac{1}{M}\sum_{i=1}^{M}U_{t}^{i}$, i.e., the average
across different machines at time $t$. We denote the stochastic gradient of
machine $i$ at time $t$ with $g_{t}^{i}:=\nabla
f_{i}^{j_{t}}(U_{t}^{i}U_{t}^{i\top})U_{t}^{i}=\nabla
g_{i}^{j_{t}}(U_{t}^{i})$, and the average of stochastic gradients across
machines with $g_{t}=\frac{1}{M}\sum_{i=1}^{M}g_{t}^{i}$. Finally, we denote
$\mathbb{E}\left[g_{t}\right]=\bar{g_{t}}.$

We first decompose the distance of $D^{2}(\hat{U}_{t+1},U^{\star})$:

$\displaystyle D^{2}(\hat{U}_{t+1},U^{\star})$
$\displaystyle=\min_{R\in\mathcal{O}}\|\hat{U}_{t+1}{-}U^{\star}R\|_{F}^{2}\leqslant\|\hat{U}_{t+1}{-}U^{\star}R^{\star}\|_{F}^{2}$
$\displaystyle=\|\hat{U}_{t}-U^{\star}R^{\star}-\eta_{t}\bar{g}_{t}\|_{F}^{2}+\eta_{t}^{2}\|\bar{g}_{t}-g_{t}\|_{F}^{2}$
$\displaystyle\quad+2\eta_{t}\langle\hat{U}_{t}-U^{\star}R^{\star}-\eta_{t}\bar{g}_{t},\bar{g}_{t}-g_{t}\rangle.$
(14)

The first term in (14) can be further decomposed to:

$\displaystyle\|\hat{U}_{t}-U^{\star}R^{\star}\|_{F}^{2}+\eta_{t}^{2}\|\bar{g}_{t}\|_{F}^{2}-2\eta_{t}\langle\hat{U}_{t}-U^{\star}R^{\star},\bar{g}_{t}\rangle.$
(15)

We bound the second and the third terms of (15) separately. For the second
term, by Jensen’s inequality, we have:

$\displaystyle\|\bar{g}_{t}\|_{F}^{2}$
$\displaystyle=\big{\|}\tfrac{1}{M}\sum_{i=1}^{M}\nabla
g_{i}(U_{t}^{i})\big{\|}_{F}^{2}\leqslant\tfrac{1}{M}\sum_{i=1}^{M}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}.$ (16)

For the third term, we decompose further to have:

$\displaystyle\big{\langle}\hat{U}_{t}{-}U^{\star}R^{\star},\tfrac{1}{M}\sum_{i=1}^{M}\nabla
g_{i}(U_{t}^{i})\big{\rangle}$
$\displaystyle=\tfrac{1}{M}\sum_{i=1}^{M}\big{\langle}\hat{U}_{t}{-}U_{t}^{i},\nabla
g_{i}(U_{t}^{i})\big{\rangle}$
$\displaystyle+\tfrac{1}{M}\sum_{i=1}^{M}\big{\langle}U_{t}^{i}-U^{\star}R^{\star},\nabla
g_{i}(U_{t}^{i})\big{\rangle}.$ (17)

We again bound the two terms in (17) separately. Using $\langle
A,B\rangle\geqslant-\tfrac{\delta}{2}\|A\|_{F}^{2}-\tfrac{1}{2\delta}\|B\|_{F}^{2},$
The first term admits:

$\displaystyle\tfrac{1}{M}\sum_{i=1}^{M}\big{\langle}\hat{U}_{t}-U_{t}^{i},\nabla
g_{i}(U_{t}^{i})\big{\rangle}$
$\displaystyle{\stackrel{{\scriptstyle}}{{\geqslant}}\tfrac{1}{M}\sum_{i=1}^{M}\big{(}-\tfrac{\delta}{2}\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}-\tfrac{1}{2\delta}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}\big{)}}$

By Lemma 1, the second term in (17) admits:

$\displaystyle\tfrac{1}{M}\sum_{i=1}^{M}\big{\langle}U_{t}^{i}-U^{\star}R^{\star},\nabla
g_{i}(U_{t}^{i})\big{\rangle}$
$\displaystyle\geqslant\tfrac{1}{M}\sum_{i=1}^{M}\tfrac{2\eta_{t}}{3}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+\tfrac{1}{M}\sum_{i=1}^{M}\tfrac{3\mu\cdot\sigma_{r}(X^{\star})}{20}D^{2}(U_{t}^{i},U^{\star})$
$\displaystyle\geqslant\tfrac{1}{M}\sum_{i=1}^{M}\tfrac{2\eta_{t}}{3}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+\tfrac{3\mu\cdot\sigma_{r}(X^{\star})}{20}D^{2}(\hat{U}_{t},U^{\star}),$

where we used convexity of $D^{2}(\cdot,\cdot)$ in the second inequality.

Combining above two bounds into (17), we have

$\displaystyle\langle\bar{g}_{t},\hat{U}_{t}-U^{\star}R^{\star}\rangle$
$\displaystyle\geqslant\tfrac{1}{M}\sum_{i=1}^{M}\big{(}-\tfrac{\delta}{2}\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}-\tfrac{1}{2\delta}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}\big{)}$
$\displaystyle+\tfrac{1}{M}\sum_{i=1}^{M}\tfrac{2\eta_{t}}{3}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+\tfrac{3\mu\cdot\sigma_{r}(X^{\star})}{20}D^{2}(\hat{U}_{t},U^{\star}).$
(18)

Substituting (16) and (18) into (15), we have

$\displaystyle\|\hat{U}_{t}-U^{\star}R^{\star}-\eta_{t}\bar{g}_{t}\|_{F}^{2}$
$\displaystyle\quad\leqslant\left(1-\eta_{t}\cdot\tfrac{3\mu}{10}\sigma_{r}(X^{\star})\right)\|\hat{U}_{t}-U^{\star}R^{\star}\|_{F}^{2}$
$\displaystyle\quad\quad+\tfrac{1}{M}\sum_{i=1}^{M}\big{[}\big{(}\tfrac{\eta_{t}}{\delta}-\tfrac{\eta_{t}^{2}}{3}\big{)}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+\eta_{t}\delta\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}\big{]}$
$\displaystyle\hskip
2.84526pt{\stackrel{{\scriptstyle\delta=4/\eta_{t}}}{{=}}\left(1-\eta_{t}\alpha\right)\|\hat{U}_{t}{-}U^{\star}R^{\star}\|_{F}^{2}}$
$\displaystyle\quad\quad{+\tfrac{1}{M}\sum_{i=1}^{M}\big{[}\eta_{t}\big{(}\tfrac{\eta_{t}}{4}-\tfrac{\eta_{t}}{3}\big{)}\|\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}+4\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}\big{]}}$
$\displaystyle\quad{\stackrel{{\scriptstyle}}{{\leqslant}}\left(1-\eta_{t}\alpha\right)\|\hat{U}_{t}{-}U^{\star}R^{\star}\|_{F}^{2}+\tfrac{4}{M}\sum_{i=1}^{M}\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2},}$
(19)

where in the equality we defined
$\alpha:=\tfrac{3\mu}{10}\sigma_{r}(X^{\star}),$ and in the last inequality we
used that $\tfrac{\eta_{t}}{4}-\tfrac{\eta_{t}}{3}<0$.

Substituting (19) into (14) and taking expectations conditional on previous
iterates, and using $\mathbb{E}\left[g_{t}\right]=\bar{g}_{t}$, we get

$\displaystyle\mathbb{E}[D^{2}(\hat{U}_{t+1},U^{\star})]=\|\hat{U}_{t}{-}U^{\star}R^{\star}{-}\eta_{t}\bar{g}_{t}\|_{F}^{2}+\eta_{t}^{2}\mathbb{E}[\|\bar{g}_{t}{-}g_{t}\|_{F}^{2}]$
$\displaystyle{\stackrel{{\scriptstyle\eqref{eq:main-rec-first-term-
bound}}}{{\leqslant}}\left(1-\eta_{t}\alpha\right)\|\hat{U}_{t}-U^{\star}R^{\star}\|_{F}^{2}}$
$\displaystyle{\quad\quad+\tfrac{4}{M}\sum_{i=1}^{M}\mathbb{E}\big{[}\|\hat{U}_{t}-U_{t}^{i}\|_{F}^{2}\big{]}+\eta_{t}^{2}\mathbb{E}[\|\bar{g}_{t}-g_{t}\|_{F}^{2}]}$
$\displaystyle{\stackrel{{\scriptstyle\eqref{eq:avg-local-
diff}}}{{\leqslant}}\left(1-\eta_{t}\alpha\right)\|\hat{U}_{t}-U^{\star}R^{\star}\|_{F}^{2}}$
$\displaystyle{\quad\quad+4\eta_{t}^{2}(h-1)^{2}G^{2}+\eta_{t}^{2}\mathbb{E}[\|\bar{g}_{t}-g_{t}\|_{F}^{2}].}$

where the last inequality is by Lemma 3. We further have:

$\displaystyle\mathbb{E}[\|\bar{g}_{t}-g_{t}\|_{F}^{2}]$
$\displaystyle=\mathbb{E}\Big{[}\Big{\|}\tfrac{1}{M}\sum_{i=1}^{M}\big{(}\nabla
g_{i}^{j_{t}}(U_{t}^{i})-\nabla
g_{i}(U_{t}^{i})\big{)}\Big{\|}_{F}^{2}\Big{]}$
$\displaystyle\stackrel{{\scriptstyle}}{{\leqslant}}\tfrac{1}{M^{2}}\sum_{i=1}^{M}\mathbb{E}[\|\nabla
g_{i}^{j_{t}}(U_{t}^{i})-\nabla
g_{i}(U_{t}^{i})\|_{F}^{2}]\stackrel{{\scriptstyle\eqref{A2:eq:stoc-grad-var-
bound}}}{{\leqslant}}\tfrac{\sigma^{2}}{M},$

where we used
$\text{Var}(\sum_{m=1}^{M}X_{m})=\sum_{m=1}^{M}\text{Var}(X_{m})$ for
independent random variables.

We now arrive at the iteration invariant bound:

$\displaystyle\mathbb{E}[D^{2}(\hat{U}_{t+1},U^{\star})]$
$\displaystyle\leqslant\left(1-\eta_{t}\alpha\right)D^{2}(\hat{U}_{t},U^{\star})+\eta_{t}^{2}\big{(}4(h-1)^{2}G^{2}+\tfrac{\sigma^{2}}{M}\big{)}.$
(20)

Lastly, we unfold (20) for $T$ iterations, and using
$\sum_{t=0}^{T}(1-\eta_{t}\alpha)^{t}\leqslant\sum_{t=0}^{\infty}(1-\eta_{t}\alpha)^{t}=\frac{1}{\eta_{t}\alpha},$
we obtain

$\displaystyle\mathbb{E}[D^{2}(\hat{U}_{T+1},U^{\star})]$
$\displaystyle\leqslant\left(1-\eta_{t}\alpha\right)^{T+1}D^{2}(\hat{U}_{0},U^{\star})$
$\displaystyle\quad+\eta_{t}\left(\tfrac{4(h-1)^{2}G^{2}}{\alpha}+\tfrac{\sigma^{2}}{M\alpha}\right),$

which completes the proof. ∎

######  Theorem 4 (Local sub-linear convergence with diminishing step size).

Let Assumptions 1, 2, and the initialization condition of Lemma 1 hold.
Moreover, let $\eta_{t}=\frac{2}{\alpha(t+2)}$ for $t\in[0:T]$ and
$\max_{p}|t_{p}-t_{p+1}|\leqslant h$. Then, the output of Algorithm 1 has the
following property:

$\begin{split}\mathbb{E}\big{[}D^{2}(\hat{U}_{T+1},U^{\star})\big{]}&\leqslant\tfrac{4C}{\alpha(T+3)},\end{split}$
(21)

where $X^{\star}$ is the optimum such that $\text{rank}(X^{\star})=r$,
$U^{\star}$ is such that $X^{\star}=U^{\star}U^{\star\top}$, and
$\alpha=\tfrac{3\mu}{10}\sigma_{r}(X^{\star})$ and
$C=4(h-1)^{2}G^{2}+\tfrac{\sigma^{2}}{M}$ are global constants.

###### Proof.

We claim the following, and prove by induction:

$\displaystyle
D^{2}(\hat{U}_{t},U^{\star})\leqslant\tfrac{4C}{\alpha^{2}(t+2)},\quad\text{with}\quad\eta_{t}=\tfrac{2}{\alpha(t+2)}.$
(22)

We start from the iteration invariant bound in (20):

$\displaystyle\mathbb{E}\big{[}D^{2}(\hat{U}_{t+1},U^{\star})\big{]}$
$\displaystyle\leqslant\left(1-\eta_{t}\alpha\right)D^{2}(\hat{U}_{t},U^{\star})+\eta_{t}^{2}\cdot
C.$

For the base case $t=0$, we have

$\displaystyle\mathbb{E}\big{[}D^{2}(\hat{U}_{1},U^{\star})\big{]}\leqslant(1-\eta_{0}\alpha)D^{2}(\hat{U}_{0},U^{\star})+\eta_{0}^{2}\cdot
C$
$\displaystyle=\left(1-\tfrac{1}{\alpha}\cdot\alpha\right)D^{2}(\hat{U}_{0},U^{\star})+\tfrac{C}{\alpha^{2}}=\tfrac{C}{\alpha^{2}}\leqslant\tfrac{4C}{3\alpha^{2}}.$

Now, we proceed to the inductive step. Assuming (22) holds for time step $t$,
we want to prove the same holds for the time step $t+1$. Starting from (20)
again, we have:

$\displaystyle\mathbb{E}\big{[}D^{2}(\hat{U}_{t+1},U^{\star})\big{]}\leqslant\left(1-\eta_{t}\alpha\right)D^{2}(\hat{U}_{t},U^{\star})+\eta_{t}^{2}\cdot
C$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:induction-
hypothesis}}}{{\leqslant}}\left(1-\tfrac{2}{t+2}\right)\cdot\tfrac{4C}{\alpha^{2}(t+2)}+\tfrac{4C}{\alpha^{2}(t+2)^{2}}=4C\cdot\tfrac{t+1}{\alpha^{2}(t+2)}\cdot\tfrac{1}{t+2}$
$\displaystyle\leqslant
4C\cdot\tfrac{t+2}{\alpha^{2}(t+3)}\cdot\tfrac{1}{t+2}=\tfrac{4C}{\alpha^{2}(t+3)},$

where in the last inequality we used the fact that
$\frac{t+1}{t+2}\leqslant\frac{t+2}{t+3}.$ This completes the proof. ∎

##  V Numerical Results

We use Local SFGD to reconstruct the Greenberger-Horne-Zeilinger (GHZ) state,
using simulated measurement data from Qiskit. GHZ state is known as _maximally
entangled_ quantum state [6], meaning it exhibits the maximal inter-particle
correlation, which does not exist in classical mechanics. We are interested
in: $(i)$ how the number of local steps affect the accuracy defined as
$\varepsilon=\|\hat{U}_{t}\hat{U}_{t}^{\top}-\rho^{\star}_{\text{ghz}}\|_{F}^{2}$,
where $\rho^{\star}_{\text{ghz}}=U^{\star}U^{\star\dagger}$ is the true
density matrix for the GHZ state; and $(ii)$ the scalability of the
distributed setup for various number of classical machines $M$.

In Fig. 1 (Top), we first fix the number of machines $M=10$ and the number of
total synchronization steps to be $100$, and vary the number of local
iterations between two synchronization steps, i.e.,
$h\in\\{1,10,25,50,100,200\\}.$ We use constant step size $\eta=1$ for all
$h$. Increasing $h$, i.e., each distributed machine performing more local
iterations, leads to faster convergence in terms of the synchronization steps.
Notably, the speed up gets marginal: e.g., there is not much difference
between $h=100$ and $h=200$, indicating there is an “optimal” $h$ that leads
to the biggest reduction in the number of synchronization steps. Further, one
can notice that higher $h$ leads to slightly worse final accuracy—this is
consistent with (12) in Theorem 2, where the first variance term that depends
on $G^{2}$ disappears with $h=1$. Finally, note that $\varepsilon$ does not
decrease below certain level due to the inherent finite sampling error of
quantum measurements [29].

In Fig. 1 (Bottom), we plot the number of synchronization steps to reach
$\varepsilon\leqslant 0.05,$ while fixing $h=20$. We vary the number of
workers $M\in\\{5,10,15,20\\}$, where each machine gets $200$ measurements.
There is a significant speed up from $M=5$ to $M=15$, while for $M=20,$ it
took one more syncrhonization step compared to $M=15,$ which is likely due to
the stochasticity of SFGD within each machine.

Figure 1: Top: Convergence speed as a function of number of synchronization
steps $t_{p}$ for various number of local iterations. Bottom: number of
synchronization steps to reach $\varepsilon\leqslant 0.05$ as a function of
number of workers $M$. The batch size $b=50$ is used for all cases.

##  VI Conclusion and Future Work

In this work, we introduced a distributed problem set up for QST as an
instance of a general distributed optimization problem with PSD/low-rank
constraints. We proposed the Local SFGD, a distributed non-convex algorithm
that utilizes local steps at each distributed worker to estimate the low-rank
factor of a density matrix. We proved the local convergence of Local SFGD for
restricted strongly convex/smooth objectives, which can be of independent
interest. For future work, extension to the heterogeneous data case as well as
the decentralized case with various topologies can be investigated.

## References

  * [1] Jee Woo Park, Zoe Z Yan, Huanqian Loh, Sebastian A Will, and Martin W Zwierlein,  “Second-scale nuclear spin coherence time of ultracold 23Na40K molecules,”  Science, vol. 357, pp. 372–375, 2017.
  * [2] John Preskill,  “Quantum computing in the NISQ era and beyond,”  Quantum, vol. 2, pp. 79, 2018.
  * [3] Joseph B Altepeter, Evan R Jeffrey, and Paul G Kwiat,  “Photonic state tomography,”  Advances in Atomic, Molecular, and Optical Physics, vol. 52, pp. 105–159, 2005.
  * [4] Yi-Kai Liu,  “Universal low-rank matrix recovery from Pauli measurements,”  Advances in NeurIPS, vol. 24, 2011.
  * [5] Robert Clifford Gunning and Hugo Rossi,  Analytic functions of several complex variables, vol. 368,  American Mathematical Soc., 2009.
  * [6] Yajuan Zhao, Rui Zhang, Wenlan Chen, Xiang-Bin Wang, and Jiazhong Hu,  “Creation of Greenberger-Horne-Zeilinger states with thousands of atoms by entanglement amplification,”  npj Quantum Information, vol. 7, no. 1, pp. 1–6, 2021.
  * [7] David Gross, Yi-Kai Liu, Steven T Flammia, Stephen Becker, and Jens Eisert,  “Quantum state tomography via compressed sensing,”  Physical review letters, vol. 105, no. 15, pp. 150401, 2010.
  * [8] Anastasios Kyrillidis, Amir Kalev, Dohyung Park, Srinadh Bhojanapalli, Constantine Caramanis, and Sujay Sanghavi,  “Provable compressed sensing quantum state tomography via non-convex methods,”  npj Quantum Information, vol. 4, no. 1, pp. 1–7, 2018.
  * [9] Junhyung Lyle Kim, George Kollias, Amir Kalev, Ken X Wei, and Anastasios Kyrillidis,  “Fast quantum state reconstruction via accelerated non-convex programming,”  arXiv:2104.07006, 2021.
  * [10] Stephen Tu, Ross Boczar, Max Simchowitz, Mahdi Soltanolkotabi, and Ben Recht,  “Low-rank solutions of linear matrix equations via Procrustes flow,”  in ICML, 2016, pp. 964–973.
  * [11] Jinshan Zeng, Ke Ma, and Yuan Yao,  “On Global Linear Convergence in Stochastic Nonconvex Optimization for Semidefinite Programming,”  IEEE Transactions on Signal Processing, pp. 4261–4275, 2019.
  * [12] Martín Abadi et al.,  “Tensorflow: Large-scale machine learning on heterogeneous distributed systems,”  arXiv:1603.04467, 2016.
  * [13] Priya Goyal, Piotr Dollár, Ross Girshick, Pieter Noordhuis, Lukasz Wesolowski, Aapo Kyrola, Andrew Tulloch, Yangqing Jia, and Kaiming He,  “Accurate, large minibatch SGD: Training ImageNet in 1 hour,”  arXiv preprint arXiv:1706.02677, 2017.
  * [14] Guanghui Lan, Soomin Lee, and Yi Zhou,  “Communication-efficient algorithms for decentralized and stochastic optimization,”  Mathematical Programming, vol. 180, no. 1, pp. 237–284, 2020.
  * [15] Dan Alistarh, Demjan Grubic, Jerry Li, Ryota Tomioka, and Milan Vojnovic,  “QSGD: Communication-efficient SGD via gradient quantization and encoding,”  Advances in NeurIPS, vol. 30, 2017.
  * [16] Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste,  “Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks,”  Journal of Machine Learning Research, vol. 22, no. 241, pp. 1–124, 2021.
  * [17] Sebastian U Stich,  “Local sgd converges fast and communicates little,”  arXiv preprint arXiv:1805.09767, 2018.
  * [18] Farzin Haddadpour, Mohammad Mahdi Kamani, Mehrdad Mahdavi, and Viveck Cadambe,  “Local SGD with periodic averaging: Tighter analysis and adaptive synchronization,”  NeurIPS, vol. 32, 2019.
  * [19] Blake Woodworth, Kumar Kshitij Patel, et al.,  “Is local SGD better than minibatch SGD?,”  in ICML, 2020, pp. 10334–10343.
  * [20] Ahmed Khaled, Konstantin Mishchenko, and Peter Richtarik,  “Tighter Theory for Local SGD on Identical and Heterogeneous Data,”  in AISTATS. June 2020, pp. 4519–4529, PMLR,  ISSN: 2640-3498.
  * [21] Jianyu Wang and Gauri Joshi,  “Cooperative SGD: A unified framework for the design and analysis of communication-efficient SGD algorithms,”  arXiv preprint arXiv:1808.07576, 2018.
  * [22] Fan Zhou and Guojing Cong,  “On the convergence properties of a $k$-step averaging stochastic gradient descent algorithm for nonconvex optimization,”  arXiv preprint arXiv:1708.01012, 2017.
  * [23] Dongruo Zhou, Yuan Cao, and Quanquan Gu,  “Accelerated Factored Gradient Descent for Low-Rank Matrix Factorization,”  AISTATS, vol. 108, pp. 10, 2020.
  * [24] Léon Bottou, Frank E Curtis, and Jorge Nocedal,  “Optimization methods for large-scale machine learning,”  Siam Review, vol. 60, no. 2, pp. 223–311, 2018.
  * [25] Lam Nguyen, Phuong Ha Nguyen, Marten Dijk, Richtárik, et al.,  “SGD and Hogwild! convergence without the bounded gradients assumption,”  in ICML. PMLR, 2018, pp. 3750–3758.
  * [26] Srinadh Bhojanapalli, Anastasios Kyrillidis, and Sujay Sanghavi,  “Dropping convexity for faster semi-definite optimization,”  in Conference on Learning Theory, 2016, pp. 530–582.
  * [27] E. Candes and T. Tao,  “Near-optimal signal recovery from random projections: Universal encoding strategies?,”  IEEE transactions on information theory, vol. 52, no. 12, pp. 5406–5425, 2006.
  * [28] Benjamin Recht, Maryam Fazel, and Pablo A Parrilo,  “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,”  SIAM review, vol. 52, no. 3, pp. 471–501, 2010.
  * [29] Ophelia Crawford, Barnaby van Straaten, Daochen Wang, Thomas Parks, Earl Campbell, and Stephen Brierley,  “Efficient quantum measurement of pauli operators in the presence of finite sampling error,”  Quantum, vol. 5, pp. 385, 2021.