Exploring the framework of assemblage moment matrices and its applications in device-independent characterizations

In a recent work [Phys. Rev. Lett. 116, 240401 (2016)], a framework known by the name of"assemblage moment matrices"(AMMs) has been introduced for the device-independent quantification of quantum steerability and measurement incompatibility. In other words, even with no assumption made on the preparation device nor the measurement devices, one can make use of this framework to certify, directly from the observed data, the aforementioned quantum features. Here, we further explore the framework of AMM and provide improved device-independent bounds on the generalized robustness of entanglement, the incompatibility robustness and the incompatibility weight. We compare the tightness of our device-independent bounds against those obtained from other approaches. Along the way, we also provide an analytic form for the generalized robustness of entanglement for an arbitrary two-qudit isotropic state. When considering a Bell-type experiment in a tri- or more-partite scenario, we further show that the framework of AMM provides a natural way to characterize a superset to the set of quantum correlations, namely, one which also allows post-quantum steering.

I. INTRODUCTION
By using a Bell-nonlocal [1,2] resource, such as an entangled pure quantum state, one can generate correlations between measurement outcomes which do not obey the principle of local causality [3], beating our intuitive understanding of nature. To date, convincing experimental demonstrations of Bell-nonlocality (hereafter abbreviated as nonlocality) have been achieved in a number of different physical systems (see, e.g., Refs. [4][5][6][7]).
Operationally, nonlocality enables one to perform some tasks that are not achievable in classical physics, including quantum cryptography [8], randomness generation [9,10], reduction of communication complexity [11] etc. For example, using nonlocal correlations, the task of quantum key distribution [8] can be achieved [12] even when one assumes nothing about the shared quantum resource or the measurement apparatuses. Since then, several quantum information tasks have been proposed within this black-box paradigm (see [13][14][15] and references therein) -forming a discipline that has come to be known as device-independent (DI) quantum information.
Another peculiar feature offered by quantum theory is steering [16] -the fact that one can remotely steer the set of conditional quantum states (called an assemblage [17]) accessible by a distant party by locally measuring a shared entangled state. This intriguing phe- * shin-liang.chen@mpq.mpg.de † costantino.budroni@oeaw.ac.at ‡ ycliang@mail.ncku.edu.tw § yuehnan@mail.ncku.edu.tw nomenon was revisited in 2007 by Wiseman, Jones, and Doherty [18]. In turn, their mathematical formulation forms the basis of a very active field of research (see, e.g., Refs. [19][20][21][22][23] and references therein) and has given rise to the so-called one-sided DI quantum information [24]. To exhibit nonlocality or to demonstrate the steerability of a quantum state, it is necessary to employ incompatible measurements [25]. In particular, among existing formulations of such measurements [26][27][28], any measurements that are incompatible-in the sense of being nonjointly-measurable [29]-can always be used [30,31] to demonstrate the steerability of some quantum states. In fact, the incompatibility robustness [32]-a quantifier for measurement incompatibility-has even been shown to be lower bounded [33,34] by the steering robustness [21] -a quantifier for quantum steerability.
In the context of DI quantum information, a moment matrix, i.e., a matrix composed of a set of expectation values of observables, is known to play a very important role. In particular, the hierarchy of moment matrices due to Navascués, Pironio, and Acín (NPA) [35] not only has provided the only known effective characterization (more precisely, approximation) of the quantum set, but also has found applications in DI entanglement detection [36,37], quantification [33,38,39], dimensionwitnessing [13,40,41], self-testing [42,43] etc. Similarly, some other variants [17,44] of the NPA hierarchy have also found applications in the context of one-sided DI quantum information. In Appendix A, we summarize in Table A some of the hierarchy of moment matrices that have been considered in (one-sided) DI quantum information.
Inspired by the moment matrices considered in Refs. [17,38], a framework known by the name of as-semblage moment matrices (AMMs) was proposed in Ref. [33]. As opposed to previous considerations, a distinctive feature of AMM is that the moment matrices considered consist of expectation values only for subnormalized quantum states (specifically, the assemblage induced in a steering experiment). This unique feature makes AMM a very natural framework for the DI quantification of steerability, and consequently the DI quantification of measurement incompatibility as well as the DI quantification of entanglement robustness, and its usefulness in certain quantum information tasks.
In this paper, we further explore the relevance of AMM for DI characterizations. We begin in Sec. II by reviewing the concept of moment matrices considered in DI quantum information. Then, we recall from Ref. [33] the framework of AMM in Sec. III. After that, we discuss the applications of AMM in DI quantum information, specifically DI characterizations. In Sec. VI, we conclude with a summary results and outline some possibilities for future research.

II. MOMENT MATRICES WITHIN THE DEVICE-INDEPENDENT PARADIGM
Moment matrices, i.e., matrices of expectation values of certain observables, were first discussed in a DI setting by NPA in Ref. [35]. For our purposes, however, it would be more convenient to think about these matrices as the result of some local, complete-positive (CP) maps acting on the underlying density matrix, as discussed in Ref. [38]. To this end, consider two local CP maps Λ A and Λ B acting, respectively, on Alice's and Bob's system (ρ A and ρ B ): Together, when applied to a quantum state ρ AB , these local CP maps give rise to a matrix χ of expectation values which is a function of ρ AB , as well as the choice of {A i } and {B j }.
Consider now a bipartite Bell experiment where Alice (Bob) can freely choose to perform any of the n x (n y ) measurements, each giving n a (n b ) possible outcomes. In quantum theory, these measurement are described by positive-operator-valued measures (POVMs). Let {E A a|x } x,a and {E B b|y } y,b respectively denote the collection of POVM elements (also known as a measurement assemblage [21]) associated with Alice's and Bob's measurements, and let 1 1 be the identity operator. Then, if we let {A i } ({B j }) be the set of operators obtained by taking all ℓ-fold products of operators from , the corresponding moment matrix, cf. Eq. (2), is said [38] to be a moment matrix of local level ℓ (see also Ref. [45]). Note that for all ℓ ≥ 1, one can find in the corresponding moment matrix χ (ℓ) expectation values that are (at most) first order in E A a|x , E B b|y . From Born's rule, one finds that they correspond to the joint probability of Alice (Bob) observing outcome a (b) conditioned on she (he) performing the x-th (y-th) measurement, i.e., P (a, b|x, y)  [38,46,47] of superset approxi-mationsQ (ℓ) to the set of legitimate quantum correlations (denoted by Q) can be obtained by solving a hierarchy of semidefinite programs, each associated with a moment matrix of local level ℓ. Moreover, the hier-archyQ (1) Q (2) ... Q provably converges to Q, i.e.,Q (ℓ→∞) → Q (see also [46,47]). In performing this algorithmic characterization, since any POVM can be realized as a projective measurement (embedded in higher-dimensional Hilbert space [48]), without loss of generality one can thus set the uncharacterized {E a|x } a and {E b|y } b to be projectors for all x and y, such that E a|x E a ′ |x = δ a,a ′ E a|x and E b|y E b ′ |y = δ b,b ′ E b|y . In addition, one can further assume that each u i is a real number; see [38] for the detailed reasonings behind these simplifications. In Table I, we provide a summary of the various elements of χ (ℓ) DI in relation to the operators whose expectation values are to be evaluated.
In the DI paradigm explained above, all preparation devices and measurement devices are treated as uncharacterized (black) boxes. In contrast, consider now a situation where the measurements devices of one party, say, Bob, are fully characterized. Then, for every outcome a that Alice obtains when she performs the x-th measurement, Bob can in principle perform quantum state tomography to determine the corresponding quantum statê ρ a|x prepared on his end.
In quantum theory, if the shared quantum state is ρ AB and Alice's measurement assemblage is given by , thenρ a|x is simply the normalized version of the conditional state where tr A (.) refers to a partial trace over Alice's Hilbert space. Explicitly, if we denote by P (a|x) = tr(ρ a|x ), thenρ a|x = ρ a|x /P (a|x). Following Ref. [17], we refer to the set of conditional quantum states {ρ a|x } a,x ({ρ a|x } in short) as an assemblage.
In certain cases, instead of the usual quantum mechanical description, the preparation of an assemblage {ρ a|x } can be understood via a semiclassical model. Specifically, following Ref. [18], we say that an assemblage {ρ a|x } admits a local-hidden-state (LHS) model if there exists legitimate probability distributions P (λ), P (a|x, λ), and normalized quantum statesσ λ such that i.e., the observed assemblage is an average of quantum statesσ λ distributed to Bob over the common-cause distribution P (λ) and the local response function P (a|x, λ) on Alice's end. In this case, it is conventional to refer to the assemblage as being unsteerable. Otherwise, an assemblage {ρ a|x } that cannot be decomposed in the form of Eq. (5) is said to be steerable, as Alice can apparently steer the ensemble of quantum states at Bob's end with her choice of local measurements. There are several ways to quantify the degree of steerability of any given assemblage {ρ a|x }, e.g., the steerable weight [20], the steering robustness [21], the relative entropy of steering [22,49], the optimal steering fraction [23], consistent trace-distance measure [50] etc. In this paper, we would focus predominantly on the steering robustness (SR), defined [21] as the minimum (unnormalized) weight associated with another assemblage {τ a|x } so that its mixture with {ρ a|x } is unsteerable, i.e., 1 where D(a|x, λ) = δ a,λx , λ = (λ 1 , . . . , λ nx ), and σ λ is a subnormalized quantum state [σ λ = P (λ)σ λ , cf. Eq. (5)]. In the above formulation, we have made use of the fact that, in determining the existence of a decomposition in the form of Eq. (5), it suffices to consider deterministic P (a|x, λ) in the form just described.
A prominent advantage of SR is that, as with steerable weight [20], it can be efficiently computed as a semidefinite program (SDP) [by setting (1+t)σ λ as ρ λ in Eq. (6)]: From the dual of this SDP (see Ref. [51]), one finds that SR actually coincides with the optimal steering fraction, a steering monotone (based on optimal steering inequalities) introduced in Ref. [23]. Finally, as remarked by Piani and Watrous [21], SR can be given an operational meaning in terms of the (relative) success probability of some quantum information tasks (more on this below).

B. The framework of assemblage moment matrices
In a DI setting, single-partite probability distributions P (a|x), P (b|y) alone cannot be used to provide nontrivial characterizations of the underlying devices. This is because for one to arrive at any nontrivial statement, the observed correlation P obs must also violate a Bell inequality [13,14]. Since single-partite probability distributions alone do not reveal any correlation between the measurement outcomes of distant parties, they cannot possibly violate any Bell inequalities. Following this reasoning, it may seem the case that moment matrices associated with single-partite density matrices are also useless for DI characterizations.
While this intuition is true for normalized singlepartite density matrices, the same cannot be said when it comes to an assemblage, which consists only of subnormalized density matrices that arise in a steering experiment. Specifically, for each combination of outcome a and setting x, applying the local CP map of Eq. (1) to the conditional state ρ a|x , cf. Eq. (4), gives rise to a matrix of expectation values: where {B i } are again operators formed from the product of {1 1} ∪ {E B b|y } y,b . When the set {B i } involves operators that are at most ℓ-fold product of Bob's POVM elements, the collection of matrices in Eq. (8) are said [33] to be the assemblage moment matrices (AMMs) of level ℓ, and we denote each of them by Indeed, as with the moment matrices introduced in Sec. II, all entries of P obs can be identified with entries in these single-partite moment matrices. For example, by using Eq. (4) in Eq. (8) and choosing an entry in As an explicit example, consider the ℓ = 1 AMMs with Although the known data in these AMMs are P obs , we shall write ] to emphasize that the underlying moment matrices are induced by an assemblage.
Eq. (8) we have that for each a and x: (9) For DI characterizations, we then write this matrix (for a fixed value of a and x) as: where we have made use of the simplification mentioned above and expressed the experimentally inaccessible expectation value as: with u (a,x) i ∈ R (see Ref. [38]).

IV. DEVICE-INDEPENDENT APPLICATIONS
Having recalled from Ref. [33] the AMM framework, we are now in a position to further explore the framework for DI characterizations.

A. Quantification of steerability
As was already noted in our previous work [33], a DI lower bound on SR forms the basis of a couple of DI applications based on the AMM framework. For completeness and for comparison with the improved lower bound that we shall present in Sec. IV D, we now explain how a DI lower bound on SR can be obtained by relaxing the optimization problem given in Eq. (7), as was proposed in Ref. [33].
To this end, let us emphasize once again that in the DI paradigm, one does not assume any knowledge (e.g., the Hilbert space dimension) of quantum states ρ λ and ρ a|x . However, if the constraints of Eq. (7) hold, it must be the case that even upon the application of the local CP map given in Eq. (8), the constraints-which demand the positivity of certain matrices-would still hold. At the same time, notice that each tr (ρ λ ) appearing in the objective function of Eq. (7) can still be identified as a specific entry, denoted by χ (ℓ) DI [ρ λ ] tr in the AMM. For example, in the AMM given in Eq. (9), the trace of the underlying matrix ρ a|x is given by the upper-left entry of the matrix. Putting all these together, we thus see that a DI lower bound on SR can be obtained by solving the following SDP: As explained above, Eq. (12b) and Eq. (12c) follow by applying the CP map of Eq. (8) to the constraints of Eq. (7). However, by themselves, physical constraints (including normalization, positivity and consistency) associated with the assemblage {ρ a|x } may be violated and thus have to be separately enforced in Eq. (12d) and (12e). Empirical observation enters at the level of observed correlation in Eq. (12f), i.e., by matching entries in the AMM with the empirical data summarized in P obs . Instead of Eq. (12f), a (weaker) lower bound can also be obtained by imposing an equality constraint of the form a,b,x,y β x,y a,b P (a, b|x, y) =Î β whereÎ β is the observed value of a certain Bell function specified by real coefficients β x,y a,b . Moreover, notice that if we have access to the observed probabilities P obs , the condition a χ (ℓ) DI [ρ a|x ] tr = 1 is automatically satisfied, as it amounts to the condition a P (a|x) = 1. On the other hand, this is not the case if we have access only to the Bell functionÎ β .
Importantly, the constraints of Eq. (12d) and Eq. (12e) do not necessarily single out {ρ a|x } as the underlying assemblage; neither do Eq. (12b) and Eq. (12c) entail the constraints of Eq. (7). The above optimization problem is thus a relaxation of that given in Eq. (7). For concreteness, let us denote the optimum of Eq. (12) by SR A→B DI,ℓ (P obs ) and that obtained for some observed Bell violation as SR A→B DI,ℓ (Î), then for all ℓ ≥ 1, thus giving the desired DI lower bound on SR({ρ a|x }) [see Sec. IV C 2 and also Ref. [34] for alternative approaches for bounding SR({ρ a|x }).] B. Quantification of the advantage of quantum states in subchannel discriminations From Eq. (13), one can also quantitatively estimate the usefulness of certain steerable quantum states in a kind of subchannel discrimination problem (see Ref. [21] and references therein). To this end, letΛ = a Λ a be a quantum channel (a trace-preserving CP map) that can be decomposed into a collection of subchannels {Λ a } a , i.e., a family of CP maps Λ a that are each trace nonincreasing for all input states ρ. Following Ref. [21], we refer to this collection of subchannels as a quantum instrument I = {Λ a } a . An example of I consists in performing measurement on the input state with "a" labeling the measurement outcome.
In its primitive form, a subchannel discrimination problem concerns the following task: input a quantum state ρ into the channelΛ and determine, for each trial, the actual evolution (described by Λ a ) that ρ undergoes by performing a measurement on Λ a [ρ]. For an input quantum state ρ, if we denote by {G a } a the POVM associated with the measurement on the output of the channel, then the probability of correctly identifying the subchannel Λ a is given by For any given quantum instrument I, the maximal probability of correctly identifying the subchannel is then obtained by maximizing the above expression over the input state ρ and the POVM {G a ′ } a ′ , i.e., where we use NE to signify "no entanglement" in the above guessing probability expression. In Refs. [21,52], the authors considered a situation where the input to the channel is a part of an entangled state ρ AB (B is the part that enters the channel) and where a measurement on the output I A ⊗ Λ B a [ρ AB ] is allowed. Suppose now that the final measurement is restricted to be separable across A and B, but allowed to be coordinated by one-way classical communication [21] (one-way LOCC) from B to A, i.e., taking the form of where Then, it was shown [21] that for any steerable quantum state ρ AB , there always exists an instrument I = {Λ a } a such that the corresponding guessing probability-after optimizing over measurements of the form given in Eq. (16)-exceeds p NE (I).
More precisely, let {G a ′ } a ′ take the form of Eq. (16). Then, for the initial state ρ AB , the corresponding guessing probability (after optimizing over such measurements) is The advantage of a steerable state ρ AB compared to unentangled resources in the subchannel discrimination task can then be quantified via the ratio of their success probabilities. In Ref. [21], this ratio was shown to be closely related to the SR A→B (ρ AB ), the steering robustness of the given quantum state ρ AB , defined as: Explicitly, since [21] sup and we can provide a DI lower bound on SR({ρ a|x }) via Eq. (13), it follows from Eq. (18) that we can also estimate in a DI manner the advantage of the measured state over unentangled resources for the task of subchannel discrimination.

C. Quantification of entanglement
The possibility to lower bound the entanglement of an underlying state in a DI setting was first demonstratedusing negativity [53] as the entanglement measure-in Ref. [38]. Subsequently, in Ref. [54], this possibility was extended to include the linear entropy of entanglement. In this subsection, we discuss how such a quantification can be achieved also for the generalized robustness of entanglement [55,56] defined as:

Via the approach of AMM
To obtain a DI lower bound on ER, we first remind that the set of unsteerable states (either from A to B, or from B to A) is a strict superset to the set of separable states. Hence, it is evident from Eq. (20) that (see also Ref. [21]) (21) It then immediately follows from Eq. (13) and Eq. (18) that for any assemblage on Bob's side {ρ a|x }, any assemblage on Alice's side {ρ b|y }, or any correlation P obs associated with these assemblages observed in a Bell experiment: which give the desired DI lower bounds on ER(ρ AB ).

Via the approach of nonlocal robustness
In Ref. [34], Cavalcanti and Skrzypczyk introduced, for any given correlation {P (a, b|x, y)}, a quantifier for nonlocality by the name of nonlocal robustness: where L and Q are, respectively, the sets of Bell-local and quantum correlations. Moreover, they [34] showed that the nonlocal robustness NR({P (a, b|x, y)}) for any correlation associated with an assemblage is a lower bound on the corresponding steering robustness, i.e., Hence, by using the first inequality of Eq. (21), we see that a DI lower bound on ER(ρ AB ) can also be obtained by computing NR(P obs ).

Via an MBLHG-based [38] approach
For comparison, let us mention here also the possibility for bounding ER(ρ AB ) based on the approach of Moroder et al. [38], abbreviated as MBLHG (see Sec. II). The idea is to first relax the separability constraint of Eq. (20) by the positive-partial-transposition constraint [57,58], thereby making the optimum of the following SDP, i.e., a lower bound on ER(ρ AB ); here, we use O TA to denote the partial transposition of operator O with respect to the Hilbert space of A. For a two-qubit state or a qubitqutrit state ρ AB , the result of Horodecki et al. [58] implies that the ER(ρ AB ) computed from Eq. (25) is tight. Next, by applying the local mapping of Eq. (1) to the linear matrix inequality constraints of Eq. (25), we obtain a further relaxation of Eq. (20)-and hence also a DI lower bound on ER(ρ AB )-by solving the following SDP: Note also that the second line of constraints on χ[ρ AB ] stems from the fact that we now no longer assume anything about the underlying state ρ AB , but only constraints of the form of Eq. (12f). Hereafter, we denote the optimum of Eq. (26) by ER DI,ℓ (P obs ).
In Table III, we summarize how a DI lower bound on ER(ρ AB ) can be obtained using the three approaches explained above.

Some explicit examples
To gain some insight on the tightness of the DI bounds provided by the aforementioned approaches, consider, for example, the isotropic states [59]: |i |i is the d-dimensional maximally entangled state, and 1 1 d 2 is the two-qudit maximally mixed state. It is known that these states are entangled if and only if v d > 1 d+1 . In Appendix D, we show that the generalized robustness of entanglement for these states are:  Table III). The best lower bounds obtainable for each approach are shown in Fig. 1 (for the lower bounds obtained for each approach in each Bell scenario, see Fig. 4). For visibilities less than v 2 ≈ 0.9314, the lower bounds obtained from the approach of AMM [33] and that of Certifiable DI lower bounds on the generalized robustness of entanglement (ER) for two-qubit isotropic states ρI,2(v2) based on various Bell-inequality-violating correlations P obtained from these states using the three approaches discussed in Sec. IV C (see text and Fig. 4 for further details). Bounds obtained from the approach of MBLHG [38], AMM [33], and of Ref. [34] are marked, respectively, using triangles (▽), squares ( ) , and crosses (+). For completeness, the actual value of ER[ρI,2(v2)] for each given value of visibility v2, cf. Eq. (28), is also included as a (red) solid line.
Ref. [34] seem to fit well with the expression ( √ 2v 2 − 1)( √ 2 − 1). But for greater values of v 2 , especially for v 2 0.9321, the AMM-based lower bounds appear to be somewhat tighter, and appear to fit nicely with the expression 2v 2 − √ 3. On the other hand, it is also clear from the Figure that the lower bounds ER DI,ℓ offered by the MBLHG-bsaed approach-which are well-represented by the expression 3 √ 2v2−1 √ 2−1 -considerably outperform the lower bounds obtained from the other two approaches.
As a second example, we consider the d = 3 case of Eq. (27) and the correlations leading to the optimal quantum violation of the I 2233 -Bell inequality [62] by these states. Our results are shown in Fig. 2. Again, as with the case shown in Fig, 2, the AMM approach appears to offer a somewhat tighter lower bounds than that of Ref. [34]. Also, the MBLHG-based approach again appears to give a much better lower bound on ER[ρ I (V 3 )] than the other two approaches. 3 In general, since the correlations P employed for a particular value of v 2 ∈ [ 1

FIG. 2. Certifiable DI lower bounds on the generalized robustness of entanglement (ER) for two-qutrit isotropic states
ρI,3(v3) based on the maximum I2233-Bell-inequality-violating correlations P obtained from these states using the three approaches discussed in Sec. IV C. Legends of the figure follow that of Fig. 1. ℓ in the legend denotes the level of the SDP hierarchy involved in the computation; a * is included as a superscript of ℓ whenever the next level of the hierarchy, ℓ + 1, gives the same SDP bound (within a numerical precision of the order of 10 −6 .

D. Quantification of measurement incompatibility
A collection of measurement, i.e., a measurement assemblage [21] {E a|x } a,x with a denoting the output and x the input, is said to be incompatible (not jointlymeasurable) whenever it cannot be written in the form where G λ 0, λ G λ = 1 1, and D(a|x, λ) can be chosen, without loss of generality, as D(a|x, λ) = δ a,λx [cf. Eq. (6) and the text thereafter]. In other words, a measurement assemblage is incompatible if there does not exist a joint measurement {G λ } providing all the outcome probabilities for any input.
The use of incompatible measurement is necessary to observe both nonlocality [25] and steering [30,31]. Moreover, steering and incompatibilty problems can be mapped from one into another [32], thus suggesting a measure of incompatibility, the incompatibility robustness (IR) introduced in Ref. [32]. In analogy to the steering robustness, IR may be computed by solving the fol-lowing SDP: where d is the dimension of {E a|x }. In Ref. [33], it has been proven that the steering robustness of a given assemblage {ρ a|x } is a lower bound on the incompatibility robustness of the steering equivalent observables [32] 4 which, in turn, is a lower bound on the incompatibility robustness of {E a|x }, namely The corresponding DI quantifier has then been discussed in Ref. [33]. An analogous observation has been made in Ref. [34], where Cavalcanti and Skrzypczyk also gave a lower bound on the degree of incompatibility, quantified by the incompatibility robustness of Alice's measurement assemblage {E A a|x } in a DI manner. In their work, they first introduced a modified quantifier of steerability, called the consistent steering robustness, defined as: {τ a|x } is a valid assemblage, Compared with Eq. (6), the consistent steering robustness needs more constraints, i.e., a τ a|x = a ρ a|x for all x. The above problem can also be formulated as the following SDP [by settingσ λ = (1 + t)σ λ and noting the non-negativity of τ a|x ]: Following an argument analogous to those in Ref. [33], one can straightforwardly prove that SR c ({ρ a|x }) = IR({B a|x }) for the steering equivalent observables {B a|x }. In fact, by a direct inspection of Eqs. (30) and (33) , and the fact that a ρ a|x = ρ B . To show the inverse transformation,it is sufficient to use the inverse of the above mappings.
In order to provide a DI lower bound on SR c ({ρ a|x }, the authors of Ref. [34] introduced a nonlocality quantifier [for a given correlation P] named consistent nonlocal robustness NR c (P): i.e., it calculates the minimum noise one has to mix into {P (a, b|x, y)} to make the mixture become local. {Q(a, b|x, y)} ∈ Q denotes {Q(a, b|x, y)} that has a quantum realization, cf. Eq. (3), and the last set of constrains requires the equivalence between the marginals of {P (a, b|x, y)} and {Q(a, b|x, y)}, similar to the case of the consistent steering robustness [see the last line of Eq. (32)]. Since the quantum set Q is not easily characterized, one can rather consider a supersetQ (ℓ) of Q by using the ℓ-th level of NPA hierarchy. In this way, one obtains a lower bound on NR c ({P (a, b|x, y)}) by solving the following SDP, which is reformulated from Eq. (34) [by setting q(λ) = 1+r r P (λ)]: Using the above quantifiers, Cavalcanti and Skrzypczyk proved [34] (36) which allows one to estimate the degree of incompatibility of Alice's measurements from the observed data P obs , i.e., in a DI manner.
Here, we would like to compare our method of lowerbounding incompatibility robustness with Eq. (36) by considering the example in Ref. [34]. That is, Alice and Bob share a pure partially entangled state |φ = cos θ|00 + sin θ|11 θ ∈ (0, π/4]. (37) For this state, optimal measurements for Alice and Bob giving the maximal violation of the Bell-Clauser-Horne (CH) inequality [65] are known analytically (see, e.g., Ref. [61]). One can, then, estimate the DI lower bounds on the incompatibility robustness of Alice's and Bob's measurements by using the above different approaches. The results are plotted in Fig. 3, together with our improved bound SR c,A→B DI,ℓ , that will be introduced below. With some attention, one observes a small but noticeable gap (of the order of 10 −3 or less) between SR DI,ℓ and NR c ℓ for some value of θ, even though we already employed the 5th level of AMM in our computation of SR DI,ℓ (while the computation of NR c was achieved using the 2nd level of the NPA hierarchy).
Such a gap may be explained by the fact that SR DI,ℓ does not take into account of the consistency condition a τ a|x = a ρ a|x , present in some form in NR c , which provides a better lower bound to IR. To improve our bound, we apply the the AMMs approach to SR c . Then, This optimization problem, however, is not in the form of an SDP since the third line contains quadratic constraints in the free variables. To circumvent this complication, we can relax the original problem by keeping, instead, only a subset of the original constraints, i.e., entries where i, j are those corresponding to [χ  Table IV summarizes the various approaches discussed above for the DI quantification of measurement incompatibility. From Fig. 3, we can see that SR c,A→B DI,ℓ provides a much better bound (in some instances, even tight bounds) on IR compared to SR DI,ℓ and NR c ℓ . On the other hand, it is also clear from the plots that, in these instances, SR DI,ℓ already provides a tight bound on the underlying SR.

V. MULTIPARTITE GENERALIZATION AND POST-QUANTUM STEERING
Evidently, the framework of AMM introduced in Sec. III can be generalized to a scenario with more than two parties. Below, we discuss this specifically for the tripartite scenario and explain how this leads to novel insights on the set of correlations characterized by the framework of AMM.
A. Steering in the tripartite scenario Following Ref. [66], let us consider a tripartite Belltype experiment where only Charlie has access to trusted (i.e., well-characterized) measurement devices. If we denote the shared quantum state by ρ ABC , the local POVM acting on Charlie's subsystem as E C c|z , then the analog of Eq. (3) reads as: while that of Eq. (4) reads as: It is straightforward to see from Eq. (42) that the assemblage {ρ C ab|xy } a,b,x,y (hereafter abbreviated as {ρ C ab|xy }) satisfy the positivity constraints and some no-signalinglike consistency constraints, i.e., ρ C ab|xy 0 ∀ a, b, x, y, tr As with the bipartite case, the assemblage {ρ C ab|xy } is said to admit an LHS model from A and B to C if there exists a collection of normalized quantum states {σ λ }, probability distribution P (λ), response functions P (a|x, λ) and P (b|y, λ) such that ρ C ab|xy = λ P (a|x, λ)P (b|y, λ)P (λ)σ λ for all a, b, x, y. Otherwise, the assemblage is said to be steerable from A and B to C.

B. AMMs in a tripartite scenario
To generalize the AMM framework to the aforementioned steering scenario, consider the analog of Eq. (1) that acts on the Hilbert space of Charlie's system ρ C : where {|i } ({|n }) are orthonormal bases vectors for the output (input) Hilbert spaceC (C) and C i , C j are some operators acting on C.
Specifically, for each combination of outcome a, b and setting x, y, applying the local CP map of Eq. (44) to the conditional state ρ C ab|xy gives rise to a matrix of expectation values: ab|xy ], is not fully determined. Following analogous procedure as that detailed in Sec. III B, one finds that the elements of the χ (ℓ) ab|xy ] fall under two categories: observable correlation (i.e., conditional probabilities) P obs (abc|xyz) 5 and unknown variables.
As an example, consider the steering scenario with binary input and output on Charlie's side such that Then, for all a, b, x, y, the first-level AMMs take the form of P obs (ab|xy) P obs (ab1|xy1) P obs (ab1|xy2) P obs (ab1|xy1) P obs (ab1|xy1) u abxy 1 P obs (ab1|xy2) u abxy where we have made use of the simplification mentioned in Sec. III and expressed the experimentally inaccessible expectation value as: with u abxy v ∈ R.
C. Correlations characterized by the AMM framework and post-quantum steering In Ref. [33], it was left as an open problem whether the set of correlations characterized by the AMM framework converges to the set of quantum distributions, i.e., the set of P that satisfy Born's rule. In this section, we show that in the tripartite scenario, the set of P allowed by demanding the positivity of AMMs-even in the limit of ℓ → ∞-generally cannot lead to the set of P that can be written in the form of Eq. (41).
To this end, we recall from Ref. [67] that there exists assemblage {ρ C ab|xy } satisfying Eq. (43) but not Eq. (42) for any ρ ABC and any local POVM {E A a|x }, {E B b|y }. The authors of Ref. [67] dubbed this phenomenon post-quantum steering. A simple example of this kind is given by 1} andρ is an arbitrary, but normalized density operator. Since the resulting marginal distribution P (a, b|x, y) is exactly that of a Popescu-Rohrlich box [68], we see that this assemblage cannot have a quantum realization. Now, note from our discussion in Sec. V B that if we start from an assemblage satisfying Eq. (43), the resulting AMMs are always positive semidefinite, and hence are compatible with the physical requirements imposed on AMMs. However, as mentioned above, there exists assemblage {ρ C ab|xy } satisfying Eq. (43) but which is not quantum realizable. We thus see that the AMM framework in the tripartite scenario, as described in Sec. V B, can, at best, lead to a characterization of the set of postquantum-steerable correlations, i.e., a superset of correlations satisfying Eq. (41) that also include, e.g., nonsignaling, but stronger-than-quantum marginal distributions between A and B.
On the other hand, it follows from the results of Refs. [69,70] that the phenomenon of post-quantum steering cannot occur in the bipartite scenario. Thus the problem of whether the set of correlations characterized by the AMM framework leads to the set of quantum distributions remains open in the bipartite scenario. Likewise, if one considers AMMs in a tripartite scenario based on one party steering the remaining two parties, the above argument does not apply either. As such, the problem whether one recovers-in the asymptotic limitthe quantum set, cf. Eq. (41), using the AMM framework remains open.

VI. CONCLUDING REMARKS
In this work, we have further explored and developed the AMM framework introduced in Ref. [33]. To begin with, we flashed out the details on how a DI bound on steering robustness (SR) provided by the AMM framework allows us to estimate the usefulness of an entangled state in the kind of subchannel discrimination problem discussed in Ref. [21].
We then went on to compare the DI bound on the generalized robustness of entanglement provided by the AMM framework against that given by the approach of Cavalcanti and Skrzypczyk [34]. Within our computational limit, the bounds of AMM appear to be slightly tighter than (or at least as good) those from the latter approach. In the process, we also offered another mean to bound the generalized robustness of entanglement from the data alone via the approach of Moroder et al. [38]. This last set of DI bounds turned out to be much stronger than that offered by the other two approaches. In these comparisons, we considered the two-qudit isotropic states where we also evaluated their generalized robustness of entanglement explicitly (see Appendix D).
Next, we compared the DI bound on the incompatibility robustness (IR) given by the AMM framework against that of Ref. [34]. In this case, the DI bounds offered by the AMM approach-based on bounding SR-do not perform as well compared with those of Ref. [34], which are based on bounding the underlying consistent steering robustness. Motivated by this difference, we then provided an alternative way to lower bound-in a DI manner-the consistent steering robustness via the AMM framework. This turned out to provide-as compared with the approaches just mentioned-much tighter (and in some instances even tight) DI bounds on the underlying IR. Even then, let us note that, in general, a tight DI bound on the underlying IR does guarantee the possibility to self-test the underlying measurements, as exemplified by the results of Ref. [71]. On a related note, we demonstrated in Appendix B how the AMM framework can be used to provide a DI lower bound on the steerable weight, and hence the incompatibility weight-another measure of incompatibility between different measurements.
We also briefly explored the framework in the tripartite scenario. This led to the observation that the AMM framework generally does not characterize the set of quantum correlations, but rather the set of correlations where the phenomenon of post-quantum steering is allowed. However, the problem of whether the set of correlations characterized by the AMM framework converges to the quantum set in the bipartite scenario, or in a multipartite scenario where one party tries to steer the states of the remaining parties remains unsolved.

Appendix A: Comparison of different moment-matrix approaches
The following table provides a comparison of different moment-matrix approaches. TABLE V. Different moment-matrix hierarchies in the context of (partially) device-independent (DI) characterizations. The first column labels the name of the author (or their initials if there is more than one author). The second column shows how the ℓ-th level construction of the ℓ-th level of each moment matrix (for simplicity, we provide an exemplification of the construction assuming two, or otherwise the minimal of parties where the hierarchy is applicable). Symbols like E A a|x ∈ L(HA) denotes the a-th POVM element of A's x-th measurement, while a POVM element like E AB a|x ∈ L(HAB) acts on the global Hilbert space of A and B. For brevity, bi (yi) etc. is introduced to denote one of the measurement outcomes (settings). M k x denotes the x-th observable of party-k.

Hierarchies
Moment matrix construction Accessible data Examples of applications • DI lower bound on negativity, Hilbert space dim.
• DI lower bound on steerability, measurement In this section we would like to show how the AMMs approach can be used to estimate the degree of steerability measured by the steerable weight [20]. Given an assemblage, one can always represent it as a (possibly trivial) convex mixture of a steerable assemblage {ρ S a|x } and an unsteerable assemblage {ρ US a|x }: The steerable weight (SW) is the minimal weight associated with ρ S a|x among all possible convex decompositions of {ρ a|x } according to Eq. (B1), i.e., SW = min(1 − µ). Importantly, SW that can be computed by solving the SDP [20]: To obtain a lower bound on SW({ρ a|x }) via the AMMs approach, we follow essentially the same steps in the derivation of Eq. (12) to obtain the following SDP: P (a, b|x, y) = P obs (a, b|x, y) ∀ a, b, x, y. (B3f) The lower bound on SW({ρ a|x }) obtained by solving Eq. (B3) will be denoted by SW A→B DI,ℓ (P obs ). As with steering robustness, SW provides [34] a lower bound on a measure of the incompatibility between measurements, called "incompatibility weight" [72] (IW). In analogy to Eq. (B1), IW is defined as the minimal weight associated with the non-jointly-mesurable assemblage in all possible convex decomposition of the given assemblage into a component that is jointly-measurable and one that is not. Since IW({E A a|x }) ≥ SW({ρ a|x }), SW A→B DI,ℓ (P obs ) provides, via the analog of Eq. (13) for SW, a DI lower bound on IW of the underlying measurement assemblage. In fact, as with the case of bounding IR by SR, our DI bounds on IW can be strengthened by introducing the additional constraints given by Eq. (39) in Eq. (B3). We will denote the corresponding DI bounds by SW c,A→B DI,ℓ (P obs ).
As an example, we may use the correlations P detailed in the caption of Fig. 3 to estimate the IW of Bob's measurement assemblage as a function of θ. Our numerical results show that, whenever the state is entangled (i.e., θ = 0), the DI bounds SW c,A→B DI,ℓ (P obs ) give essentially the value 1, which is also the incompatibility weights of the underlying measurement assemblages. based on various Bell-inequality-violating correlations P obtained from the two-qubit isotropic states ρI,2(v2) using the three approaches discussed in Sec. IV C. Specifically, these correlations P are obtained from Eq. (3) for ρI,2(v2) and measurements leading to their optimal CHSH-inequality violation (marked with "2222" in the legend and solid line in the plot), optimal I3322-Bell-inequality violation (marked with "3322" in the legend and dashed line in the plot), and the optimal elegant-Bell-inequality-violation (marked with "4322" in the legend and dotted line in the plot) of these states. Bounds obtained from the approach of MBLHG [38], AMM [33], and of Ref. [34] are marked, respectively, using triangles (▽), squares ( ) , and crosses (+). For completeness, the actual value of ER[ρI,2(v2)] for each given value of visibility v2, as given in Eq. (28), is also included as a (red) solid line. ℓ in the legend denotes the level of the SDP hierarchy involved in the computation; a * is included as a superscript of ℓ whenever the next level of the hierarchy, ℓ + 1, gives the same SDP bound (within a numerical precision of the order of 10 −6 ).

Bounds based on Bell-inequality violations
Instead of solving Eq. (26), we have also computed a relaxation thereof where we fixed only the value of speθ ≈ 0.75 • , we see a deviation of the order of 10 −4 . Alternative, tight DI lower bounds on IW for the same set of correlations can also be obtained via the consistent nonlocal weight introduced by Cavalcanti and Skrzypczyk [34].
where E xy := a,b=0,1 (−1) a+b P (a, b|x, y), our numerical results suggest 7 the following tight lower bound: 8 On the other hand, for the elegant Bell inequality, where δ (,) is the Kronecker delta function, we have instead the following tight lower bound on ER: where S E is the observed value of the elegant Bell inequality violation. For the I 3322 Bell inequality [see Eq. (19) of Ref. [62] for its explicit form], our numerical results up to level ℓ local = 3 are shown in Fig. 5. With the highest-level relaxation that we have considered, the minimal value of ER compatible an I 3322 violation in the interval of [0, 0.25] appears to be linear (up to a numerical precision of 10 −3 ); this lower bound can again be saturated by considering a two-qubit isotropic state in conjunction with its maximal quantum violation of the I 3322 inequality. However, we do not know if the nonlinear part of the curve where the Bell inequality is violated beyond 0.25 can be saturated. In general, the fact that these DI lower bounds are saturated by the two-qubit isotropic state means that the results for ER DI shown in Fig. 4 can also be obtained by fixing the value of the observed Bell inequality, instead of considering the full set of probability distributions P.  26) and a consideration of the hierarchy defined in Ref. [38].
is attained by choosing u d such that ρ I,d (u d ) is separable and is furthest away from ρ I,d (v d ) among all the separable ρ I,d (u d ). In other words, the optimization problem of Eq. (D1) is solved by setting u d = − 1 d 2 −1 . Equating the resulting mixture ω AB with an isotropic state that is barely separable, i.e., ρ I,d ( 1 d+1 ) gives: For a separable ρ I,d (v d ), its generalized robustness of entanglement is easily seen to be ER[ρ I,d (v d )] = 0. We thus arrive at the desired analytic expression of


extfieldsofstudy:
Mathematics
Computer Science

s2fieldsofstudy:
Computer Science
Mathematics

text:
Conditional Sparse $\ell_p$-norm Regression With Optimal Probability

We consider the following conditional linear regression problem: the task is to identify both (i) a $k$-DNF condition $c$ and (ii) a linear rule $f$ such that the probability of $c$ is (approximately) at least some given bound $\mu$, and $f$ minimizes the $\ell_p$ loss of predicting the target $z$ in the distribution of examples conditioned on $c$. Thus, the task is to identify a portion of the distribution on which a linear rule can provide a good fit. Algorithms for this task are useful in cases where simple, learnable rules only accurately model portions of the distribution. The prior state-of-the-art for such algorithms could only guarantee finding a condition of probability $\Omega(\mu/n^k)$ when a condition of probability $\mu$ exists, and achieved an $O(n^k)$-approximation to the target loss, where $n$ is the number of Boolean attributes. Here, we give efficient algorithms for solving this task with a condition $c$ that nearly matches the probability of the ideal condition, while also improving the approximation to the target loss. We also give an algorithm for finding a $k$-DNF reference class for prediction at a given query point, that obtains a sparse regression fit that has loss within $O(n^k)$ of optimal among all sparse regression parameters and sufficiently large $k$-DNF reference classes containing the query point.

Introduction
In areas such as advertising, it is common to break a population into segments, on account of a belief that the population as a whole is heterogeneous, and thus better modeled as separate subpopulations. A natural question is how we can identify such subpopulations. A related question arises in personalized medicine. Namely, we need existing cases to apply data-driven methods, but how can we identify cases that should be grouped together in order to obtain more accurate models?
In this work, we consider a problem of this sort recently formalized by Juba 2017: we are given a data set where each example has both a real-valued vector together with a Boolean-valued vector associated with it. Let Y denote the real-valued vector for some example, and let X denote the vector of Boolean attributes. We are also given some target prediction Z of interest, and we wish to identify a k-DNF c such that c( X) specifies a subpopulation on which we can find a linear predictor for Z, a, Y , that achieves small loss. In particular, it may be that neither a, nor any other linear predictor achieves small loss over the entire data distribution, and we simply wish to find the subset of the distribution, determined by our k-DNF c evaluated on the Boolean attributes of our data, on which accurate regression is possible. We refer to this subset of the distribution as a "segment". Juba proposed two algorithms for different families of loss functions, one for the sup norm (that only applies to O(1)-sparse regression) and one for the usual ℓ p norms. But, the algorithm for the ℓ p norm has a major weakness, namely that it can only promise to identify a condition c ′ ( X) that has a probability polynomially smaller than optimal. This is especially problematic if the segment had relatively small probability to begin with: it may be that there is no such condition with enough examples in the data for adequate generalization. By contrast, the algorithm for the sup norm recovers a condition c ′ with probability essentially the same as the optimal condition, but of course the sup norm is an undesirable norm to use for regression, as it is maximally sensitive to outliers. In general, ℓ p regression is more tolerant to outliers as p decreases, and the sup norm is roughly "p = ∞." In this work, in Theorem 10, we show how to give an algorithm for O(1)-sparse ℓ p regression with a two-fold improvement over the previous ℓ p regression algorithm (mentioned above) when the coefficient vector is sparse: 1. The new algorithm recovers a condition c ′ with probability essentially matching the optimal condition. 2. We also obtain a smaller blow-up of the loss relative to the optimal regression parameters a and optimal condition c, reducing the degree of this polynomial factor. More concretely, Juba's previous algorithm identifies k-DNF conditions, as does ours. This algorithm only identified a condition of probability Ω(µ/n k ) when a condition with probability µ exists, and only achieved loss bounded by O(ǫn k ) when it is possible to achieve loss ǫ. We improve the probability of the condition to (1 − η)µ for any desired η > 0, and also reduce the bound on the loss toÕ(ǫn k/2 ), or further toÕ(ǫg log log(n k )) if the condition has only g "terms" (ANDs of literals), in Corollary 12. This latter algorithm furthermore features a smaller sample complexity, that only depends logarithmically on the number n of Boolean attributes. We include a synthetic data experiment demonstrating that the latter algorithm can successfully recover more terms of a planted solution than Juba's sup norm algorithm.
We also give an algorithm for the closely related reference class regression problem in Theorem 14. In this problem, we are given a query point x * (along with values for the real attributes y * ), and we wish to estimate the corresponding z * . In order to compute this estimate for a single point, we find some k-DNF c ′ such that c ′ ( x * ) = 1 so x * is in the support of the distribution conditioned on c ′ ( X), and such that there are (sparse) regression parameters a ′ for c ′ that (nearly) match the optimal loss for sparse regression under any condition c with c( x * ) = 1. In this case, the condition c ′ is a reference class containing x * , obtaining the tightest possible fit among k-DNF classes containing x * . We think of this as describing a collection of "similar cases" for use in estimating the target value z * from the attributes y * . For this problem, we guarantee that we find a condition c ′ and sparse regression parameters a ′ that achieve loss within O(n k ) of optimal.
Our algorithms are based on sparsifiers for linear systems, for example as obtained for the ℓ 2 norm by Batson et al. 2012, and as obtained for non-Euclidean ℓ p norms by Cohen and Peng 2015 based on Lewis weights [Lewis, 1978]. The strategy is similar to Juba's sup norm algorithm: there, by enumerating the vertices of the polytopes obtained on various subsets of the data, he is able to obtain a list of all possible candidates for the optimal estimates of the sup-norm regression coefficients. Here, we can similarly obtain such a list of estimates for the ℓ p norms by enumerating the possible sparsified linear systems. Sparsifiers are often used to accelerate algorithms, e.g., to run in time in terms of the number of nonzero entries rather than the size of the overall input, though here we use that the size of a sparsified representation is small in a rather different way. Namely, the small representation allows for a feasible enumeration of possible candidates for the regression coefficients.
We note that this means that we also obtain an algorithm for O(1)-sparse ℓ p -norm regression in the list-learning model [Balcan et al., 2008, Charikar et al., 2017. Although Charikar et al. in particular are able to solve a large family of problems in the list-learning model, they observe that their technique can only obtain trivial estimates for regression problems. We show that once we have such a list, algorithms for the conditional distribution search problem [Juba, 2016, Zhang et al., 2017, Juba et al., 2018 can be generically used to extract a description of a good event c for the given candidate linear rule. Similar algorithms can also then be used to find a good reference class.

Relationship to other work
The conditional linear regression problem is most closely related to the problem of prediction with a "reject" option. In this model, a prediction algorithm has the option to abstain from making a prediction. In particular, El-Yaniv and Weiner 2012 considered linear regression in such a model. However, like most of the work in this area, their approach is based on scoring the "confidence" of a prediction function, and only making a prediction when the confidence is sufficiently high. This does not necessarily yield a nice description of the region in which the predictor will actually make predictions. On the other hand, Cortes et al. 2016 consider an alternative variant of learning with rejection that does obtain such nice descriptions; but, in their model, they assume that abstention comes at a fixed, known cost, and they simply optimize the overall loss when their space of possible prediction values includes this fixed-cost "abstain" option. By contrast, we can explicitly consider various rates of coverage or loss.
Our work is also related to the field of robust statistics [Huber, 1981, Rousseeuw andLeroy, 1987]: in this area, one seeks to mitigate the effect of outliers. That is, roughly, we suppose that a small fraction of the data has been corrupted, and we wish for our models and inferences to be relatively unaffected by this corrupted data. The difference between robust statistics and the conditional regression problem is that in conditional regression we are willing to ignore most of the data. We only wish to find a small fraction, described by some rule, on which we can obtain a good regression fit.
Another work that similarly focuses on learning for small subsets of the data is the the work by Charikar et al. 2017 in the list-learning model of Balcan et al. 2008. In that work, one seeks to find a setting of some parameters that (nearly) minimizes a given loss function on an unknown small subset of the data. Since it is generally impossible to produce a single set of parameters to solve this task -many different small subsets could have different choices of optimal parameters -the objective is to produce a small list of parameters that contains a near-optimal solution for the unknown subset somewhere in it. As we mentioned above, our approach actually gives a list-learning algorithm for O(1)-sparse ℓ p regression; but, we show furthermore how to identify a conditional distribution such that some regression fit in the list obtains a good fit, thus also solving the conditional linear regression task. Indeed, in general, we find that algorithms for solving the list-learning task for regression yield algorithms for conditional linear regression. The actual technique of Charikar et al. does not solve our problem since their approximation guarantee essentially always admits the zero vector as a valid approximation. Thus, they also do not obtain an algorithm for list-learning (sparse) regression either.
Another task similar to list-learning of linear regression is the problem of finding dense linear relationships, as solved by RANSAC [Fischler and Bolles, 1981]. But, these techniques only work in constant dimension. By contrast, although we are seeking sparse regression rules, this is a sparse fit in high dimension. As with list-learning, in contrast to the conditional regression problem, such algorithms do not provide a description of the points fitting the dense linear relationship.
Finally, yet another task similar to list-learning for regression is fitting linear mixed models [McCulloch andSearle, 2001, Jiang, 2007]. In this approach, one seeks to explain all (or almost all) of the data as a mixture of several linear rules. The guarantee here is incomparable to ours: in contrast to list-learning, the linear mixed model simply needs a list of linear rules that accounts for nearly all of the data; it does not need to find a list that accounts for all possible sufficiently large subsets of the data. So, there is no guarantee that any of the mixture components represent an approximation to the regression parameters corresponding to an event of interest. On the other hand, in linear mixed models, one does not need to give any description at all of which points should lie in which of the mixture components. In applications, one usually assigns points to the linear rule that gives it the smallest residual, but this may be less useful for predicting the values for new points.

Preliminaries
We now formally define the problems we consider in this work, and recall the relevant background.

The conditional regression and distribution search tasks
Formally, we focus on the following task: Definition 1 (Conditional ℓ p -norm Regression) Conditional ℓ p -norm linear regression is the following task. Suppose we are given access to i.i.d. examples drawn from a joint distribution over If a * has at most s nonzero components, and we requireˆ a likewise has at most s nonzero components, then this is the conditional s-sparse ℓ p -norm regression task.
As stressed by Juba 2017, the restriction to k-DNF conditions is not arbitrary. If we could identify conditions that capture arbitrary conjunctions, even for one-dimensional regression, this would yield PAC-learning algorithms for general DNFs. In addition to this being an unexpected breakthrough, recent work by Daniely and Shalev-Shwartz 2016 shows that such algorithms would imply new algorithms for random k-SAT, falsifying a slight strengthening of Feige's hypothesis [Feige, 2002]. We thus regard it as unlikely that any algorithm can hope to find conditions of this kind. Of the classes of Boolean functions that do not contain arbitrary conjunctions, k-DNFs are the most natural large class, and hence are the focus for this model.
The difficulty of the problem lies in the fact that initially we are given neither the condition nor the linear predictor. Naturally, if we are told what the relevant subset is, we can just use standard methods for linear regression to obtain the linear predictor; conversely, if we are given the linear predictor, then we can use algorithms for the conditional distribution search (learning abduction) task (introduced by Juba 2016, recalled below) to identify a condition on which the linear predictor has small error. Our final algorithm actually uses this connection, by considering a list of candidates for the linear predictors to use for labeling the data, and choosing the linear predictor that yields a condition that selects a large subset. In particular, we will use algorithms for the following weighted variant of the conditional distribution search task: Definition 2 (Conditional Distribution Search) Weighted conditional distribution search is the following problem. Suppose we are given access to examples drawn i.i.d. from a distribution over This task is closely related to agnostic conditional distribution search, which is the special case where the weights only take values 0 or 1. The current state of the art for agnostic conditional distribution search is an algorithm given by Zhang et al. 2017, achieving anÕ( √ n k )-approximation to the optimal error. That work built on an earlier algorithm due to Peleg 2007. Peleg's work already showed how to extend his original algorithm to a weighted variant of the problem, and we observe that an analogous modification of the algorithm used by Zhang et al. will obtain an algorithm for the more general weighted conditional distribution search problem we are considering here: Theorem 3 (Peleg 2007, Zhang et al. 2017) There is a polynomial-time algorithm for weighted conditional distribution search achieving anÕ(n k/2 (log b+ log 1/η + log log 1/δ))-approximation with probability 1 − δ given m = Θ b 3 µǫη 2 (n k + log 1 δ ) examples.
We note that it is possible to obtain much stronger guarantees when we are seeking a small formula. Juba et al. 2018 present an algorithm that, when there is a k-DNF with g terms achieving error ǫ, uses only m =Õ( gk µǫη 2 log n δ ) examples and obtains a k-DNF with probability at least (1−η)µ and errorÕ(ǫg log m). Although this is stated for the unweighted case (i.e., ǫ is a probability), it is easy to verify that since our loss is nonnegative and bounded, by rescaling the losses to lie in [0, 1], we can obtain an analogous guarantee for the weighted case: Theorem 4 (Juba et al. 2018) There is a polynomial-time algorithm for weighted conditional distribution search when the condition has g terms using m =Õ( bgk µǫη 2 log n δ ) examples achieving añ O(g log m)-approximation using a k-DNF withÕ(g log m) terms with probability 1 − δ.

Reference class regression
Using a similar approach, we will also solve a related problem, selecting a best k-DNF "reference class" for regression. In this task, we are not merely seeking some k-DNF event of probability µ on which the conditional loss is small. Rather, we are given some specific observed Boolean attribute values x * , and we wish to find a k-DNF condition c that is satisfied by x * solving the previous task. That is, c should have probability at least µ and our sparse regression fit has small conditional loss, conditioned on c. Naturally, the motivation here is that we have some specific point ( x * , y * ) for which we are seeking to predict z * , and so we are looking for a "reference class" c such that we can get the tightest possible regression estimate of z * from y * ; to do so, we need to take µ large enough that we have enough data to get a high-confidence estimate, and we need x * to lie in the support of the conditional distribution for which we are computing this estimate.
Definition 5 (Reference Class ℓ p -Regression) Reference class ℓ p -norm regression is the following task. We are given a query point x * ∈ {0, 1} n , target density µ ∈ (0, 1), ideal loss bound ǫ 0 > 0 approximation parameter η > 0, confidence parameter δ > 0, and access to i.i.d. examples drawn from a joint distribution over We wish to findˆ a ∈ R d with ˆ a 2 ≤ b and a reference class k-DNFĉ such that with probability and k-DNFs c * such that c * ( x * ) = 1 and Pr[c * ( X)] ≥ µ. If we also require bothˆ a and a * to have at most s nonzero components, then this is the reference class s-sparse ℓ p -norm regression task.
The selection and use of such reference classes for estimation goes back to work by Reichenbach 1949. Various refinements of this approach were proposed by Kyburg 1974 andPollock 1990, e.g., to choose the estimate provided by the highest-accuracy reference class that is consistent with the most specific reference class containing the point of interest x * . Our approach is not compatible with these proposals, as they essentially disallow the use of the kind of disjunctive classes that are our exclusive focus. Along the lines we noted earlier, it is unlikely that there exist efficient algorithms for selecting reference classes that capture arbitrary conjunctions, so k-DNFs are essentially the most expressive class for which we can hope to solve this task. Bacchus et al. 1996 give a nice discussion of other unintended shortcomings of disallowing disjunctions. A concrete example discussed by Bacchus et al. is the genetic disease Tay-Sachs. Tay-Sachs only occurs in two very specific, distinct populations: Eastern European Jews and French Canadians. Thus, a study of Tay-Sachs should consider a reference class at least partially defined by a disjunction over membership in these two populations.

ℓ p sparsification
Our approach is based on techniques for extracting low-dimensional sketches of small subspaces in high dimensions. The usual ℓ 2 norm uses much simpler underlying techniques, and we describe it first. The extension to ℓ p norms for p = 2 is obtained via Lewis weights [Lewis, 1978].

Euclidean sparsifiers
The kind of sketches we need originate in the work of Batson et al. 2012. Specifically, it will be convenient to start from the following variant due to Boutsidis et al. 2014: Lemma 6 (BSS weights [Boutsidis et al., 2014] where λ j denotes the jth largest eigenvalue. In particular, taking r = t/γ 2 for some γ ∈ (0, 1), we obtain that for the Proof: Observe that since the each ith row of [u] has unit norm, it must have an entry u i,j * that is at least 1/ √ d in magnitude. By the above argument, where notice in particular, the ith row of [s] ⊤ contributes at least ( √ s i u i,j * ) 2 ≥ s i /d to the norm.

Sparsifiers for non-Euclidean norms
It is possible to obtain an analogue of the BSS weights for p = 2 using techniques based on Lewis weights [Lewis, 1978]. Lewis weights are a general way to reduce problems involving ℓ p norms to analogous ℓ 2 computations. Cohen and Peng 2015 applied this to sparsification to obtain the following family of sparsifiers:

Theorem 8 (ℓ p weights [Cohen and Peng, 2015]) Given a d × t matrix [u]
there exists a set of r(p, t, γ) weights s 1 , . . . , s r such that for the d × r matrix [s] which has as its ith column s i e i , where r(p, t, γ) is asymptotically bounded as in Table 1. p Required dimension r p = 1 t log t γ 2 1 < p < 2 1 γ 2 t log(t/γ) log 2 log(t/γ) p = 2 t/γ 2 p > 2 log 1/γ γ 5 t p/2 log t Cohen and Peng also show how to construct the sparsifiers for a given matrix efficiently, but we won't be able to make use of this, since we will be searching for the sparsifier for an unknown subset of the rows.
We furthermore obtain an analogue of Lemma 7 for the ℓ p weights, using essentially the same argument: Lemma 9 Suppose the rows of [u] are orthonormal. Then the matrix [s] obtained by Theorem 8 has entries of magnitude at most (1 + γ) √ d.

The weighting algorithms for conditional and reference class regression
The results from sparsification of linear systems tell us that we can estimate the loss on the subset of the data selected by the unknown condition by computing the loss on an appropriate "sketch," a weighted average of the losses on a small subset of the data. The dimensions in Table 1 give the size of these sparsified systems, i.e., the number of examples from the unknown subset we need to use to estimate the loss over the whole subset. The key point is that since the predictors are sparse, the dimension is small; and, since we are willing to accept a constant-factor approximation to the loss, a small number of points suffice. Therefore, it is feasible to enumerate these small tuples of points to obtain a list of candidate sets of points for use in the sketch. We also need to enumerate the weights for these points, but since we have also argued that the weights are bounded and we are (again) willing to tolerate a constant-factor approximation to the overall expression, there is also a small list of possible approximate weights for the points. The collection of points, together with the weights, gives a candidate for what might be an appropriate sketch for the empirical loss on the unknown subset. We can use each such candidate approximation for the loss to recover a candidate for the linear predictor. Thus, we obtain a list of candidate linear predictors that we can use to label our data as described above. More precisely, our algorithm is as shown in Algorithm 1. In the following, let Π d 1 ,...,ds denote the projection to coordinates d 1 , d 2 , . . . , d s .
In our proof of this theorem, we will find it convenient to use the Rademacher generalization bounds for linear predictors (note that x → |x| p is pb p−1 -Lipschitz on [−b, b]): Theorem 11 ( Bartlett andMendelson 2002, Kakade et al. 2009 b], and any δ ∈ (0, 1), let L p ( a) denote E[| a, Y − Z| p ], and for an an i.i.d. sample of size m letL p ( a) be the empirical loss 1 m m j=1 | a, y (j) − z (j) | p . We then have that with probability 1 − δ for all a with a 2 ≤ b, Note that although this bound is stated in terms of the ℓ 2 norm of the attribute and parameter vectors y and a, we can obtain a bound in terms of the dimension s of the sparse rule if we are given a bound B on the magnitude of the entries: b ≤ √ sB.
Proof of Theorem 10: Given that we are directly checking the empirical ℓ p loss before returning a and c, for the quoted number of examples m it is immediate by a union bound over the iterations that any a and c we return are satisfactory with probability 1 − δ. All that needs to be shown is that the algorithm will find a pair that passes this final check.
By Theorem 11, we note that it suffices to have b 2p γ 2p ǫ 2p (2pb + 2 ln(12/δ)) 2 examples from the distribution conditioned on the unknown k-DNF event c * to obtain that the ℓ p loss of each candidate for a is estimated to within an additive γǫ with probability 1 − δ/3. By Hoeffding's inequality, therefore when we draw m 0 examples, there is a sufficiently large subset satisfying c * with probability 1 − δ/3.
We let [u] be an orthonormal basis for span{(Π d 1 ,...,ds y (j) , z (j) ) : c * ( x (j) ) = 1, j ≤ m 0 } and invoke Lemma 6 for ℓ 2 or Theorem 8 for p = 2. In either case, there is some set of weights s 1 , . . . , s r 0 for a subset of r 0 coordinates j 1 , . . . , j r 0 such that for any v in the column span of [u], [s] ⊤ [u] v has ℓ p norm that is a 1 ± γ-approximation to the ℓ p norm of [u] v. In particular, for any a, observing v = [y] a − z is in the column span of [u] by construction, we obtain Now, we observe that we can discard weights (and dimensions) from s 1 , . . . , s r 0 of magnitude smaller than γ/r 1/p 0 , since for any unit vector v, the contribution of such entries to [s] ⊤ [u] v p p (recalling there are at most r 0 nonzero entries) is at most γ p . So we may assume the r ≤ r 0 remaining weights all have magnitude at least γ/r 1/p . Furthermore, if we round each weight to the nearest power of (1 + γ), this only changes [s] ⊤ [u] v| p p by an additional (1 ± γ) factor. Finally, we note that since (Π d 1 ,...,ds y (j) , z (j) ) has dimension s+1, Lemmas 7 and 9 guarantee that the magnitude is also at most (1 + γ) √ s + 1. Thus it indeed suffices to find the powers (q 1 , . . . , q r ) for our r examples j 1 , . . . , j r such that (1 + γ) q ℓ is within (1 + γ) of s ℓ , and the resulting set of weights will approximate the ℓ p -norm of every v in the column span to within a 1 + 3γ-factor. Now, when the loop in Algorithm 1 considers (i) the dimensions d * 1 , . . . , d * s contained in the optimal s-sparse regression rule a * (ii) the set of examples j * 1 , . . . , j * r used for the sparse approximation for these coordinates and (iii) the appropriate weights (1 + γ) q * 1 , . . . , (1 + γ) q * r , the algorithm will obtain a vector a that achieves a (1 + 3γ)-approximation to the empirical ℓ p -loss of a * on the same s coordinates. It then follows from Theorem 3 that with probability at least 1 − δ/3 over the data, WtCond will in turn return to us a k-DNF c with probability (1 − η)µ that selects a subset of the data on which a achieves an αǫ =Õ((1 + γ) √ n k (log b + log 1/η + log log 1/δ)ǫ) approximation to the empirical ℓ p loss of a * on c * . This choice of a and c passes the final check and is thus sufficient.
By simply plugging in the algorithm from Theorem 4 for WtCond, we can obtain the following improvement when the desired k-DNF condition is small: Corollary 12 For any constant s and γ > 0, r as given in Table 1 for t = (s + 1), and m = Θ (1+γ)bgk µǫη 2 log n δ + s log d + r log m 0 log(γ 1/p s/r) γ + m 0 examples, if there is a g-term k-DNF solution to the conditional s-sparse ℓ p -regression task, then the modified Algorithm 1 runs in polynomial time and solves the task with aÕ(g log m)-term k-DNF with α =Õ((1 + γ)ǫg log m).
Note that this guarantee is particularly strong in the case where the k-DNF would be small enough to be reasonably interpretable by a human user.
The extension to reference class ℓ p -norm regression proceeds by replacing the weighted condition search algorithm with a variant of the tolerant elimination algorithm from Juba 2016, given in Algorithm 2.
Proof: For convenience, let N ≤ (1 + log 1+η b/ǫ 0 )(1 + log 1+η 1/µ 0 ) denote the total number of iterations. Consider first what happens when the loop considers the largest µ ≤ Pr[c * ( X)]/(1 + η) and the smallest ǫ that is at least (1 + η) 2 ǫ * . On this iteration, for each term T of c * , we observe that So, since W is bounded by b, by a Chernoff bound, 1 . Since this is in turn at most ǫµ, T will be included in c on this iteration.
But similarly, for T not in c with E[W · T ( X)] > (1 + η)µǫ, the Chernoff bound also yields that j:T (x (j) )=1 w (j) ≥ ǫµm with probability 1 − δ/2 n ≤k N . By a union bound over all T ∈ c * and T not in c * with such large error, we see that with probability at least 1 − δ/2N , all of the terms of c * are included in c and only terms with E[W · T ( X)] ≤ (1 + η)µǫ are included in c. So, Furthermore, by yet another application of a Chernoff bound, c * is true of at least µm examples with probability at least 1 − δ/2N . Thus, with probability 1 − δ/N , after this iterationĉ is set to some k-DNF andǫ ≤ (1 + η) 2 max{ǫ * , ǫ 0 }. Now, furthermore, on every iteration, we see more generally that with probability 1 − δ/N , only . Thus, with probability 1 − δ overall, since we found above thatǫ ≤ (1 + η) 2 max{ǫ * , ǫ 0 }, we return a k-DNFĉ as claimed.

Experimental evaluation
To evaluate our algorithm's performance in practice, we two kinds of experiments: one using synthetic data with a planted solution, and another using some standard benchmark data sets from the LIBSVM repository [Chang and Lin, 2011].
Synthetic data. To generate the synthetic data sets, we first chose a random 2-DNF over Boolean attributes by sampling terms uniformly at random. We also fixed two out of d coordinates at random and sampled parameters a from a mean 0, variance σ 2 Gaussian for our target regression fit. For the actual data, we then sampled Boolean attributes X that, with 25% probability are a uniformly random satisfying assignment to the DNF and with 75% probability are a uniformly random falsifying assignment. We sampled Y from a standard Gaussian in R d . Finally, for the examples where X was a satisfying assignment, we let Z be given by the linear rule a, Y + ν where ν is a mean 0, variance σ 2 Gaussian; otherwise, if X was a falsifying assignment, Z was set to ν * where ν * is a standard mean 0, variance 1 Gaussian. (Parameters in Table 2.) Thus, conditioned on the planted formula, there is a small-error regression fit a-the expected error of a is |ν|. Off of this planted formula, the expected error for the optimal prediction z = 0 is |ν * | = 1. For these simple data sets, we made several modifications to simplify and accelerate the algorithm. First, we simply fixed all of the weights (of the form (1 + γ) q i ) to 1, i.e., q i = 0 for all i, since the leverage scores for Gaussian data (and indeed, most natural data sets) are usually small. Second, we set m 0 = 200 for the 1000-example data set and m 0 = 500 for the 5000-example data set. Third, we used the version of the algorithm for small DNFs considered in Corollary 12, that uses algorithm of Juba et al. 2018 (recalled in Theorem 4) for WtCond.
With these modifications, we ran Algorithm 1 for ℓ 2 regression with s = 2, γ = 1, and µ = 0.22 on the 1000-example data set. The algorithm was able to identify the planted 2-DNF together with the two components used in the regression rule in this case, and thus also we approximately identify the linear predictor. The mean squared error of the approximate predictor returned by our algorithm was 0.0149 on the planted 2-DNF (c.f. the ideal linear fit has variance-0.01 Gaussian noise added to it). However, we note that once the planted 2-DNF has been identified, a good regression fit can be found by any number of methods. We also ran Algorithm 1 for ℓ 2 regression with s = 2, γ = 1, and µ = 0.2465 on the 5000-example data set. The algorithm was able to identify a subset of 11 out of the 16 planted k-DNF terms, picking up about 98% of the planted k-DNF data. The mean squared error of the approximate predictor was 0.4221.
As a baseline, we also used the algorithm from Juba 2017 for conditional sparse sup-norm regression on our synthetic data. In order to obtain reasonable performance, we needed to make some modifications: we modified the algorithm along the lines of Algorithm 1, to only search over tuples of a random m 0 = 200 examples for the 1000-example data set and m 0 = 500 for the 5000-example data set in producing its candidate regression parameters. We also chose to take the parameters that achieved the smallest residuals under the sup norm, among those that obtained a condition satisfied by the desired fraction of the data. Again, we used s = 2, and for the norm bound we used a threshold of ǫ = 0.24 for the 1000-example data set and 1.1 for the 5000-example data set. (We considered a range of possible thresholds; when ǫ ≈ 1.5 on the σ 2 = 0.1 data, the sup norm algorithm starts adding terms outside the planted DNF, significantly increasing its error.) For the 1000-example data set, this algorithm was able to identify the planted 2-DNF. However, for the 5000-example data set, it only found 10 out 16 terms of the planted 2-DNF, picking up about 96% of the planted k-DNF data.
Real-world benchmarks. We also compared Algorithm 1 against selective regressors [El-Yaniv and Wiener, 2012] on some of the LIBSVM [Chang and Lin, 2011] regression data sets that were used to evaluate that work. We split each of those datasets into a training set (1/3 of the data) and test set (2/3 of the data). We generated Boolean attributes by choosing different binary splits (median, quartile) on the numerical features, excluding the target attribute. We then run our algorithm 1 on the training set to obtain a 2-DNF. Next, we filtered both the training and test data satisfying the 2-DNF, using the first one to train a new regression fit on the selected subset and the other as a holdout to estimate the error of this resulting regression fit. We compared this to the test error for selective regression. Risk-Coverage (RC) curves of the results are shown in Figure 1. We outperform the baseline on the Boston housing dataset, and we generally achieved lower error than the baseline for lower coverage (less than 0.5) and somewhat higher error for higher coverage for the other three data sets (Body fat, Space, and Cpusmall). Note that in addition, we obtained a 2-DNF that describes the subpopulation on all four, which is the main advantage of our method. Since selective regression first tries to fit the entire data set, and then chooses a subset where that fixed predictor does well, it is to be expected that it may miss a small subset where a different predictor can do better, but that its freedom to abstain on a somewhat arbitrary subset may give it an advantage at high coverage.

Directions for future work
There are several natural open problems. First, although sparsity is desirable, our exponential dependence of the running time (or list size) on the sparsity is problematic. (The sup norm regression algorithm [Juba, 2017] also suffered this deficiency.) Is it possible to avoid this? Second, our algorithm for reference class regression has an O(n k ) blow-up of the loss, as compared toÕ(n k/2 ) for conditional regression. Can we achieve a similar approximation factor for reference class regression? Finally, we still do not know how close to optimal this blow-up of the loss is; in particular, we do not have any lower bounds. Note that this is a computational and not a statistical issue, since we can obtain uniform convergence over all k-DNF conditions.