Paper deep dive
Partial Identification under Causal Orders by Linear Programming
Eric Rossetto, Alessandro Antonucci
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 8/29/2026, 4:38:27 AM
Summary
This paper introduces a method for partial identification of counterfactual queries without requiring a fully specified causal graph. By leveraging the partial topological ordering implied by the counterfactual query itself, the authors reduce the identification task to a linear programming problem. They prove the tightness of the resulting bounds and demonstrate the method's utility through case studies, generalizing the Tian and Pearl (2000) framework.
Entities (10)
Relation Signals (7)
Alessandro Antonucci → affiliatedwith → SUPSI
confidence 95% · Affiliation: Scuola Universitaria Professionale della Svizzera Italiana (SUPSI)
Alessandro Antonucci → affiliatedwith → IDSIA
confidence 95% · Affiliation: Istituto Dalle Molle di Studi sull’Intelligenza Artificiale (IDSIA)
Eric Rossetto → affiliatedwith → SUPSI
confidence 95% · Affiliation: Scuola Universitaria Professionale della Svizzera Italiana (SUPSI)
Eric Rossetto → affiliatedwith → IDSIA
confidence 95% · Affiliation: Istituto Dalle Molle di Studi sull’Intelligenza Artificiale (IDSIA)
Counterfactual Queries → reducedto → Linear Programming
confidence 95% · enables an explicit query parametrisation reducing the identification task to a linear program
Partial Identification under Causal Orders by Linear Programming → generalizes → Tian and Pearl (2000)
confidence 90% · Our work can be viewed as a generalisation of the classical bounding framework of Tian and Pearl (2000)
Causal Graph → notrequiredfor → Partial Identification under Causal Orders by Linear Programming
confidence 90% · challenge this requirement by leveraging structural assumptions that are inherently implied by the query itself
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Non-parametric (partial) identification of counterfactual queries typically relies on a fully specified causal graph. Motivated by settings with incomplete domain knowledge, we challenge this requirement by leveraging structural assumptions that are inherently implied by the query itself. We show that any counterfactual inquiry induces a, mostly partial, topological ordering over relevant variables, which, in turn, enables an explicit query parametrisation reducing the identification task to a linear program. This allows bounding arbitrary counterfactual and nested counterfactual queries. Our work can be viewed as a generalisation of the classical bounding framework of Tian and Pearl (2000), originally developed for probabilities of causation. We also prove the \emph{tightness} of our bounds by constructing structural causal models that attain the bounds whilst being compatible with both the observed data and the query-implied order. To assess both the generality and practical utility of the proposed bounding procedure, we revisit several case studies from the literature, demonstrating how the derived bounds can be used to yield informative insights even in the absence of an input causal graph.
Tags
Links
- Source: https://arxiv.org/abs/2608.24427v1
- Canonical: https://arxiv.org/abs/2608.24427v1
Trouble viewing inline? Open PDF directly →
Full Text
74,805 characters extracted from source content.
Expand or collapse full text
Partial Identification under Causal Orders by Linear Programming Eric Rossetto and Alessandro Antonucci Affiliation: Istituto Dalle Molle di Studi sull’Intelligenza Artificiale (IDSIA) Affiliation: Scuola Universitaria Professionale della Svizzera Italiana (SUPSI) Affiliation: Lugano, Switzerland Affiliation: eric.rossetto, alessandro.antonucci@supsi.ch Abstract Non-parametric (partial) identification of counterfactual queries typically relies on a fully specified causal graph. Motivated by settings with incomplete domain knowledge, we challenge this requirement by leveraging structural assumptions that are inherently implied by the query itself. We show that any counterfactual inquiry induces a, mostly partial, topological ordering over relevant variables, which, in turn, enables an explicit query parametrisation reducing the identification task to a linear program. This allows bounding arbitrary counterfactual and nested counterfactual queries. Our work can be viewed as a generalisation of the classical bounding framework of Tian and Pearl (2000), originally developed for probabilities of causation. We also prove the tightness of our bounds by constructing structural causal models that attain the bounds whilst being compatible with both the observed data and the query-implied order. To assess both the generality and practical utility of the proposed bounding procedure, we revisit several case studies from the literature, demonstrating how the derived bounds can be used to yield informative insights even in the absence of an input causal graph. A Preprint Keywords Structural Causal Models ⋅· Partial Identifiability ⋅· Linear Programming 1 Introduction Causal reasoning deepens our understanding of the effects of policies, actions, and treatments in domains that require accountable decision-making, such as algorithmic fairness, personalised medicine, economics, and the broader social sciences. Quantifying such causal inquiries demands a rigorous mathematical foundation. The work of Pearl (2009) formalised the systems under study using structural causal models (SCMs), which elegantly encode prior causal beliefs via directed graphs (Pearl, 1995). Ideally, these structural assumptions allow causal queries to be uniquely computed—formally identified—from the available set of observations. In many practical settings, however, the available assumptions are insufficient for point identification, necessitating partial identification: computing a range within which the target quantity is guaranteed to lie. Early work by Manski (1990) derived closed-form bounds for simple treatment effect estimation without relying on structural assumptions. Similarly, Tian and Pearl (2000) formulated a linear program to compute sharp bounds on the probabilities of causation (Pearl, 1999) from experimental and observational data under either general or mild assumptions. Whilst weak structural assumptions typically yield substantially wide bounds, assuming a fully specified causal graph might produce tighter intervals. Graph-based reductions to linear programs were first developed for the instrumental variable setting (Balke and Pearl, 1994a; Balke and Pearl, 1994b) and later extended to generalised instrumental constraints (Sachs et al., 2023). For general graphical assumptions, partial identification is reduced to a polynomial program via canonicalisation, though the exponential growth of the resulting exogenous state space ultimately necessitates approximate methods (Zhang et al., 2022; Duarte et al., 2024). We revisit partial identification in the spirit of Tian and Pearl (2000), whose bounding method did not require a specification of a causal graph; instead, their proposal relied on a direct parametrisation of the joint counterfactual distribution. To extend this framework to arbitrary queries, we leverage the minimal structural assumptions inherently implied by the query itself, reducing the identification task to a linear program without imposing additional structural commitments. Our specific contributions are as follows: (i) we reduce partial identification under a total order of variables to a linear program, capable of handling arbitrary counterfactual, including nested, queries; (i) we prove that the resulting bounds are tight by constructing witnessing structural causal models that attain them; (i) we prove that the order-dependent optimisation task is invariant under both the marginalisation of unqueried variables and the specific choice of total order (compatible with the query); and (iv) we show that the linear program admits a lifted reformulation by exploiting symmetries induced by the observations and the query’s logical evaluations. The presented results define an automated framework that enables the bounding of arbitrary queries without the need for explicit graphical assumptions. It is important to emphasise that graph-based methods generally yield narrower bounds by exploiting conditional independencies; in contrast, we pose ourselves in a regime of incomplete domain knowledge, abstaining from such graphical commitments. After introducing the necessary background in Sect. 2, Sect. 3 reduces partial identification to a linear programme, which Sect. 4 extends to accommodate conditional queries, nested counterfactuals, and structural constraints. Sect. 5 demonstrates practical case studies, followed by conclusions in Sect. 6; proofs, additional experiments and discussions are in Apps. A, B and C, respectively. 2 Background Notation. Let us first review some necessary background on causality. We denote variables by capital letters, whilst small and calligraphic letters are used instead for the states and the sample spaces. Thus, x∈x is a state of X. We consider discrete variables only. Bold is used for sets of variables. The indicator function that is one when its argument is true and zero otherwise is denoted as ⟦⋅⟧ · , whilst |⋅||·| is the cardinality of the set in the argument. Structural Causal Models (SCMs). Following Pearl (2009), an SCM ℳ≔⟨,,ℱ,M V, U,F, P()⟩P( U) is such that V is its set of endogenous variables, U the set of exogenous variables distributed according to P()P( U), and ℱF a set fVV∈\f_V\_V∈ V of structural equations (SEs). The SE of each V∈V∈ V determines v=fV(paV,V)v=f_V(pa_V, u_V), where paV∈PaV⊆∖Vpa_V _V V \V\ and V⊆ U_V U are, respectively, the endogenous and the exogenous inputs of the SE. We may view an SE as a collection of deterministic mechanisms, or mappings, from input configurations to output states. If P()P( U) is unavailable, M≔⟨,,ℱ⟩M V, U,F is termed partially specified SCM (PSCM). Both SCM ℳM and PSCM M induce a directed graph G whose nodes represent variables and an edge (X,Y)(X,Y) exists iff X∈PaYX _Y. We restrict our analysis to recursive SCMs, i.e., such that G is acyclic. We adopt a kinship terminology (e.g., parents) to present graphical relations between variables. The notation MM_M is finally used for the set of all the SCMs that share the SEs’ signature, i.e., parent and child variables, with M. Interventions. In an SCM ℳM, P()P( U) and ℱF induce a joint observational distribution Pℳ()P_M( V). An SCM natural regime can be modified by acting on its SEs. For an arbitrary set ⊆ X V, an intervention yields a new model ℳM_ x (i.e., a sub-model) by replacing the SEs of X with constant assignments ← X← x, whilst keeping P()P( U) and all other SEs unchanged. The resulting, interventional, distribution is denoted as Pℳ(∖)P_M_ x( V X). For any Y∈Y∈ V, the potential response Y()Y_ x( u) denotes the value of Y in ℳM_ x given = U= u. A counterfactual variable (or potential outcome) Y_ x is the variable induced by Y()Y_ x( u) when ∼P() u P( U). Counterfactuals. Following the logical formalisation of a PSCM by Halpern (2000), we define an atomic proposition as the assignment Y=yY_ x=y (for short, y_ x), where Y∈,⊆Y∈ V, X V and y∈y . A counterfactual event γ is a conjunction of k atomic propositions, γ≡y(1)(1)∧⋯∧y(k)(k)γ≡ y^(1)_ x^(1) … y^(k)_ x^(k). A counterfactual query Pℳ(γ)P_M(γ) asks for the probability of γ with respect to an SCM ℳM. This requires the simultaneous computation of potential responses of multiple sub-models ℳ(i)i=1k\M_ x^(i)\_i=1^k, coupled by the shared P()P( U), i.e., Pℳ(γ)=∑∈⟦⋀i=1kY(i)(i)()=y(i)⟧⋅P().P_M(γ)= _ u∈ U _i=1^kY^(i)_ x^(i)( u)=y^(i) · P( u)\,. (1) The evaluation accumulates over the exogenous instantiations u satisfying the formula γ. We denote such an entailment as (ℳ,)⊧γ(M, u) γ. When it is clear from the context which SCM we refer to, we omit specifying it and rewrite the right-hand side of Eq. (1) as ∑⟦⊧γ⟧P() _ u u γ P( u). If all the subscripts agree on x, then Eq. (1) can be retrieved from the interventional distribution PℳP_M_ x. In general, this is not the case and we eventually characterise a counterfactual distribution. Finally, we denote as γ⊆ V_γ V the endogenous, no matter whether queried or intervened, variables in γ. Nested Counterfactuals. We can further generalise the intervention for the counterfactual variable YxY_x to allow representing settings where the variable X is set to behave as another variable, say XzX_z. A random variable Y in such system is represented with a potential outcome of the form YXz=yY_X_z=y, which is called a nested counterfactual. This expresses the statement “Y becomes y when we set X to whatever value it would have taken if we had set Z to z”. Queries involving nested counterfactuals can be rewritten into a marginalisation over conjunctions of standard, atomic counterfactual events (Correa et al., 2021), i.e., (YXz=y)≡⋁x∈(Yx=y∧Xz=x). (Y_X_z=y )≡ _x (Y_x=y X_z=x )\,. (2) (Partial) Identifiability. Eq. (1) allows to compute counterfactual queries in an SCM ℳM. Yet, the exogenous distribution P()P( U) is rarely available, and one must cope with the corresponding PSCM M and a dataset D of endogenous observations, used to estimate an empirical distribution P^() P( V). In general, there are multiple SCMs in MM_M, that could have generated P^() P( V) and yield different values to a query over γ. In such a partially identifiable setting, we can only compute bounds to the query, i.e., the minimum: minℳ∈MPℳ(γ)s.t.Pℳ()=P^(), M _M \,\,P_M(γ) .t. P_M( V)= P( V)\,, (3) and the maximum, which is analogously defined. In the remainder of the paper, just for the sake of brevity, we focus on the minimisation problem, noticing that the maximisation counterpart can be handled analogously. Those optimisation tasks have been shown to be NP-hard (Zaffalon et al., 2024). Nevertheless, approximate techniques were proposed to estimate bounds efficiently (Zhang et al., 2022; Zaffalon et al., 2024; Duarte et al., 2024; Bjøru et al., 2025). Standard (point) identifiability emerges as a special case when the two bounds coincide. Queries and Partial Orders. Although potential outcome semantics allow us to consider a variable YxY_x for any Y,X∈Y,X∈ V and x∈x , an intervention on X can alter the realisation of Y only if there exists a directed path from X to Y in the causal graph induced by the assumed SCM ℳM over V. This is formalised as follows. Proposition 1 (Tian, 2002). If Y,X∈Y,X∈ V are distinct endogenous variables of an SCM ℳM such that X∉Anc(Y)X (Y), then, ∀ u, Yx()=Y()Y_x( u)=Y( u), and thus Pℳ(Yx)=Pℳ(Y)P_M(Y_x)=P_M(Y). Thus, a counterfactual variable YxY_x can be non-trivial only if X precedes Y in the underlying causal graph. In this paper, we restrict attention to non-trivial queries: consequently, when the graph is not available, a query over γ can be also used to induce a collection of precedence constraints, corresponding to a partial order ≺γ _γ on γ V_γ. A “cyclic” event such as yx∧xyy_x x_y would be ill-posed as it induces incompatible precedence constraints. Such cases would naturally call for non-recursive models (see, for instance, Cozman et al., 2025 and Bongers et al., 2021), which fall outside the scope of this paper. 3 Bounding Framework The main goal of this paper is to relax the optimisation task in Eq. (3) to cope with situations where the underlying PSCM M is unavailable. We start by assuming that the endogenous variables are sorted by a total order σ, corresponding to ≔(V(1),…,V(k)) V (V^(1),…,V^(k)). In this setting, we can replace in Eq. (3) the set MM_M of SCMs compatible with the PSCM M with the set σM_σ of all SCMs whose causal graph is consistent with σ, i.e., minℳ∈σPℳ(γ)s.t.Pℳ()=P^(). M _σ P_M(γ) .t. P_M( V)= P( V)\,. (4) For the moment, we consider γ such that γ= V_γ= V, i.e., all endogenous variables appear in the query. A more general account of the case γ⊂ V_γ⊂ V is addressed later. Response Representation. Unlike Eq. (3), the optimisation task in Eq. (4) refers to a collection of SCMs possibly based on different PSCMs. Yet, given only σ and V, we consider a complete canonical PSCM over V denoted as MσM_σ and such that: (i) PaV(i)Pa_V^(i) coincides with all the predecessors of V(i)V^(i) according to σ and they are denoted as ≺i V_ i; (i) there is a single exogenous variable U, which is a parent of all the endogenous variables; (i) the SEs of each V∈V∈ V and the cardinality of U are such that MσM_σ is canonical in the sense of Zhang et al. (2022), this basically meaning that the exogenous states of U are indexing all possible structural relations between the endogenous variables and their endogenous parents (see App. C for a more detailed account). The optimisation task in Eq. (4) can be equivalently achieved by considering the complete canonical PSCM MσM_σ only, as stated by the following result. Theorem 1. Eq. (3) with M=MσM=M_σ gives the same minimum of Eq. (4), i.e., minℳ∈σPℳ(γ)=minℳ∈MσPℳ(γ), _M _σP_M(γ)= _M _M_σP_M(γ)\,, (5) where Pℳ()=P^()P_M( V)= P( V) is required on both sides. The above result reduces the optimisation in Eq. (4)—which considers a collection of structurally heterogenous SCMs—to span only those SCMs sharing the same causal graph as MσM_σ. This reduction provides a parametrisation of the task, where the optimisation variables are associated with the exogenous probabilities as in Eq. (1). To avoid explicitly modelling these latent variables, we build upon the approach proposed by Balke and Pearl (1994a); Balke and Pearl (1994b). Here, we consider all the potential outcomes that MσM_σ is allowed to generate as to enumerate all possible endogenous mappings from parents to child. This is achieved by considering a response signature σ S_σ for MσM_σ, corresponding to the following collection of potential outcomes: σ≔V≺i(i)≺i∈≺i=1,…,k. S_σ \V^(i)_ v_ i \_ v_ i∈ V_ i^i=1,…,k\,. (6) The following result shows that the probability of any counterfactual query can be expressed not only as a linear combination of joint probabilities over the exogenous variables, as in Eq. (1), but also as a linear combination of joint probabilities over σ S_σ. Theorem 2. A query Pℳ(γ)P_M(γ) in an SCM ℳ∈σM _σ can be written as a linear combination of the probabilities of the response signature, i.e., Pℳ(γ)=∑σ∈σ⟦σ⊧γ⟧⋅P(σ).P_M(γ)= _ s_σ∈ S_σ s_σ γ · P( s_σ)\,. (7) By viewing observational queries as a special case of counterfactual queries in which interventions are performed on the empty set, we also obtain the following additional result. Corollary 1. Any observational probability Pℳ()P_M( v) in an SCM ℳ∈σM _σ can be written as a linear combination of the probabilities of the response signature, i.e., Pℳ()=∑σ∈σ⟦σ⊧⟧⋅P(σ).P_M( v)= _ s_σ∈ S_σ s_σ v · P( s_σ)\,. (8) Linear Programming (LP) Formulation. Thm. 2 and Cor. 1 allow us to address the optimisation task in Eq. (4) via the signature joint states. The resulting linearity of the constraints and the objective function allows us to rewrite the bounding task as a LP task: min∑σ∈σ⟦σ⊧γ⟧⋅qσs.t.∑σ∈σ⟦σ⊧⟧⋅qσ=P^(),∀∈;∑σ∈σqσ=1;qσ≥0,∀σ∈σ. split & _ s_σ∈ S_σ s_σ γ · q_ s_σ\\ s.t. & _ s_σ∈ S_σ s_σ v · q_ s_σ= P( v)\,,∀ v∈ V\;; _ s_σ∈ S_σq_ s_σ=1\;; q_ s_σ≥ 0\,,\,∀ s_σ∈ S_σ\,. split (9) with semantics qσ=P(σ)q_ s_σ=P( s_σ) for the optimisation variables. The following result establishes that the above LP task solves the relaxed optimisation over SCMs consistent with σ by yielding tight identification bounds. Theorem 3. The LP task in Eq. (9) gives the solution of Eq. (4). Furthermore, there exists an SCM ℳ∈σM _σ such that Pℳ(γ)P_M(γ) coincides with the solution of the LP. Signature State Partitioning. The latent space reduction to finite response-function classes (Balke and Pearl, 1994a; Duarte et al., 2024) enables us to neatly solve Eq. (4) through Eq. (9). However, a further reduction can be achieved by observing that the LP in Eq. (9) depends exclusively on the Boolean evaluations of the signatures against the observations v and the counterfactual event γ. Signatures states that evaluate identically against these observations and events can be grouped into equivalence classes σ∈σ≔σ/∼γ t_σ∈ T_σ S_σ/ _γ, where: σ∼γσ′⇔((σ⊧γ)≡(σ′⊧γ))∧(∃∈s.t.(σ⊧)∧(σ′⊧)). s_σ _γ s _σ (( s_σ γ)≡( s _σ γ) ) (∃ v∈ V\,\,s.t.\,\,( s_σ v) ( s _σ v) )\,. (10) Each class σ∈σ t_σ∈ T_σ inherits the binary evaluations from the signatures it includes, hence ⟦σ⊧⟧ t_σ v and ⟦σ⊧γ⟧ t_σ γ . By defining aggregated masses as variables, qσ≔∑σ∈σqσq_ t_σ _ s_σ∈ t_σq_ s_σ, and applying Thm. 2 and Cor. 1, the LP in Eq. (9) reduces to a lifted LP task (Kersting et al., 2017): min∑σ∈σ⟦σ⊧γ⟧⋅qσs.t.∑σ∈σ⟦σ⊧⟧⋅qσ=P^(),∀∈;∑σ∈σqσ=1;qσ≥0,∀σ∈σ. split & _ t_σ∈ T_σ t_σ γ · q_ t_σ\\ s.t. & _ t_σ∈ T_σ t_σ v · q_ t_σ= P( v)\,,∀ v∈ V\;; _ t_σ∈ T_σq_ t_σ=1\;; q_ t_σ≥ 0\,,∀ t_σ∈ T_σ\,. split (11) The equivalence mapping improves the tractability of the LP in Eq. (9) whilst maintaining both soundness and tightness, as it follows from the following result. Theorem 4. The lifted LP task in Eq. (11) gives the solution of Eq. (4). Furthermore, there exists an SCM ℳ∈σM _σ such that Pℳ(γ)P_M(γ) coincides with the solution of the LP. Marginalising Non-Queried Variables. So far, we assumed that all the endogenous variables to appear in the query, i.e., γ= V_γ= V. We now extend our results to the more general case in which γ⊆ V_γ V. To this end, let σγM_σ V_γ denote the set of SCMs defined only over γ V_γ and compatible with σ. Observe that σ=σM_σ V=M_σ and, in general, σγ⊆σM_σ V_γ _σ, where, for notational ease, we assumed that the endogenous variables might not appear explicitly in an SCM, thus being uniformly distributed. The following result holds. Theorem 5. The solution of Eq. (4) coincides with: minℳ∈σγPℳ(γ)s.t.Pℳ(γ)=P^(γ). M _σ V_γ P_M(γ) .t. P_M( V_γ)= P( V_γ)\,. (12) where P^(γ) P( V_γ) is obtained from P^() P( V) by marginalisation. As a consequence of the above result, when γ⊂ V_γ⊂ V, we marginalise the empirical distribution before applying the LP formulation to solve Eq. (12) and hence Eq. (4). Coping with Partial Orders. The optimisation task in Eq. (4) assumes a total order σ over the endogenous variables. However, as discussed in Sect. 2, a counterfactual query over γ might induce only a partial order ≺γ _γ. Nevertheless, a general formulation is readily obtained by considering the bound over all linear extensions Σ(≺γ) ( _γ): minσ∈Σ(≺γ)minℳ∈σγPℳ(γ)s.t.Pℳ(γ)=P^(γ). σ∈ ( _γ) \, M V_γ_σ P_M(γ) .t. P_M( V_γ)= P( V_γ)\,. (13) When γ contains no interventions, ≺γ _γ imposes no ordering constraints, so that all possible total orders must be considered, and we have |Σ(≺γ)|=|γ|!| ( _γ)|=| V_γ|!. Interventions restrict this search space by ruling out inconsistent orderings; nonetheless, the number of linear extensions can still scale exponentially with respect to |γ|| V_γ|. In spite of this, the following result establishes that evaluating an arbitrary linear extension is sufficient, thus eliminating the need for a brute-force optimisation over all the orders in Σ(≺γ) ( _γ). Theorem 6. Consider a counterfactual query over γ inducing a partial order ≺γ _γ. Any total order σ∈Σ(≺γ)σ∈ ( _γ) gives the same minimum in Eq. (12) and, hence, the solution of Eq. (13). As a consequence of Thm. 6, we are not required to assume or enforce any specific total order. The choice of a particular σ serves only as a possible parametrisation of the task. Complexity. Thm. 3 enables us to solve the optimisation problem in Eq. (4) using a linear solver that implements the LP defined in Eq. (9). It is straightforward to verify that, if the number of variables in the query is k≔|γ|k | V_γ| (see Thm. 5) and n≔maxi=1k|(i)|n _i=1^k|V^(i)| denotes the maximum endogenous cardinality, the number of potential outcomes comprising the signature σ S_σ in Eq. (6) is bounded by ∑i=0k−1ni=nk−1n−1 _i=0^k-1n^i= n^k-1n-1. The number of joint signature states |σ|| S_σ|, which dictates the number of optimisation variables in the naïve LP of Eq. (9), is therefore double exponential, namely (nnk−1)O(n^n^k-1). In contrast, by Thm. 4, the size of the quotient space |σ|| T_σ|, which governs the size of the lifted LP in Eq. (11), is bounded by (∏i=1k|(i)|)=(nk)O ( _i=1^k|V^(i)| )=O(n^k). Thus, to partially identify any query given the empirical marginal P^(γ) P( V_γ), we only need to solve Eq. (11), whose optimisation space scales linearly with the size of the sample space of γ V_γ. 4 Extensions In the current section, we present some important extensions of the LP mapping introduced in the previous section. Conditional Queries. To bound a conditional query P(γ∣δ)≔P(γ∧δ)/P(δ)P(γ δ) P(γ δ)/P(δ), it suffices to note that both terms decompose linearly over σ S_σ via Thm. 2 and Cor. 1. This reformulates the bounding task as a linear-fractional programme, which is reduced to a standard LP via the Charnes-Cooper transformation (Charnes and Cooper, 1962). We incorporate evidence by refining the equivalence relation in Eq. (10) to evaluate both γ∧δγ δ and δ: σ∼γ|δσ′⇔σ∼γ∧δσ′∧((σ⊧δ)≡(σ′⊧δ)). s_σ _γ δ s _σ s_σ _γ δ s _σ (( s_σ δ)≡( s _σ δ) )\,. (14) Letting bψ≔⟦σ⊧ψ⟧b_ψ s_σ ψ , the logical constraint γ∧δ⇒δγ δ δ enforces bγ∧δ≤bδ∈0,1b_γ δ≤ b_δ∈\0,1\. This restricts evaluations to three valid Boolean states: (bγ∧δ,bδ)∈(0,0),(0,1),(1,1)(b_γ δ,b_δ)∈\(0,0),(0,1),(1,1)\. The optimisation space retains the same asymptotic size discussed in Sect. 3. Nested Queries. The LP mapping in Eq. (9), and by extension the lifted LP in Eq. (11), accommodates nested counterfactuals. Although nested variables do not appear directly in the response signature σ S_σ, their evaluation remains a linear optimisation task. For a single nested counterfactual, it is sufficient to consider the probabilistic version of Eq. (2), i.e., Pℳ(YXz=y)=∑x∈Pℳ(Yx=y,Xz=x).P_M(Y_X_z=y)= _x P_M(Y_x=y,X_z=x)\,. (15) In the general case, it is sufficient to perform such un-nesting on all the nested terms to obtain a linear function of standard counterfactual queries and thus a linear objective function. Experimental Data. Beyond observational data, it is straightforward to incorporate experimental distributions obtained directly from randomised controlled trials. By Thm. 2, experimental probabilities translate directly into linear constraints on the signature probabilities in Eq. (9). To preserve the validity of the reduction in Eq. (11), the equivalence relation is refined by further requiring that equivalent signatures agree on the Boolean evaluation of the relevant interventional events. For instance, given an experimental probability P(yx)P(y_x), the relation in Eq. (10) is augmented with the logical condition (σ⊧yx)≡(σ′⊧yx)( s_σ y_x)≡( s _σ y_x). Tracking these additional constraints partitions the signature space more finely, thereby scaling the quotient space by a factor determined by the number of interventional events considered. Furthermore, when multiple studies or heterogeneous data sources are available, one can seamlessly incorporate all evidence by appending the constraints induced by each trial in the same LP. Such a procedure natively performs automated data fusion, analogously to the approach developed by Zaffalon et al. (2023) for PSCMs. Monotonicity and Weak Exogeneity. For ordinal variables, monotonicity (Angrist et al., 1996; Manski, 1997) posits a directional effect (e.g., Yx1≥Yx0Y_x_1≥ Y_x_0 if x1>x0x_1>x_0). Weak exogeneity (Tian and Pearl, 2000) requires that a variable is unaffected by unobserved confounders. Both assumptions translate expert knowledge into logical restrictions on the response signature σ S_σ: monotonicity excludes states violating the inequality by fixing relevant optimisation variables to zero; exogeneity requires that causal effects are identified from data, i.e., P(yx)=P(y∣x)P(y_x)=P(y x), functionally mirroring the embedding of experimental data. XXYYUU (a) ZZXXYYUU (b) XXWWYYUU (c) Figure 1: Causal diagrams considered in the examples. (a) Confounded treatment-outcome pair. (b) Instrumental Variable (IV) setting, where Z serves as an instrument for X. (c) Mediation model with a direct effect and confounding between the mediator W and the outcome. Grey nodes denote unobserved exogenous variables. 5 Examples Let us validate the application of our LP-based bounding scheme across a range of scenarios. We compare our bounds against analytical formulae and results from the literature.11 1 A Python implementation of our method is available at w.github.com/Erhtric/copid. The experimental setup employed the COIN-OR CBC open-source solver via the PuLP modeller (w.coin-or.org). Probabilities of Causation. As a first illustrative application of our method, we derive bounds for the well-known probability of necessity and sufficiency (PNS, Pearl, 1999). This query refers to the counterfactual event γ≔(YX=0=0)∧(YX=1=1)γ (Y_X=0=0) (Y_X=1=1), which is intended to capture the causal relationship between two Boolean endogenous variables, X and Y. The only total order σ compatible with ≺γ _γ is σ:=(X,Y)σ:=(X,Y). The associated response signature is therefore σ=X,YX=0,YX=1 S_σ=\X,Y_X=0,Y_X=1\. Accordingly, PNS bounds are obtained by solving LPs with eight variables, corresponding to the joint states of σ S_σ; and four linear constraints induced by the observational distribution P^(X,Y) P(X,Y). The resulting solutions coincide with the classical closed-form bounds derived by Tian and Pearl (2000). The equivalence relation in Eq. (10) prunes the parameter space by identifying empty equivalence classes. Specifically, the candidate classes defined by the Boolean tuples ⟨(X=x,Y=y),⟦γ⟧⟩ (X=x,Y=y), γ for (x,y)∈(1,0),(0,1)(x,y)∈\(1,0),(0,1)\ contain no valid signatures. This is because an observational state of (X=1,Y=0)(X=1,Y=0) forces YX=1=0Y_X=1=0 via consistency (Galles and Pearl, 1998), directly contradicting the requirement in γ that YX=1=1Y_X=1=1. We similarly proceed for (X=0,Y=1)(X=0,Y=1). The tightness (referred to as sharpness in the original work) of these bounds follows from Thm. 3: one can construct two complete canonical SCMs under which the query attains the lower and upper bounds (see also the constructive argument in the proof). Finally, by considering the extension to conditional queries discussed in Sect. 4, analogous results can be readily obtained for the other probabilities of causation (PN and PS, Pearl, 1999). Moreover, unlike the analytical results of Tian and Pearl (2000), which are restricted to binary variables, our approach naturally generalises to the non-binary case. By simply adjusting the cardinality of the variables in the response signature σ S_σ, we can derive tight bounds for probabilities of causation involving multi-valued treatments and outcomes (Sun et al., 2025) without requiring new algebraic derivations. From this perspective, our work is a genuine generalisation of the analysis in Tian and Pearl (2000), with a clear semantics, and allowing to bound arbitrary counterfactual queries, whether nested or not. Average Causal Effect (ACE). The ACE measures the average change due to treatment: ACEx,x′(y)≔P(Yx=y)−P(Yx′=y).ACE_x,x (y) P(Y_x=y)-P(Y_x =y)\,. (16) For an outcome Y∈y,y′Y∈\y,y \ and treatment ||=3|X|=3, the unique order σ=(X,Y)σ=(X,Y) induces an LP with 24 signature variables and six empirical constraints in Eq.(9). The resulting bounds, [P(x,y)+P(x′,y′)−1,1−P(x,y′)−P(x′,y)][P(x,y)+P(x ,y )-1,1-P(x,y )-P(x ,y)] (Balke and Pearl, 1997), recover those of Sachs et al. (2023, Sect. 6.1) whose causal assumptions are in Fig. 1(a). Crucially, our method achieves this without explicitly parametrising a fully specified causal graph nor the SEs. A different account of the ACE bounding can be found in App. B. Controlled Direct Effect (CDE). Assumption Set CDEx0,x1w0(y1)CDE^w_0_x_0,x_1(y_1) CDEx0,x1w1(y1)CDE^w_1_x_0,x_1(y_1) |σ|| S_σ| |σ|| T_σ| Lower Upper Lower Upper No additional assumptions Cai et al. (2008) (Fig. 1(c)) -0.200 0.385 -0.781 0.634 – – This paper -0.599 0.694 -0.891 0.817 128 20 Monotonicity Cai et al. (2008) 0.000 0.036 0.000 0.582 – – This paper 0.000 0.519 0.000 0.791 36 13 Weak exogeneity This paper -0.200 0.385 -0.782 0.633 128 64 Weak exogeneity + monotonicity This paper 0.000 0.035 0.000 0.580 36 24 Table 1: Comparison of CDE bounds on the (binary) Lipid data (Cai et al., 2008, Tab. 1) against analytical bounds under varying assumptions. In this experiment, Monotonicity forces positive effects of X on W, X on Y, and W on Y; whilst, weak exogeneity imposes no-interaction between X and W on Y. Consider a setting in which the effect of a treatment X on an outcome Y is mediated by a variable W. In this context, the ACE in Eq. (16) can be refined into a CDE to measure the change due to treatment at different mediator strata: CDEx0,x1w(y)≔P(Yx1,w=y)−P(Yx0,w=y).CDE^w_x_0,x_1(y) P(Y_x_1,w=y)-P(Y_x_0,w=y)\,. (17) The query induces a partial order, namely X,W≺Y\X,W\ Y. By Thm. 6, we can choose an arbitrary order among its set of linear extensions, together with an empirical distribution P^(X,W,Y) P(X,W,Y), to solve Eq. (11). We validate the method against the problem proposed by Cai et al. (2008) of bounding Eq. (17) under the assumptions encoded in Fig. 1(c). In such a setting, we first consider all binary variables and a randomised treatment (i.e., no confounders between X and Y), examining data from the Coronary Primary Prevention Trial (Cai et al., 2008). A summary of the computed bounds is reported in Tab. 1. One can observe that when no assumptions are made, or when solely monotonicity is assumed, the minimal structure inherently implied by the query yields valid, but wide, intervals. Notably, the tight bounds reported by Cai et al. (2008) rely on a fully specified causal graph. To mirror their structural assumptions within our framework, one must incorporate a weak exogeneity constraint—also referred to as the no-interaction assumption—into the signature space σ S_σ. Specifically, enforcing the independence Y,W⟂X\Y,W\ \!\!\! X recovers Cai et al. (2008) analytical bounds. As a special instance, one can also require Y⟂X,WY \!\!\! \X,W\ to collapse the bound to a point estimate, since P(yx,w)=P(y∣x,w)P(y_x,w)=P(y x,w). As expected, introducing these independence constraints expands the number of optimisation parameters |σ|| T_σ|, whilst introducing monotonicity constraints shrinks the space (cf. Tab. 1). A similar discussion holds for categorical variables, where the LP reduction can be systematically constructed and compared against the equations presented in Cai et al. (2008). Natural Direct Effect (NDE). Lower Upper |σ|| S_σ| |σ|| T_σ| Assumptions for NDEx0,x1(y1)NDE_x_0,x_1(y_1) Query-induced Order -0.691 0.688 294,912 60 Y⟂XY \!\!\! X -0.503 0.497 294,912 96 Y,W⟂X\Y,W\ \!\!\! X -0.503 0.497 294,912 528 Y⟂X,WY \!\!\! \X,W\ -0.628 0.679 294,912 89,172 Y⟂X,WY \!\!\! \X,W\ & W⟂XW \!\!\! X -0.489 0.637 294,912 294,912 Query P(Yx1,W=w0=y1,Wx0=w0)P(Y_x_1,W=w_0=y_1,W_x_0=w_0) 0.000 0.610 294,912 37 P(Yx1,W=w1=y1,Wx0=w1)P(Y_x_1,W=w_1=y_1,W_x_0=w_1) 0.000 0.493 294,912 37 P(Yx1,W=w2=y1,Wx0=w2)P(Y_x_1,W=w_2=y_1,W_x_0=w_2) 0.000 0.365 294,912 37 P(Yx1,W=w3=y1,Wx0=w3)P(Y_x_1,W=w_3=y_1,W_x_0=w_3) 0.000 0.415 294,912 37 P(Yx1,W=w4=y1,Wx0=w4)P(Y_x_1,W=w_4=y_1,W_x_0=w_4) 0.000 0.357 294,912 37 P(Yx1,W=w5=y1,Wx0=w5)P(Y_x_1,W=w_5=y_1,W_x_0=w_5) 0.000 0.391 294,912 37 P(Yx0=y1)P(Y_x_0=y_1) 0.312 0.691 8 6 Table 2: UC Berkeley dataset: NDE bounds under various weak exogeneity assumptions (top) and individual unnested query terms under query-induced assumptions (bottom). The CDE in Eq. (17) refers to a fixed state of the mediator. By contrast, the NDE evaluates the effect of a treatment when the mediator is set to the value it would have attained in the absence of treatment. This leads to the nested counterfactual query: NDEx0,x1(y1)≔P(Yx1,Wx0=y1)−P(Yx0=y1).NDE_x_0,x_1(y_1) P(Y_x_1,W_x_0=y_1)-P(Y_x_0=y_1)\,. (18) As a matter of fact, the nested syntactical structure of the query induces a unique complete order X≺W≺YX W Y. We proceed to apply the un-nesting procedure in Eq. (15), which, in turn, allows one to rewrite nested counterfactuals like in Eq. (18) as sums of standard (non-nested) counterfactual queries. This transformation preserves the linear structure of the objective function in our LP formulation and one can consider the PSCM whose endogenous structure is like in Fig. 1(c) but with one exogenous U parent of all variables. As discussed previously, by Thm. 4, the optimisation space’s size is bounded by O(||⋅||⋅||)O(|X|·|W|·|Y|) (cf. Tab. 2), making our approach computationally feasible for modern LP solvers. To the best of our knowledge, no general bounding techniques are currently available for queries of this kind without imposing additional structural assumptions about the underlying causal graph. Berkeley Admissions Dataset. For a numerical illustration, we consider the popular UC Berkeley admissions dataset (Bickel et al., 1975). In this dataset, the binary variables X and Y represent gender (with x1x_1 indicating female) and admission outcome (with y1y_1 denoting acceptance). The mediator W is a categorical variable with six states corresponding to the department to which the applicant applied. We omit monotonicity assumptions, as deterministic effects are tricky to justify in social contexts. Tab. 2 summarises the total NDE bounds under varying weak exogeneity assumptions, alongside the decomposed counterfactual terms from Eq. (15). In the context of university admissions, enforcing Y⟂XY \!\!\! X assumes no confounding jointly affecting gender and admission. The stricter condition Y,W⟂X\Y,W\ \!\!\! X treats gender as if it were completely randomised, assuming department choice does not depend on gender. Conversely, Y⟂X,WY \!\!\! \X,W\ dictates that the admission committee’s decision is entirely independent of unobserved background factors influencing either the applicant’s gender or department preference. It is worth noting that when all possible exogeneity assumptions are employed, the equivalence reduction from Sect. 2 fails to compress the optimisation space, as enforcing mutual independence across all response mechanisms requires tracking each signature state individually. Whilst these assumptions narrow the bounds, at the cost of expanding the optimisation space |σ|| T_σ|, they risk misrepresenting complex dynamics. Unlike bespoke analytical methods that require graphical constraints to bound nested counterfactuals, our LP formulation automatically evaluates the NDE using only the minimal assumptions inherent to the query itself. For reference, under the standard fairness model (Plečko and Bareinboim, 2024), the NDE evaluates to a point estimate of 0.043. Concluding a lack of discriminative policies based on this point estimate might be completely inaccurate, as the stark contrast with our wider, assumption-free intervals demonstrates how such conclusions rely heavily on rigid structural assumptions. Furthermore, it is interesting to notice that directly bounding the full NDE yields substantially tighter intervals than combining separately bounded decomposed components via interval arithmetic. 6 Conclusions We presented a linear programming framework for partially identifying arbitrary, including nested, counterfactual queries, relying exclusively on structural assumptions inherent to the query itself. In doing so, we successfully extended the classical bounding results of Tian and Pearl (2000) for probabilities of causation. We proved that partial identification under these minimal assumptions yields tight intervals, and established their invariance with respect to the marginalisation of the variables not appearing in the query. Rather than requiring an input causal graph, an SCM attaining the bounds emerges constructively as a by-product. Furthermore, exploiting query symmetries drastically reduces the optimisation space, rendering the LP tractable. We acknowledge an inherent trade-off in the proposed method: reducing the reliance on prior domain knowledge, and explicit conditional independencies, produces intervals that are naturally wider than those derived from a fully specified causal graph. Nevertheless, as demonstrated in our case studies, this framework serves as a general and computationally efficient method to bound arbitrary inquiries under minimal assumptions. Future work could investigate extending this framework to cyclic causal structures, examining the impact of different valid causal orders, and investigate finer quotient-space reductions to further enhance scalability. Acknowledgements The research of Eric Rossetto was funded by Swiss Post AG as part of a Doctoral Studies Grant on AI and Robustness. The authors thank the anonymous reviewers for their constructive comments, and in particular the first reviewer for insightful remarks on the underlying assumptions of our approach and its possible extension to cyclic models. References Angrist et al. (1996) Joshua D. Angrist, Guido W. Imbens, and Donald B. Rubin. Identification of causal effects using instrumental variables. Journal of the American Statistical Association, 91(434):444–455, 1996. Balke and Pearl (1994a) Alexander Balke and Judea Pearl. Counterfactual probabilities: Computational methods, bounds and applications. In Proceedings of the Tenth International Conference on Uncertainty in Artificial Intelligence, pages 46–54, 1994a. Balke and Pearl (1994b) Alexander Balke and Judea Pearl. Probabilistic evaluation of counterfactual queries. In Probabilistic and Causal Inference, volume Proceedings of the Twelfth National Conference on Artificial Intelligence of AAAI ’94, pages 230–237. American Association for Artificial Intelligence, 1994b. Balke and Pearl (1997) Alexander Balke and Judea Pearl. Bounds on treatment effects from studies with imperfect compliance. Journal of the American Statistical Association, 92(439):1171–1176, 1997. Bickel et al. (1975) Peter Bickel, Eugene Hammel, and William O’Connell. Sex bias in graduate admissions: data from Berkeley. Science, 187(4175):398–404, 1975. Bjøru et al. (2025) Anna Rodum Bjøru, Rafael Cabañas, Helge Langseth, and Antonio Salmerón. Divide and conquer for causal computation. International Journal of Approximate Reasoning, 186:109520, 2025. Bongers et al. (2021) Stephan Bongers, Patrick Forré, Jonas Peters, and Joris M Mooij. Foundations of structural causal models with cycles and latent variables. The Annals of Statistics, 49(5):2885–2915, 2021. Cai et al. (2008) Zhihong Cai, Manabu Kuroki, Judea Pearl, and Jin Tian. Bounds on direct effects in the presence of confounded intermediate variables. Biometrics, 64(3):695–701, 2008. Charnes and Cooper (1962) Abraham Charnes and William W. Cooper. Programming with linear fractional functionals. Naval Research Logistics Quarterly, 9(3-4):181–186, 1962. Correa et al. (2021) Juan D. Correa, Sanghack Lee, and Elias Bareinboim. Nested counterfactual identification from arbitrary surrogate experiments. In Proceedings of the 35th International Conference on Neural Information Processing Systems. Curran Associates Inc., 2021. Cozman et al. (2025) Fabio G. Cozman, Radu Marinescu, Junkyu Lee, Alexander Gray, and Denis D. Mauá. Dealing with cycles in graph-based probabilistic models: the case of logical credal networks. In Proceedings of the Fourteenth International Symposium on Imprecise Probabilities: Theories and Applications, volume 290 of Proceedings of Machine Learning Research, pages 93–102. PMLR, 15–18 Jul 2025. Duarte et al. (2024) Guilherme Duarte, Noam Finkelstein, Dean Knox, Jonathan Mummolo, and Ilya Shpitser. An automated approach to causal inference in discrete settings. Journal of the American Statistical Association, 119(547):1778–1793, 2024. Evans (2016) Robin J. Evans. Graphs for margins of Bayesian networks. Scandinavian Journal of Statistics, 43(3):625–648, 2016. Galles and Pearl (1998) David Galles and Judea Pearl. An Axiomatic Characterization of Causal Counterfactuals. Foundations of Science, 3(1):151–182, 1998. Halpern (2000) Joseph Y. Halpern. Axiomatizing Causal Reasoning. Journal of Artificial Intelligence Research, 12:317–337, 2000. Kawakami et al. (2024) Yuta Kawakami, Manabu Kuroki, and Jin Tian. Probabilities of causation for continuous and vector variables. In Proceedings of the Fortieth Conference on Uncertainty in Artificial Intelligence, volume 244 of Proceedings of Machine Learning Research, pages 1901–1921. PMLR, 15–19 Jul 2024. Kersting et al. (2017) Kristian Kersting, Martin Mladenov, and Pavel Tokmakov. Relational linear programming. Artificial Intelligence, 244:188–216, 2017. Manski (1990) Charles F. Manski. Nonparametric bounds on treatment effects. The American Economic Review, 80(2):319–323, 1990. Manski (1997) Charles F. Manski. Monotone Treatment Response. Econometrica, 65(6):1311, 1997. Pearl (1995) Judea Pearl. Causal diagrams for empirical research. Biometrika, 82(4):669–688, 1995. Pearl (1999) Judea Pearl. Probabilities of causation: Three counterfactual interpretations and their identification. Synthese, 121(1):93–149, 1999. Pearl (2009) Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2009. Plečko and Bareinboim (2024) Drago Plečko and Elias Bareinboim. Causal Fairness Analysis: A Causal Toolkit for Fair Machine Learning. Foundations and Trends® in Machine Learning, 17(3):304–589, 2024. Sachs et al. (2023) Michael C. Sachs, Gustav Jonzon, Arvid Sjölander, and Erin E. Gabriel. A General Method for Deriving Tight Symbolic Bounds on Causal Effects. Journal of Computational and Graphical Statistics, 32(2):567–576, 2023. Shpitser and Pearl (2007) Ilya Shpitser and Judea Pearl. What Counterfactuals Can Be Tested. 2007. Sun et al. (2025) Hanmei Sun, Chengfeng Shi, and Qiang Zhao. Bounding the probability of causation under ordinal outcomes. Communications in Statistics-Theory and Methods, 54(24):8121–8132, 2025. Tian (2002) Jin Tian. Studies in causal reasoning and learning. PhD thesis, University of California, Los Angeles, 2002. Tian and Pearl (2000) Jin Tian and Judea Pearl. Probabilities of causation: Bounds and identification. Annals of Mathematics and Artificial Intelligence, 28(1-4):287–313, 2000. Zaffalon et al. (2023) Marco Zaffalon, Alessandro Antonucci, Rafael Cabañas, and David Huber. Approximating counterfactual bounds while fusing observational, biased and randomised data sources. International Journal of Approximate Reasoning, 162:109023, 2023. Zaffalon et al. (2024) Marco Zaffalon, Alessandro Antonucci, Rafael Cabañas, David Huber, and Dario Azzimonti. Efficient computation of counterfactual bounds. International Journal of Approximate Reasoning, 171:109111, 2024. Zhang et al. (2022) Junzhe Zhang, Jin Tian, and Elias Bareinboim. Partial counterfactual identification from observational and experimental data. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 26548–26558, 2022. Appendix A Proofs Proof of Thm. 1. Let us separately prove the non-strict inequalities between the two sides of Eq. (5). By definition, MσM_σ is based on a causal graph compatible with σ. This implies Mσ⊇σM_M_σ _σ and hence the ≤ inequality. In order to prove the opposite inequality, let ℳ∗M^* be an SCM that solves Eq. (4). By definition of σM_σ, the causal graph of ℳ∗M^* is consistent with the total order σ. Therefore, for each V(i)∈V^(i)∈ V, the endogenous parents of V(i)V^(i) in ℳ∗M^* form a subset of its predecessors ≺i V_ i under σ. We can enlarge each SE of V(i)V^(i) by including as inputs all variables in ≺i V_ i that are not already endogenous parents of V(i)V^(i), and by similarly incorporating all exogenous variables not already among its exogenous parents. In this way, the collection U of exogenous variables can be treated as a single aggregated variable U. Consequently, ℳ∗M^* may be regarded as an SCM that is consistent with the same causal graph as the PSCM MσM_σ. Since MσM_σ is canonical, there exists an equivalent SCM in MσM_M_σ. This implies the ≥ inequality. ∎ Lemma 1. For any atomic counterfactual proposition V=vV_ x=v and σ∈σ s_σ∈ S_σ, the truth value of the entailment σ⊧(V=v) s_σ (V_ x=v) is determined. Furthermore, if γ is a conjunction of such propositions, the truth value of σ⊧γ s_σ γ is also determined. Proof. Consider an arbitrary intervention x. We proceed by induction on the topological order σ to show that the potential response of every variable V(i)V^(i) is uniquely determined by σ s_σ and x. If i=1i=1, the variable V(1)V^(1) has no endogenous parents. If V(1)∈V^(1)∈ X, its value is fixed to the assignment specified in x. If V(1)∉V^(1)∉ X, its value is uniquely determined by the realization encoded in σ s_σ. Assume now that for some i>1i>1, the interventional values of all preceding variables ≺i V_ i are fixed to (≺i)( v_ i)_ x. We determine the value of V(i)V^(i) as follows: 1. If V(i)∈V^(i)∈ X, its value is determined by the assignment in x. 2. If V(i)∉V^(i)∉ X, its value is uniquely determined by σ s_σ given the values of its parents, which are fixed to (≺i)( v_ i)_ x by the induction hypothesis. Since ||<∞| V|<∞, every counterfactual variable V_ x evaluates to a unique constant under σ s_σ. Consequently, the atomic entailment is well-defined. By extension, the entailment for a conjunction γ=⋀j=1k(Y(j)(j)=y(j))γ= _j=1^k(Y^(j)_ x^(j)=y^(j)) is determined by: (σ⊧γ)⇔⋀j=1k(σ⊧Y(j)(j)=y(j)).( s_σ γ) _j=1^k ( s_σ Y^(j)_ x^(j)=y^(j) )\,. (19) ∎ Lemma 2. For any arbitrary SCM ℳ∈MσM _M_σ, there exists a bijection h:→σh:U→ S_σ between its exogenous domain U and the response signature domain σ S_σ. Proof. Let h:→σh:U→ S_σ map each exogenous state u to its induced signature σ s_σ. Because SEs are deterministic, each u induces exactly one signature, making h a well-defined function. We establish that h is a bijection. By definition, σ S_σ comprises all valid configurations of potential outcomes which MσM_σ generates. Because ℳ∈MσM _M_σ, it follows a canonical specification (Zhang et al., 2022, Def. 2.3); its exogenous domain U explicitly indexes the Cartesian product of every possible deterministic function mapping from endogenous parents to children. Thus, for any arbitrary signature σ∈σ s_σ∈ S_σ, there is at least one state u∈u that generates it. We conclude that h is surjective. By the same canonical specification of ℳM, the cardinality of U is exactly equal to the cardinality of the Cartesian product of all possible functional mappings—which coincides with |σ|| S_σ| by construction in Eq. (6). Because h is a surjective function mapping between two finite sets of identical cardinality, it must necessarily be injective. Therefore, every state u corresponds uniquely to a signature σ s_σ. ∎ Proof of Thm. 2. By Thm. 1, without loss of generality, we can restrict our analysis to an arbitrary SCM ℳ∈MσM _M_σ. By Lem. 2, there exists a bijection h mapping each exogenous state u∈u to a unique response signature σ∈σ s_σ∈ S_σ. Consequently, the logical equivalence (u⊧γ)⇔(σ⊧γ)(u γ) ( s_σ γ) holds true for σ=h(u) s_σ=h(u). The left-hand side is well-defined by the soundness axiomatisation of SCMs (Halpern, 2000), and the right-hand side is determined by Lem. 1. Because h is bijective, the probability measure over U directly translates to σ S_σ, yielding P(U=u)=P(σ=σ)P(U=u)=P( S_σ= s_σ). Substituting u with σ s_σ in Eq. (1), we obtain: Pℳ(γ)=∑σ∈σ⟦σ⊧γ⟧⋅P(σ=σ).P_M(γ)= _ s_σ∈ S_σ s_σ γ · P( S_σ= s_σ)\,. ∎ Proof of Thm. 3. Let L∗L^* be the minimum of the optimisation problem in Eq. (4), and let L be the minimum of the LP task defined in Eq. (9). To prove the theorem, we first show that L≤L∗L≤ L^*, and then L≥L∗L≥ L^*. Assume the optimal value L∗L^* is achieved by some SCM ℳ∗∈σM^* _σ for query Pℳ∗(γ)P_M^*(γ). By Thm. 1, we can consider ℳ∗M^* to be an SCM parametrising the canonical PSCM MσM_σ. By construction of σ S_σ, the exogenous distribution of ℳ∗M^* maps to a valid joint probability distribution ∗q^* over the states σ∈σ s_σ∈ S_σ. By Thm. 2, the objective query can be expressed as a linear combination of the probabilities associated with such states. Furthermore, since ℳ∗M^* is a feasible solution to Eq. (4), it must satisfy the observational constraint Pℳ∗()=P^()P_M^*( V)= P( V). Correspondingly, ∗q^* satisfies the LP’s constraint imposed by observational data as by Cor. 1. Consequently, ∗q^* constitutes a feasible solution within the constraints of the LP formulation in Eq. (9). Because the LP identifies the global minimum L over all feasible distributions, it must be that L≤L∗L≤ L^*. Conversely, let ∗q^* be the optimal solution that minimises the LP task in Eq. (9), such that the objective evaluates to L. Consider a specific SCM ℳ∈σM _σ; our argument will equivalently consider ℳ∈MσM _M_σ (after Thm. 1). We define its exogenous distribution by assigning Pℳ(u)=qh(u)∗P_M(u)=q^*_h(u), where h is the bijection established in Lem. 2. Additionally, because the LP enforces the observational constraints, ℳM matches the empirical distribution, then it is also a feasible candidate for the optimisation task in Eq. (4). Since L∗L^* represents the minimum over all feasible SCMs in σM_σ, and we have found a specific feasible SCM ℳM that achieves L, it must hold that L∗≤L^*≤ L. Since L≤L∗L≤ L^* and L∗≤L^*≤ L, we conclude that L=L∗L=L^*. Furthermore, the explicit construction in the second half of the proof demonstrates that there exists an SCM ℳ∈MσM _M_σ that precisely attains the LP bounds. ∎ Proof of Thm. 4. Let QLPQ_LP and QLLPQ_LLP denote the feasible regions of the LP in Eq. (9) and of the lifted LP in Eq. (11), respectively. For any feasible solution ∗∈QLP q^*∈ Q_LP, by definition of the equivalence relation in Eq. (10), qσ=∑σ∈σqσ∗q_ t_σ= _ s_σ∈ t_σq^*_ s_σ trivially satisfies the constraints of Eq. (11) whilst preserving the objective value, as the Boolean evaluations are constant within a class σ∈σ t_σ∈ T_σ. Conversely, consider any feasible solution σ∗∈QLLP q^*_ t_σ∈ Q_LLP. Let kσ≔|σ|k_ t_σ | t_σ| be the number of signatures σ s_σ contained within a given equivalence class σ t_σ. Note that by construction every class σ t_σ is non-empty, hence kσ>0k_ t_σ>0. We construct a valid distribution in the original space by assigning a uniform mass qσ=qσ∗/kσq_ s_σ=q^*_ t_σ/k_ t_σ to every signature σ∈σ s_σ∈ t_σ, for all σ∈σ t_σ∈ T_σ. Because the logical constraints evaluate identically across σ t_σ, summing over this uniformly distributed mass recovers the constraints and the objective value of the lifted LP, satisfying Eq. (9). Since these mappings strictly preserve both feasibility and the objective value in both directions, the extrema evaluated over QLLPQ_LLP must exactly match those over QLPQ_LP. Following the constructive argument detailed in the proof of Thm. 3, there exists a witness SCM that attains these bounds, ensuring the tightness is inherited by the lifted LP in Eq. (11). ∎ Proof of Thm. 5. Let ℳ∗∈σM^* _σ and ℳ~∗∈σγ M^* _σ V_γ denote the SCMs giving the optima in, respectively, Eq. (4) and Eq. (12). Let also L≔Pℳ∗(γ)L P_M^*(γ) and L′≔Pℳ~∗(γ)L P_ M^*(γ). Following, for instance, Bongers et al. (2021), we can marginalise out from ℳ∗M^* the variables in ∖γ V V_γ and obtain a well-defined SCM over γ V_γ. Moreover, by construction, we have Pℳ∗()=P^()P_M^*( V)= P( V), and this implies, in the marginal SCM, Pℳ∗(γ)=P^(γ)P_M^*( V_γ)= P( V_γ). Overall, we have an SCM in σγM_σ V_γ giving L, and this proves L≥L′L≥ L . To prove the inverse inequality, let P(′∣γ):=P^()/P^(γ)P( V V_γ):= P( V)/ P( V_γ), where ′:=∖γ V := V V_γ. We augment the SCM ℳ~∗ M^* with the variables in ′ V . By setting all the variables in γ V_γ as endogenous parents of those in ′ V , we can easily obtain P()=P^()P( V)= P( V). Thus we have an SCM in σM_σ such that P()=P^()P( V)= P( V) and giving L′L . This proves L≤L′L≤ L , and hence the thesis. ∎ Proof of Thm. 6. Let σ∈Σ(≺γ)σ∈ ( _γ) be an arbitrary total order compatible with ≺γ _γ, and consider the lifted LP in Eq. (11) constructed over σ. Each equivalence class σ∈σ t_σ∈ T_σ is uniquely characterised by an observational assignment ∈ v∈ V and a Boolean query evaluation b=⟦σ⊧γ⟧∈0,1b= t_σ γ ∈\0,1\. An equivalence class indexed by ⟨,b⟩ v,b is non-empty (i.e., contains at least one signature ∈σ s∈ S_σ) if and only if v and b do not contradict the consistency axiom (Galles and Pearl, 1998). Specifically, for any atomic proposition Yx=y∗Y_x=y^* in γ, if v assigns X=xX=x, consistency forces Yx=vYY_x=v_Y. Hence, if vY≠y∗v_Y≠ y^*, consistency prevents γ from being satisfied, ruling out b=1b=1; conversely, if v forces every atomic term in γ to hold, consistency rules out b=0b=0. Crucially, whether a pair ⟨,b⟩ v,b satisfies the consistency axiom depends solely on the assignment v and the syntax of γ, independently of the topological order σ. Moreover, because the canonical specification contains all possible functional responses, every pair ⟨,b⟩ v,b that is compatible with consistency is realised by at least one signature ∈σ s∈ S_σ: the observational response is fixed to v, while the counterfactual responses under unobserved interventions can be set to satisfy γ (yielding b=1b=1) or violate at least one atom of γ (yielding b=0b=0). We can conclude that the quotient space σ T_σ of non-empty equivalence classes is invariant to σ. Indexing the lifted optimisation variables as q⟨,b⟩q_ v,b , the LP objective function is ∑q⟨,1⟩ _ vq_ v,1 and the empirical constraints are q⟨,0⟩+q⟨,1⟩=P^()q_ v,0 +q_ v,1 = P( v) for each ∈ v∈ V. Because the optimisation variables, constraints, and objective function are entirely decoupled from σ, the resulting LP is identical for all linear extensions σ∈Σ(≺γ)σ∈ ( _γ). Solving the LP for any arbitrary compatible σ therefore yields the same bounds as Eq. (12) and solves the global optimisation task in Eq. (13). ∎ Appendix B Additional Experiments ACE with Imperfect Compliance. Method Lower Upper Vitamin A Balke and Pearl (1997) -0.195 0.005 This paper -0.587 0.412 Coronary Balke and Pearl (1997) 0.262 0.868 This paper -0.685 0.971 Table 3: Comparison of ACE bounds. Lower and upper bounds for the ACE grouped by Study. We consider a model with three Boolean variables. In addition to the outcome Y, we distinguish between treatment exposure X and treatment assignment Z. To evaluate the causal effect of the treatment when assigned, we focus on the query ACEX=0,X=1(Y=1)ACE_X=0,X=1(Y=1), defined as in Eq. (16). Observational information is available only in the form of a conditional empirical distribution P^(X,Y|Z) P(X,Y|Z). In our framework, consistency with the empirical distribution is expressed as: ∑σ∈σ⟦σ⊧(x,y,z)⟧⋅qσ∑σ∈σ⟦σ⊧z⟧⋅qσ=P^(x,y|z), _ s_σ∈ S_σ s_σ (x,y,z) · q_ s_σ _ s_σ∈ S_σ s_σ z · q_ s_σ= P(x,y|z)\,, (20) which yields linear constraints over the optimisation variables qσq_ s_σ, thus preserving the LP structure of our method. Although the query involves only X and Y, the conditional nature of the empirical distribution prevents a marginalisation of Z. Consequently, the response signature must be based on all three Boolean variables. By Thm. 6 we choose an arbitrary admissible order σ such that X≺YX Y. We compare our bounds with those obtained from the LP formulation introduced by Balke and Pearl (1997) and later generalised by Sachs et al. (2023). That approach requires an explicit causal graph. In particular, the model in Fig. 1(b) assumes that Z affects Y only through X, and that Z is independent of an unobserved confounder U affecting both X and Y. Tab. 3 reports the resulting comparison for two observational studies: the vitamin A supplementation study in northern Sumatra and the coronary primary prevention trial. Our bounds remain informative, but are substantially wider than those of Balke and Pearl (1997). This difference reflects the fact that our framework relies only on the order restrictions induced by the query, whereas the graphical model in Fig. 1(b) imposes additional structure, most notably the exclusion of a direct effect from Z to Y. To model these effects with such a granularity one would either impose additional (non-linear) constraints or eventually require a causal graph, but this lies outside manuscript’s scope. By Thm. 3, the more extreme values admitted by our procedure are attainable under SCMs compatible with our weaker assumptions, including models in which Z may directly affect Y. As pointed out by Balke and Pearl (1994a), bounds can be further tightened if additional assumptions are made about subject behaviours: imposing a monotonic constraint, e.g. a positive response on treatment, Yx1,z≥Yx0,zY_x_1,z≥ Y_x_0,z for every z∈z . As a consequence, negative lower bounds on the effects are raised to zero. Nevertheless, even when combined with additional assumptions such as weak exogeneity, these restrictions are still insufficient to recover the analytical bounds implied by the graphical model, as discussed in Sect. 5. PNS Generalisation for Non-Binary Variables 44556677889910101111121213131414151510−110^-110010^010110^110210^2Cardinality nnComputation Time (seconds)GPNS(2)GPNS(3)GPNS(4) Figure 2: End-to-end computation time for GPNS(r)GPNS(r) bounds, on a semi-logarithmic scale, as a function of the variable cardinalities (n=||=||n=|X|=|Y|). To assess the practical scalability of the lifted LP in Eq. (11), we measured the end-to-end time required to construct the signature equivalence classes, aggregate the lifted variables, and solve the resulting LP. The benchmark query is a multi-hypothetical generalisation of PNS for ordinal treatment and outcome variables, in the spirit of Kawakami et al. (2024). Let X and Y be ordinal treatment and effect variables, respectively, defined over finite domains =x(0),…,x(n−1)X=\x^(0),…,x^(n-1)\ and =y(0),…,y(n−1)Y=\y^(0),…,y^(n-1)\ where n is an experimental parameter. The GPNS of order r (where 2≤r≤n2≤ r≤ n) evaluates the joint probability of strict compliance across r distinct interventions, defined as: GPNS(r)≔P(⋀i=0r−1Yx(i)=y(i)).GPNS(r) P ( _i=0^r-1Y_x^(i)=y^(i) )\,. (21) For binary X and Y, GPNS(2)GPNS(2) reduces to the standard PNS event. For each n∈4,…,15n∈\4,…,15\, we randomly parametrise a Bayesian network with the endogenous structure in Fig. 1(a), estimated P^(X,Y) P(X,Y) from simulated observational data, and computed bounds for GPNS(r)GPNS(r) with r∈2,3,4r∈\2,3,4\. Figure 2 reports the resulting end-to-end computation times. For this benchmark, γ=X,Y V_γ=\X,Y\, so the lifted LP has at most ||||=n2|X||Y|=n^2 variables, as derived in Sect. 3. Results are consistent with the polynomial scaling of the lifted formulation, although we must consider the overhead introduced by the symmetry detection, LP construction, and numerical optimisation. As expected, solving the LP remains tractable for standard solvers as the domain expands. Appendix C Canonical Specifications XXWWYYUXU_XUWU_WUYU_Y (a) XXWWYYUXU_XUWYU_WY (b) XXWWYYUU (c) Figure 3: Examples of canonical PSCM specifications compatible with the endogenous order σ≔(X,W,Y)σ (X,W,Y). (a) Markovian specification. (b) Specification with one exogenous variable shared by W and Y. (c) Specification with one exogenous variable shared by all endogenous variables. Grey nodes denote exogenous variables. We expand the discussion of canonicalisation of SCMs and PSCMs briefly introduced in Sect. 3. We recall here that the canonical PSCM used in the proof of Thm. 1 builds upon the theoretical contributions of Zhang et al. (2022) and Balke and Pearl (1994a). Generally speaking, an SCM takes its structural equations (SEs) as given. When these are unknown and the endogenous variables take on finitely many values, the literature shows that one may still specify a PSCM from the causal graph alone by adopting a canonical representation. For an endogenous variable V∈V∈ V with parents PaVPa_V, a state of its exogenous parent UVU_V acts as a selector for a deterministic mapping (or response function) from the domain of PaVPa_V to V. There are exactly |||PaV||V|^|Pa_V| such possible functions, where |PaV|≔∏Z∈PaV|||Pa_V| _Z _V|Z|. In the absence of unobserved confounding, as in Markovian models, each UVU_V is independent of all other exogenous variables and can be replaced without loss of generality by a discrete canonical variable of finite cardinality |V|=|||PaV||U_V|=|V|^|Pa_V|. When latent confounding is present, as in semi-Markovian models, unobserved variables may jointly affect multiple endogenous variables. The graph is then partitioned into confounded components (or c-components) (Tian, 2002), defined as the maximal sets of endogenous variables connected by paths of latent confounders (see confounded paths in Shpitser and Pearl, 2007). In a canonical specification, all latent confounding within each c-component ⊆ C V is replaced without loss of generality by a single discrete exogenous variable U_ C that jointly selects the response functions for every variable in C, whilst the exogenous variables associated with distinct c-components remain mutually independent. Because the state space of U_ C must encompass all joint assignments of response functions, its canonical space U_ C indexes all possible combinations of individual mappings, yielding ||=∏V∈|||PaV||U_ C|= _V∈ C|V|^|Pa_V|. Canonicalisation guarantees that any discrete SCM (or PSCM) admits a discrete equivalent inducing an identical set of observational and counterfactual distributions (Evans, 2016; Duarte et al., 2024). To ground the ideas, we proceed by means of an example. Consider ≔X,W,Y V \X,W,Y\ with binary domains ==0,1X=Y=\0,1\ and ternary domain =0,1,2W=\0,1,2\. We assume the causal order X≺W≺YX W Y. The three panels in Fig. 3 share the same endogenous structure compatible with this order. A canonical specification of these causal diagrams is obtained by enumerating the deterministic response functions for each endogenous variable. First, consider the unconfounded setting in Fig. 3(a), where every variable forms a singleton c-component. Since PaX=∅Pa_X= , there are two constant mappings, giving |X|=2|U_X|=2. For W, whose parent is X, the mechanism w←fW(x,uW)w← f_W(x,u_W) assigns a state in W to each x∈x ; hence, |W|=||||=32=9|U_W|=|W|^|X|=3^2=9. Finally, the mapping y←fY(x,w,uY)y← f_Y(x,w,u_Y) assigns a value in Y to every configuration (x,w)∈×(x,w) ×W, yielding |Y|=||||⋅||=22⋅3=64|U_Y|=|Y|^|X|·|W|=2^2· 3=64. Next, consider Fig. 3(b), where unobserved confounding between W and Y induces the partition into c-components X\X\ and W,Y\W,Y\. The canonical confounder UWYU_WY jointly indexes the mappings of both W and Y, yielding |WY|=|W|⋅|Y|=9⋅64=576|U_WY|=|U_W|·|U_Y|=9· 64=576, whilst |X|=2|U_X|=2. Lastly, in the fully confounded model of Fig. 3(c), all endogenous variables belong to a single c-component = C= V. A single global exogenous parent U governs all mechanisms, with cardinality ||=|X|⋅|W|⋅|Y|=2⋅9⋅64=1,152|U|=|U_X|·|U_W|·|U_Y|=2· 9· 64=1,152. The canonical PSCM constructed in the proof of Thm. 1 matches this full-confounding specification.