Device-independent point estimation from finite data and its application to device-independent property estimation

TL;DR

This paper introduces practical, statistically sound methods for estimating quantum distributions in device-independent scenarios from finite data, enabling reliable property estimation like entanglement without prior Bell inequality knowledge.

Contribution

It develops the first physical, unique, and computationally feasible point estimation techniques for quantum distributions in device-independent settings from finite data.

Findings

Estimates are physical, unique, and computationally tractable.

Methods enable trustworthy entanglement estimation from finite data.

No prior Bell inequality knowledge needed for property estimation.

Abstract

The device-independent approach to physics is one where conclusions are drawn directly from the observed correlations between measurement outcomes. In quantum information, this approach allows one to make strong statements about the properties of the underlying systems or devices solely via the observation of Bell-inequality-violating correlations. However, since one can only perform a {\em finite number} of experimental trials, statistical fluctuations necessarily accompany any estimation of these correlations. Consequently, an important gap remains between the many theoretical tools developed for the asymptotic scenario and the experimentally obtained raw data. In particular, a physical and concurrently practical way to estimate the underlying quantum distribution has so far remained elusive. Here, we show that the natural analogs of the maximum-likelihood estimation technique and the least-square-error estimation technique in the device-independent context result in point estimates of the true distribution that are physical, unique, computationally tractable and consistent. They thus serve as sound algorithmic tools allowing one to bridge the aforementioned gap. As an application, we demonstrate how such estimates of the underlying quantum distribution can be used to provide, in certain cases, trustworthy estimates of the amount of entanglement present in the measured system. In stark contrast to existing approaches to device-independent parameter estimations, our estimation does not require the prior knowledge of {\em any} Bell inequality tailored for the specific property and the specific distribution of interest.

Figures

Tables

Table 1. Table 1: A summary of the various regularization methods discussed and some of their key properties. Note that all these considered methods are invariant under relabelings, in the sense that for any relabeling operation M 𝑀 M , we have M ​ P → Reg ​ ( f → ) = P → Reg ​ ( M ​ f → ) 𝑀 subscript → 𝑃 Reg → 𝑓 subscript → 𝑃 Reg 𝑀 → 𝑓 M\vec{P}_{\text{\tiny Reg}}(\vec{f})=\vec{P}_{\text{\tiny Reg}}(M\vec{f}) . In the event that the output of the regularization method is unique, we indicate in the last column if the method provides (in general) unbiased estimation of Bell violation.

Method	Output set $ℛ$	Optimization	Unique?	Unbiased?
Projection	$\overline{𝒩}$	-
NNA1	$𝒩$	LP	✗	-
LSN		SOCP		✗
NNA∞		LP	✗	-
NQA1	$𝒬$	SDP	✗	-
LS	(approximated by ${\widetilde{𝒬}}_{\ell }\supseteq 𝒬$)	SDP		✗
NQA∞		SDP	✗	-
MLN	$𝒩$	CP		✗
ML	$𝒬$	CP		✗
	(approximated by ${\widetilde{𝒬}}_{\ell }\supseteq 𝒬$)

Equations

$$\begin{array}{rl}P(a|x,y)\equiv \sum _{b}P(a,b|x,y)& =P(a|x,y^{\prime }),\forall a,x,y,y^{\prime },\\ P(b|x,y)\equiv \sum _{a}P(a,b|x,y)& =P(b|x^{\prime },y),\forall b,x,x^{\prime },y.\end{array}$$

$$D_{\text{KL}}\left(\vec{f}||\vec{P}\right)=\sum _{a,b,x,y}f(x,y)f(a,b|x,y){\log }_2\left[\frac{f(a,b|x,y)}{P(a,b|x,y)}\right],$$

$$\sum _{a,b}P(a,b|x,y)$$

$$P(a,b|x,y)$$

$$\begin{array}{rl}P(a|x,y)\equiv \sum _{b}P(a,b|x,y)& =P(a|x,y^{\prime }),\forall a,x,y,y^{\prime },\\ P(b|x,y)\equiv \sum _{a}P(a,b|x,y)& =P(b|x^{\prime },y),\forall b,x,x^{\prime },y.\end{array}$$

$$\underset{\vec{P}\in \overline{\mathcal{N}}}{\min }{\Vert \vec{P}-\vec{f}\Vert }_2={\Vert {\vec{P}}_{\Pi }(\vec{f})-\vec{f}\Vert }_2$$

$$\vec{f}={\vec{f}}_{\overline{\mathcal{N}}}+{\vec{f}}_{\text{SI}},$$

$${\vec{P}}_{\Pi }(\vec{f})=\Pi \vec{f}$$

$${\Pi }_{abxy,a^{\prime }{b}^{\prime }{x}^{\prime }{y}^{\prime }}=\text{mathbbm}{1}_{16}-\frac{1}{16}\sum _{(i,j,k,l)\in \mathcal{I}}{i}^{a+a^{\prime }}{j}^{b+b^{\prime }}{k}^{x+x^{\prime }}{l}^{y+y^{\prime }},$$

$$⟨A_x⟩=\sum _{A=\pm 1}A\text{leavevmode}P(A|x)=\sum _{a=0,1}(-1{)}^{a}P(a|x)$$

$$⟨A_x{B}_y⟩=\sum _{A,B}A\text{leavevmode}B\text{leavevmode}P(a,b|x,y)=\sum _{a,b}(-1{)}^{a+b}P(a,b|x,y).$$

$$\begin{split}P(a,b|x,y)=\frac{1}{4}\big{[}1&+(-1)^{a}\left<A_{x}\right>+(-1)^{b}\left<B_{y}\right>\\ &+(-1)^{a+b}\left<A_{x}B_{y}\right>\big{]}\end{split}$$

$${\tilde{P}}_{\text{A}}(a|x)=\frac{1}{2}\sum _{b,y}P(a,b|x,y)$$

$${\tilde{P}}_{\text{A}}(a|x)=\frac{1}{m^{n-1}}\sum _{b,y,c,z...}P(a,b,c,...|x,y,z...),$$

$$⟨A_x^i⟩=\sum _a{c}_{ia}\tilde{P}(a|x),0\le i\le k-2,$$

$$⟨A_x^i{B}_y^j⟩=\sum _{a,b}{c}_{ia}{c}_{jb}\tilde{P}(a,b|x,y),0\le i,j\le k-2.$$

$$\!\!\!\!\vec{f}=\left[\begin{array}[]{c|c}f(a,b|0,0)&f(a,b|0,1)\\ \hline\cr f(a,b|1,0)&f(a,b|1,1)\\ \end{array}\right]=\tfrac{1}{10}\left[\begin{array}[]{cc|cc}3&0&7&0\\ 1&6&1&2\\ \hline\cr 5&1&1&6\\ 1&3&3&0\\ \end{array}\right]$$

$$\vec{P}_{\Pi}(\vec{f})\!=\!\left[\begin{array}[]{c|c}P(a,b|0,0)&P(a,b|0,1)\\ \hline\cr P(a,b|1,0)&P(a,b|1,1)\\ \end{array}\right]\!=\!\tfrac{1}{40}\!\!\left[\begin{array}[]{cc|cc}18&2&20&0\\ 2&18&4&16\\ \hline\cr 19&7&7&19\\ 1&13&17&-3\\ \end{array}\right],$$

$${\vec{P}}_{\mathrm{LS}}(\vec{f})=\underset{\vec{P}\in \mathcal{Q}}{\mathrm{argmin}}||\vec{f}-\vec{P}|{|}_2.$$

$${\vec{P}}_{\mathrm{LS}}(\vec{f})$$

$$=\underset{\vec{P}\in {\tilde{\mathcal{Q}}}_{\ell }}{\mathrm{argmin}}\Vert \vec{P}-{\vec{f}}_{\overline{}\mathcal{N}}{\Vert }_2=\underset{\vec{P}\in {\tilde{\mathcal{Q}}}_{\ell }}{\mathrm{argmin}}\Vert \vec{P}-{\vec{P}}_{\Pi }(\vec{f}){\Vert }_2$$

$$\begin{array}{rl}& \underset{\vec{P}\in {\tilde{\mathcal{Q}}}_{\ell }}{\mathrm{argmin}}s\\ & \text{s.t.}\left(\begin{array}{cc}s\text{mathbbm}1& \vec{f}-\vec{P}\\ {\vec{f}}^{\text{T}}-{\vec{P}}^{\text{T}}& s\end{array}\right)⪰0,\end{array}$$

$$D_{\text{KL}}\left(\vec{v}||\vec{w}\right)=\sum _i{v}_i{\log }_2\left(v_i/w_i\right),$$

$$\begin{array}{rl}& \underset{\vec{P}\in \mathcal{C}}{\min }{D}_{\text{KL}}\left(\vec{f}||\vec{P}\right)\\ =& \underset{\vec{P}\in \mathcal{C}}{\min }\sum _{a,b,x,y}f(x,y)f(a,b|x,y){\log }_2\left[\frac{f(a,b|x,y)}{P(a,b|x,y)}\right],\\ =& \kappa +\underset{\vec{P}\in \mathcal{C}}{\min }-\frac{1}{N}\sum _{a,b,x,y}{N}_{a,b,x,y}{\log }_{2}P(a,b|x,y),\\ =& \kappa -\frac{1}{N}\underset{\vec{P}\in \mathcal{C}}{\max }\sum _{a,b,x,y}{\log }_{2}P(a,b|x,y{)}^{N_{a,b,x,y}},\\ =& \kappa -\frac{1}{N}\underset{\vec{P}\in \mathcal{C}}{\max }{\log }_2\prod _{a,b,x,y}P(a,b|x,y{)}^{N_{a,b,x,y}},\end{array}$$

$$\underset{\vec{P}\in {\tilde{\mathcal{Q}}}_{\ell }}{\mathrm{argmin}}\sum _{a,b,x,y}f(x,y)f(a,b|x,y){\log }_2\left[\frac{f(a,b|x,y)}{P(a,b|x,y)}\right].$$

$$\min$$

$$A\vec{x}=\vec{b},$$

$$\vec{x}\in K,$$

$$K_{\text{exp}}=\{(u,v,w)\text{leavevmode}|\text{leavevmode}v{e}^{u/v}\le w,v\ge 0\}.$$

$${\vec{P}}_{\text{ML}}(\vec{f})=\mathrm{argmin}$$

$$\text{s.t.}\chi$$

$$\mathrm{tr}[F_{abxy}\chi ]$$

$$\mathrm{tr}[G_k\chi ]$$

$${\vec{P}}_{\text{ML}}(\vec{f})=\mathrm{argmax}$$

Full text

I Introduction

The proper analysis of empirical data is an indispensable part in the development of both science and technologies. In quantum information, for instance, the careful preparation followed by the proper characterization of quantum systems (which includes estimating reliably the prepared quantum state or the confirmation that it possesses certain desired properties) is often the first step to many quantum information processing protocols.

In practice, however, the execution of this preliminary task is far from trivial. For example, systematic uncertainties arising from various imperfections in the setup may compromise the reliability of the estimate Rosset:2012 ; Moroder2013b ; vanEnk:2013 . Moreover, unavoidable statistical fluctuations result in situations where the ideal, theoretical description become inapplicable. Hence, quantum state estimation Hradil1997 using real data is a daunting task Blume-Kohout2010 ; arXiv:1202.5270 ; Christandl2012 ; Sugiyama2013 ; Shang2013 ; Faist2016 where there remains an ongoing debate on the preferred approach (see, e.g., arXiv:1405.5350:brief ; Schwemmer2015 and references therein).

Interestingly, the first of these problems can be circumvented, to some extent, by the so-called device-independent approach Scarani2012 ; Brunner:RMP . There, the nature of the devices employed is deduced directly from the measurement statistics Mayers1998 ; Mayers2004 , without relying on any assumption about the devices’ detailed functioning or the associated Hilbert space dimension. Consequently, robust characterizations of quantum systems and instruments are now in principle possible with minimal assumptions. Likewise, the distribution of shared secret keys Ekert1991 ; Acin2007 ; Vazirani2014 and the generation of random bits—secured by the laws of physics ColbeckPhD ; Pironio:2010aa ; Colbeck2011 —are now possibilities at our disposal.

Crucially, in order to make nontrivial statements from the empirical data, the latter have to be Bell-inequality Bell1964 violating, which cannot arise from local measurements on a separable quantum state Werner1989 . Indeed, the extent of the observed Bell-inequality violation can be used to provide an estimate, e.g., of the amount of shared entanglement Moroder2013 ; Toth2015 , or even the incompatibility between the measurements Chen2016 ; Cavalcanti2016 employed. Nonetheless, as with the case where the measurement devices are fully characterized, the underlying distribution—which serves as the analog of a quantum state in this black-box setting—contains all the available information and thus generally provides a much better estimate of the system’s properties Zhang2011 ; Zhang2013 ; Bancal2014 ; Nieto2014 .

Indeed, various theoretical tools taking into account the full quantum distribution have been developed for device-independent characterizations: from the nature of the (multipartite) entanglement Bancal2011 ; LiangPRL2015 ; Baccari2017 present to their quantification Moroder2013 ; Toth2015 , from the steerability Wiseman2007 of the underlying state to the incompatibility of the measurements employed Chen2016 ; Cavalcanti2016 , and from the minimal compatible Hilbert space dimension Brunner2008 ; NavascuesPRX ; NavascuesPRL2015 to the self-testing Mayers1998 ; Mayers2004 of the quantum apparatus Yang2014 ; BancalPRA2015 . Stemming from the algorithmic characterization of the set of quantum distributions due to Navascués-Pironio-Acín (NPA) Navascues2007 ; Navascues2008a , they share the common assumption that the estimated distributions satisfy the physically motivated conditions of nonsignaling Popescu1994 ; Barrett2005 .

However, raw distributions estimated from the relative frequencies of experimental outcomes—due to statistical fluctuations—generically do not satisfy these conditions. As such, none of the aforementioned tools can be directly applied to experimentally observed statistics. In other words, while the device-independent approach offers an elegant solution to overcome the problem of mistrusting the measurement devices, there remains an important gap between the theoretical tools developed for such purposes and the actual data available from any Bell experiment.

For the very specific problem of device-independent randomness certification and quantum key distribution, techniques based on hypothesis testing have been shown, respectively, in Pironio:2010aa ; Pironio2013 ; arXiv:1611.00352:brief ; arXiv:1702.05178:brief ; Knill2017 and in Dupuis:1607.01796:brief ; Rotem2017 to be applicable even in the presence of finite statistics. These approaches are, however, very problem specific, and it is not yet known how to generalize them for the general problem of device-independent characterizations. Here, we consider an alternative approach inspired by estimation theory, which consists in constructing a point estimate for the underlying quantum distribution from the observed frequencies. In particular, we show that the natural analogs of two physical estimators employed in usual quantum state tomography, namely, maximum-likelihood (ML) estimation and least-square-error estimation, also serve as sound estimators in the device-independent context, thereby allowing us to regularize these raw data and obtain a direct estimation of the properties of interest through the respective theoretical techniques.

Although there have been attempts to perform regularizations for device-independent property estimations Bancal2014 ; Schwarz2016 and for the quantification of nonlocality Bernhard2014 , these proposals turn out to suffer from the drawback of generating estimates that can be either nonphysical or nonunique. In contrast, our methods are provably free from such problems. Armed with these point estimates of the underlying distribution, a device-independent estimation of the property of interest then follows naturally by applying the algorithmic tools mentioned above (see Fig. 1).

II Preliminaries

The starting point of device-independent estimations is a Bell experiment. Consider the simplest Bell scenario where Alice and Bob each randomly performs two possible measurements (labeled, respectively, by $x,y\in\{0,1\}$), and where each measurement gives binary outcomes (labeled, respectively, by $a,b\in\{0,1\}$). Generalization of our discussion to other finite Bell scenarios is obvious from the context. The correlations between their measurement outcomes can be summarized by a vector of joint conditional probability distributions $\vec{P}=\{P(a,b|x,y)\}_{a,b,x,y}\in\mathbb{R}^{16}$.

Denote by $\rho$ the state shared by Alice and Bob, and by $M^{\text{A}}_{a|x}$ ($M^{\text{B}}_{b|y}$) the positive-operator-valued-measure elements associated with their measurements. Born’s rule dictates that for all $a,b,x,y$, the conditional probability distributions read as $P(a,b|x,y)\stackrel{{\scriptstyle\mathcal{Q}}}{{=}}\operatorname{tr}\left(\rho\,M^{\text{A}}_{a|x}\otimes M^{\text{B}}_{b|y}\right)$ where the positivity and the normalization of probabilities demand that $M^{\text{A}}_{a|x},M^{\text{B}}_{b|y}\succeq 0$ (matrix positivity) and $\sum_{a}M^{\text{A}}_{a|x}=\mathbb{I}_{A},\sum_{b}M^{\text{B}}_{b|y}=\mathbb{I}_{B}$ with $\mathbb{I}_{A},\mathbb{I}_{B}$ being identity operators. Throughout, we use $\mathcal{Q}$ to denote the set of quantum distributions, i.e., the collection of $\vec{P}$ that follows from Born’s rule.

Importantly, quantum distributions satisfy the nonsignaling conditions Popescu1994 ; Barrett2005 , i.e., their marginal distributions are independent of the measurement choice of the distant party:

[TABLE]

In an experiment, the underlying quantum distribution $P(a,b|x,y)$ is often estimated by computing the relative frequency, i.e., $P(a,b|x,y)\approx f(a,b|x,y)=\tfrac{N_{a,b,x,y}}{N_{x,y}}$ where $N_{a,b,x,y}$ is the number of coincidences registered for the combination of outcomes and settings $(a,b,x,y)$ while $N_{x,y}=\sum_{ab}N_{a,b,x,y}$ is the total number of trials pertaining to the measurement choice $(x,y)$.

Of course, in the asymptotic limit of a large number of trials, i.e., when $N_{\text{\tiny trials}}=\min_{x,y}N_{x,y}\rightarrow\infty$, the difference between the true distribution $\vec{P}$ and the relative frequency $\vec{f}=\{f(a,b|x,y)\}_{a,b,x,y}$ vanishes. In practice, as $\max_{x,y}N_{x,y}$ is necessarily finite, not only is this difference nonzero but $\vec{f}$ typically also violates the weaker requirement of the nonsignaling conditions [see Eq. (1)]. As mentioned above, this discrepancy between theory and practice immediately renders many of the tools developed for device-independent characterizations inapplicable.

III Regularization methods (estimators)

To overcome this mismatch, one may project the observed frequency $\vec{f}$ onto an affine subspace $\overline{\mathcal{N}}$ of $\mathbb{R}^{16}$ which contains only $\vec{P}$’s that satisfy Eq. (1). For example, if one demands that this projection (via the corresponding projector $\Pi$) commutes Renou2017 with all possible permutations of the labels for parties, settings, and outcomes (e.g., $a=0\leftrightarrow a=1$), then it happens to be equivalent to finding the unique minimizer of the least-square-error problem: \vec{P}_{\mbox{\tiny\Pi}}(\vec{f})=\operatorname*{argmin}_{\vec{P}\in\overline{\mathcal{N}}}||\vec{f}-\vec{P}||_{2}, where $\|\mathbf{.}\|_{p}$ denotes the $p$-norm; the regularization invoked in Bancal2014 is precisely an application of such a projection (see Appendix B).

Albeit intuitive and straightforward, such a projection suffers from the serious drawback that it may give “negative probabilities” (see Appendix B.4 for an explicit example). Indeed, the possibility of giving rise to an unphysical estimate is a problem that such a projection shares with the linear inversion technique employed in standard quantum state tomography (see, e.g., Schwemmer2015 ). Moreover, even when \vec{P}_{\mbox{\tiny\Pi}}(\vec{f}) represents a legitimate probability vector, it may well be outside the quantum set $\mathcal{Q}$. To overcome these issues, one is naturally led to the least-square (LS) estimator in the device-independent context, i.e., $\vec{P}_{{\rm LS}}({\vec{f}})=\operatorname*{argmin}_{\vec{P}\in\mathcal{Q}}||\vec{f}-\vec{P}||_{2}$.

We thus see that various estimators (see also Schwarz2016 ; Bernhard2014 and Appendix H) map $\vec{f}$ to a regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ that is non-negative and normalized and which satisfies the nonsignaling conditions. However, for a regularization method to be relevant for subsequent property estimation, it is also convenient that $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ is uniquely determined by $\vec{f}$. In particular, for a given $\vec{f}$, a nonunique estimator may give rise to $\vec{P}_{\text{\tiny Reg}}(\vec{f})$’s with drastically different properties, e.g., some being Bell-inequality violating (and therefore implying some nontrivial features of the underlying system) and some not [which renders that particular $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ useless for device-independent property estimation]. An ambiguity then arises: which of these estimates should we rely on for subsequent property estimation? A possibility would be to consider the worst case over all such estimates, but this clearly complicates the property estimation as one would now need to consider an entire solution set $\{\vec{P}_{\text{\tiny Reg}}(\vec{f})\}$ (the characterization of which is generally nontrivial). This makes evident the inconveniences of 1-norm estimators, and more generally nonunique estimators in the present context. The regularization procedure previously considered in Schwarz2016 ; Bernhard2014 —both being 1-norm estimators—precisely suffers from this nonuniqueness drawback.

In this regard, note that a regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ obtained from minimizing a strictly convex function $g$ over a convex set (such as $\mathcal{Q}$) is provably unique, and is determined by $\vec{f}$ and $g$ (see, e.g., Theorem 8.3 of Beck2018 ). Using this observation, we show in Appendix F.1 that the aforementioned LS estimator is unique. Similarly, the device-independent analog of the ML estimator Banaszek1999 is provably unique (see Appendix F.1 for a proof). To this end, consider the Kullback-Leibler (KL) divergence vanDam2005 ; Acin2005 ; Kullback1951 (i.e., the relative entropy Cover:Book ) from some $\vec{P}\in\mathcal{Q}$ to $\vec{f}$:

[TABLE]

where $f(x,y)$ is the relative frequency of choosing the measurement settings labeled by $(x,y)$.

The quantity $D_{\text{\tiny KL}}\left(\vec{f}||\vec{P}\right)$ can be seen as a measure of “statistical closeness” vanDam2005 ; Acin2005 between $\vec{f}$ and $\vec{P}$. Indeed, its minimization over $\vec{P}\in\mathcal{Q}$ is equivalent Cover:Book ; vanDam2005 to maximizing the likelihood of producing the observed frequency by $\vec{P}\in\mathcal{Q}$ (we provide in Appendix D.1 a proof adapted to the present context). The unique minimizer of $D_{\text{\tiny KL}}\left(\vec{f}||\vec{P}\right)$ over $\vec{P}\in\mathcal{Q}$, i.e., $\vec{P}_{\text{\tiny ML}}({\vec{f}})=\operatorname*{argmin}_{\vec{P}\in{\mathcal{Q}}}D_{\text{\tiny KL}}\left(\vec{f}||\vec{P}\right)$, therefore serves as the equivalent of the ML estimator in the device-independent context. Hereafter, we focus predominantly on this operationally well-motivated estimator. For further details of the LS estimator and some other plausible regularization methods, see, respectively, Appendix C and Appendix H.

As it stands, since there is no known exact characterization of $\mathcal{Q}$ using only finite resources, the ML estimator cannot be computed exactly. Nonetheless, via a converging hierarchy of semidefinite programs (SDPs) Navascues2007 ; Navascues2008a ; Doherty2008 ; Moroder2013 , one can in principle obtain an arbitrary good approximation to $\mathcal{Q}$. To fix ideas, we hereafter focus on employing the hierarchy $\tilde{\mathcal{Q}}_{\ell}$ of approximations to $\mathcal{Q}$ discussed in Moroder2013 and Vallins2017 . The lowest level of this hierarchy $\tilde{\mathcal{Q}}_{1}$ gives a decent outer approximation of $\mathcal{Q}$ known as the almost-quantum set Navascues2015 . In general, $\tilde{\mathcal{Q}}_{\ell}\subseteq\tilde{\mathcal{Q}}_{\ell-1}$ for all $\ell\geq 2$ and $\lim_{\ell\to\infty}\tilde{\mathcal{Q}}_{\ell}\to\mathcal{Q}$. For any fixed $\ell$, although the nonlinear optimization problem $\operatorname*{argmin}_{\vec{P}\in{\mathcal{Q}}}D_{\text{\tiny KL}}\left(\vec{f}||\vec{P}\right)$ does not appear to be a semidefinite program, we show in Appendix D.2 that it belongs to a more general class of convex optimization problems Boyd2004Book — an exponential conic program. A minimization of the KL divergence with NPA constraints is thus also efficiently solvable on a computer with a numerical precision of $10^{-6}$ or better.

IV Notable properties of point estimates

The uniqueness of $\vec{P}_{\text{\tiny ML}}({\vec{f}})$ and the nonnegativity of the KL divergence ensure that our estimators are consistent Shao2003Book , in the sense that they provide an estimate that converges to the true distribution $\vec{P}_{\mathcal{Q}}$ in the asymptotic limit of $N_{\text{\tiny trials}}\to\infty$. This can be seen by noting that in the asymptotic limit, $\vec{f}\to\vec{P}_{\mathcal{Q}}$, the nonnegativity of the KL divergence then implies that the unique minimizer of $D_{\text{\tiny KL}}\left(\vec{f}||\vec{P}\right)$ over $\vec{P}\in\mathcal{Q}$ is necessarily given by $\vec{P}=\vec{P}_{\mathcal{Q}}$. [Likewise, the LS estimator $\vec{P}_{{\rm LS}}({\vec{f}})$ is provably consistent.]

In practice, one would evidently be more interested in how these methods fare for finite values of $N_{\text{\tiny trials}}$. To gain insights into this, we carry out extensive numerical simulations by (i) picking some ideal $\vec{P}_{\mathcal{Q}}$, (ii) numerically simulating the outcomes of a Bell experiment according to $\vec{P}_{\mathcal{Q}}$ and computing the relative frequency $\vec{f}$, (iii) computing the point estimate $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ and calculating various quantities of interest, and (iv) repeating steps i-iii $10^{4}$ times for $N_{\text{\tiny trials}}=10^{2},10^{3},\ldots,10^{10}$ for all $x,y$. For simplicity, we take $N_{x,y}=N_{\text{\tiny trials}}$, i.e., a constant independent of $x,y$ [this amounts to setting $f(x,y)$ as a constant in Eq. (2)].

Our numerical results in Appendix E.2 suggest that in general, the difference between $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ and the ideal distribution $\vec{P}_{\mathcal{Q}}$, quantified, e.g., via $||\vec{P}_{\text{\tiny Reg}}(\vec{f})-\vec{P}_{\mathcal{Q}}||_{1}$ (or other $p$ norms), diminishes, as with $||\vec{f}-\vec{P}_{\mathcal{Q}}||_{1}$, at a rate proportional to $\tfrac{1}{\sqrt{N_{\text{\tiny trials}}}}$. Similar convergence is also observed for $D_{\text{\tiny KL}}\left(\vec{P}_{\text{\tiny ML}}({\vec{f}})||\vec{P}_{\mathcal{Q}}\right)$. Moreover, although we have only employed an outer approximation to $\mathcal{Q}$ in the regularization step iii, as our example below illustrates, the regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ can already be used to perform reasonable device-independent property estimations.

V Application to device-independent estimations

As a concrete example of such property estimations via regularization (see Fig. 1), consider $\vec{P}_{\mathcal{Q}}=\vec{P}_{\mathcal{Q}}^{\tiny\tau_{1.25}}$, a quantum distribution considered in a recent Bell test Christensen2015 (see Appendix E.1 for details). Device-independent estimations of the underlying negativity Vidal2002 (a well-known entanglement measure) based on ideal quantum distributions are known to be possible Moroder2013 . Here, we illustrate how such an estimation can be realized for finite data through the regularization of $\vec{f}$. To facilitate comparison, we plot in Fig. 2 the average negativity N$(\rho)$ of the underlying state $\rho$ estimated from the regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ (via the SDP described in Moroder2013 ) against that deducible from the amount of Clauser-Horne-Shimony-Holt CHSH (CHSH) Bell-inequality violation $\mathcal{S}_{\text{\tiny CHSH}}$ Moroder2013 : N$(\rho)\geq\tfrac{\mathcal{S}_{\text{\tiny CHSH}}-2}{4\sqrt{2}-4}$ for $\mathcal{S}_{\text{\tiny CHSH}}\in[2,2\sqrt{2}]$.

A few features of this comparison are worth noting. First, since the negativity estimated directly from the CHSH Bell-inequality violation $\mathcal{S}_{\text{\tiny CHSH}}$ of $\vec{f}$ depends linearly on this violation, the mean value of the negativity estimated hardly depends on $N_{\text{\tiny trials}}$, and is suboptimal. In fact, even in the asymptotic limit, a negativity estimation based on $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ and $\mathcal{S}_{\text{\tiny CHSH}}$ is suboptimal. In contrast, the mean value of the negativity estimated from the regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ rapidly converges to the true value as $N_{\text{\tiny trials}}$ increases; already at $N_{\text{\tiny trials}}=10^{4}$, this mean value only differs from the true value by a few percents. Second, note that our negativity estimations based on $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ clearly systematically underestimate the amount of negativity present, which is in strong contrast with the results presented in Schwemmer2015 for the non-device-independent scenario using both least-square and maximum-likelihood estimators. (For further examples of underestimation using $\vec{f}$ sampled from other quantum distributions, see Appendix F.2.)

Of course, instead of the CHSH Bell inequality, one could hope to improve the negativity estimation by considering a Bell-inequality the quantum violation of which is optimized for the negativity estimation of $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ [see Eq. (29) and Appendix G for details]. Such an optimized device-independent witness may be obtained, e.g., by feeding the DI algorithm with some regularized distribution sampled from $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$, as indicated in Fig. 1. In practice, however, the relative frequencies sampled from $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ turn out to give—independent of $N_{\text{\tiny trials}}$—about half the time, a Bell violation more than that allowed by quantum theory, thereby rendering negativity estimation from this Bell violation impossible in all these cases.

VI Discussion

The device-independent state estimation problem is the core of all state estimation problems from finite data, as it addresses the generic problem of matching empirical data (subjected to statistical fluctuations) with the ideal, theoretical description given by Born’s rule. The recent demonstration of loophole-free Bell tests hensen_loophole-free_2015 ; shalm_strong_2015 ; giustina_significant-loophole-free_2015 ; rosenfeld_event-ready_2017 has made it clear that the development of reliable techniques for the device-independent estimation of underlying properties from finite data not only is of fundamental interest but also would play an indispensable role in the next generation of quantum information protocols. In fact, although our focus is on a fully device-independent setting, the insights obtained thereof are also relevant in analogous problems in a partially device-independent scenario, such as those incurred in a quantum steering experiment Cavalcanti:2015aa .

In this paper, we have provided the device-independent analog of the maximum-likelihood and the least-square estimators and shown that they are physical, computationally tractable, and unique. These features render them ideal for bridging the device-independent tools developed for ideal quantum distributions and the experimentally obtained raw data. Generalizing the arguments in Schwemmer2015 , however, it can be shown that all quantum estimators are necessarily biased, as summarized in the following proposition (see Appendix F.2 for a proof and the corresponding definition of being strictly nonorthogonal).

Proposition 1**.**

Let $\mathcal{C}$ be a closed convex subset of the nonsignaling polytope $\mathcal{N}$ (such as $\mathcal{Q}$) with two strictly nonorthogonal extreme points. Any point estimator constrained to give $\vec{P}_{\text{\tiny Reg}}(\vec{f})\in\mathcal{C}$ is necessarily biased.

As stressed above, one of the goals of performing a device-independent state estimation is to obtain therefrom various device-independent parameter (property) estimates (see Fig. 1). Given Proposition 1, one may expect that any such parameter estimates are also biased. Indeed, for data sampled from a distribution that is not Bell-inequality violating, but which is sufficiently near to the boundary of the local polytope, an overestimation of the corresponding negativity is to be expected.111Evidently, due to statistical fluctuations, some of the $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ would lie inside the local polytope, giving a zero negativity value while some other $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ would violate a Bell inequality, giving a strictly positive negativity value. Their average is thus positive, thereby resulting in an overestimation. On the other hand, our numerical results show that for data sampled from some extremal quantum distribution, we observe, instead, an underestimation of the underlying negativity (see Fig. 2) and/or Bell violation (see Fig. 4). Admittedly, for such point estimation to be useful, a more thorough investigation is needed in order to determine when such device-independent estimations are trustworthy (in the sense of not leading to an overestimation).

As one can perform analogous property estimation directly from the observed Bell-inequality violation, our approach of estimation by regularization may seem redundant at first glance. However, as we illustrate in the negativity example, the quality of an estimate obtained from Bell-inequality violation depends heavily on the choice of the inequality and it is a priori not always obvious which Bell inequality (the violation of which is used as a device-independent witness) will provide an optimal estimate. In contrast, our approach yields an optimized Bell-like inequality for witnessing the desired property as a byproduct. Moreover, due to the signaling nature of the relative frequencies, even physically equivalent Collins2004 Bell inequalities may give rise to different estimates Renou2017 , thereby resulting in further ambiguities.

We now briefly comment on some other possibilities for future research. Implicit in our discussion is the assumption that the experimental trials are independent and identically distributed. While this is an often adopted assumption in the non-device-independent setup, its justification in a device-independent setting is far from trivial. A natural line of research thus consists of relaxing this assumption while maintaining the possibility to perform a reliable state estimation.

In close connection to this is the problem of establishing a confidence region: how does one generalize the tools presented here to construct a region of estimates in accordance to, say, some pre-defined likelihood arXiv:1202.5270 ; Christandl2012 ; Sugiyama2013 ; Shang2013 ; Faist2016 ? For some device-independent tasks, specific techniques Zhang2011 ; Zhang2013 ; Pironio2013 ; arXiv:1611.00352:brief ; arXiv:1702.05178:brief ; Dupuis:1607.01796:brief ; Rotem2017 for dealing with finite statistics (possibly with the inclusion of confidence regions) have been developed, but general techniques for establishing confidence regions associated with generic quantum properties are still lacking. To appreciate the importance of constructing these confidence regions, recall from our example that the Bell violation given by the observed relative frequency—due to statistical fluctuations—may give rise to a value beyond that allowed by quantum theory, thus rendering this Bell violation useless for the estimate of a quantum parameter. To this end, we remark that the work of Wills:Unpublished suggests that the hypothesis-testing technique of Zhang2011 , together with the numerical technique developed here, can indeed be used to provide a confidence region for general device-independent estimates. Addressing these questions, however, clearly goes beyond the scope of the present paper and is something that we plan to take up in the sequel to this paper.

Acknowledgements.

We are grateful to Bänz Bessire, Daniel Cavalcanti, Flavien Hirsch, Sacha Schwarz and Paul Skrzypczyk for useful discussions, and to Boris Bourdoncle, Jedrzej Kaniewski, Lukas Knips, as well as a few anonymous referees for useful comments on an earlier version of this paper. This work is supported by the Ministry of Science and Technology, Taiwan (Grant No. 104-2112-M-006-021-MY3), the National Center for Theoretical Sciences of Taiwan, the Swiss National Science Foundations [through the NCCR QSIT as well as Grants No. PP00P2-150579 and No. P2GEP2$\_$162060 (Early Postdoc.Mobility)], the Ontario Research Fund, and the Natural Sciences and Engineering Research Council of Canada.

Appendix A Notations and definitions

Throughout, we label the measurement settings (inputs) by $x,y,z,\ldots$ and the corresponding measurement outcomes by $a,b,c,\ldots$ where each of these labels are elements from some finite sets. The correlations between their measurement outcomes are succinctly summarized by the vector of joint conditional distributions $\vec{P}:=\{P(a,b,c,\ldots|x,y,z,\ldots)\}$. For simplicity, we will focus our discussion on the bipartite cases.

We denote by $\mathcal{P}$ the set of all legitimate probability distributions, i.e., those that obey:

[TABLE]

Moreover, we denote by $\mathcal{N}$ the subset of $\mathcal{P}$ which satisfies also the nonsignaling conditions Popescu1994 ; Barrett2005 :

[TABLE]

It is worth noting that both $\mathcal{P}$ and $\mathcal{N}$ Barrett2005 are convex polytopes, i.e., convex sets having a finite number of extreme points, and are conventionally referred to, respectively, as the signaling and the nonsignaling polytope. Here, we describe both $\mathcal{P}$ and $\mathcal{N}$ in their H-representation using linear equations and inequalities.

Definition 1**.**

The nonsignaling affine space $\overline{\mathcal{N}}\supset\mathcal{N}$ is the smallest-dimensional affine space containing the set $\mathcal{N}$ and is given by the distributions satisfying Eqs. (3a) and (1).

Appendix B Further details about the projection method

Here, we provide some further details about the projection mentioned in the main text.

B.1 Equivalent definitions

The projection method can be defined in three equivalent ways.

It is the minimizer $\vec{P}_{\Pi}(\vec{f})$ of the following optimization problem:

[TABLE] 2. 2.

It is the nonsignaling part $\vec{f}_{\overline{\mathcal{N}}}\in\overline{\mathcal{N}}$, i.e., the first component of the decomposition

[TABLE]

where $\vec{f}_{\text{SI}}$ is the signaling component of $\vec{f}$ and is orthogonal to all vectors in the affine subspace $\overline{\mathcal{N}}$. 3. 3.

It is the result

[TABLE]

of the projection onto $\overline{\mathcal{N}}$ by the linear operator $\Pi$, where $\Pi$ is uniquely determined by the set of commutation relations $\Pi M=M\Pi$ with $M$ being any permutation matrix corresponding to relabeling Collins2004 of outputs, inputs and parties.

For the bipartite Bell scenario with binary inputs and outputs, their equivalence can be shown using the decomposition given in (Renou2017, , Prop. in Sec. 3). For more general Bell scenarios, a proof of their equivalence analogously follows from group representation theory, and will be made available in Rosset:Unpublished .

B.2 Explicit form of the projection matrix in the simplest Bell scenario

Using the notation of Renou2017 , the projection operator $\Pi$ in the bipartite Bell scenario with binary inputs and outputs admits the explicit form

[TABLE]

where $\mathbbm{1}_{16}$ is the 16$\times$16 identity matrix, the sum is carried out over the four quadruplets $\mathcal{I}=\{(+1,-1,-1,\pm 1),(-1,+1,\pm 1,-1)\}$, the rows of $\Pi$ are indexed by $(a,b,x,y)$, while its columns are indexed by $(a^{\prime},b^{\prime},x^{\prime},y^{\prime})$.

B.3 An algorithm for performing the projection for the more general Bell scenarios

In general, in (bipartite) Bell scenarios where parties have binary outputs, it is customary to write $A=(-1)^{a}$ and $B=(-1)^{b}$, and compute the expectation values (correlators Bancal2012 ):

[TABLE]

and similarly for $\left<B_{y}\right>$, while

[TABLE]

Together, the $\left<A_{x}\right>$, $\left<B_{y}\right>$ and $\left<A_{x}B_{y}\right>$ represent a parametrization of the nonsignaling subspace. In particular, the above relations can be inverted to give:

[TABLE]

However, when the distribution $P(a,b|x,y)$ is signaling, the marginals $P(a|x,y)$, $P(b|x,y)$ depend on both inputs $(x,y)$. As shown in Renou2017 , in the case of binary inputs, the (nonsignaling) correlators $\left<A_{x}\right>$ should be computed as $\left<A_{x}\right>=\sum_{a}(-1)^{a}\tilde{P}_{\text{A}}(a|x)$ where $\tilde{P}_{\text{A}}(a|x)$ is given by:

[TABLE]

i.e., averaged uniformly over $y=0,1$. This is the only choice that keeps the projection invariant under permutations of $y$. We make the same choice to compute $\left<B_{y}\right>$. Now, the correlators $\left<A_{x}\right>$, $\left<B_{y}\right>$ and $\left<A_{x}B_{y}\right>$ correspond to a unique distribution $P_{\Pi}(a,b|x,y)$ in the nonsignaling subspace, cf. Eq. (11), which is the result of the projection. This is the essence of the regularization method employed in Bancal2014 .

The construction above extends to an arbitrary number of parties, inputs, and outputs, specified by the tuple $(n,m,k)$. We give a sketch below of this generalization, which satisfies criteria (5)–(7) (a proper proof will be discussed in a future work Rosset:Unpublished ).

The generalization to additional parties and inputs is simple. Write, for example:

[TABLE]

and the same for other marginal distributions $\tilde{P}_{\text{...}}$, by averaging uniformly over all inputs not fixed by the indices of the marginals $\tilde{P}_{\text{...}}$. Then, use these $\tilde{P}_{\text{...}}$ to compute the correlators according to the straightforward multipartite generalizations of Eqs. (9) and (10).

To cater for scenarios with nonbinary outputs, the correlators have to be generalized. We use the framework proposed in Bancal2010 , adapting slightly the notation to the present paper:

[TABLE]

where $c_{ia}=k\leavevmode\nobreak\ \delta_{ia}-1$ and we omitted the coefficient for $i=k-1$ as it is linearly dependent on the others. Note that $\left<A^{0}_{x}\right>=\left<A_{x}\right>$ in the case of binary outcomes ($k=2$). Multipartite correlators are written, for example:

[TABLE]

Starting from a signaling distribution $P(a,b,\ldots|x,y,\ldots)$, we compute the generalized correlators using the averaged marginals of Eq. (13). Then we interpret the resulting correlators as a description of the nonsignaling distribution $P_{\Pi}(a,b,\ldots|x,y,\ldots)$. The whole process is a projection: all operations are linear, and the averaging [e.g., $P(a|x,y)\rightarrow\tilde{P}(a|x)$] is injective for nonsignaling distributions. Importantly, the projection defined by this algorithm commutes with relabelings and thus corresponds to the three equivalent definitions given in Appendix B.1.

B.4 An explicit example showing that the output of the projection method may be nonphysical

Although easy to compute, the projection method may give rise to coefficients of $\vec{P}_{\Pi}(\vec{f})$ that are negative. Indeed, the output space of the projection method is the nonsignaling affine space $\overline{\mathcal{N}}$. We now give an explicit example to illustrate this fact. Consider some relative frequency $\vec{f}$ given in the compact matrix representation:

[TABLE]

where the entries in each block are arranged such that the value of $a$ ($b$) increases downward (rightward). By applying the projection matrix $\Pi$ given in Eq. (8) to this signaling distribution, one obtains:

[TABLE]

which is easily seen to satisfy the nonsignaling conditions of Eq. (1). However, this $\vec{P}_{\Pi}(\vec{f})$ is evidently nonphysical as its entry for $x=y=a=b=1$ is negative.

Appendix C Further details about the device-independent least-square method

Here, we present the details of the device-independent analog of the least-square tomography method. Formally, the method amounts to finding the unique minimizer of the least-square problem:

[TABLE]

C.1 Equivalence to performing a projection and minimization of the 2-norm distance from $\vec{P}_{\Pi}(\vec{f})$ to $\mathcal{Q}$

We now prove that the above optimization is equivalent to first performing the projection method, followed by performing a minimization of the 2-norm distance from $\vec{P}_{\Pi}(\vec{f})$ to $\mathcal{Q}$. For convenience, we shall prove this equivalence, instead, for any converging superset relaxation Navascues2008a ; Doherty2008 ; Moroder2013 $\tilde{\mathcal{Q}}_{\ell}$ of the quantum set $\mathcal{Q}$. The desired equivalence then follows from the fact that $\lim_{\ell\to\infty}\tilde{\mathcal{Q}}_{\ell}=\mathcal{Q}$.

Lemma 1**.**

Given the relative frequency $\vec{f}$, the least-square estimator (when $\tilde{\mathcal{Q}}_{\ell}$ is used to approximate $\mathcal{Q}$) satisfies $\vec{P}_{{\rm LS}}({\vec{f}})=\operatorname*{argmin}_{\vec{P}\in\tilde{\mathcal{Q}}_{\ell}}\|\vec{P}-\vec{P}_{\Pi}(\vec{f})\|_{2}$.

Proof.

First, note from Appendix B.1 that any given relative frequency $\vec{f}$ can be decomposed as $\vec{f}=\vec{f}_{\overline{}\mathcal{N}}+\vec{f}_{\rm SI}$, where $\vec{f}_{\overline{}\mathcal{N}}\in\overline{\mathcal{N}}$ is orthogonal to $\vec{f}_{\rm SI}$. Similarly, for any $\vec{P}\in\tilde{\mathcal{Q}}_{\ell}\subset\mathcal{N}\subset\overline{}\mathcal{N}$, we must have $(\vec{P}-\vec{f}_{\overline{}\mathcal{N}})\in\overline{}\mathcal{N}$, which is orthogonal to $\vec{f}_{\rm SI}$. It then follows from the definition of the least-square method that

[TABLE]

where the second to last equality follows from the orthogonality of $\vec{P}-\vec{f}_{\overline{}\mathcal{N}}$ and $\vec{f}_{\rm SI}$ and the fact that $\|\vec{f}_{\rm SI}\|^{2}$ is a constant in the 2-norm minimization, while the last equality follows from the equivalence between the second and the third definition of the projection method. ∎

Note that although $\vec{P}_{\Pi}(\vec{f})$ is not necessarily in $\tilde{\mathcal{Q}}_{\ell}$, if it happens that $\vec{P}_{\Pi}(\vec{f})\in\tilde{\mathcal{Q}}_{\ell}$ then the equivalent regularization shown in Appendix C implies that $\vec{P}_{{\rm LS}}({\vec{f}})=\vec{P}_{\Pi}(\vec{f})$.

C.2 Formulation as a semidefinite program

When the quantum set $\mathcal{Q}$ is approximated by a superset relaxation $\tilde{\mathcal{Q}}_{\ell}$ that admits a semidefinite programming Boyd2004Book characterization, the optimization problem of Eq. (17) can also be solved as an SDP. To this end, note that the SDP characterization of $\tilde{\mathcal{Q}}_{\ell}$ is achieved in terms of some moment matrix $\chi$ that contains all the entries of $\vec{P}$ as some of its matrix elements.

Using the characterization of positive semidefinite matrices via their Schur complements (see, e.g., Theorem 7.7 of Horn1985Book ), we can then reformulate Eq. (17) (with $\tilde{\mathcal{Q}}_{\ell}$ approximating $\mathcal{Q}$) as:

[TABLE]

where $\mathbbm{1}$ is the identity matrix having the same dimension as the column vector $\vec{f}$, and $\vec{f}^{\text{\tiny T}}$ is the transpose of $\vec{f}$. Evidently, we see from Eq. (18) that the equivalent optimization problem of Eq. (17) now involves only an objection function and matrix inequality constraints that are linear in all their optimization variables: $s$, $\vec{P}$ and some other entries of $\chi$ (that cannot be estimated from experimental data). Thus, as claimed, the minimization of the 2-norm of $\vec{f}-\vec{P}$ over $\vec{P}\in\tilde{\mathcal{Q}}_{\ell}$ can indeed be cast as an SDP.

Appendix D Further details about the Kullback-Leibler divergence and the corresponding regularization method

The Kullback-Leibler (KL) divergence:

[TABLE]

is conventionally defined for unconditional probability distributions $v_{i}$/ $w_{i}$ [such as $P(a,b,x,y)$ and $f(a,b,x,y)$]. In the explicit examples studied, we fixed $P(x,y)=f(x,y)=1/|\mathcal{X}\times\mathcal{Y}|$. This allows us to keep our definition of the ML method valid when $f(x,y)$, sampled from $P(x,y)=$ constant, has itself statistical fluctuations, as reflected in our definition of the KL divergence given in Eq. (2) in the main text.

Note that while the KL divergence is a statistical distance, it is not a metric as it is asymmetric and violates the triangle inequality. To appreciate the relevance of this asymmetry, see vanDam2005 .

D.1 Connection to maximum likelihood

The equivalence between the minimization of the KL divergence over some $\vec{P}\in\mathcal{C}$ (for some set $\mathcal{C}$) to $\vec{f}$ and the maximization of the likelihood of generating $\vec{f}$ from $\vec{P}$ can be seen as follows:

[TABLE]

where $\kappa:=\sum_{a,b,x,y}f(x,y)f(a,b|x,y)\log_{2}f(a,b|x,y)$ is a constant of the optimization, $N:=\sum_{x,y,a,b}N_{a,b,x,y}$, and in the second equality we have used the definition of the relative frequency $\vec{f}$ and the fact that $f(x,y)=\tfrac{N_{x,y}}{N}$. In the last line of Eq. (20), the argument of the maximization is the log likelihood of observing $N_{a,b,x,y}$ times the event labeled by $(x,y,a,b)$ with probability $P(a,b|x,y)$. Hence, we see that the minimization of the KL divergence is equivalent to maximizing the likelihood of generating $\vec{f}$ given $P(a,b|x,y)$.

D.2 Formulation as a conic program

Here, we briefly explain how the device-independent ML method can be formulated and solved as a conic program (CP) with an exponential cone. Recall from the main text that for any given relative frequency $\vec{f}$, the ML estimator (with the quantum set approximated by $\tilde{\mathcal{Q}}_{\ell}$) works by solving:

[TABLE]

A conic program takes the canonical form of:

[TABLE]

where $K$ is a convex cone, such as an exponential cone:

[TABLE]

After discarding the constant term and folding the factors in $f(a,b,x,y)=f(x,y)f(a,b|x,y)$, the minimizer of Eq. (21) can be obtained by solving:

[TABLE]

where $\chi$ is the moment matrix associated with $\tilde{\mathcal{Q}}_{\ell}$, while $F_{abxy}$ and $G_{k}$ encode the equality constraints associated with the structure of this moment matrix. This problem has the conic form of Eq. (D.2):

[TABLE]

where the constraint (25b) is that for a positive semidefinite cone and the constraint (25c) is that for copies of the exponential cone (23) (with dummy variables $v_{abxy}=1$ for all $a,b,x,y$). Thus, we see that the problem of Eq. (21) is indeed an exponential conic program.

Appendix E Details of numerical investigations

E.1 Explicit form of the quantum distribution $\vec{P}_{\mathcal{Q}}$ considered

Here, we provide the explicit form of the ideal quantum distributions $\vec{P}_{\mathcal{Q}}$ employed in our numerical studies and an explicit quantum strategy realizing each of these correlations. First, we consider $\vec{P}_{\mathcal{Q}}=\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ with entries given by $\tfrac{1}{4}+(-1)^{a+b+xy}\tfrac{\sqrt{2}}{8}$. $\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ is known to violates maximally the CHSH CHSH Bell inequality:

[TABLE]

up to the limit of 2$\sqrt{2}$ allowed by quantum theory. $\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ can be realized by both parties locally measuring $\cos\tfrac{3\pi}{8}\sigma_{z}+\sin\tfrac{3\pi}{8}\sigma_{x}$ and $\cos\tfrac{7\pi}{8}\sigma_{z}+\sin\tfrac{7\pi}{8}\sigma_{x}$ for, respectively, input 0 and 1 on the shared state $|\Psi^{+}\rangle=\tfrac{1}{\sqrt{2}}(|01\rangle+|10\rangle)$.

Next, we consider $\vec{P}_{\mathcal{Q}}^{\text{\tiny 90\%CHSH}}$, which consists of a mixture of $\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ with the uniformly random distribution $\vec{P}_{\mathbb{I}}=\tfrac{1}{4}$. This allows us to gain some insight on how the regularization methods fare for noisy quantum distributions, which are more readily accessible in the laboratory. Explicitly, the entries of this distribution are:

[TABLE]

which can be realized by performing the measurements mentioned above on the mixed state $\rho=\tfrac{9}{10}|\Psi^{+}\rangle\!\langle\Psi^{+}|+\tfrac{1}{10}\tfrac{\mathbb{I}}{4}$ where $\tfrac{\mathbb{I}}{4}$ is the maximally mixed two-qubit state.

For the purpose of device-independent property estimations, it is known Moroder2013 that $\tilde{\mathcal{Q}}_{1}$ generally does not provide a tight estimate of, e.g., the amount of negativity Vidal2002 present  in the system. An example of this is given by the quantum distribution $\vec{P}_{\mathcal{Q}}=\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ which arises from both parties locally measuring the state $|\psi\rangle\simeq 0.91\,|00\rangle+0.42\,|11\rangle$ in the basis of $\vec{n}_{0}\cdot\vec{\sigma}$ and $\vec{n}_{1}\cdot\vec{\sigma}$, where $\vec{n}_{0}\simeq(0.26,-0.97)$ and $\vec{n}_{1}\simeq-(0.87,0,0.49)$ are Bloch vectors. Explicitly, using the matrix representation given just below Eq. (16), $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ reads as:

[TABLE]

where $\alpha_{00}\simeq 0.00$, $\beta_{00}=\gamma_{00}\simeq 0.01$, $\alpha_{01}=\alpha_{10}\simeq 0.00$, $\beta_{01}=\gamma_{10}\simeq 0.01$, $\beta_{10}=\gamma_{01}\simeq 0.21$, $\alpha_{11}\simeq 0.05$, $\beta_{11}=\gamma_{11}\simeq 0.16$ and $\epsilon_{xy}=1-\alpha_{xy}-\beta_{xy}-\gamma_{xy}$ for all $x,y\in\{0,1\}$. Note that $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ can be used to demonstrate more nonlocality with less entanglement Junge2011 ; Liang2011 ; Vidick2011 , as was achieved in Christensen2015 . In particular, being on the boundary of $\mathcal{Q}$, $\vec{P}_{\mathcal{Q}}^{\tau_{1.25}}$ maximally violates the $\tau=\tfrac{5}{4}$ version of the Bell inequality from Liang2011 :

[TABLE]

which is provably Christensen2015 satisfied by all finite-dimensional maximally entangled states whenever $\tfrac{1}{\sqrt{2}}+\tfrac{1}{2}\leq\tau\leq\tfrac{3}{2}$.

Finally, we also consider the correlation $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$ of Putz:NJP :

[TABLE]

which can be realized with both parties locally measuring $\sigma_{x}$ and $\sigma_{z}$ for, respectively, input 0 and 1 on the shared state $|\Psi\rangle=\tfrac{1}{\sqrt{3}}(|01\rangle+|10\rangle-|11\rangle)$. $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$ can be used to demonstrate the Hardy paradox Hardy1993 , as well as a violation of measurement-dependent locality Puetz2014 (MDL), e.g., via the following MDL inequality:

[TABLE]

where $h>l>0$. Note also that as opposed to $\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ which lies on the boundary of $\mathcal{Q}$, but strictly inside $\mathcal{N}$, $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$ lies on the boundary of both $\mathcal{Q}$ and $\mathcal{N}$. Experimental realizations of a correlation analogous to $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$ have been achieved in Aktas2015 ; Putz2016 .

E.2 Rate of convergence to the true distribution

We provide in Fig. 3 the plots of the mean value of the 1-norm deviation ${||\vec{P}_{\text{\tiny Reg}}(\vec{f})-\vec{P}_{\mathcal{Q}}||_{1}}$, between the regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ and the various $\vec{P}_{\mathcal{Q}}$ discussed above as a function of the number of trials $N_{\text{\tiny trials}}=10^{2},10^{3},\ldots,10^{10}$ for the regularization methods discussed in the main text and two additional regularization methods discussed in Appendix H. For ease of comparison, we also include in each of these figures the corresponding plot for $\vec{f}$.

Notice that from some basic numerical fitting, one finds that for all these methods, the mean value of ${||\vec{P}_{\text{\tiny Reg}}(\vec{f})-\vec{P}_{\mathcal{Q}}||_{1}}$, as with the mean value of ${||\vec{f}-\vec{P}_{\mathcal{Q}}||_{1}}$ diminishes at a rate of $\tfrac{1}{\sqrt{N_{\text{\tiny trials}}}}$. In addition, since the 1-norm upper bounds all other $p$-norms with $p$ being an integer greater than or equal to 2, our results of 1-norm deviation also upper bound the deviation when measured in terms of other $p$ norms.

E.3 Bias and mean squared errors of estimates

To gain further insight into the bias and the mean squared errors of various regularization methods, we provide in Fig. 4 and Fig. 5 our simulation results for the mean Bell-inequality violation and the corresponding mean squared error based on the regularized distributions as a function of the number of trials $N_{\text{\tiny trials}}=10^{2},10^{3},\ldots,10^{10}$. Our results clearly suggest that the bias in the Bell value obtained from $\vec{P}_{\Pi}(\vec{f})$ is essentially negligible,222As the projection method involves a linear transformation of $\vec{f}$, its bias is in theory identically zero. The nonzero bias that we observe in this case arises from the fact that our numerical simulations involve only a finite number of samples ($10^{4}$). whereas that obtained from $\vec{P}_{{\rm LS}}({\vec{f}})$ and $\vec{P}_{\text{\tiny ML}}({\vec{f}})$—for extremal $\vec{P}_{\mathcal{Q}}$— systematically underestimates (on average) the true value. However, as can be seen from the corresponding insets, such underestimations rapidly shrink with $N_{\text{\tiny trials}}$, diminishing at a rate of the order of $\tfrac{1}{\sqrt{N_{\text{\tiny trials}}}}$. For the case of nonextremal $\vec{P}_{\mathcal{Q}}$, such as $\vec{P}_{\mathcal{Q}}^{\text{\tiny 90\%CHSH}}$, we see that the bias present is essentially of the same order as that given by the projection method, which is basically negligible already for small $N_{\text{\tiny trials}}$.

Similarly, as can be seen from Fig. 4, the mean squared error rapidly decreases with $N_{\text{\tiny trials}}$ at a rate of the order of $\tfrac{1}{N_{\text{\tiny trials}}}$ in all the cases investigated. In particular, it is worth noting that for the case of $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$, the mean squared error present for the ML method is approximately three orders of magnitude less than all those given by the other methods. This superiority of the ML method over the others is, to some extent, anticipated from the fact that the KL divergence is superior as a statistical distance over, e.g., the total variation distance in terms of discriminating probability distribution that contains zero entries, such as $\vec{P}_{\mathcal{Q}}^{\text{\tiny MDL}}$ (see page 28 of the preprint version of vanDam2005 for a discussion).

In general, the bias nature of the few physical point estimators considered in Fig. 4 can be rigorously shown. See Appendix F.2 for a proof.

Appendix F Proof of certain properties of the estimators

F.1 Uniqueness of estimators

We give here a proof that the output of certain regularization methods $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ is determined uniquely by $\vec{f}$. To this end, we first recall from Theorem 8.3 of Beck2018 that the minimizer of a strictly convex function over a convex set is unique. To see that the LS method provides a unique estimate, let us first note that the Euclidean norm squared $(\|\vec{x}\|_{2})^{2}$ is strictly convex in $\vec{x}$, since its Hessian is two times the identity matrix. Moreover, $\min_{\vec{P}\in\tilde{\mathcal{Q}}_{\ell}}(\|\vec{f}-\vec{P}\|_{2})^{2}$ and $\min_{\vec{P}\in\tilde{\mathcal{Q}}_{\ell}}\|\vec{f}-\vec{P}\|_{2}$ share exactly the same set of minimizer(s). Thus, the output of the LS method with $\vec{P}_{\text{\tiny Reg}}\in\mathcal{C}$ for any convex set $\mathcal{C}\subseteq\overline{\mathcal{N}}$ is necessarily unique. Likewise, since $-\log_{2}(x)$ is a strictly convex function of $x$, the KL divergence from $\vec{P}\in\mathcal{C}$ to $\vec{f}$ is also strictly convex. In other words, the output of the regularization via both the LS method and the ML method is unique. In the main text, we have focused on $\mathcal{C}$ being some superset relaxation $\tilde{\mathcal{Q}}_{\ell}$ of $\mathcal{Q}$, but this uniqueness clearly applies also to the nonsignaling polytope $\mathcal{N}$, thus allowing one to show the uniqueness of the LS and the ML estimators even when the target set is $\mathcal{N}$ (see Appendix H).

F.2 Bias of estimators

We give here a proof of the biased nature of the physical estimator provided by the ML method and the LS method. To this end, we will first introduce the following definition for identifying correlations that may give rise to the same set of relative frequencies.

Definition 2**.**

Let $\vec{P}_{1}$ and $\vec{P}_{2}$ be two correlation vectors. We say that $\vec{P}_{1}$ and $\vec{P}_{2}$ are strictly nonorthogonal if for all input combinations, there is at least one combination of outcomes where the corresponding probability distributions of both $\vec{P}_{1}$ and $\vec{P}_{2}$ are nonvanishing.

Explicitly, in the bipartite scenario, the nonnegativity of probabilities means that $\vec{P}_{1}$ and $\vec{P}_{2}$ are strictly nonorthogonal if and only if $\sum_{a,b}P_{1}(a,b|x,y)P_{2}(a,b|x,y)>0$ for all $x,y$. The biased nature of $\vec{P}_{\text{\tiny ML}}({\vec{f}})$ and $\vec{P}_{{\rm LS}}({\vec{f}})$ is then an immediate consequence of Proposition 1, which we prove as follows.

Proof.

We will give a proof by contradiction. The steps follow closely those given in the proof of the main proposition of Schwemmer2015 . For simplicity, we give the proof in a finite bipartite Bell scenario. The generalization to more complicated but still finite Bell scenarios follows analogously. Suppose a Bell experiment is carried out with $\vec{P}_{j}$ governing the underlying joint distribution of measurement outcomes. If the experiment involves only a finite number of trials $N_{\text{\tiny trials}}$, and $\vec{P}_{j}$ is nondeterministic, one ends up with finite (nonsingleton) possibilities of relative frequencies $\vec{f}_{i}$ indexed by $i$. Let $\mathcal{F}_{j}=\{\vec{f}_{i}|q_{\vec{P}_{j}}(\vec{f}_{i})>0\}$ be the set of all such relative frequencies obtainable from $\vec{P}_{j}$, where $q_{\vec{P}_{j}}(\vec{f}_{i})\geq 0$ is the probability of observing the relative frequency $\vec{f}_{i}$ given that the underlying distribution is $\vec{P}_{j}$.

By assumption, one can find two strictly nonorthogonal extremal distributions of $\mathcal{C}$, say, $\vec{P}_{1}$ and $\vec{P}_{2}$, with $\vec{P}_{1}\neq\vec{P}_{2}$. Due to statistical fluctuations and the fact that $\vec{P}_{1}$ and $\vec{P}_{2}$ are strictly nonorthogonal, the set of relative frequencies simultaneously obtainable by both $\vec{P}_{1}$ and $\vec{P}_{2}$ is not empty, i.e., $\mathcal{F}_{1}\bigcap\mathcal{F}_{2}\neq\{\}$. Now, suppose that a regularization procedure gives an unbiased estimator. The expected distribution returned by the regularization scheme is then exactly the true distribution:

[TABLE]

where $\sum_{j}q_{\vec{P}_{j}}(\vec{f})=1$. In particular, since $\vec{P}_{1}\in\mathcal{C}$, we have

[TABLE]

By assumption, $\vec{P}_{1}$ is extremal in $\mathcal{C}$ and $\vec{P}_{\text{\tiny Reg}}(\vec{f})\in\mathcal{C}$ for all $\vec{f}\in\mathcal{F}_{1}$. Thus, the unbiased nature of $\vec{P}_{\text{\tiny Reg}}(\vec{f})$, see Eq. (33), implies that $\vec{P}_{\text{\tiny Reg}}(\vec{f})=\vec{P}_{1}$ for all $\vec{f}\in\mathcal{F}_{1}$.

Exactly the same line of reasonings can be applied to $\vec{P}_{2}$ to give $\vec{P}_{\text{\tiny Reg}}(\vec{f})=\vec{P}_{2}$ for all $\vec{f}\in\mathcal{F}_{2}$. But as argued above, there exists $\vec{f}=\vec{f}_{c}$ such that $\vec{f}_{c}\in\mathcal{F}_{1}$ and $\vec{f}_{c}\in\mathcal{F}_{2}$. This implies that $\vec{P}_{1}=\vec{P}_{\text{\tiny Reg}}(\vec{f}_{c})=\vec{P}_{2}$, which contradicts our assumption that $\vec{P}_{1}\neq\vec{P}_{2}$. Hence, given the premise that there exist extremal $\vec{P}_{1}$, $\vec{P}_{2}\in\mathcal{C}$ that are strictly nonorthogonal, the estimator $\vec{P}_{\text{\tiny Reg}}(\vec{f})$, which is constrained to be a member of $\mathcal{C}$, must be biased.

∎

To complete the proof that $\vec{P}_{\text{\tiny ML}}({\vec{f}})$ and $\vec{P}_{{\rm LS}}({\vec{f}})$ are biased, it suffices to choose $\vec{P}_{1}=\vec{P}_{\mathcal{Q}}^{\text{\tiny CHSH}}$ and $\vec{P}_{2}$ as the deterministic point with $P_{2}(a,b|x,y)=\delta_{a,x}\delta_{b,y}$. Indeed, both these quantum distributions are known to be extremal in $\mathcal{C}=\mathcal{Q}$ (as well as any of its relaxation $\tilde{\mathcal{Q}}_{\ell}$). Evidently, since no two extreme points of the nonsignaling polytope $\mathcal{N}$ in the simplest Bell scenario are strictly nonorthogonal, Proposition 1 cannot be used to show that a regularization method having $\mathcal{N}$ as its target is biased. Nonetheless, we conjecture that such estimators (which include all the other unique estimators given in Table 1) must also be biased.

Appendix G Device-independent negativity estimation, optimized witnesses and Bell inequality

In the scenarios we studied, the sets $\tilde{\mathcal{Q}}_{1}$, $\tilde{\mathcal{Q}}_{2}$, …, $\tilde{\mathcal{Q}}_{4}$ correspond to outer approximations of the quantum set $\mathcal{Q}$ using moment matrices of (local) level $\ell=1,2,\ldots$ (see Moroder2013 ; Vallins2017 for details). As shown in Moroder2013 , the same relaxations provide a way to lower bound the amount of entanglement present in a quantum system. We discuss here a geometrical formulation of their method. Let $\mathcal{Q}^{\leq\nu}$ be the set of all correlations $\vec{P}_{\mathcal{Q}}$ obtained from Born’s rule using a state $\rho$ of maximal negativity $\nu$. Following Moroder2013 , we write $\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}\supseteq\mathcal{Q}^{\leq\nu}$ as the corresponding semidefinite relaxation of level $\ell$. Then, given correlations $\vec{P}_{\mathcal{Q}}$ and approximation level $\ell$, a lower bound on the negativity of $\rho$ can be obtained by computing the supremum $\nu$ such that $\vec{P}_{\mathcal{Q}}\not\in\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}$.

What are the requirements on the distributions used as input to the negativity estimation algorithm? We first remark that $\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}$ is a subset of the nonsignaling set, so the negativity estimation algorithm can never be performed on the relative frequency $\vec{f}$, but should instead employ one of the regularized distributions $\vec{P}_{\text{\tiny Reg}}(\vec{f})$. Secondly, the set $\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}$ is a subset of $\tilde{\mathcal{Q}}_{\ell}$, so that $\vec{f}$ should be regularized to a target set $\tilde{\mathcal{Q}}_{\ell^{\prime}}$ with $\ell^{\prime}\geq\ell$. Otherwise, we run the risk of having $\vec{P}_{\text{\tiny Reg}}\in\tilde{\mathcal{Q}}_{\ell^{\prime}}\setminus\tilde{\mathcal{Q}}_{\ell}$ and the semidefinite program will turn out to be infeasible. In practice, this discrepancy could happen even by regularizing to $\ell^{\prime}=\ell$ due to insufficient numerical precision, in which case a small amount of white noise can be added to restore feasibility. For both $\ell=1,2$, this typical fraction of white noise is found to be of the order of 10*-8* or less.

In our negativity estimation shown in Fig. 2 of the main text, we have employed $\tilde{\mathcal{Q}}_{2}$ as our approximation to the quantum set $\mathcal{Q}$. This may seem unnecessary as $\tilde{\mathcal{Q}}_{1}$ Navascues2015 is known to be a pretty good approximation of the quantum set. However, as in the case of obtaining a negativity bound from the CHSH Bell-inequality violation Moroder2013 , one finds that $\tilde{\mathcal{Q}}_{1}$ generally does not provide a tight negativity bound. This becomes evident by plotting $\tilde{\mathcal{Q}}_{\ell}$ and $\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}$ around the distribution $\vec{P}_{\mathcal{Q}}^{\tiny\tau_{1.25}}$, as shown in Fig. 6. There, we see that although $\vec{P}_{\mathcal{Q}}^{\tiny\tau_{1.25}}$ lies visually on the boundary of both $\tilde{\mathcal{Q}}_{1}$ and $\tilde{\mathcal{Q}}_{2}$, the approximations $\tilde{\mathcal{Q}}_{\ell}^{\leq\nu}$ of the bounded negativity sets differ significantly, to the point that, in the considered slice, the approximation level $\ell=1$ is unable to certify any point with a negativity higher than $0.3$.

Let us also clarify here the connection between the optimized device-independent witness obtained by running a DI algorithm and an ordinary Bell inequality. By solving an SDP analogous to that given in Moroder2013 with a regularized distribution $\vec{P}_{\text{\tiny Reg}}(\vec{f})$, the dual of the SDP gives an optimized device-independent witness (a Bell-like inequality) of the form:333Extracting this witness follows a procedure analogous to that given in Appendix A of Schwarz2016 .

[TABLE]

where $\alpha\geq 0$ and $\mathcal{S}_{\alpha}$ are some fixed numbers, and $\sum_{a,b,x,y}\beta^{xy}_{ab}P_{\text{\tiny Reg}}(a,b|x,y)>\mathcal{S}_{\alpha}$.

Inequality (34) is an (optimized) device-independent negativity witness in the sense that for any (nonsignaling) $\vec{P}$, if it gives rise to the left-hand-side of Eq. (34) a value that is within the limit allowed by quantum theory and is greater than $\mathcal{S}_{\alpha}$, then whatever quantum state gives rise to $\vec{P}$ must have a negativity greater than $\alpha$. Thus, assuming that $\vec{P}_{\text{\tiny Reg}}(\vec{f})$ is a reasonable estimate of the underlying distribution, we see that it gives a negativity estimate of the underlying state that is at least $\alpha$.

Alternatively,  Moroder2013 provides a way to estimate negativity, in a device-independent manner, by starting from a given Bell-inequality violation: given any vector of coefficients $\beta_{ab}^{xy}$, one can compute, for each value $\mathcal{S}=\sum_{abxy}\beta_{ab}^{xy}P(a,b|x,y)$ a corresponding lower bound $\alpha$ on the negativity. For example, when one uses the coefficients of the CHSH inequality as $\beta_{ab}^{xy}$, we obtain the negativity bound mentioned in the main text: N$(\rho)\geq\tfrac{\mathcal{S}_{\text{\tiny CHSH}}-2}{4\sqrt{2}-4}$ for $\mathcal{S}_{\text{\tiny CHSH}}\in[2,2\sqrt{2}]$, which gives a monotonous relation between $\alpha$ and $\mathcal{S}$. Note that only when $\mathcal{S}_{\text{\tiny CHSH}}$ is greater than the local bound 2 one can hope to obtain a nonzero lower bound on negativity. Also, $\mathcal{S}_{\text{\tiny CHSH}}>2\sqrt{2}$ is not quantum realizable and thus cannot be used for the device-independent estimation of negativity.

In the same spirit, any vector of coefficients $\beta_{ab}^{xy}$ can be converted to a family of witnesses. In this process, one can even reuse the optimized Bell inequality given in  (34), providing a way to interpret values of $\mathcal{S}\geq\mathcal{S}_{\alpha}$ for which the coefficients $\beta_{ab}^{xy}$ were originally computed.

Appendix H Some other plausible regularization methods and their properties

Here, we discuss a few other plausible regularization methods with a target set $\mathcal{C}\in\{\mathcal{Q},\mathcal{N}\}$. Properties of these methods are summarized in Table 1.

H.1 Nearest quantum or nonsignaling approximation via $p$-norms

Obviously, one can regularize a given $\vec{f}$ to $\mathcal{Q}$ (resp. $\mathcal{N}$) by determining the nearest quantum (resp. nonsignaling) approximation NQA (resp. NNA) of $\vec{f}$ to $\mathcal{Q}$ (resp. $\mathcal{N}$) with a metric induced by any of the $p$-norms:

[TABLE]

When $p=2$, NQA2 gives the LS method described in Appendix C. By the same token, we shall refer to NNA2 of Eq. (35b) as LSN. We can rewrite it as the following second-order cone program (SOCP) Boyd2004Book :

[TABLE]

where inequalities in the last line are positivity constraints used to define the nonsignaling polytope. Hence, the LSN regularization can be also efficiently computed using an SOCP solver. Note that the inequality constraints in the last line of Eq. (36) can be replaced by imposing the nonsignaling constraints of Eq. (1). The regularization method of LSN is then evidently a least-square minimization problem with linear equality constraints. This regularization method has previously been implemented in Cavalcanti:2015aa as part of their data analysis.

For the case of $p=1$ (previously considered in Schwarz2016 ) or $p=\infty$, the optimization problem of Eqs. (35a) and  (35b) can be cast, respectively, as a semidefinite program and a linear program. For both values of $p$, one can find easily an example of $\vec{f}$ where some \vec{P}_{\text{\tiny NQA{}{p}}}(\vec{f}) [\vec{P}_{\text{\tiny NNA{}{p}}}(\vec{f})] is Bell-inequality violating but some other is not. For explicit examples, see the Supplemental Material MatFile . Similarly, the variant of NNA1 employed in Bernhard2014 is known to give nonunique estimators Schwarz:Private .

H.2 Minimizing the KL divergence to $\mathcal{N}$

As with the LSN method, one can also consider minimizing the KL divergence from the nonsignaling polytope $\mathcal{N}$ to some given relative frequency $\vec{f}$

[TABLE]

We shall refer to the corresponding regularization method as the MLN method. Note that, as with the ML method, performing the MLN regularization method amounts to solving a conic program, which can be achieved using, e.g., the SCS solver.

Although not explicitly discussed as a regularization method, the MLN method has been employed in Zhang2011 and was noted to help in the analysis of the hypothesis testing of local causality, as discussed in Appendix 2 of Zhang2011 [see Eqs. (A1) and (A2) therein]. Conceivably, the ML method may help in the analysis of Bell tests further.