# Semi-device-dependent blind quantum tomography

Ingo Roth<EMAIL_ADDRESS>Dahlem Center for Complex Quantum Systems,
Freie Universität Berlin, Germany  Jadwiga Wilkens  Dahlem Center for Complex
Quantum Systems, Freie Universität Berlin, Germany  Dominik Hangleiter  Dahlem
Center for Complex Quantum Systems, Freie Universität Berlin, Germany  Jens
Eisert  Dahlem Center for Complex Quantum Systems, Freie Universität Berlin,
Germany  Mathematics and Computer Science, Freie Universität Berlin, Germany
Helmholtz-Zentrum Berlin für Materialien und Energie, 14109 Berlin, Germany

###### Abstract

Extracting tomographic information about quantum states is a crucial task in
the quest towards devising high-precision quantum devices. Current schemes
typically require measurement devices for tomography that are a priori
calibrated to a high precision. Ironically, the accuracy of the measurement
calibration is fundamentally limited by the accuracy of state preparation,
establishing a vicious cycle. Here, we prove that this cycle can be broken and
the fundamental dependence on the measurement devices significantly relaxed.
We show that exploiting the natural low-rank structure of quantum states of
interest suffices to arrive at a highly scalable blind tomography scheme with
a classically efficient post-processing algorithm. We further improve the
efficiency of our scheme by making use of the sparse structure of the
calibrations. This is achieved by relaxing the blind quantum tomography
problem to the task of de-mixing a sparse sum of low-rank quantum states.
Building on techniques from model-based compressed sensing, we prove that the
proposed algorithm recovers a low-rank quantum state and the calibration
provided that the measurement model exhibits a restricted isometry property.
For generic measurements, we show that our algorithm requires a close-to-
optimal number measurement settings for solving the blind tomography task.
Complementing these conceptual and mathematical insights, we numerically
demonstrate that blind quantum tomography is possible by exploiting low-rank
assumptions in a practical setting inspired by an implementation of trapped
ions using constrained alternating optimization.

###### Contents

  1. I Introduction
    1. I.1 Provable blind tomography via sparse de-mixing
    2. I.2 Practical blind tomography
    3. I.3 Related work and applications in signal processing
  2. II Quantum state tomography with imperfect Pauli correlation measurements
  3. III Formal problem definition
  4. IV Algorithm
    1. IV.1 Hard-thresholding algorithms: Ease and hardness of the projection
    2. IV.2 Relaxing the blind tomography problem: sparse de-mixing
    3. IV.3 Details of the SDT algorithm
    4. IV.4 Blind tomography via alternating least-square optimization
  5. V Recovery guarantees
  6. VI Numerical results
    1. VI.1 $\operatorname{GUE}$ measurements
    2. VI.2 Sub-sampled random Pauli measurements
    3. VI.3 Pauli measurements with coherent single-qubit errors
  7. VII Summary and outlook
  8. A Hardness of projection
  9. B Convergence proof
  10. C RIP guarantee for Hermitian random matrices

##  I Introduction

The development of quantum technologies is arguably one of the most vivid
scientific endeavours of current times. This development is faced with a
daunting challenge: To achieve the promising advantages of those technologies
one must engineer individual quantum components with an enormous precision.
The main limiting factors in implementing the many existing theoretical
proposals for exciting applications in quantum computing today are the
achievable noise levels and scalability of the components.

From an engineering perspective, improving such _noisy intermediate scale
quantum (NISQ) devices_ Preskill (2018); Acin  _et al._ (2017) requires
advanced and flexible _diagnostic techniques_ to extract actionable advice on
how to improve the device in the engineering cycle. Such diagnostic schemes
must meet tight practical constraints in terms of their complexity as well as
the required accuracy of the used devices. One of the most basic diagnostic
tasks is the extraction of tomographic information about quantum states from
experimentally measured data. Indeed, at the heart of every quantum
computation is the preparation of a quantum state. Quantum state tomography
can therefore provide valuable information for improving quantum devices
beyond a mere benchmarking of their correct functioning Eisert  _et al._
(2019).

However, in any such endeavour one encounters the following fundamental
challenge: In order to arrive at an accurate state estimate, most tomography
schemes rely on measurement devices that are calibrated to a very high
precision. At the same time, a precise and detailed characterization of a
measurement device requires an accurate state preparation. But improving the
accuracy of the state preparation using tomographic information was our goal
to begin with. We are trapped in a vicious cycle. This vicious cycle, depicted
in Figure 1, constitutes a fundamental obstacle to the improvement of quantum
devices.

Figure 1:  In the quest to engineer high fidelity quantum technologies one
encounters a vicious cycle: Extracting _actionable advice_ to correct for
error in the state preparation requires accurate quantum state estimation. The
accuracy of a state estimate crucially relies on the precise calibration of
the measurement device. But the calibration can ultimately only be tested and
improved if high fidelity quantum states are provided.

A make-or-break question is therefore: Is there any hope to break this cycle
at all? In other words, can one perform quantum state tomography _blindly_ ,
that is, without full knowledge of the measurement to begin with? More
specifically, can one simultaneously infer a quantum state and learn certain
unknown calibration parameters of the measurement device in a _self-
calibrating tomography scheme_ Brańczyk  _et al._ (2012)? A simple parameter
count indicates that this is typically impossible by just measuring a set of
mutually orthogonal observables: While an arbitrary quantum state in a
$d$-dimensional Hilbert space is characterized by $d^{2}-1$ many real
parameters, at the same time, the number of linearly independent measurements
in this space implies that we can learn at most $d^{2}$ independent
parameters. This leaves room for a single additional (calibration) parameter
only. This prohibits even slightly relaxing the requirement of a complete and
accurate characterization of the measurement device towards a partially
uncalibrated device. Tomography of an arbitrary quantum state is therefore
typically intrinsically measurement _device-dependent_ in this sense.

In this work, we break this vicious cycle. We observe that the above naïve
parameter count misses the point that in reasonably controlled quantum
devices, commonly encountered quantum states exhibit a natural structure: they
are close to being pure. We leverage this natural property to prove that one
can simultaneously learn an unknown calibration and a low-rank quantum state.
We thus arrive at what we coin a _semi-device-dependent_ scheme in which the
dependence on the measurement apparatus is significantly softened.

In order to achieve this goal, we formulate the _blind tomography problem_ as
the recovery task of a highly structured signal. This allows us to exploit and
further develop a powerful formal machinery from modern signal processing to
propose a scalable self-calibrating state tomography scheme that comes with
mathematically guarantees. We use a general model of the measurement device
which applies in a variety of relevant experiments: we model the measurements
as depending linearly on the unknown parameters of the possible calibration
errors. Indeed, in many situations the daunting uncertainty about the device
calibration is small and can be approximated as a linear deviation from an
empirically known calibration baseline.

Our scheme makes a trade-off between the dependence on the measurement device
and the state preparation device explicit and allow to optimally exploit this
dependence in a practical scheme. It is an intriguing feature of our results
that while structural assumptions on the quantum state to be learned typically
allow for _more efficient_ solutions Gross  _et al._ (2010); Flammia  _et al._
(2012); Kalev  _et al._ (2015); Riofrio  _et al._ (2017); Steffens  _et al._
(2017), here, structural assumptions allows one to solve a task in settings
where it _could not be solved at all_ in the absence of this assumption.

Going further, we exploit yet another structure to significantly extend the
realm of applicability and efficiency of our scheme, namely the _sparsity_ of
the calibration. Physically, this structural property amounts to the
assumption that only a small number out of the many possible calibration
errors has occurred in the specific experiment. In our scheme we therefore
simultaneously exploit the low-rank structure of the quantum state and
sparsity of the calibration coefficients to overcome the vicious tomography
cycle and provide rigorous guarantees with a favourable scaling in terms of
both the system dimension and the number of calibration errors.

###  I.1 Provable blind tomography via sparse de-mixing

Let us be slightly more formal in order to give an overview over the technical
methods and contributions of this work. In mathematical terms, the blind
tomography task that we solve is to infer a vector $\xi$ of $n$ calibration
parameters and a rank-$r$ quantum state $\rho$ from data of the form

$y=\mathcal{B}_{\xi}(\rho)=\mathcal{A}(\xi\otimes\rho)$ (1)

where $\mathcal{B}:\xi,\rho\mapsto\mathcal{B}_{\xi}(\rho)$ is a bi-linear map
describing the measurement model. The measured data $y$ might for example be
estimates for the expectation values of observables or probabilities of POVM
elements. For the time being, we ignore the error of the estimates induced by
finite statistics. It is convenient to regard the data as associated to a
structured linear estimation problem: we can equivalently model the
measurement map as a linear map $\mathcal{A}$ acting on $\xi\otimes\rho$.

Such structured linear inverse problems are studied in the mathematical
discipline of model-based _compressed sensing_ Baraniuk  _et al._ (2010);
Foucart and Rauhut (2013), where efficient algorithms with analytical
performance guarantees have been developed. A work horse of compressed sensing
that most rapidly solve the relevant inverse problems are so-called _iterative
hard-thresholding (IHT)_ algorithms Blumensath and Davies (2008).

As a first result of this work, we establish that the key step of an IHT
algorithm that solves the blind tomography problem is $\mathsf{NP}$-hard. To
overcome this obstacle, we propose an IHT algorithm that solves a slightly
relaxed version of the blind tomography problem: the task of de-mixing a sum
of $n$ different low-rank quantum states $\rho_{i}$, i.e., data of the form

$y=\mathcal{A}\left(\sum_{i=1}^{n}\xi_{i}e_{i}\otimes\rho_{i}\right),$ (2)

where $\\{e_{i}\\}_{i=1}^{n}$ denotes the standard orthonormal basis. An
efficient IHT algorithm for the de-mixing problem of low-rank matrices was
developed and analyzed in Ref. Strohmer and Wei (2019). This algorithm can be
readily adapted to our problem.

But relaxing the blind tomography problem to the de-mixing problem
artificially introduces an overhead in the number of unknown degrees of
freedom of the problem scaling as $2drn$, and in particular linearly with the
number of calibration parameters in the model. This leads to an unfavourable
situation in a two-fold manner: First, determining a large number of
calibration parameters also requires a large number of measurement settings as
the cost per calibration parameter scales with the dimension $d$ of the
quantum system. Second, a necessary condition for a well-posed blind de-mixing
problem of rank-$r$ with a maximal number of $d^{2}$ linearly independent
measurements of the form (2) is that there are more linearly independent
measurements than real parameters, i.e., $2rnd\leq d^{2}$. This means that the
simultaneous determination of a certain number of calibration parameters $n$
can in principle only work for sufficiently large system dimension $d$ in many
situations. This causes severe constraints in the achievable self-calibration
for small system sizes.

We argue that an additional well-motivated structural assumption can render
the blind tomography much more broadly applicable. Our argument is based on
the observation that the problem of determining an accurate estimate of the
quantum state in the blind setting involves solving two distinct sub-problems:
first, one needs to determine which ones of many potential error models of the
measurement contribute. Second, one needs to estimate the calibration
parameters of these models. Generically, there are many potential models that
parametrize, for instance, the deviation of every imperfect implementation of
a fixed measurement setting from its ideal implementation.

In this case, the first problem becomes combinatorially costly since many
distinct measurement settings need to be simultaneously calibrated. In
contrast, in our approach, it is straightforward to solve both tasks
simultaneously and even avoid a combinatorical overhead using the built-in
relaxations of compressed sensing. To this end, we observe that allowing for
many potential errors with associated calibration parameters only a small
number $s$ of which contribute amounts to assuming that the calibration vector
$\xi$ is $s$-sparse, i.e., it has only $s$ non-vanishing entries. Of course,
we do not assume that we know the support of the vector $\xi$. This falls
naturally into the framework of structured signal recovery. To summarize: we
observe data generated by linear measurements acting on $\xi\otimes\rho$ where
$\xi$ is an $s$-sparse vector and $\rho$ is a rank $r$ quantum state.

We are now faced with the recovery problem of de-mixing a _sparse sum_ of
different _low-rank quantum states_. We show that the projection onto this
structure can be efficiently calculated using hierarchical thresholding Roth
_et al._ (2016a) and therefore circumvents our $\mathsf{NP}$-hardness result.
We derive the corresponding iterative hard-thresholding algorithm and prove
that it successfully recovers the states $\rho_{i}$ and the sparse vector
$\xi$ provided that the measurement map $\mathcal{A}$ acts isometrically on
sparse sums of low-rank states. We further show that generic measurement
ensembles with $m$ different measurement settings exhibit this restricted
isometry property provided that $m$ scales at least as $srd+s\log n$. Thus, we
find that our algorithm solves the blind tomography problem with an overhead
in the required number of measurements that scales linearly in $s$ as compared
to the number of degrees of freedom in the problem given by $rd+s$. In
particular, the number of potential calibration models $n$ enters only
logarithmically in the measurement complexity of the scheme. This renders the
scheme highly scalable in $n$ providing flexibility in the modelling of
systematic measurement errors or calibration corrections. Furthermore, it
leaves sufficiently many linearly independent parameters to allow one to infer
a couple of calibration parameters already for comparably small system sizes.
We demonstrate the performance of the algorithm for the physically relevant
case of measuring Pauli operators that are locally mixed with the unknown
calibration parameters.

###  I.2 Practical blind tomography

Going beyond working out the theoretical guarantees, we numerically
demonstrate the functioning of the scheme and the mindset behind it.
Specifically, we show that the iterative hard-thresholding algorithm solves
the blind tomography problem from much fewer samples than competing methods
from generic (Gaussian) measurements as well as sub-sampled random Pauli
measurements. We then take the theoretical model to the practical testbed and
turn to a realistic model of measurement errors given by a coherent over-
rotation along some axis. Those measurements have significantly more
structure. We observe that the measurement structure together with the
sparsity constraints causes the SDT algorithm to frequently get stuck at
objective variables with an incorrect support. For this reason, we also study
the performance of a more pragmatically minded optimization strategy, namely,
constrained alternating minimization that does not require the relaxation to
the de-mixing problem. We numerically demonstrate that the blind tomography
problem in a realistic setting can be solved using this adapted algorithmic
approach. Thereby, we show that exploiting the low-rank structures of quantum
states allows to perform tomography blindly in realistic calibration and
measurement models. Therefore, we expect our approach to be directly
applicable in a variety of experimental settings that are practically relevant
in the quantum technologies.

###  I.3 Related work and applications in signal processing

DIDDsemiDIsemiDD Self- testing Pawłowski and Brunner (2011); Liang  _et al._
(2011); Li  _et al._ (2011, 2012) Blind tomography Standard tomography Figure
2:  Illustration of the spectrum between fully device-independent (DI) and
fully device-dependent (DD) quantum system characterization methods such as
self-testing and standard tomography, respectively. Semi-device-independent
(semi-DI) methods relax the stringent requirements of full device-
independence. Self-calibrating tomography relaxes the assumptions on the
calibration of the measurement device and therefore exemplifies a semi-device-
dependent (semi-DD) scheme. The blind tomography scheme presented here is an
example of such a semi-device-dependent scheme.

In our semi-device-dependent, self-calibrating scheme we aim at softening the
requirements of fully device-dependent schemes that crucially rely on perfect
measurement apparata. Coming from the opposite end, in quantum communications
introducing mild assumptions such as bounds on the system dimension Gallego
_et al._ (2010), one can weaken the impractical stringency of full device
independence to semi-device independence Pawłowski and Brunner (2011); Liang
_et al._ (2011); Li  _et al._ (2011, 2012). Device-independent and device-
dependent approaches can, thus, be seen as the extreme ends of an axis that
quantifies the amount of assumptions on the measurement device, see Figure 2
for an illustration. Hand in hand with reducing the amount of assumptions and
gaining robustness to imperfections, the amount of novel information that can
be gained is dramatically reduced. Semi-device-independent schemes move away
from the requirements of full device-independence towards more practical
settings but are still extremely demanding in terms of the required resources
and acceptance criteria. Our semi-device-dependent tomography scheme lies in
the opposite regime. It slightly relaxes the assumptions on the precision of
the measurement apparatus but still extracts tomographic information.

Self-calibrating tomography schemes have been previously proposed in specific
contexts using different methods and assumptions as a leverage to break the
vicious blind tomography cycle. In Ref. Mogilevtsev (2010) it has been argued
that single photon detectors can be simultaneously calibrated during state
tomography under the assumptions that the state is squeezed extending the
mindset of Ref. Mogilevtsev  _et al._ (2009); see also Ref. Mogilevtsev  _et
al._ (2012) for a more extensive discussion of potential classes of states and
the recent Ref. Sim  _et al._ (2019) for error bars in this context. Ref.
Brańczyk  _et al._ (2012) has reported experimental demonstration of
simultaneously reconstructing a quantum state together with certain unknown
unitary rotations associated to the measurement device via maximum likelihood
estimation in a linear optics setting. Here, the term ‘self-calibrating
tomography’ was coined. From a practical perspective, this work is perhaps
closest in mindset to the current work. Complementing and going significantly
beyond this work, here we prove that such a self-calibrating approach works
under very mild and natural structural assumptions and give rigorous
guarantees.

Another general-purpose framework is the _Gram matrix completion_ proposed in
Refs. Stark (2012, 2014). Here, a correlation matrix encoding information
about the measurement, the state and the measured data is completed from a
subset of known indices. So-called data-pattern tomography Řeháček  _et al._
(2010) avoids the calibration of the measurement device by comparing the data
to previously determined signatures of well-controlled reference states such
as coherent states Motka  _et al._ (2014). More conceptually speaking, schemes
incorporating model selection for quantum state tomography can be also viewed
as self-calibrating Ferrie (2014).

Another set of approaches focuses on characterizing entire gate sets or their
respective unitary errors self-consistently from the observed statistics when
applying different sequences Merkel  _et al._ (2013); Blume-Kohout  _et al._
(2013); Greenbaum (2015); Blume-Kohout  _et al._ (2017); Cerfontaine  _et al._
(2019). These methods typically rely on a certain design of the measurement
sequences and cost-intensive classical post-processing. In contrast, our model
of the measurement devices as linearly depending on a set of calibration
parameters is much simpler and requires far less resources.

Our work builds-on and further develops compressed-sensing techniques for the
tomography of quantum devices. Previous compressed-sensing schemes for quantum
tomography reduce the effort in the data acquisition while still ensuring an
efficient classical post-processing Gross  _et al._ (2010); Gross (2011); Liu
(2011); Kalev  _et al._ (2015); Kueng (2015); Kabanava  _et al._ (2016). These
schemes come with theoretical guarantees and have as well been successfully
employed in experiments Steffens  _et al._ (2017); Riofrio  _et al._ (2017);
Shabani  _et al._ (2011). The practical applicability of compressed sensing
tomography schemes rests on their robustness and stability against various
imperfections of the experimental setup. Small deviations from the compressive
model assumption and additive errors to the measurement outcomes, e.g. induced
by finite statistics, are reflected in a proportional and only slightly
enhanced estimation error. Still, the schemes rely on measurement devices that
are calibrated to very high precision with the notable exception of
compressive tomography schemes for quantum processes that use randomized
benchmarking data Kimmel and Liu (2017); Roth  _et al._ (2018). Here, we relax
this requirement using a semi-device-dependent approach. In distinction, in
the previous schemes low-rank assumptions were considered to reduce the
complexity of a tomography scheme, giving rise to an important _quantitative
improvement_. Here, those assumptions are expected to often make blind
tomography possible in the first place and therefore permit even a
_qualitative improvement_ over the known schemes.

Recovery problems of the form (1) or the related de-mixing problem (2) also
arise in other disciplines. For example, these problems appear in future
mobile communication scenarios with the promise to yield much more scalable
protocols with respect to the number of served devices Wunder  _et al._
(2015). More specifically, our work can be applied in order to extend the
_internet-of-things setup_ described in Ref. Strohmer and Wei (2019) in case
one additionally wants to exploit the sporadic (sparse) user activity of
machine-type messaging. Furthermore, our work identifies yet another set of
hierarchical signal structures that allow for an efficient projection: It
extends the work on compressed sensing with hierarchically sparse signals of a
subset of the authors to low-rank matrices Roth  _et al._ (2016b, a).

The remainder of this work is organized as follows. In the subsequent Section
II, we give a detailed description of a concrete experimental setup that
motivates our mathematical formulation of the blind tomography problem. In
Section III, we provide the formal definitions of the blind tomography problem
and introduce the notation used in the subsequent parts of the work. The
details of the sparse demixing algorithm and its variant based on alternating
optimisation are derived in Section IV. On the way, we establish the
$\mathsf{NP}$-hardness of the projection associated to the original blind
tomography problem. The theorems guaranteeing the performance of the sparse
demixing algorithm are explained in Section V. The corresponding proofs are
given in the appendix. Finally, numerical simulations of the algorithms
performance and its application to practical use cases are shown in Section VI
before we conclude with an outlook in Section VII.

##  II Quantum state tomography with imperfect Pauli correlation measurements

So far, our description of the measurement scheme has been fairly abstract. In
the following, we describe a concrete scenario in which our formalism applies.
Consider an ion trap experiment preparing a multi-qubit quantum state $\rho$.
We perform Pauli correlation measurements, i.e., we estimate $m$ expectation
values of the form

$y_{0}^{(k)}=\mathcal{A}_{0}(\rho)^{(k)}=\operatorname{Tr}\left[\rho\left(W^{(k)}_{1}\otimes
W^{(k)}_{2}\otimes\cdots\otimes W^{(k)}_{l}\right)\right],$ (3)

where $W^{(k)}_{j}\in\\{X,Y,Z,\operatorname{Id}\\}$ is a Pauli matrix acting
on the $j$th qubit and $k\in[m]\coloneqq\\{1,2,\ldots,m\\}$. We refer to
$\mathcal{A}_{0}:\mathbb{C}^{d\times d}\to\mathbb{R}^{m}$ as the measurement
map or sampling operator 111 Note that one might actually implement projective
measurements in the multi-qubit Pauli basis as done, e.g., in Refs. Riofrio
_et al._ (2017); Steffens  _et al._ (2017). While such projective measurements
contain more information than the Pauli correlation measurement, we restrict
ourselves to Pauli expectation values both for the sake of simplicity and to
remain in a setting for which theoretical guarantees can be proven Gross  _et
al._ (2010); Liu (2011)..

In many experimental setups, it is natural to implement measurements of a
certain Pauli observable – in the case of ion traps Pauli $Z$ – while the
other Pauli observables require more effort. A measurement of any other Pauli
observable – in the case of ion traps Pauli $X$ and Pauli $Y$ – can then be
implemented by applying a suitable sequence of unitary gates prior to the
measurement. For example, using addressed laser pulses of different duration
one can implement rotations around different axes and thus implement the
Hadamard gate $H$ as well as the phase gate $S$. In this way, one can realize
measurements in the $X=HZH$ and $Y=SHZHS^{\dagger}$ basis.

But each application of an additional gate may come with a coherent error _in
addition_ to the native error associated with the measurement itself. In this
way, we end up with different systematic errors for different Pauli
observables parametrized by the angles $\theta,\varphi$ of a coherent error
given by $\mathrm{e}^{\mathrm{i}\theta X}\mathrm{e}^{\mathrm{i}\varphi Z}$.
This gives rise to some probability of actually measuring the expectation
value of another local Pauli matrix than the targeted one. For example,
consider a coherent error given by a (small) rotation around the $Z$-axis as
given by $\mathrm{e}^{i\varphi Z}$. The faulty implementation of the Hadamard
gate is then given by $\tilde{H}=\mathrm{e}^{i\varphi Z}H$. Of course, the
native $Z$-measurement is untouched by this coherent error, since no unitary
rotation precedes this measurement. However, instead of $Y$ one now actually
measures
$\tilde{Y}=S\tilde{H}Z\tilde{H}^{\dagger}S^{\dagger}=\cos(2\varphi)Y+\sin(2\varphi)X$.
At the same time $X$ remains undisturbed.

More generally, we can introduce calibration parameters $\xi_{W\to\tilde{W}}$
measuring the strength of the error that replaces a certain target Pauli
matrix $W$ by $\tilde{W}$. For instance, in the above example those parameters
are given by $\xi_{Y\to Y}=\cos(2\varphi)$, $\xi_{Y\to X}=\sin(2\varphi)$ and
$\xi_{Z\to Z}=\xi_{X\to X}=1$. For simplicity, we assume that these
calibration parameters are identical for different qubit registers. Assuming
that errors are not too large, the calibration parameters for the target
measurement fulfil $\xi_{W\to W}\approx 1$ or all $W\in\\{X,Y,Z\\}$. This
leaves us with six independent calibration parameters corresponding to the
cross-contributions. To construct the measurement map $\mathcal{A}{}$, we
start from the definition of the target measurement $\mathcal{A}_{0}$ in (3).
From $\mathcal{A}_{0}$ we can derive calibration measurement components
$\mathcal{A}_{W\to\tilde{W}}$ appearing with the coefficient
$\xi_{W\to\tilde{W}}$ by replacing all appearances of the Pauli matrix $W$ in
the definition of $\mathcal{A}_{0}$ with $\tilde{W}$. If $W$ appears in a
multi-qubit Pauli observable several times the resulting observable is the sum
of all Pauli observables generated by replacing only one of the $W$ by
$\tilde{W}$, assuming that the coherent errors are small so that the higher-
order terms can be neglected. For example, a faulty realization of the
observable $ZYZZY$ is now given by $\xi_{Y\rightarrow Y}ZYZZY+\xi_{Y\to
X}(ZXZZY+ZYZZX)$.

Altogether, to linear order in the calibration parameters
$\xi_{W\to\tilde{W}}$ with $W\neq\tilde{W}$ we end up constructing a
description of the effective faulty measurement by

$y=\xi_{0}\mathcal{A}_{0}(\rho)+\sum_{W\neq\tilde{W}\in\\{X,Y,Z\\}}\xi_{W\to\tilde{W}}\mathcal{A}_{W\to\tilde{W}}(\rho),$
(4)

which can be written as linear map $\mathcal{A}$ action on $\xi\otimes\rho$
with $\xi=[\xi_{0},\xi_{X\to Y},\xi_{X\to Z},\ldots,\xi_{Z\to Y}]^{T}$. By
assumption, we set $\xi_{0}=1$.

In this measurement model the sparsity assumption is justified if unitary
errors in a certain coordinate plane are dominant compared to others thus
singling out certain types of calibration measurement components. Importantly,
we do not assume that we know which corrections are dominant (i.e., the
support of $\xi$) _a priori_. The measurement model also exemplifies a setting
in which one is ultimately limited to measuring a maximal set of $d^{2}$
observables. Thus, blind tomography becomes only possible exploiting structure
assumptions if one does not allow for different ways of implementing the same
measurement that yield different calibration corrections without introducing
too many new calibration parameters.

##  III Formal problem definition

Motivated by this example, we set out to provide a formal definition of the
blind tomography problem and the related sparse-de-mixing problem. The
notation and terminology introduced in this section allows us to formulate a
general signal-processing framework using which the blind tomography can be
provably solved. Both, the blind tomography and the sparse de-mixing, problems
are linear inverse problems that feature a combination of different
compressive structures. These are smaller sets of linear vector spaces and it
will be convenient to introduce some notation to refer to these sets. The
prototypical example is the _set of $s$-sparse real vectors_

$\Sigma_{s}^{n}\coloneqq\\{\xi\in\mathbb{R}^{n}\mid|\operatorname{supp}\xi|\leq
s\\}\subset\mathbb{R}^{n},$ (5)

which is defined by the support $\operatorname{supp}\xi$ of a vector $\xi$,
i.e., the index set of the non-vanishing entries of $\xi$, having cardinality
smaller or equal than $s$. The set of $s$-sparse vectors is not a vector space
itself but the union of $\binom{n}{s}$ $s$-dimensional subspaces.

In the realm of quantum mechanics, the non-commutative analog of sparse
vectors, namely low-rank matrices, is important. We denote the _set of complex
rank $r$ matrices_ by

$\mathbb{C}^{d\times d}_{r}\coloneqq\\{x\in\mathbb{C}^{d\times
d}\mid\operatorname{rank}x\leq r\\}.$ (6)

Since we are dealing with quantum states we will restrict our attention to the
set $\mathcal{D}^{d}\subset\mathbb{C}^{d\times d}$ of trace-normalized,
positive semidefinite matrices, i.e., $\rho\geq 0$ and
$\operatorname{Tr}\rho=1$ for all $\rho\in\mathcal{D}^{d}$. Our results can be
straightforwardly generalized to general matrices without these constraints.
We denote the set of rank-$r$ quantum states as
$\mathcal{D}^{d}_{r}=\mathcal{D}^{d}\cap\mathbb{C}^{d\times d}_{r}$. In
particular, $\mathcal{D}^{d}_{1}$ is the set of pure quantum states. In order
to solve the blind tomography problem we need to simultaneously recover an
$s$-sparse real vector $\xi$ and a rank-$r$ quantum state $\rho$. It is
convenient to regard both $\xi$ and $\rho$ as a combined signal
$X=\xi\otimes\rho$ and model the measurement including its dependence on the
calibration parameter as a linear map $\mathcal{A}$ acting on $X$. Considering
such linear maps instead of bi-linear maps is sometimes refered to as
‘lifting’ in the compressed sensing literature Ahmed  _et al._ (2014). For a
physicist, ‘lifting’ is also the natural isomorphism at the heart of the
density matrix formulation of quantum mechanics. The signal $X$ is highly
structured as it is a tensor product of a sparse vector and a low-rank quantum
state. We denote the set of all potential signals as

$\Omega^{n,d}_{s,r}\coloneqq\\{\xi\otimes x\mid\xi\in\Sigma^{n}_{s},\
x\in\mathcal{D}^{d}_{r}\\}\subset\mathbb{C}^{nd\times d}.$ (7)

One can regard a signal $X\in\Omega^{n,d}_{s,r}$ as an $nd\times d$ matrix
consisting of $n$ blocks of size $d\times d$ stacked on top of each other as
depicted in Figure 3, where each $d\times d$ block is proportional to the same
quantum state $\rho$ and only $s$ of the blocks are non-vanishing.

$\mathbb{C}^{nd\times d}$

= $\Omega^{n,d}_{s,r}$

+…+=$\hat{\Omega}^{n,d}_{s,r}$

Figure 3:  The signal sets of the blind tomography and sparse de-mixing
problem can be regarded as subsets of $\mathbb{C}^{nd\times d}$, i.e.,
matrices consisting of $n$ blocks of $d\times d$. For a blind tomography
signal in $\Omega^{n,d}_{s,r}$, only $s$ out of the $n$ blocks are non-zero
and are proportional to the same rank $r$ matrix. In contrast, a signal of the
sparse de-mixing problem in $\hat{\Omega}^{n,d}_{s,r}$ comprises $s$ non-
vanishing blocks with potentially different rank $r$ matrices.

We are now equipped to concisely state the problem we would like to study.

######  Problem 1 (Blind tomography).

Let $\mathcal{A}:\mathbb{C}^{nd\times d}\to\mathbb{R}^{m}$ be a linear map.
Given data $y=\mathcal{A}(X)\in\mathbb{R}^{m}$ and the linear map
$\mathcal{A}$, recover $X$ under the assumption that

$X\in\Omega^{n,d}_{s,r}$

Our approach to algorithmically solving the blind tomography problem makes use
of a proxy problem: we relax it to signals that are a bit less restrictively
structured

$\hat{\Omega}^{n,d}_{s,r}\coloneqq\left\\{\sum_{i=1}^{n}\xi_{i}e_{i}\otimes
x_{i}\ \middle|\ \xi\in\Sigma^{n}_{s},\ x_{i}\in\mathcal{D}^{d}_{r}\ \forall
i\in[n]\right\\}.$ (8)

Both sets $\hat{\Omega}^{n,d}_{s,r}$ and $\Omega^{n,d}_{s,r}$ are subsets of
$\mathbb{C}^{nd\times d}$. The difference between them as illustrated in
Figure 3 is the following: While for $X\in\Omega^{n,d}_{s,r}$ all $d\times d$
blocks are proportional to the same quantum state $x$, we allow the $d\times
d$ blocks of $\hat{X}\in\Omega^{n,d}_{s,r}$ to be proportional to _different
quantum states_ $x_{i}$. Again only $s$ out of the $n$ blocks of $\hat{X}$ are
non-vanishing. Analogously to Problem 1, we define the linear inverse problem
associated with $\hat{\Omega}$.

######  Problem 2 (Sparse de-mixing).

Let $\mathcal{A}:\mathbb{C}^{nd\times d}\to\mathbb{R}^{m}$ be a linear map.
Given data $y=\mathcal{A}(X)\in\mathbb{R}^{m}$ and the linear map
$\mathcal{A}$, recover $X$ under the assumption that

$X\in\hat{\Omega}^{n,d}_{s,r}.$

The observed data of the sparse de-mixing problem can be equivalently
described as

$y=\sum_{k=1}^{n}\xi_{k}\mathcal{A}_{k}(x_{k}),$ (9)

where we have split up $X$ into trace-normalized $d\times d$ blocks $x_{k}$
and their norm $\xi_{k}$ according to the definition of
$\hat{\Omega}^{n,d}_{s,r}$. Correspondingly, we can decompose the linear map
$\mathcal{A}$ into the set of linear maps $\\{\mathcal{A}_{k}\\}_{k=1}^{n}$
where each $\mathcal{A}_{k}$ acts only on the $k$-th $d\times d$ block of $X$.
From this reformulation it becomes clear that the problem amounts to
reconstructing a set of low-rank signals $\\{x_{k}\\}_{k}$ from observing its
sparse mixture under linear maps, hence the name sparse de-mixing.

For both the blind-tomography and the sparse-de-mixing problem, we
alternatively write each of the $n$ linear maps $\mathcal{A}_{k}$ in terms of
$m$ observables

$\\{A_{k}^{(i)}\in\mathbb{C}^{d\times
d}\mid(A_{k}^{(i)})^{\dagger}=A_{k}^{(i)}\\}_{i=1}^{m}$ (10)

via

$\mathcal{A}_{k}(x_{k})^{(i)}=\langle A_{k}^{(i)},x_{k}\rangle$ (11)

with the Hilbert-Schmidt inner product $\langle
X,Y\rangle=\operatorname{Tr}(X^{\dagger}Y)$.

Note that as long as we consider Hermitian matrices for the measurement
$A_{k}^{(i)}$ and signals $x_{i}$, we end up with an real data vector
$y\in\mathbb{R}^{m}$. For applications other than quantum tomography it is
straightforward to adopt our proofs and results to real signals or complex-
valued measurement maps. Furthermore, for the sake of simplicity we have
formulated both recovery problems without noise. More generally, the data can
be assumed to be of the form $y=\mathcal{A}(X)+\epsilon$ where $\epsilon$
denotes additive, e.g. statistical, noise.

In the following, we will also make use of the inner product of vectors
$x,y\in\mathbb{R}^{n}$ defined as $\langle x,y\rangle=\sum_{i}x_{i}y_{i}$,
their $\ell_{2}$-norm $\|x\|_{\ell_{2}}=\sqrt{\langle x,x\rangle}$ and the
Frobenius norm a matrix $X\in\mathbb{C}^{d_{1},d_{2}}$ induced by the Hilbert-
Schmidt inner product $\|X\|_{F}=\sqrt{\langle X,X\rangle}$.

##  IV Algorithm

We now turn to the technical derivation of our algorithm for the blind quantum
tomography and the sparse de-mixing problem. Our algorithm builds on
primitives developed in the field of compressed sensing. In particular, we
generalize the hard thresholding algorithm to accommodate the structural
assumptions of both problems. As a first step, we establish the hardness of
direct thresholding approaches to the blind tomography problem before stating
a tractable algorithm for the sparse de-mixing problem.

Let us be more precise: the blind quantum tomography problem requires
different assumptions on two levels. First, we want the signal to be a tensor
product $\xi\otimes\rho$, i.e., of rank one. Second, both tensor factors are
assumed to be structured. Concretely, we assume $\xi$ to be $s$-sparse and
$\rho$ to be of rank $r$. We are therefore faced with low-rank structures on
two separate levels: first, the block-structured signal as given by the tensor
product of calibration vector and quantum state has unit rank. Second, by
assumption the target quantum states, i.e., the individual blocks of the
signal, have low rank.

It has been observed in the compressed sensing literature that multi-level
structures with structured tensor components can be notoriously hard to
reconstruct. One prototypical example of this is combined sparsity and low-
rankness in the sense that the signal is the tensor product of two sparse
vectors, i.e., $X=\xi\otimes\tau$ with $\xi,\tau\in\Sigma^{n}_{s}$. This
problem is already very similar to the blind tomography problem where one of
the sparse vectors is replaced by a low-rank matrix, the quantum state.

The obstacle arising from such structures can be understood from a different
perspective present in the compressed sensing literature that is related to
different algorithmic approaches. The perhaps most prominent approach in
compressed sensing is the convex relaxation of structure-promoting
regularizers yielding efficient convex optimization programs. Minimizing the
$\ell_{1}$-norm or the Schatten-$1$-norm is known to solve linear inversion
problems involving sparse or low-rank vectors efficiently and sampling
optimal, respectively. However, simply combining both regularizers in a convex
fashion does not yield a sampling-optimal reconstruction of problems that
feature both structures anymore Oymak  _et al._ (2015).

###  IV.1 Hard-thresholding algorithms: Ease and hardness of the projection

Another algorithmic approach used in compressed sensing are so-called _hard
thresholding algorithms_ such as CoSAMP, IHT or HTP Needell and Tropp (2008);
Blumensath and Davies (2008); Foucart (2011); see also the textbook Foucart
and Rauhut (2013) for an introduction. These are typically iterative
procedures that minimize the deviation from the linear constraints in some way
or other, e.g. by gradient descent, and in each iteration project onto the
structure of the signal. For many compressed sensing problems this is possible
because even though recovery problems, such as

$\operatorname*{minimize}_{\xi}\|\mathcal{A}(\xi)-y\|_{\ell_{2}}^{2}\quad\text{subject
to $\xi\in\Sigma^{n}_{s}$},$ (12)

are $\mathsf{NP}$-hard Magdon-Ismail (2017), the related projection

$P_{\Sigma_{s}^{n}}(\tau)\coloneqq\operatorname*{arg\,min}_{\xi\in\Sigma_{s}^{n}}\|\xi-\tau\|_{\ell_{2}}$
(13)

can be computed efficiently. For the given example of projecting onto
$s$-sparse vectors, this solution is given by the hard-thresholding operation
defined as follows: Let $\Sigma_{\text{max}}$ be the set of indices of the $s$
absolutely largest entries of $\tau$. Then,

$(P_{\Sigma_{s}^{n}}(\tau))_{i}=\begin{cases}\tau_{i}&\text{for
$i\in\Sigma_{\text{max}}$}\\\ 0&\text{otherwise}.\end{cases}$ (14)

In words, one keeps the largest entries of $\tau$ and replaces the other
entries by zero. Analogously, the projection of Hermitian matrices onto low-
rank matrices can be efficiently calculated by calculating the eigenvalue
decomposition and applying $P_{\Sigma_{r}^{d}}$ to the eigenvalue vector. Let
$X\in\mathbb{C}^{d\times d}$ be a Hermitian matrix with eigenvalue
decomposition $X=U\operatorname{diag}(\lambda)U^{\dagger}$. We define the
projection onto positive semi-definite low-rank matrices as

$P_{\mathcal{D}_{r}^{d}}(X)=U\operatorname{diag}(P_{\Sigma_{r}^{d}}(\lambda|_{\geq
0}))U^{\dagger},$ (15)

where $\lambda|_{\geq 0}$ denotes the restriction of $\lambda$ to its non-
negative entries.

In hard-thresholding algorithms, the problems associated with simultaneously
exploiting sparse and low-rank structures are manifest in the computational
hardness of computing the respective projections. For the case of unit rank
matrix with sparse singular vectors, calculating the projection is the so-
called _sparse PCA_ problem, i.e., given a matrix $A\in\mathbb{R}^{n\times n}$

$\operatorname*{minimize}_{\xi,\tau\in\Sigma^{n}_{s}}\
\|A-\xi\otimes\tau\|_{F}.$ (16)

Indeed, exactly solving this problem in the worst case is $\mathsf{NP}$-hard
by a trivial reduction to the CLIQUE problem Magdon-Ismail (2017). But it
turns out that the hardness is much worse: one can even make average-case
hardness statements based on conjectures regarding the hardness of the
_planted clique problem_ Berthet and Rigollet (2013a, b); Brennan and Bresler
(2019). Moreover, the sparsePCA problem remains just as hard even when one
merely asks for an approximation up to a _constant relative error_ Chan  _et
al._ (2016); Magdon-Ismail (2017).

As the first technical result of this work, we show that also the projection
onto $\Omega^{n,d}_{s,r}$ is an $\mathsf{NP}$-hard problem by reducing it to
the sparse PCA problem.

######  Theorem 3 (Hardness of constrained minimization).

There exists no polynomial time algorithm that calculates

$\operatorname*{minimize}\quad\|A-X\|_{F}\quad\text{subject to
$X\in\Omega^{n,d}_{s,r}$},$ (17)

for all $A\in\mathbb{C}^{nd\times d}$ unless ${\mathsf{P}}={\mathsf{NP}}$.
This still holds for $s=n$.

The details of the proof are given in Appendix A. This hardness result
provides a strong indication that a straightforward adaptation of compressed
sensing techniques is not feasible. In this work, our way out of this is to
sacrifice sampling optimality of the algorithm for a lower runtime and being
able to prove analytical performance guarantees. Alternating minimization
approaches that make the factorization explicit is also a viable way forward.
Via provide a detailed description of such an algorithm in Section IV.4. But
proving global recovery guarantees for non-convex algorithms typically becomes
much more involved.

###  IV.2 Relaxing the blind tomography problem: sparse de-mixing

In fact, the bi-sparse and low-rank structure can be relaxed to a simple
hierarchical sparsity constraint Wunder  _et al._ (2018a); Foucart  _et al._
(2019). A vector $\xi\in\mathbb{C}^{Nn}$ consisting of $N$ blocks of size $n$
is called $(s,\sigma)$-hierarchically sparse if it has at most $s$ blocks with
non-vanishing entries, that themself are $\sigma$-sparse Sprechmann  _et al._
(2010); Friedman  _et al._ (2010); Sprechmann  _et al._ (2011); Simon  _et
al._ (2013). For this structure a hard-thresholding algorithm together with
theoretical recovery guarantees has been derived in Refs. Roth  _et al._
(2016a, b, 2018); Groß  _et al._ (2020). It has been applied in different
contexts Wunder  _et al._ (2017, 2018b, 2018c, 2019) including _sparse blind
deconvolution_ Wunder  _et al._ (2018a) which features the combined low-rank,
sparse structure.

Here, we make use of this approach to solve the blind quantum tomography
problem as formalized in Problem 1. At the heart of our approach is the
insight that the projection onto $\hat{\Omega}^{n,d}_{s,r}$ can be efficiently
computed since the $n$ $d\times d$ blocks may be different. This allows one to
combine the projection onto $\Sigma^{n}_{s}$ and the projection onto
$\mathcal{D}^{d}_{r}$: First, the low-rank projection
$P_{\mathcal{D}^{d}_{r}}$ is applied to each of the $d\times d$ blocks of the
input matrix $X$. Subsequently, the sparse projection operator is applied by
setting the $n-s$ smallest blocks in Frobenius norm to zero. The resulting
algorithm is summarized as Algorithm 1.

Algorithm 1 Projection onto $\hat{\Omega}^{n,d}_{s,r}$

1:$X\in\mathbb{C}^{nd\times d}.$

2:$Y=0$

3:for $k\in[n]$ do

4: $n_{k}=\|P_{\mathcal{D}^{d}_{r}}(x_{k})\|_{F}$

5:end for

6:$W=\operatorname{supp}P_{\Sigma^{n}_{s}}(n)$

7:$Y_{W}=X_{W}$.

8:$Y$ is projection of $X$ onto $\hat{\Omega}^{n,d}_{s,r}$.

The computational cost of the projection onto $\hat{\Omega}^{n,d}_{s,r}$ is
dominated by the eigenvalue decomposition required to compute the low-rank
approximation $\mathcal{D}^{d}_{r}$ of each block. Computing the full
eigenvalue decomposition of the $d\times d$ blocks requires computation time
of $O(d^{3})$ using, e.g., Householder reflections Golub and van Loan (1989).
Since we are only interested in the dominant $r\ll d$ eigenvalues, the effort
can be reduced to $O(rdw)$ using the Lanczos algorithm, where $w$ is the
average number of non-zero elements in a row of a block Golub and van Loan
(1989). Using randomized techniques one might be able to further reduce the
computational costs Halko  _et al._ (2011). The calculation of the Frobenius
norms contributes $O(nd^{2})$ flops. The largest blocks can be selected using
the quick-select algorithm Hoare (1961) in $O(n)$. Note that the low-rank
projections and Frobenius norms of all blocks can also be performed in
parallel without any modification to the algorithm.

Equipped with an efficient projection for $\hat{\Omega}_{s,r}^{n,d}$, we can
construct a structured iterative gradient descent algorithm. This is a variant
of the IHT algorithm, that was originally developed for sparse vectors
Blumensath and Davies (2008). The IHT algorithm is a projective gradient
descent algorithm that iteratively alternates gradient steps to optimize the
$\ell_{2}$-norm deviation between the data and a projection onto the
constraint set. The resulting recovery algorithm for the _sparse de-mixing
(SDT)_ problem is stated as Algorithm 2.

Algorithm 2 SDT-algorithm

1:Data $y$, measurement $\mathcal{A}$, sparsity $s$ and rank $r$ of signal

2:Initialize $X^{0}=0$.

3:repeat

4: Calculate step-widths $\mu^{l}$

5:
$X^{l+1}=P_{\hat{\Omega}^{n,d}_{s,r}}\left(X^{l}+\operatorname{diag}(\mu^{l})P_{\mathcal{T}_{X^{l}}}\left(\mathcal{A}^{\dagger}\left(y-\mathcal{A}(X^{l})\right)\right)\right)$

6:until stopping criterion is met at $l=l^{\ast}$

7:Recovered signal $X^{l^{\ast}}$

The SDT algorithm is closely related to the IHT algorithm for de-mixing low-
rank matrices that was developed in Ref. Strohmer and Wei (2019). We will
refer to this algorithm as the _DT algorithm_. The main difference between our
SDT and the DT algorithm of Ref. Strohmer and Wei (2019) is that the latter
does not make the additional sparsity assumptions on the signal. For this
reason, the SDT algorithm differs in the projection
$P_{\hat{\Omega}^{n,d}_{s,r}}$ that additionally applies the projection
$P_{\Sigma^{n}_{s}}$ selecting the $s$ dominant blocks. In fact, in the
special case of considering non-sparse signals in $\hat{\Omega}^{n,d}_{n,r}$
the SDT algorithm coincides with the DT algorithm.

###  IV.3 Details of the SDT algorithm

To be fully self-contained, let us now go through the individual steps of the
SDT algorithm and specify the relevant details. Every iteration of the
algorithm starts with the computation of
$G^{l}=\mathcal{A}^{\dagger}(y-\mathcal{A}(X^{l}))$, the gradient for the
$\ell_{2}$-norm deviation
$f(X)=\frac{1}{2}\|y-\mathcal{A}(X)\|_{\ell_{2}}^{2}$ evaluated at $X^{l}$.
The algorithm subsequently employs a modification from Ref. Wei  _et al._
(2016) in calculating the steepest gradient inspired by geometrical
optimization techniques which leads to a faster convergence Absil  _et al._
(2009); Vandereycken (2013): The set of rank $r$ matrices is an embedded
differential manifold in the linear vector space of all matrices. Thus, a
direction on this embedded manifold is characterized by a tangent vector on
the manifold. While this geometry straight-forwardly generalizes to the set of
$nd\times d$ matrices with rank $r$ blocks, due to sparsity constraint
$\hat{\Omega}^{n,d}_{s,r}$ fails to be a differential manifold. Nonetheless,
we can make use of tangent vectors as ‘natural’ search directions in our
optimization problem for the non-vanishing blocks of $X^{l}$ that are conform
with a fixed rank constraints.

The tangent space of rank $r$ matrices at point $x$ is given by the set of
matrices that share the same column or row space $x$ Absil  _et al._ (2009).
Correspondingly, the tangent space projection of a non-vanishing block of $X$
can be defined as follows: Let $x_{k}=U_{k}\Lambda_{k}U^{\dagger}_{k}$ be the
eigenvalue decomposition of the $k$-th block of $X$ with $\Lambda_{k}$ the
diagonal matrix with eigenvalues in decreasing order. Further, let
$U^{(r)}_{k}$ denote the restriction of $U_{k}$ to its first $r$ columns
corresponding to the range of $x_{k}$. Then, the tangent space projection
acting on $g_{k}$ the $k$-th block of $G$ is given by

$P_{\mathcal{T}_{X}}(G)_{k}=g_{k}-(\operatorname{Id}-P_{U})g_{k}(1-P_{U}),$
(18)

with $(P_{U})_{k}=U^{(r)}_{k}(U^{(r)}_{k})^{\dagger}$. The entire tangent-
space projection $P_{\mathcal{T}_{X}}(G)$ is defined by acting trivially on
the blocks of $G$ corresponding to vanishing blocks of $X$ and as the
projection (18) otherwise.

As we prove below in generic situations the SDT algorithm converges for a
constant step-width set to $\mu_{l}=1$ and even without using the tangent
space projection. Empirically, a faster convergence is achieved with the
tangent space projection and using the following prescription for the step-
width calculation: From the projected gradient
$G^{l}_{P}=P_{\mathcal{T}_{X}}(G^{l})$ in the $l$-th iteration we then
calculate the algorithms step width for each block individually as

$\mu_{k}^{l}=\frac{\|(G_{P}^{l})_{k}\|_{F}^{2}}{\|\mathcal{A}((G_{P}^{l})_{k})\|_{\ell_{2}}^{2}}$
(19)

and multiply each block by the corresponding $\mu_{k}^{l}$. In order to have a
compact notation, we introduce the diagonal matrix
$\operatorname{diag}(\mu^{l})=\operatorname{diag}(\mu^{1}_{l},\ldots\mu_{1}^{l},\mu_{2}^{l},\ldots,\mu_{2}^{l},\ldots,\mu_{n}^{l})$
where each step width is repeated $d$ times. The new state of the algorithm,
$X^{l+1}$, is given by the projection of the result of a gradient step with
step width $\mu^{l}$ onto the set $\hat{\Omega}_{s,r}^{n,d}$.

Finally, we have to specify a stopping criterion at which the loop of the
algorithm is exited. We terminate the algorithms if the objective function is
below a specified threshold, i.e.,

$\frac{\|y-\mathcal{A}(X^{l})\|_{\ell_{2}}}{\|y\|_{\ell_{2}}}\leq\gamma_{\text{break}}$
(20)

or a maximal number of iteration is reached. If the data vector $y$ has
additive noise, $\gamma_{\text{break}}$ has to be chosen to be larger than the
expected norm of the noise. To be less relying on expectations on the noise
levels, one can alternatively make use of criteria on the gradient and step
width or test for oscillating patterns in the identified support.

###  IV.4 Blind tomography via alternating least-square optimization

A more direct algorithmic approach to the blind tomography problem is to use a
constrained alternating least square (ALS) optimization. In ALS optimization,
one performs a constrained optimization of the objective function

$\displaystyle
f_{\text{ALS}}(\xi,\rho)=\frac{1}{2}\|y-\mathcal{A}(\xi,\rho)\|_{\ell_{2}}^{2},$
(21)

with respect to one of the two variables while regarding the respective other
variable as constant in an alternating fashion, see Algorithm 3.

Algorithm 3 ALS-BT algorithm

1:Data $y$, measurement $\mathcal{A}$, sparsity $s$ and rank $r$ of signal

2:Initialize $\rho^{0}$.

3:repeat

4: $\quad\xi^{l}=\operatorname*{arg\,min}_{\xi\in\Sigma^{n}_{s}}\
f_{\text{ALS}}(\xi,\rho^{l-1})$

5: $\rho^{l}=\operatorname*{arg\,min}_{\rho\in\mathbb{C}^{d\times d}_{r}}\
f_{\text{ALS}}(\xi^{l},\rho)$

6:until stopping criterion is met at $l=l^{\ast}$

7:Recovered signal $\rho^{l^{\ast}}$, $\xi^{l^{\ast}}$

We perform the optimization over $\Sigma^{n}_{s}$, Algorithm 3 Step 4, using
the standard IHT algorithm for sparse vector recovery. Note that calculating
the linear measurement map for $\xi$ given a fixed $\rho$ simply involves
evaluating all calibration measurement blocks individually, i.e. calculating
$\mathcal{A}_{i}(\rho)$ for all $i\in[N]$. Analogously, the low-rank
optimization over $\mathbb{C}^{d\times d}_{r}$, Algorithm 3 Step 5, can be
performed with iterative hard-thresholding on the manifold of low-rank
matrices. A detailed description of a suitable algorithmic implementation is
given by Algorithm 2 in the special case of a single matrix block, i.e.
$N,s=1$. Computing the corresponding linear map acting on $\rho$ for fixed
$\xi$ amounts to summing up the individual measurement blocks weighted by
their corresponding calibration coefficient.

The ALS optimization requires an initialization with a suitable $\rho^{0}$ in
order to evaluate the first objective function for optimizing $\xi$. One
method that we found viable is to randomly draw a rank-$r$ state using Haar-
random eigenvectors. Note that, in general, constrained ALS optimization can
be highly sensitive to the chosen initialization. For this reason, depending
on the measurement map and calibration model, alternative initialization
strategies might become necessary. As break-off criteria we can again use a
bound on the objective function as in (20) and an allowed maximal number of
iterations.

##  V Recovery guarantees

We now prove that for certain simple measurement ensembles, the SDT algorithm
converges to the optimal solution before we numerically demonstrate its
performance in the following section. More precisely, following the outline of
model-based compressed sensing Baraniuk  _et al._ (2010); Blumensath and
Davies (2009), the SDT-algorithm can be accompanied with recovery guarantees
based on a _restricted isometry property_ (RIP) of the measurement ensemble
that is custom-tailored to the structure at hand. Intuitively, it seems clear
that a measurement map should at least in principle allow for solving the
associated linear inverse problem uniquely if it acts as an isometry on
signals from the constraint set. So-called RIP constants formalize this
intuition:

######  Definition 4 ($\hat{\Omega}_{s,r}$-RIP).

Given a linear map $\mathcal{A}:\mathbb{C}^{nd^{2}}\to\mathbb{C}^{m}$, we
denote by $\delta_{s,r}$ the smallest $\delta\geq 0$ such that

$(1-\delta)\|x\|_{F}^{2}\leq\|\mathcal{A}(x)\|_{\ell_{2}}^{2}\leq(1+\delta)\|x\|_{F}^{2}$
(22)

for all $x\in\hat{\Omega}_{s,r}$.

The constant $\delta_{s,r}$ measures how much the action of $\mathcal{A}$ when
restricted to elements of $\hat{\Omega}_{s,r}$ deviates from that of an
isometry. Correspondingly, if $\delta_{s,r}$ is sufficiently small we expect
this to be sufficient to ensure that the restricted action of $\mathcal{A}$
becomes invertible. In fact, if a measurement map has a sufficiently small RIP
constant one can prove the convergence of projective gradient descent
algorithms to the correct solution of the structured linear inverse problem.
For the sake of simplicity, we analyze the SDT-algorithm omitting the tangent
space projection and also assuming a constant step widths $\mu^{l}=1$. In
numerically simulations we observe that making use of the tangent space
projection and a more sophisticated heuristic for the step width yields faster
convergence and better recovery performance. But the RIP assumption is in fact
strong enough to already for this simpler algorithmic variant ensure that the
following theorem holds:

######  Theorem 5 (Recovery guarantee).

Let $\mathcal{A}:\mathbb{C}^{nd\times d}\to\mathbb{C}^{m}$ be a linear map and
suppose that the following RIP condition for $\mathcal{A}$ holds

$\delta_{3s,3r}<\frac{1}{2}.$ (23)

Then, for $X\in\hat{\Omega}_{s,r}$, the sequence $(X^{l})$ defined by the SDT-
Algorithm (Algorithm 2) with $\mu^{l}=1$ and
$P_{\mathcal{T}_{X^{l}}}=\operatorname{Id}$ with $y=\mathcal{A}(X)$ satisfies,
for any $l\geq 0$,

$\left\|X^{l}-X\right\|_{F}\leq\gamma^{l}\left\|X^{0}-X\right\|_{F},$ (24)

where $\gamma=2\delta_{3s,3r}<1$

We establish that the SDT-Algorithm converges to the correct solution of the
sparse de-mixing problem at a rate that is upper bounded by the RIP constant
$\delta_{3s,3r}$ of the measurement map. The right hand side of the RIP
condition (23) is not expected to be optimal. Typically, one can at least
improve the bound to $\frac{1}{\sqrt{3}}$ with a slightly more complicated
argument Foucart and Rauhut (2013). Since we are interested in the parametric
scalings here, we choose to present a simpler argument at the cost of slightly
worse constants. Furthermore, the statement of Theorem 5 does not account for
statistical noise or potential mild violation of the signal constraints. Again
following standard techniques Foucart and Rauhut (2013), we expect that a more
complicated noise- and model-robust version of Theorem 5 can be derived. But
for the current analysis, we are content with the significantly simpler
version.

The pressing next question is, of course, which measurement ensembles actually
exhibit the required RIP. Interestingly, it is notoriously hard to give
deterministic constructions of measurement maps that are sample optimal and
feature the RIP. In fact, already for the RIP for $s$-sparse vectors there are
no sample optimal deterministic measurement maps known to date Foucart and
Rauhut (2013). To further complicate the state of affairs, it is also known to
be $\mathsf{NP}$-hard to check whether a given measurement map exhibits the
$s$-sparse RIP with RIP constant small than a given $\delta$ Bandeira  _et
al._ (2013).

For this reason, the field of compressed sensing uses probabilistic
constructions to arrive at provably sampling optimal measurement maps. Using a
random ensemble of measurement maps of sampling optimal dimension one
establishes that with high probability a randomly drawn instance will exhibit
the RIP property. In other words, one proves that the originally hard linear
inverse problem typically becomes easy for a certain measurement ensemble. The
arguably simplest measurement ensemble consists of observables given by i.i.d.
chosen random Gaussian matrices. In our setting a fully Gaussian measurement
map can be constructed from a set of $\\{A_{i}\in\mathbb{R}^{nd\times
d}\\}_{i=1}^{m}$ of $m$ Gaussian matrices with entries draws i.i.d. from the
normal distribution $\mathcal{N}(0,1)$ and defining
$y^{(l)}=\operatorname{Tr}(A_{i}X)$.

As a toy model for quantum tomography it is more natural to consider
observables drawn from a random ensemble of Hermitian matrices such as the
Gaussian unitary ensemble ($\operatorname{GUE}$). Operationally, we define the
$\operatorname{GUE}$ by drawing a matrix $X$ with complex Gaussian entries,
$X_{k,l}\sim\mathcal{N}(0,1)+\mathrm{i}\mathcal{N}(0,1)$, and subsequently
projecting $X$ onto Hermitian matrices using
$P_{\Yboxdim{4pt}\,\yng(2)}{}:X\mapsto\frac{1}{2}(X+X^{\dagger})$. For
measurement maps from $\operatorname{GUE}$ we prove the following statement:

######  Theorem 6 ($\hat{\Omega}^{n,d}_{s,r}$-RIP for random Hermitian
matrices.).

Let $\\{{A}^{(k)}_{i}\\}_{i=1,k=1}^{n,m}$ be a set of Hermitian matrices drawn
i.i.d. from the $\operatorname{GUE}$. Let $\mathcal{A}$ be the measurement
operator defined by $\\{{A}^{(k)}_{i}\\}_{i=1,k=1}^{n,m}$ via Eqs. (9) and
(11). Then $\frac{1}{\sqrt{m}}\mathcal{A}$ satisfies the
$\hat{\Omega}^{n,d}_{s,r}$-RIP with parameter $\delta_{s,r}$ with probability
at least $1-\tau$ provided that

$m\geq\frac{C}{\delta_{s,r}^{2}}\left[s\ln\frac{\mathrm{e}n}{s}+(2d+1)rs\ln\frac{c}{\delta}+\ln\frac{2}{\tau}\right]$
(25)

for sufficiently large numerical constants $C,c>0$.

The proof of the theorem is provided in Appendix C. Based on the result for
random Hermitian measurement maps we now discuss the asymptotic scaling of the
measurement complexity of our approach to the blind tomography problem and the
sparse de-mixing problem. First, the derived measurement complexity (25) is in
accordance with the degrees of freedom of signal
$X\in\hat{\Omega}^{n,d}_{s,r}$. The second term of $O(drs)$ corresponds to the
number of degrees of freedom specifying the $s$ rank-$r$ matrices of dimension
$d$. The first term of $O(s\ln n)$ is the minimal sampling complexity in $s$
for learning the $s$ non-trivial entries and their support Foucart and Rauhut
(2013). Second, in analogy, we expect the optimal number of measurements for
the blind tomography problem, i.e., reconstructing signals in
$\Omega_{s,r}^{n,d}$ instead of $\hat{\Omega}_{s,r}^{n,d}$, to scale as
$O(s\ln n+dr)$. Hence, having a provably efficient algorithms capable of
solving the blind tomography as well as the sparse de-mixing problem comes at
the cost of an increase in the sampling complexity by an additional factor of
$s$ in the second term of the sampling complexity. Most importantly, invoking
the sparsity assumption on the calibration vector $\xi$ allows us to get away
without a linear increase $n$ of the number of calibration parameters. Thus,
the overhead in measurement complexity of our approach to the blind tomography
problem is relatively mild.

In fact, the measurement complexity derived for Gaussian measurements can
often be used as a guideline for the sampling complexity of other measurement
ensembles that are also sufficiently unstructured. However, the proof
techniques for model-based compressed sensing that exploit the combination of
different structures are not easily translatable to other measurement
ensembles. An exception are measurement ensembles that feature a structure
that is sufficiently aligned with the signal structure such as the one
exploited in Ref. Roth  _et al._ (2016a) for hierarchically sparse signals. We
leave the study of more involved measurement ensembles to future work.

##  VI Numerical results

The analytical results of the previous section provide worst-case bounds on
the asymptotic scaling for a class of idealized, unstructured measurements. In
order to benchmark and assess the non-asymptotic performance of compressed
sensing algorithms in practice, however, numerical simulations are
indispensable. In a first step we therefore perform numerical simulations for
the idealized measurement model as given by random $\operatorname{GUE}$
matrices, comparing the performance of our algorithm to related established
algorithms that do not entirely exploit the structure of the problem. In a
second step, we compare the SDT Algorithm 2 with standard CS tomography in a
blind tomography setting involving measurements of Pauli correlators, cmp.
(3). To do so we randomly draw subsets of the possible Pauli measurements as
possible calibrations $\mathcal{A}_{i}$ of the measurement apparatus. Finally,
we demonstrate the feasibility of blind tomography under structure assumptions
in the realistic measurement and calibration setting involving single-qubit
coherent errors described in Section II. To this end, we employ the Algorithm
3 that performs alternating constrained optimisation. The algorithms and the
scripts producing the plots have been implemented in Python and will be made
available under Ref. Wilkens  _et al._ (2020).

Figure 4:  The figure displays the recovery rate for the SDT, DT and informed
DT algorithm for different number of observables $m$ for GUE measurements.
Each point is averaged over $50$ random measurement and signal instances with
$r=1$, $d=16$, $n=10$ and $s=3$. A signal is considered successfully recovered
if its Frobenius norm deviation from the original signal is smaller than
$10^{-3}$. One observes nearly coinciding recovery performances for the
informed DT and the SDT algorithm. In comparison, the DT algorithm requires
significantly more observables for recovery.

###  VI.1 $\operatorname{GUE}$ measurements

The SDT algorithm goes beyond existing IHT algorithms for the de-mixing
problem of low-rank matrices in that it additionally allows to exploit a
sparse mixtures. We demonstrate that this yields a drastic and practically
important improvement in the number of measurement required for the
reconstruction.

To this end, we draw signal instances $X=\xi\otimes\rho$ at random from
$\Omega^{n,d}_{s,r}$. We use four qubit pure states
$\rho=\left.\left|{\psi}\right.\right\rangle\\!\\!\left.\left\langle{\psi}\right.\right|$
with $r=1$ and $d=16$, where $\left.\left|{\psi}\right.\right\rangle$ is drawn
uniformly (Haar) random from the complex $\ell_{2}$-norm sphere. The
calibration vector $\xi\in\mathbb{R}^{n}$ with $n=10$ has a support of size
$s=3$ drawn uniformly from the set of all $\binom{n}{s}$ possible supports.
The non-vanishing entries of $\xi$ are normal distributed with unit variance.
The measurements are drawn at random from the $\operatorname{GUE}$ ensemble as
defined above with a varying number of observables $m$.

The closest competitor to the SDT algorithm is the related algorithm of Ref.
Strohmer and Wei (2019). The algorithm of Ref. Strohmer and Wei (2019)
coincides with the special case of the SDT algorithm where we use the
projection on to $\hat{\Omega}^{n,d}_{n,r}$ with $s=n$ ignoring the sparsity
in the block structure. We will refer to this algorithm as the _DT algorithm_.
We can also give the DT algorithm the ‘unfair’ advantage of restricting the
problem to the correct block support of the signal from the beginning. We will
refer to this variant as the _informed DT algorithm_.

Figure 4 shows the recovery rate for the SDT algorithm, the DT algorithm and
its informed variant for different $m$. Each point is average over $50$ random
signal and measurement instances. We consider a signal as successfully
recovered if the distance of the algorithm’s output to the original signal is
smaller than $10^{-3}$ in Frobenius norm. The algorithm terminated if either
the stopping criterion (20) with $\gamma_{\text{break}}=10^{-5}$ is met or
after a maximal number of 600 iteration. We observe that if one of the
algorithm successfully recovers a signal it typically meets the stopping
criterion after less than 100 iterations.

The curves for all three algorithm in Figure 4 display a sharp phase
transition from a regime where the number of measurement is too small to
recover any signal to a regime of reliable recovery. While the phase
transition for the SDT algorithm appears in a similar regime to the informed
DT algorithm, the DT algorithm requires considerably more samples in order to
recover the signal instances.

We conclude that the sparsity of the calibration parameters can be exploited
by the SDT algorithm to considerably reduce the required number of
measurements. Even more so, this does not require many more sampling points as
compared to an algorithm which is given _a priori_ knowledge which errors were
present, that is, the block support of the signal. This shows that the SDT
algorithm solves the de-mixing and blind tomography task in a highly efficient
way and scalable. Finally, the number of possible erroneous measurements
$\mathcal{A}_{i}$ can be scaled up at a very low cost in terms of required
measurement settings.

Figure 5:  The figure displays the trace norm reconstruction error for the SDT
compared to the standard tomography algorithm for different number of
observables $m$ for sub-sampled random Pauli measurements. Each point is
averaged over $30$ random measurement and signal instances with $r=1$, $d=8$,
$n=10$ and $s=3$. The inline figure shows the mean $\ell_{2}$-norm
reconstruction error of the calibration coefficients for the SDT algorithm.
Figure 6:  The figure displays the trace norm reconstruction error for the SDT
compared to the standard tomography algorithm for different number of
observables $m$ for sub-sampled random Pauli measurements. Each point is
averaged over $30$ random measurement and signal instances with $r=1$, $d=8$,
$n=10$ and $s=4$. The inline figure shows the mean $\ell_{2}$-norm
reconstruction error of the calibration coefficients for the SDT algorithm.

###  VI.2 Sub-sampled random Pauli measurements

For the application in characterizing quantum devices, it is key to compare
the recovery performance of the SDT algorithm with standard low-rank quantum
tomography algorithms. To this end note that the SDT algorithm restricted to
$n,s=1$ is also a state-of-the-art algorithm for standard low-rank state
tomography without the on-the-fly calibration. Thus, we will make use of this
implementation of conventional low-rank state tomography in the following.

We draw signal instances as before but using three-qubit states,
$s\in\\{3,4\\}$ and altering the model for the calibration parameter: We set
the first entry of $\xi$ to $\xi_{0}=1$. The support of the remaining entries
is drawn uniformly at random. The non-vanishing entries are then i.i.d. taken
from the normal distribution rescaled by a factor of $1/10$. This mimics a
setting where we have a dominant target measurement and a couple of small
systematic deviation from a known set of candidates. The target measurements
as well as the systematic deviations are uniformly sub-sampled Pauli
observables. Thus, $\mathcal{A}_{0}$ till $\mathcal{A}_{n}$ have the form of
(3) with differently i.i.d. selected Pauli-observables uniformly selected from
$\\{\operatorname{Id},X,Y,Z\\}$. We simulate statistical noise using $10^{8}$
samples per expectation value in order to realistically limit the resolution
of the SDT algorithm.

We simultaneously perform recoveries with the SDT-algorithm using the entire
measurement matrix including the calibration measurement components and the
SDT-algorithm using only the target measurement $\mathcal{A}_{0}$ as in a
conventional tomography setting.

The resulting trace distance of the state estimate, i.e., the trace-normalized
first block of $X$, from the original $\rho$ is shown for different number of
measurements in Figure 5 and Figure 6 for different sparsities $s=3$ and
$s=4$, respectively. The curves indicate the median over the depicted 30
sample points per value of $m$. The inline plot of boths Figures further show
the $\ell_{2}$-norm deviation of the reconstructed calibration parameters and
the original $\xi$.

One observes that the conventional low-rank tomography becomes more accurate
with an increasing number of measurement but is asymptotically still bounded
from below by the systematic error induced by the calibration on the order of
$10^{-1}$. This agrees with the order of magnitude of variance of the
calibration coefficients. In contrast, the SDT-algorithm while performing
slightly worse in a regime of insufficient measurements outperforms the
conventional algorithm for a moderate number of samples and is ultimately only
limited by the statistical noise. However, in the parameter regime under
investigation there are even for large number of samples $m>150$ a small
number (well below 10%) of instance where SDT only reaches an accuracy
comparable to standard tomography. In these instances we find that the support
for the calibration measurement components was incorrectly identified. For
$s=4$ we furthermore observe one pathological instance of SDT for $m=240$ that
is worse in recovery than standard tomography is in this regime. For $s=4$ the
phase transition of SDT appears for a slightly larger values of $m$ compared
to $s=3$. The curves for the reconstruction error of the quantum state
approximately coincide with the curves for the error in the calibration
parameter. We conclude that for a sufficient number of measurement settings,
the SDT algorithm almost always performs a significantly more acurate state
reconstruction and simultaneously extracts the calibration parameters. The
precision is ultimately only limited by the statistical error in the
estimation of the expectation values.

Figure 7: The figure displays the trace norm reconstruction error for the ALS
compared to the standard tomography algorithm for different number of
observables $m$ for Pauli measurements with coherent single-qubit errors. Each
point is averaged over $50$ random measurement and signal instances with
$r=1$, $d=16$, $n=7$ and $s=2$. The inline figure shows the mean
$\ell_{2}$-norm reconstruction error of the calibration coefficients for the
SDT algorithm. Figure 8: The figure displays the trace norm reconstruction
error for the ALS compared to the standard tomography algorithm for different
number of observables $m$ for Pauli measurements with coherent single-qubit
errors. Each point is averaged over $30$ random measurement and signal
instances with $r=1$, $d=16$, $n=7$ and $s=3$. The inline figure shows the
mean $\ell_{2}$-norm reconstruction error of the calibration coefficients for
the SDT algorithm.

###  VI.3 Pauli measurements with coherent single-qubit errors

We now come back to the concrete realistic scenario described in Section II.
There we derived the calibration measurement model originating from coherent
errors in the gates that implement the single-qubit measurements.

For the numerical simulations, we draw a set of $m$ Pauli observables
uniformly at random as the target measurement. Subsequently, we introduce six
calibration blocks such that every observable in the set $\\{X,Y,Z\\}$ is
swapped with another Pauli observable in $\\{X,Y,Z\\}$ in a specific block. We
generate data $y$ for given states and calibration parameters using the linear
calibration measurement model without noise as induced by finite statistics.

We find that in the parameter regimes that are easily amenable to numerical
studies on desktop hardware the SDT algorithm is not capable of successfully
reconstructing the states when the calibration parameters for the corrections
are considerably smaller than the leading order measurement. To thoroughly
understand this limitation, in the following, we briefly report the
performance of the SDT algorithm on different sub-tasks related to the
recovery problem.

First, we choose $d=16$ and $n=s=1$ such that only a single block, either the
ideal measurement or one of the correction blocks, is used to generate the
signal from a random pure state ($r=1$). We observe that the SDT algorithm is
able to recover the signals in this standard tomography problem. This
indicates that also the calibration blocks individually allow for tomographic
reconstruction of low-rank states. Second, the SDT algorithm can discriminate
between different mixtures of the six correction blocks. To demonstrate this,
we ignore the ideal measurement and employ only the correction blocks to
generate the signal. We set the active calibration coefficients to one. Thus,
$n=6$, $s\leq n$ and $\xi_{i}=1$ for $i$ active. We observe that given a
sufficient number of measurement settings the SDT algorithm correctly
reconstructs pure states in this measurement setting. The same findings hold
true if the target measurement is again considered as long as the active
calibration coefficients are set to $1$. We observe successful reconstructions
of unit rank states for $n=7$ and $s\in{1,2,3}$.

A more natural setting however would typically have calibration coefficients
that are considerably smaller than the ideal measurement. This justifies the
linear expansion for the measurement model in the first place. If we choose,
e.g., $\xi_{i}=1/10$ for the indices $i$ of active blocks, we were unable to
identify a parameter regime on desktop hardware where the SDT algorithm can
successfully recover the majority of instances of pure states. We observe that
if the SDT algorithm settles on an objective variable with an incorrect block
support in the first few iteration, it is not able to subsequently run into
objective variables with a different block support in most instances.

Despite the negative result for the SDT algorithm in the most realistic
setting, the general mindset to exploit structure (low-rankness) to allow
quantum state tomography in a blind fashion is fruitful using a slightly
different algorithmic strategy.

To this end, we use the constrained _alternating least square_ (ALS) algorithm
described in Section IV.4. We set the first calibration coefficient
corresponding to the ideal measurement to one. The support of the remaining
active calibration coefficients is drawn uniformly at random and their value
are i.i.d. drawn from a shifted normal distribution with standard deviation
$0.05$ and mean value $0.2$. We use Haar random pure states, $r=1$ of a four-
qubit system, $d=16$, as the target states.

The algorithm is initialized with a Haar-randomly drawn pure state. We allow
for a maximal number of $1000$ iterations of the algorithm or terminate if the
criterion (20) with $\gamma_{\text{break}}=10^{-5}$ is met. Furthermore, if
the stopping criterion is not met after $50$ iterations, we re-initialize the
algorithm with a new random pure state. We allowed for a maximal number of
$10$ or $20$ re-initializations for $s=2$ and $s=3$, respectively. We observe
that in case of successful recovery typically at most $3$ re-initializations
are required with most instances already correctly converging from the initial
state.

As in the previous section, we compare the recovery performance of the ALS
with the standard low-rank tomography algorithm. The trace-norm error and
calibration error for different numbers of measurement settings for $s=2$ and
$s=3$ are displayed in Figure 7 and 8, respectively. We observe that, as
expected, the reconstruction error of standard low-rank tomography is again
lower-bounded by a scale set by the magnitude of the calibration parameters.
In contrast, with an only slightly larger number of measurement settings, the
constrained ALS algorithm is capable of recovering the states and the
calibration parameter with an accuracy that is improved by orders of magnitude
and in the noiseless scenario only limited by the algorithms stopping
criterion. Compared to recovery performance of the SDT algorithm we observe an
even sharper phase transition to the regime of recovery.

##  VII Summary and outlook

In this work, we have shown that the natural assumption of low-rankness allows
one to perform self-calibrating quantum state tomography. Relaxing the blind
tomography problem to a sparse de-mixing problem has allowed us to develop an
efficient classical post-processing algorithm, the SDT-algorithm, that is
theoretically guaranteed to recover both the quantum state and the device
calibration under a restricted isometry condition of the measurement model. We
have demonstrated the necessity of relaxing the blind tomography problem
within the framework of hard-thresholding algorithms by establishing the
$\mathsf{NP}$-hardness of the projection onto the set consisting of the outer
products of vectors and fixed-rank matrices. Introducing a sparsity assumption
on the calibration coefficients ensures that the reconstruction scheme can
already be applied for fairly small system dimension. We have explicitly
proven that a Gaussian random measurement model meets the required restricted
isometry condition with a close-to-optimal measurement complexity in $O(s\ln
n+drs)$. Furthermore, we have numerically demonstrated an improved performance
of the SDT-algorithm for random instances of measurement models compared to
previously proposed non-sparse de-mixing algorithms and standard low-rank
state tomography. While these generic measurement and calibration models
allows us to derive analytical guarantees, it is fair to argue that these
models might at best capture some aspects of actual experimental
implementations. A potential starting point for extending recovery guarantees
to more realistic settings is the generalization of our results to random
Pauli measurements as considered in Sec. VI Liu (2011) together with the
coherence measures and structured measurement guarantees developed in the
context of hierarchically spares signals Sprechmann  _et al._ (2011); Roth
_et al._ (2016a, 2018); Groß  _et al._ (2020).

To complement our conceptually and rigorously minded work with a more
pragmatic approach, we have additionally developed and implemented a
structure-exploiting blind tomography algorithm based on alternating
optimization. We have numerically demonstrated that the alternating algorithm
is able to perform self-calibrating low-rank tomography in a realistic
measurement and calibration model that is well-motivated by gate
implementations in ion traps. These numerical simulations indicate that the
approach to the blind tomography problem developed here might be well-suited
to improve tomographic diagnostics in current experiments. Ultimately, the
recovery performance of the proposed algorithms has to be evaluated on
experimental data. It is the hope that this work contributes to establishing a
new mindset in quantum system identification and specifically tomographic
recovery in which no component used has to be precisely known, but still under
physically meaningful structural assumptions, a mindset here referred to as
being semi-device-dependent.

##  Appendix A Hardness of projection

As a starting point we state the sparsePCA problem.

######  Problem 7 (SparsePCA).

_Input:_ Symmetric matrix $A\in\mathbb{R}^{n\times n}$, sparsity $s$, positive
real number $a>0$. _Question:_ Does there exist an $s$-sparse unit vector
$v\in\mathbb{R}^{n}$ with $v^{T}Av\geq a$?

It has been folklore for quite some time that the sparse PCA problem is
$\mathsf{NP}$-hard. A formal proof can be found in Ref. Magdon-Ismail (2017),
where the CLIQUE problem is encoded into instances of sparsePCA. From the
hardness of sparsePCA it follows that there does not exist a polynomial time
algorithm for the projection onto the set of symmetric, unit rank matrices
with sparse eigenvectors. Formally, we have:

######  Proposition 8 (Hardness of projection onto the set of symmetric, unit
rank matrices with sparse eigenvectors).

Given a matrix $A\in\mathbb{R}^{d\times n}$ and $s,\sigma\in\mathbb{N}$, there
exist no polynomial time algorithm that calculates

$\begin{split}&\operatorname*{minimize}\quad\|A-vw^{T}\|_{F},\\\
&\text{subject to $v\in\Sigma^{d}_{\sigma}$,
$w\in\Sigma^{n}_{s}$}.\end{split}$ (26)

This still holds for $\sigma=d$.

###### Proof.

It turns out to be sufficient to only consider the case where $\sigma=d$,
i.e., only one of the factors is required to be sparse. It is straightforward
to see that solving the problem with both vectors being sparse allows one to
solve the projection with only one sparse vector: Define

$A=\begin{pmatrix}\mathbf{0}_{d-\sigma,n}\\\ A^{\prime}\end{pmatrix}$ (27)

with $\mathbf{0}_{a,b}$ being an $a\times b$ matrix filled with zeros. It then
holds that

$\min_{v\in\Sigma^{d}_{\sigma},w\in\Sigma^{n}_{s}}\|A-vw^{T}\|_{F}=\min_{v^{\prime}\in\mathbb{C}^{\sigma},w\in\Sigma^{n}_{s}}\|A^{\prime}-v^{\prime}w^{T}\|_{F}.$
(28)

We now embed the sparsePCA problem. To do so we first make the normalization
of the vectors $v,w$ in the optimization problem explicit to it to a
maximization problem over normalized vectors:

$\begin{split}&\min_{v\in\mathbb{R}^{\sigma},w\in\Sigma^{n}_{s}}\|A-vw^{T}\|^{2}_{F}\\\
=&\min_{\lambda\in\mathbb{R},v\in\mathbb{R}^{\sigma}\cap
B_{\ell_{2}}^{\sigma},w\in\Sigma^{n}_{s}\cap B_{\ell_{2}}^{n}}\|A-\lambda
vw^{T}\|^{2}_{F}\end{split}$ (29)

with $B_{\ell_{2}}^{n}=\\{v\in\mathbb{R}^{n}\mid\|v\|_{\ell_{2}}\leq 1\\}$ the
$\ell_{2}$-norm ball. Solving the optimization problem over $\lambda$ yields

$\displaystyle\min_{\lambda\in\mathbb{R}}\|A-\lambda vw^{T}\|^{2}_{F}$ (30)
$\displaystyle=\min_{\lambda\in\mathbb{R}}\left\\{\|A\|_{F}^{2}+\lambda^{2}\|v\|_{\ell_{2}}^{2}\|w\|_{\ell_{2}}^{2}-2\lambda\langle
w,Av\rangle\right\\}$ (31)
$\displaystyle=\|A\|_{F}^{2}-\min_{v\in\mathbb{R}^{\sigma}\cap
B_{\ell_{2}}^{\sigma},w\in\Sigma^{n}_{s}\cap B_{\ell_{2}}^{n}}\langle
w,Av\rangle^{2}.$ (32)

Since $A$ is fixed we conclude that the optimization problem (29) is
equivalent to

$\begin{split}&\qquad\mathrm{maximize}\ |\langle w,Av\rangle|\ \\\
&\text{subject to $v\in\mathbb{R}^{\sigma}\cap B_{\ell_{2}}^{\sigma}$,
$w\in\Sigma^{n}_{s}\cap B_{\ell_{2}}^{n}$}.\end{split}$ (33)

Furthermore, using the Cauchy-Schwarz inequality we find that

$\max_{v\in\mathbb{R}^{\sigma}\cap
B_{\ell_{2}}^{\sigma},w\in\Sigma^{n}_{s}\cap B_{\ell_{2}}^{n}}|\langle
v,Aw\rangle|=\max_{w\in\Sigma^{n}_{s}\cap B_{\ell_{2}}^{n}}\|Aw\|_{\ell_{2}}.$
(34)

Now consider an instance of the sparsePCA problem with a symmetric input
matrix $B\in\mathbb{R}^{n\times n}$, sparsity $s$ and $a>0$. W.l.o.g. we can
assume that $B$ is a positive matrix since solving the sparsePCA problem for
the $B-\min{}\\{0,\lambda_{\min{}}(B)\\}\operatorname{Id}$ shifted by the
smallest eigenvalue $\lambda_{\min{}}(B)$ of $B$ and $a$ shifted
correspondingly, allows one to solve the sparsePCA problem for $B$. For a
positive matrix $B$ we find a factorization $B=A^{T}A$. Hence, deciding
whether the maximum over all $w\in\Sigma_{s}^{n}$ of $w^{T}Bw$ is larger than
$a$ is solved by calculating the maximum of $\|Aw\|_{\ell_{2}}^{2}=w^{T}Bw$.
This completes the reduction. ∎

We are now prepared to tackle our related problem: the projection onto
$\Omega^{n,d}_{s,r}$. We have the following statement:

######  Theorem 9 (Hardness of constrained minimization).

There exist no polynomial time algorithm that calculates for all
$A\in\mathbb{C}^{n\times n}$:

$\operatorname*{minimize}\quad\|A-X\|_{F}\quad\text{subject to
$X\in\Omega^{n,d}_{s,r}$},$ (35)

unless ${\mathsf{P}}{}={\mathsf{NP}}{}$. This still holds for $s=n$.

We note that our result for exactly computing straightforwardly generalizes to
the case of approximating the target function up to constant relative error
using results on the approximatability of sparsePCA Chan  _et al._ (2016).

###### Proof.

Suppose the existed an efficient algorithm that determines the objective value
of the projection (35). To encode the sparsePCA problem, we choose an instance
of $A$ as follows: Let $A^{\prime}\in\mathbb{R}^{n\times d}$ be a matrix and
let $A^{\prime}_{i}$ denote the $i$-th row of $A$. Let $e_{i}$ be the basis
vectors $(e_{i})_{j}=\delta_{i,j}$, with $\delta_{i,j}$ the Kronecker symbol.
We choose $A=\sum_{i=1}^{n}e_{i}\otimes\operatorname{diag}(A^{\prime}_{i})$,
where $\operatorname{diag}(A^{\prime}_{i})$ denotes the diagonal matrix with
the $i$-th row of $A^{\prime}$ on its diagonal. Furthermore, we define
$a^{\prime}=(A^{\prime}_{1},\ldots,A^{\prime}_{n})\in\mathbb{R}^{nd}$ to be
the vector arising by concatenating all rows of $A$. By definition an
$X\in\Omega^{n,d}_{s,r}$ can be decomposed as $X=\xi\otimes\rho$ with
$\xi\in\Sigma^{n}_{s}$ and $\rho\in H^{d}_{r}$. Let
$\rho=U\operatorname{diag}(\lambda)U^{\dagger}$ the eigenvalue decomposition
of $\rho$ with a suitable unitary $U\in U(n)$ and $\lambda$ the vector of its
eigenvalues. Then, we can rewrite

$\displaystyle\|A$
$\displaystyle-\xi\otimes\rho\|_{2}^{2}=\sum_{i=1}^{n}\|\operatorname{diag}(A^{\prime}_{i})-\xi_{i}\rho)\|_{2}^{2}$
$\displaystyle=\|\operatorname{diag}(a^{\prime})-(\operatorname{Id}_{n}\otimes
U)\operatorname{diag}(\xi\otimes\lambda)(\operatorname{Id}_{n}\otimes
U^{\dagger})\|_{2}^{2}$
$\displaystyle=\|\operatorname{diag}(a^{\prime})(\operatorname{Id}_{n}\otimes
U)-(\operatorname{Id}_{n}\otimes
U)\operatorname{diag}(\xi\otimes\lambda)\|_{2}^{2}$
$\displaystyle=\sum_{i,j=1}^{nd}|A^{\prime}_{i}-(\xi\otimes\lambda)_{j}|^{2}|(\operatorname{Id}_{n}\otimes
U)_{i,j}|^{2},$

where we have used the unitary invariance of the $\ell_{2}$-norm in the third
step. We can introduce the doubly stochastic matrix $W$ with entries
$W_{k,l}=|U_{k,l}|^{2}$ and relax the optimization to

$\displaystyle\min_{\xi\in\Sigma^{n}_{s},\rho\in
H^{d}_{r}}\|A-\xi\otimes\rho\|^{2}_{2}$ (36)
$\displaystyle\quad\leq\min_{W\in\mathrm{DS}^{n},\,\xi\in\Sigma_{n}^{s},\,\lambda\in\Sigma_{r}^{d}}\sum_{i,j=1}^{nd}|A^{\prime}_{i}-(\xi\otimes\lambda)_{j}|^{2}(\operatorname{Id}_{n}\otimes
W)_{i,j},$

where $W$ is optimized over all doubly stochastic matrices
$\mathrm{DS}^{d}\subset\mathbb{C}^{d\times d}$. For $\sigma\in S_{d}$, a
permutation of the the symbols in $[d]$, we denote the corresponding
permutation matrix by $\Pi_{\sigma}:\mathbb{C}^{d}\to\mathbb{C}^{d}$,
$\xi\mapsto\Pi_{\sigma}\xi$ with $(\Pi_{\sigma}\xi)_{i}=\xi_{\sigma(i)}$. By
Birkhoff’s theorem, see e.g., Ref. (Bhatia, 1997, Theorem II.2.3), the set of
extremal points of the convex set of doubly stochastic matrices
$\mathrm{DS}^{d}$ are the permutation matrices
$\Pi_{S_{d}}=\\{\Pi_{\sigma}\mid\sigma\in S_{d}\\}$.

Since the optimum is, hence, attained for a permutation matrix
$W=\Pi_{\sigma}$ and $U_{i,j}=(\Pi_{\sigma})^{1/2}_{i,j}=(\Pi_{\sigma})_{i,j}$
is a unitary matrix, the inequality (36) is saturated. Therefore, we conclude
that

$\displaystyle\min_{\xi\in\Sigma^{n}_{s},\rho\in
H^{y}_{r}}\|A^{\prime}-\xi\otimes\rho\|^{2}_{2}$ (37)
$\displaystyle\quad=\min_{\xi\in\Sigma^{n}_{s},\lambda\in\Sigma^{d}_{r},\sigma\in
S_{d}}\|a^{\prime}-\xi\otimes\Pi_{\sigma}\lambda\|_{2}^{2}$ (38)
$\displaystyle\quad=\min_{\xi\in\Sigma^{n}_{s},\lambda\in\Sigma^{d}_{r}}\|a^{\prime}-\xi\otimes\lambda\|_{2}^{2}$
(39)
$\displaystyle\quad=\min_{\xi\in\Sigma^{n}_{s},\lambda\in\Sigma^{d}_{r}}\|A^{\prime}-\xi\lambda^{T}\|_{2}^{2}.$
(40)

Thus, an algorithm calculating the projection onto $\Omega^{n,d}_{s,r}$ for
the matrix $A$ chosen here solves the sparsePCA problem for $A^{\prime}$. We
conclude that there exists no polynomial time algorithm for the problem. ∎

##  Appendix B Convergence proof

In this section we provide the proof of Theorem 5. We first introduce a bit
more notation. Consider $X\in\Omega^{n,d}_{s,r}$. By definition, it can be
written as $X=\sum_{i=1}^{n}\xi_{i}e_{i}\otimes x_{i}$ with
$\xi\in\Sigma^{n}_{s}$ and $x_{i}\in\mathcal{D}_{r}^{d}$ for all $i$. Let
$Q_{i}$ be the projector onto the range of $x_{i}$. Furthermore, we set
$Q_{i}=0$ for all $i$ not in the support of $\xi$. Slightly overloading our
notation, we define the projection of every ‘block’ onto the range of the
corresponding ‘block’ of $X$ as $\mathcal{P}_{\hat{\Omega}(X)}(Y)\coloneqq
P_{\hat{\Omega}(X)}YP_{\hat{\Omega}(X)}$ with

$P_{\hat{\Omega}(X)}\coloneqq\operatorname{diag}(Q_{1},\ldots,Q_{n}).$ (41)

Note that the projection simultanously projects onto the “block-wise support”
of $X$.

It is common and useful to rewrite the RIP inequalities such as in Definition
4 as an equivalent spectral condition of restrictions of
$\mathcal{A}^{\dagger}\mathcal{A}$.

######  Proposition 10.

Let $X\in\hat{\Omega}^{n,d}_{s,r}$ and $\mathcal{A}:\mathbb{C}^{nd\times
d}\to\mathbb{R}^{m}$ a linear map. Then the following two statements are
equivalent:

  1. a.)

$\left\|\mathcal{P}_{\hat{\Omega}(X)}\circ(\operatorname{Id}-\mathcal{A}^{\dagger}\circ\mathcal{A})\circ\mathcal{P}_{\hat{\Omega}(X)}\right\|_{\infty}\leq\delta$.

  2. b.)

For all $Y\in\operatorname{range}{\mathcal{P}_{\hat{\Omega}(X)}}$ it holds
that

$(1-\delta)\left\|Y\right\|_{F}^{2}\leq\left\|\mathcal{A}(Y)\right\|_{F}^{2}\leq(1+\delta)\left\|Y\right\|_{F}^{2}.$
(42)

###### Proof.

The inequality

$\displaystyle\delta$
$\displaystyle\geq\left\|\mathcal{P}_{\hat{\Omega}(X)}\circ(\operatorname{Id}-\mathcal{A}^{\dagger}\circ\mathcal{A})\circ\mathcal{P}_{\hat{\Omega}(X)}\right\|_{\infty}$
(43)
$\displaystyle=\max_{Y\in\operatorname{range}{\mathcal{P}_{\hat{\Omega}(X)}}}\frac{|\langle
Y,(\operatorname{Id}-\mathcal{A}^{\dagger}\circ\mathcal{A})Y\rangle|}{\left\|Y\right\|_{F}^{2}}$
(44)

holds if and only if for all
$Y\in\operatorname{range}{\mathcal{P}_{\hat{\Omega}(X)}}$

$\delta\left\|Y\right\|_{F}^{2}\geq|\left\|Y\right\|_{F}^{2}-\left\|\mathcal{A}(Y)\right\|_{F}^{2}|.$
(45)

The last bound is equivalent to (42). ∎

We will now prove the recovery guarantee, Theorem 5. The derivation of
recovery guarantees for the IHT algorithm follows largely the same blue print
developed in the original IHT proposal for sparse vectors Blumensath and
Davies (2008), see also Ref. Foucart and Rauhut (2013) for a detailed
description of the proof. Here, we are in addition in the comfortable position
that Ref. Strohmer and Wei (2019) already fleshed out the details of the
recovery proof for an IHT algorithm for de-mixing low-rank matrices. However,
in order to accommodate a non-trivial choice of the step width the proof of
Ref. Strohmer and Wei (2019) yields a slightly weaker result than what can be
shown by a simpler argument for a fixed step width. Thus, we give a slightly
simpler proof that carefully adapts the one given in Ref. Strohmer and Wei
(2019) to account for the additional sparsity constraint and uses a slightly
more concise notation.

###### Proof of Theorem 5.

Let $X\in\hat{\Omega}_{s,r}^{n,d}$ be the matrix to be recovered. Let $X^{l}$
denote the $l$-th iterate of the vector of matrices in the SDT-Algorithm
(Algorithm 2). Since the algorithm always involves a projection step onto
$\hat{\Omega}_{s,r}^{n,d}$ the $l$-th iterate $X^{l}$ is in
$\hat{\Omega}_{s,r}^{n,d}$. Furthermore, we observe that
$X+X^{l}+X^{l+1}\in\hat{\Omega}_{3s,3r}^{n,d}$. For convenience, we denote the
projection onto the (“block-wise”) joint range and support of $X$, $X^{l}$ and
$X^{l+1}$ simply by
$\mathcal{P}^{l}\coloneqq\mathcal{P}_{\hat{\Omega}(X+X^{l}+X^{l+1})}$ and its
orthogonal complement by $\mathcal{P}^{l}_{\perp}$. It is crucial for the
proof to bound norm deviations restricted to the range of $\mathcal{P}^{l}$ as
this eventually allows us to apply a RIP bound.

We want to show the convergence of the iterates of the algorithm $X^{l}$ to
the correct solution $X$. In other words, we want to derive a bound of the
form

$\left\|X^{l+1}-X\right\|_{F}\leq\gamma\left\|X^{l}-X\right\|_{F}$ (46)

with constant $\gamma<1$. Note that by the theorem’s assumption we set the
step width to $\mu^{l}=1$ in the following and omit the tangent space
projection $P_{\mathcal{T}_{X^{l}}}$.

We first derive the following consequence of the thresholding operation: Let
$G^{l}\coloneqq\mathcal{A}^{\dagger}(y-\mathcal{A}(X^{l}))=\mathcal{A}^{\dagger}\circ\mathcal{A}(X-X^{l})$.
By the definition of $X^{l+1}$ as the best approximation to $X^{l}+G^{l}$ in
$\hat{\Omega}_{s,r}^{n,d}$ it holds that

$\begin{split}&\left\|X^{l+1}-\left[X^{l}+G^{l}\right]\right\|_{F}\\\
&\quad\leq\left\|X-\left[X^{l}+G^{l}\right]\right\|_{F}.\end{split}$ (47)

Since the parts of both sides of the inequality that are not in the kernel of
$\mathcal{P}^{l}_{\perp}$ conincides, we get the same inequality also for the
with $\mathcal{P}^{l}$ inserted

$\begin{split}&\left\|X^{l+1}-\left[X^{l}+\mathcal{P}^{l}(G^{l})\right]\right\|_{F}\\\
&\quad\leq\left\|X-\left[X^{l}+\mathcal{P}^{l}(G^{l})\right]\right\|_{F}.\end{split}$
(48)

With the help of this inequality, we can bound

$\begin{split}\left\|X^{l+1}-X\right\|_{F}&\leq\left\|X^{l+1}-\left[X^{l}+\mathcal{P}^{l}(G^{l})\right]\right\|_{F}\\\
&\quad+\left\|X-\left[X^{l}+\mathcal{P}^{l}(G^{l})\right]\right\|_{F}\\\ &\leq
2\left\|X-\left[X^{l}+\mathcal{P}^{l}(G^{l})\right]\right\|_{F}\\\
&=2\left\|\mathcal{M}_{1}(X^{l}-X)\right\|_{F}\\\ &\leq
2\left\|\mathcal{M}_{1}\right\|_{\infty}\left\|X^{l}-X\right\|_{F},\end{split}$
(49)

where in the last step we used the definition of $G^{l}$, the fact that
$\mathcal{P}^{l}$ acts trivially on $X^{l}-X$ and defined
$\mathcal{M}_{1}\coloneqq\mathcal{P}^{l}\circ(\operatorname{Id}-\mathcal{A}^{\dagger}\circ\mathcal{A})\circ\mathcal{P}^{l}$.
To arrive at the theorem’s assertion, we now bound the spectral norm of
$\mathcal{M}_{1}$ using the RIP property of $\mathcal{A}$ and Proposition 10:

$\left\|\mathcal{M}_{1}\right\|_{\infty}=\left\|\mathcal{P}^{l}\circ(\operatorname{Id}-\mathcal{A}^{\dagger}\circ\mathcal{A})\circ\mathcal{P}^{l}\right\|_{\infty}\leq\delta_{3s,3r}$
(50)

since the range of $\mathcal{P}^{l}$ is in $\Omega^{n,d}_{3s,3r}$. Using (50)
in (49) completes the proof. ∎

##  Appendix C RIP guarantee for Hermitian random matrices

In this section we provide the proof of Theorem 6 that establishes the RIP
condition for measurement matrices consisting of Hermitian matrices i.i.d.
drawn from the $\operatorname{GUE}$. Establishing RIP conditions for Gaussian
matrices for a set of structured signals typically proceeds in two steps: One
first derives a strong concentration result for a single signal in the set
using standard concentration of measure. Second, one takes the union bound
over the signal set with the help of an $\epsilon$-covering net construction
to arrive at the uniform statement of RIP. We can readily adapt this strategy
also to $\operatorname{GUE}$.

For the first step, we derive a Gaussian-type concentration result, modifying
a standard line of arguments for our example, see, e.g., Ref. Strohmer and Wei
(2019). The result is summarized as the following lemma:

######  Lemma 11 (Gaussian-type concentration).

Let $X\in\hat{\Omega}^{n,d}_{s,r}$. Let $\\{{A}^{(k)}_{i}\\}_{i=1,k=1}^{n,m}$
be a set of Hermitian matrices drawn i.i.d. from the $\operatorname{GUE}$ and
$\mathcal{A}$ be the measurement operator defined by
$\\{{A}^{(k)}_{i}\\}_{i=1,k=1}^{n,m}$ via Eqs. (9) and (11). Then, for
$0<\delta<1$

$(1-\delta)\left\|X\right\|_{F}^{2}\leq\frac{1}{m}\left\|\mathcal{A}(X)\right\|_{\ell_{2}}^{2}\leq(1+\delta)\left\|X\right\|_{F}^{2}$
(51)

with probability of at least $1-2\mathrm{e}^{-m\delta^{2}/C_{\delta}}$ and
constant $C_{\delta}\geq 40$.

Our proof essentially follows the argument of Ref. Strohmer and Wei (2019) for
Gaussian measurements and then exploits that the Hermitian blocks of the
signal $X\in\hat{\Omega}^{n,d}_{s,r}$ only overlap with the Hermitian part of
the Gaussian measurement matrix.

###### Proof.

Let $X\in\hat{\Omega}^{n,d}_{s,r}$ and denote its $n$ $d\times d$ blocks by
$x_{i}$. Consider a set $\\{B_{i}^{(k)}\in\mathbb{C}^{d\times
d}\\}_{k,i=1}^{m,n}$ of $m\cdot n$ $d\times d$ matrices with entries
independently drawn from the complex-valued normal distribution. Let
$A_{i}^{(k)}\coloneqq P_{\Yboxdim{4pt}\,\yng(2)}{}B_{i}^{(k)}$ be
corresponding matrices drawn from the $\operatorname{GUE}$ and $\mathcal{A}$
the corresponding measurement map. Since all blocks $x_{i}$ are Hermitian, we
have

$\begin{split}\mathcal{A}(X)^{(k)}&=\sum_{i=1}^{n}\langle
A_{i}^{(k)},x_{i}\rangle=\sum_{i=1}^{n}\langle
P_{\Yboxdim{4pt}\,\yng(2)}{}B_{i}^{(k)},x_{i}\rangle\\\
&=\sum_{i=1}^{n}\operatorname{Re}\\{\langle
B_{i}^{(k)},x_{i}\rangle\\}.\end{split}$ (52)

Since all entries of $B_{i}^{(k)}$ are i.i.d. complex normal random variables
and $x_{i}$ is Hermitian, $\operatorname{Re}\\{\langle
B_{i}^{(k)},x_{i}\rangle\\}$ are i.i.d. real random variables from the
distribution $\mathcal{N}(0,\|x_{i}\|^{2}_{F})$ for all $i$ and $k$. We
conclude that all entries $y_{k}=\mathcal{A}(X)^{(k)}$ of $\mathcal{A}(X)$ are
Gaussian distributed with variance
$\sigma^{2}=\sum_{i}\|x_{i}\|^{2}_{F}=\|X\|^{2}_{F}$ and have even moments
$\mathbb{E}[{y_{k}}^{2t}]=2^{-t}t!\binom{2t}{t}\sigma^{2t}$ (Foucart and
Rauhut, 2013, Corollary 7.7). Correspondingly, the squared entries are sub-
exponential random variables with mean $\mathbb{E}[y_{k}^{2}]=\sigma^{2}$. We
denote the associated centered sub-exponential variable as

$z_{k}\coloneqq y_{k}^{2}-\sigma^{2}.$ (53)

The moments of $z_{k}$ are bounded by

$\mathbb{E}[|z_{k}|^{t}]\leq
2^{t}\mathbb{E}[|y_{k}|^{2t}]=t!\binom{2t}{t}\sigma^{2t},$ (54)

where the first inequality follows from the triangle and Jensen’s inequality.
The binomial can be upper bounded using Stirlings formula (Foucart and Rauhut,
2013, (C.13)) by $\binom{2t}{t}=4^{t}r_{t}/\sqrt{\pi t}$ with
$r_{t}\leq\mathrm{e}^{1/(24t)}$. Thus, for $t\geq 2$ we have
$\mathbb{E}[|z_{k}|^{t}]\leq t!R^{t-2}\Sigma^{2}/2$ with $R=4\sigma^{2}$ and
$\Sigma^{2}=\sqrt{2/\pi}\mathrm{e}^{1/48}16\sigma^{4}\leq 0.815\cdot
16\sigma^{4}$. Controlling the moments of $z_{k}$ for $t\geq 2$, we can apply
the Bernstein inequality (Foucart and Rauhut, 2013, Theorem 7.30) and bound
the probability that $\|\mathcal{A}(X)\|_{\ell_{2}}^{2}$ varies by more than
$\Delta>0$ from its expectation value

$\displaystyle\begin{split}&\operatorname{\mathbb{P}}\left[\left|\frac{1}{m}\left\|\mathcal{A}(X)\right\|_{\ell_{2}}^{2}-\|X\|_{F}^{2}\right|\geq\Delta\right]\\\
&\quad\quad=\operatorname{\mathbb{P}}\left[\left|\sum_{k=1}^{m}z_{k}\right|\geq
m\Delta\right]\\\ &\quad\quad\leq
2\exp\left[-\frac{m\Delta^{2}/2}{\Sigma^{2}+R\Delta}\right]\\\ &\quad\quad\leq
2\exp\left[\frac{-m\Delta^{2}}{32\|X\|_{F}^{4}+8\|X\|_{F}^{2}\Delta}\right].\end{split}$
(55)

Let $\Delta=\delta\|X\|_{F}^{2}$ for some $0<\delta<1$. Then we can rewrite
the tail bound (LABEL:eq:bernsteinbound) as

$\displaystyle\mathbb{P}\left[\left|\frac{1}{m}\left\|\mathcal{A}(X)\right\|_{\ell_{2}}^{2}-\|X\|_{F}^{2}\right|\geq\delta\|X\|_{F}^{2}\right]\leq
2\exp\left[-\frac{m\delta^{2}}{C_{\delta}}\right]$ (56)

with a constant $C_{\delta}\geq 40$. Hence, the condition

$\displaystyle(1-\delta)\|X\|_{F}^{2}\leq\frac{1}{m}\left\|\mathcal{A}(X)\right\|_{\ell_{2}}^{2}\leq(1+\delta)\|X\|_{F}^{2}$
(57)

holds with probability at least $1-2\mathrm{e}^{-m\delta^{2}/C_{\delta}}$. ∎

Note that by the homogeneity of the RIP condition it suffices to restrict
ourselves to normalized elements of $\hat{\Omega}^{n,d}_{s,r}$ in the proof of
Theorem 6. In the following, we will therefore focus on the set

$\bar{\Omega}^{n,d}_{s,r}\coloneqq\\{X\in\hat{\Omega}^{n,d}_{s,r}\mid\|X\|_{F}^{2}=1\\}.$
(58)

To take a union bound over the set $\bar{\Omega}^{n,d}_{s,r}$ we need to bound
the size of an $\epsilon$-net that covers the set $\bar{\Omega}^{n,d}_{s,r}$.
An $\epsilon$-net $\mathcal{S}$ covering a set of matrices
$\mathcal{M}\subset\mathbb{C}^{nd\times d}$ is a finite subset of
$\mathcal{M}$ such that for all $X\in\mathcal{M}$ there exists
$\bar{X}\in\mathcal{S}$ such that $\|X-\bar{X}\|_{F}\leq\epsilon$. Our
construction generalizes the construction of Ref. Strohmer and Wei (2019).
Therein, a covering net for the set of normalized block-wise low-rank matrices
$\bar{\Omega}^{n,d}_{n,r}$ was derived. We summarize the statement given in
Ref. Strohmer and Wei (2019) in the following lemma without giving a proof.

######  Lemma 12 (Strohmer and Wei (2019)).

For $\bar{\Omega}^{n,d}_{n,r}$ there exists an $\epsilon$-covering net
$\mathcal{S}^{n,d}_{r}$ with cardinality bounded by $(9/\epsilon)^{(2d+1)nr}$.

The proof of Lemma 12 basically lifts the result of an $\epsilon$-net for low-
rank matrices of Ref. Candes and Plan (2011) to the set
$\bar{\Omega}^{n,d}_{n,r}$ using the triangle inequality.

We can combine multiple $\epsilon$-nets for $\bar{\Omega}^{s,d}_{s,r}$ to
construct an $\epsilon$-covering net for the set $\bar{\Omega}^{n,d}_{s,r}$ of
block-sparse matrix vectors with low-rank blocks. The bound on the cardinality
of the resulting $\epsilon$-covering net is given in the following lemma:

######  Lemma 13 (Bound on the cardinality of a covering net).

For $\bar{\Omega}^{n,d}_{s,r}$ there exists an $\epsilon$-covering net
$\mathcal{S}^{n,d}_{s,r}$ of cardinality bounded by
$\binom{n}{s}(9/\epsilon)^{(2d+1)sr}$. Furthermore, for each
$X=[X_{1},\ldots,X_{n}]\in\bar{\Omega}^{n,d}_{s,r}$ there exists
$\bar{X}=[\bar{X}_{1},\ldots,\bar{X}_{n}]\in\mathcal{S}^{n,d}_{s,r}$ such that
$\|X-\bar{X}\|_{F}\leq\epsilon$ and $\|\bar{X}_{k}\|_{F}=0$ for all $k$ for
which $\|X_{k}\|_{F}=0$.

###### Proof.

Let $\Gamma\subset[n]$ with $|\Gamma|\leq s$, i.e., the indices of the support
of an $s$-sparse vector. The set

$\bar{\Omega}^{\Gamma}_{r}\coloneqq\left\\{\sum_{i\in\Gamma}\xi_{i}e_{i}\otimes
x_{i}\ \middle|\ \xi_{i}\in\mathbb{R},\ x_{i}\in\mathcal{D}^{d}_{r}\ \forall
i\right\\}\subset\bar{\Omega}^{n,d}_{s,r}$ (59)

shall consist of all elements of $\bar{\Omega}^{n,d}_{s,r}$ which have non-
vanishing blocks only supported on $\Gamma$. To each element of
$\bar{\Omega}^{\Gamma}_{r}$, we can associate an element of
$\bar{\Omega}^{s,d}_{s,r}$ by omitting the vanishing blocks in the matrix
vector and vice versa. By virtue of Lemma 12 we thus know that
$\bar{\Omega}^{\Gamma}_{r}$ has a covering net $\mathcal{S}^{\Gamma}_{r}$ of
cardinality bounded by $(9/\epsilon)^{(2d+1)sr}$.

We can decompose the entire set $\bar{\Omega}^{n,d}_{s,r}$ as

$\bar{\Omega}^{n,d}_{s,r}=\bigcup_{\Gamma\subset[n],|\Gamma|\leq
s}\bar{\Omega}^{\Gamma}_{r},$ (60)

and thus, the set

$\mathcal{S}^{n,d}_{s,r}=\bigcup_{\Gamma\subset[n],|\Gamma|\leq
s}\mathcal{S}^{\Gamma}_{r}$ (61)

is an $\epsilon$-covering net for $\bar{\Omega}^{n,d}_{s,r}$. The union is
taken over $\binom{n}{s}$ different sets. Thus, the cardinality of
$\mathcal{S}^{n,d}_{s,r}$ is upper bounded by
$\binom{n}{s}(9/\epsilon)^{(2d+1)sr}$. The second statement follows by
construction. ∎

We are now in the position to prove Theorem 6.

###### Proof of Theorem 6.

The proof proceeds in two steps. First, we prove the RIP for elements of the
$\epsilon$-covering net $\mathcal{S}^{n,d}_{s,r}$ of
$\bar{\Omega}^{n,d}_{s,r}$. To do so, we combine the concentration result of
Lemma 11 and the union bound of Lemma 13 to establish uniform concentration.
In a second step, following Ref. Strohmer and Wei (2019), we then use the
definition of an $\epsilon$-covering net to show that for elements
$X\in\bar{\Omega}^{n,d}_{s,r}$ that are close enough to an element of the net,
the RIP condition still holds.

_Step 1:_ Taking the union bound over the $\epsilon$-net
$\mathcal{S}^{n,d}_{s,r}$ constructed in Lemma 13 and using the result of
Lemma 11 in the form of (56) with constant $C_{\delta}\geq 40$ we get

$\begin{split}&\operatorname{\mathbb{P}}\left(\max_{\begin{subarray}{c}X\in\mathcal{S}^{n,d}_{s,r}\end{subarray}}\left|\frac{1}{m}\|\mathcal{A}(X)\|_{\ell_{2}}^{2}-\|X\|_{F}^{2}\right|\quad\geq\delta/2\right)\\\
&\quad\leq
2|\mathcal{S}^{n,d}_{s,r}|\mathrm{e}^{-m\delta^{2}/(4C_{\delta})}\\\
&\quad\leq
2\binom{n}{s}\left(\frac{9}{\epsilon}\right)^{(2d+1)sr}e^{-m\delta^{2}/(4C_{\delta})}.\end{split}$
(62)

The aim is to find a lower bound for the number of measurements $m$ for which
the probability (62) small. To this end, we rewrite

$\displaystyle\begin{split}&2\binom{n}{s}\left(\frac{9}{\epsilon}\right)^{(2d+1)sr}e^{-m\delta^{2}/(4C_{\delta})}\\\
&\quad\leq
2\exp\left[s\ln\frac{\mathrm{e}n}{s}+(2d+1)sr\ln\frac{9}{\epsilon}-\frac{m\delta^{2}}{4C_{\delta}}\right]\\\
&\quad\leq\tau,\end{split}$ (63)

using $\binom{n}{s}\leq\left(\frac{en}{s}\right)^{s}$ (Foucart and Rauhut,
2013, Lemma C.5). The latter inequality becomes true under the condition that

$m\geq\frac{4C_{\delta}}{\delta^{2}}\left[s\ln\frac{\mathrm{e}n}{s}+(2d+1)sr\ln\frac{9}{\epsilon}+\ln\frac{2}{\tau}\right].$
(64)

Assuming that (64) holds, we have established the RIP condition for the
$\epsilon$-net $\mathcal{S}^{n,d}_{s,r}$, i.e., for all vectors
$\overline{X}\in\mathcal{S}^{n,d}_{s,r}$ it holds that

$(1-\delta/2)\|\overline{X}\|_{F}^{2}\leq\|\mathcal{A}(\overline{X})\|_{\ell_{2}}^{2}\leq(1+\delta/2)\|\overline{X}\|_{F}^{2}$
(65)

with probability at least $1-\tau$.

_Step 2:_ Let us now transfer the RIP of $\mathcal{S}^{n,d}_{s,r}$ to the
entire set $\bar{\Omega}^{n,d}_{s,r}$ while keeping the error under control.
To this end, we choose the net parameter $\epsilon$ as
$\frac{\delta}{4\sqrt{2}}$. By definition of an $\epsilon$-net, for elements
$X\in\bar{\Omega}^{n,d}_{s,r}$, there exists an element
$\overline{X}\in\mathcal{S}^{n,d}_{s,r}$ such that

$\|X-\overline{X}\|_{F}\leq\frac{\delta}{4\sqrt{2}}.$ (66)

To prove the RIP for the set $\bar{\Omega}^{n,d}_{s,r}$ we need to bound
$\|\mathcal{A}(X)\|_{F}$ from above and below.

We start with the upper bound, making use of Eq. (65):

$\begin{split}\|\mathcal{A}(X)\|_{\ell_{2}}&\leq\|\mathcal{A}(\overline{X})\|_{\ell_{2}}+\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}\\\
&\leq
1+\frac{\delta}{2}+\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}.\end{split}$
(67)

Now $\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}$ has to be bounded from above.
We use that by the second statement of Lemma 13 the block supports of $X$ and
$\overline{X}$ coincide. Therefore, $X-\overline{X}$ has also $s$ non-
vanishing blocks that have rank of atmost $2r$. We can, thus, decompose
$X-\overline{X}=B+C$ in terms of orthogonal matrices
$B,C\in\hat{\Omega}^{n,d}_{s,r}$ that obey $\langle B,C\rangle=0$. In
particular, $B$ and $C$ have the same block support as $X$. Let us define

$\kappa_{s,r}\coloneqq\sup_{X\in\bar{\Omega}^{n,d}_{s,r}}\|\mathcal{A}(X)\|_{\ell_{2}}.$
(68)

Then we get using homogeniety

$\begin{split}\|\mathcal{A}&(X-\overline{X})\|_{\ell_{2}}\leq\|\mathcal{A}(B)\|_{\ell_{2}}+\|\mathcal{A}(C)\|_{\ell_{2}}\\\
&\leq\kappa_{s,r}(\|B\|_{F}+\|C\|_{F})\leq\sqrt{2}\,\kappa_{s,r}\sqrt{\|B\|_{F}^{2}+\|C\|_{F}^{2}}\\\
&=\sqrt{2}\,\kappa_{s,r}\|X-\bar{X}\|_{F},\end{split}$ (69)

where the last step makes use of the orthogonality of $B$ and $C$. Together
with (66) it follows that

$\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}\leq\frac{\delta\cdot\kappa_{s,r}}{4}.$
(70)

It remains to derive an upper bound for $\kappa_{s,r}$. To this end, we use
that, by definition, $\kappa_{s,r}$ is the best upper bound of the left hand
side of (67). Inserting (70) into the right hand side of (67), we find the
condition

$\kappa_{s,r}\leq 1+\frac{\delta}{2}+\frac{\delta\cdot\kappa_{s,r}}{4}.$ (71)

Solving for $\kappa_{s,r}$, Eq. (71) implies for $0<\delta<1$

$\begin{split}\kappa_{s,r}\leq\frac{1+\delta/2}{1-\delta/4}\leq
1+\delta.\end{split}$ (72)

Altogether, this yields the desired upper bound

$\|\mathcal{A}(X)\|_{\ell_{2}}\leq
1+\frac{3}{4}\delta+\frac{\delta^{2}}{4}\leq 1+\delta,$ (73)

for $\delta<1$. The lower bound is analogously obtained by combining the
inequality

$\displaystyle\|\mathcal{A}(X)\|_{\ell_{2}}$
$\displaystyle\geq\|\mathcal{A}(\overline{X})\|_{\ell_{2}}-\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}$
(74) $\displaystyle\geq 1-\delta/2-\|\mathcal{A}(X-\overline{X})\|_{\ell_{2}}$
(75)

with (70) (72) to arrive at

$\|\mathcal{A}(X)\|_{\ell_{2}}\geq 1-\delta/2-\delta(1+\delta)/4\geq
1-\delta.$ (76)

With the choice of $\epsilon$, we can rewrite the condition (64) on $m$ as

$m\geq\frac{C}{\delta^{2}}\left[s\ln\frac{\mathrm{e}n}{s}+(2d+1)sr\ln\frac{c}{\delta}+\ln\frac{2}{\tau}\right]$
(77)

with constants $C\geq 4C_{\delta}\geq 160$ and $c\geq 36\sqrt{2}\geq 51$. This
completes the proof.
∎

## Acknowledgements

We thank David Gross, Steven T. Flammia, Christian Krumnow, Robin Harper, Yi-
Kai Liu, and Carlos A. Riofrio for inspiring discussions and helpful comments.
Furthermore, we are grateful to Alireza Seif and Nobert Linke for valuable
comments on realistic error models and Peter Jung for making us aware of Ref.
Strohmer and Wei (2019). We are grateful to Susane Calegari for kindly
providing drawings used in Figure 1. This work has been supported by the DFG,
specifically SPP1798 CoSIP, but also EI 519/9-1, EI 519/7-1, and CRC 183. It
has also received funding from the Templeton Foundation and from the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No 817482 (PASQuanS).

## References

  * Preskill (2018) J. Preskill,  _Quantum computing in the NISQ era and beyond_ , Quantum 2, 79 (2018).
  * Acin  _et al._ (2017) A. Acin, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm,  _The European quantum technologies roadmap_ , (2017), arXiv:1712.03773 [quant-ph].
  * Eisert  _et al._ (2019) J. Eisert, D. Hangleiter, N. Walk, I. Roth, D. Markham, R. Parekh, U. Chabaud, and E. Kashefi,  _Quantum certification and benchmarking_ , ArXiv e-prints (2019), arXiv:1910.06343 [quant-ph].
  * Brańczyk  _et al._ (2012) A. M. Brańczyk, D. H. Mahler, L. A. Rozema, A. Darabi, A. M. Steinberg, and D. F. V. James,  _Self-calibrating quantum state tomography_ , New J. Phys. 14, 085003 (2012).
  * Gross  _et al._ (2010) D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert,  _Quantum state tomography via compressed sensing_ , Phys. Rev. Lett. 105, 150401 (2010).
  * Flammia  _et al._ (2012) S. T. Flammia, D. Gross, Y.-K. Liu, and J. Eisert,  _Quantum tomography via compressed sensing: error bounds, sample complexity and efficient estimators_ , New J. Phys. 14, 095022 (2012), arXiv:1205.2300 [quant-ph].
  * Kalev  _et al._ (2015) A. Kalev, R. L. Kosut, and I. H. Deutsch,  _Quantum tomography protocols with positivity are compressed sensing protocols_ , npj Quantum Information 1, 15018 (2015), arXiv:1502.00536.
  * Riofrio  _et al._ (2017) C. A. Riofrio, D. Gross, S. T. Flammia, T. Monz, D. Nigg, R. Blatt, and J. Eisert,  _Experimental quantum compressed sensing for a seven-qubit system_ , Nature Comm. 8, 15305 (2017).
  * Steffens  _et al._ (2017) A. Steffens, C. A. Riofrío, W. McCutcheon, I. Roth, B. A. Bell, A. McMillan, M. S. Tame, J. G. Rarity, and J. Eisert,  _Experimentally exploring compressed sensing quantum tomography_ , Quant. Sc. Tech. 2, 025005 (2017).
  * Baraniuk  _et al._ (2010) R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde,  _Model-based compressive sensing_ , IEEE Trans. Inf. Th. 56, 1982 (2010).
  * Foucart and Rauhut (2013) S. Foucart and H. Rauhut,  _A mathematical introduction to compressive sensing_ (Springer, 2013).
  * Blumensath and Davies (2008) T. Blumensath and M. E. Davies,  _Iterative thresholding for sparse approximations_ , J. Four. An. App. 14, 629 (2008).
  * Strohmer and Wei (2019) T. Strohmer and K. Wei,  _Painless breakups-efficient demixing of low rank matrices_ , J. Four. Ana. App. 25, 1 (2019).
  * Roth  _et al._ (2016a) I. Roth, M. Kliesch, A. Flinth, G. Wunder, and J. Eisert,  _Reliable recovery of hierarchically sparse signals for Gaussian and Kronecker product measurements_ , Arxiv e-prints (2016a), accepted in IEEE Trans. Signal Process., arXiv:1612.07806.
  * Pawłowski and Brunner (2011) M. Pawłowski and N. Brunner,  _Semi-device-independent security of one-way quantum key distribution_ , Phys. Rev. A 84, 010302 (2011).
  * Liang  _et al._ (2011) Y.-C. Liang, T. Vértesi, and N. Brunner,  _Semi-device-independent bounds on entanglement_ , Phys. Rev. A 83, 022108 (2011).
  * Li  _et al._ (2011) H.-W. Li, Z.-Q. Yin, Y.-C. Wu, X.-B. Zou, S. Wang, W. Chen, G.-C. Guo, and Z.-F. Han,  _Semi-device-independent random-number expansion without entanglement_ , Phys. Rev. A 84, 034301 (2011).
  * Li  _et al._ (2012) H.-W. Li, M. Pawłowski, Z.-Q. Yin, G.-C. Guo, and Z.-F. Han,  _Semi-device-independent randomness certification using $n\rightarrow 1$ quantum random access codes_, Phys. Rev. A 85, 052308 (2012).
  * Gallego  _et al._ (2010) R. Gallego, N. Brunner, C. Hadley, and A. Acín,  _Device-independent tests of classical and quantum dimensions_ , Phys. Rev. Lett. 105, 230501 (2010).
  * Mogilevtsev (2010) D. Mogilevtsev,  _Calibration of single-photon detectors using quantum statistics_ , Phys. Rev. A 82, 021807 (2010).
  * Mogilevtsev  _et al._ (2009) D. Mogilevtsev, J. Řeháček, and Z. Hradil,  _Relative tomography of an unknown quantum state_ , Phys. Rev. A 79, 020101 (2009).
  * Mogilevtsev  _et al._ (2012) D. Mogilevtsev, J. Řeháček, and Z. Hradil,  _Self-calibration for self-consistent tomography_ , New J. Phys. 14, 095001 (2012).
  * Sim  _et al._ (2019) J. Y. Sim, J. Shang, H. K. Ng, and B.-G. Englert,  _Proper error bars for self-calibrating quantum tomography_ , Arxiv e-prints (2019), arXiv:1904.11202.
  * Stark (2012) C. Stark,  _Simultaneous estimation of dimension, states and measurements: Computation of representative density matrices and POVMs_ , Arxiv e-prints (2012), arXiv:1210.1105.
  * Stark (2014) C. Stark,  _Self-consistent tomography of the state-measurement Gram matrix_ , Phys. Rev. A 89, 052109 (2014), arXiv:1209.5737.
  * Řeháček  _et al._ (2010) J. Řeháček, D. Mogilevtsev, and Z. Hradil,  _Operational tomography: Fitting of data patterns_ , Phys. Rev. Lett. 105, 010402 (2010).
  * Motka  _et al._ (2014) L. Motka, B. Stoklasa, J. Rehacek, Z. Hradil, V. Karasek, D. Mogilevtsev, G. Harder, C. Silberhorn, and L. L. Sánchez-Soto,  _Efficient algorithm for optimizing data-pattern tomography_ , Phys. Rev. A 89, 054102 (2014).
  * Ferrie (2014) C. Ferrie,  _Quantum model averaging_ , New J. Phys. 16, 093035 (2014).
  * Merkel  _et al._ (2013) S. T. Merkel, J. M. Gambetta, J. A. Smolin, S. Poletto, A. D. Córcoles, B. R. Johnson, C. A. Ryan, and M. Steffen,  _Self-consistent quantum process tomography_ , Phys. Rev. A 87, 062119 (2013), arXiv:1211.0322 [quant-ph].
  * Blume-Kohout  _et al._ (2013) R. Blume-Kohout, J. King Gamble, E. Nielsen, J. Mizrahi, J. D. Sterk, and P. Maunz,  _Robust, self-consistent, closed-form tomography of quantum logic gates on a trapped ion qubit_ , ArXiv e-prints (2013), arXiv:1310.4492 [quant-ph].
  * Greenbaum (2015) D. Greenbaum,  _Introduction to quantum gate set tomography_ , arXiv:1509.02921 [quant-ph].
  * Blume-Kohout  _et al._ (2017) R. Blume-Kohout, J. K. Gamble, E. Nielsen, K. Rudinger, J. Mizrahi, K. Fortier, and P. Maunz,  _Demonstration of qubit operations below a rigorous fault tolerance threshold with gate set tomography_ , Nature Comm. 8, 14485 (2017), arXiv:1605.07674 [quant-ph].
  * Cerfontaine  _et al._ (2019) P. Cerfontaine, R. Otten, and H. Bluhm,  _Self-consistent calibration of quantum gate sets_ , Arxiv e-prints (2019), arXiv:1906.00950.
  * Gross (2011) D. Gross,  _Recovering low-rank matrices from few coefficients in any basis_ , IEEE Trans. Inf. Th. 57, 1548 (2011), arXiv:0910.1879 [cs.IT].
  * Liu (2011) Y.-K. Liu,  _Universal low-rank matrix recovery from Pauli measurements_ , Adv. Neural Inf. Process. Syst. 24, 1638 (2011), 1103.2816.
  * Kueng (2015) R. Kueng,  __Low rank matrix recovery from few orthonormal basis measurements_ , _in  _Sampling Theory and Applications (SampTA), 2015 International Conference on_ (2015) pp. 402–406.
  * Kabanava  _et al._ (2016) M. Kabanava, R. Kueng, H. Rauhut, and U. Terstiege,  _Stable low-rank matrix recovery via null space properties_ , Inf. Inf. 5, 405 (2016).
  * Shabani  _et al._ (2011) A. Shabani, R. L. Kosut, M. Mohseni, H. Rabitz, M. A. Broome, M. P. Almeida, A. Fedrizzi, and A. G. White,  _Efficient measurement of quantum dynamics via compressive sensing_ , Phys. Rev. Lett. 106, 100401 (2011).
  * Kimmel and Liu (2017) S. Kimmel and Y. K. Liu,  __Phase retrieval using unitary 2-designs_ , _in  _2017 International Conference on Sampling Theory and Applications (SampTA)_ (2017) pp. 345–349, arXiv:1510.08887.
  * Roth  _et al._ (2018) I. Roth, R. Kueng, S. Kimmel, Y.-K. Liu, D. Gross, J. Eisert, and M. Kliesch,  _Recovering quantum gates from few average gate fidelities_ , Phys. Rev. Lett. 121, 170502 (2018), arXiv:1803.00572 [quant-ph].
  * Wunder  _et al._ (2015) G. Wunder, H. Boche, T. Strohmer, and P. Jung,  _Sparse signal processing concepts for efficient 5G system design_ , IEEE Access 3, 195 (2015).
  * Roth  _et al._ (2016b) I. Roth, M. Kliesch, G. Wunder, and J. Eisert,  __Reliable recovery of hierarchically sparse signals_ , _in  _Proceedings of the third “international traveling workshop on interactions between sparse models and technology” (iTWIST’16)_ (2016) arXiv:1609.04167 [cs.NA].
  * Note (1) Note that one might actually implement projective measurements in the multi-qubit Pauli basis as done, e.g., in Refs. Riofrio  _et al._ (2017); Steffens  _et al._ (2017). While such projective measurements contain more information than the Pauli correlation measurement, we restrict ourselves to Pauli expectation values both for the sake of simplicity and to remain in a setting for which theoretical guarantees can be proven Gross  _et al._ (2010); Liu (2011).
  * Ahmed  _et al._ (2014) A. Ahmed, B. Recht, and J. Romberg,  _Blind deconvolution using convex programming_ , IEEE Trans. Inf. Th. 60, 1711 (2014).
  * Oymak  _et al._ (2015) S. Oymak, A. Jalali, M. Fazel, Y. C. Eldar, and B. Hassibi,  _Simultaneously structured models with application to sparse and low-rank matrices_ , IEEE Trans. Inf. Th. 61, 2886 (2015).
  * Needell and Tropp (2008) D. Needell and J. A. Tropp,  _CoSaMP: Iterative signal recovery from incomplete and inaccurate samples_ , Appl. Comp. Harm. An. (2008).
  * Foucart (2011) S. Foucart,  _Hard thresholding pursuit: An algorithm for compressive sensing_ , SIAM J. Num. An. 49, 2543 (2011).
  * Magdon-Ismail (2017) M. Magdon-Ismail,  _NP-hardness and inapproximability of sparse PCA_ , Inf. Proc. Lett. 126, 35 (2017).
  * Berthet and Rigollet (2013a) Q. Berthet and P. Rigollet,  __Complexity theoretic lower bounds for sparse principal component detection_ , _in  _Conference on Learning Theory_ (2013) pp. 1046–1066.
  * Berthet and Rigollet (2013b) Q. Berthet and P. Rigollet,  _Optimal detection of sparse principal components in high dimension_ , Ann. Statist. 41, 1780 (2013b).
  * Brennan and Bresler (2019) M. Brennan and G. Bresler,  _Optimal average-case reductions to sparse PCA: From weak assumptions to strong hardness_ , Arxiv e-prints (2019), arXiv:1902.07380.
  * Chan  _et al._ (2016) S. O. Chan, D. Papailliopoulos, and A. Rubinstein,  __On the approximability of sparse PCA_ , _in  _PMLR_, Vol. 49 (2016) pp. 623–646, arXiv:1507.05950.
  * Wunder  _et al._ (2018a) G. Wunder, I. Roth, R. Fritschek, B. Groß, and J. Eisert,  __Secure massive IoT using hierarchical fast blind deconvolution_ , _in  _2018 IEEE Wireless Communications and Networking Conference Workshops, WCNC 2018 Workshops, Barcelona, Spain, April 15-18, 2018_ (2018) pp. 119–124, arXiv:1801.09628.
  * Foucart  _et al._ (2019) S. Foucart, R. Gribonval, L. Jacques, and H. Rauhut,  _Jointly low-rank and bisparse recovery: Questions and partial answers_ , arXiv e-prints (2019), arXiv:1902.04731 [math.NA].
  * Sprechmann  _et al._ (2010) P. Sprechmann, I. Ramirez, G. Sapiro, and Y. Eldar,  __Collaborative hierarchical sparse modeling_ , _in  _2010 44th Annual Conference on Information Sciences and Systems (CISS)_ (2010) pp. 1–6.
  * Friedman  _et al._ (2010) J. Friedman, T. Hastie, and R. Tibshirani,  _A note on the group lasso and a sparse group lasso_ , Arxiv e-prints (2010), arXiv:1001.0736.
  * Sprechmann  _et al._ (2011) P. Sprechmann, I. Ramirez, G. Sapiro, and Y. C. Eldar,  _C-HiLasso: A collaborative hierarchical sparse modeling framework_ , IEEE Trans. Sig. Proc. 59, 4183 (2011).
  * Simon  _et al._ (2013) N. Simon, J. Friedman, T. Hastie, and R. Tibshirani,  _A sparse-group Lasso_ , J. Comp. Graph. Stat. 22, 231 (2013).
  * Roth  _et al._ (2018) I. Roth, A. Flinth, R. Kueng, J. Eisert, and G. Wunder,  __Hierarchical restricted isometry property for Kronecker product measurements_ , _in  _2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ (2018) pp. 632–638.
  * Groß  _et al._ (2020) B. Groß, A. Flinth, I. Roth, J. Eisert, and G. Wunder,  _Hierarchical sparse recovery from hierarchically structured measurements_ , ArXiv e-prints (2020), arXiv:2005.10379.
  * Wunder  _et al._ (2017) G. Wunder, I. Roth, R. Fritschek, and J. Eisert,  __HiHTP: A custom-tailored hierarchical sparse detector for massive MTC_ , _in  _2017 51st Asilomar Conference on Signals, Systems, and Computers_ (2017) pp. 1929–1934.
  * Wunder  _et al._ (2018b) G. Wunder, I. Roth, R. Fritschek, and J. Eisert,  _Performance of hierarchical sparse detectors for massive MTC_ , Arxiv e-print (2018b), arXiv:1806.02754.
  * Wunder  _et al._ (2018c) G. Wunder, I. Roth, M. Barzegar, A. Flinth, S. Haghighatshoar, G. Caire, and G. Kutyniok,  __Hierarchical sparse channel estimation for massive MIMO_ , _in  _WSA 2018; 22nd International ITG Workshop on Smart Antennas_ (VDE, 2018) pp. 1–8.
  * Wunder  _et al._ (2019) G. Wunder, S. Stefanatos, A. Flinth, I. Roth, and G. Caire,  _Low-overhead hierarchically-sparse channel estimation for multiuser wideband massive MIMO_ , IEEE Transactions on Wireless Communications 18, 2186 (2019).
  * Golub and van Loan (1989) G. H. Golub and C. F. van Loan,  _Matrix computations_ (The Johns Hopkins University Press, 1989).
  * Halko  _et al._ (2011) N. Halko, P.-G. Martinsson, and J. A. Tropp,  _Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions_ , SIAM Rev. 53, 217 (2011).
  * Hoare (1961) C. A. R. Hoare,  _Algorithm 65: Find_ , Commun. ACM 4, 321 (1961).
  * Wei  _et al._ (2016) K. Wei, J.-F. Cai, T. F. Chan, and S. Leung,  _Guarantees of Riemannian optimization for low rank matrix recovery_ , SIAM J. Mat. An. App. 37, 1198 (2016).
  * Absil  _et al._ (2009) P.-A. Absil, R. Mahony, and R. Sepulchre,  _Optimization algorithms on matrix manifolds_ (Princeton University Press, 2009).
  * Vandereycken (2013) B. Vandereycken,  _Low-rank matrix completion by Riemannian optimization_ , SIAM J. Opt. 23, 1214 (2013).
  * Blumensath and Davies (2009) T. Blumensath and M. E. Davies,  _Sampling theorems for signals from the union of finite-dimensional linear subspaces_ , IEEE Trans. Inf. Theory 55, 1872 (2009).
  * Bandeira  _et al._ (2013) A. S. Bandeira, E. Dobriban, D. G. Mixon, and W. F. Sawin,  _Certifying the restricted isometry property is hard_ , IEEE Trans. Inf. Th. 59, 3448 (2013).
  * Wilkens  _et al._ (2020) J. Wilkens, D. Hangleiter, and I. Roth, (2020), Gitlab repository at https://gitlab.com/wilkensJ/blind-quantum-tomography.
  * Bhatia (1997) R. Bhatia,  _Matrix analysis_, Graduate texts in mathematics No. 169 (Springer, New York, 1997).