# A universal scheme for robust self-testing in the prepare-and-measure
scenario

Nikolai Miklin<EMAIL_ADDRESS>Institute of Theoretical Physics and
Astrophysics, National Quantum Information Center, Faculty of Mathematics,
Physics and Informatics, University of Gdansk, 80-306 Gdańsk, Poland
International Centre for Theory of Quantum Technologies (ICTQT), University of
Gdansk, 80-308 Gdańsk, Poland  Michał Oszmaniec<EMAIL_ADDRESS>Center
for Theoretical Physics, Polish Academy of Sciences, Al. Lotników 32/46,
02-668 Warszawa, Poland

###### Abstract

We consider the problem of certification of arbitrary ensembles of pure states
and projective measurements solely from the experimental statistics in the
prepare-and-measure scenario assuming the upper bound on the dimension of the
Hilbert space. To this aim we propose a universal and intuitive scheme based
on establishing perfect correlations between target states and suitably-chosen
projective measurements. The method works in all finite dimensions and allows
for robust certification of the overlaps between arbitrary preparation states
and between the corresponding measurement operators. Finally, we prove that
for qubits our technique can be used to robustly self-test arbitrary
configurations of pure quantum states and projective measurements. These
results pave the way towards practical application of the prepare-and-measure
paradigm to certification of quantum devices.

##  1 Introduction

Quantum devices are becoming more and more complex and the possibilities of
their precise control and manipulation keep increasing. Recently reported
demonstration of quantum computational advantage by Google [1] is only an
intermediate milestone and quantum technologies have a potential real-life
applications in fields such as quantum sensing [2], simulation of quantum
systems [3], efficient computation [4] and machine learning [5, 6].

With the increasing complexity of quantum systems, there is a growing need for
certification and verification of their performance. This task is usually
realized via the combination of quantum tomography and various benchmarking
schemes (see [7] for a recent review). However, these methods, despite being
powerful and universally applicable, depend on the assumptions about the inner
workings of quantum systems, such as perfect measurements or uncorrelated and
independent errors. In contrast to these approaches _self-testing_ is a method
which aims at proving the uniqueness of the implemented states or measurements
based solely on the observed statistics and under minimal physical
assumptions.

The paradigm of self-testing was first introduced in the context of quantum
cryptography [8], with the aim to obtain trust in cryptographic devices (see
[9] for a recent review). Initially, it was applied to correlations observed
in the Bell scenario [10] (see e.g., [8, 11, 12, 13]). The most know result in
this area is certification of the singlet state in the case of maximal
violation of the Clauser-Horne-Shimony-Holt Bell inequality [14].

With the growing number of results on self-testing, more and more attention is
being drawn to prepare-and-measure scenarios, that are more experimentally
appealing as compared to the one of Bell (see e.g., [15, 16, 17, 18]).
Therein, one does not need to ensure space-like separation of the measurement
events by two parties; in contrast, one party, Alice, communicates some of her
states to Bob, who measures them in some basis of his choice. In order to get
meaningful certification results further assumptions are needed. In the most
commonly studied _semi-device-independent_ (SDI) scenario [19], one assumes
that the dimension of the quantum system used for transmitting information is
bounded from above. There exist, however, alternative approaches based on
other constraints like minimal overlap [16], mean energy constraint [20] or
entropy constraint [21].

Let us briefly recap on what has been done so far in the area of SDI self-
testing. First, self-testing results were proven for mutually unbiased bases
(MUBs) in dimension $2$, for both state ensembles and measurements [22]. This
was further generalised to SDI certification of mutual unbiasedness of pairs
of bases in an arbitrary dimension in Ref. [23]. Methods for self-testing of
extremal qubit positive operator-valued measures (POVMs) were proposed in [24,
25] and further extended to symmetric-informationally complete (SIC) POVMs
[26]. Importantly, all of the above results either rely on numerical
approaches, for general state preparations and POVMs, or work only for special
scenarios that exhibit many symmetries.

In this work we propose a simple analytical method allowing to certify
overlaps between preparations of arbitrary pure states and arbitrary
projective measurements in _qudit_ systems. The scheme relies on establishing
perfect correlations between preparation states and outcomes of suitably-
chosen projective measurements. The method is universally applicable and
robust to experimental noise. We prove that for qubits our SDI certification
method can be used to obtain a robust self-testing result for arbitrary
preparations of pure qubit states and corresponding projective measurements.
While for higher dimensions we do not show self-testing, our scheme allows for
SDI certification of numerous inherently quantum properties. The examples
include, but are not limited to: _arbitrary_ overlaps between any number of
measurement bases, MUB conditions, information-completeness of measurements,
and SIC relations among measurement effects.

We believe that our findings greatly extend the applicability of the paradigm
of self-testing in the SDI setting. They will likely have application for
certification of near-term quantum computers [7], especially since our scheme
is not sensitive to state-preparation and measurement errors, which is one of
the major problems in certification of quantum devices [27]. We also expect
possible cryptographic applications as our setup is very similar to the one of
textbook quantum key distribution schemes [28, 29]. The perfect correlations
can be utilized for generation of the secret key, while the rest can be used
to estimate the security. Thus, our methods can be directly applied for
certification of quantum devices implementing protocols such as BB84 [28],
which is normally achieved by introducing additional preparation states or
measurement bases [30].

Before we proceed, let us first fix some of the notations used in this paper.
Let $X$ be a linear operator acting on a finite-dimensional Hilbert space
$\mathcal{H}$. Throughout the paper we use $\left|\left|\hskip 1.0ptX\hskip
1.0pt\right|\right|$, $\left|\left|\hskip 1.0ptX\hskip
1.0pt\right|\right|_{F}$ to denote operator norm and Frobenius norm of $X$. We
will also use $|\mathbf{n}|$ to denote the Euclidean norm of
$\mathbf{n}\in\mathbb{R}^{3}$, and $[n]$ to denote an $n$-element set
$\\{1,2,\dots,n\\}$.

##  2 Description of the scenario

We consider the prepare-and-measure scenario in which in each run of the
experiment Alice prepares a quantum $d$-level system in one of the states from
a finite set of preparations $\varrho_{a}^{x}$, for which we use two indexes
$x\in[n]$ and $a\in[d]$. Subsequently, Bob performs a measurement on this
state with a finite choice of measurement settings $y\in[n]$ having the
possible outcomes $b\in[d]$. We assume that measurement process performed by
Bob is described by quantum mechanics and furthermore that the parties do not
communicate in any other way and do not have access to any entangled states or
shared randomness [31] (the role of this assumption is discussed in detail
later in the text). This implies that the observed statistics $p(b|a,x,y)$ are
given by the Born rule i.e.,
$p(b|a,x,y)={\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})$, where
$\mathbf{M}^{y}=(M^{y}_{1},M^{y}_{2},\dots,M^{y}_{d})$ is a quantum
measurement (POVM) performed by Bob upon the choice of the setting $y$. The
goal of SDI certification is then to identify the preparation states and the
measurements, or their properties, based solely on the observed statistics
$p(b|a,x,y)$ _assuming_ the upper bound on the dimension $d$ and the validity
of Born’s rule. We say that certain states $\varrho^{x}_{a}$ and measurements
$\mathbf{M}^{y}$ can be _self-tested_ if the observed statistics specify these
objects uniquely up to a unitary transformation and, perhaps, a global
transposition.

Figure 1: The idea of the certification scheme for the overlaps between
states. Alice chooses her inputs $a,x$ and Bob his input $y$. Alice sends
states $\varrho^{x}_{a}$ to Bob, who produces his outcome $b$. After
establishing perfect correlations between preparation states and measurements
for $y=x$, the rest of the statistics can be used to compute the overlaps
between the preparation states. This can be used to directly certify presence
of coherences between Alice’s input states and between effects of measurements
implemented by Bob.

##  3 Certification of overlaps

We start with a presentation of our scheme for certifying pairwise overlaps
between _pure qudit_ preparation states and between the corresponding
_projective_ measurements. By “corresponding" we mean that in our scheme we
set these states and measurements to be “equal" in a sense that
$\varrho^{x}_{a}=(\varrho^{x}_{a})^{2}=M^{x}_{a}$ for all $a\in[d]$ and
$x\in[n]$. In what follows we will refer to these objects as target pure
states and target measurements respectively. Their experimental counterparts
we denote as $\tilde{\varrho}_{a}^{x}$, and $\mathbf{\tilde{M}}^{y}$
respectively. By “experimental" we mean any states and measurements defined
over the same Hilbert space as the target ones and that reproduce the
_observed statistics_ i.e.,
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$.
Clearly, we _do not_ assume that the experimental states and measurement have
to be “equal".

The idea of our certification scheme is very intuitive, yet powerful (see Fig.
1). Assume that Alice and Bob prepared their devices in a way that
$\tilde{p}(b|a,x,y)=1$, whenever $y=x$ and $b=a$. In other words, outcomes of
Bob’s measurement are _perfectly correlated_ with the preparations of Alice
(whenever $x=y$). Since the quantum dimension is upper bounded by $d$, we can
easily conclude that $\tilde{\varrho}^{x}_{a}=\tilde{M}^{x}_{a}$, for all
$a\in[d]$ and $x\in[n]$. Clearly, after these perfect correlations are
established, the “cross-terms", can be used to compute the overlaps between
the preparation states and between measurement operators:
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{y}_{b})={\mathrm{tr}}(\tilde{M}^{x}_{a}\tilde{M}^{y}_{b})$.
Therefore, if the experimental statistics $\tilde{p}(b|a,x,y)$ match the
target statistics $p(b|a,x,y)$, we can certify that overlaps between
experimental states match those of the target states. The same holds for the
corresponding measurement operators.

Our method can be also applied when experimental statistics do not match
exactly the target ones.

######  Theorem 1 (Robust SDI certification of overlaps).

Consider pure target qudit preparation states $\varrho^{x}_{a}$ and target
projective measurements $\mathbf{M}^{y}$, where $a\in[d]$ and $x,y\in[n]$.
Assume that $\varrho^{x}_{a}=M^{x}_{a}$ for all $a,x$ and furthermore that
experimental states $\tilde{\varrho}^{x}_{a}$ and measurements
$\mathbf{\tilde{M}}^{y}$ act on Hilbert space of dimension at most $d$ and
generate statistics
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$
such that
$|\tilde{p}(b|a,x,y)-{\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})|\leq\varepsilon$,
for all $a,b,x,y$. Then, input states $\tilde{\varrho}^{x}_{a}$ are almost
pure and measurements $\mathbf{\tilde{M}}^{y}$ are almost projective in the
sense that

$\displaystyle\text{for all $x$\ \ }\sum_{a=1}^{d}\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip 1.0pt\right|\right|\geq d(1-2\varepsilon)\
,$ (1) $\displaystyle\text{for all $y$\ \ }\sum_{b=1}^{d}\left|\left|\hskip
1.0pt\tilde{M}^{y}_{b}\hskip 1.0pt\right|\right|\geq d(1-\varepsilon)\ .$ (2)

Moreover, for all $x\neq x^{\prime}$, $a\neq a^{\prime}$, $y\neq y^{\prime}$,
and $b\neq b^{\prime}$, we have

$\displaystyle|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|\leq\varepsilon+\sqrt{2\varepsilon+d^{2}\varepsilon^{2}},$
(3)
$\displaystyle|{\mathrm{tr}}(\tilde{M}^{y}_{b}\tilde{M}^{y^{\prime}}_{b^{\prime}})-{\mathrm{tr}}(M^{y}_{b}M^{y^{\prime}}_{b^{\prime}})|\leq\varepsilon+(1+d\varepsilon)\sqrt{2\varepsilon+d^{2}\varepsilon^{2}}\
.$

###### Proof.

We start with a straightforward proof of Eq.(2). Since $\left|\left|\hskip
1.0pt\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|\geq{\mathrm{tr}}(\tilde{M}^{x}_{a}\varrho)$ for any state
$\varrho$, the relation
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{x}_{a})\geq 1-\varepsilon$,
valid for all $a,x$, implies $\left|\left|\hskip 1.0pt\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|\geq 1-\varepsilon$, for all $a,x$, and hence for all $x$
we have $\sum_{a=1}^{d}\left|\left|\hskip 1.0pt\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|\geq d-d\varepsilon$.

To prove the bounds on the norms of experimental states in Eq. (1), we fix a
setting $x$ and use a decomposition:

$\tilde{M}^{x}_{a}=\lambda_{a}|\phi_{a}\rangle\langle\phi_{a}|+\text{Rest}_{a}\
,$ (4)

where $\lambda_{a}$ is the largest eigenvalue of $\tilde{M}^{x}_{a}$,
$|\phi_{a}\rangle\langle\phi_{a}|$ is the corresponding eigenvector and,
$\text{Rest}_{a}$ is a positive-semi-definite operator satisfying
$\langle\phi_{a}|\text{Rest}_{a}|\phi_{a}\rangle=0$ (to keep the derivations
legible, we omit the index $x$). Using the fact that observed statistics are
$\varepsilon$ close to the target ones we get:

$\displaystyle
1-\varepsilon\leq{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{x}_{a})=\lambda_{a}\langle\phi_{a}|\tilde{\varrho}^{x}_{a}|\phi_{a}\rangle+{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\text{Rest}_{a})\leq\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip
1.0pt\right|\right|+{\mathrm{tr}}(\text{Rest}_{a})\ .$ (5)

By taking the sum over $a$ we obtain $\sum_{a=1}^{d}\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip 1.0pt\right|\right|\geq
d(1-\varepsilon)-\sum_{a=1}^{d}{\mathrm{tr}}(\text{Rest}_{a})$. We can now use
the identity $d=\sum_{a=1}^{d}{\mathrm{tr}}(\tilde{M}^{x}_{a})$, which follows
form that fact that operators $M^{x}_{a}$ form a POVM in $\mathbb{C}^{d}$.
This equality allows us to give a bound:

$d=\sum_{a=1}^{d}\lambda_{a}+\sum_{a=1}^{d}{\mathrm{tr}}(\text{Rest}_{a})\leq
d(1-\varepsilon)+\sum_{a=1}^{d}{\mathrm{tr}}(\text{Rest}_{a})\ ,$ (6)

which is equivalent to $\sum_{a=1}^{d}{\mathrm{tr}}(\text{Rest}_{a})\leq
d\varepsilon$. Inserting this to $\sum_{a=1}^{d}\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip 1.0pt\right|\right|\geq
d(1-\varepsilon)-\sum_{a=1}^{d}{\mathrm{tr}}(\text{Rest}_{a})$ we obtain Eq.
(1).

We now proceed to the proof of Eqs. (3). We start with the one for the
overlaps between preparation states. The proof is given by the following
sequence of inequalities which hold for every $a\neq a^{\prime}$, $x\neq
x^{\prime}$:

$\displaystyle|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|\leq|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\tilde{M}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})|+|{\mathrm{tr}}(\tilde{M}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(M^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|$
$\displaystyle\leq\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|+\varepsilon\leq\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|_{F}+\varepsilon=\sqrt{{\mathrm{tr}}((\tilde{\varrho}^{x}_{a})^{2})+{\mathrm{tr}}((\tilde{M}^{x}_{a})^{2})-2{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{x}_{a})}+\varepsilon$
$\displaystyle\leq\sqrt{1+(1+d^{2}\varepsilon^{2})-2(1-\varepsilon)}=\varepsilon+\sqrt{2\varepsilon+d^{2}\varepsilon^{2}}.$

All of the above inequalities are pretty straightforward apart from
${\mathrm{tr}}((\tilde{M}^{x}_{a})^{2})\leq 1+d^{2}\varepsilon^{2}$ that we
prove below. For this we again use the partial spectral decomposition form Eq.
(4) (without writing the superscript $x$ as above) which implies:

$\displaystyle{\mathrm{tr}}((\tilde{M}^{x}_{a})^{2})=\lambda_{a}^{2}+{\mathrm{tr}}(\text{Rest}_{a}^{2})\leq
1+({\mathrm{tr}}(\text{Rest}_{a}))^{2}\leq 1+d^{2}\varepsilon^{2}.$ (7)

The proof for the overlaps of POVM effects is very similar to the one for
states with the only difference being the following inequality:

$\displaystyle|{\mathrm{tr}}(\tilde{M}^{x}_{a}\tilde{M}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\tilde{M}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})|\leq(1+d\varepsilon)\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|,$ (8)

that is used in the second step in the proof in Eq. (3). One can easily verify
the validity of this inequality by writing once more the decomposition
$\tilde{M}^{x^{\prime}}_{a^{\prime}}=\lambda_{a}|\phi_{a^{\prime}}\rangle\langle\phi_{a^{\prime}}|+\text{Rest}_{a^{\prime}}$
and remembering that ${\mathrm{tr}}(\text{Rest}_{a})\leq d\varepsilon$. ∎

The above result states that if the statistics observed in our certification
scheme vary just a little bit form the target ones, the overlaps of the
experimental states are also close to the overlaps between target states (and
analogously for measurements). To our best knowledge analogous results have
been previously known only for very special symmetric target states and
measurements forming MUBs [22, 23].

In Appendix B (Lemma 1) we improve the above bounds for the case of qubit
systems. Moreover, in Appendix A we prove, by giving explicit examples, that
the bounds in Eq. (3) are tight in the first orders in $\sqrt{\varepsilon}$
and $d$.

##  4 Self-testing of qubits

We now show that certification of overlaps allows to prove robust self-testing
result for _arbitrary_ pure qubit preparations and projective measurements
appearing in our certification scheme.

######  Theorem 2 (Robust self-testing of qubit systems).

Consider target pure qubit states $\varrho^{x}_{a}$ and projective
measurements $\mathbf{M}^{y}$, where $a=1,2$, $x,y\in[n]$, and
$\varrho^{x}_{a}=M^{x}_{a}$ for all $a,x$. Assume that experimental qubit
states $\tilde{\varrho}^{x}_{a}$ and measurements $\mathbf{\tilde{M}}^{y}$
generate statistics
$\tilde{p}(b|a,x,y)=\mathrm{tr}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$
such that
$|\tilde{p}(b|a,x,y)-{\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})|\leq\varepsilon$,
for all $a,b=1,2$ and $x,y\in[n]$. Then, there exist $\varepsilon_{0}$ such
that for $\varepsilon\leq\varepsilon_{0}$ there exist a qubit unitary $U$ such
that

$\displaystyle\frac{1}{2n}\sum_{a,x}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger}\varrho^{x}_{a})\geq
1-f(\varepsilon)\ ,$
$\displaystyle\frac{1}{2n}\sum_{b,y}{\mathrm{tr}}(U(\tilde{M}^{y}_{b})^{(T)}U^{\dagger}M^{y}_{b})\geq
1-g(\varepsilon)\ ,$ (9)

where $(\cdot)^{(T)}$ is the transposition with respect to a fixed basis in
$\mathbb{C}^{2}$ that may have to be applied to all experimental states and
measurements at the same time. Moreover, functions
$f,g:[0,\varepsilon_{0})\rightarrow\mathbb{R}_{+}$ depend solely on the target
states and measurements and, for small $\varepsilon$, have the asymptotics
$f(\varepsilon)\propto\varepsilon,\ g(\varepsilon)\propto\varepsilon$.

In the above we used the fidelity
$F(\varrho,\sigma)={\mathrm{tr}}(\varrho\sigma)$ to indicate the closeness
between _rotated_ experimental states and target pure states (and analogously
for measurements), following existing literature. The case $\varepsilon=0$
corresponds to ideal reconstruction of target states and measurements after
applying a suitable isometry. We remark that we allow only unitary operations
(and possible transposition), as opposed to general channels [22, 25], to be
applied to the experimental states in order to approximate the target states
_as well as possible_.

In Appendix B we give a formal version of the above result and its proof.
Moreover, we present there robustness bounds expressed in terms of the trace
distance and its analogue for measurements [32]. We remark that the functions
$f,g$ become unbounded once Bloch vectors of target qubit states become
singular, i.e., once the target states are close to being aligned in a space
of smaller dimension.

######  Remark.

Certification of overlaps between pure states in general does not allow for
their self-testing in higher-dimensional systems. This is e.g., due to the
existence of unitary inequivalent sets of SIC-POVMs for $d=3$ [33] and MUBs
for $d=4$ [34] (even if we allow for complex conjugation).

###### Sketch of the proof.

We give a full proof for the ideal case ($\varepsilon=0$) below. From Theorem
1 it follows that for all $x,x^{\prime},a,a^{\prime}$ we have
$\tilde{\varrho^{x}_{a}}=\tilde{M}^{x}_{a}$ and
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})={\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})$.
Using the Bloch representation,
$\varrho^{x}_{a}=\frac{1}{2}(\mathbbm{1}+\mathbf{n}^{x}_{a}\cdot\bm{\sigma})$,
we can conclude that also
$\mathbf{\tilde{n}}^{x}_{a}\cdot\mathbf{\tilde{n}}^{x^{\prime}}_{a^{\prime}}=\mathbf{n}^{x}_{a}\cdot\mathbf{n}^{x^{\prime}}_{a^{\prime}}$,
where $\mathbf{n}^{x}_{a}$ and $\tilde{\mathbf{n}}^{x}_{a}$ are Bloch vectors
of $\varrho^{x}_{a}$ and $\tilde{\varrho}^{x}_{a}$ respectively.

Assume now that the vectors
$\mathbf{n}^{1}_{1},\mathbf{n}^{2}_{1},\mathbf{n}^{3}_{1}$ are linear
independent. Let $O$ be the linear transformation defined by
$O\mathbf{n}^{x}_{1}=\mathbf{\tilde{n}}^{x}_{1}$, $x=1,2,3$, and let $L$ be
the matrix whose rows are the vectors $\mathbf{n}^{x}_{1}$, $x=1,2,3$. Then,
we have $LO^{T}OL^{T}=LL^{T}$, and consequently, since $L$ is invertible by
the construction, $O^{T}O=\mathbbm{1}_{3}$, i.e., $O$ is an orthogonal
transformation in $\mathbb{R}^{3}$. It is well-known [35] that if $\det(O)=1$,
there exist a unitary matrix $U$ such that
${\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})U^{\dagger}\varrho^{x}_{a})=1$ for
$x=1,2,3$ and $a=1$. By our assumption all remaining states $\varrho^{x}_{a}$
can be decomposed in the basis
$\\{\mathbbm{1},\varrho^{1}_{1},\varrho^{2}_{1},\varrho^{3}_{1}\\}$, with the
coefficients depending solely on the overlaps
${\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})$. The, the
same $U$ that maps $\tilde{\varrho}^{x}_{1}$ to $\varrho^{x}_{1}$, for
$x=1,2,3$, also connects the remaining pairs of states. Finally, if
$\det(O)=-1$, the transformation $O$ corresponds to application of the
transposition in the standard basis of $\mathbb{C}^{2}$ followed by
application of a unitary operation $U$, determined by $O$ [36].

For the general case $\varepsilon>0$ we consider Cholesky factorisations of
Gram matrices, denoted by $\Gamma$ and $\tilde{\Gamma}$, of Bloch vectors
$\mathbf{n}^{1}_{1},\mathbf{n}^{2}_{1},\mathbf{n}^{3}_{1}$ and their
experimental counterparts respectively. Theorem 1 is then used to bound
$\left|\left|\hskip 1.0pt\Gamma-\tilde{\Gamma}\hskip 1.0pt\right|\right|_{F}$
and utilize results of Ref. [37] to gauge how much the Cholesky decompositions
of $\Gamma$ and $\tilde{\Gamma}$ differ in the Frobenius norm. The latter can
be connected to the average fidelity between the selected target and
experimental states. The robustness for the remaining states follows from the
fact that they can be decomposed (as operators) using the three initially
chosen target states and the identity.

If vectors from the set $\\{{\mathbf{n}^{x}_{a}\\}}$ span two dimensional
space then the same arguments for both $\varepsilon=0$ and $\varepsilon>0$ can
be repeated for the considered subspace. Importantly, in this case additional
transposition is not necessary. ∎

Figure 2: Lower bounds for the average fidelity between the experimental and
target states. Part (a) presents results for $n=2,3$ qubit MUBs. Part (b)
presents results for two qubit bases for different degree of bias $\alpha$.
Value $\alpha=0$ corresponds to two MUBs while $\alpha=1$ gives two identical
bases.

##  5 Examples

We now apply the quantitative formulation of Theorem 2 to lower-bound average
fidelities for different configurations of target states as a function of the
allowed error $\varepsilon$. In Figure 2 we present results for eigenstates of
$n=2$ and $n=3$ Pauli matrices (i.e., states forming $n=2$ and $n=3$ qubit
MUBs), and states belonging to two biased projective measurements satisfying
${\mathrm{tr}}{(\varrho^{1}_{a}\varrho^{2}_{a})}=\frac{1+\alpha}{2}$, $a=1,2$
where we take $\alpha\in[0,1]$.

For $n=2$ MUBs we compare our results with Ref. [22] that aimed at self-
testing of qubit MUBs. The results of Ref. [22] give the upper bound on
average fidelity equal $0.75$ for the deviation of $\simeq 0.1$ in the figure
of merit. In our scheme this happens for $\varepsilon\simeq 0.033$ as shown on
Fig. 2. To test the versatility of our scheme, we also applied it to $n=3$
MUBs, trine and tetrahedral configurations of qubit states. Quantitative
results concerning these examples are listed in Table 1, while detailed
derivations are given in Appendix D. For trine and tetrahedral configurations
the robustness is obtained via the so-called Procrustes problem, described in
Appendix C, based on Ref. [38]. For cases other than two MUBs we cannot make
any comparison with the existing literature since, to our best knowledge,
these case have not been studied previously.

Configuration | $\varepsilon_{0}$ | $C$
---|---|---
2 MUBs | $\approx 0.062$ | $\frac{7}{2}+\sqrt{2}$
3 MUBs | $\approx 0.030$ | $6$
2 biased bases | $\lesssim\frac{4-3\alpha-\sqrt{7-6\alpha}}{18}$ | $2+\left(1+\frac{\sqrt{1+\alpha}}{\sqrt{2}(1-\alpha)}\right)^{2}$
Trine | $\approx 0.058$ | $\frac{19}{3}$
Tetrahedron | $\approx 0.037$ | $10$
Table 1:  Results of quantitative variant of Theorem 2 applied to different
configurations of target quantum states. The threshold $\varepsilon_{0}$ sets
the maximal noise level which is tolerated by our scheme. The constant $C$ is
defined via the relation
$f(\varepsilon)\stackrel{{\scriptstyle\varepsilon\rightarrow
0}}{{\approx}}C\varepsilon$, where $1-f$ is the lower bound on the average
fidelity from Eq. (2).

##  6 Shared randomness

Throughout the article we assumed that preparation and measurement devices are
uncorrelated, which is often true in practice. However, one may also consider
a situation in which Alice and Bob share a random variable $\lambda$, that
they can use to coordinate their actions in a given round of the experiment.
The most general statistics that can be generated in such a scenario can be
expressed as $p(b|a,x,y)=\int d\lambda
p(\lambda){\mathrm{tr}}(\varrho^{x}_{a}(\lambda)M^{y}_{b}(\lambda))$, where
$p(\lambda)$ denotes the probability distribution of $\lambda$.

Interestingly, the presence of shared randomness makes our main results, as
stated by Theorems 1 and 2, inapplicable. The easiest way to see it is to
consider the following simple example of $n=d=2$, one bit of shared randomness
$\lambda\in\\{1,2\\}$, and $\tilde{\varrho}^{1,2}_{a}=|a\rangle\langle a|$,
$\tilde{M}^{1,2}_{b}=|b\rangle\langle b|$. This clearly satisfies the
requirement on $\tilde{p}(a|a,x,x)=1$. Now, Alice’s and Bob’s devices can
decide to “flip" their preparations and measurement operators whenever $x=2$,
$y=2$ and $\lambda=2$ (note that $\lambda$ can be distributed in an arbitrary
way). This procedure does not affect the correlations $\tilde{p}(a|a,x,y=x)$,
but it can be used to set $\tilde{p}(a|a,1,2)$ to an arbitrary value.

Additionally, we believe that in the presence of shared randomness, the usual
notion of self-testing in the SDI setting has to be reconsidered.
Specifically, for _any_ quantum realisation giving the statistics
$p(b|a,x,y)={\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})$, one can always consider
a strategy in which with probability $p(\lambda)$ Alice prepares states
$\varrho^{x}_{a}(\lambda)=U_{\lambda}\varrho^{x}_{a}U^{\dagger}_{\lambda}$ and
Bob implements measurements
$M^{y}_{b}(\lambda)=U_{\lambda}M^{y}_{b}U^{\dagger}_{\lambda}$, where
$U_{\lambda}$ is some arbitrary unitary transformation depending on $\lambda$.
Clearly, such strategy reproduces the original statistics, and makes it
impossible to find a single unitary that connects the target and experimental
states (or measurements).

While, we do not consider the assumption on no shared randomness strong, below
we propose a modification of our SDI certification scheme that allows to
certify overlaps between arbitrary pure states even in the presence of shared
randomness. The idea is to introduce, for every pair of non-orthogonal states,
a suitable _intermediate state_ [39], that enforces fixed overlaps between
experimental states in every round of the experiment.

Recall that in our certification scheme we consider pure target qudit states
$\varrho^{x}_{a}$ and target projective measurements $\mathbf{M}^{y}$, where
$a\in[d]$ and $x,y\in[n]$. We extend this scheme by introducing an additional
intermediate target state $\varrho_{z}$ for every two pairs of Alice’s input
$(x,a)$ , $(x^{\prime},a^{\prime})$, where $a\neq a^{\prime}$ and $x\neq
x^{\prime}$. The state $\varrho_{z}$ is chosen as the unique state satisfying:

${\mathrm{tr}}\left(\varrho_{z}(\varrho^{x}_{a}+\varrho^{x^{\prime}}_{a^{\prime}})\right)=1+{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})^{\frac{1}{2}}\
.$ (10)

Let $\tilde{p}(b|z,y)$ denote the experimental statistics, corresponding to
inputs $z,y$ and outcome $b$. In the ideal scenario, in which experimental
statistics satisfy assumptions of Theorem 1 with $\varepsilon=0$, assume
additionally that
$\tilde{p}(a|z,x)+\tilde{p}(a^{\prime}|z,x^{\prime})=1+\sqrt{\tilde{p}(a^{\prime}|x,a,x^{\prime})}\
.$ Below we prove that under the above assumptions for all $\lambda$ we have
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda))={\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})$.
Since in the ideal case overlaps between preparation states _do not depend_ on
$\lambda$ and match the target value, our modified certification scheme works
also in the presence of shared randomness. In the noisy scenario one can
generally expect to be able to upper-bound the probability (over the choice of
$\lambda$) of the deviations from perfect correlations between states and
measurements, as a function of the error $\varepsilon$. We postpone the
discussion of this interesting point to the future work.

###### Proof of soundness of certification in the presence of shared
randomness.

Consider a general physical realisation of experimental statistics
$\tilde{p}(b|x,a,y)$, $\tilde{p}(b|z,y)$ via classically correlated
preparations and measurements on $\mathbb{C}^{d}$:

$\displaystyle\tilde{p}(b|x,a,y)$ $\displaystyle=\int d\lambda
p(\lambda){\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{M}^{y}_{b}(\lambda))\
,$ (11) $\displaystyle\tilde{p}(b|z,y)$ $\displaystyle=\int d\lambda
p(\lambda){\mathrm{tr}}(\tilde{\varrho}_{z}(\lambda)\tilde{M}^{y}_{b}(\lambda))\
,$ (12)

where $x,y\in[n]$, $a,b\in[d]$ and the variable (input) $z$ labels elements in
the set of unordered pairs $\\{{(x,a),(x^{\prime},a^{\prime})\\}}$, where
$x\neq x^{\prime}$ and $a\neq a^{\prime}$. Assume now that the above
statistics match _exactly_ the ones required by our scheme and are compatible
with Eq. (10) in the sense that:

$\displaystyle\tilde{p}(b|x,a,y)$ $\displaystyle=$
$\displaystyle{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{y}_{b})\ ,$ (13)
$\displaystyle\tilde{p}(a|z,x)+\tilde{p}(a^{\prime}|z,x^{\prime})$
$\displaystyle=$ $\displaystyle 1+\sqrt{\tilde{p}(a^{\prime}|x,a,x^{\prime})}\
.$ (14)

In what follows we prove that if the above constraints are satisfied, then for
all111Technically speaking this equality holds for _almost all_ $\lambda$.:

$\displaystyle\tilde{\varrho}^{x}_{a}(\lambda)$ $\displaystyle=$
$\displaystyle\tilde{M}^{x}_{a}(\lambda)\ ,$ (15)
$\displaystyle{\mathrm{tr}}(\varrho^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})$
$\displaystyle=$
$\displaystyle{\mathrm{tr}}(\varrho^{x}_{a}(\lambda)\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda))\
,$ (16)

where $x,x^{\prime}\in[n]$, $a,a^{\prime}\in[d]$. The above equation means
that the overlaps between preparation sates (and measurements) do not depend
on on the value of the shared random variable $\lambda$.

The proof of Eq. (15) is straightforward. Namely, form Eq. (13) if follows
that:

$\int d\lambda
p(\lambda){\mathrm{tr}}\left(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{M}^{x}_{a}(\lambda)\right)=1\
.$ (17)

Since
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{M}^{x}_{a}(\lambda))\leq
1$, we get that for all $\lambda$ we have
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{M}^{x}_{a}(\lambda))=1$.
Since this reasoning can be repeated for all $a\in[d]$ (for the fixed value of
$x$) and since operators $\tilde{M}^{x}_{a}$ form a POVM on $\mathbb{C}^{d}$ ,
we finally get Eq. (15).

The proof of Eq. (16) is more involved and relies on both Eqs. (13) and (14).
Specifically, from Eq. (13) and the already established identity
$\tilde{\varrho}^{x}_{a}=\tilde{M}^{x}_{a}$, it follows that Eq. (14) is
equivalent to:

$\int d\lambda
p(\lambda){\mathrm{tr}}\left(\tilde{\varrho}_{z}(\lambda)[\tilde{\varrho}^{x}_{a}(\lambda)+\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda)]\right)=1+\sqrt{{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})}\
,$ (18)

where moreover:

$\int d\lambda
p(\lambda){\mathrm{tr}}\left(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda)\right)={\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})\
.$ (19)

Using Bloch representation of qubit states (states
$\tilde{\varrho}^{x}_{a}(\lambda),\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda)$
are pure and hence span a two dimensional subspace of $\mathbb{C}^{d}$) it is
straightforward to obtain the bound:

${\mathrm{tr}}\left(\tilde{\varrho}_{z}(\lambda)[\tilde{\varrho}^{x}_{a}(\lambda)+\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda)]\right)\leq
1+\sqrt{{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda))}\
.$ (20)

After setting
$g(\lambda)=\sqrt{{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}(\lambda)\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}(\lambda))}$
and using Eq. (19) we obtain:

$\int d\lambda p(\lambda)g(\lambda)\geq\sqrt{\int d\lambda
p(\lambda)g(\lambda)^{2}}\ .$ (21)

Using Chauchy-Schwartz inequality for the left hand side of the above
inequality we finally get:

$\int d\lambda p(\lambda)g(\lambda)=\sqrt{\int d\lambda
p(\lambda)g(\lambda)^{2}}\ ,$ (22)

which is equivalent to saying that the variance of the random variable
$g(\lambda)$ vanishes. Therefore, $g(\lambda)=\alpha$, where $\alpha$ is some
numerical constant. We conclude the proof by noting that from Eq. (19) it
follows that
$\alpha^{2}={\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})$
which implies Eq. (16). ∎

##  7 Extension of the scheme to SDI characterization of rank-1 non-projective
measurements

In our work we focus on self-testing of arbitrary configurations of
preparation states and projective measurements. In this section we discuss
possible extension of our results to certification of properties of general
rank-1 measurements in arbitrary dimension $d$. To that end, we introduce
additional states coming from suitably-chosen orthonormal bases per every
effect of the measurement we want to certify. With the help of these extra
states (and corresponding projective measurements) we can e.g., certify the
following properties of POVMs in arbitrary dimension $d$: (i) information
completeness, (ii) extremality of rank-1 POVMs, and (iii) SIC-property.

The main idea of the extension is the following. Suppose we would like to
self-test an $m$-outcome POVM $\mathbf{N}=(N_{1},\dots,N_{m})$ acting on a
Hilbert space of dimension $d<m$ with each effect $N_{b}$ being of rank $1$
i.e., $N_{b}=\alpha_{b}\sigma_{b}$, for some pure states $\sigma_{b}$ and
positive scalars $\alpha_{b}$. We consider $m$ sets of preparation states
$\\{\varrho^{x}_{a}\\}$, $a\in[d]$, ${x\in[m]}$ and the corresponding
projective measurements $\mathbf{M}^{y}$, $y\in[m]$. Assume now that in the
idealized case we observe statistics satisfying:

$\displaystyle{\mathrm{tr}}(\varrho^{b}_{a}N_{b})=0,\quad\forall
b\in[m],\,a\neq 1,$
$\displaystyle{\mathrm{tr}}(\varrho^{x}_{a}M^{x}_{b})=\delta_{a,b},\quad\forall
x\in[m]\ .$ (23)

From the second set of conditions we conclude that $\\{\varrho^{x}_{a}\\}_{a}$
forms, as before, a basis of pure states for each $x$. From the first
condition in Eq. (7) we get that $\mathbf{N}$ has to be a rank-1 POVM and,
moreover, that $\varrho^{b}_{1}=\sigma_{b}$ ($N_{b}$ is orthogonal to all of
the states except the one with $a=1$). Consequently, we have
${\mathrm{tr}}(\varrho^{b}_{1}N_{b})=\alpha_{b}$ for all $b\in[m]$. This fact
together with the measured overlaps
${\mathrm{tr}}(\varrho^{b}_{1}\varrho^{b^{\prime}}_{1})={\mathrm{tr}}(\varrho^{b}_{1}M^{b^{\prime}}_{1})$
gives the information about the overlaps ${\mathrm{tr}}(N_{b}N_{b^{\prime}})$
between effects of the POVM $\mathbf{N}$.

The knowledge of all overlaps ${\mathrm{tr}}(N_{b}N_{b^{\prime}})$ (including
$b^{\prime}=b$) and the fact that $\mathbf{N}$ is a rank-1 POVM allow to
certify many interesting properties of $\mathbf{N}$. First, as long as $m\leq
d^{2}$ we can certify extremality of $\mathbf{N}$. This follows from the fact
that rank-1 POVMs are extremal if and only if their effects are linearly
independent (see e.g., Ref. [40]). Linear independence can be directly
inferred from the experimentally accessible moment matrix
$\Gamma_{b,b^{\prime}}={\mathrm{tr}}(N_{b}N_{b^{\prime}})$. Specifically,
$\Gamma$ is non-singular if and only if operators $\\{N_{b}\\}_{b}$ are linear
independent. Now, if $m=d^{2}$ we can use the same reasoning to infer
information completeness of $\mathbf{N}$. Finally, from the moment matrix
$\Gamma$ we can directly verify the SIC property i.e., the condition
${\mathrm{tr}}(N_{b}N_{b^{\prime}})=\frac{1+\delta_{b,b^{\prime}}}{d^{2}(d+1)},\forall
b,b^{\prime}$.

##  8 Discussion

We have presented a systematic analytical scheme for noise-resilient
certification of overlaps between arbitrary configurations of pure quantum
states and rank-$1$ projective measurements in the prepare-and-measure
scenario. For qubits our scheme can be used to robustly self-test general
ensembles of pure quantum sates and the corresponding projective measurements.
We believe that these findings pave the way towards systematic certification
of general quantum systems in the semi-device-independent paradigm. This is
supported by the concrete qubit results from Table 1 and by the universality
of our protocol for certifying overlaps of quantum sates in arbitrary
dimension.

There are many exciting research directions that stem from our work. One of
them is a systematic SDI characterization of generalized measurements beyond
qubits. Another research direction is the robustness analysis of our scheme in
the presence of shared randomness. An important contribution to the filed of
self-testing would be a generalization of our self-testing results for qubit
to higher-dimensional systems. Finally, it is interesting to combine our
methods with recent findings relating state discrimination games and resource
theories of measurements and channels [41, 42, 43].

We would also like to draw reader’s attention to a related work [44].

## Acknowledgements

We thank Marcin Pawłowski, Jedrzej Kaniewski and Tanmoy Biswas for interesting
discussions and comments. NM acknowledges the financial support by First TEAM
Grant No. 2016-1/5. We acknowledge partial support by the Foundation for
Polish Science (IRAP project, ICTQT, contract no. 2018/MAB/5, co-financed by
EU within Smart Growth Operational Programme). We acknowledge the support from
International Research Agendas Programme of the Foundation for Polish Science
funded by the European Regional Development Fund. MO acknowledges the
financial support by TEAM-NET project (contract no.
POIR.04.04.00-00-17C1/18-00).

## Appendix

In this Appendix we provide technical details that were omitted in the main
text. First, in Appendix A we prove statements regarding saturation of the
bounds in Theorem 1. Second, in Appendix B we formulate and prove a
quantitative version of Theorem 2, which concerns robustness of our self-
testing protocol expressed in terms of average fidelities. We also prove there
Theorem 3 that gives analogous results expressed in terms of trace distance.
In Appendix C we provide alternative robustness analysis based on orthogonal
Procrustes problem. In Appendix D we provide technical details about the
considered examples in Table 1. Finally, in Appendix E we give proofs of two
auxiliary Lemmas needed in the proof of the technical version of Theorem 2.

##  Appendix A Explicit form of states and measurements saturating the bounds
in Theorem 1

We start by giving an example of four states and two qubit measurements
($n=2$, $d=2$), for which the bound in Eq. (3) on the deviation of the
overlaps scales like $\sqrt{\varepsilon}$. Their explicit form is given below:

$\displaystyle\varrho^{1}_{1}=|\tilde{0}\rangle\langle\tilde{0}|,\quad|\tilde{0}\rangle=\sqrt{1-\varepsilon}|0\rangle+\sqrt{\varepsilon}|1\rangle,\quad\varrho^{2}_{1}=|+\rangle\langle+|,$
(24)

where $|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$, and as required by
our scheme $\varrho^{x}_{2}=\mathbbm{1}-\varrho^{x}_{1},\;x=1,2$ and
$M^{x}_{a}=\varrho^{x}_{a}$, for all $x,a$. Let the experimental states and
measurements be the following:

$\displaystyle\tilde{\varrho}^{1}_{1}=|0\rangle\langle
0|,\quad\tilde{\varrho}^{1}_{2}=|1\rangle\langle
1|,\quad\tilde{\varrho}^{2}_{1}=|+\rangle\langle+|,\quad\tilde{\varrho}^{2}_{2}=|-\rangle\langle-|,$
$\displaystyle\tilde{M}^{1}_{1}=|\tilde{0}\rangle\langle\tilde{0}|,\quad\tilde{M}^{2}_{1}=|\tilde{+}\rangle\langle\tilde{+}|,\quad|\tilde{+}\rangle=\sqrt{1-\varepsilon}|+\rangle+\sqrt{\varepsilon}|-\rangle\
.$ (25)

First, we need to show that the experimental statistics is deviated from the
target one by at most $\varepsilon$, i.e.,
$|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{y}_{b})|\leq\varepsilon,\forall
a,b,x,y$. Since $\tilde{\varrho}^{x}_{1}+\tilde{\varrho}^{x}_{2}=\mathbbm{1}$
for both $x=1,2$, it is sufficient to consider the cases of $a=1,b=1$.
Calculating the statistics we can see that
${\mathrm{tr}}(\tilde{\varrho}^{1}_{1}\tilde{M}^{1}_{1})={\mathrm{tr}}(\tilde{\varrho}^{2}_{1}\tilde{M}^{2}_{1})=1-\varepsilon$,
and
${\mathrm{tr}}(\tilde{\varrho}^{1}_{1}\tilde{M}^{2}_{1})={\mathrm{tr}}(\tilde{\varrho}^{2}_{1}\tilde{M}^{1}_{1})=\frac{1}{2}+\sqrt{\varepsilon-\varepsilon^{2}}$
and the latter is just the same as the target probability,
${\mathrm{tr}}(\varrho^{1}_{1}\varrho^{2}_{1})$. Finally, calculation of the
overlap between the states yields:

$\displaystyle|{\mathrm{tr}}(\varrho^{1}_{1}\varrho^{2}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{1}_{1}\tilde{\varrho}^{2}_{1})|=\sqrt{\varepsilon-\varepsilon^{2}},$
(26)

which confirms the scaling of the first order of $\sqrt{\varepsilon}$.

Now, let us show that the bounds from Eq. (3) are also tight in the first
order of $d$. We still consider the case of $n=2$, but now $d$ is arbitrary.
Let us consider the following POVM:

$M_{1}=|1\rangle\langle 1|+\varepsilon(d-1)|+\rangle\langle+|\ ,\
M_{a}=(|a\rangle-\delta|+\rangle)(\langle a|-\delta\langle+|)\ a=2,3,\dots,d\
,$ (27)

where $\delta=\frac{1}{\sqrt{d-1}}+\sqrt{\frac{1}{d-1}-\varepsilon}$,
$\\{|a\rangle\\}_{a=1}^{d}$ is the computational basis of $\mathbb{C}^{d}$ and
$|+\rangle=\frac{1}{\sqrt{d-1}}\sum_{a=2}^{d}|a\rangle$ is the "maximally
coherent" state in the subspace spanned by $\\{|i\rangle\\}_{i=2}^{d}$. In the
proof of Theorem 1 the quantity
$|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|$
is upper-bounded by $\varepsilon+\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|$, for
which there exist a state $\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}$ and an
effect $\tilde{M}^{x^{\prime}}_{a^{\prime}}$ reaching the bound. Now if we
take $\mathbf{\tilde{M}}^{x}$ to be the POVM we just introduced, and
$\tilde{\varrho}^{x}_{1}=|1\rangle\langle 1|$,
$\tilde{\varrho}^{x}_{a}=\frac{M_{a}}{{\mathrm{tr}}(M_{a})}$, $a=2,3,\dots,d$,
the conditions of the Theorem 1 will be satisfied and the resulting bound will
be $d\varepsilon$.

##  Appendix B Quantitative statement of the robust self-testing

In this part we first formulate Theorem 3 in which we present robustness of
our self-testing scheme in terms of trace distance222Recall that the trace
distance is defined via
$\mathrm{d}_{\mathrm{tr}}(\sigma,\varrho)=\frac{1}{2}\left|\left|\hskip
1.0pt\sigma-\varrho\hskip
1.0pt\right|\right|_{1}=\frac{1}{2}{\mathrm{tr}}\sqrt{(\sigma-\varrho)^{\dagger}(\tilde{\varrho}-\varrho)}$
and it has a neat operational interpretation in term of optimal success
probability $p$ of distinguishing $\varrho$ and $\sigma$ via the most general
quantum measurements:
$p=\frac{1}{2}(1+\mathrm{d}_{\mathrm{tr}}(\sigma,\varrho))$.  for states and
operator norm for measurements. Then, we give a quantitative statement of
Theorem 2, that concerns robustness of our self-testing scheme in terms of
average fidelities. We then proceed with proofs of both results.

In what follows we will need the following definition.

######  Definition 1.

Let $\\{\mathbf{n}_{i}\\}_{i\in[n]}$ be a set of $n$ vectors in
$\mathbb{R}^{l}$. Matrix $\Gamma\in\mathbb{R}^{n\times n}$ is called the Gram
matrix of the set $\\{\mathbf{n}_{i}\\}_{i\in[n]}$, if its elements are given
by $\Gamma_{i,j}=\mathbf{n}_{i}\cdot\mathbf{n}_{j}$, $i,j\in[n]$.

######  Theorem 3 (Robust self-testing for qubits via trace distance and
operator norm).

Consider pure target qubit preparation states $\varrho^{x}_{a}$ and target
projective measurements $\mathbf{M}^{y}$, where $a=1,2$ and $x,y\in[n]$.
Assume that $\varrho^{x}_{a}=M^{x}_{a}$ for all $a,x$ and, furthermore, that
experimental states $\tilde{\varrho}^{x}_{a}$ and measurements
$\mathbf{\tilde{M}}^{y}$ act on Hilbert space of dimension at most $d$ and
generate statistics
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$
such that
$|\tilde{p}(b|a,x,y)-{\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})|\leq\varepsilon$,
for all $a,b,x,y$.

Let $k\in\\{{2,3\\}}$ be the cardinality of the maximal set of linearly
independent Bloch vectors of states $\varrho^{x}_{1}$. Fix a set $S\subset[n]$
of $k$ linearly independent vectors $\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$,
construct their Gram matrix $\Gamma_{S}$, and let $L_{S}$ be a Cholesky factor
of $\Gamma_{S}$ (i.e., $\Gamma_{S}=L_{S}L^{T}_{S}$ and $L_{S}$ is lower
triangular). For every $x\in[n]\setminus S$ and both $a=1,2$ let
$\mathbf{c}^{x,a}_{S}$, denote the coefficients of decomposition of
$\mathbf{n}^{x}_{a}$ as a linear combination of
$\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$. Finally, let us define three auxiliary
functions

$\displaystyle
F_{k}(\varepsilon)\coloneqq\sqrt{\varepsilon}\sqrt{4k(k-1)}\sqrt{1+2\sqrt{\varepsilon}+\frac{k+3}{k-1}\varepsilon}\
,$ $\displaystyle O_{k}(\varepsilon)\coloneqq
2((k-1)\sqrt{\varepsilon}+(k+1)\varepsilon)\ ,$ (28) $\displaystyle
E_{S,k}(\varepsilon)\coloneqq\frac{1}{2\sqrt{2}}\frac{\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|F_{k}(\varepsilon)}{\sqrt{1-\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|O_{k}(\varepsilon)}}\min\left[\frac{\left|\left|\hskip
1.0pt\Gamma_{S}\hskip 1.0pt\right|\right|}{\left|\left|\hskip
1.0pt\Gamma_{S}\hskip 1.0pt\right|\right|_{F}},\frac{\left|\left|\hskip
1.0ptL_{S}\hskip 1.0pt\right|\right|}{\sqrt{k}}\right]\ .$

Then, there is a region of $\varepsilon\in[0,\varepsilon_{0})$ determined by

$\displaystyle\left|\left|\hskip 1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|O_{k}(\varepsilon)\leq 1,$ (29)

for which there exist a qubit unitary matrix $U$ such that

$\displaystyle\frac{1}{k}\sum_{x\in
S}\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger},\varrho^{x}_{1})\leq
E_{S,k}(\varepsilon)\ ,$ (30a) $\displaystyle\frac{1}{k}\sum_{x\in
S}\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{2})^{(T)}U^{\dagger},\varrho^{x}_{2})\leq
E_{S,k}(\varepsilon)+2\sqrt{\varepsilon}\ ,$ (30b)
$\displaystyle\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger},\varrho^{x}_{a})\leq\sqrt{k}{|\mathbf{c}^{a,x}_{S}|E_{S,k}(\varepsilon)}+\frac{\sqrt{k}}{2}\left(|\mathbf{c}^{a,x}_{S}|+\frac{\sqrt{k}}{k-1}\right)\frac{\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|O_{k}(\varepsilon)}{1-\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip 1.0pt\right|\right|O_{k}(\varepsilon)},$ (30c)
$\displaystyle\text{for}\;x\notin S,\ a=1,2\ ,$
$\displaystyle\frac{1}{k}\sum_{y\in S}\left|\left|\hskip
1.0ptU(\tilde{M}^{y}_{1})^{(T)}U^{\dagger}-M^{y}_{1}\hskip
1.0pt\right|\right|\leq E_{S,k}(\varepsilon)+\sqrt{\varepsilon}\ ,$ (30d)
$\displaystyle\left|\left|\hskip
1.0ptU(\tilde{M}^{y}_{1})^{(T)}U^{\dagger}-M^{y}_{1}\hskip
1.0pt\right|\right|\leq\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{y}_{1})^{(T)}U^{\dagger},\varrho^{y}_{1})+\sqrt{\varepsilon},\quad\text{for}\;y\notin
S\ ,$ (30e)

where $(\cdot)^{(T)}$ is the transposition (with respect a fixed basis in
$\mathbb{C}^{2}$) that may have to be applied to all experimental states and
measurements simultaneously.

The reason for such a formulation of the theorem lays in the proof techniques
that were used to derive it. Namely, we use results on the stability of the
Cholesky factorization. Since the result employed by us (Ref. [37]) is valid
for positive-definite matrices we cannot apply it to the Gram matrix of all of
the vectors $\mathbf{n}^{x}_{a}$ due to their linear dependence. As a result,
our bounds depend on the particular choice of the set $S$ of states whose
Bloch vectors are linearly independent. In what follows, without the loss of
generality we take $S=\\{1,2\\}$ for $k=2$ and $S=\\{1,2,3\\}$ for $k=3$ in
the proof. In practice, however, one should consider all subsets $S$ of states
that give non-singular Gram matrices having in mind a specific application.

_Outline of the proof.–_ In what follows we will consider a particular subset
$\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$ of $k$ linearly independent Bloch vectors.
For simplify we will omit the subscript $S$ whenever possible. The main idea
of the proof is to use a Cholesky factors of $k\times k$ the Gram matrices
$\Gamma$ and $\tilde{\Gamma}$ of target ($\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$)
and experimental ($\\{\mathbf{\tilde{n}}^{x}_{1}\\}_{x\in S}$) Bloch vectors
respectively. We then make use of the result of Ref. [37] on the stability of
Cholesky factorization which, in simple terms, states that if $\Gamma$ and
$\tilde{\Gamma}$ are close to each other, so are their Cholesky factors $L$
and $\tilde{L}$. Specifically, this result, which we quote bellow (see Theorem
[Sun 1991]), sets an upper bound on the Frobenius norm of $\Delta
L=\tilde{L}-L$ in terms of the Frobenius norm of the perturbation
$\Delta\Gamma=\tilde{\Gamma}-\Gamma$. The Frobeniuous norm of the perturbation
$\Delta L$ can be connected to the trace distance between states
$\varrho^{x}_{1}$ and the rotated states $\tilde{\varrho}^{x}_{1}$ in the
selected subset $S$. On the other hand, the bound on $\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|_{F}$ can be estimated from our
assumption: $|p(b|a,x,y)-\tilde{p}(b|a,x,y)|\leq\varepsilon$. This bound
follows directly from the results of Theorem 1, but here we use an improved
qubit stability bounds, which are given by Lemma 1. This, in turn, leads to
stronger estimates for the norms of $\Delta\Gamma$ in Lemma 2 bellow.
Combining all these results produces bounds in Eq. (30a).

In the next part of the proof we determine the trace distance between states
$\varrho^{x}_{2}$ and $\tilde{\varrho}^{x}_{2}$ for $x\in S$ based solely on
the fact that $\varrho^{x}_{2}=\mathbbm{1}-\varrho^{x}_{1}$ and
$\tilde{\varrho}^{x}_{2}$ are _close to_
$\mathbbm{1}-\tilde{\varrho}^{x}_{1}$. This gives the bounds in Eq. (30b). For
states $x\notin S$ we use the fact that they can be decomposed in the basis of
the states in $S$. Since this linear decomposition can be different for target
and experimental states, we use the result on stability of linear systems (see
Theorem [Higham 2002], which we also quote below). The bounds for the states
$x\notin S$ are given by Eq. (30c). Finally, we connect the bounds for the
distance between the target and experimental measurements and the distance
between the corresponding states resulting in Eq. (30d) and (30e).

######  Theorem 2 (Quantitative formulation of robust self-testing for
qubits).

Consider pure target qubit preparation states $\varrho^{x}_{a}$ and target
projective measurements $\mathbf{M}^{y}$, where $a=1,2$ and $x,y\in[n]$.
Assume that $\varrho^{x}_{a}=M^{x}_{a}$ for all $a,x$ and, furthermore, that
experimental states $\tilde{\varrho}^{x}_{a}$ and measurements
$\mathbf{\tilde{M}}^{y}$ act on Hilbert space of dimension at most $d$ and
generate statistics
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$
such that
$|\tilde{p}(b|a,x,y)-{\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})|\leq\varepsilon$,
for all $a,b,x,y$.

Let $k\in\\{{2,3\\}}$ be the cardinality of the maximal set of linearly
independent Bloch vectors of states $\varrho^{x}_{1}$. Fix a set $S\subset[n]$
of $k$ linearly independent vectors $\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$,
construct their Gram matrix $\Gamma_{S}$, and let $L_{S}$ be a Cholesky factor
of $\Gamma_{S}$ (i.e., $\Gamma_{S}=L_{S}L^{T}_{S}$ and $L_{S}$ is lower
triangular). For every $x\in[n]\setminus S$ and both $a=1,2$ let
$\mathbf{c}^{x,a}_{S}$, denote the coefficients of decomposition of
$\mathbf{n}^{x}_{a}$ as a linear combination of
$\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$. Finally, let us define three auxiliary
functions

$\displaystyle
F_{k}(\varepsilon)\coloneqq\sqrt{\varepsilon}\sqrt{4k(k-1)}\sqrt{1+2\sqrt{\varepsilon}+\frac{k+3}{k-1}\varepsilon}\
,$ $\displaystyle O_{k}(\varepsilon)\coloneqq
2((k-1)\sqrt{\varepsilon}+(k+1)\varepsilon)\ ,$ $\displaystyle
E_{S,k}(\varepsilon)\coloneqq\frac{1}{2\sqrt{2}}\frac{\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|F_{k}(\varepsilon)}{\sqrt{1-\left|\left|\hskip
1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|O_{k}(\varepsilon)}}\min\left[\frac{\left|\left|\hskip
1.0pt\Gamma_{S}\hskip 1.0pt\right|\right|}{\left|\left|\hskip
1.0pt\Gamma_{S}\hskip 1.0pt\right|\right|_{F}},\frac{\left|\left|\hskip
1.0ptL_{S}\hskip 1.0pt\right|\right|}{\sqrt{k}}\right]\ .$ (31)

Then, there is a region of $\varepsilon\in[0,\varepsilon_{0})$ determined by

$\displaystyle\left|\left|\hskip 1.0pt\Gamma_{S}^{-1}\hskip
1.0pt\right|\right|O_{k}(\varepsilon)\leq 1,$ (32)

for which there exist a qubit unitary matrix $U$ such that

$\displaystyle\frac{1}{k}\sum_{x\in
S}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger}\varrho^{x}_{1})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-E_{S,k}(\varepsilon)^{2}\
,$ (33a) $\displaystyle\frac{1}{k}\sum_{x\in
S}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{2})^{(T)}U^{\dagger}\varrho^{x}_{2})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-(E_{S,k}(\varepsilon)+2\sqrt{\varepsilon})^{2}\
,$ (33b)
$\displaystyle{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger}\varrho^{x}_{a})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-\left(\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger},\varrho^{x}_{a})\right)^{2},\quad\text{for}\;x\notin
S,\ a=1,2\ ,$ (33c) $\displaystyle\frac{1}{k}\sum_{y\in
S}{\mathrm{tr}}(U(\tilde{M}^{y}_{1})^{(T)}U^{\dagger}M^{y}_{1})\geq
1-\frac{5}{2}\varepsilon-\left(\sqrt{2}E_{S,k}(\varepsilon)+\sqrt{\varepsilon}\right)^{2}\
,$ (33d)
$\displaystyle{\mathrm{tr}}(U(\tilde{M}^{y}_{1})^{(T)}U^{\dagger}M^{y}_{1})\geq
1-\varepsilon-\left(\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{y}_{1})^{(T)}U^{\dagger},\varrho^{y}_{1})+\sqrt{\varepsilon}\right)^{2},\quad\text{for}\;y\notin
S\ ,$ (33e)

where $(\cdot)^{(T)}$ is the transposition (with respect a fixed basis in
$\mathbb{C}^{2}$) that may have to be applied to all experimental states and
measurements simultaneously. Note that in formulas (33c) and (30e) we have
used trace distance
$\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger},\varrho^{x}_{a})$
in order to simplify the resulting formulas. The bound on this quantity is
given in Eq. (30c).

######  Remark.

Results of Theorem 2 are obtained using a similar reasoning to that given in
the outline of the proof of Theorem 3. However, the proof steps are
supplemented by bounds connecting fidelity to the trace distance (for states)
and operator norm (for measurement operators). The proof of Theorem 2 is given
at the end of this part of the Appendix.

Before we proceed we state two theorems from the literature that are needed
for our proof. First, we repeat the statement of the Theorem 1.4 from Ref.
[37] for real-valued matrices. Second, we state a result concerning stability
of systems of linear equation, which we borrow from Ref. [45] (Theorem 7.2).
We changed the notation used in these papers in accordance with our paper.

######  Theorem (Sun 1991 - Stability of Cholesky factorisation).

Let $\Gamma$ be an $k\times k$ positive definite matrix and $\Gamma=LL^{T}$
its Cholesky factorization. If $\Delta\Gamma$ is an $k\times k$ symmetric
matrix satisfying

$\displaystyle\left|\left|\hskip 1.0pt\Gamma^{-1}\hskip
1.0pt\right|\right|\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|\leq 1,$ (34)

then there is a unique Cholesky factorization

$\displaystyle\Gamma+\Delta\Gamma=(L+\Delta L)(L+\Delta L)^{T},$

and

$\displaystyle\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}\leq\frac{\left|\left|\hskip 1.0pt\Gamma^{-1}\hskip
1.0pt\right|\right|\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|_{F}}{\sqrt{2(1-\left|\left|\hskip 1.0pt\Gamma^{-1}\hskip
1.0pt\right|\right|\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|})}\min\left[\frac{\left|\left|\hskip 1.0ptL\hskip
1.0pt\right|\right|_{F}\left|\left|\hskip 1.0pt\Gamma\hskip
1.0pt\right|\right|}{\left|\left|\hskip 1.0pt\Gamma\hskip
1.0pt\right|\right|_{F}},\left|\left|\hskip 1.0ptL\hskip
1.0pt\right|\right|\right].$ (35)

######  Remark.

This theorem dates back to 1991, and, of course there have been attempts to
improve this result. However, the recent review [46] on this topic suggests
that the bound given in the Theorem above is the most appropriate in our case
(Remark 3.2, Ref. [46]).

######  Theorem (Higham 2002 - Stability of systems of linear equations).

Let $\mathbf{c}$ be a solution of a system of linear equations
$\Gamma\mathbf{c}=\mathbf{g}$, where $\mathbf{c},\mathbf{g}\in\mathbb{R}^{k}$
and $\Gamma\in\mathrm{Mat}_{k\times k}(\mathbb{R})$. Let now
$\mathbf{\tilde{c}}$ be a solution to
$(\Gamma+\Delta\Gamma)\mathbf{\tilde{c}}=\mathbf{g}+\mathbf{\Delta g}$, where
$\Delta\Gamma\in\mathrm{Mat}_{k\times k}(\mathbb{R})$, $\mathbf{\Delta
g}\in\mathbb{R}^{k}$. Assume that there exits $\delta^{\prime}>0$ such that
$\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|\leq\delta^{\prime}\left|\left|\hskip 1.0ptE\hskip
1.0pt\right|\right|$ and $|\mathbf{\Delta g}|\leq\delta^{\prime}|\mathbf{f}|$
(for some $E\in\mathrm{Mat}_{k\times k}(\mathbb{R})$,
$\mathbf{f}\in\mathbb{R}^{k}$), and that $\delta^{\prime}\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip 1.0ptE\hskip
1.0pt\right|\right|\leq 1$. Then, we have

$\displaystyle\frac{|\mathbf{c}-\mathbf{\tilde{c}}|}{|\mathbf{c}\,|}\leq\frac{\delta^{\prime}}{1-\delta^{\prime}\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip 1.0ptE\hskip
1.0pt\right|\right|}\left(\frac{\left|\left|\hskip 1.0pt\Gamma^{-1}\hskip
1.0pt\right|\right||\mathbf{f}|}{|\mathbf{c}\,|}+\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip 1.0ptE\hskip
1.0pt\right|\right|\right)\ .$ (36)

######  Lemma 1.

Under the conditions of Theorem 1 qubit states and measurements
$\\{\varrho^{x}_{a}\\}_{a,x}$, $\\{\mathbf{M}^{y}\\}_{y}$ and
$\\{\tilde{\varrho}^{x}_{a}\\}_{a,x}$, $\\{\mathbf{\tilde{M}}^{y}\\}_{y}$
satisfy

$\displaystyle\left|\left|\hskip 1.0pt\tilde{\varrho}^{x}_{a}\hskip
1.0pt\right|\right|\geq\frac{1-2\varepsilon}{1-\varepsilon},\;\forall a,x,$
$\displaystyle|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|\leq\varepsilon+\sqrt{\varepsilon},$
$\displaystyle|{\mathrm{tr}}(\tilde{M}^{y}_{b}\tilde{M}^{y^{\prime}}_{b^{\prime}})-{\mathrm{tr}}(M^{y}_{b}M^{y^{\prime}}_{b^{\prime}})|\leq\varepsilon+(1+\varepsilon)\sqrt{\varepsilon},$

$\forall x\neq x^{\prime},a\neq a^{\prime}$ and $\forall y\neq
y^{\prime},b\neq b^{\prime}$ respectively, and $\varepsilon\leq\frac{1}{3}$.

######  Lemma 2.

Let $\Gamma$ be the Gram matrix of Bloch vectors of target states
$\\{\varrho^{x}_{1}\\}_{x\in S}$. Likewise, let $\tilde{\Gamma}$ be the Gram
matrix of Bloch vectors of experimental states
$\\{\tilde{\varrho}^{x}_{1}\\}_{x\in S}$. Let
$|p(b|a,x,y)-\tilde{p}(b|a,x,y)|\leq\varepsilon$ for all $b,a,x,y$ and let
$k=|S|$ be the cardinality of the set $S$. Then, we have the following bounds
on norms of $\Delta\Gamma=\tilde{\Gamma}-\Gamma$

$\displaystyle\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|_{F}\leq F_{k}(\varepsilon)\ ,$
$\displaystyle\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|\leq O_{k}(\varepsilon)\ ,$ (37)

where

$\displaystyle
F_{k}(\varepsilon)=\sqrt{\varepsilon}\sqrt{4k(k-1)}\sqrt{1+2\sqrt{\varepsilon}+\frac{k+3}{k-1}\varepsilon}\
,$ $\displaystyle
O_{k}(\varepsilon)=2((k-1)\sqrt{\varepsilon}+(k+1)\varepsilon)\ .$

###### Proof of Theorem 3.

Let $\tilde{\Gamma}$ be the Gram matrix of Bloch vectors
$\\{\mathbf{\tilde{n}}^{x}_{1}\\}_{x\in S}$ of the experimental states for
$x\in S$. Let us for now assume that the Cholesky factorization
$\tilde{\Gamma}=\tilde{L}\tilde{L}^{T}$ exists and let $\Delta L=\tilde{L}-L$.

We start our _first part of the proof_ regarding the states $x\in S$, $a=1$ by
connecting the square of the norm $\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}$ with trace distance between the states
$\\{\tilde{\varrho}^{x}_{1}\\}_{x\in S}$ and $\\{\varrho^{x}_{1}\\}_{x\in S}$.
From the sketch of the proof of Theorem 3 it follows that the Bloch vectors
$\\{\mathbf{n}^{x}_{1}\\}_{x\in S}$ of the target states and the vectors
$\\{\mathbf{l}^{x}\\}_{x\in S}$, transpose of which are rows of $L$, are
connected via an orthogonal transformation, which we denote as $O$, i.e.,
$\mathbf{l}^{x}=O\mathbf{n}^{x}_{1}$, $x\in S$. Similarly, for the Bloch
vectors $\\{\mathbf{\tilde{n}}^{x}_{1}\\}_{x\in S}$ of the experimental states
and the vectors $\\{\mathbf{\tilde{l}}^{x}\\}_{x\in S}$ that form $\tilde{L}$
there exists an orthogonal transformation $\tilde{O}$, such that
$\mathbf{\tilde{l}}^{x}=\tilde{O}\mathbf{\tilde{n}}^{x}_{1}$, $x\in S$. Let us
define states
$\tau^{x}=\frac{1}{2}(\mathbbm{1}+\mathbf{l}^{x}\cdot\bm{\sigma})$, and
$\tilde{\tau}^{x}=\frac{1}{2}(\mathbbm{1}+\mathbf{\tilde{l}}^{x}\cdot\bm{\sigma})$,
$x\in S$. We know that there exist unitary transformations $V$ and $\tilde{V}$
such that: $\tau^{x}=V(\varrho^{x}_{1})^{(T)}V^{\dagger}$, and
$\tilde{\tau}^{x}=\tilde{V}(\tilde{\varrho}^{x}_{1})^{(T)}\tilde{V}^{\dagger}$
for $x\in S$ (the optional transposition $(\cdot)^{(T)}$ is reserved for the
cases $\det(O)=-1$ and $\det(\tilde{O})=-1$). We have the following sequence
of equalities valid for $x\in S$:

$\displaystyle\mathrm{d}_{\mathrm{tr}}(\tilde{\tau}^{x},\tau^{x})$
$\displaystyle=$
$\displaystyle\mathrm{d}_{\mathrm{tr}}(\tilde{V}(\tilde{\varrho}^{x}_{1})^{(T)}\tilde{V}^{\dagger},V(\varrho^{x}_{1})^{(T)}V^{\dagger})=\mathrm{d}_{\mathrm{tr}}(V^{\dagger}\tilde{V}(\tilde{\varrho}^{x}_{1})^{(T)}\tilde{V}^{\dagger}V,(\varrho^{x}_{1})^{(T)})$
(38) $\displaystyle=$
$\displaystyle\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger},\varrho^{x}_{1})\
,$

where we used the invariance of trace distancee under unitary evolution and
transposition, and also denoted the resulting unitary $V^{\dagger}\tilde{V}$
as $U$. It is clear that one transposition in the final formula in Eq. (38) is
enough as the case $\det(\tilde{O})=\det(O)=-1$ is equivalent to having
$\det(\tilde{O})=\det(O)=1$ and changing $U$ to $U^{T}$.

We need to show now that this unitary $U$ satisfies the claim of Theorem 2 for
the states $x\in S$. For qubits, trace distance can be expressed directly via
Eucleedian distance between Bloch vectors:

$\displaystyle\mathrm{d}_{\mathrm{tr}}(\tilde{\tau}^{x},\tau^{x})=\frac{1}{2}|\mathbf{\tilde{l}}^{x}-\mathbf{l}^{x}|\
.$ (39)

Using this fact and Eq. (38), it is straightforward to derive the following
upper bound:

$\displaystyle\frac{1}{k}\sum_{x=1}^{k}\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger},\varrho^{x}_{1})\leq\frac{1}{2\sqrt{k}}\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|_{F}\ .$ (40)

We now use Theorem [Sun 1991] (Ref. [37]) to upper-bound $\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|_{F}$. Specifically, we apply it to the
Gram matrix of the Bloch vectors of states $\\{\varrho^{x}_{1}\\}_{x\in S}$.
We have $\left|\left|\hskip 1.0ptL\hskip 1.0pt\right|\right|_{F}=\sqrt{k}$,
since the target states are assumed to be pure. The direct substitution of the
bound in Eq. (35) into Eq. (40) gives the bound in Eq. (30a), and also the
condition in Eq. (29) when the Theorem 3 (and also Theorem 2) applies.

In the beginning of the proof we assumed that the Cholesky factorization
$\tilde{\Gamma}=\tilde{L}\tilde{L}^{T}$ exists. The condition in Eq. (29)
gives the sufficient condition for this to hold. We conclude this part of the
proof by noting that a similar statement in terms of fidelity (Theorem 2, Eq.
(33a)) follows from Eq. (40) by connecting fidelity to the trace distance as
described in Eq. (54).

In _the second part of the proof_ we derive upper bounds on the trace
distances between the states corresponding to $x\in S$ and $a=2$. As we will
see, those can be connected to the bounds for the states with $x\in S$ and
$a=1$. Indeed, we can write the following for every $x\in S$:

$\displaystyle
2\mathrm{d}_{\mathrm{tr}}(U\tilde{\varrho}^{x}_{2}U^{\dagger},\varrho^{x}_{2})$
$\displaystyle=\left|\left|\hskip
1.0pt\varrho^{x}_{2}-U\tilde{\varrho}^{x}_{2}U^{\dagger}\hskip
1.0pt\right|\right|_{1}\leq\left|\left|\hskip
1.0pt\mathbbm{1}-\varrho^{x}_{1}-U(\mathbbm{1}-\tilde{\varrho}^{x}_{1})U^{\dagger}\hskip
1.0pt\right|\right|_{1}+\left|\left|\hskip
1.0ptU(\mathbbm{1}-\tilde{\varrho}^{x}_{1}-\tilde{\varrho}^{x}_{2})U^{\dagger}\hskip
1.0pt\right|\right|_{1}$ $\displaystyle\leq\left|\left|\hskip
1.0pt\varrho^{x}_{1}-U\tilde{\varrho}^{x}_{1}U^{\dagger}\hskip
1.0pt\right|\right|_{1}+\left|\left|\hskip
1.0pt\mathbbm{1}-\tilde{\varrho}^{x}_{1}+\tilde{M}^{x}_{2}\hskip
1.0pt\right|\right|_{1}-\left|\left|\hskip
1.0pt\tilde{M}^{x}_{2}-\tilde{\varrho}^{x}_{2}\hskip 1.0pt\right|\right|_{1}$
$\displaystyle=2\mathrm{d}_{\mathrm{tr}}(U\tilde{\varrho}^{x}_{1}U^{\dagger},\varrho^{x}_{1})+\sum_{a=1,2}\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|_{1}.$

Exactly the same reasoning can be applied to upper bound the trace distance
$\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{2})^{T}U^{\dagger},\varrho^{x}_{2})$,
in which case the transposition will propagate to
$\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{T}U^{\dagger},\varrho^{x}_{1})$
but would not affect the term $\sum_{a=1,2}\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|_{1}$
(as we can take $(\tilde{M}^{x}_{a})^{T}$).

To finish the calculations we need to upper-bound the second summand in Eq.
(B). We do it by writing $\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|_{1}\leq 2\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|\leq
2\sqrt{\varepsilon}$, $\forall a,x$, where the last inequality is proven in
Lemma 1 (see Eq. (77)). From here, it is easy to obtain the bound in Eq.
(30b). The reason why we do not simply apply the same reasoning to the states
corresponding to $a=2$ and $x\in S$ as for the states with $a=1$ is because we
want the isometry $U$ to be the same for all of the states.

In _the third part of the proof_ we derive the bounds for the preparations
states corresponding to $x\notin S$ and both $a=1,2$. We start by reminding
ourselves that the set of states $\\{\varrho^{x}_{1}\\}_{x\in S}$ is assumed
to be tomographically complete. Is is equivalent to the assumption that the
vectors $\\{\mathbf{l}^{x}\\}_{x\in S}$, defined above, are linearly
independent and span $\mathbb{R}^{3}$ for $k=3$, or the considered subspace
for $k=2$. The condition in Eq. (29) of the Theorem 3 (also stated in Eq.
(34)) ensures that the same holds for the vectors
$\\{\mathbf{\tilde{l}}^{x}\\}_{x\in S}$. If so, let us expand the Bloch
vectors of the states $\varrho^{x^{\prime}}_{a^{\prime}}$ and
$U\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}U^{\dagger}$, $x^{\prime}\notin S$
in terms of $\\{\mathbf{l}^{x}\\}_{x\in S}$ and
$\\{\mathbf{\tilde{l}}^{x}\\}_{x\in S}$ respectively. Let us denote the
coefficients of these linear expansions as
$\\{c^{a^{\prime},x^{\prime}}_{x}\\}_{x\in S}$ and
$\\{\tilde{c}^{a^{\prime},x^{\prime}}_{x}\\}_{x\in S}$ respectively. These
coefficients will, of course, depend on $x^{\prime}$ and $a^{\prime}$ as well
as on the choice of the set $S$. However, for simplicity of the derivations we
are going omit the subscripts $x^{\prime}$ and $a^{\prime}$ until we present
the final result. We will also assume, without loss of generality, that
$S=\\{1,2,3\\}$ (or $S=\\{1,2\\}$ for $k=2$).

It is clear that the coefficients $\\{c_{x}\\}_{x\in S}$ and
$\\{\tilde{c}_{x}\\}_{x\in S}$ satisfy the following respective systems of
linear equations:

$\displaystyle\Gamma\mathbf{c}=\mathbf{g},\quad\tilde{\Gamma}\mathbf{\tilde{c}}=\mathbf{\tilde{g}},$
(42)

where $\mathbf{c}=(c_{1},c_{2},c_{3})$,
$\mathbf{g}=(2{\mathrm{tr}}(\varrho^{x^{\prime}}_{a^{\prime}}\varrho^{1}_{1})-1,2{\mathrm{tr}}(\varrho^{x^{\prime}}_{a^{\prime}}\varrho^{2}_{1})-1,2{\mathrm{tr}}(\varrho^{x^{\prime}}_{a^{\prime}}\varrho^{3}_{1})-1)$
and analogously
$\mathbf{\tilde{c}}=(\tilde{c}_{1},\tilde{c}_{2},\tilde{c}_{3})$,
$\mathbf{g}=(2{\mathrm{tr}}(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}\tilde{\varrho}^{1}_{1})-1,2{\mathrm{tr}}(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}\tilde{\varrho}^{2}_{1})-1,2{\mathrm{tr}}(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}\tilde{\varrho}^{3}_{1})-1)$
whenever the cardinality $k$ of the set $S$ is $3$. For $k=2$ and
$\mathbf{c},\mathbf{\tilde{c}},\mathbf{g},\mathbf{\tilde{g}}\in\mathbb{R}^{2}$
their definition is analogous. We again omitted the subscripts
$a^{\prime},x^{\prime}$ for the vectors $\mathbf{g}$ and $\mathbf{\tilde{g}}$
for simplicity. The matrices $\Gamma$ and $\tilde{\Gamma}$ are still the
moment matrices for the states $x\in S$, $a=1$.

In complete analogy to Eq. (39) we can write:

$\displaystyle\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})^{(T)}U^{\dagger},\varrho^{x^{\prime}}_{a^{\prime}})=\frac{1}{2}\Big{|}\sum_{x\in
S}(\tilde{c}_{x}\mathbf{\tilde{l}}^{x}-c_{x}\mathbf{l}^{x})\Big{|}.$ (43)

We then can upper-bound the latter norm as follows:

$\displaystyle\Big{|}\sum_{x\in
S}(\tilde{c}_{x}\mathbf{\tilde{l}}^{x}-c_{x}\mathbf{l}^{x})\Big{|}\leq\Big{|}\sum_{x\in
S}c_{x}(\mathbf{\tilde{l}}^{x}-\mathbf{l}^{x})\Big{|}+\Big{|}\sum_{x\in
S}(\tilde{c}_{x}-c_{x})\mathbf{\tilde{l}}^{x}\Big{|}\leq|\mathbf{c}\,|\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|_{F}+\sum_{x\in
S}|\tilde{c}_{x}-c_{x}||\mathbf{\tilde{l}}^{x}|.$

In the equation above we used the following relation:

$\displaystyle|\sum_{x\in
S}c_{x}(\mathbf{\tilde{l}}^{x}-\mathbf{l}^{x})|=\sqrt{\sum_{i=1}^{3}\left(\sum_{x\in
S}c_{x}(\tilde{l}^{x}_{i}-l^{x}_{i})\right)^{2}}\leq\sqrt{\sum_{x^{\prime}\in
S}c^{2}_{x^{\prime}}\sum_{i=1}^{3}\sum_{x\in
S}(\tilde{l}^{x}_{i}-l^{x}_{i})^{2}}=|\mathbf{c}\,|\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|_{F}.$ (45)

Also, since the norms $|\mathbf{\tilde{l}}^{x}|$ can be upper-bounded by $1$
for all $x$, the resulting upper bound can be simplified further to be
$|\mathbf{c}\,|\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}+\sum_{x\in S}|\tilde{c}_{x}-c_{x}|$.

To estimate the deviation $|\tilde{c}_{x}-c_{x}|$ we apply Theorem [Higham
2002] by taking $\Delta\Gamma=\tilde{\Gamma}-\Gamma$ and $\mathbf{\Delta
g}=\mathbf{\tilde{g}}-\mathbf{g}$. The bound on the operator norm
$\left|\left|\hskip 1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|$ is given by
Lemma 2. At the same time $|\mathbf{\Delta g}|=2\sqrt{\sum_{x\in
S}({\mathrm{tr}}(\varrho^{x^{\prime}}_{a^{\prime}}\varrho^{x}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}}\tilde{\varrho}^{x}_{1}))^{2}}\leq
2\sqrt{k}(\sqrt{\varepsilon}+\varepsilon)$. If we take
$\delta^{\prime}=2((k-1)\sqrt{\varepsilon}+(k+1)\varepsilon)$, some matrix $E$
with $\left|\left|\hskip 1.0ptE\hskip 1.0pt\right|\right|=1$, vector
$\mathbf{f}$ such that $|\mathbf{f}|=\frac{\sqrt{k}}{k-1}$, we satisfy the
conditions of the Theorem [Higham 2002] and get the following bound:

$\displaystyle|\mathbf{c}-\mathbf{\tilde{c}}|\leq\left(|\mathbf{c}|+\frac{\sqrt{k}}{k-1}\right)\frac{\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|}{1-\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|}\ ,$ (46)

As the final step we need to connect the bound in Eq. (B) with the one in Eq.
(46) by the relation between $1$-norm and the Euclidean norm in
$\mathbb{R}^{k}$, which effectively adds a factor of $\sqrt{k}$ to the bound
in Eq. (46). Combining everything together we obtain the following:

$\displaystyle\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})^{(T)}U^{\dagger},\varrho^{x^{\prime}}_{a^{\prime}})\leq\frac{1}{2}|\mathbf{c}\,|\left|\left|\hskip
1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}+\frac{\sqrt{k}}{2}\left(|\mathbf{c}|+\frac{\sqrt{k}}{k-1}\right)\frac{\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|}{1-\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|},$ $\displaystyle x^{\prime}\notin
S,a^{\prime}=1,2.$ (47)

The above bound is given in Eq. (30c) in terms of the quantity
$E_{S,k}(\varepsilon)$ (Eq. (28)), where we took $\left|\left|\hskip
1.0pt\Delta L_{S}\hskip 1.0pt\right|\right|_{F}=2\sqrt{k}E_{S,k}(\varepsilon)$
and brought back all the necessary subscripts. Again, using the relation in
Eq. (55) we can derive bounds on the fidelity between the states
$\varrho^{x}_{a}$ and $\tilde{\varrho}^{x}_{a}$ for $x\notin S$.

In _the fourth, final, part of the proof_ we derive the bounds for the
measurements. We do it by connecting the distance between the experimental
measurements and the distance between the corresponding states. Indeed, we can
write the following:

$\displaystyle\left|\left|\hskip
1.0ptU(\tilde{M}^{y}_{1})^{(T)}U^{\dagger}-M^{y}_{1}\hskip
1.0pt\right|\right|\leq\left|\left|\hskip
1.0ptU(\tilde{\varrho}^{y}_{1})^{(T)}U^{\dagger}-\varrho^{y}_{1}\hskip
1.0pt\right|\right|+\left|\left|\hskip
1.0pt\tilde{\varrho}^{y}_{1}-\tilde{M}^{y}_{1}\hskip
1.0pt\right|\right|,\quad\forall y,$ (48)

where we used the triangle inequality and the invariance of the norm
$\left|\left|\hskip 1.0pt\tilde{\varrho}^{y}_{1}-\tilde{M}^{y}_{1}\hskip
1.0pt\right|\right|$ under unitary transformations. Remembering that
$\left|\left|\hskip 1.0pt\tilde{\varrho}^{y}_{1}-\tilde{M}^{y}_{1}\hskip
1.0pt\right|\right|\leq\sqrt{\varepsilon}$ (Lemma 1, Eq. (77)) and
$\left|\left|\hskip
1.0ptU(\tilde{\varrho}^{y}_{1})^{(T)}U^{\dagger}-\varrho^{y}_{1}\hskip
1.0pt\right|\right|=\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{y}_{1})^{(T)}U^{\dagger},\varrho^{y}_{1})$,
we produce the bounds in Eqs. (30d,30e). ∎

Finally, we proceed to the proof of Theorem 2, which gives a quantitative
robustness analysis expressed in terms of average fidelity.

###### Proof of Theorem 2.

Let us first state a useful identity between Frobenius distance and the
fidelity, valid for an arbitrary pure states $\varrho$ and a Hermitian
operator $X$:

$\left|\left|\hskip 1.0ptX-\varrho\hskip
1.0pt\right|\right|^{2}_{F}=1+{\mathrm{tr}}(X^{2})-2{\mathrm{tr}}(X\varrho)\
.$ (49)

We follow exactly the same steps that were given in the proof of Theorem 3. We
assume that the reader is familiar with the notation introduced there. First,
in order to derive the bound in Eq. (33a) we repeat the reasoning preceding
Eq. (38). Analogously to the trace distance discussed there, we use the
invariance of the fidelity under unitary transformations and transposition
which gives us for all $x\in S=[k]$:

${\mathrm{tr}}(\tilde{\tau}^{x}\tau^{x})={\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger}\varrho^{x}_{1})\
,$ (50)

where $U=V^{\dagger}\tilde{V}$. Using the above formula and standard algebra
involving Pauli matrices we obtain:

$\frac{1}{k}\sum_{x=1}^{k}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger}\varrho^{x}_{1})=\frac{1}{2}+\frac{1}{2k}\sum_{x=1}^{k}\mathbf{l}^{x}\cdot\mathbf{\tilde{l}}^{x}\
.$ (51)

On the other hand, we have the following identity:

$\displaystyle\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|^{2}_{F}=\sum_{x=1}^{k}|\mathbf{l}^{x}|^{2}+\sum_{x=1}^{k}|\mathbf{\tilde{l}}^{x}|^{2}-2\sum_{x=1}^{k}\mathbf{l}^{x}\cdot\mathbf{\tilde{l}}^{x}\
,$ (52)

where $|\cdot|$ is a standard Euclidean norm in $\mathbb{R}^{3}$. Using the
identity
$|\mathbf{\tilde{l}}^{x}|^{2}=2{\mathrm{tr}}((\tilde{\varrho}^{x}_{1})^{2})-1$
and the fact that target states are pure (which implies $|\mathbf{l}^{x}|=1$)
we obtain:

$\displaystyle\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|^{2}_{F}=2\sum_{x=1}^{k}{\mathrm{tr}}((\tilde{\varrho}^{x}_{1})^{2})-2\sum_{x=1}^{k}\mathbf{l}^{x}\cdot\mathbf{\tilde{l}}^{x}\
.$ (53)

Inserting the above into Eq. (51) and using the fact that
${\mathrm{tr}}{(\tilde{\varrho}^{x}_{1})^{2}}\geq
1-\frac{2\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}$ (this follows
straightforwardly from Lemma 2), we finally obtain:

$\frac{1}{k}\sum_{x=1}^{k}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger}\varrho^{x}_{1})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-\frac{\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|^{2}_{F}}{4k}\ .$ (54)

We complete the proof of Eq. (33a) by again employing, exactly as before,
Theorem [Sun 1991] (Ref. [37]) in order to upper-bound $\left|\left|\hskip
1.0pt\Delta L\hskip 1.0pt\right|\right|_{F}$.

Now, to derive Eq. (33b), which gives the bounds for the fidelity of states
for $x\in[k]$ and $a=2$, we can use the following inequality which can be
derived from Eq. (49) and from the relation
${\mathrm{tr}}{(\tilde{\varrho}^{x}_{1})^{2}}\geq
1-\frac{2\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}$:

$\displaystyle{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger}\varrho^{x}_{a})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-(\mathrm{d}_{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{a})^{(T)}U^{\dagger},\varrho^{x}_{a}))^{2}.$
(55)

From Eqs. (40,B) it follows that:

$\displaystyle\frac{1}{k}\sum_{x=1}^{k}(\mathrm{d}_{\mathrm{tr}}(\tilde{\varrho}^{x}_{2})^{(T)}U^{\dagger},\varrho^{x}_{2}))^{2}\leq\frac{1}{k}\sum_{x=1}^{k}(\mathrm{d}_{\mathrm{tr}}(\tilde{\varrho}^{x}_{1})^{(T)}U^{\dagger},\varrho^{x}_{1})+2\sqrt{\varepsilon})^{2}$
(56) $\displaystyle\leq\frac{\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|^{2}_{F}}{4k}+\frac{2\left|\left|\hskip 1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}}{\sqrt{k}}\sqrt{\varepsilon}+4\varepsilon=\left(\frac{\left|\left|\hskip
1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}}{2\sqrt{k}}+2\sqrt{\varepsilon}\right)^{2},$

which gives the desired bound form Eq. (33b):

$\displaystyle\frac{1}{k}\sum_{x=1}^{k}{\mathrm{tr}}(U(\tilde{\varrho}^{x}_{2})^{(T)}U^{\dagger}\varrho^{x}_{2})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-\left(\frac{\left|\left|\hskip
1.0pt\Delta L\hskip
1.0pt\right|\right|_{F}}{2\sqrt{k}}+2\sqrt{\varepsilon}\right)^{2}\ .$ (57)

The proof of the remaining bounds in Eqs. (33c,33d,33e) is straightforward and
follows directly form the formula in Eq. (49). In particular, in order to
prove Eq. (33c) we set $X=\tilde{\varrho}^{x}_{a}$ and
$\varrho=\varrho^{x}_{a}$ in Eq. (49), use the inequality
${\mathrm{tr}}{(\tilde{\varrho}^{x}_{a})^{2}}\geq
1-\frac{2\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}$ and the relation
$\left|\left|\hskip 1.0pt\sigma-\varrho\hskip
1.0pt\right|\right|^{2}_{F}=2\mathrm{d}_{\mathrm{tr}}(\varrho,\sigma)^{2}$,
with latter being true for arbitrary qubits states $\varrho$ and $\sigma$.
Derivation of Eq. (33d) is completely analogous to the one of Eq. (33b).
Finally, for Eq. (33e), the reasoning is again analogous to Eq. (33c), where,
after setting $X=\tilde{M}^{y}_{1}$ and $\varrho=\varrho^{y}_{1}$, it is
additionally necessary to use a simple lower bound:
${\mathrm{tr}}((M^{y}_{1})^{2})\geq(1-\varepsilon)^{2}$, which follows from
the conditions of Theorem 2. ∎

##  Appendix C Alternative bounds from Procrustes

######  Lemma 3 (Robust self-testing for qubits from Procrustes).

Consider pure target qubit preparation states $\varrho^{x}_{a}$ and target
projective measurements $\mathbf{M}^{y}$, where $a=1,2$ and $x,y\in[n]$.
Assume that $\varrho^{x}_{a}=M^{x}_{a}$ for all $a,x$ and, furthermore, that
experimental states $\tilde{\varrho}^{x}_{a}$ and measurements
$\mathbf{\tilde{M}}^{y}$ act on Hilbert space of dimension at most $d$ and
generate statistics
$\tilde{p}(b|a,x,y)={\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{y}_{b})$
such that
$|\tilde{p}(b|a,x,y)-{\mathrm{tr}}(\varrho^{x}_{a}M^{y}_{b})|\leq\varepsilon$,
for all $a,b,x,y$.

Let $\\{\varrho_{i}\\}_{i=1}^{m}$ be a subset of $m$ considered states among
$\varrho^{x}_{a}$. Let $L$ be a matrix whose rows are the Bloch vectors of
states $\varrho_{i}$, $i\in[m]$, and let $k\in\\{{2,3\\}}$ be its rank ($m\geq
k$). Assume, without loss of generality, that for $k=2$ the third component of
the Bloch vectors is $0$. In that case, truncate $L$ to the first two columns.
Let us define two auxiliary functions:

$\displaystyle P_{m}(\varepsilon,L)=\begin{cases}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)+\min\left[\frac{\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)}{\sqrt{1-\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|^{2}F_{m}(\varepsilon)}},\sqrt{\sqrt{k}F_{m}(\varepsilon)}\right],&\text{if}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip 1.0pt\right|\right|\sqrt{F_{m}(\varepsilon)}<1\\\
\left|\left|\hskip 1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)+\sqrt{\sqrt{k}F_{m}(\varepsilon)}&\text{otherwise,}\end{cases}$
(58)

and

$\displaystyle
F_{m}(\varepsilon)=\sqrt{4m(m-1)\varepsilon\left(1+2\sqrt{\varepsilon}+\frac{m+3}{m-1}\varepsilon\right)},$
(59)

where $L^{\ddagger}=(LL^{T})^{-1}L^{T}$. There exist a unitary matrix $U$ such
that:

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}{\mathrm{tr}}(\varrho_{i}U\tilde{\varrho}_{i}U^{\dagger})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-\frac{1}{4m}P^{2}_{m}(\varepsilon,L).$
(60)

In this section we derive alternative bounds for the fidelity between
preparation states that follows from the bounds on the so-called orthogonal
Procrustes problem [47]. The problem itself can be formulated as follows.
Given two sets of vectors $\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{m}$
and $\mathbf{y}_{1},\mathbf{y}_{2},\dots,\mathbf{y}_{m}$ in $\mathbb{R}^{d}$
find an orthogonal transformation $O\in\mathrm{O}(d)$ in $\mathbb{R}^{d}$ that
minimizes $\sum_{i=1}^{m}|\mathbf{x}_{i}-O\mathbf{y}_{i}|$. This problem has a
clear relevance to our task. Indeed, if we take $\mathbf{x}_{i}$ to be the
Bloch vectors of the target qubit preparation states and $\mathbf{y}_{i}$ the
Bloch vectors of the experimental states, then minimization over
$\mathrm{O}(3)$ is the same as the problem of finding a unitary transformation
that connects those qubit states. In Ref. [38] the bounds on Procrustes
problem were derived. We give formulation of Theorem 1 from Ref. [38] below,
where we change the notation according to our problem.

######  Theorem (Arias-Castro et.al. 2020 - A perturbation bound for
Procrustes).

Given two tall matrices $L$ and $\tilde{L}$ of same size with $L$ having full
rank, and set $\delta^{2}=\left|\left|\hskip
1.0pt\tilde{L}\tilde{L}^{T}-LL^{T}\hskip 1.0pt\right|\right|_{F}$. Then we
have

$\displaystyle\min_{O\in\mathrm{O}(k)}\left|\left|\hskip
1.0ptL-O\tilde{L}\hskip
1.0pt\right|\right|_{F}\leq\begin{cases}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|\delta^{2}+\min\left[\frac{\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|\delta^{2}}{\sqrt{1-\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|^{2}\delta^{2}}},k^{\frac{1}{4}}\delta\right],&\text{if}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip 1.0pt\right|\right|\delta<1\\\ \left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|\delta^{2}+k^{\frac{1}{4}}\delta&\text{otherwise.}\end{cases}$
(61)

In the formulation of the above theorem $L^{\ddagger}$ stands for the
Moore–Penrose inverse, which can be defined as
$L^{\ddagger}=(LL^{T})^{-1}L^{T}$ for tall matrices of full rank.

###### Proof of Lemma 3.

. Unitizing the results of the Theorem [Arias-Castro et.al. 2020] is rather
straightforward. Let $L$ be a matrix which rows are the Bloch vectors of all
target preparation states $\\{\varrho^{x}_{a}\\}_{a,x}$. Let $\tilde{L}$ in
turn be the matrix of Bloch vectors of the experimental states
$\\{\tilde{\varrho}^{x}_{a}\\}_{a,x}$. The matrices $LL^{T}$ and
$\tilde{L}\tilde{L}^{T}$ are then, of course, the full Gram matrices $\Gamma$
and $\tilde{\Gamma}$. By “full" we mean that now we do not select a subset $S$
of linearly independent vectors among the Bloch vectors of $\varrho^{x}_{a}$.
Instead, $\Gamma$ and $\tilde{\Gamma}$ are formed by all the considered
states, which can still be a subset of the $2n$ states ($x\in[n]$, $a\in[2]$).
We will be using $m$ to denote the number of the considered states, and a
simple one-indexed set $\\{\varrho_{i}\\}_{i}$ to denote the states
themselves. Estimating $\delta^{2}$ from the formulation of the Theorem
[Arias-Castro et.al. 2020] is a direct application of the bound on
$\left|\left|\hskip 1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|_{F}$ from
Lemma 2, where now instead of $k$ one should put the number $m$ of the
considered preparation states.

As for the left-hand side of Eq. (61), we can write the following:

$\displaystyle\left|\left|\hskip 1.0ptL-O\tilde{L}\hskip
1.0pt\right|\right|_{F}^{2}={\mathrm{tr}}(LL^{T})+{\mathrm{tr}}(\tilde{L}\tilde{L}^{T})-2{\mathrm{tr}}(L^{T}O\tilde{L})$
$\displaystyle\geq
m+m-m\frac{4\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-4\sum_{i=1}^{m}{\mathrm{tr}}(\varrho_{i}U\tilde{\varrho}_{i}U^{\dagger})+2m$
(62)

where we use the following identity:

$\displaystyle{\mathrm{tr}}(L^{T}O\tilde{L})=\sum_{i=1}^{m}\mathbf{n}_{i}\cdot
O\mathbf{\tilde{n}}_{i}=2\sum_{i=1}^{m}{\mathrm{tr}}(\varrho_{i}U\tilde{\varrho}_{i}U^{\dagger})-m,$
(63)

and $U$ is the unitary transformation in $\mathrm{SU}(2)$ corresponding to the
orthogonal transformation $O$ of the Bloch vectors. The bound on the average
fidelity between $m$ target preparation states and the corresponding
experimental states is then simply:

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}{\mathrm{tr}}(\varrho_{i}U\tilde{\varrho}_{i}U^{\dagger})\geq
1-\frac{\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}-\frac{1}{4m}P^{2}_{m}(\varepsilon,L)$
(64)

where

$\displaystyle P_{m}(\varepsilon,L)=\begin{cases}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)+\min\left[\frac{\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)}{\sqrt{1-\left|\left|\hskip
1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|^{2}F_{m}(\varepsilon)}},\sqrt{\sqrt{k}F_{m}(\varepsilon)}\right],&\text{if}\left|\left|\hskip
1.0ptL^{\ddagger}\hskip 1.0pt\right|\right|\sqrt{F_{m}(\varepsilon)}<1\\\
\left|\left|\hskip 1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|F_{m}(\varepsilon)+\sqrt{\sqrt{k}F_{m}(\varepsilon)}&\text{otherwise,}\end{cases}$
(65)

and

$\displaystyle
F_{m}(\varepsilon)=\sqrt{4m(m-1)\varepsilon\left(1+2\sqrt{\varepsilon}+\frac{m+3}{m-1}\varepsilon\right)}.$
(66)

∎

##  Appendix D Examples

Below we provide detailed derivations of the results presented in Table 1.

The first example concerns $n=2,3$ MUBs in $d=2$. Since for MUBs
${\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})=\frac{1}{2}$,
for $x\neq x^{\prime},\forall a,a^{\prime}$, it follows that $\Gamma$ is an
identity matrix in $\mathbb{R}^{n}$ ($n\in\\{2,3\\}$). Hence, in Theorem 2
(see Appendix B) we should take $\left|\left|\hskip 1.0pt\Gamma_{S}\hskip
1.0pt\right|\right|=\left|\left|\hskip 1.0pt\Gamma^{-1}_{S}\hskip
1.0pt\right|\right|=\left|\left|\hskip 1.0ptL_{S}\hskip
1.0pt\right|\right|=1$, and $\left|\left|\hskip 1.0pt\Gamma\hskip
1.0pt\right|\right|_{F}=\sqrt{n}$. The resulting bound is the average between
expressions in Eq. (30a) and Eq. (30b) with the function
$E_{S,k}(\varepsilon)$ being simply:

$\displaystyle
E_{S,k}(\varepsilon)=\frac{1}{2\sqrt{2n}}\frac{F_{k}(\varepsilon)}{\sqrt{1-O_{k}(\varepsilon)}}.$
(67)

The leading linear term is given in Table 1 for both $n=2,3$.

The second example is a little less straightforward. From the condition
${\mathrm{tr}}(\varrho^{1}_{1}\varrho^{2}_{1})=\frac{1+\alpha}{2}$,
$\alpha\in(-1,1)$ we obtain that $\Gamma=\left(\begin{smallmatrix}1&\alpha\\\
\alpha&1\end{smallmatrix}\right)$, and hence $\left|\left|\hskip
1.0pt\Gamma\hskip 1.0pt\right|\right|=1+|\alpha|$, $\left|\left|\hskip
1.0pt\Gamma^{-1}\hskip 1.0pt\right|\right|=\frac{1}{1-|\alpha|}$, and
$\left|\left|\hskip 1.0pt\Gamma\hskip
1.0pt\right|\right|_{F}=\sqrt{2+2\alpha^{2}}$. The output $L$ of the Cholesky
factorization is $L=\left(\begin{smallmatrix}1&0\\\
\alpha&\sqrt{1-\alpha^{2}}\end{smallmatrix}\right)$, which leads to
$\left|\left|\hskip 1.0ptL\hskip 1.0pt\right|\right|=\sqrt{1+|a|}$. This also
determines the minimum in Eq. (28) to be $\frac{\sqrt{1+|a|}}{\sqrt{2}}$.
Plugging this values in Eq. (28) gives:

$\displaystyle
E_{\\{1,2\\},2}(\varepsilon)=\frac{1}{4}\frac{\sqrt{1+|\alpha|}F_{k}(\varepsilon)}{\sqrt{1-|\alpha|-O_{k}(\varepsilon)}}.$
(68)

The final bound is again the average of the bounds in Eq. (30a) and Eq. (30b).
The first order in $\varepsilon$ for this bound is given in Table 1. The
applicability of the above bound is determined by the inequality
$1-|\alpha|-O_{2}(\varepsilon)\geq 0$. Importantly, the latter condition gives
nonempty region $\varepsilon\in[0,\varepsilon_{0})$ whenever $|\alpha|>0$.

The third example is a trine ensemble of states
$(\varrho^{1}_{1},\varrho^{2}_{2},\varrho^{3}_{1})$, with
$\varrho^{x}_{1}=\frac{\mathbbm{1}}{2}+\frac{1}{2}\mathbf{n}_{x}\cdot\bm{\sigma}$,
$x=1,2,3$, and where $\mathbf{n}_{1}=(1,0,0)$,
$\mathbf{n}_{2}=\left(-\frac{1}{2},\frac{\sqrt{3}}{2},0\right)$, and
$\mathbf{n}_{3}=\left(-\frac{1}{2},-\frac{\sqrt{3}}{2},0\right)$. For this
configuration of the preparation states the alternative robustness analysis
via Procrustes (see Appendix C) gives better bounds. Given the vectors
$\mathbf{n}_{i}$, $i=1,2,3$ we can directly compute that $\left|\left|\hskip
1.0ptL^{\ddagger}\hskip 1.0pt\right|\right|=\sqrt{\frac{2}{3}}$. Inserting
this value to Lemma 3 produces the results given in Table 1.

The fourth example is the tetrahedron, with $\mathbf{n}_{1}=(0,0,1)$,
$\mathbf{n}_{2}=\left(\sqrt{\frac{8}{9}},0,-\frac{1}{3}\right)$,
$\mathbf{n}_{3}=\left(-\sqrt{\frac{2}{9}},\sqrt{\frac{2}{3}},-\frac{1}{3}\right)$,
$\mathbf{n}_{4}=\left(-\sqrt{\frac{2}{9}},-\sqrt{\frac{2}{3}},-\frac{1}{3}\right)$,
and
$\varrho^{x}_{1}=\frac{\mathbbm{1}}{2}+\frac{1}{2}\mathbf{n}_{x}\cdot\bm{\sigma}$,
$x=1,2,3,4$ as before. In this case, we also employ the bounds from
Procrustes. For the above configuration of states we have that
$\left|\left|\hskip 1.0ptL^{\ddagger}\hskip
1.0pt\right|\right|=\frac{\sqrt{3}}{2}$, which leads to the results in Table
1.

##  Appendix E Proofs of auxiliary results

###### Proof of Lemma 1.

First of all, we improve the bound on the norm of each of the experimental
states. From $\left|\left|\hskip 1.0pt\tilde{M}^{y}_{b}\hskip
1.0pt\right|\right|\geq 1-\varepsilon$, for $d=2$ it follows immediately that
the second (second largest) eigenvalue of each of the effects
$\tilde{M}^{y}_{b}$ cannot exceed $\varepsilon$. Hence, we can conclude that
${\mathrm{tr}}(\tilde{M}^{y}_{b})\leq 1+\varepsilon,\forall y,b$.

Secondly, we can improve the bound on the norm of each of the experimental
states. For that let us write the spectral decomposition of each experimental
state and POVM effect:

$\displaystyle\tilde{\varrho}^{x}_{a}$
$\displaystyle=\eta(\mathbbm{1}-|\psi\rangle\langle\psi|)+(1-\eta)|\psi\rangle\langle\psi|=\eta\mathbbm{1}+(1-2\eta)|\psi\rangle\langle\psi|$
(69) $\displaystyle\tilde{M}^{x}_{a}$
$\displaystyle=\lambda_{0}|\phi\rangle\langle\phi|+\lambda_{1}(\mathbbm{1}-|\phi\rangle\langle\phi|),$

where we assume that $\lambda_{0}\geq\lambda_{1}$. We omitted the indices
$x,a$ for $\eta$, $\lambda_{0},\lambda_{1}$ and $\psi$,$\phi$ for simplicity.
We can assume, without loss of generality, that $\eta\geq\frac{1}{2}$, and let
us also assume for now that $|\langle\phi|\psi\rangle|^{2}\leq\frac{1}{2}$.
From the condition
${\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{M}^{x}_{a})\geq 1-\varepsilon$,
it then follows that:

$\displaystyle\eta(\lambda_{0}+\lambda_{1})+(1-2\eta)(\lambda_{0}|\langle\phi|\psi\rangle|^{2}+\lambda_{1}(1-|\langle\phi|\psi\rangle|^{2}))\geq
1-\varepsilon,$ (70)

from where we can obtain a lower bound on $\eta$, namely:

$\displaystyle\eta\geq\frac{1}{1-2|\langle\phi|\psi\rangle|^{2}}\left(\frac{1-\varepsilon-\lambda_{1}}{\lambda_{0}-\lambda_{1}}-\frac{1}{2}\right)+\frac{1}{2}.$
(71)

The expression on the right-hand side of the above inequality is maximal when
$|\langle\phi|\psi\rangle|^{2}=0$ and $\lambda_{0}=1$,
$\lambda_{1}=\varepsilon$, which returns the bound $\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip
1.0pt\right|\right|\geq\frac{1-2\varepsilon}{1-\varepsilon}$.

Now let us return to our assumption
$|\langle\phi|\psi\rangle|^{2}\leq\frac{1}{2}$, for which the above bound is
valid. We can upper-bound $\eta$ by $1$ in Eq. (70), which returns nontrivial
upper bound on $|\langle\phi|\psi\rangle|^{2}$ that happens to be
$1-\frac{1-2\varepsilon}{1-\varepsilon}$. This function is below $\frac{1}{2}$
for $\varepsilon\leq\frac{1}{3}$, i.e., for $\varepsilon\leq\frac{1}{3}$ our
newly-derived bound on $\left|\left|\hskip 1.0pt\tilde{\varrho}^{x}_{a}\hskip
1.0pt\right|\right|$ is valid. The region $\varepsilon\in[0,\frac{1}{3}]$ is
significantly larger than the resulting region in which our self-testing
argument are valid, so this assumption does not affect our final results.

Let us now try to improve the bounds for the overlaps. In the proof of Theorem
1 (see Eq. (3)) we have already established that:

$\displaystyle|{\mathrm{tr}}(\tilde{\varrho}^{x}_{a}\tilde{\varrho}^{x^{\prime}}_{a^{\prime}})-{\mathrm{tr}}(\varrho^{x}_{a}\varrho^{x^{\prime}}_{a^{\prime}})|\leq\varepsilon+\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|,\quad\forall x\neq x^{\prime},\forall a,a^{\prime}.$ (72)

Let us now refine the bound on $\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip 1.0pt\right|\right|$. For
simplicity, we present this result below for a pair of operators
$\varrho=(1-\eta)\mathbbm{1}+(2\eta-1)|\psi\rangle\langle\psi|$ and
$M=\lambda_{1}\mathbbm{1}+(\lambda_{0}-\lambda_{1})|\phi\rangle\langle\phi|$,
that satisfy the conditions ${\mathrm{tr}}(\varrho M)\geq 1-\varepsilon$, and
$1-\varepsilon\leq\lambda_{0}\leq 1$, $0\leq\lambda_{1}\leq\varepsilon$. To
compute the norm, we look for the eigenvector $|\xi\rangle$ of the operator
$(\varrho-M)$, i.e., $(\varrho-M)|\xi\rangle=\Lambda|\xi\rangle$,
corresponding to the maximal eigenvalue. We have the following quadratic
equation for the eigenvalue $\Lambda$:

$\displaystyle\Lambda^{2}-\Lambda(1-\lambda_{0}-\lambda_{1})+\eta(1-\eta)-\lambda_{1}+\lambda_{0}\lambda_{1}-\eta(\lambda_{0}-\lambda_{1})+(2\eta-1)(\lambda_{0}-\lambda_{1})|\langle\phi|\psi\rangle|^{2}=0.$

The sum of the roots of this equation is equal to $1-\lambda_{0}-\lambda_{1}$.
Since we know that $|1-\lambda_{0}-\lambda_{1}|\leq\varepsilon$, then either
both roots are of the same sign, in which case the largest eigenvalue could
only be $\varepsilon$, or they are of the opposite sign. In the latter case,
it is evident that the absolute values of both roots are maximal whenever the
free term in Eq. (E) is minimal. We then can upper-bound the solutions to Eq.
(E) by lower-bounding the free term using the condition ${\mathrm{tr}}(\varrho
M)\geq 1-\varepsilon$. Indeed, from ${\mathrm{tr}}(\varrho M)\geq
1-\varepsilon$ we infer immediately that:

$\displaystyle{\mathrm{tr}}(\varrho
M)=\lambda_{0}-\eta(\lambda_{0}-\lambda_{1})+(2\eta-1)(\lambda_{0}-\lambda_{1})|\langle\phi|\psi\rangle|^{2}\geq
1-\varepsilon,$ (74)

and thus we reduce Eq. (E) to:

$\displaystyle\Lambda^{2}-\Lambda(1-\lambda_{0}-\lambda_{1})+\eta(1-\eta)-\lambda_{0}-\lambda_{1}+\lambda_{0}\lambda_{1}+1-\varepsilon=0.$
(75)

Using the same argument we can set $\eta=1$, because $\eta(1-\eta)\geq 0$. The
positive root of Eq. (75) is equal to:

$\displaystyle\Lambda=\frac{1-\lambda_{0}-\lambda_{1}}{2}+\sqrt{\left(\frac{1-\lambda_{0}-\lambda_{1}}{2}\right)^{2}-(1-\lambda_{0})(1-\lambda_{1})+\varepsilon}.$
(76)

It is easy to check that the above expression does not have any local maxima
w.r.t. $\lambda_{0}$,$\lambda_{1}$ on the domain
$1-\varepsilon\leq\lambda_{0}\leq 1$, $0\leq\lambda_{1}\leq\varepsilon$,
whenever $\varepsilon<\frac{1}{2}$, which we assume to be the case. Thus, we
conclude that the maximal value of $\Lambda$ corresponds to the boundary of
the region of $(\lambda_{0},\lambda_{1})$. By considering this boundary we
find that this maximal value corresponds to the case of $\lambda_{0}=1$ and
$\lambda_{1}=0$ which yields $\Lambda=\sqrt{\varepsilon}$. From the above
argument we finally conclude that:

$\displaystyle\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}-\tilde{M}^{x}_{a}\hskip
1.0pt\right|\right|\leq\sqrt{\varepsilon},\quad\forall a,$ (77)

which completes the proof for the state overlaps. From Eq. (77) and the fact
that ${\mathrm{tr}}(\tilde{M}^{y}_{b})\leq 1+\varepsilon$, it is easy to
obtain the improved bound on the overlaps between measurement effect. ∎

###### Proof of Lemma 2.

Let us start by deriving the bound on the Frobenius norm $\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|_{F}$:

$\displaystyle\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|^{2}_{F}=\sum_{x=1}^{k}(|\Gamma_{x,x}-\tilde{\Gamma}_{x,x}|^{2}+\sum_{x^{\prime}\neq
x}|\Gamma_{x,x^{\prime}}-\tilde{\Gamma}_{x,x^{\prime}}|^{2})$
$\displaystyle=4\sum_{x=1}^{k}\Big{(}(1-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1})^{2})^{2}+\sum_{x^{\prime}\neq
x}|{\mathrm{tr}}(\varrho^{x}_{1}\varrho^{x^{\prime}}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1}\tilde{\varrho}^{x^{\prime}}_{1})|^{2}\Big{)}.$
(78)

We have already established the bound on
$|{\mathrm{tr}}(\varrho^{x}_{1}\varrho^{x^{\prime}}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1}\tilde{\varrho}^{x^{\prime}}_{1})|$
in Lemma 1. The bound on $(1-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1})^{2})^{2}$
can be obtained from the bound on the norm of each $\tilde{\varrho}^{x}_{1}$.
Namely, from the condition $\left|\left|\hskip
1.0pt\tilde{\varrho}^{x}_{a}\hskip
1.0pt\right|\right|\geq\frac{1-2\varepsilon}{1-\varepsilon}$ we can
immediately conclude that:

$\displaystyle
1-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1})^{2}\leq\frac{2\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}.$
(79)

From here, it is easy to get to the final bound on $\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|^{2}_{F}$, which reads:

$\displaystyle\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|^{2}_{F}\leq
4k(k-1)(\varepsilon+\sqrt{\varepsilon})^{2}+\frac{16k\varepsilon^{2}(1-2\varepsilon)^{2}}{(1-\varepsilon)^{4}}\leq
4k(k-1)\varepsilon\left(1+2\sqrt{\varepsilon}+\frac{k+3}{k-1}\varepsilon\right),$
(80)

where we made some approximations to simplify the result.

Now, let us derive the bound on $\left|\left|\hskip 1.0pt\Delta\Gamma\hskip
1.0pt\right|\right|$. In principle, we know that $\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|\leq\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|_{F}$, but we can derive a better
bound based on the fact that the diagonal entries of $\Delta\Gamma$ are much
less than the off-diagonal entries.

First of all, due to the triangle inequality, we can write $\left|\left|\hskip
1.0pt\Delta\Gamma\hskip 1.0pt\right|\right|\leq\left|\left|\hskip
1.0pt\mathrm{diag}(\Delta\Gamma)\hskip 1.0pt\right|\right|+\left|\left|\hskip
1.0pt\mathrm{offdiag}(\Delta\Gamma)\hskip 1.0pt\right|\right|$, where we split
$\Delta\Gamma$ on the diagonal and off-diagonal parts. The first term,
$\left|\left|\hskip 1.0pt\mathrm{diag}(\Delta\Gamma)\hskip
1.0pt\right|\right|$, can be easily bounded as follows:

$\displaystyle\left|\left|\hskip 1.0pt\mathrm{diag}(\Delta\Gamma)\hskip
1.0pt\right|\right|=2\max_{x}(1-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1})^{2})\leq\frac{4\varepsilon(1-2\varepsilon)}{(1-\varepsilon)^{2}}\leq
4\varepsilon.$ (81)

As for the off-diagonal part, we give the proof for two cases $k=2$ and $k=3$
separately. For $k=2$, $\left|\left|\hskip
1.0pt\mathrm{offdiag}(\Delta\Gamma)\hskip
1.0pt\right|\right|=2|{\mathrm{tr}}(\varrho^{1}_{1}\varrho^{2}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{1}_{1}\tilde{\varrho}^{2}_{1})|\leq
2\sqrt{\varepsilon}+2\varepsilon$, where we used the results of Lemma 1. As
for $k=3$, we will need some intermediate result, namely the following
relation:

$\displaystyle\left|\left|\hskip 1.0ptA\hskip
1.0pt\right|\right|\leq\sqrt{\frac{k-1}{k}}\left|\left|\hskip 1.0ptA\hskip
1.0pt\right|\right|_{F},$ (82)

where $k$ is the size of the matrix $A$ with ${\mathrm{tr}}(A)=0$. We give a
proof of this below.

Let us assume that $\\{\lambda_{i}\\}_{i=1}^{k}$ are the eigenvalues of matrix
$A$, hence we know that $\sum_{i=1}^{k}\lambda_{i}=0.$ Let $\lambda_{1}$ be
the largest eigenvalue, i.e., the operator norm of $A$, if $\lambda_{1}\geq
0$. If we wish to maximize the Frobenius norm of $A$ for fixed $\lambda_{1}$,
the following lower-bound has to be satisfied:

$\displaystyle\left|\left|\hskip 1.0ptA\hskip
1.0pt\right|\right|_{F}^{2}=\lambda_{1}^{2}+\sum_{i=2}^{k}\lambda_{i}^{2}\geq\lambda_{1}^{2}+\frac{1}{k-1}\left(\sum_{i=2}^{k}|\lambda_{i}|\right)^{2}\geq\lambda_{1}^{2}+\frac{1}{k-1}\left(\sum_{i=2}^{k}\lambda_{i}\right)^{2}=\lambda_{1}^{2}\frac{k}{k-1},$
(83)

which proves the bound in Eq. (82). Using the above result, we obtain the
following bound:

$\displaystyle\left|\left|\hskip 1.0pt\mathrm{offdiag}(\Delta\Gamma)\hskip
1.0pt\right|\right|\leq\sqrt{\frac{2}{3}}\left|\left|\hskip
1.0pt\mathrm{offdiag}(\Delta\Gamma)\hskip
1.0pt\right|\right|_{F}=2\sqrt{\frac{2}{3}}\sqrt{\sum_{x\neq
x^{\prime}}|{\mathrm{tr}}(\varrho^{x}_{1}\varrho^{x^{\prime}}_{1})-{\mathrm{tr}}(\tilde{\varrho}^{x}_{1}\tilde{\varrho}^{x^{\prime}}_{1})|^{2}}\leq
4(\sqrt{\varepsilon}+\varepsilon),$

which completes out proof. ∎

## References

  * Arute et al. [2019] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando G. S. L. Brandao, David A. Buell, Brian Burkett, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins, William Courtney, Andrew Dunsworth, Edward Farhi, Brooks Foxen, Austin Fowler, Craig Gidney, Marissa Giustina, Rob Graff, Keith Guerin, Steve Habegger, Matthew P. Harrigan, Michael J. Hartmann, Alan Ho, Markus Hoffmann, Trent Huang, Travis S. Humble, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang, Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov, Sergey Knysh, Alexander Korotkov, Fedor Kostritsa, David Land huis, Mike Lindmark, Erik Lucero, Dmitry Lyakh, Salvatore Mandrà, Jarrod R. McClean, Matthew McEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen, Masoud Mohseni, Josh Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Murphy Yuezhen Niu, Eric Ostby, Andre Petukhov, John C. Platt, Chris Quintana, Eleanor G. Rieffel, Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim Smelyanskiy, Kevin J. Sung, Matthew D. Trevithick, Amit Vainsencher, Benjamin Villalonga, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Hartmut Neven, and John M. Martinis.  Quantum supremacy using a programmable superconducting processor.  _Nature_ , 574(7779):505–510, Oct 2019.  doi: 10.1038/s41586-019-1666-5.
  * Degen et al. [2017] C. L. Degen, F. Reinhard, and P. Cappellaro.  Quantum sensing.  _Rev. Mod. Phys._ , 89:035002, Jul 2017.  doi: 10.1103/RevModPhys.89.035002.  URL https://link.aps.org/doi/10.1103/RevModPhys.89.035002.
  * Georgescu et al. [2014] I. M. Georgescu, S. Ashhab, and Franco Nori.  Quantum simulation.  _Rev. Mod. Phys._ , 86:153–185, Mar 2014.  doi: 10.1103/RevModPhys.86.153.  URL https://link.aps.org/doi/10.1103/RevModPhys.86.153.
  * Montanaro [2016] Ashley Montanaro.  Quantum algorithms: an overview.  _npj Quantum Information_ , 2:15023, Jan 2016.  doi: 10.1038/npjqi.2015.23.
  * Biamonte et al. [2017] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  _Nature_ , 549(7671):195–202, Sep 2017.  doi: 10.1038/nature23474.
  * Dunjko and Briegel [2018] Vedran Dunjko and Hans J Briegel.  Machine learning & artificial intelligence in the quantum domain: a review of recent progress.  _Reports on Progress in Physics_ , 81(7):074001, jun 2018.  doi: 10.1088/1361-6633/aab406.  URL https://doi.org/10.1088%2F1361-6633%2Faab406.
  * 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_ , art. arXiv:1910.06343, Oct 2019.
  * Mayers and Yao [2003] Dominic Mayers and Andrew Yao.  Self testing quantum apparatus.  _arXiv preprint quant-ph/0307205_ , 2003.
  * Šupić and Bowles [2019] Ivan Šupić and Joseph Bowles.  Self-testing of quantum systems: a review.  _arXiv preprint arXiv:1904.10042_ , 2019.
  * Bell [1964] John S Bell.  On the einstein podolsky rosen paradox.  _Physics Physique Fizika_ , 1(3):195, 1964.
  * Chen et al. [2016] Shin-Liang Chen, Costantino Budroni, Yeong-Cherng Liang, and Yueh-Nan Chen.  Natural framework for device-independent quantification of quantum steerability, measurement incompatibility, and self-testing.  _Physical Review Letters_ , 116(24):240401, 2016\.
  * Coladangelo et al. [2017] Andrea Coladangelo, Koon Tong Goh, and Valerio Scarani.  All pure bipartite entangled states can be self-tested.  _Nature communications_ , 8:15485, 2017.
  * Bowles et al. [2018] Joseph Bowles, Ivan Šupić, Daniel Cavalcanti, and Antonio Acín.  Device-independent entanglement certification of all entangled states.  _Physical Review Letters_ , 121(18):180503, 2018\.
  * Popescu and Rohrlich [1992] Sandu Popescu and Daniel Rohrlich.  Which states violate Bell’s inequality maximally?  _Physics Letters A_ , 169(6):411–414, 1992.
  * Ahrens et al. [2014] Johan Ahrens, Piotr Badziąg, Marcin Pawłowski, Marek Żukowski, and Mohamed Bourennane.  Experimental tests of classical and quantum dimensionality.  _Phys. Rev. Lett._ , 112:140401, Apr 2014.  doi: 10.1103/PhysRevLett.112.140401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.112.140401.
  * Brask et al. [2017] Jonatan Bohr Brask, Anthony Martin, William Esposito, Raphael Houlmann, Joseph Bowles, Hugo Zbinden, and Nicolas Brunner.  Megahertz-rate semi-device-independent quantum random number generators based on unambiguous state discrimination.  _Phys. Rev. Applied_ , 7:054018, May 2017.  doi: 10.1103/PhysRevApplied.7.054018.  URL https://link.aps.org/doi/10.1103/PhysRevApplied.7.054018.
  * Aguilar et al. [2018] Edgar A. Aguilar, Máté Farkas, Daniel Martínez, Matías Alvarado, Jaime Cariñe, Guilherme B. Xavier, Johanna F. Barra, Gustavo Cañas, Marcin Pawłowski, and Gustavo Lima.  Certifying an irreducible 1024-dimensional photonic state using refined dimension witnesses.  _Phys. Rev. Lett._ , 120:230503, Jun 2018.  doi: 10.1103/PhysRevLett.120.230503.  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.230503.
  * Anwer et al. [2020] Hammad Anwer, Sadiq Muhammad, Walid Cherifi, Nikolai Miklin, Armin Tavakoli, and Mohamed Bourennane.  Experimental characterization of unsharp qubit observables and sequential measurement incompatibility via quantum random access codes, Aug 2020\.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.080403.
  * Pawłowski and Brunner [2011] Marcin Pawłowski and Nicolas Brunner.  Semi-device-independent security of one-way quantum key distribution.  _Phys. Rev. A_ , 84:010302, Jul 2011.  doi: 10.1103/PhysRevA.84.010302.  URL https://link.aps.org/doi/10.1103/PhysRevA.84.010302.
  * Van Himbeeck et al. [2017] Thomas Van Himbeeck, Erik Woodhead, Nicolas J. Cerf, Raúl García-Patrón, and Stefano Pironio.  Semi-device-independent framework based on natural physical assumptions.  _Quantum_ , 1:33, November 2017.  ISSN 2521-327X.  doi: 10.22331/q-2017-11-18-33.  URL https://doi.org/10.22331/q-2017-11-18-33.
  * Chaves et al. [2015] Rafael Chaves, Jonatan Bohr Brask, and Nicolas Brunner.  Device-independent tests of entropy.  _Phys. Rev. Lett._ , 115:110501, Sep 2015.  doi: 10.1103/PhysRevLett.115.110501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.115.110501.
  * Tavakoli et al. [2018a] Armin Tavakoli, Jędrzej Kaniewski, Tamás Vértesi, Denis Rosset, and Nicolas Brunner.  Self-testing quantum states and measurements in the prepare-and-measure scenario.  _Physical Review A_ , 98(6):062307, 2018a.
  * Farkas and Kaniewski [2019] Máté Farkas and Jędrzej Kaniewski.  Self-testing mutually unbiased bases in the prepare-and-measure scenario.  _Physical Review A_ , 99(3):032316, 2019.
  * Mironowicz and Pawłowski [2019] Piotr Mironowicz and Marcin Pawłowski.  Experimentally feasible semi-device-independent certification of four-outcome positive-operator-valued measurements.  _Phys. Rev. A_ , 100:030301, Sep 2019.  doi: 10.1103/PhysRevA.100.030301.  URL https://link.aps.org/doi/10.1103/PhysRevA.100.030301.
  * Tavakoli et al. [2018b] Armin Tavakoli, Massimiliano Smania, Tamás Vértesi, Nicolas Brunner, and Mohamed Bourennane.  Self-testing non-projective quantum measurements.  _arXiv preprint arXiv:1811.12712_ , 2018b.
  * Tavakoli et al. [2019] Armin Tavakoli, Denis Rosset, and Marc-Olivier Renou.  Enabling computation of correlation bounds for finite-dimensional quantum systems via symmetrization.  _Phys. Rev. Lett._ , 122:070501, Feb 2019.  doi: 10.1103/PhysRevLett.122.070501.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.070501.
  * Blume-Kohout et al. [2017] Robin Blume-Kohout, John King Gamble, Erik Nielsen, Kenneth Rudinger, Jonathan Mizrahi, Kevin Fortier, and Peter Maunz.  Demonstration of qubit operations below a rigorous fault tolerance threshold with gate set tomography.  _Nature Communications_ , 8, Feb 2017.  doi: 10.1038/ncomms14485.
  * Bennett and Brassard [1984] Charles H Bennett and Gilles Brassard.  Quantum cryptography: Public-key distribution and coin tossing.  In _Proceedings of IEEE International Conference on Computers, Systems and Signal Processing_ , pages 175–179. Bangalore, India, IEEE Press, 1984\.
  * Bennett [1992] Charles H Bennett.  Quantum cryptography using any two nonorthogonal states.  _Phys. Rev. Lett._ , 68:3121–3124, May 1992.  doi: 10.1103/PhysRevLett.68.3121.  URL https://link.aps.org/doi/10.1103/PhysRevLett.68.3121.
  * Woodhead et al. [2013] Erik Woodhead, Charles Ci Wen Lim, and Stefano Pironio.  Semi-device-independent QKD based on BB84 and a CHSH-type estimation.  In _Theory of Quantum Computation, Communication, and Cryptography_ , pages 107–115, Berlin, Heidelberg, 2013. Springer Berlin Heidelberg.  ISBN 978-3-642-35656-8.
  * de Vicente [2017] Julio I. de Vicente.  Shared randomness and device-independent dimension witnessing.  _Phys. Rev. A_ , 95:012340, Jan 2017.  doi: 10.1103/PhysRevA.95.012340.  URL https://link.aps.org/doi/10.1103/PhysRevA.95.012340.
  * Navascués and Popescu [2014] Miguel Navascués and Sandu Popescu.  How energy conservation limits our measurements.  _Phys. Rev. Lett._ , 112:140502, Apr 2014.  doi: 10.1103/PhysRevLett.112.140502.  URL https://link.aps.org/doi/10.1103/PhysRevLett.112.140502.
  * Renes et al. [2004] Joseph M. Renes, Robin Blume-Kohout, A. J. Scott, and Carlton M. Caves.  Symmetric informationally complete quantum measurements.  _Journal of Mathematical Physics_ , 45(6):2171–2180, 2004.  doi: 10.1063/1.1737053.  URL https://doi.org/10.1063/1.1737053.
  * Brierley et al. [2009] Stephen Brierley, Stefan Weigert, and Ingemar Bengtsson.  All mutually unbiased bases in dimensions two to five, 2009.
  * Nielsen and Chuang [2011] Michael A. Nielsen and Isaac L. Chuang.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_.  Cambridge University Press, USA, 10th edition, 2011.  ISBN 1107002176.
  * Bargmann [1964] V. Bargmann.  Note on Wigner’s theorem on symmetry operations.  _Journal of Mathematical Physics_ , 5(7):862–868, 1964.  doi: 10.1063/1.1704188.  URL https://doi.org/10.1063/1.1704188.
  * Sun [1991] Ji-Guang Sun.  Perturbation bounds for the Cholesky and QR factorizations.  _BIT Numerical Mathematics_ , 31(2):341–352, Jun 1991.  ISSN 1572-9125.  doi: 10.1007/BF01931293.  URL https://doi.org/10.1007/BF01931293.
  * Arias-Castro et al. [2020] Ery Arias-Castro, Adel Javanmard, and Bruno Pelletier.  Perturbation bounds for Procrustes, classical scaling, and trilateration, with applications to manifold learning.  _Journal of Machine Learning Research_ , 21:15–1, 2020\.
  * Bechmann-Pasquinucci and Gisin [2003] H. Bechmann-Pasquinucci and N. Gisin.  Intermediate states in quantum cryptography and Bell inequalities.  _Phys. Rev. A_ , 67:062310, Jun 2003.  doi: 10.1103/PhysRevA.67.062310.  URL https://link.aps.org/doi/10.1103/PhysRevA.67.062310.
  * D’Ariano et al. [2005] Giacomo Mauro D’Ariano, Paoloplacido Lo Presti, and Paolo Perinotti.  Classical randomness in quantum measurements.  _Journal of Physics A: Mathematical and General_ , 38(26):5979, 2005.
  * Oszmaniec and Biswas [2019] Michał Oszmaniec and Tanmoy Biswas.  Operational relevance of resource theories of quantum measurements.  _Quantum_ , 3:133, April 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-04-26-133.  URL https://doi.org/10.22331/q-2019-04-26-133.
  * Takagi and Regula [2019] Ryuji Takagi and Bartosz Regula.  General resource theories in quantum mechanics and beyond: Operational characterization via discrimination tasks.  _Phys. Rev. X_ , 9:031053, Sep 2019.  doi: 10.1103/PhysRevX.9.031053.  URL https://link.aps.org/doi/10.1103/PhysRevX.9.031053.
  * Uola et al. [2019] Roope Uola, Tristan Kraft, Jiangwei Shang, Xiao-Dong Yu, and Otfried Gühne.  Quantifying quantum resources with conic programming.  _Phys. Rev. Lett._ , 122:130404, Apr 2019.  doi: 10.1103/PhysRevLett.122.130404.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.130404.
  * Tavakoli [2020] Armin Tavakoli.  Semi-device-independent certification of independent quantum state and measurement devices.  _arXiv preprint arXiv:2003.03859_ , 2020.
  * Higham [2002] Nicholas J. Higham.  _Accuracy and Stability of Numerical Algorithms_.  Society for Industrial and Applied Mathematics, second edition, 2002.  doi: 10.1137/1.9780898718027.  URL https://epubs.siam.org/doi/abs/10.1137/1.9780898718027.
  * Chang and Stehlé [2010] X.-W. Chang and D. Stehlé.  Rigorous perturbation bounds of some matrix factorizations.  _SIAM Journal on Matrix Analysis and Applications_ , 31(5):2841–2859, 2010.  doi: 10.1137/090778535.  URL https://doi.org/10.1137/090778535.
  * Hurley and Cattell [1962] John R Hurley and Raymond B Cattell.  The Procrustes program: Producing direct rotation to test a hypothesized factor structure.  _Behavioral science_ , 7(2):258–262, 1962.