Introduction
Alarms in industrial and scientific systems demand decisions before their causes are clear. Immediately after an alarm, incident data are scarce and causal mechanisms may have changed. Ranking variables by observed shift is unreliable: a large shift may be downstream or have little effect on the alarm score, whereas a modest upstream shift may propagate through the system and strongly affect the score.
The diagnostic question is therefore how changes in the data-generating process propagated through the system to move a supplied alarm score, rather than merely which variables shifted. We seek an additive allocation of this contrast to structural-function and exogenous-law changes. Producing it soon enough to guide intervention requires three coupled capabilities. First, the report must accommodate a changing SCM and resolve parent-dependent mechanisms from exogenous inputs; otherwise opposing effects at one node can cancel and obscure which kind of repair could help. Second, to compute the report quickly, its Shapley allocation must avoid factorial permutation averaging and exponential coalition enumeration. Third, to work from few incident samples, it must leverage the reference regime before enough incident data accumulate to relearn the SCM. Resolution, computation, and estimation thus jointly determine whether attribution can guide intervention immediately after an alarm.
Existing methods solve only subsets of this immediate-diagnosis problem. Feature-attribution methods explain a fixed score map through observed inputs [16 ] ; event-level causal attribution explains an individual prediction or outlier within one model [11 , 6 , 20 ] ; node-wise population attribution replaces a whole causal conditional [4 , 35 , 25 ] ; and shift discovery locates changed components [12 , 7 , 28 ] . The first two hold generating mechanisms fixed. Node-wise population reports compare regimes, but they can require substantial target-regime data to estimate conditional changes or hybrid values and can merge opposing function/law effects. Shift discovery stops before allocating a supplied score contrast after propagation.
Locating where the system changed and determining what moved the alarm score are therefore distinct tasks; none of these approaches provides the resolved, tractable, and scarce-data attribution needed here.
Our starting point is a unified Shapley game over hybrid SCM states, indexed by source, target, player granularity, and score. Each coalition takes selected components from the incident regime and the rest from the reference regime, then evaluates the score after causal propagation. The primary population report uses structural functions and exogenous laws as separate players, so their propagated contributions remain visible. The same framework defines realized-event reports, with an event–population coherence identity linking their expected attributions to the population report. Because the allocation is relative to aligned endpoint SCMs, physical interpretation still requires system-specific semantics.
To make this estimand available in the short incident window, we specialize to ordered linear SCMs and quadratic scores. Every hybrid expectation then depends only on causally propagated first and second moments. This includes common PCA- and ICA-based monitoring scores, quadratic scores formed from Mahalanobis or graphical-lasso precision matrices, and squared z z z [24 , 14 , 17 , 8 ] . This second-order closure resolves both remaining bottlenecks: a causal-order moment recursion, shared backward-sensitivity sweep, and exact quadrature compute the full Shapley vector exactly up to roundoff in polynomial time, while a well-observed reference SCM anchors cross-fitted estimates of incident mechanism deviations and residual moments. The resulting player-specific error bounds yield a uniform rate, sign recovery, and top-k k k
The resulting method is S4 {}^{4} 4
Exact polynomial-time computation. The moment program computes the full function/law vector, up to roundoff, in O ( d 4 ) O(d^{4}) O ( d 4 ) O ( d 2 ) O(d^{2}) O ( d 2 )
Scarce-incident estimation. Reference-anchored cross-fitting yields player-specific error bounds and guarantees for sign and top-k k k
Unified Causal Score Attribution
Source, target, players, and score. The primary attribution problem has four objects: (i) a normal-regime source world , represented by an SCM M ( 0 ) M^{(0)} M ( 0 ) target world , represented by M ( 1 ) M^{(1)} M ( 1 ) U \mathcal{U} U τ \tau τ T T T ( M ( 0 ) , M ( 1 ) , U , T ) (M^{(0)}, M^{(1)}, \mathcal{U}, T) ( M ( 0 ) , M ( 1 ) , U , T )
Let X ⊆ R d \mathcal{X} \subseteq\mathbb{R}^{d} X ⊆ R d x = ( x j ) j = 1 d ∈ X x = (x_j)_{j=1}^{d} \in\mathcal{X} x = ( x j ) j = 1 d ∈ X e = ( e j ) j = 1 d e = (e_j)_{j=1}^{d} e = ( e j ) j = 1 d r ∈ { 0 , 1 } r \in\{0,1\} r ∈ { 0 , 1 } F ( r ) = ( F j ( r ) ) j = 1 d F^{(r)} = (F_j^{(r)})_{j=1}^{d} F ( r ) = ( F j ( r ) ) j = 1 d Q ( r ) = ⨂ j = 1 d Q j ( r ) Q^{(r)} = \bigotimes_{j=1}^{d} Q_j^{(r)} Q ( r ) = ⨂ j = 1 d Q j ( r ) M ( r ) = M ( F ( r ) , Q ( r ) ) M^{(r)} = M(F^{(r)}, Q^{(r)}) M ( r ) = M ( F ( r ) , Q ( r ) )
x j = F j ( r ) ( x < j , e j ) , e j ∼ Q j ( r ) . x_j = F_j^{(r)}(x_{<j}, e_j), \qquad e_j \sim Q_j^{(r)}. x j = F j ( r ) ( x < j , e j ) , e j ∼ Q j ( r ) . The exogenous coordinates are mutually independent, and x = G F ( r ) ( e ) x = G_{F^{(r)}}(e) x = G F ( r ) ( e ) τ : X → R \tau: \mathcal{X} \to\mathbb{R} τ : X → R
T ( M ( F , Q ) ) : = E e ∼ Q [ τ ( G F ( e ) ) ] . (1) T(M(F,Q)) := \mathbb{E}_{e \sim Q}[\tau(G_F(e))]. \tag*{(1)} T ( M ( F , Q )) := E e ∼ Q [ τ ( G F ( e ))] . (1) The primary population question attributes T ( M ( 1 ) ) − T ( M ( 0 ) ) T(M^{(1)}) - T(M^{(0)}) T ( M ( 1 ) ) − T ( M ( 0 ) )
SCM-state players to attribute. Let A = { F j , Q j : j ∈ [ d ] } \mathcal{A} = \{F_j, Q_j : j \in[d]\} A = { F j , Q j : j ∈ [ d ]} F j F_j F j F j ( 0 ) F_j^{(0)} F j ( 0 ) F j ( 1 ) F_j^{(1)} F j ( 1 ) Q j Q_j Q j Q j ( 0 ) Q_j^{(0)} Q j ( 0 ) Q j ( 1 ) Q_j^{(1)} Q j ( 1 ) u u u A \mathcal{A} A U \mathcal{U} U u u u node-local if u ⊆ { F j , Q j } u \subseteq\{F_j,Q_j\} u ⊆ { F j , Q j } j j j U = U F ∪ Q \mathcal{U} = \mathcal{U}_{F\cup Q} U = U F ∪ Q
U F ∪ Q = { { F j } : j ∈ [ d ] } ∪ { { Q j } : j ∈ [ d ] } . \mathcal{U}_{F\cup Q} = \{\{F_j\}:j \in[d]\} \cup\{\{Q_j\}:j \in[d]\}. U F ∪ Q = {{ F j } : j ∈ [ d ]} ∪ {{ Q j } : j ∈ [ d ]} . Its 2 d 2d 2 d U ( F , Q ) = { { F j , Q j } : j ∈ [ d ] } \mathcal{U}_{(F,Q)} = \{\{F_j,Q_j\}:j \in[d]\} U ( F , Q ) = {{ F j , Q j } : j ∈ [ d ]} U Q = { { Q j } : j ∈ [ d ] } \mathcal{U}_Q = \{\{Q_j\}:j \in[d]\} U Q = {{ Q j } : j ∈ [ d ]}
A coalition is U ∈ 2 U U \in2^{\mathcal{U}} U ∈ 2 U U ⊆ U U \subseteq\mathcal{U} U ⊆ U u ∈ U u \in\mathcal{U} u ∈ U U U U M U = M ( F U , Q U ) M^U = M(F^U,Q^U) M U = M ( F U , Q U ) M ( 1 ) M^{(1)} M ( 1 ) M ( 0 ) M^{(0)} M ( 0 )
F j U = F j 1 { F j ∈ ⋃ w ∈ U w } , Q j U = Q j 1 { Q j ∈ ⋃ w ∈ U w } . F_j^U = F_j^{1\{\!F_j\in\bigcup_{w\in U}w\!\}}, \qquad Q_j^U = Q_j^{1\{\!Q_j\in\bigcup_{w\in U}w\!\}}. F j U = F j 1 { F j ∈ ⋃ w ∈ U w } , Q j U = Q j 1 { Q j ∈ ⋃ w ∈ U w } . Each coalition is therefore an SCM replacement followed by causal propagation, not feature masking. Switching the label pair { F j , Q j } \{F_j,Q_j\} { F j , Q j }
Causal Shapley estimand. The causal question is how much each prespecified SCM-state player contributes to the score contrast between the source and target worlds. Each coalition U U U v : 2 U → R v:2^{\mathcal{U}}\to\mathbb{R} v : 2 U → R v ( U ) = T ( M U ) v(U)=T(M^U) v ( U ) = T ( M U ) u u u [29 ] . Denote it by Φ u ( M ( 0 ) , M ( 1 ) , U , T ) \Phi_u(M^{(0)},M^{(1)},\mathcal{U},T) Φ u ( M ( 0 ) , M ( 1 ) , U , T ) Φ u \Phi_u Φ u
Φ u : = ∑ U ⊆ U ∖ { u } ∣ U ∣ ! ( ∣ U ∣ − ∣ U ∣ − 1 ) ! ∣ U ∣ ! { T ( M U ∪ { u } ) − T ( M U ) } . (2) \Phi_u := \sum_{U\subseteq\mathcal{U}\setminus\{u\}} \frac{|U|!(|\mathcal{U}|-|U|-1)!}{|\mathcal{U}|!}\left\{T(M^{U\cup\{u\}})-T(M^U)\right\}. \tag*{(2)} Φ u := U ⊆ U ∖ { u } ∑ ∣ U ∣ ! ∣ U ∣ ! ( ∣ U ∣ − ∣ U ∣ − 1 )! { T ( M U ∪ { u } ) − T ( M U ) } . (2) For singleton players, we suppress braces in attribution subscripts, writing Φ F j : = Φ { F j } \Phi_{F_j}:=\Phi_{\{F_j\}} Φ F j := Φ { F j } Q j Q_j Q j U = U F ∪ Q \mathcal{U}=\mathcal{U}_{F\cup Q} U = U F ∪ Q ( Φ F 1 , … , Φ F d , Φ Q 1 , … , Φ Q d ) (\Phi_{F_1},\ldots,\Phi_{F_d},\Phi_{Q_1},\ldots,\Phi_{Q_d}) ( Φ F 1 , … , Φ F d , Φ Q 1 , … , Φ Q d )
Report variants and granularity. For an event x x x e ( r ) ( x ; F ( r ) ) e^{(r)}(x;F^{(r)}) e ( r ) ( x ; F ( r ) ) r r r δ e ( r ) \delta_{e^{(r)}} δ e ( r ) Φ u \Phi_u Φ u Φ u s h i f t ( e ( 1 ) ) \Phi_u^{\mathrm{shift}}(e^{(1)}) Φ u shift ( e ( 1 ) ) Φ u f i x ( e ( 0 ) ) \Phi_u^{\mathrm{fix}}(e^{(0)}) Φ u fix ( e ( 0 ) ) F ( 0 ) F^{(0)} F ( 0 ) F ( 1 ) F^{(1)} F ( 1 ) [22 ] .
Causal replacement semantics. Activating F j F_j F j Q j Q_j Q j
Proposition 2.1 (Efficiency and event–population coherence). Let Φ u \Phi_u Φ u Φ u s h i f t ( e ) \Phi_u^{\mathrm{shift}}(e) Φ u shift ( e ) U \mathcal{U} U M ( 1 ) M^{(1)} M ( 1 ) e ∼ Q ( 1 ) e \sim Q^{(1)} e ∼ Q ( 1 )
∑ u ∈ U Φ u = T ( M ( 1 ) ) − T ( M ( 0 ) ) , (3) \sum_{u\in\mathcal{U}} \Phi_u = T(M^{(1)}) - T(M^{(0)}), \tag*{(3)} u ∈ U ∑ Φ u = T ( M ( 1 ) ) − T ( M ( 0 ) ) , (3) E e ∼ Q ( 1 ) [ Φ u s h i f t ( e ) ] = Φ u for every u ∈ U . (4) \mathbb{E}_{e\sim Q^{(1)}}\left[\Phi_u^{\mathrm{shift}}(e)\right] = \Phi_u \qquad\text{for every } u\in\mathcal{U}. \tag*{(4)} E e ∼ Q ( 1 ) [ Φ u shift ( e ) ] = Φ u for every u ∈ U . (4) The efficiency identity makes the population diagnosis complete for the selected score contrast: its component attributions account for the entire mean-score change. Event–population coherence makes the two causal estimands directly comparable: averaging shifted-event attributions over incident events recovers each population attribution exactly, provided the endpoints, players, and replacement semantics agree. The proof is in Supplement Sec. C.1; the fixed-SCM event corollary is in Sec. C.3.
Relation to causal attribution and shift discovery. The source–target–player–score framework makes the distinctions from related methods explicit. Feature and fixed-SCM event attribution explain inputs, disturbances, or interventions under one supplied model, rather than changes between regime SCMs [16 , 6 , 11 , 31 , 20 ] . do-Shapley attributes interventional outcome changes to variable-intervention players within a fixed SCM, not function/law replacements between regime SCMs [13 , 34 ] . Unit-level change attribution decomposes individual outcome differences, whereas causal change attribution (CCA) and multiply-robust CCA (MRCCA) replace node-wise generating distributions; neither provides function/law-resolved population alarm-score attribution [5 , 4 , 35 , 25 ] . Mechanism-shift and difference-DAG methods detect shifted mechanisms or structural differences, or leverage shifts to recover causal structure, but do not allocate a supplied score contrast after causal propagation [12 , 23 , 33 , 18 , 7 ] . Source-to-target regression transfers predictive models but does not define or compute this causal attribution estimand [1 , 32 , 15 ] . S 4 S^4 S 4
Why function/law resolution matters. Separating function and law effects matters for intervention and repair because the two can favor different actions at the same node. In a controlled oracle construction, a target t t t ( Φ F t , Φ Q t ) = ( − β , + α ) (\Phi_{F_t}, \Phi_{Q_t}) = (-\beta, +\alpha) ( Φ F t , Φ Q t ) = ( − β , + α ) δ \delta δ α > β > δ > α − β > 0 \alpha> \beta> \delta> \alpha- \beta> 0 α > β > δ > α − β > 0 [28 ] , but does not attribute their propagated contributions to a supplied score contrast; the two tasks are complementary.
Exact Second-Order Causal Attribution
Second-order specialization. To support immediate response, we specialize the Section 2 estimand to ordered linear SCMs and quadratic scores (polynomials of degree at most two). Their population coalition values depend only on first and second exogenous moments. The estimand remains Φ U ( M ( 0 ) , M ( 1 ) , U , T ) \Phi_{\mathcal{U}}(M^{(0)}, M^{(1)}, \mathcal{U}, T) Φ U ( M ( 0 ) , M ( 1 ) , U , T ) M ( r ) = M ( B ( r ) , Q ( r ) ) M^{(r)} = M(B^{(r)}, Q^{(r)}) M ( r ) = M ( B ( r ) , Q ( r ) ) r ∈ { 0 , 1 } r \in\{0,1\} r ∈ { 0 , 1 }
x = ( B ( r ) ) T x + e , e ∼ Q ( r ) . (5) x = (B^{(r)})^{\mathsf{T}}x + e,\qquad e \sim Q^{(r)}. \tag*{(5)} x = ( B ( r ) ) T x + e , e ∼ Q ( r ) . (5) where B ( r ) B^{(r)} B ( r ) F j F_j F j B : , j ( r ) B^{(r)}_{:,j} B : , j ( r ) Q j Q_j Q j e j e_j e j μ j ( r ) = E Q ( r ) [ e j ] \mu_j^{(r)} = \mathbb{E}_{Q^{(r)}}[e_j] μ j ( r ) = E Q ( r ) [ e j ] ω j ( r ) = Var Q ( r ) ( e j ) \omega_j^{(r)} = \operatorname{Var}_{Q^{(r)}}(e_j) ω j ( r ) = Var Q ( r ) ( e j )
A ( r ) = ( I − B ( r ) ) − 1 , x = ( A ( r ) ) T e . (6) A^{(r)} = (I - B^{(r)})^{-1},\qquad x = (A^{(r)})^{\mathsf{T}}e. \tag*{(6)} A ( r ) = ( I − B ( r ) ) − 1 , x = ( A ( r ) ) T e . (6) The supplied anomaly score is
τ ( x ) = x T W x + z T x . (7) \tau(x) = x^{\mathsf{T}}Wx + z^{\mathsf{T}}x. \tag*{(7)} τ ( x ) = x T W x + z T x . (7) We omit additive constants, which cancel from all coalition marginals and endpoint contrasts, and take W = W T W = W^{\mathsf{T}} W = W T ( W + W T ) / 2 (W + W^{\mathsf{T}})/2 ( W + W T ) /2 T 2 T^2 T 2 I 2 I^2 I 2 z z z x j x_j x j Diag ( v ) \operatorname{Diag}(v) Diag ( v ) v v v diag ( M ) \operatorname{diag}(M) diag ( M ) M M M Ω = Diag ( ω ) \Omega= \operatorname{Diag}(\omega) Ω = Diag ( ω ) V = A W A T V = AWA^{\mathsf{T}} V = A W A T ℓ = A z \ell= Az ℓ = A z x = A T e x = A^{\mathsf{T}}e x = A T e
τ ( A T e ) = e T V e + ℓ T e . (8) \tau(A^{\mathsf{T}}e) = e^{\mathsf{T}}Ve + \ell^{\mathsf{T}}e. \tag*{(8)} τ ( A T e ) = e T V e + ℓ T e . (8) Taking expectations gives the closed-form target functional:
T ( M ( B , Q ) ) = tr ( V Ω ) + μ T V μ + ℓ T μ . (9) T(M(B,Q)) = \operatorname{tr}(V\Omega) + \mu^{\mathsf{T}}V\mu+ \ell^{\mathsf{T}}\mu. \tag*{(9)} T ( M ( B , Q )) = tr ( V Ω ) + μ T V μ + ℓ T μ . (9) We also write this same value as T ( B , μ , ω ) T(B,\mu,\omega) T ( B , μ , ω )
Hybrid values. We use the primary function/law-resolved granularity U = U F ∪ U Q \mathcal{U} = \mathcal{U}_F \mathbin{\cup} \mathcal{U}_Q U = U F ∪ U Q U ⊆ U U \subseteq\mathcal{U} U ⊆ U T ( M U ) = T ( B U , μ U , ω U ) T(M^U) = T(B^U,\mu^U,\omega^U) T ( M U ) = T ( B U , μ U , ω U )
B : , j U = B : , j ( 1 ) 1 { F j ∈ ⋃ w ∈ U w } + B : , j ( 0 ) 1 { F j ∉ ⋃ w ∈ U w } , B^U_{:,j} = B^{(1)}_{:,j}\mathbf{1}\left\{F_j \in\bigcup_{w\in U}w\right\} + B^{(0)}_{:,j}\mathbf{1}\left\{F_j \notin\bigcup_{w\in U}w\right\}, B : , j U = B : , j ( 1 ) 1 { F j ∈ w ∈ U ⋃ w } + B : , j ( 0 ) 1 { F j ∈ / w ∈ U ⋃ w } , μ j U = μ j ( 1 ) 1 { Q j ∈ ⋃ w ∈ U w } + μ j ( 0 ) 1 { Q j ∉ ⋃ w ∈ U w } , \mu^U_j = \mu^{(1)}_j\mathbf{1}\left\{Q_j \in\bigcup_{w\in U}w\right\} + \mu^{(0)}_j\mathbf{1}\left\{Q_j \notin\bigcup_{w\in U}w\right\}, μ j U = μ j ( 1 ) 1 { Q j ∈ w ∈ U ⋃ w } + μ j ( 0 ) 1 { Q j ∈ / w ∈ U ⋃ w } , ω j U = ω j ( 1 ) 1 { Q j ∈ ⋃ w ∈ U w } + ω j ( 0 ) 1 { Q j ∉ ⋃ w ∈ U w } . (10) \omega^U_j = \omega^{(1)}_j\mathbf{1}\left\{Q_j \in\bigcup_{w\in U}w\right\} + \omega^{(0)}_j\mathbf{1}\left\{Q_j \notin\bigcup_{w\in U}w\right\}. \tag*{(10)} ω j U = ω j ( 1 ) 1 { Q j ∈ w ∈ U ⋃ w } + ω j ( 0 ) 1 { Q j ∈ / w ∈ U ⋃ w } . (10) The game is simulation-free, but direct Shapley evaluation still requires all K ! K! K ! 2 K 2^K 2 K K = ∣ U ∣ K = |\mathcal{U}| K = ∣ U ∣ R p ⊆ U R_p \subseteq\mathcal{U} R p ⊆ U u ∈ U u \in\mathcal{U} u ∈ U p u p_u p u U ↦ T ( M U ) U \mapsto T(M^U) U ↦ T ( M U ) [21 ]
T ~ ( p ) : = E R p [ T ( M R p ) ] = ∑ U ⊆ U T ( M U ) Pr ( R p = U ) , (11) \widetilde{T}(p) := \mathbb{E}_{R_p}\left[T(M^{R_p})\right] = \sum_{U\subseteq\mathcal{U}} T(M^U)\Pr(R_p = U), \tag*{(11)} T ( p ) := E R p [ T ( M R p ) ] = U ⊆ U ∑ T ( M U ) Pr ( R p = U ) , (11) for p = ( p u ) u ∈ U ∈ [ 0 , 1 ] K p = (p_u)_{u\in\mathcal{U}} \in[0,1]^K p = ( p u ) u ∈ U ∈ [ 0 , 1 ] K [29 , 21 ] ), but the extension and its derivatives must still be evaluated without constructing the exponential coalition table.
Moment recursion. Draw R p R_p R p x x x x ˉ ( p ) = E [ x ] \bar{x}(p)=\mathbb{E}[x] x ˉ ( p ) = E [ x ] S ( p ) = E [ x x ⊤ ] S(p)=\mathbb{E}[xx^{\top}] S ( p ) = E [ x x ⊤ ] T ~ \widetilde{T} T b j ( r ) = B < j , j ( r ) ∈ R j − 1 b_j^{(r)}=B^{(r)}_{<j,j}\in\mathbb{R}^{j-1} b j ( r ) = B < j , j ( r ) ∈ R j − 1 Z F , j = 1 { F j ∈ ⋃ w ∈ R p w } Z_{F,j}=\mathbf{1}\{F_j\in\bigcup_{w\in R_p}w\} Z F , j = 1 { F j ∈ ⋃ w ∈ R p w } Z Q , j = 1 { Q j ∈ ⋃ w ∈ R p w } Z_{Q,j}=\mathbf{1}\{Q_j\in\bigcup_{w\in R_p}w\} Z Q , j = 1 { Q j ∈ ⋃ w ∈ R p w }
b ˉ j ( p ) = E [ b j Z F , j ] , μ ˉ j ( p ) = E [ μ j Z Q , j ] , Γ j ( p ) = E [ b j Z F , j ( b j Z F , j ) ⊤ ] , c j ( p ) = E [ b j Z F , j μ j Z Q , j ] , q j ( p ) = E [ ( μ j Z Q , j ) 2 + ω j Z Q , j ] . (12) \begin{aligned}
\bar{b}_j(p)&=\mathbb{E}[b_j^{Z_{F,j}}], & \bar{\mu}_j(p)&=\mathbb{E}[\mu_j^{Z_{Q,j}}],\\
\Gamma_j(p)&=\mathbb{E}[b_j^{Z_{F,j}}(b_j^{Z_{F,j}})^{\top}], & c_j(p)&=\mathbb{E}[b_j^{Z_{F,j}}\mu_j^{Z_{Q,j}}],\\
q_j(p)&=\mathbb{E}[(\mu_j^{Z_{Q,j}})^2+\omega_j^{Z_{Q,j}}]. \tag*{(12)}
\end{aligned} b ˉ j ( p ) Γ j ( p ) q j ( p ) = E [ b j Z F , j ] , = E [ b j Z F , j ( b j Z F , j ) ⊤ ] , = E [( μ j Z Q , j ) 2 + ω j Z Q , j ] . μ ˉ j ( p ) c j ( p ) = E [ μ j Z Q , j ] , = E [ b j Z F , j μ j Z Q , j ] , (12) Because the SCM is ordered, the required moments can be propagated node by node.
Theorem 3.1 (Exact second-order moment recursion). Assume the ordered source and target states have independent exogenous coordinates with finite second moments, U U U x ˉ ( p ) \bar{x}(p) x ˉ ( p ) S ( p ) S(p) S ( p )
x ˉ j ( p ) = b ˉ j ( p ) ⊤ x ˉ < j ( p ) + μ ˉ j ( p ) , S < j , j ( p ) = S < j , < j ( p ) b ˉ j ( p ) + x ˉ < j ( p ) μ ˉ j ( p ) , S j j ( p ) = tr { Γ j ( p ) S < j , < j ( p ) } + 2 c j ( p ) ⊤ x ˉ < j ( p ) + q j ( p ) . (13) \begin{aligned}
\bar{x}_j(p)&=\bar{b}_j(p)^{\top}\bar{x}_{<j}(p)+\bar{\mu}_j(p),\\
S_{<j,j}(p)&=S_{<j,<j}(p)\bar{b}_j(p)+\bar{x}_{<j}(p)\bar{\mu}_j(p),\\
S_{jj}(p)&=\operatorname{tr}\{\Gamma_j(p)S_{<j,<j}(p)\}\\
&\quad+2c_j(p)^{\top}\bar{x}_{<j}(p)+q_j(p). \tag*{(13)}
\end{aligned} x ˉ j ( p ) S < j , j ( p ) S j j ( p ) = b ˉ j ( p ) ⊤ x ˉ < j ( p ) + μ ˉ j ( p ) , = S < j , < j ( p ) b ˉ j ( p ) + x ˉ < j ( p ) μ ˉ j ( p ) , = tr { Γ j ( p ) S < j , < j ( p )} + 2 c j ( p ) ⊤ x ˉ < j ( p ) + q j ( p ) . (13) Then, for every u ∈ U u\in U u ∈ U
Φ u ( M ( 0 ) , M ( 1 ) , U , T ) = ∫ 0 1 ∂ T ~ ( p ) ∂ p u ∣ p = t 1 d t , (14) \Phi_u(M^{(0)},M^{(1)},U,T)=\int_0^1\left.\frac{\partial\widetilde{T}(p)}{\partial p_u}\right|_{p=t\mathbf{1}}\,dt, \tag*{(14)} Φ u ( M ( 0 ) , M ( 1 ) , U , T ) = ∫ 0 1 ∂ p u ∂ T ( p ) p = t 1 d t , (14) where
T ~ ( p ) = tr { W S ( p ) } + z ⊤ x ˉ ( p ) . (15) \widetilde{T}(p)=\operatorname{tr}\{WS(p)\}+z^{\top}\bar{x}(p). \tag*{(15)} T ( p ) = tr { W S ( p )} + z ⊤ x ˉ ( p ) . (15) These equations compute every Shapley value without enumerating coalitions.
The key is second-order closure: random function/law replacements enter through the local moments in (12) , and the causal recursion propagates exactly the moments required by the score rather than approximating the coalition game.
Corollary 3.2 (Exact all-component computation). For the primary function/law-resolved game, one forward evaluation of the recursion in (13) followed by one shared backward-sensitivity sweep returns ∇ p T ~ ( p ) \nabla_p\widetilde{T}(p) ∇ p T ( p ) K = 2 d K=2d K = 2 d O ( d 3 ) O(d^3) O ( d 3 ) O ( d 2 ) O(d^2) O ( d 2 ) d d d O ( K d 3 ) = O ( d 4 ) O(Kd^3)=O(d^4) O ( K d 3 ) = O ( d 4 ) O ( d 2 ) O(d^2) O ( d 2 )
At each of the d d d 2 d 2d 2 d 2 d − 1 2d-1 2 d − 1 d d d [10 ] . Supplement Sec. C.2 and Algorithm 1 give the reverse equations, proof, and full procedure. The resulting speedup preserves the resolved estimand: tractability does not require pairing function and law players.
Event corollary. Fixing F F F [6 ] : it shares factual-noise replacement but uses E M [ τ ~ ( X ) ] \mathbb{E}_M[\widetilde{\tau}(X)] E M [ τ ( X )] F F F
Sample-Efficient Plug-in Estimation
Given the supplied score, causal order or parent sets, and player set, we estimate the Section 3 target by plugging sample estimates of B ( r ) , μ ( r ) , ω ( r ) B^{(r)}, \mu^{(r)}, \omega^{(r)} B ( r ) , μ ( r ) , ω ( r ) τ \tau τ
Identification for the linear plug-in. Under (5), a correct common order (or parent sets) and nonsingular parent covariances identify B : , j ( r ) B^{(r)}_{:,j} B : , j ( r ) Q j ( r ) Q^{(r)}_j Q j ( r )
Incident adaptation from a reference SCM. To estimate B ( 1 ) B^{(1)} B ( 1 ) n 1 n_1 n 1 Q Q Q X ( r ) X^{(r)} X ( r ) b ~ k \widetilde b_k b k k k k
δ ^ k ∈ argmin δ { 1 2 n 1 ∥ y k − X : , < k ( 1 ) δ ∥ 2 2 + λ 1 ∥ δ ∥ 1 } , \widehat{\delta}_k \in\operatorname*{argmin}_{\delta}\left\{\frac{1}{2n_1}\left\|y_k-X^{(1)}_{:,<k}\delta\right\|_2^2+\lambda_1\|\delta\|_1\right\}, δ k ∈ δ argmin { 2 n 1 1 y k − X : , < k ( 1 ) δ 2 2 + λ 1 ∥ δ ∥ 1 } , y k = X : , k ( 1 ) − X : , < k ( 1 ) b ~ k , b ^ k = b ~ k + δ ^ k . (16) y_k=X^{(1)}_{:,k}-X^{(1)}_{:,<k}\widetilde b_k,\qquad\widehat b_k=\widetilde b_k+\widehat\delta_k. \tag*{(16)} y k = X : , k ( 1 ) − X : , < k ( 1 ) b k , b k = b k + δ k . (16) For downstream plug-ins, set B ^ ( 0 ) : = B ~ \widehat B^{(0)}:=\widetilde B B ( 0 ) := B B ^ ( 1 ) : = B ^ δ \widehat B^{(1)}:=\widehat B^{\delta} B ( 1 ) := B δ Z k ( r ) = X : , < k ( r ) Z_k^{(r)}=X^{(r)}_{:,<k} Z k ( r ) = X : , < k ( r ) C ( S ; 3 ) = { θ : ∥ θ S c ∥ 1 ≤ 3 ∥ θ S ∥ 1 } C(S;3)=\{\theta:\|\theta_{S^c}\|_1\leq3\|\theta_S\|_1\} C ( S ; 3 ) = { θ : ∥ θ S c ∥ 1 ≤ 3∥ θ S ∥ 1 } S k ( 0 ) S_k^{(0)} S k ( 0 ) ∥ ⋅ ∥ o p \|\cdot\|_{\mathrm{op}} ∥ ⋅ ∥ op ∥ ⋅ ∥ F \|\cdot\|_{\mathrm{F}} ∥ ⋅ ∥ F ∥ ⋅ ∥ 1 , 1 \|\cdot\|_{1,1} ∥ ⋅ ∥ 1 , 1 ℓ 1 \ell_1 ℓ 1 O ( ⋅ ) O(\cdot) O ( ⋅ ) o ( ⋅ ) o(\cdot) o ( ⋅ ) O P ( ⋅ ) O_P(\cdot) O P ( ⋅ ) o P ( ⋅ ) o_P(\cdot) o P ( ⋅ ) σ r \sigma_r σ r
Assumption 4.1 (Reference mechanism fit (M0)). Reference rows are independent. With S k ( 0 ) = supp ( b k ( 0 ) ) S_k^{(0)}=\operatorname{supp}(b_k^{(0)}) S k ( 0 ) = supp ( b k ( 0 ) ) ∑ k ∣ S k ( 0 ) ∣ = s 0 \sum_k|S_k^{(0)}|=s_0 ∑ k ∣ S k ( 0 ) ∣ = s 0 κ 0 > 0 \kappa_0>0 κ 0 > 0 n 0 − 1 ∥ Z k ( 0 ) θ ∥ 2 2 ≥ κ 0 ∥ θ ∥ 2 2 n_0^{-1}\|Z_k^{(0)}\theta\|_2^2\geq\kappa_0\|\theta\|_2^2 n 0 − 1 ∥ Z k ( 0 ) θ ∥ 2 2 ≥ κ 0 ∥ θ ∥ 2 2 k ∈ [ d ] k\in[d] k ∈ [ d ] θ ∈ C ( S k ( 0 ) ; 3 ) \theta\in C(S_k^{(0)};3) θ ∈ C ( S k ( 0 ) ; 3 )
Assumption 4.2 (Incident transfer (M1)). Incident rows are mutually independent and independent of the reference sample. For common κ 1 > 0 \kappa_1>0 κ 1 > 0 κ 2 ≥ 0 \kappa_2\geq0 κ 2 ≥ 0 n 1 − 1 ∥ Z k ( 1 ) θ ∥ 2 2 ≥ κ 1 ∥ θ ∥ 2 2 − κ 2 ( log d / n 1 ) ∥ θ ∥ 1 2 n_1^{-1}\|Z_k^{(1)}\theta\|_2^2\geq\kappa_1\|\theta\|_2^2-\kappa_2(\log d/n_1)\|\theta\|_1^2 n 1 − 1 ∥ Z k ( 1 ) θ ∥ 2 2 ≥ κ 1 ∥ θ ∥ 2 2 − κ 2 ( log d / n 1 ) ∥ θ ∥ 1 2 k ∈ [ d ] k\in[d] k ∈ [ d ] θ ∈ R k − 1 \theta\in\mathbb{R}^{k-1} θ ∈ R k − 1
For the residual-moment plug-in below, we additionally use the following scaling condition.
Assumption 4.3 (Residual-moment scaling (R)). Uniformly over r ∈ { 0 , 1 } r\in\{0,1\} r ∈ { 0 , 1 } σ r = O ( 1 ) \sigma_r=O(1) σ r = O ( 1 ) ∥ A ( r ) ∥ o p = O ( 1 ) \|A^{(r)}\|_{\mathrm{op}}=O(1) ∥ A ( r ) ∥ op = O ( 1 )
Theorem 4.4 (Radius-based SCM Delta-Lasso transfer). Under the column-normalization and sub-Gaussian conditions above and Assumptions 4.1 and 4.2, let Δ = B ( 1 ) − B ( 0 ) \Delta=B^{(1)}-B^{(0)} Δ = B ( 1 ) − B ( 0 ) λ r ≳ σ r log d / n r \lambda_r\gtrsim\sigma_r\sqrt{\log d/n_r} λ r ≳ σ r log d / n r
R Δ : = ∥ Δ ∥ 1 , 1 + 12 λ 0 s 0 κ 0 . (17) R_\Delta:=\|\Delta\|_{1,1}+\frac{12\lambda_0s_0}{\kappa_0}. \tag*{(17)} R Δ := ∥Δ ∥ 1 , 1 + κ 0 12 λ 0 s 0 . (17) For any fixed c > 0 c>0 c > 0 κ 2 R Δ log d / n 1 ≤ c κ 1 \kappa_2R_\Delta\sqrt{\log d/n_1}\leq c\kappa_1 κ 2 R Δ log d / n 1 ≤ c κ 1
∥ B ^ δ − B ( 1 ) ∥ F 2 ≤ C R Δ σ 1 κ 1 log d n 1 , (18) \|\widehat B^\delta-B^{(1)}\|_{\mathrm{F}}^2\leq CR_\Delta\frac{\sigma_1}{\kappa_1}\sqrt{\frac{\log d}{n_1}}, \tag*{(18)} ∥ B δ − B ( 1 ) ∥ F 2 ≤ C R Δ κ 1 σ 1 n 1 log d , (18) where C C C n 0 n_0 n 0 n 1 n_1 n 1
When transfer helps. Let s 1 = ∥ B ( 1 ) ∥ 0 s_1 = \lVert B^{(1)}\rVert_0 s 1 = ∥ B ( 1 ) ∥ 0 κ > 0 \kappa> 0 κ > 0 O P ( ( σ 1 / κ ) s 1 log d / n 1 ) O_{\mathbb{P}}((\sigma_1/\kappa)\sqrt{s_1\log d/n_1}) O P (( σ 1 / κ ) s 1 log d / n 1 ) R Δ R_{\Delta} R Δ o ( ( σ 1 κ 1 / κ 2 ) s 1 log d / n 1 ) o((\sigma_1\kappa_1/\kappa^2)s_1\sqrt{\log d/n_1}) o (( σ 1 κ 1 / κ 2 ) s 1 log d / n 1 ) s 0 = O ( s 1 ) s_0 = O(s_1) s 0 = O ( s 1 ) n 0 / n 1 → ∞ n_0/n_1 \to\infty n 0 / n 1 → ∞ ∥ Δ ∥ 1 , 1 = o ( s 1 log d / n 1 ) \lVert\Delta\rVert_{1,1} = o(s_1\sqrt{\log d/n_1}) ∥ Δ ∥ 1 , 1 = o ( s 1 log d / n 1 ) n 0 ≫ n 1 n_0 \gg n_1 n 0 ≫ n 1 ∥ Δ ∥ 1 , 1 \lVert\Delta\rVert_{1,1} ∥ Δ ∥ 1 , 1
Held-out residual moments. For a regime-r r r ℓ \ell ℓ I r , ℓ \mathcal{I}_{r,\ell} I r , ℓ m r , ℓ = ∣ I r , ℓ ∣ m_{r,\ell} = \lvert\mathcal{I}_{r,\ell}\rvert m r , ℓ = ∣ I r , ℓ ∣ B ^ − ℓ ( r ) \widehat{B}^{(r)}_{-\ell} B − ℓ ( r ) ℓ \ell ℓ
e ^ i ( r ) = ( I − B ^ − ℓ ( r ) T ) x i ( r ) , i ∈ I r , \widehat{e}^{(r)}_i = (I-\widehat{B}^{(r)\mathsf{T}}_{-\ell})x^{(r)}_i,\qquad i\in\mathcal{I}_r, e i ( r ) = ( I − B − ℓ ( r ) T ) x i ( r ) , i ∈ I r , ( μ ^ j ( r ) , ω ^ j ( r ) ) = m r − 1 ∑ i ∈ I r ( e ^ i j ( r ) , ( e ^ i j ( r ) − μ ^ j ( r ) ) 2 ) . (19) (\widehat{\mu}^{(r)}_j,\widehat{\omega}^{(r)}_j)=m_r^{-1}\sum_{i\in\mathcal{I}_r}(\widehat{e}^{(r)}_{ij},(\widehat{e}^{(r)}_{ij}-\widehat{\mu}^{(r)}_j)^2). \tag*{(19)} ( μ j ( r ) , ω j ( r ) ) = m r − 1 i ∈ I r ∑ ( e ij ( r ) , ( e ij ( r ) − μ j ( r ) ) 2 ) . (19) Excluding I r , ℓ \mathcal{I}_{r,\ell} I r , ℓ B ^ − ℓ ( r ) \widehat{B}^{(r)}_{-\ell} B − ℓ ( r )
Plug-in attribution and ranking. For a node-local game, replace the endpoints in (10) to form B ^ U , μ ^ U , ω ^ U \widehat{B}^{U},\widehat{\mu}^{U},\widehat{\omega}^{U} B U , μ U , ω U
T ^ ( U ) = T ( B ^ U , μ ^ U , ω ^ U ) . (20) \widehat{T}(U)=T(\widehat{B}^{U},\widehat{\mu}^{U},\widehat{\omega}^{U}). \tag*{(20)} T ( U ) = T ( B U , μ U , ω U ) . (20) Let Φ ^ = ( Φ ^ u ) u ∈ U \widehat{\Phi}=(\widehat{\Phi}_u)_{u\in\mathcal{U}} Φ = ( Φ u ) u ∈ U (2) to U ↦ T ^ ( U ) U\mapsto\widehat{T}(U) U ↦ T ( U ) Θ = ( θ ( 0 ) , θ ( 1 ) ) \Theta=(\theta^{(0)},\theta^{(1)}) Θ = ( θ ( 0 ) , θ ( 1 ) ) θ ( r ) = ( B ( r ) , μ ( r ) , ω ( r ) ) \theta^{(r)}=(B^{(r)},\mu^{(r)},\omega^{(r)}) θ ( r ) = ( B ( r ) , μ ( r ) , ω ( r ) ) Θ ^ \widehat{\Theta} Θ Φ u ( Θ ) \Phi_u(\Theta) Φ u ( Θ ) Θ s = Θ + s ( Θ ^ − Θ ) \Theta_s=\Theta+s(\widehat{\Theta}-\Theta) Θ s = Θ + s ( Θ − Θ ) s ∈ [ 0 , 1 ] s\in[0,1] s ∈ [ 0 , 1 ] B \mathcal{B} B Π B \Pi_{\mathcal{B}} Π B η B ( r ) = ∥ B ^ ( r ) − B ( r ) ∥ F \eta^{(r)}_B=\lVert\widehat{B}^{(r)}-B^{(r)}\rVert_{\mathrm{F}} η B ( r ) = ∥ B ( r ) − B ( r ) ∥ F η μ ( r ) = ∥ μ ^ ( r ) − μ ( r ) ∥ ∞ \eta^{(r)}_\mu=\lVert\widehat{\mu}^{(r)}-\mu^{(r)}\rVert_{\infty} η μ ( r ) = ∥ μ ( r ) − μ ( r ) ∥ ∞ η ω ( r ) = ∥ ω ^ ( r ) − ω ( r ) ∥ ∞ \eta^{(r)}_\omega=\lVert\widehat{\omega}^{(r)}-\omega^{(r)}\rVert_{\infty} η ω ( r ) = ∥ ω ( r ) − ω ( r ) ∥ ∞ u u u
L u , B ( r ) : = sup s ∈ [ 0 , 1 ] ∥ Π B ∇ B ( r ) Φ u ( Θ s ) ∥ F , L^{(r)}_{u,B}:=\sup_{s\in[0,1]}\lVert\Pi_{\mathcal{B}}\nabla_{B^{(r)}}\Phi_u(\Theta_s)\rVert_{\mathrm{F}}, L u , B ( r ) := s ∈ [ 0 , 1 ] sup ∥ Π B ∇ B ( r ) Φ u ( Θ s ) ∥ F , L u , μ ( r ) : = sup s ∈ [ 0 , 1 ] ∥ ∇ μ ( r ) Φ u ( Θ s ) ∥ 1 , L^{(r)}_{u,\mu}:=\sup_{s\in[0,1]}\lVert\nabla_{\mu^{(r)}}\Phi_u(\Theta_s)\rVert_1, L u , μ ( r ) := s ∈ [ 0 , 1 ] sup ∥ ∇ μ ( r ) Φ u ( Θ s ) ∥ 1 , L u , ω ( r ) : = sup s ∈ [ 0 , 1 ] ∥ ∇ ω ( r ) Φ u ( Θ s ) ∥ 1 . (21) L^{(r)}_{u,\omega}:=\sup_{s\in[0,1]}\lVert\nabla_{\omega^{(r)}}\Phi_u(\Theta_s)\rVert_1. \tag*{(21)} L u , ω ( r ) := s ∈ [ 0 , 1 ] sup ∥ ∇ ω ( r ) Φ u ( Θ s ) ∥ 1 . (21) Uniform error rate. For deterministic ρ r \rho_r ρ r η B ( r ) = O P ( ρ r ) \eta^{(r)}_B=O_{\mathbb{P}}(\rho_r) η B ( r ) = O P ( ρ r ) ν r = log ( 2 d ) / m r \nu_r=\sqrt{\log(2d)/m_r} ν r = log ( 2 d ) / m r M r = ∥ E x ( r ) ∥ 2 M_r=\lVert\mathbb{E}x^{(r)}\rVert_2 M r = ∥ E x ( r ) ∥ 2 q μ , r = ν r + ( M r + ν r ) ρ r q_{\mu,r}=\nu_r+(M_r+\nu_r)\rho_r q μ , r = ν r + ( M r + ν r ) ρ r q ω , r = ν r + ρ r 2 q_{\omega,r}=\nu_r+\rho_r^2 q ω , r = ν r + ρ r 2 a ∗ = max U ⊆ U ∥ ( I − B U ) − 1 ∥ o p a_*=\max_{U\subseteq\mathcal{U}}\lVert(I-B^U)^{-1}\rVert_{\mathrm{op}} a ∗ = max U ⊆ U ∥( I − B U ) − 1 ∥ op W ∗ = ∥ W ∥ o p W_*=\lVert W\rVert_{\mathrm{op}} W ∗ = ∥ W ∥ op z ∗ = ∥ z ∥ 2 z_*=\lVert z\rVert_2 z ∗ = ∥ z ∥ 2 d d d
m ∗ = 1 + ∑ r = 0 1 { ∥ μ ( r ) ∥ 2 + d q μ , r } , m_* = 1+\sum_{r=0}^{1}\{\lVert\mu^{(r)}\rVert_2+\sqrt{d}\,q_{\mu,r}\}, m ∗ = 1 + r = 0 ∑ 1 {∥ μ ( r ) ∥ 2 + d q μ , r } , ω ∗ = 1 + max r { ∥ ω ( r ) ∥ ∞ + q ω , r } , \omega_* = 1+\max_r\{\lVert\omega^{(r)}\rVert_\infty+q_{\omega,r}\}, ω ∗ = 1 + r max {∥ ω ( r ) ∥ ∞ + q ω , r } , H B = 16 a ∗ 3 ( d ω ∗ + m ∗ 2 ) W ∗ + 4 a ∗ 2 m ∗ z ∗ , H_B=16a_*^3(\sqrt{d}\,\omega_*+m_*^2)W_*+4a_*^2m_*z_*, H B = 16 a ∗ 3 ( d ω ∗ + m ∗ 2 ) W ∗ + 4 a ∗ 2 m ∗ z ∗ , H μ = d { 8 a ∗ 2 W ∗ m ∗ + 2 a ∗ z ∗ } , H ω = 4 d a ∗ 2 W ∗ . H_\mu=\sqrt{d}\{8a_*^2W_*m_*+2a_*z_*\},\qquad H_\omega=4da_*^2W_*. H μ = d { 8 a ∗ 2 W ∗ m ∗ + 2 a ∗ z ∗ } , H ω = 4 d a ∗ 2 W ∗ . Theorem 4.5 (Plug-in attribution error). For any fitted endpoints with slopes in B \mathcal{B} B
∣ Φ ^ u − Φ u ∣ ≤ ∑ r = 0 1 { L u , B ( r ) η B ( r ) + L u , μ ( r ) η μ ( r ) + L u , ω ( r ) η ω ( r ) } . (22) \left|\widehat{\Phi}_u-\Phi_u\right| \le\sum_{r=0}^{1}\left\{L_{u,B}^{(r)}\eta_{B}^{(r)}+L_{u,\mu}^{(r)}\eta_{\mu}^{(r)}+L_{u,\omega}^{(r)}\eta_{\omega}^{(r)}\right\}. \tag*{(22)} Φ u − Φ u ≤ r = 0 ∑ 1 { L u , B ( r ) η B ( r ) + L u , μ ( r ) η μ ( r ) + L u , ω ( r ) η ω ( r ) } . (22) For the out-of-sample construction above, under the held-out-sample independence and sub-Gaussian conditions above and Assumption 4.3, if max r ν r = o ( 1 ) \max_r \nu_r=o(1) max r ν r = o ( 1 ) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 ) a_*(\rho_0+\rho_1)=o(1) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 ) ∥ Φ ^ − Φ ∥ ∞ = O P ( R d ) \lVert\widehat{\Phi}-\Phi\rVert_{\infty}=O_{\mathbb{P}}(R_d) ∥ Φ − Φ ∥ ∞ = O P ( R d )
R d : = ∑ r = 0 1 { H B β ρ r + H μ q μ , r + H ω q ω , r } . (23) R_d:=\sum_{r=0}^{1}\left\{H_B\beta\rho_r+H_\mu q_{\mu,r}+H_\omega q_{\omega,r}\right\}. \tag*{(23)} R d := r = 0 ∑ 1 { H B β ρ r + H μ q μ , r + H ω q ω , r } . (23) For a fixed number of folds, averaging their foldwise bounds preserves the corresponding order.
The theorem converts endpoint rates into attribution accuracy. For the estimators above, take ρ 0 = ( σ 0 / κ 0 ) s 0 log d / n 0 \rho_0=(\sigma_0/\kappa_0)\sqrt{s_0\log d/n_0} ρ 0 = ( σ 0 / κ 0 ) s 0 log d / n 0 ρ 1 = [ R Δ ( σ 1 / κ 1 ) log d / n 1 ] 1 / 2 \rho_1=[R_{\Delta}(\sigma_1/\kappa_1)\sqrt{\log d/n_1}]^{1/2} ρ 1 = [ R Δ ( σ 1 / κ 1 ) log d / n 1 ] 1/2 R d = o ( 1 ) R_d=o(1) R d = o ( 1 ) a ∗ a_* a ∗ a ∗ , W ∗ = O ( 1 ) a_*,W_*=O(1) a ∗ , W ∗ = O ( 1 ) H ω q ω , r H_\omega q_{\omega,r} H ω q ω , r O P ( d log d / m r + d ρ r 2 ) O_{\mathbb{P}}(d\sqrt{\log d/m_r}+d\rho_r^2) O P ( d log d / m r + d ρ r 2 ) m r ≫ d 2 log d m_r\gg d^2\log d m r ≫ d 2 log d
Corollary 4.6 (Top-k k k Under the conditions of Theorem 4.5, let S k S_k S k k k k g k : = min u ∈ S k ∣ Φ u ∣ − max v ∉ S k ∣ Φ v ∣ g_k:=\min_{u\in S_k}|\Phi_u|-\max_{v\notin S_k}|\Phi_v| g k := min u ∈ S k ∣ Φ u ∣ − max v ∈ / S k ∣ Φ v ∣ g k > 0 g_k>0 g k > 0 R d = o ( g k ) R_d=o(g_k) R d = o ( g k ) k k k u ∈ U u\in\mathcal{U} u ∈ U R d = o ( ∣ Φ u ∣ ) R_d=o(|\Phi_u|) R d = o ( ∣ Φ u ∣ )
Thus the leading function/law attributions are ranked correctly once their population gap dominates the uniform estimation error.
Experiments
We test exact computation, scarce-incident function/law recovery, and model-relative measured-intervention diagnosis. Supplement Secs. D.1–D.8 give protocols, uncertainty summaries, stress tests, and reproducibility details.
Exact computation. All methods query the same oracle for true-endpoint coalition values. We use DoWhy-GCM for the DW comparators [3 ] ; Figure 1 compares runtime and relative ℓ 2 \ell_2 ℓ 2 4 ^4 4 100 × 100\times 100 × K ! K! K ! 2 K 2^K 2 K K = 20 K=20 K = 20 4 ^4 4 K = 1000 K=1000 K = 1000 10 3 10^3 1 0 3
Scarce-incident attribution benchmark. We use 30 ordered Gaussian d = 100 d=100 d = 100 4 ^4 4 Figure 2 reports F F F Q Q Q F + Q F+Q F + Q ∣ Δ B ∣ |\Delta B| ∣Δ B ∣
Figure 2. Scarce-incident recovery over 30 ordered d = 100 d = 100 d = 100 ± \pm ± F F F Q Q Q F + Q F + Q F + Q F + Q F + Q F + Q 4 ^{4} 4 F F F
The advantage is clearest for function recovery: at the largest incident size, recall is 0.687, versus 0.460 for ∣ Δ B ∣ |\Delta B| ∣Δ B ∣ 4 ^4 4 F + Q F+Q F + Q Q Q Q F + Q F+Q F + Q
is comparable. Supplement Sec. D.4 gives the protocol and a broader benchmark including fixed-F F F
Measured interventions. Across 16 22-channel Causal Chamber Light-Tunnel interventions [9 ] , we test whether a measured target-law component separates from propagated fitted-response components of a reference-fitted PCA T 2 + S P E T^2 + \mathrm{SPE} T 2 + SPE F F F
S 4 S^4 S 4 F F F Q Q Q F F F Figure 3 ). Within the fitted ordered linear SCM, this allocation separates the target-law component from downstream relinearization. Supplement Sec. D.6 gives the full protocol and split sensitivity; Sec. D.7 extends the stress test to two other datasets.
Figure 3. Causal Chamber all-channel score. (a) Target rank (1 best; worst of five splits; m/s: severity). (b) Target-Q Q Q i r 3 ir_3 i r 3 F F F
Conclusion and Limitations
S 4 S^4 S 4
Appendix organization. Appendix A gives the framework scope, closest-work reductions, and hierarchy-preserving Owen extension. Appendix B gives pseudocode and implementation details. Appendix C proves the general identities, the linear/quadratic computation results, and the plug-in estimation guarantees used in Sections 2–4. Appendix D records the experimental protocols and boundary diagnostics needed to reproduce and interpret the main claims.
Extensions and Closest-Work Reductions
This appendix develops two aspects compressed in the main paper. First, it states the scope of the source–target–player–score construction and maps closely related methods into that common notation. It then gives an Owen-value extension that preserves paired-node hierarchy while retaining separate function/law allocations.
Framework scope and closest-work reductions
The Shapley definition in Eq. (2) applies to any finite scalar functional T T T T ( M ) = E M [ τ ( X ) ] T(M)=\mathbb{E}_{M}[\tau(X)] T ( M ) = E M [ τ ( X )]
A function/law-resolved report is conditional on aligned variables and exogenous coordinates across its endpoints; unrestricted observational laws do not identify a unique physical function/noise decomposition. The endpoints may be supplied or estimated under additional assumptions. In the ordered linear specialization, a correct common order (or parent sets), nonsingular parent covariances, and independent disturbances identify the slopes and residual laws and justify product-law hybrids. There is no separate structural intercept: F j F_j F j Q j Q_j Q j
Any attribution report is determined by four choices: source M ( 0 ) M^{(0)} M ( 0 ) M ( 1 ) M^{(1)} M ( 1 ) U \mathcal{U} U T T T Table 1 translates every framework into this paper’s source–target–player–score notation. For attribution methods, “source” and “target” mean the empty- and grand-coalition endpoints M ∅ M^{\varnothing} M ∅ M U M^{\mathcal{U}} M U T ( M U ) − T ( M ∅ ) T(M^{\mathcal{U}})-T(M^{\varnothing}) T ( M U ) − T ( M ∅ )
Table 1. Three representative causal score-attribution questions in the source–target–player–score framework.
Reading the map. Population CCA/MRCCA ([4 , 25 ] ) is the closest population-attribution framework: it becomes the same coalition game when each function/law pair is treated as one node player. Fixed-SCM outlier RCA ([6 ] ) shares a factual-noise replacement orientation, but attributes a different event-level tail functional. Unit-level change attribution ([5 ] ) instead switches observed values, and do-Shapley ([34 ] ) switches variable interventions within one SCM. These distinctions explain why the first two relations merit explicit mappings below, while the remaining rows are adjacent causal questions rather than reductions of the population alarm-score estimand.
Paired-node population reduction. Suppose the endpoints are Markovian, have aligned variables, and share a common acyclic parent universe. If each pair { F j , Q j } \{F_j,Q_j\} { F j , Q j } A ⊆ R A\subseteq\mathbb{R} A ⊆ R
C j ( r ) ( A ∣ x < j ) : = Q j ( r ) ( { e j : F j ( r ) ( x < j , e j ) ∈ A } ) . C_j^{(r)}(A\mid x_{<j}) := Q_j^{(r)}\left(\left\{e_j:F_j^{(r)}(x_{<j},e_j)\in A\right\}\right). C j ( r ) ( A ∣ x < j ) := Q j ( r ) ( { e j : F j ( r ) ( x < j , e j ) ∈ A } ) . Specifically,
P U ( d x ) = ∏ j = 1 d C j 1 { { F j , Q j } ∈ U } ( d x j ∣ x < j ) . P^{\mathcal{U}}(dx)=\prod_{j=1}^{d}C_j^{\mathbf{1}_{\{\{F_j,Q_j\}\in\mathcal{U}\}}}(dx_j\mid x_{<j}). P U ( d x ) = j = 1 ∏ d C j 1 {{ F j , Q j } ∈ U } ( d x j ∣ x < j ) . Table 1 (PDF p. 13). Closest frameworks translated into this paper’s notation. CCA/MRCCA neither identifies nor separately switches the function/law representation of an induced node conditional. For outlier RCA, A x = { τ ( X ) ≥ τ ( x ) } A_x=\{\tau(X)\ge\tau(x)\} A x = { τ ( X ) ≥ τ ( x )} S M ( 0 ) T ( x ) = − log P M ( 0 ) ( A x ) S_{M^{(0)}}^T(x)=-\log P_{M^{(0)}}(A_x) S M ( 0 ) T ( x ) = − log P M ( 0 ) ( A x ) δ x ( r ) = ⨂ j = 1 d δ x j ( r ) \delta_{x}^{(r)}=\bigotimes_{j=1}^{d}\delta_{x_j}^{(r)} δ x ( r ) = ⨂ j = 1 d δ x j ( r ) X j X_j X j
Thus, for T ( M ) = E M [ τ ( X ) ] T(M)=\mathbb{E}_M[\tau(X)] T ( M ) = E M [ τ ( X )] 4 {}^{4} 4 τ ( X ) \tau(X) τ ( X ) [4 , 25 ] ), and so has the same Shapley values. This reduction does not imply Ψ j = Φ F j + Φ Q j \Psi_j=\Phi_{F_j}+\Phi_{Q_j} Ψ j = Φ F j + Φ Q j
Fixed-SCM event mapping. The native fixed-SCM outlier game (Budhathoki et al. [6 ] ) indexes a coalition I I I δ e ( 0 ) ( x ) \delta_{e^{(0)}(x)} δ e ( 0 ) ( x ) Q ( 0 ) Q^{(0)} Q ( 0 ) S = [ d ] ∖ I S=[d]\setminus I S = [ d ] ∖ I
w ( S ) = log P M ( F ( 0 ) , δ e S ( 0 ) ( x ) ⊗ Q − S ( 0 ) ) ( A x ) . w(S)=\log P_{M(F^{(0)},\delta_{e_S^{(0)}(x)}\otimes Q_{-S}^{(0)})}(A_x). w ( S ) = log P M ( F ( 0 ) , δ e S ( 0 ) ( x ) ⊗ Q − S ( 0 ) ) ( A x ) . If v ( I ) v(I) v ( I ) w ( S ) = − v ( [ d ] ∖ S ) w(S)=-v([d]\setminus S) w ( S ) = − v ([ d ] ∖ S ) Table 1 correctly writes its source and target as Q ( 0 ) Q^{(0)} Q ( 0 ) δ e ( 0 ) ( x ) \delta_{e^{(0)}(x)} δ e ( 0 ) ( x ) T x T_x T x S M ( 0 ) T ( x ) S_{M^{(0)}}^T(x) S M ( 0 ) T ( x )
FANS and nonlinear function/law dissection. FANS shows that function/noise dissection need not be confined to additive models [28 ] . It considers the following nonlinear, non-additive class, with a differentiable noise inverse for each fixed h h h
F j ( r ) ( p a j , e j ) = g j ( h j ( r ) ( p a j ) , e j ) , F_j^{(r)}(pa_j,e_j)=g_j\left(h_j^{(r)}(pa_j),e_j\right), F j ( r ) ( p a j , e j ) = g j ( h j ( r ) ( p a j ) , e j ) , e j ↦ g j ( h , e j ) is a differentiable bijection. e_j\mapsto g_j(h,e_j)\text{ is a differentiable bijection.} e j ↦ g j ( h , e j ) is a differentiable bijection. where the link g j g_j g j g j ( h j ( 1 ) ( p a j ) , e j ( 0 ) ) g_j\left(h_j^{(1)}(pa_j),e_j^{(0)}\right) g j ( h j ( 1 ) ( p a j ) , e j ( 0 ) ) p a j pa_j p a j g j ( h j ( 0 ) ( p a j ) , ψ ( e j ( 0 ) ) ) g_j\left(h_j^{(0)}(pa_j),\psi(e_j^{(0)})\right) g j ( h j ( 0 ) ( p a j ) , ψ ( e j ( 0 ) ) ) ψ \psi ψ
Hierarchy-preserving Owen attribution
Let N = { F j , Q j : j ∈ [ d ] } N=\{F_j,Q_j:j\in[d]\} N = { F j , Q j : j ∈ [ d ]} v ( A ) = T ( M A ) v(A)=T(M^A) v ( A ) = T ( M A ) N N N P = { P 1 , … , P d } \mathcal{P}=\{P_1,\ldots,P_d\} P = { P 1 , … , P d } P j = { F j , Q j } P_j=\{F_j,Q_j\} P j = { F j , Q j } R ⊆ [ d ] R\subseteq[d] R ⊆ [ d ] P ( R ) = ⋃ k ∈ R P k \mathcal{P}(R)=\bigcup_{k\in R}P_k P ( R ) = ⋃ k ∈ R P k
w d ( R ) = ∣ R ∣ ! ( d − ∣ R ∣ − 1 ) ! d ! . w_d(R)=\frac{|R|!(d-|R|-1)!}{d!}. w d ( R ) = d ! ∣ R ∣ ! ( d − ∣ R ∣ − 1 )! . The Owen value [22 ] averages marginal contributions over orders in which the two players of each block remain contiguous. Since every block here has size two, its allocations have the explicit form
Ω F j = 1 2 ∑ R ⊆ [ d ] ∖ { j } w d ( R ) [ v ( P ( R ) ∪ { F j } ) − v ( P ( R ) ) + v ( P ( R ) ∪ P j ) − v ( P ( R ) ∪ { Q j } ) ] . (24) \Omega_{F_j}=\frac{1}{2}\sum_{R\subseteq[d]\setminus\{j\}}w_d(R)\left[v\left(\mathcal{P}(R)\cup\{F_j\}\right)-v\left(\mathcal{P}(R)\right)+v\left(\mathcal{P}(R)\cup P_j\right)-v\left(\mathcal{P}(R)\cup\{Q_j\}\right)\right]. \tag*{(24)} Ω F j = 2 1 R ⊆ [ d ] ∖ { j } ∑ w d ( R ) [ v ( P ( R ) ∪ { F j } ) − v ( P ( R ) ) + v ( P ( R ) ∪ P j ) − v ( P ( R ) ∪ { Q j } ) ] . (24) Ω Q j = 1 2 ∑ R ⊆ [ d ] ∖ { j } w d ( R ) [ v ( P ( R ) ∪ { Q j } ) − v ( P ( R ) ) + v ( P ( R ) ∪ P j ) − v ( P ( R ) ∪ { F j } ) ] . (25) \Omega_{Q_j}=\frac{1}{2}\sum_{R\subseteq[d]\setminus\{j\}}w_d(R)\left[v\left(\mathcal{P}(R)\cup\{Q_j\}\right)-v\left(\mathcal{P}(R)\right)+v\left(\mathcal{P}(R)\cup P_j\right)-v\left(\mathcal{P}(R)\cup\{F_j\}\right)\right]. \tag*{(25)} Ω Q j = 2 1 R ⊆ [ d ] ∖ { j } ∑ w d ( R ) [ v ( P ( R ) ∪ { Q j } ) − v ( P ( R ) ) + v ( P ( R ) ∪ P j ) − v ( P ( R ) ∪ { F j } ) ] . (25) Let v P ( R ) = v ( P ( R ) ) v^{\mathcal{P}}(R)=v\left(\mathcal{P}(R)\right) v P ( R ) = v ( P ( R ) ) Ψ j \Psi_j Ψ j
Ω F j + Ω Q j = Ψ j , ∑ j = 1 d ( Ω F j + Ω Q j ) = v ( N ) − v ( ∅ ) . (26) \Omega_{F_j}+\Omega_{Q_j}=\Psi_j,\qquad\sum_{j=1}^{d}\left(\Omega_{F_j}+\Omega_{Q_j}\right)=v(N)-v(\varnothing). \tag*{(26)} Ω F j + Ω Q j = Ψ j , j = 1 ∑ d ( Ω F j + Ω Q j ) = v ( N ) − v ( ∅ ) . (26) For exact computation in the ordered-linear/quadratic-score specialization, give node j j j r s ∈ { 00 , 10 , 01 , 11 } rs\in\{00,10,01,11\} r s ∈ { 00 , 10 , 01 , 11 } H ( π ) H(\boldsymbol{\pi}) H ( π ) π 00 , j = 1 − t \pi_{00,j}=1-t π 00 , j = 1 − t π 11 , j = t \pi_{11,j}=t π 11 , j = t π 10 , j = π 01 , j = 0 \pi_{10,j}=\pi_{01,j}=0 π 10 , j = π 01 , j = 0 g j r s ( t ) = ∂ H / ∂ π r s , j g_j^{rs}(t)=\partial H/\partial\pi_{rs,j} g j r s ( t ) = ∂ H / ∂ π r s , j D p a i r , j = g j 11 − g j 00 D_{\mathrm{pair},j}=g_j^{11}-g_j^{00} D pair , j = g j 11 − g j 00 D c r o s s , j = g j 10 − g j 01 D_{\mathrm{cross},j}=g_j^{10}-g_j^{01} D cross , j = g j 10 − g j 01 d − 1 d-1 d − 1
Ω F j = 1 2 ∫ 0 1 { D p a i r , j ( t ) + D c r o s s , j ( t ) } d t , (27) \Omega_{F_j}=\frac{1}{2}\int_0^1\left\{D_{\mathrm{pair},j}(t)+D_{\mathrm{cross},j}(t)\right\}\,dt, \tag*{(27)} Ω F j = 2 1 ∫ 0 1 { D pair , j ( t ) + D cross , j ( t ) } d t , (27) Ω Q j = 1 2 ∫ 0 1 { D p a i r , j ( t ) − D c r o s s , j ( t ) } d t . \Omega_{Q_j}=\frac{1}{2}\int_0^1\left\{D_{\mathrm{pair},j}(t)-D_{\mathrm{cross},j}(t)\right\}\,dt. Ω Q j = 2 1 ∫ 0 1 { D pair , j ( t ) − D cross , j ( t ) } d t . Each difference is a Bernstein polynomial of degree at most d − 1 d-1 d − 1 ⌈ d / 2 ⌉ \lceil d/2\rceil ⌈ d /2 ⌉ O ( d 3 ) O(d^3) O ( d 3 ) O ( d 2 ) O(d^2) O ( d 2 ) O ( d 4 ) O(d^4) O ( d 4 ) O ( d 2 ) O(d^2) O ( d 2 )
Linearity in coalition values also preserves event–population coherence for these prespecified coarse F j / Q j F_j/Q_j F j / Q j
Algorithm 1 Exact second-order attribution for the resolved function/law game
Input: Endpoint moments { b j ( r ) , μ j ( r ) , ω j ( r ) } j , r \{b_j^{(r)},\mu_j^{(r)},\omega_j^{(r)}\}_{j,r} { b j ( r ) , μ j ( r ) , ω j ( r ) } j , r W W W z z z
Output: Resolved attribution vector Φ \Phi Φ
Gauss–Legendre quadrature and endpoint setup: Set n q = d n_q=d n q = d G ∈ R n q × n q G\in\mathbb{R}^{n_q\times n_q} G ∈ R n q × n q G k k = 0 G_{kk}=0 G k k = 0 k ∈ [ n q ] k\in[n_q] k ∈ [ n q ] G k , k + 1 = G k + 1 , k = k / 4 k 2 − 1 G_{k,k+1}=G_{k+1,k}=k/\sqrt{4k^2-1} G k , k + 1 = G k + 1 , k = k / 4 k 2 − 1 k ∈ [ n q − 1 ] k\in[n_q-1] k ∈ [ n q − 1 ]
Compute G = V Diag ( ξ 1 , … , ξ n q ) V ⊤ G=V\operatorname{Diag}(\xi_1,\ldots,\xi_{n_q})V^\top G = V Diag ( ξ 1 , … , ξ n q ) V ⊤ t ℓ ← ( 1 + ξ ℓ ) / 2 t_\ell\leftarrow(1+\xi_\ell)/2 t ℓ ← ( 1 + ξ ℓ ) /2 w ℓ ← V 1 ℓ 2 w_\ell\leftarrow V_{1\ell}^2 w ℓ ← V 1 ℓ 2
Initialize Φ ← 0 \Phi\leftarrow0 Φ ← 0
for j = 1 , … , d j=1,\ldots,d j = 1 , … , d
Δ b j ← b j ( 1 ) − b j ( 0 ) \Delta b_j\leftarrow b_j^{(1)}-b_j^{(0)} Δ b j ← b j ( 1 ) − b j ( 0 ) Δ μ j ← μ j ( 1 ) − μ j ( 0 ) \Delta\mu_j\leftarrow\mu_j^{(1)}-\mu_j^{(0)} Δ μ j ← μ j ( 1 ) − μ j ( 0 )
Δ Γ j ← b j ( 1 ) ( b j ( 1 ) ) ⊤ − b j ( 0 ) ( b j ( 0 ) ) ⊤ \Delta\Gamma_j\leftarrow b_j^{(1)}(b_j^{(1)})^\top-b_j^{(0)}(b_j^{(0)})^\top Δ Γ j ← b j ( 1 ) ( b j ( 1 ) ) ⊤ − b j ( 0 ) ( b j ( 0 ) ) ⊤ Δ q j ← ( ( μ j ( 1 ) ) 2 + ω j ( 1 ) ) − ( ( μ j ( 0 ) ) 2 − ω j ( 0 ) ) \Delta q_j\leftarrow((\mu_j^{(1)})^2+\omega_j^{(1)})-((\mu_j^{(0)})^2-\omega_j^{(0)}) Δ q j ← (( μ j ( 1 ) ) 2 + ω j ( 1 ) ) − (( μ j ( 0 ) ) 2 − ω j ( 0 ) )
end for
Forward–reverse evaluation at quadrature nodes:
for ℓ = 1 , … , n q \ell=1,\ldots,n_q ℓ = 1 , … , n q
t ← t ℓ t\leftarrow t_\ell t ← t ℓ p ← 1 p\leftarrow\mathbf{1} p ← 1
for j = 1 , … , d j=1,\ldots,d j = 1 , … , d
b ˉ j ← ( 1 − t ) b j ( 0 ) + t b j ( 1 ) \bar b_j\leftarrow(1-t)b_j^{(0)}+tb_j^{(1)} b ˉ j ← ( 1 − t ) b j ( 0 ) + t b j ( 1 ) μ ˉ j ← ( 1 − t ) μ j ( 0 ) + t μ j ( 1 ) \bar\mu_j\leftarrow(1-t)\mu_j^{(0)}+t\mu_j^{(1)} μ ˉ j ← ( 1 − t ) μ j ( 0 ) + t μ j ( 1 )
Γ j ← ( 1 − t ) b j ( 0 ) ( b j ( 0 ) ) ⊤ + t b j ( 1 ) ( b j ( 1 ) ) ⊤ \Gamma_j\leftarrow(1-t)b_j^{(0)}(b_j^{(0)})^\top+tb_j^{(1)}(b_j^{(1)})^\top Γ j ← ( 1 − t ) b j ( 0 ) ( b j ( 0 ) ) ⊤ + t b j ( 1 ) ( b j ( 1 ) ) ⊤
q j ← ( 1 − t ) { ( μ j ( 0 ) ) 2 + ω j ( 0 ) } + t { ( μ j ( 1 ) ) 2 + ω j ( 1 ) } q_j\leftarrow(1-t)\{(\mu_j^{(0)})^2+\omega_j^{(0)}\}+t\{(\mu_j^{(1)})^2+\omega_j^{(1)}\} q j ← ( 1 − t ) {( μ j ( 0 ) ) 2 + ω j ( 0 ) } + t {( μ j ( 1 ) ) 2 + ω j ( 1 ) } c j ← b ˉ j μ ˉ j c_j\leftarrow\bar b_j\bar\mu_j c j ← b ˉ j μ ˉ j
end for
Run the forward causal-moment recursion in (13) to obtain x ˉ ( p ) \bar x(p) x ˉ ( p ) S ( p ) S(p) S ( p )
Initialize the shared reverse sweep with sensitivities a ← z a\leftarrow z a ← z H ← W H\leftarrow W H ← W
for j = d , d − 1 , … , 1 j=d,d-1,\ldots,1 j = d , d − 1 , … , 1
α j ← a j \alpha_j\leftarrow a_j α j ← a j ρ j ← H < j , j + H j , < j ⊤ \rho_j\leftarrow H_{<j,j}+H^\top_{j,<j} ρ j ← H < j , j + H j , < j ⊤ γ j ← H j j \gamma_j\leftarrow H_{jj} γ j ← H j j
g F j ← ( α j x ˉ < j + S < j , < j ρ j + 2 γ j μ ˉ j x ˉ < j ) ⊤ Δ b j + γ j ⟨ Δ Γ j , S < j , < j ⟩ F g_{Fj}\leftarrow(\alpha_j\bar x_{<j}+S_{<j,<j}\rho_j+2\gamma_j\bar\mu_j\bar x_{<j})^\top\Delta b_j+\gamma_j\langle\Delta\Gamma_j,S_{<j,<j}\rangle_F g F j ← ( α j x ˉ < j + S < j , < j ρ j + 2 γ j μ ˉ j x ˉ < j ) ⊤ Δ b j + γ j ⟨ Δ Γ j , S < j , < j ⟩ F
g Q j ← ( α j + ρ j ⊤ x ˉ < j + 2 γ j μ ˉ j x ˉ < j ) Δ μ j + γ j Δ q j g_{Qj}\leftarrow(\alpha_j+\rho_j^\top\bar x_{<j}+2\gamma_j\bar\mu_j\bar x_{<j})\Delta\mu_j+\gamma_j\Delta q_j g Q j ← ( α j + ρ j ⊤ x ˉ < j + 2 γ j μ ˉ j x ˉ < j ) Δ μ j + γ j Δ q j
a < j + = α j b ˉ j + ρ j μ ˉ j + 2 γ j μ ˉ j b ˉ j a_{<j}+=\alpha_j\bar b_j+\rho_j\bar\mu_j+2\gamma_j\bar\mu_j\bar b_j a < j + = α j b ˉ j + ρ j μ ˉ j + 2 γ j μ ˉ j b ˉ j
H < j , < j + = ρ j b ˉ j ⊤ + γ j Γ j H_{<j,<j}+=\rho_j\bar b_j^\top+\gamma_j\Gamma_j H < j , < j + = ρ j b ˉ j ⊤ + γ j Γ j
end for
Accumulate Φ + = w ℓ g \Phi+=w_\ell g Φ + = w ℓ g
end for
return Φ \Phi Φ
Algorithms and Implementation Details
This section turns the computation and estimation results into two executable procedures. Algorithm 1 takes two SCM endpoint moment sets and a fixed quadratic score, propagates the interpolated moments in causal order, and obtains every resolved Shapley value through one shared reverse sweep at each Gauss–Legendre point. The reverse sweep differentiates the moment program; it neither reverses causal edges nor defines a reverse SCM. Supplement Sec. C.2 gives the corresponding recursions and proof. Each quadrature point costs O ( d 3 ) O(d^3) O ( d 3 ) O ( d 2 ) O(d^2) O ( d 2 ) K = 2 d K=2d K = 2 d d d d
Algorithm 2 implements the scarce-incident estimator. It fits the reference endpoint once, learns incident deviations only on training folds, evaluates residual moments on held-out rows, and sends each fold-matched endpoint pair to the exact routine before averaging. Supplement Secs. C.4–C.7 give the estimation conditions and guarantees.
Preprocessing is fitted on reference rows and applied unchanged to both regimes. Centering or intercepts may be used within training folds, but held-out residual location remains part of Q Q Q F F F Algorithm 2 Sample-efficient reference-anchored cross-fitted plug-in attribution
Input: Reference data D ( 0 ) \mathcal{D}^{(0)} D ( 0 ) D ( 1 ) \mathcal{D}^{(1)} D ( 1 ) W , z W,z W , z
Output: Fold-averaged attribution Φ ^ \widehat{\Phi} Φ
Split reference rows into fit and moment sets ( D f i t ( 0 ) , I 0 ) (\mathcal{D}^{(0)}_{\mathrm{fit}}, I_0) ( D fit ( 0 ) , I 0 ) W , z W,z W , z D f i t ( 0 ) \mathcal{D}^{(0)}_{\mathrm{fit}} D fit ( 0 )
for k = 1 , … , d k=1,\ldots,d k = 1 , … , d
Tune λ k ( 0 ) \lambda_k^{(0)} λ k ( 0 ) b ^ k ( 0 ) ∈ arg min b { ∣ D f i t ( 0 ) ∣ − 1 ∥ X f i t , k ( 0 ) , c − X f i t , < k ( 0 ) , c b ∥ 2 2 / 2 + λ k ( 0 ) ∥ b ∥ 1 } \widehat{b}_k^{(0)} \in\arg\min_b\{|\mathcal{D}^{(0)}_{\mathrm{fit}}|^{-1}\|X_{\mathrm{fit},k}^{(0),c}-X_{\mathrm{fit},<k}^{(0),c}b\|_2^2/2+\lambda_k^{(0)}\|b\|_1\} b k ( 0 ) ∈ arg min b { ∣ D fit ( 0 ) ∣ − 1 ∥ X fit , k ( 0 ) , c − X fit , < k ( 0 ) , c b ∥ 2 2 /2 + λ k ( 0 ) ∥ b ∥ 1 }
end for
On I 0 I_0 I 0 e ^ i ( 0 ) = ( I d − B ^ ( 0 ) ⊤ ) x ~ i ( 0 ) \widehat{e}_i^{(0)}=(I_d-\widehat{B}^{(0)\top})\widetilde{x}_i^{(0)} e i ( 0 ) = ( I d − B ( 0 ) ⊤ ) x i ( 0 ) ( μ ^ ( 0 ) , ω ^ ( 0 ) ) (\widehat{\mu}^{(0)},\widehat{\omega}^{(0)}) ( μ ( 0 ) , ω ( 0 ) )
Form incident holdout folds I 1 , … , I L I_1,\ldots,I_L I 1 , … , I L J ℓ = [ n 1 ] ∖ I ℓ J_\ell=[n_1]\setminus I_\ell J ℓ = [ n 1 ] ∖ I ℓ
for ℓ = 1 , … , L \ell=1,\ldots,L ℓ = 1 , … , L
for k = 1 , … , d k=1,\ldots,d k = 1 , … , d
y ℓ k ← X ℓ , J ℓ , k ( 1 ) , c − X ℓ , J ℓ , < k ( 1 ) , c b ^ k ( 0 ) y_{\ell k}\leftarrow X_{\ell,J_\ell,k}^{(1),c}-X_{\ell,J_\ell,<k}^{(1),c}\widehat{b}_k^{(0)} y ℓ k ← X ℓ , J ℓ , k ( 1 ) , c − X ℓ , J ℓ , < k ( 1 ) , c b k ( 0 )
Tune λ ℓ k \lambda_{\ell k} λ ℓ k J ℓ J_\ell J ℓ δ ^ ℓ k ∈ arg min δ { ∣ J ℓ ∣ − 1 ∥ y ℓ k − X ℓ , J ℓ , < k ( 1 ) , c δ ∥ 2 2 / 2 + λ ℓ k ∥ δ ∥ 1 } \widehat{\delta}_{\ell k}\in\arg\min_\delta\{|J_\ell|^{-1}\|y_{\ell k}-X_{\ell,J_\ell,<k}^{(1),c}\delta\|_2^2/2+\lambda_{\ell k}\|\delta\|_1\} δ ℓ k ∈ arg min δ { ∣ J ℓ ∣ − 1 ∥ y ℓ k − X ℓ , J ℓ , < k ( 1 ) , c δ ∥ 2 2 /2 + λ ℓ k ∥ δ ∥ 1 } b ^ ℓ k ( 1 ) ← b ^ k ( 0 ) + δ ^ ℓ k \widehat{b}_{\ell k}^{(1)}\leftarrow\widehat{b}_k^{(0)}+\widehat{\delta}_{\ell k} b ℓ k ( 1 ) ← b k ( 0 ) + δ ℓ k
end for
On I ℓ I_\ell I ℓ B ^ ( 1 ) \widehat{B}^{(1)} B ( 1 ) ( μ ^ ℓ ( 1 ) , ω ^ ℓ ( 1 ) ) (\widehat{\mu}_\ell^{(1)},\widehat{\omega}_\ell^{(1)}) ( μ ℓ ( 1 ) , ω ℓ ( 1 ) )
Apply Algorithm 1 to the reference and fold-ℓ \ell ℓ Φ ^ ℓ \widehat{\Phi}_\ell Φ ℓ
end for
return L − 1 ∑ ℓ = 1 L Φ ^ ℓ L^{-1}\sum_{\ell=1}^{L}\widehat{\Phi}_\ell L − 1 ∑ ℓ = 1 L Φ ℓ
uses a separate reference score-fitting split or requires joint control of score estimation. The theorem uses a penalty satisfying the score event; experiments tune inside each training fold only.
Proofs and Technical Lemmas
This appendix collects the proofs used by the main paper. The order follows the reader’s pipeline: efficiency and coherence, exact computation, its fixed-SCM event specialization, high-dimensional definitions, Delta-Lasso, held-out residual moments, and finally plug-in attribution and ranking.
Proof of Proposition 2.1
Proof of Proposition 2.1. By construction, M ∅ = M ( 0 ) M^{\varnothing}=M^{(0)} M ∅ = M ( 0 ) M U = M ( 1 ) M^{\mathcal{U}}=M^{(1)} M U = M ( 1 ) (2) gives
∑ u ∈ U Φ u = T ( M ( 1 ) ) − T ( M ( 0 ) ) . \sum_{u\in\mathcal{U}}\Phi_u=T(M^{(1)})-T(M^{(0)}). u ∈ U ∑ Φ u = T ( M ( 1 ) ) − T ( M ( 0 ) ) . For coherence, fix a coalition U ⊆ U U\subseteq\mathcal{U} U ⊆ U e ∼ Q ( 1 ) e\sim Q^{(1)} e ∼ Q ( 1 ) M δ e U M_{\delta_e}^{U} M δ e U M ( F ( 1 ) , δ e ) M(F^{(1)},\delta_e) M ( F ( 1 ) , δ e ) M Q ( 1 ) U M_{Q^{(1)}}^{U} M Q ( 1 ) U M ( F ( 1 ) , Q ( 1 ) ) M(F^{(1)},Q^{(1)}) M ( F ( 1 ) , Q ( 1 ) ) M δ e U M_{\delta_e}^{U} M δ e U e ∼ Q ( 1 ) e\sim Q^{(1)} e ∼ Q ( 1 ) Q j ( 1 ) Q_j^{(1)} Q j ( 1 ) Q j ( 0 ) Q_j^{(0)} Q j ( 0 )
E e ∼ Q ( 1 ) [ T ( M δ e U ) ] = T ( M Q ( 1 ) U ) . \mathbb{E}_{e\sim Q^{(1)}}\left[T(M_{\delta_e}^{U})\right]=T(M_{Q^{(1)}}^{U}). E e ∼ Q ( 1 ) [ T ( M δ e U ) ] = T ( M Q ( 1 ) U ) . Each Shapley value is a finite linear combination of these coalition values. Taking expectations term by term in (2) proves E e ∼ Q ( 1 ) [ Φ u s h i f t ( e ) ] = Φ u \mathbb{E}_{e\sim Q^{(1)}}[\Phi_u^{\mathrm{shift}}(e)]=\Phi_u E e ∼ Q ( 1 ) [ Φ u shift ( e )] = Φ u u u u
Exact moment-recursion SCM-replacement Shapley: proof of Theorem 3.1
Proof of Theorem 3.1. Use the notation of the theorem. For node j j j b j R p b_j^{R_p} b j R p e j R p e_j^{R_p} e j R p R p R_p R p μ j R p = E [ e j R p ∣ R p ] \mu_j^{R_p} = \mathbb{E}[e_j^{R_p} \mid R_p] μ j R p = E [ e j R p ∣ R p ] ω j R p = Var ( e j R p ∣ R p ) \omega_j^{R_p} = \operatorname{Var}(e_j^{R_p} \mid R_p) ω j R p = Var ( e j R p ∣ R p ) R p R_p R p
x j = ( b j R p ) ⊤ x < j + e j R p . x_j = (b_j^{R_p})^\top x_{<j} + e_j^{R_p}. x j = ( b j R p ) ⊤ x < j + e j R p . Because modules are node-local and independently switched in the multilinear extension, the current node’s switching indicators are independent of x < j x_{<j} x < j
E [ x j ] = E [ ( b j R p ) ⊤ ] E [ x < j ] + E [ e j R p ] = b ˉ j ⊤ x ˉ < j + μ ˉ j . \mathbb{E}[x_j] = \mathbb{E}[(b_j^{R_p})^\top]\mathbb{E}[x_{<j}] + \mathbb{E}[e_j^{R_p}] = \bar{b}_j^\top\bar{x}_{<j} + \bar{\mu}_j. E [ x j ] = E [( b j R p ) ⊤ ] E [ x < j ] + E [ e j R p ] = b ˉ j ⊤ x ˉ < j + μ ˉ j . For the raw cross second moment,
E [ x < j x j ] = E [ x < j x < j ⊤ ] b ˉ j + E [ x < j ] μ ˉ j , \mathbb{E}[x_{<j}x_j] = \mathbb{E}[x_{<j}x_{<j}^\top]\bar{b}_j + \mathbb{E}[x_{<j}]\bar{\mu}_j, E [ x < j x j ] = E [ x < j x < j ⊤ ] b ˉ j + E [ x < j ] μ ˉ j , which is the second line of (13). For the diagonal raw second moment,
E [ x j 2 ] = E [ ( b j R p ) ⊤ x < j x < j ⊤ b j R p ] + 2 E [ μ j R p ( b j R p ) ⊤ ] x ˉ < j + E [ ( e j R p ) 2 ] = tr ( Γ j S < j , < j ) + 2 c j ⊤ x ˉ < j + q j . \begin{aligned}
\mathbb{E}[x_j^2] = \mathbb{E}\left[(b_j^{R_p})^\top x_{<j}x_{<j}^\top b_j^{R_p}\right]
+ 2\mathbb{E}\left[\mu_j^{R_p}(b_j^{R_p})^\top\right]\bar{x}_{<j} + \mathbb{E}[(e_j^{R_p})^2] \\
&= \operatorname{tr}(\Gamma_j S_{<j,<j}) + 2c_j^\top\bar{x}_{<j} + q_j.
\end{aligned} E [ x j 2 ] = E [ ( b j R p ) ⊤ x < j x < j ⊤ b j R p ] + 2 E [ μ j R p ( b j R p ) ⊤ ] x ˉ < j + E [( e j R p ) 2 ] = tr ( Γ j S < j , < j ) + 2 c j ⊤ x ˉ < j + q j . This proves the recursion in (13) for x ˉ ( p ) \bar{x}(p) x ˉ ( p ) S ( p ) S(p) S ( p ) τ ( x ) = x ⊤ W x + z ⊤ x \tau(x) = x^\top W x + z^\top x τ ( x ) = x ⊤ W x + z ⊤ x
T ~ ( p ) = E R p [ T ( M R p ) ] = E [ τ ( x ) ] = tr ( W S ( p ) ) + z ⊤ x ˉ ( p ) , \begin{aligned}
\widetilde{T}(p) = \mathbb{E}_{R_p}\left[T(M^{R_p})\right] \\
&= \mathbb{E}[\tau(x)] \\
&= \operatorname{tr}(W S(p)) + z^\top\bar{x}(p),
\end{aligned} T ( p ) = E R p [ T ( M R p ) ] = E [ τ ( x )] = tr ( W S ( p )) + z ⊤ x ˉ ( p ) , which is (15).
It remains to connect this program to the Shapley weights. Differentiating (11) directly with respect to p u p_u p u u u u
∂ T ~ ( p ) ∂ p u = ∑ U ⊆ U ∖ { u } { T ( M U ∪ { u } ) − T ( M U ) } ∏ v ∈ U p v ∏ v ∈ U ∖ ( U ∪ { u } ) ( 1 − p v ) . \frac{\partial\widetilde{T}(p)}{\partial p_u}
=
\sum_{U\subseteq\mathcal{U}\setminus\{u\}}
\left\{T(M^{U\cup\{u\}})-T(M^U)\right\}
\prod_{v\in U}p_v
\prod_{v\in\mathcal{U}\setminus(U\cup\{u\})}(1-p_v). ∂ p u ∂ T ( p ) = U ⊆ U ∖ { u } ∑ { T ( M U ∪ { u } ) − T ( M U ) } v ∈ U ∏ p v v ∈ U ∖ ( U ∪ { u }) ∏ ( 1 − p v ) . Consequently, with K = ∣ U ∣ K = |\mathcal{U}| K = ∣ U ∣
∂ T ~ ( t 1 ) ∂ p u = ∑ U ⊆ U ∖ { u } { T ( M U ∪ { u } ) − T ( M U ) } × t ∣ U ∣ ( 1 − t ) K − ∣ U ∣ − 1 . \frac{\partial\widetilde{T}(t\mathbf{1})}{\partial p_u}
=
\sum_{U\subseteq\mathcal{U}\setminus\{u\}}
\left\{T(M^{U\cup\{u\}})-T(M^U)\right\}
\times t^{|U|}(1-t)^{K-|U|-1}. ∂ p u ∂ T ( t 1 ) = U ⊆ U ∖ { u } ∑ { T ( M U ∪ { u } ) − T ( M U ) } × t ∣ U ∣ ( 1 − t ) K − ∣ U ∣ − 1 . Integrating term by term over t ∈ [ 0 , 1 ] t \in[0,1] t ∈ [ 0 , 1 ]
∫ 0 1 t ∣ U ∣ ( 1 − t ) K − ∣ U ∣ − 1 d t = ∣ U ∣ ! ( K − ∣ U ∣ − 1 ) ! K ! , \int_0^1 t^{|U|}(1-t)^{K-|U|-1}\,dt
=
\frac{|U|!(K-|U|-1)!}{K!}, ∫ 0 1 t ∣ U ∣ ( 1 − t ) K − ∣ U ∣ − 1 d t = K ! ∣ U ∣ ! ( K − ∣ U ∣ − 1 )! , recovers (2) and proves (14) directly [29 , 21 ] . The displayed derivative is a polynomial in t t t K − 1 K-1 K − 1 m m m 2 m − 1 2m-1 2 m − 1 m = ⌈ K / 2 ⌉ m = \lceil K/2\rceil m = ⌈ K /2 ⌉ [10 ] .
Proof of Corollary 3.2. We give the backward-sensitivity equations for the separated function/law game; paired modules use the same sweep with the mixed local moment c j c_j c j p p p α j \alpha_j α j ρ j \rho_j ρ j γ j \gamma_j γ j x ˉ j \bar{x}_j x ˉ j S j j S_{jj} S j j b j = b ˉ j b_j=\bar{b}_j b j = b ˉ j μ j = μ ˉ j \mu_j=\bar{\mu}_j μ j = μ ˉ j Γ j = E [ b j ( r ) ( b j ( r ) ) ⊤ ] \Gamma_j=\mathbb{E}[b_j^{(r)}(b_j^{(r)})^\top] Γ j = E [ b j ( r ) ( b j ( r ) ) ⊤ ]
Δ b j = b j ( 1 ) − b j ( 0 ) , Δ μ j = μ j ( 1 ) − μ j ( 0 ) , \Delta b_j=b_j^{(1)}-b_j^{(0)},\qquad\Delta\mu_j=\mu_j^{(1)}-\mu_j^{(0)}, Δ b j = b j ( 1 ) − b j ( 0 ) , Δ μ j = μ j ( 1 ) − μ j ( 0 ) , Δ Γ j = b j ( 1 ) ( b j ( 1 ) ) ⊤ − b j ( 0 ) ( b j ( 0 ) ) ⊤ , \Delta\Gamma_j=b_j^{(1)}(b_j^{(1)})^\top-b_j^{(0)}(b_j^{(0)})^\top, Δ Γ j = b j ( 1 ) ( b j ( 1 ) ) ⊤ − b j ( 0 ) ( b j ( 0 ) ) ⊤ , q j ( r ) = ( μ j ( r ) ) 2 + ω j ( r ) , Δ q j = q j ( 1 ) − q j ( 0 ) . q_j^{(r)}=(\mu_j^{(r)})^2+\omega_j^{(r)},\qquad\Delta q_j=q_j^{(1)}-q_j^{(0)}. q j ( r ) = ( μ j ( r ) ) 2 + ω j ( r ) , Δ q j = q j ( 1 ) − q j ( 0 ) . Differentiating the three forward assignments at node j j j
∂ T ~ ∂ p F j = ( α j x ˉ < j + S < j , < j ⊤ ρ j + 2 γ j μ j x ˉ < j ) ⊤ Δ b j + γ j ⟨ Δ Γ j , S < j , < j ⟩ . (28) \frac{\partial\widetilde{T}}{\partial p_{F_j}}=\left(\alpha_j\bar{x}_{<j}+S_{<j,<j}^{\top}\rho_j+2\gamma_j\mu_j\bar{x}_{<j}\right)^{\top}\Delta b_j+\gamma_j\langle\Delta\Gamma_j,S_{<j,<j}\rangle. \tag*{(28)} ∂ p F j ∂ T = ( α j x ˉ < j + S < j , < j ⊤ ρ j + 2 γ j μ j x ˉ < j ) ⊤ Δ b j + γ j ⟨ Δ Γ j , S < j , < j ⟩ . (28) ∂ T ~ ∂ p Q j = ( α j + ρ j ⊤ x ˉ < j + 2 γ j b j ⊤ x ˉ < j ) Δ μ j + γ j Δ q j . (29) \frac{\partial\widetilde{T}}{\partial p_{Q_j}}=\left(\alpha_j+\rho_j^{\top}\bar{x}_{<j}+2\gamma_j b_j^{\top}\bar{x}_{<j}\right)\Delta\mu_j+\gamma_j\Delta q_j. \tag*{(29)} ∂ p Q j ∂ T = ( α j + ρ j ⊤ x ˉ < j + 2 γ j b j ⊤ x ˉ < j ) Δ μ j + γ j Δ q j . (29) Writing a a a H H H a = z a=z a = z H = W H=W H = W
a < j + = α j b j + ρ j μ j + 2 γ j μ j b j , a_{<j}+=\alpha_jb_j+\rho_j\mu_j+2\gamma_j\mu_jb_j, a < j + = α j b j + ρ j μ j + 2 γ j μ j b j , H < j , < j + = ρ j b ˉ j ⊤ + γ j Γ ˉ j . (30) H_{<j,<j}+=\rho_j\bar{b}_j^{\top}+\gamma_j\bar{\Gamma}_j. \tag*{(30)} H < j , < j + = ρ j b ˉ j ⊤ + γ j Γ ˉ j . (30) The two triangle contributions are combined in ρ j \rho_j ρ j 2 d 2d 2 d
The forward and backward matrix–vector/inner-product work summed over nodes is O ( d 3 ) O(d^3) O ( d 3 ) O ( d 2 ) O(d^2) O ( d 2 ) ⌈ K / 2 ⌉ = O ( K ) \lceil K/2\rceil=O(K) ⌈ K /2 ⌉ = O ( K ) O ( K d 3 ) O(Kd^3) O ( K d 3 ) K = 2 d K=2d K = 2 d O ( d 4 ) O(d^4) O ( d 4 ) O ( K 2 d 3 ) = O ( d 5 ) O(K^2d^3)=O(d^5) O ( K 2 d 3 ) = O ( d 5 ) K = d K=d K = d
Fixed-SCM event specialization
Corollary C.1 (Fixed-SCM event attribution). Fix an ordered linear SCM M ( B , Q ) M(B,Q) M ( B , Q ) μ = E Q [ e ] \mu=\mathbb{E}_{Q}[e] μ = E Q [ e ] Ω = Diag ( ω ) \Omega=\operatorname{Diag}(\omega) Ω = Diag ( ω ) A = ( I − B ) − 1 A=(I-B)^{-1} A = ( I − B ) − 1 V = A W A ⊤ V=AWA^{\top} V = A W A ⊤ ℓ = A z \ell=Az ℓ = A z x x x e = ( I − B ⊤ ) x e=(I-B^{\top})x e = ( I − B ⊤ ) x Q Q Q M ( B , Q ) M(B,Q) M ( B , Q ) M ( B , δ e ) M(B,\delta_e) M ( B , δ e )
Φ f i x ( e ) = Diag ( e − μ ) { V ( e + μ ) + ℓ } − diag ( V Ω ) . (31) \Phi^{\mathrm{fix}}(e)=\operatorname{Diag}(e-\mu)\{V(e+\mu)+\ell\}-\operatorname{diag}(V\Omega). \tag*{(31)} Φ fix ( e ) = Diag ( e − μ ) { V ( e + μ ) + ℓ } − diag ( V Ω ) . (31) It satisfies 1 ⊤ Φ f i x ( e ) = τ ( x ) − E X ∼ M ( B , Q ) [ τ ( X ) ] \mathbf{1}^{\top}\Phi^{\mathrm{fix}}(e)=\tau(x)-\mathbb{E}_{X\sim M(B,Q)}[\tau(X)] 1 ⊤ Φ fix ( e ) = τ ( x ) − E X ∼ M ( B , Q ) [ τ ( X )] Embedding the endpoints in the primary function/law-resolved player set adds only null F F F
Proof. This is the fixed-mechanism specialization of Theorem 3.1. Write Δ = e − μ \Delta=e-\mu Δ = e − μ j j j
m j ( p ) = μ j + p j Δ j , m_j(p)=\mu_j+p_j\Delta_j, m j ( p ) = μ j + p j Δ j , q j ( p ) = μ j 2 + ω j + p j ( e j 2 − μ j 2 − ω j ) . q_j(p)=\mu_j^2+\omega_j+p_j(e_j^2-\mu_j^2-\omega_j). q j ( p ) = μ j 2 + ω j + p j ( e j 2 − μ j 2 − ω j ) . Independence gives cross moments m i ( p ) m j ( p ) m_i(p)m_j(p) m i ( p ) m j ( p ) i ≠ j i\ne j i = j
T ~ ( p ) = ∑ j V j j q j ( p ) + ∑ i ≠ j V i j m i ( p ) m j ( p ) + ℓ ⊤ m ( p ) . \widetilde{T}(p)=\sum_j V_{jj}q_j(p)+\sum_{i\ne j}V_{ij}m_i(p)m_j(p)+\ell^{\top}m(p). T ( p ) = j ∑ V j j q j ( p ) + i = j ∑ V ij m i ( p ) m j ( p ) + ℓ ⊤ m ( p ) . Using V = V ⊤ V=V^{\top} V = V ⊤ p = t 1 p=t1 p = t 1
∂ T ~ ( t 1 ) ∂ p j = V j j ( e j 2 − μ j 2 − ω j ) + 2 Δ j ∑ k ≠ j V j k ( μ k + t Δ k ) + ℓ j Δ j . \frac{\partial\widetilde{T}(t1)}{\partial p_j}=V_{jj}(e_j^2-\mu_j^2-\omega_j)+2\Delta_j\sum_{k\ne j}V_{jk}(\mu_k+t\Delta_k)+\ell_j\Delta_j. ∂ p j ∂ T ( t 1 ) = V j j ( e j 2 − μ j 2 − ω j ) + 2 Δ j k = j ∑ V j k ( μ k + t Δ k ) + ℓ j Δ j . Integration over t ∈ [ 0 , 1 ] t\in[0,1] t ∈ [ 0 , 1 ]
( Φ f i x ) j = Δ j { V ( e + μ ) + ℓ } j − V j j ω j , (\Phi^{\mathrm{fix}})_j=\Delta_j\{V(e+\mu)+\ell\}_j-V_{jj}\omega_j, ( Φ fix ) j = Δ j { V ( e + μ ) + ℓ } j − V j j ω j , which is (C-4) because diag ( V Ω ) j = V j j ω j \operatorname{diag}(V\Omega)_j=V_{jj}\omega_j diag ( V Ω ) j = V j j ω j
High-dimensional conditions and auxiliary estimators
Identification versus estimation. The nonsingular parent covariance invoked in the main-text identification paragraph is a population condition, not invertibility of an empirical high-dimensional covariance. With known parent sets, it is needed only on the parent subvector. With only a known order and a parameter class of nodewise sparsity at most s ≥ 1 s\ge1 s ≥ 1
inf δ ≠ 0 ∥ δ ∥ 0 ≤ 2 s δ ⊤ Cov ( x < j ( r ) ) δ ∥ δ ∥ 2 2 > 0 for every r and j > 1. \inf_{\substack{\delta\ne0\\ \|\delta\|_0\le2s}}\frac{\delta^{\top}\operatorname{Cov}(x^{(r)}_{<j})\delta}{\|\delta\|_2^2}>0\quad\text{for every }r\text{ and }j>1. δ = 0 ∥ δ ∥ 0 ≤ 2 s inf ∥ δ ∥ 2 2 δ ⊤ Cov ( x < j ( r ) ) δ > 0 for every r and j > 1. Indeed, the difference of two s s s 2 s 2s 2 s
Notation. Nodewise slope regressions use centered or intercept-residualized designs, denoted by X ( r ) X^{(r)} X ( r ) k k k Z k ( r ) = X : , 1 : ( k − 1 ) ( r ) Z_k^{(r)}=X^{(r)}_{:,1:(k-1)} Z k ( r ) = X : , 1 : ( k − 1 ) ( r ) S k ( r ) = supp ( b k ( r ) ) S_k^{(r)}=\operatorname{supp}(b_k^{(r)}) S k ( r ) = supp ( b k ( r ) ) s r , k = ∣ S k ( r ) ∣ s_{r,k}=|S_k^{(r)}| s r , k = ∣ S k ( r ) ∣ s r = ∑ k s r , k = ∥ B ( r ) ∥ 0 s_r=\sum_k s_{r,k}=\|B^{(r)}\|_0 s r = ∑ k s r , k = ∥ B ( r ) ∥ 0 S ⊆ [ k − 1 ] S\subseteq[k-1] S ⊆ [ k − 1 ] S c = [ k − 1 ] ∖ S S^c=[k-1]\setminus S S c = [ k − 1 ] ∖ S k k k ∥ M ∥ 1 , 1 = ∑ i , j ∣ M i j ∣ \|M\|_{1,1}=\sum_{i,j}|M_{ij}| ∥ M ∥ 1 , 1 = ∑ i , j ∣ M ij ∣ Δ = B ( 1 ) − B ( 0 ) \Delta=B^{(1)}-B^{(0)} Δ = B ( 1 ) − B ( 0 )
Sub-Gaussian noise, residualization, and the joint score event. For a scalar random variable Y Y Y X X X
∥ Y ∥ ψ 2 : = inf { c > 0 : E exp ( Y 2 / c 2 ) ≤ 2 } , \|Y\|_{\psi_2}:=\inf\{c>0:\mathbb{E}\exp(Y^2/c^2)\le2\}, ∥ Y ∥ ψ 2 := inf { c > 0 : E exp ( Y 2 / c 2 ) ≤ 2 } , ∥ X ∥ ψ 2 : = sup ∥ v ∥ 2 = 1 ∥ v ⊤ X ∥ ψ 2 . \|X\|_{\psi_2}:=\sup_{\|v\|_2=1}\|v^{\top}X\|_{\psi_2}. ∥ X ∥ ψ 2 := ∥ v ∥ 2 = 1 sup ∥ v ⊤ X ∥ ψ 2 . Write ε k ( r ) = ( ε 1 k ( r ) , … , ε n r k ( r ) ) ⊤ \varepsilon_k^{(r)}=(\varepsilon_{1k}^{(r)},\ldots,\varepsilon_{n_rk}^{(r)})^{\top} ε k ( r ) = ( ε 1 k ( r ) , … , ε n r k ( r ) ) ⊤ e k ( r ) − μ k ( r ) e_k^{(r)}-\mu_k^{(r)} e k ( r ) − μ k ( r ) P ( r ) P^{(r)} P ( r ) P ( r ) 1 = 0 P^{(r)}1=0 P ( r ) 1 = 0 Z k ( r ) = P ( r ) Z k , r a w ( r ) Z_k^{(r)}=P^{(r)}Z_{k,\mathrm{raw}}^{(r)} Z k ( r ) = P ( r ) Z k , raw ( r ) ε k ( r ) \varepsilon_k^{(r)} ε k ( r ) i ∈ [ n r ] i\in[n_r] i ∈ [ n r ] k ∈ [ d ] k\in[d] k ∈ [ d ] t ∈ R t\in\mathbb{R} t ∈ R
E [ exp { t ε i k ( r ) } ∣ Z k , r a w ( r ) ] ≤ exp ( σ r 2 t 2 / 2 ) . (32) \mathbb{E}\left[\exp\{t\varepsilon_{ik}^{(r)}\}\mid Z_{k,\mathrm{raw}}^{(r)}\right]\le\exp(\sigma_r^2t^2/2). \tag*{(32)} E [ exp { t ε ik ( r ) } ∣ Z k , raw ( r ) ] ≤ exp ( σ r 2 t 2 /2 ) . (32) Column normalization means
max k ∈ [ d ] max j < k ∥ Z k , : , j ( r ) ∥ 2 n r ≤ 1. (33) \max_{k\in[d]}\max_{j<k}\frac{\|Z_{k,:,j}^{(r)}\|_2}{\sqrt{n_r}}\le1. \tag*{(33)} k ∈ [ d ] max j < k max n r ∥ Z k , : , j ( r ) ∥ 2 ≤ 1. (33) The residualized error ε ˉ k ( r ) : = P ( r ) ε k ( r ) \bar{\varepsilon}_k^{(r)} := P^{(r)}\varepsilon_k^{(r)} ε ˉ k ( r ) := P ( r ) ε k ( r )
Z k ( r ) ⊤ ε ˉ k ( r ) = Z k ( r ) ⊤ ε k ( r ) . Z_k^{(r)\top}\bar{\varepsilon}_k^{(r)} = Z_k^{(r)\top}\varepsilon_k^{(r)}. Z k ( r ) ⊤ ε ˉ k ( r ) = Z k ( r ) ⊤ ε k ( r ) . For regime r r r
λ r = 2 σ r 2 log ( d 2 / δ r ) n r , (34) \lambda_r = 2\sigma_r\sqrt{\frac{2\log(d^2/\delta_r)}{n_r}}, \tag*{(34)} λ r = 2 σ r n r 2 log ( d 2 / δ r ) , (34) or any larger order-equivalent value, where 0 < δ r < 1 0 < \delta_r < 1 0 < δ r < 1
max k ∈ [ d ] max j < k ∣ n r − 1 Z k , : , j ( r ) ⊤ ε k ( r ) ∣ ≤ λ r / 2 (35) \max_{k\in[d]}\max_{j<k}\left|n_r^{-1}Z_{k,:,j}^{(r)\top}\varepsilon_k^{(r)}\right| \le\lambda_r/2 \tag*{(35)} k ∈ [ d ] max j < k max n r − 1 Z k , : , j ( r ) ⊤ ε k ( r ) ≤ λ r /2 (35) with probability at least 1 − δ r 1-\delta_r 1 − δ r δ r → 0 \delta_r \to0 δ r → 0 log ( 1 / δ r ) = O ( log d ) \log(1/\delta_r) = O(\log d) log ( 1/ δ r ) = O ( log d ) λ r ≍ σ r log d / n r \lambda_r \asymp\sigma_r\sqrt{\log d/n_r} λ r ≍ σ r log d / n r
Restricted curvature. For S ⊆ [ k − 1 ] S \subseteq[k-1] S ⊆ [ k − 1 ]
C ( S ; 3 ) : = { θ ∈ R k − 1 : ∥ θ S c ∥ 1 ≤ 3 ∥ θ S ∥ 1 } . (36) \mathcal{C}(S;3) := \{\theta\in\mathbb{R}^{k-1}: \|\theta_{S^c}\|_1 \le3\|\theta_S\|_1\}. \tag*{(36)} C ( S ; 3 ) := { θ ∈ R k − 1 : ∥ θ S c ∥ 1 ≤ 3∥ θ S ∥ 1 } . (36) The cone restricted-eigenvalue condition used for sparse reference and incident-only fits is Bickel et al. [2 ]
1 n ∥ Z θ ∥ 2 2 ≥ κ ∥ θ ∥ 2 2 for every θ ∈ C ( S ; 3 ) . (37) \frac{1}{n}\|Z\theta\|_2^2 \ge\kappa\|\theta\|_2^2 \quad\text{for every }\theta\in\mathcal{C}(S;3). \tag*{(37)} n 1 ∥ Z θ ∥ 2 2 ≥ κ ∥ θ ∥ 2 2 for every θ ∈ C ( S ; 3 ) . (37) Our κ \kappa κ ∥ Z θ ∥ 2 / n \|Z\theta\|_2/\sqrt{n} ∥ Z θ ∥ 2 / n [19 ] :
1 n ∥ Z θ ∥ 2 2 ≥ κ 1 ∥ θ ∥ 2 2 − κ 2 log d n ∥ θ ∥ 1 2 for every θ ∈ R k − 1 . (38) \frac{1}{n}\|Z\theta\|_2^2 \ge\kappa_1\|\theta\|_2^2-\kappa_2\frac{\log d}{n}\|\theta\|_1^2 \quad\text{for every }\theta\in\mathbb{R}^{k-1}. \tag*{(38)} n 1 ∥ Z θ ∥ 2 2 ≥ κ 1 ∥ θ ∥ 2 2 − κ 2 n log d ∥ θ ∥ 1 2 for every θ ∈ R k − 1 . (38) We call (C-10) RE and (C-11) RSC below.
These conditions are assumptions, not consequences of SCM linearity. They hold with high probability for Gaussian random designs under well-conditioned population covariance, typically once n ≳ s log d n \gtrsim s\log d n ≳ s log d [26 ] ; Negahban et al. [19 ] . A union bound over O ( d 2 ) O(d^2) O ( d 2 ) log d \log d log d
Auxiliary Lasso estimators. The pre-incident reference estimator used by Delta-Lasso is
b ~ k ∈ arg min b { 1 2 n 0 ∥ X : , k ( 0 ) − Z k ( 0 ) b ∥ 2 2 + λ 0 ∥ b ∥ 1 } . (39) \widetilde{b}_k \in\underset{b}{\arg\min}\left\{\frac{1}{2n_0}\|\mathbf{X}_{:,k}^{(0)}-\mathbf{Z}_k^{(0)}b\|_2^2+\lambda_0\|b\|_1\right\}. \tag*{(39)} b k ∈ b arg min { 2 n 0 1 ∥ X : , k ( 0 ) − Z k ( 0 ) b ∥ 2 2 + λ 0 ∥ b ∥ 1 } . (39) For the incident-only comparison in Lemma C.2, define
b ^ k n a i v e ∈ arg min b { 1 2 n 1 ∥ X : , k ( 1 ) − Z k ( 1 ) b ∥ 2 2 + λ 1 ∥ b ∥ 1 } . (40) \widehat{b}_k^{\mathrm{naive}} \in\underset{b}{\arg\min}\left\{\frac{1}{2n_1}\|\mathbf{X}_{:,k}^{(1)}-\mathbf{Z}_k^{(1)}b\|_2^2+\lambda_1\|b\|_1\right\}. \tag*{(40)} b k naive ∈ b arg min { 2 n 1 1 ∥ X : , k ( 1 ) − Z k ( 1 ) b ∥ 2 2 + λ 1 ∥ b ∥ 1 } . (40) Reference, Delta-Lasso, and incident-only bounds
Lemma C.2 (Standard nodewise Lasso bounds). Under the joint score event and the nodewise RE condition,
∥ B ~ − B ( 0 ) ∥ 1 , 1 ≤ 12 λ 0 s 0 κ 0 , (41) \lVert\widetilde{B} - B^{(0)} \rVert_{1,1} \le\frac{12\lambda_0 s_0}{\kappa_0}, \tag*{(41)} ∥ B − B ( 0 ) ∥ 1 , 1 ≤ κ 0 12 λ 0 s 0 , (41) ∥ B ~ − B ( 0 ) ∥ F 2 ≤ 9 λ 0 2 s 0 κ 0 2 , (42) \lVert\widetilde{B} - B^{(0)} \rVert_{F}^{2} \le\frac{9\lambda_0^{2}s_0}{\kappa_0^{2}}, \tag*{(42)} ∥ B − B ( 0 ) ∥ F 2 ≤ κ 0 2 9 λ 0 2 s 0 , (42) ∥ B ^ n a i v e − B ( 1 ) ∥ F 2 ≤ 9 λ 1 2 s 1 κ 2 . (43) \lVert\widehat{B}^{\mathrm{naive}} - B^{(1)} \rVert_{F}^{2} \le\frac{9\lambda_1^{2}s_1}{\kappa^{2}}. \tag*{(43)} ∥ B naive − B ( 1 ) ∥ F 2 ≤ κ 2 9 λ 1 2 s 1 . (43) With (C-7), the reference and incident bounds hold with probabilities at least 1 − δ 0 1-\delta_0 1 − δ 0 1 − δ 1 1-\delta_1 1 − δ 1
Proof. We give the common argument to make the constants and simultaneous event explicit. For either sparse endpoint, suppress r , k r,k r , k h ^ = b ^ − b \widehat{h}=\widehat{b}-b h = b − b S = supp ( b ) S=\operatorname{supp}(b) S = supp ( b ) y = Z b + ε ‾ y=Zb+\overline{\varepsilon} y = Z b + ε b ^ \widehat{b} b b b b ∥ Z T ε ‾ / n ∥ ∞ ≤ λ / 2 \lVert Z^{\mathsf{T}}\overline{\varepsilon}/n\rVert_{\infty}\le\lambda/2 ∥ Z T ε / n ∥ ∞ ≤ λ /2 [2 , 19 ]
∥ Z h ^ ∥ 2 2 2 n + λ 2 ∥ h ^ S c ∥ 1 ≤ 3 λ 2 ∥ h ^ S ∥ 1 . (44) \frac{\lVert Z\widehat{h}\rVert_2^2}{2n}+\frac{\lambda}{2}\lVert\widehat{h}_{S^c}\rVert_1\le\frac{3\lambda}{2}\lVert\widehat{h}_S\rVert_1. \tag*{(44)} 2 n ∥ Z h ∥ 2 2 + 2 λ ∥ h S c ∥ 1 ≤ 2 3 λ ∥ h S ∥ 1 . (44) Hence h ^ ∈ C ( S ; 3 ) \widehat{h}\in\mathcal{C}(S;3) h ∈ C ( S ; 3 ) ∥ h ^ S ∥ 1 ≤ ∣ S ∣ ∥ h ^ ∥ 2 \lVert\widehat{h}_S\rVert_1\le\sqrt{\lvert S\rvert}\lVert\widehat{h}\rVert_2 ∥ h S ∥ 1 ≤ ∣ S ∣ ∥ h ∥ 2
∥ h ^ ∥ 2 ≤ 3 λ ∣ S ∣ κ , ∥ h ^ ∥ 1 ≤ 4 ∥ h ^ S ∥ 1 ≤ 12 λ ∣ S ∣ κ . \lVert\widehat{h}\rVert_2\le\frac{3\lambda\sqrt{\lvert S\rvert}}{\kappa},\qquad\lVert\widehat{h}\rVert_1\le4\lVert\widehat{h}_S\rVert_1\le\frac{12\lambda\lvert S\rvert}{\kappa}. ∥ h ∥ 2 ≤ κ 3 λ ∣ S ∣ , ∥ h ∥ 1 ≤ 4 ∥ h S ∥ 1 ≤ κ 12 λ ∣ S ∣ . Applying these inequalities with ( r , λ , κ , S ) = ( 0 , λ 0 , κ 0 , S k ( 0 ) ) (r,\lambda,\kappa,S)=(0,\lambda_0,\kappa_0,S_k^{(0)}) ( r , λ , κ , S ) = ( 0 , λ 0 , κ 0 , S k ( 0 ) ) ℓ 1 \ell_1 ℓ 1 ℓ 2 \ell_2 ℓ 2 k k k ( 1 , λ 1 , κ , S k ( 1 ) ) (1,\lambda_1,\kappa,S_k^{(1)}) ( 1 , λ 1 , κ , S k ( 1 ) ) ℓ 2 \ell_2 ℓ 2 □ \square □
Recall from (17) that
R Δ = ∥ Δ ∥ 1 , 1 + 12 λ 0 s 0 κ 0 . R_{\Delta}=\lVert\Delta\rVert_{1,1}+\frac{12\lambda_0s_0}{\kappa_0}. R Δ = ∥ Δ ∥ 1 , 1 + κ 0 12 λ 0 s 0 . The tolerance condition in Theorem 4.4 is κ 2 R Δ log d / n 1 ≤ c σ 1 \kappa_2R_{\Delta}\sqrt{\log d/n_1}\le c\sigma_1 κ 2 R Δ log d / n 1 ≤ c σ 1
∥ B ( 1 ) − B ~ ∥ 1 , 1 ≤ ∥ B ( 1 ) − B ( 0 ) ∥ 1 , 1 + ∥ B ~ − B ( 0 ) ∥ 1 , 1 ≤ R Δ , (45) \lVert B^{(1)}-\widetilde{B}\rVert_{1,1}\le\lVert B^{(1)}-B^{(0)}\rVert_{1,1}+\lVert\widetilde{B}-B^{(0)}\rVert_{1,1}\le R_{\Delta}, \tag*{(45)} ∥ B ( 1 ) − B ∥ 1 , 1 ≤ ∥ B ( 1 ) − B ( 0 ) ∥ 1 , 1 + ∥ B − B ( 0 ) ∥ 1 , 1 ≤ R Δ , (45) Proof of Theorem 4.4. Work on the reference event above and condition on the independent fit B ~ \widetilde{B} B k k k Δ ~ k = b k ( 1 ) − b ~ k \widetilde{\Delta}_k=b_k^{(1)}-\widetilde{b}_k Δ k = b k ( 1 ) − b k Δ k = δ ^ k − Δ ~ k \Delta_k=\widehat{\delta}_k-\widetilde{\Delta}_k Δ k = δ k − Δ k y k = Z k ( 1 ) Δ ~ k + ε ‾ k ( 1 ) y_k=Z_k^{(1)}\widetilde{\Delta}_k+\overline{\varepsilon}_k^{(1)} y k = Z k ( 1 ) Δ k + ε k ( 1 ) δ ^ k \widehat{\delta}_k δ k Δ ~ k \widetilde{\Delta}_k Δ k
1 2 n 1 ∥ ε ‾ k ( 1 ) − Z k ( 1 ) Δ ^ k ∥ 2 2 + λ 1 ∥ δ ^ k ∥ 1 ≤ 1 2 n 1 ∥ ε ‾ k ( 1 ) ∥ 2 2 + λ 1 ∥ Δ ~ k ∥ 1 . \frac{1}{2n_1}\lVert\overline{\varepsilon}_k^{(1)}-Z_k^{(1)}\widehat{\Delta}_k\rVert_2^2+\lambda_1\lVert\widehat{\delta}_k\rVert_1\le\frac{1}{2n_1}\lVert\overline{\varepsilon}_k^{(1)}\rVert_2^2+\lambda_1\lVert\widetilde{\Delta}_k\rVert_1. 2 n 1 1 ∥ ε k ( 1 ) − Z k ( 1 ) Δ k ∥ 2 2 + λ 1 ∥ δ k ∥ 1 ≤ 2 n 1 1 ∥ ε k ( 1 ) ∥ 2 2 + λ 1 ∥ Δ k ∥ 1 . After expanding and cancelling the noise norm,
1 2 n 1 ∥ Z k ( 1 ) Δ k ∥ 2 2 + λ 1 ∥ δ ^ k ∥ 1 ≤ 1 n 1 ( ε ‾ k ( 1 ) ) T Z k ( 1 ) Δ k + λ 1 ∥ Δ ~ k ∥ 1 . \frac{1}{2n_1}\lVert Z_k^{(1)}\Delta_k\rVert_2^2+\lambda_1\lVert\widehat{\delta}_k\rVert_1\le\frac{1}{n_1}(\overline{\varepsilon}_k^{(1)})^{\mathsf{T}}Z_k^{(1)}\Delta_k+\lambda_1\lVert\widetilde{\Delta}_k\rVert_1. 2 n 1 1 ∥ Z k ( 1 ) Δ k ∥ 2 2 + λ 1 ∥ δ k ∥ 1 ≤ n 1 1 ( ε k ( 1 ) ) T Z k ( 1 ) Δ k + λ 1 ∥ Δ k ∥ 1 . The incident joint score event is the corresponding simultaneous event obtained by the union bound in C.4. Hence
1 n 1 e ~ k ( 1 ) T Z k ( 1 ) Δ ^ k ≤ λ 1 2 ∥ Δ ^ k ∥ 1 . \frac{1}{n_1}\tilde{e}_k^{(1)\mathsf{T}}Z_k^{(1)}\widehat{\Delta}_k \leq\frac{\lambda_1}{2}\lVert\widehat{\Delta}_k\rVert_1. n 1 1 e ~ k ( 1 ) T Z k ( 1 ) Δ k ≤ 2 λ 1 ∥ Δ k ∥ 1 . Together with ∥ Δ ^ k ∥ 1 ≤ ∥ δ ~ k ∥ 1 + ∥ Δ ~ k ∥ 1 \lVert\widehat{\Delta}_k\rVert_1 \leq\lVert\widetilde{\delta}_k\rVert_1 + \lVert\widetilde{\Delta}_k\rVert_1 ∥ Δ k ∥ 1 ≤ ∥ δ k ∥ 1 + ∥ Δ k ∥ 1
1 2 n 1 ∥ Z k ( 1 ) Δ ^ k ∥ 2 2 + λ 1 2 ∥ δ ~ k ∥ 1 ≤ 3 λ 1 2 ∥ Δ ~ k ∥ 1 . (46) \frac{1}{2n_1}\lVert Z_k^{(1)}\widehat{\Delta}_k\rVert_2^2 + \frac{\lambda_1}{2}\lVert\widetilde{\delta}_k\rVert_1 \leq\frac{3\lambda_1}{2}\lVert\widetilde{\Delta}_k\rVert_1. \tag*{(46)} 2 n 1 1 ∥ Z k ( 1 ) Δ k ∥ 2 2 + 2 λ 1 ∥ δ k ∥ 1 ≤ 2 3 λ 1 ∥ Δ k ∥ 1 . (46) In particular,
∥ Δ ^ k ∥ 1 ≤ 4 ∥ Δ ~ k ∥ 1 , 1 n 1 ∥ Z k ( 1 ) Δ ^ k ∥ 2 2 ≤ 3 λ 1 ∥ Δ ~ k ∥ 1 . (47) \lVert\widehat{\Delta}_k\rVert_1 \leq4\lVert\widetilde{\Delta}_k\rVert_1,
\qquad
\frac{1}{n_1}\lVert Z_k^{(1)}\widehat{\Delta}_k\rVert_2^2 \leq3\lambda_1\lVert\widetilde{\Delta}_k\rVert_1. \tag*{(47)} ∥ Δ k ∥ 1 ≤ 4 ∥ Δ k ∥ 1 , n 1 1 ∥ Z k ( 1 ) Δ k ∥ 2 2 ≤ 3 λ 1 ∥ Δ k ∥ 1 . (47) Applying the incident RSC condition (C-11), as in the standard restricted-convexity analysis of regularized estimators [19 ] , gives
κ 1 ∥ Δ ^ k ∥ 2 2 ≤ 3 λ 1 ∥ Δ ~ k ∥ 1 + 16 κ 2 ∥ Δ ~ k ∥ 1 2 log d n 1 . (48) \kappa_1\lVert\widehat{\Delta}_k\rVert_2^2 \leq3\lambda_1\lVert\widetilde{\Delta}_k\rVert_1 + 16\kappa_2\lVert\widetilde{\Delta}_k\rVert_1^2\frac{\log d}{n_1}. \tag*{(48)} κ 1 ∥ Δ k ∥ 2 2 ≤ 3 λ 1 ∥ Δ k ∥ 1 + 16 κ 2 ∥ Δ k ∥ 1 2 n 1 log d . (48) Because Δ ^ k \widehat{\Delta}_k Δ k k k k B ^ δ − B ( 1 ) \widehat{B}^{\delta}-B^{(1)} B δ − B ( 1 ) k k k ∑ k ∥ Δ ~ k ∥ 1 ≤ R Δ \sum_k\lVert\widetilde{\Delta}_k\rVert_1 \leq R_{\Delta} ∑ k ∥ Δ k ∥ 1 ≤ R Δ ∑ k ∥ Δ ~ k ∥ 1 2 ≤ R Δ 2 \sum_k\lVert\widetilde{\Delta}_k\rVert_1^2 \leq R_{\Delta}^2 ∑ k ∥ Δ k ∥ 1 2 ≤ R Δ 2
∥ B ^ δ − B ( 1 ) ∥ F 2 ≤ 3 λ 1 R Δ + 16 κ 2 R Δ 2 log d / n 1 κ 1 . \lVert\widehat{B}^{\delta}-B^{(1)}\rVert_F^2 \leq\frac{3\lambda_1R_{\Delta}+16\kappa_2R_{\Delta}^2\log d/n_1}{\kappa_1}. ∥ B δ − B ( 1 ) ∥ F 2 ≤ κ 1 3 λ 1 R Δ + 16 κ 2 R Δ 2 log d / n 1 . The tolerance condition and λ 1 ≍ σ 1 log d / n 1 \lambda_1\asymp\sigma_1\sqrt{\log d/n_1} λ 1 ≍ σ 1 log d / n 1 C λ 1 R Δ C\lambda_1R_{\Delta} C λ 1 R Δ 1 − δ 0 − δ 1 1-\delta_0-\delta_1 1 − δ 0 − δ 1
The tolerance condition is implied by n 1 ≳ ( κ 2 R Δ / σ 1 ) 2 log d n_1\gtrsim(\kappa_2R_{\Delta}/\sigma_1)^2\log d n 1 ≳ ( κ 2 R Δ / σ 1 ) 2 log d ℓ 1 \ell_1 ℓ 1
Held-out residual-moment estimation
We first bound the held-out residual mean and variance for one fitted regime, then substitute the mechanism-estimation rates used in the main theorem. Fold independence is a sampling condition of the out-of-sample construction, not an additional moment restriction: conditional on the mechanism estimate, the evaluation rows are independent and follow the corresponding regime law.
Proposition C.3 (Residual mean and variance estimation). *For regime r r r χ r = ∥ A ( r ) ∥ o p σ r \chi_r=\lVert A^{(r)}\rVert_{\mathrm{op}}\sigma_r χ r = ∥ A ( r ) ∥ op σ r M r = ∥ E x ( r ) ∥ 2 M_r=\lVert\mathbb{E}x^{(r)}\rVert_2 M r = ∥ E x ( r ) ∥ 2 m r m_r m r r r r r r r η ˉ B = max j ∥ b ^ j − b j ∥ 2 \bar{\eta}_B=\max_j\lVert\widehat{b}_j-b_j\rVert_2 η ˉ B = max j ∥ b j − b j ∥ 2 u r , ζ = log ( 4 d / ζ ) / m r u_{r,\zeta}=\sqrt{\log(4d/\zeta)/m_r} u r , ζ = log ( 4 d / ζ ) / m r 0 < ζ < 1 0<\zeta<1 0 < ζ < 1 1 − ζ 1-\zeta 1 − ζ
∥ μ ^ − μ ∥ ∞ ≤ C { σ r u r , ζ + ( M r + χ r u r , ζ ) η ˉ B } , (49) \lVert\widehat{\mu}-\mu\rVert_\infty\leq C\{\sigma_ru_{r,\zeta}+(M_r+\chi_ru_{r,\zeta})\bar{\eta}_B\}, \tag*{(49)} ∥ μ − μ ∥ ∞ ≤ C { σ r u r , ζ + ( M r + χ r u r , ζ ) η ˉ B } , (49) ∥ ω ^ − ω ∥ ∞ ≤ C { σ r 2 ( u r , ζ + u r , ζ 2 ) + σ r χ r η ˉ B ( u r , ζ + u r , ζ 2 ) + χ r 2 η ˉ B 2 ( 1 + u r , ζ + u r , ζ 2 ) } . (50) \lVert\widehat{\omega}-\omega\rVert_\infty\leq C\{\sigma_r^2(u_{r,\zeta}+u_{r,\zeta}^2)+\sigma_r\chi_r\bar{\eta}_B(u_{r,\zeta}+u_{r,\zeta}^2)+\chi_r^2\bar{\eta}_B^2(1+u_{r,\zeta}+u_{r,\zeta}^2)\}. \tag*{(50)} ∥ ω − ω ∥ ∞ ≤ C { σ r 2 ( u r , ζ + u r , ζ 2 ) + σ r χ r η ˉ B ( u r , ζ + u r , ζ 2 ) + χ r 2 η ˉ B 2 ( 1 + u r , ζ + u r , ζ 2 )} . (50) Here and below, C > 0 C>0 C > 0
Proof. Condition on B ^ \widehat{B} B h j = b ^ j − b j h_j=\widehat{b}_j-b_j h j = b j − b j ξ i = x i − E x \xi_i=x_i-\mathbb{E}x ξ i = x i − E x x i , < j = E x < j + ξ i , < j x_{i,<j}=\mathbb{E}x_{<j}+\xi_{i,<j} x i , < j = E x < j + ξ i , < j e i j = μ j + ε i j e_{ij}=\mu_j+\varepsilon_{ij} e ij = μ j + ε ij ∥ ε i j ∥ ψ 2 ≤ C σ r \lVert\varepsilon_{ij}\rVert_{\psi_2}\le C\sigma_r ∥ ε ij ∥ ψ 2 ≤ C σ r x − E x = A ( r ) T ( e − μ ) x-\mathbb{E}x=A^{(r)\mathsf{T}}(e-\mu) x − E x = A ( r ) T ( e − μ )
max j ∥ ξ i , < j ∥ ψ 2 ≤ C χ r , max j ∥ E x < j ∥ 2 ≤ M r . \max_j\lVert\xi_{i,<j}\rVert_{\psi_2}\le C\chi_r,\qquad\max_j\lVert\mathbb{E}x_{<j}\rVert_2\le M_r. j max ∥ ξ i , < j ∥ ψ 2 ≤ C χ r , j max ∥ E x < j ∥ 2 ≤ M r . Moreover,
e ^ i j = μ j − h j T E x < j + ε i j + q i j , q i j = − h j T ξ i , < j . \widehat{e}_{ij}=\mu_j-h_j^{\mathsf{T}}\mathbb{E}x_{<j}+\varepsilon_{ij}+q_{ij},\qquad q_{ij}=-h_j^{\mathsf{T}}\xi_{i,<j}. e ij = μ j − h j T E x < j + ε ij + q ij , q ij = − h j T ξ i , < j . where ∥ q i j ∥ ψ 2 ≤ C χ r ∥ h j ∥ 2 \lVert q_{ij}\rVert_{\psi_2}\le C\chi_r\lVert h_j\rVert_2 ∥ q ij ∥ ψ 2 ≤ C χ r ∥ h j ∥ 2
μ ^ j − μ j = − h j T E x < j + ε ‾ j + q ‾ j . \widehat{\mu}_j-\mu_j=-h_j^{\mathsf{T}}\mathbb{E}x_{<j}+\overline{\varepsilon}_j+\overline{q}_j. μ j − μ j = − h j T E x < j + ε j + q j . Sub-Gaussian concentration gives, simultaneously over j j j ∣ ε ‾ j ∣ ≤ C σ r u r , ζ |\overline{\varepsilon}_j|\le C\sigma_r u_{r,\zeta} ∣ ε j ∣ ≤ C σ r u r , ζ ∣ q ‾ j ∣ ≤ C χ r ∥ h j ∥ 2 u r , ζ |\overline{q}_j|\le C\chi_r\lVert h_j\rVert_2u_{r,\zeta} ∣ q j ∣ ≤ C χ r ∥ h j ∥ 2 u r , ζ ∣ h j T E x < j ∣ ≤ M r ∥ h j ∥ 2 |h_j^{\mathsf{T}}\mathbb{E}x_{<j}|\le M_r\lVert h_j\rVert_2 ∣ h j T E x < j ∣ ≤ M r ∥ h j ∥ 2
The constant term − h j T E x < j -h_j^{\mathsf{T}}\mathbb{E}x_{<j} − h j T E x < j
ω ^ j − ω j = m r − 1 ∑ i ( ε i j 2 − ω j ) + 2 m r − 1 ∑ i ε i j q i j \widehat{\omega}_j-\omega_j=m_r^{-1}\sum_i(\varepsilon_{ij}^{2}-\omega_j)+2m_r^{-1}\sum_i\varepsilon_{ij}q_{ij} ω j − ω j = m r − 1 i ∑ ( ε ij 2 − ω j ) + 2 m r − 1 i ∑ ε ij q ij + m r − 1 ∑ i q i j 2 − ( ε ‾ j + q ‾ j ) 2 . {}+m_r^{-1}\sum_iq_{ij}^{2}-(\overline{\varepsilon}_j+\overline{q}_j)^2. + m r − 1 i ∑ q ij 2 − ( ε j + q j ) 2 . Exogenous independence gives E ( ε i j q i j ) = 0 \mathbb{E}(\varepsilon_{ij}q_{ij})=0 E ( ε ij q ij ) = 0 ∥ ε i j 2 ∥ ψ 1 ≤ C σ r 2 \lVert\varepsilon_{ij}^{2}\rVert_{\psi_1}\le C\sigma_r^2 ∥ ε ij 2 ∥ ψ 1 ≤ C σ r 2 ∥ ε i j q i j ∥ ψ 1 ≤ C σ r χ r ∥ h j ∥ 2 \lVert\varepsilon_{ij}q_{ij}\rVert_{\psi_1}\le C\sigma_r\chi_r\lVert h_j\rVert_2 ∥ ε ij q ij ∥ ψ 1 ≤ C σ r χ r ∥ h j ∥ 2 ∥ q i j 2 ∥ ψ 1 ≤ C χ r 2 ∥ h j ∥ 2 2 \lVert q_{ij}^{2}\rVert_{\psi_1}\le C\chi_r^2\lVert h_j\rVert_2^2 ∥ q ij 2 ∥ ψ 1 ≤ C χ r 2 ∥ h j ∥ 2 2 C σ r 2 ( u r , ζ + u r , ζ 2 ) C\sigma_r^2(u_{r,\zeta}+u_{r,\zeta}^2) C σ r 2 ( u r , ζ + u r , ζ 2 ) C σ r χ r ∥ h j ∥ 2 ( u r , ζ + u r , ζ 2 ) C\sigma_r\chi_r\lVert h_j\rVert_2(u_{r,\zeta}+u_{r,\zeta}^2) C σ r χ r ∥ h j ∥ 2 ( u r , ζ + u r , ζ 2 )
m r − 1 ∑ i q i j 2 ≤ C χ r 2 ∥ h j ∥ 2 2 ( 1 + u r , ζ + u r , ζ 2 ) . m_r^{-1}\sum_iq_{ij}^{2}\le C\chi_r^2\lVert h_j\rVert_2^2(1+u_{r,\zeta}+u_{r,\zeta}^2). m r − 1 i ∑ q ij 2 ≤ C χ r 2 ∥ h j ∥ 2 2 ( 1 + u r , ζ + u r , ζ 2 ) . Finally, ( ε ‾ j + q ‾ j ) 2 ≤ C ( σ r + χ r ∥ h j ∥ 2 ) 2 u r , ζ 2 (\overline{\varepsilon}_j+\overline{q}_j)^2\le C(\sigma_r+\chi_r\lVert h_j\rVert_2)^2u_{r,\zeta}^2 ( ε j + q j ) 2 ≤ C ( σ r + χ r ∥ h j ∥ 2 ) 2 u r , ζ 2 d d d ∥ h j ∥ 2 ≤ η ‾ B \lVert h_j\rVert_2\le\overline{\eta}_B ∥ h j ∥ 2 ≤ η B
Corollary C.4 (Regimewise held-out moment rates). Define η B ( r ) = ∥ B ^ ( r ) − B ( r ) ∥ F \eta_B^{(r)}=\lVert\widehat{B}^{(r)}-B^{(r)}\rVert_{\mathrm{F}} η B ( r ) = ∥ B ( r ) − B ( r ) ∥ F η μ ( r ) = ∥ μ ^ ( r ) − μ ( r ) ∥ ∞ \eta_\mu^{(r)}=\lVert\widehat{\mu}^{(r)}-\mu^{(r)}\rVert_\infty η μ ( r ) = ∥ μ ( r ) − μ ( r ) ∥ ∞ η ω ( r ) = ∥ ω ^ ( r ) − ω ( r ) ∥ ∞ \eta_\omega^{(r)}=\lVert\widehat{\omega}^{(r)}-\omega^{(r)}\rVert_\infty η ω ( r ) = ∥ ω ( r ) − ω ( r ) ∥ ∞ Under Proposition C.3 , suppose η B ( r ) = O P ( ρ r ) \eta_B^{(r)}=O_{\mathrm{P}}(\rho_r) η B ( r ) = O P ( ρ r ) ν r = log ( 2 d ) / m r \nu_r=\sqrt{\log(2d)/m_r} ν r = log ( 2 d ) / m r
η μ ( r ) = O P ( σ r ν r + ( M r + χ r ν r ) ρ r ) , (51) \eta_\mu^{(r)}=O_{\mathrm{P}}\bigl(\sigma_r\nu_r+(M_r+\chi_r\nu_r)\rho_r\bigr), \tag*{(51)} η μ ( r ) = O P ( σ r ν r + ( M r + χ r ν r ) ρ r ) , (51) η ω ( r ) = O P ( σ r 2 ( ν r + ν r 2 ) + σ r χ r ρ r ( ν r + ν r 2 ) + χ r 2 ρ r 2 ( 1 + ν r + ν r 2 ) ) . (52) \eta_\omega^{(r)}=O_{\mathrm{P}}\bigl(\sigma_r^2(\nu_r+\nu_r^2)+\sigma_r\chi_r\rho_r(\nu_r+\nu_r^2)
+\chi_r^2\rho_r^2(1+\nu_r+\nu_r^2)\bigr). \tag*{(52)} η ω ( r ) = O P ( σ r 2 ( ν r + ν r 2 ) + σ r χ r ρ r ( ν r + ν r 2 ) + χ r 2 ρ r 2 ( 1 + ν r + ν r 2 ) ) . (52) Under Assumption 4.3 , if ν r = o ( 1 ) \nu_r=o(1) ν r = o ( 1 ) η μ ( r ) = O P ( ν r + ( M r + ν r ) ρ r ) \eta_\mu^{(r)}=O_{\mathrm{P}}(\nu_r+(M_r+\nu_r)\rho_r) η μ ( r ) = O P ( ν r + ( M r + ν r ) ρ r ) η ω ( r ) = O P ( ν r + ρ r 2 ) \eta_\omega^{(r)}=O_{\mathrm{P}}(\nu_r+\rho_r^2) η ω ( r ) = O P ( ν r + ρ r 2 )
Proof. Every nodewise coefficient error is bounded by the matrix Frobenius error: max j ∥ b ^ j ( r ) − b j ( r ) ∥ 2 ≤ η B ( r ) \max_j\lVert\widehat{b}_j^{(r)}-b_j^{(r)}\rVert_2\le\eta_B^{(r)} max j ∥ b j ( r ) − b j ( r ) ∥ 2 ≤ η B ( r ) 0 < ζ < 1 0<\zeta<1 0 < ζ < 1 log ( 4 d / ζ ) ≤ C ζ log ( 2 d ) \log(4d/\zeta)\le C_\zeta\log(2d) log ( 4 d / ζ ) ≤ C ζ log ( 2 d ) u r , ζ ≤ C ζ ν r u_{r,\zeta}\le\sqrt{C_\zeta}\nu_r u r , ζ ≤ C ζ ν r O P O_{\mathrm{P}} O P ζ \zeta ζ σ r , χ r = O ( 1 ) \sigma_r,\chi_r=O(1) σ r , χ r = O ( 1 ) ν r = o ( 1 ) \nu_r=o(1) ν r = o ( 1 ) ρ r ν r ≤ ( ρ r 2 + ν r 2 ) / 2 \rho_r\nu_r\le(\rho_r^2+\nu_r^2)/2 ρ r ν r ≤ ( ρ r 2 + ν r 2 ) /2 ν r 2 ≤ ν r \nu_r^2\le\nu_r ν r 2 ≤ ν r 1 + ν r + ν r 2 = O ( 1 ) 1+\nu_r+\nu_r^2=O(1) 1 + ν r + ν r 2 = O ( 1 )
Plug-in attribution sensitivity and ranking
The proof separates four reusable steps: hybrid-path stability, exact endpoint gradients, pathwise norm envelopes, and deterministic rate sequences. Put
D U F = Diag ( 1 { F j ∈ ⋃ w ∈ U w } j = 1 d ) , D_U^F = \operatorname{Diag}\left(\mathbf{1}\left\{F_j \in\bigcup_{w \in U} w\right\}_{j=1}^d\right), D U F = Diag 1 { F j ∈ w ∈ U ⋃ w } j = 1 d , D U Q = Diag ( 1 { Q j ∈ ⋃ w ∈ U w } j = 1 d ) . D_U^Q = \operatorname{Diag}\left(\mathbf{1}\left\{Q_j \in\bigcup_{w \in U} w\right\}_{j=1}^d\right). D U Q = Diag 1 { Q j ∈ w ∈ U ⋃ w } j = 1 d . For r ∈ { 0 , 1 } r \in\{0,1\} r ∈ { 0 , 1 } D U F , ( 1 ) = D U F D_U^{F,(1)} = D_U^F D U F , ( 1 ) = D U F D U F , ( 0 ) = I − D U F D_U^{F,(0)} = I - D_U^F D U F , ( 0 ) = I − D U F D U Q , ( 1 ) = D U Q D_U^{Q,(1)} = D_U^Q D U Q , ( 1 ) = D U Q D U Q , ( 0 ) = I − D U Q D_U^{Q,(0)} = I - D_U^Q D U Q , ( 0 ) = I − D U Q
Proposition C.5 (Derived hybrid-path stability). Let
a ∗ = max U ⊆ U ∥ ( I − B U ) − 1 ∥ o p , a_* = \max_{U \subseteq\mathcal{U}} \left\|(I - B^U)^{-1}\right\|_{\mathrm{op}}, a ∗ = U ⊆ U max ( I − B U ) − 1 op , δ B = ∑ r = 0 1 ∥ B ^ ( r ) − B ( r ) ∥ o p . \delta_B = \sum_{r=0}^{1}\left\|\widehat{B}^{(r)} - B^{(r)}\right\|_{\mathrm{op}}. δ B = r = 0 ∑ 1 B ( r ) − B ( r ) op . For Θ s = Θ + s ( Θ ^ − Θ ) \Theta_s = \Theta+ s(\widehat{\Theta} - \Theta) Θ s = Θ + s ( Θ − Θ ) B s U B_s^U B s U U U U a ∗ δ B < 1 a_*\delta_B < 1 a ∗ δ B < 1
max U , s ∥ ( I − B s U ) − 1 ∥ o p ≤ a ∗ 1 − a ∗ δ B . (53) \max_{U,s}\left\|(I - B_s^U)^{-1}\right\|_{\mathrm{op}} \le\frac{a_*}{1-a_*\delta_B}. \tag*{(53)} U , s max ( I − B s U ) − 1 op ≤ 1 − a ∗ δ B a ∗ . (53) Consequently, the right side is at most 2 a ∗ 2a_* 2 a ∗ ∥ B ^ ( r ) − B ( r ) ∥ F = O P ( ρ r ) \|\widehat{B}^{(r)} - B^{(r)}\|_{\mathrm{F}} = O_{\mathrm{P}}(\rho_r) ∥ B ( r ) − B ( r ) ∥ F = O P ( ρ r ) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 ) a_*(\rho_0+\rho_1) = o(1) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 )
Proof. Every B U B^U B U a ∗ < ∞ a_* < \infty a ∗ < ∞ E B ( r ) = B ^ ( r ) − B ( r ) E_B^{(r)} = \widehat{B}^{(r)} - B^{(r)} E B ( r ) = B ( r ) − B ( r ) D U F D_U^F D U F
E B U : = B ^ U − B U = E B ( 0 ) ( I − D U F ) + E B ( 1 ) D U F , E_B^U := \widehat{B}^U - B^U = E_B^{(0)}(I-D_U^F) + E_B^{(1)}D_U^F, E B U := B U − B U = E B ( 0 ) ( I − D U F ) + E B ( 1 ) D U F , ∥ E B U ∥ o p ≤ δ B . \left\|E_B^U\right\|_{\mathrm{op}} \le\delta_B. E B U op ≤ δ B . Because B \mathcal{B} B B s U B_s^U B s U a ∗ δ B < 1 a_*\delta_B < 1 a ∗ δ B < 1 A U = ( I − B U ) − 1 A^U = (I-B^U)^{-1} A U = ( I − B U ) − 1
I − B s U = ( I − B U ) ( I − s A U E B U ) . I-B_s^U = (I-B^U)(I-sA^UE_B^U). I − B s U = ( I − B U ) ( I − s A U E B U ) . The Neumann bound gives ∥ ( I − s A U E B U ) − 1 ∥ o p ≤ ( 1 − a ∗ δ B ) − 1 \left\|(I-sA^UE_B^U)^{-1}\right\|_{\mathrm{op}} \le(1-a_*\delta_B)^{-1} ( I − s A U E B U ) − 1 op ≤ ( 1 − a ∗ δ B ) − 1 δ B ≤ ∑ r ∥ E B ( r ) ∥ F = O P ( ρ 0 + ρ 1 ) \delta_B \le\sum_r \|E_B^{(r)}\|_{\mathrm{F}} = O_{\mathrm{P}}(\rho_0+\rho_1) δ B ≤ ∑ r ∥ E B ( r ) ∥ F = O P ( ρ 0 + ρ 1 ) a ∗ δ B = o P ( 1 ) a_*\delta_B = o_{\mathrm{P}}(1) a ∗ δ B = o P ( 1 )
Lemma C.6 (Endpoint-gradient representation). Write θ s U = ( B s U , μ s U , ω s U ) \theta_s^U = (B_s^U,\mu_s^U,\omega_s^U) θ s U = ( B s U , μ s U , ω s U ) Θ s \Theta_s Θ s
w U = ∣ U ∣ ! ( ∣ U ∣ − ∣ U ∣ − 1 ) ! ∣ U ∣ ! . w_U = \frac{|U|!(|\mathcal{U}|-|U|-1)!}{|\mathcal{U}|!}. w U = ∣ U ∣ ! ∣ U ∣ ! ( ∣ U ∣ − ∣ U ∣ − 1 )! . At any hybrid state, with A = ( I − B ) − 1 A=(I-B)^{-1} A = ( I − B ) − 1
∇ μ T = 2 A W A ⊤ μ + A z , \nabla_\mu T = 2AWA^\top\mu+ Az, ∇ μ T = 2 A W A ⊤ μ + A z , ∇ ω T = diag ( A W A ⊤ ) , \nabla_\omega T = \operatorname{diag}(AWA^\top), ∇ ω T = diag ( A W A ⊤ ) , ∇ B T = A ⊤ { 2 ( Diag ( ω ) + μ μ ⊤ ) A W + μ z ⊤ } A ⊤ . (54) \nabla_B T = A^\top\left\{2\left(\operatorname{Diag}(\omega)+\mu\mu^\top\right)AW+\mu z^\top\right\}A^\top. \tag*{(54)} ∇ B T = A ⊤ { 2 ( Diag ( ω ) + μ μ ⊤ ) A W + μ z ⊤ } A ⊤ . (54) For r ∈ { 0 , 1 } r \in\{0,1\} r ∈ { 0 , 1 } ϑ ∈ { μ , ω } \vartheta\in\{\mu,\omega\} ϑ ∈ { μ , ω } u u u
∇ B ( r ) Φ u ( Θ s ) = ∑ U ⊆ U ∖ { u } w U { ∇ B T ( θ s U ∪ { u } ) D U ∪ { u } F , ( r ) − ∇ B T ( θ s U ) D U F , ( r ) } , \nabla_{B^{(r)}}\Phi_u(\Theta_s)=\sum_{U\subseteq\mathcal{U}\setminus\{u\}}w_U\left\{\nabla_B T(\theta_s^{U\cup\{u\}})D_{U\cup\{u\}}^{F,(r)}-\nabla_B T(\theta_s^U)D_U^{F,(r)}\right\}, ∇ B ( r ) Φ u ( Θ s ) = U ⊆ U ∖ { u } ∑ w U { ∇ B T ( θ s U ∪ { u } ) D U ∪ { u } F , ( r ) − ∇ B T ( θ s U ) D U F , ( r ) } , ∇ ϑ ( r ) Φ u ( Θ s ) = ∑ U ⊆ U ∖ { u } w U { D U ∪ { u } Q , ( r ) ∇ ϑ T ( θ s U ∪ { u } ) − D U Q , ( r ) ∇ ϑ T ( θ s U ) } . (55) \nabla_{\vartheta^{(r)}}\Phi_u(\Theta_s)=\sum_{U\subseteq\mathcal{U}\setminus\{u\}}w_U\left\{D_{U\cup\{u\}}^{Q,(r)}\nabla_\vartheta T(\theta_s^{U\cup\{u\}})-D_U^{Q,(r)}\nabla_\vartheta T(\theta_s^U)\right\}. \tag*{(55)} ∇ ϑ ( r ) Φ u ( Θ s ) = U ⊆ U ∖ { u } ∑ w U { D U ∪ { u } Q , ( r ) ∇ ϑ T ( θ s U ∪ { u } ) − D U Q , ( r ) ∇ ϑ T ( θ s U ) } . (55) Proof. Let S = Diag ( ω ) + μ μ ⊤ S=\operatorname{Diag}(\omega)+\mu\mu^\top S = Diag ( ω ) + μ μ ⊤
T ( B , μ , ω ) = tr ( S A W A ⊤ ) + μ ⊤ A z . T(B,\mu,\omega)=\operatorname{tr}(SAWA^\top)+\mu^\top Az. T ( B , μ , ω ) = tr ( S A W A ⊤ ) + μ ⊤ A z . For fixed B B B μ \mu μ ω \omega ω d A = A ( d B ) A dA=A(dB)A d A = A ( d B ) A A A A S S S W W W
d B T = ⟨ A ⊤ { 2 S A W + μ z ⊤ } A ⊤ , d B ⟩ F , d_BT=\left\langle A^\top\{2SAW+\mu z^\top\}A^\top,dB\right\rangle_F, d B T = ⟨ A ⊤ { 2 S A W + μ z ⊤ } A ⊤ , d B ⟩ F , which proves the third line.
When only endpoint r r r
d B s U = ( d B ( r ) ) D U F , ( r ) , dB_s^U=(dB^{(r)})D_U^{F,(r)}, d B s U = ( d B ( r ) ) D U F , ( r ) , d μ s U = D U Q , ( r ) d μ ( r ) , d ω s U = D U Q , ( r ) d ω ( r ) . d\mu_s^U=D_U^{Q,(r)}d\mu^{(r)},\qquad d\omega_s^U=D_U^{Q,(r)}d\omega^{(r)}. d μ s U = D U Q , ( r ) d μ ( r ) , d ω s U = D U Q , ( r ) d ω ( r ) . The Shapley value is the finite sum
Φ u ( Θ s ) = ∑ U ⊆ U ∖ { u } w U { T ( θ s U ∪ { u } ) − T ( θ s U ) } . \Phi_u(\Theta_s)=\sum_{U\subseteq\mathcal{U}\setminus\{u\}}w_U\{T(\theta_s^{U\cup\{u\}})-T(\theta_s^U)\}. Φ u ( Θ s ) = U ⊆ U ∖ { u } ∑ w U { T ( θ s U ∪ { u } ) − T ( θ s U )} . Differentiating it term by term and applying the three mask identities proves (C-28). In particular, the before/after endpoint gradients are subtracted before any norm is taken.
Equation (C-28) is the exact player-specific formula: endpoint directions common to the two marginal values cancel. The exact sensitivities can therefore be evaluated by differentiating the moment-and-quadrature program.
For deterministic rate notation, put
q μ , r = ν r + ( M r + ν r ) ρ r , q ω , r = ν r + ρ r 2 , q_{\mu,r}=\nu_r+(M_r+\nu_r)\rho_r,\qquad q_{\omega,r}=\nu_r+\rho_r^2, q μ , r = ν r + ( M r + ν r ) ρ r , q ω , r = ν r + ρ r 2 , and define
m ∗ = 1 + ∑ r = 0 1 { ∥ μ ( r ) ∥ 2 + d q μ , r } , m_* = 1+\sum_{r=0}^{1}\{\|\mu^{(r)}\|_2+\sqrt{d}\,q_{\mu,r}\}, m ∗ = 1 + r = 0 ∑ 1 { ∥ μ ( r ) ∥ 2 + d q μ , r } , ω ∗ = 1 + max r { ∥ ω ( r ) ∥ ∞ + q ω , r } , \omega_* = 1+\max_r\{\|\omega^{(r)}\|_\infty+q_{\omega,r}\}, ω ∗ = 1 + r max { ∥ ω ( r ) ∥ ∞ + q ω , r } , H B = 16 a ∗ 3 ( d ω ∗ + m ∗ 2 ) W ∗ + 4 a ∗ 2 m ∗ z ∗ , H_B=16a_*^3(\sqrt{d}\,\omega_*+m_*^2)W_*+4a_*^2m_*z_*, H B = 16 a ∗ 3 ( d ω ∗ + m ∗ 2 ) W ∗ + 4 a ∗ 2 m ∗ z ∗ , H μ = d { 8 a ∗ 2 W ∗ m ∗ + 2 a ∗ z ∗ } , H ω = 4 d a ∗ 2 W ∗ . H_\mu=\sqrt{d}\{8a_*^2W_*m_*+2a_*z_*\},\qquad H_\omega=4da_*^2W_*. H μ = d { 8 a ∗ 2 W ∗ m ∗ + 2 a ∗ z ∗ } , H ω = 4 d a ∗ 2 W ∗ . Lemma C.7 (Uniform endpoint-gradient rate). Under the conditions of Corollary C.4 in both regimes, assume Assumption 4.3, ν r = o ( 1 ) \nu_r=o(1) ν r = o ( 1 ) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 ) a_*(\rho_0+\rho_1)=o(1) a ∗ ( ρ 0 + ρ 1 ) = o ( 1 )
max u , r L u , ϑ ( r ) = O P ( H ϑ ) , ϑ ∈ { B , μ , ω } . (56) \max_{u,r}L_{u,\vartheta}^{(r)}=O_{\mathbb{P}}(H_\vartheta),\qquad\vartheta\in\{B,\mu,\omega\}. \tag*{(56)} u , r max L u , ϑ ( r ) = O P ( H ϑ ) , ϑ ∈ { B , μ , ω } . (56) Proof. Put α p a t h = sup U , s ∥ ( I − B s U ) − 1 ∥ o p \alpha_{\mathrm{path}}=\sup_{U,s}\|(I-B_s^U)^{-1}\|_{\mathrm{op}} α path = sup U , s ∥ ( I − B s U ) − 1 ∥ op m p a t h = sup U , s ∥ μ s U ∥ 2 m_{\mathrm{path}}=\sup_{U,s}\|\mu_s^U\|_2 m path = sup U , s ∥ μ s U ∥ 2 ω p a t h = sup U , s ∥ ω s U ∥ ∞ \omega_{\mathrm{path}}=\sup_{U,s}\|\omega_s^U\|_\infty ω path = sup U , s ∥ ω s U ∥ ∞
m p a t h ≤ ∑ r = 0 1 { ∥ μ ( r ) ∥ 2 + d η μ ( r ) } , m_{\mathrm{path}}\le\sum_{r=0}^{1}\{\|\mu^{(r)}\|_2+\sqrt{d}\,\eta_{\mu}^{(r)}\}, m path ≤ r = 0 ∑ 1 { ∥ μ ( r ) ∥ 2 + d η μ ( r ) } , ω p a t h ≤ max r = 0 , 1 { ∥ ω ( r ) ∥ ∞ + η ω ( r ) } . \omega_{\mathrm{path}}\le\max_{r=0,1}\{\|\omega^{(r)}\|_\infty+\eta_{\omega}^{(r)}\}. ω path ≤ r = 0 , 1 max { ∥ ω ( r ) ∥ ∞ + η ω ( r ) } . Corollary C.4 gives m p a t h = O P ( m ∗ ) m_{\mathrm{path}}=O_{\mathrm{P}}(m_*) m path = O P ( m ∗ ) ω p a t h = O P ( ω ∗ ) \omega_{\mathrm{path}}=O_{\mathrm{P}}(\omega_*) ω path = O P ( ω ∗ ) α p a t h ≤ 2 a ∗ \alpha_{\mathrm{path}}\le2a_* α path ≤ 2 a ∗
∥ Diag ( ω ) + μ μ ⊤ ∥ F ≤ d ∥ ω ∥ ∞ + ∥ μ ∥ 2 2 , \|\operatorname{Diag}(\omega)+\mu\mu^\top\|_{\mathrm{F}}\le\sqrt{d}\|\omega\|_\infty+\|\mu\|_2^2, ∥ Diag ( ω ) + μ μ ⊤ ∥ F ≤ d ∥ ω ∥ ∞ + ∥ μ ∥ 2 2 , ∥ v ∥ 1 ≤ d ∥ v ∥ 2 , \|v\|_1\le\sqrt{d}\|v\|_2, ∥ v ∥ 1 ≤ d ∥ v ∥ 2 , ∥ diag ( C ) ∥ 1 ≤ d ∥ C ∥ o p . \|\operatorname{diag}(C)\|_1\le d\|C\|_{\mathrm{op}}. ∥ diag ( C ) ∥ 1 ≤ d ∥ C ∥ op . Applying these inequalities to (C-27) bounds the three score-gradient norms by the corresponding H B H_B H B H μ H_\mu H μ H ω H_\omega H ω Π B \Pi_B Π B L L L □ \square □
Proof of Theorem 4.5. Put E ϑ ( r ) = ϑ ^ ( r ) − ϑ ( r ) E_\vartheta^{(r)}=\widehat{\vartheta}^{(r)}-\vartheta^{(r)} E ϑ ( r ) = ϑ ( r ) − ϑ ( r ) ϑ ∈ { B , μ , ω } \vartheta\in\{B,\mu,\omega\} ϑ ∈ { B , μ , ω } B \mathcal{B} B Θ s \Theta_s Θ s
Φ ^ u − Φ u = ∫ 0 1 ∑ r = 0 1 [ ⟨ ∇ B ( r ) Φ u ( Θ s ) , E B ( r ) ⟩ F + ∇ μ ( r ) Φ u ( Θ s ) ⊤ E μ ( r ) + ∇ ω ( r ) Φ u ( Θ s ) ⊤ E ω ( r ) ] d s . \widehat{\Phi}_u-\Phi_u=\int_0^1\sum_{r=0}^{1}\left[\left\langle\nabla_{B^{(r)}}\Phi_u(\Theta_s),E_B^{(r)}\right\rangle_{\mathrm{F}}+\nabla_{\mu^{(r)}}\Phi_u(\Theta_s)^\top E_\mu^{(r)}+\nabla_{\omega^{(r)}}\Phi_u(\Theta_s)^\top E_\omega^{(r)}\right]\,\mathrm{d}s. Φ u − Φ u = ∫ 0 1 r = 0 ∑ 1 [ ⟨ ∇ B ( r ) Φ u ( Θ s ) , E B ( r ) ⟩ F + ∇ μ ( r ) Φ u ( Θ s ) ⊤ E μ ( r ) + ∇ ω ( r ) Φ u ( Θ s ) ⊤ E ω ( r ) ] d s . Since E B ( r ) ∈ B E_B^{(r)}\in\mathcal{B} E B ( r ) ∈ B Π B \Pi_B Π B ℓ 1 / ℓ ∞ \ell_1/\ell_\infty ℓ 1 / ℓ ∞
For the stochastic part, Proposition C.5 gives the 2 a ∗ 2a_* 2 a ∗ q μ , r q_{\mu,r} q μ , r q ω , r q_{\omega,r} q ω , r H ϑ H_\vartheta H ϑ η B ( r ) = O P ( ρ r ) \eta_B^{(r)}=O_{\mathrm{P}}(\rho_r) η B ( r ) = O P ( ρ r ) ∥ Φ ^ − Φ ∥ ∞ = O P ( R d ) \|\widehat{\Phi}-\Phi\|_\infty=O_{\mathrm{P}}(R_d) ∥ Φ − Φ ∥ ∞ = O P ( R d ) R d R_d R d
For fold-specific games targeting the same population endpoints, applying the deterministic inequality to each fold and then the triangle inequality bounds the averaged attribution error by the average of the foldwise bounds. For a fixed number of folds, this preserves the displayed stochastic order. Each training fold must satisfy its score, RE/RSC, radius, and tolerance conditions at its own training-sample size; the folds themselves need not be independent. □ \square □
Proof of Corollary 4.6. Theorem 4.5 gives max u ∣ Φ ^ u − Φ u ∣ = O P ( R d ) \max_u|\widehat{\Phi}_u-\Phi_u|=O_{\mathrm{P}}(R_d) max u ∣ Φ u − Φ u ∣ = O P ( R d ) u ∈ S k u\in S_k u ∈ S k v ∉ S k v\notin S_k v ∈ / S k
∣ Φ ^ u ∣ − ∣ Φ ^ v ∣ ≥ ∣ Φ u ∣ − ∣ Φ v ∣ − ∣ Φ ^ u − Φ u ∣ − ∣ Φ ^ v − Φ v ∣ . |\widehat{\Phi}_u|-|\widehat{\Phi}_v|\ge|\Phi_u|-|\Phi_v|-|\widehat{\Phi}_u-\Phi_u|-|\widehat{\Phi}_v-\Phi_v|. ∣ Φ u ∣ − ∣ Φ v ∣ ≥ ∣ Φ u ∣ − ∣ Φ v ∣ − ∣ Φ u − Φ u ∣ − ∣ Φ v − Φ v ∣. Since every population boundary difference is at least g k g_k g k R d = o ( g k ) R_d=o(g_k) R d = o ( g k ) o P ( g k ) o_{\mathrm{P}}(g_k) o P ( g k ) k k k □ \square □
Evaluation Details and Reproducibility
This appendix documents the constructions, protocols, and diagnostics behind the main-paper example and experiments. Sections D.2–D.5 provide controlled evidence for resolution, computation, scarce-incident estimation, and the main modeling boundaries; Sections D.6–D.7 examine measured interventions. Unless stated otherwise, S 4 S^{4} S 4
Common protocol and metrics
For attribution magnitudes a a a G G G ∣ G ∣ = ∣ G ∩ Top ∣ G ∣ ( a ) ∣ / ∣ G ∣ |G| = |G \cap\operatorname{Top}_{|G|}(a)|/|G| ∣ G ∣ = ∣ G ∩ Top ∣ G ∣ ( a ) ∣/∣ G ∣ k k k ∣ Top k ( Φ ^ ) ∩ Top k ( Φ ) ∣ / k |\operatorname{Top}_{k}(\hat{\Phi}) \cap\operatorname{Top}_{k}(\Phi)|/k ∣ Top k ( Φ ^ ) ∩ Top k ( Φ ) ∣/ k F F F Q Q Q F + Q F + Q F + Q 2 d 2d 2 d Figure 2 evaluates the function-only ∣ Δ B ∣ |\Delta B| ∣Δ B ∣ F F F k k k
Function/law resolution and action ambiguity
Purpose and source SCMs. The construction isolates information lost by a paired-node report: whether the candidate response should restore a changed exogenous law or block an incoming mechanism. We generate 100 zero-mean, ordered linear-Gaussian source SCMs with d = 20 d = 20 d = 20 m m m R T R / m + Diag ( h ) R^{\mathsf{T}}R/m + \operatorname{Diag}(h) R T R / m + Diag ( h ) R ∈ R m × m R \in\mathbb{R}^{m \times m} R ∈ R m × m h j ∼ Unif ( 0.25 , 0.75 ) h_j \sim\operatorname{Unif}(0.25, 0.75) h j ∼ Unif ( 0.25 , 0.75 ) R T R / m R^{\mathsf{T}}R/m R T R / m h h h m m m Unif ( 0.70 , 1.30 ) \operatorname{Unif}(0.70, 1.30) Unif ( 0.70 , 1.30 ) Unif ( 0.70 , 1.30 ) \operatorname{Unif}(0.70, 1.30) Unif ( 0.70 , 1.30 ) Unif ( 0.04 , 0.12 ) \operatorname{Unif}(0.04, 0.12) Unif ( 0.04 , 0.12 )
Matched incident endpoints. Each source SCM is paired with two incident endpoints whose target attributions are ( Φ F , Φ Q ) = ( − β , + α ) (\Phi_F, \Phi_Q) = (-\beta, +\alpha) ( Φ F , Φ Q ) = ( − β , + α ) ( + α , − β ) (+\alpha, -\beta) ( + α , − β ) Φ Q = δ \Phi_Q = \delta Φ Q = δ α > β > δ > α − β > 0 \alpha> \beta> \delta> \alpha- \beta> 0 α > β > δ > α − β > 0 α − β \alpha- \beta α − β
To realize these identities, for each source SCM draw α ∼ Unif ( 1.5 , 2.5 ) \alpha\sim\operatorname{Unif}(1.5, 2.5) α ∼ Unif ( 1.5 , 2.5 ) β / α ∼ Unif ( 0.65 , 0.78 ) \beta/\alpha\sim\operatorname{Unif}(0.65, 0.78) β / α ∼ Unif ( 0.65 , 0.78 ) γ / α ∼ Unif ( 0.15 , 0.30 ) \gamma/\alpha\sim\operatorname{Unif}(0.15, 0.30) γ / α ∼ Unif ( 0.15 , 0.30 ) δ ∼ Unif ( α − β + 0.10 α , α − β + 0.20 α ) \delta\sim\operatorname{Unif}(\alpha- \beta+ 0.10\alpha, \alpha- \beta+ 0.20\alpha) δ ∼ Unif ( α − β + 0.10 α , α − β + 0.20 α ) β + γ \beta+ \gamma β + γ s t s_t s t s r s_r s r γ / ( β + γ ) \sqrt{\gamma/(\beta+ \gamma)} γ / ( β + γ ) α / s t \alpha/s_t α / s t δ / s r \delta/s_r δ / s r ( β + γ + α ) / ( β + γ ) \sqrt{(\beta+ \gamma+ \alpha)/(\beta+ \gamma)} ( β + γ + α ) / ( β + γ ) β / s t \beta/s_t β / s t β + Unif ( 0.40 α , 0.80 α ) \beta+ \operatorname{Unif}(0.40\alpha, 0.80\alpha) β + Unif ( 0.40 α , 0.80 α ) Unif ( 0.30 α , 0.60 α ) \operatorname{Unif}(0.30\alpha, 0.60\alpha) Unif ( 0.30 α , 0.60 α )
Candidate actions and evaluation. In the Q-actionable task, a rule restores either the target or distractor law to its source value. The resolved rule ranks the corresponding Q Q Q T ( M ( 1 ) ) − T ( M ( 1 ) , a ) T(M^{(1)}) - T(M^{(1),a}) T ( M ( 1 ) ) − T ( M ( 1 ) , a ) T ( M ( 1 ) , a ) − T ( M ( 0 ) ) T(M^{(1),a)} - T(M^{(0)}) T ( M ( 1 ) , a ) − T ( M ( 0 ) ) Table 2 shows that the resolved rule always selects the target law and clears the alarm, whereas the paired rule always selects the distractor and does not. Across both matched endpoints, the resolved action rule has balanced accuracy 1.0, whereas constant-Q, constant-F, and paired-only rules attain or are bounded by 0.5. This oracle construction demonstrates repair ambiguity under aggregation, not physical function/law identification or generic repair optimization.
Table 2. Q-actionable repair localization over 100 independent source problems (mean ± standard error). Score reduction is T ( M ( 1 ) ) − T ( M ( 1 ) , a ) T(M^{(1)}) - T(M^{(1),a)} T ( M ( 1 ) ) − T ( M ( 1 ) , a ) T ( M ( 1 ) , a ) − T ( M ( 0 ) ) T(M^{(1),a)} - T(M^{(0)}) T ( M ( 1 ) , a ) − T ( M ( 0 ) )
Exact computation: scalability and numerical agreement
Question and protocol. We ask whether the specialized algorithm removes the combinatorial Shapley bottleneck while preserving the same attribution. All methods target the same generated game; generic solvers receive its common coalition-value callable. The comparison therefore isolates Shapley aggregation rather than endpoint estimation. Stochastic solvers are included only to show the accuracy–runtime tradeoff. For each repetition and size, we generate a sparse ordered DAG with independent zero-mean disturbances, a diagonal quadratic score, and all K = 2 d K = 2d K = 2 d K ∈ { 10 , 20 , 40 , 100 , 200 , 400 , 1000 } K \in\{10, 20, 40, 100, 200, 400, 1000\} K ∈ { 10 , 20 , 40 , 100 , 200 , 400 , 1000 } K ∈ { 4 , 6 , 8 , 10 , 14 , 20 } K \in\{4, 6, 8, 10, 14, 20\} K ∈ { 4 , 6 , 8 , 10 , 14 , 20 } [3 ] ). Relative error is ∥ Φ ^ − Φ S 4 ∥ 2 / ∥ Φ S 4 ∥ 2 \lVert\widehat{\Phi} - \Phi^{S^4}\rVert_2/\lVert\Phi^{S^4}\rVert_2 ∥ Φ − Φ S 4 ∥ 2 / ∥ Φ S 4 ∥ 2 S 4 S^4 S 4
Results. Main Figure 1 shows that S 4 S^4 S 4 K = 1000 K = 1000 K = 1000 S 4 S^4 S 4 K = 20 K = 20 K = 20
Figure 1. Same-oracle Shapley aggregation on the full K = 2 d K = 2d K = 2 d ± \pm ± 4 ^{4} 4 K = 1000 K = 1000 K = 1000 ℓ 2 \ell_{2} ℓ 2 K = 20 K = 20 K = 20
Scarce-incident attribution recovery
Question and synthetic design. This benchmark asks whether reference anchoring recovers the leading function and law changes from scarce incident data. We generate 30 ordered Gaussian d = 100 d = 100 d = 100 n 1 ∈ { 20 , 40 , 60 , 80 , 100 } n_1 \in\{20, 40, 60, 80, 100\} n 1 ∈ { 20 , 40 , 60 , 80 , 100 } ω 0 ( j ) = Unif ( 0.8 , 1.2 ) 2 \omega^{(j)}_0 = \operatorname{Unif}(0.8, 1.2)^2 ω 0 ( j ) = Unif ( 0.8 , 1.2 ) 2 r j = min { p j , a j , max ( 1 , ⌈ 0.2 p j ⌉ ) } r_j = \min\{p_j, a_j, \max(1, \lceil0.2p_j\rceil)\} r j = min { p j , a j , max ( 1 , ⌈ 0.2 p j ⌉)} p j p_j p j a j a_j a j r j r_j r j r j r_j r j Q Q Q F + Q F+Q F + Q
Reference monitor and sample splitting. The reference size is n 0 = 10,000 n_0 = 10{,}000 n 0 = 10 , 000 μ ( 0 ) , ω ( 0 ) \mu^{(0)}, \omega^{(0)} μ ( 0 ) , ω ( 0 ) b b b
T 2 ( x b ) = ( x b − x ˉ b ) ⊤ { Σ ^ b + 10 − 3 tr ( Σ ^ b ) I / 50 } + ( x b − x ˉ b ) , T^2(x_b) = (x_b-\bar{x}_b)^\top\left\{\widehat{\Sigma}_b + 10^{-3}\operatorname{tr}(\widehat{\Sigma}_b)I/50\right\}^{+}(x_b-\bar{x}_b), T 2 ( x b ) = ( x b − x ˉ b ) ⊤ { Σ b + 1 0 − 3 tr ( Σ b ) I /50 } + ( x b − x ˉ b ) , where x ˉ b \bar{x}_b x ˉ b Σ ^ b \widehat{\Sigma}_b Σ b
For each incident sample size, five shuffled outer folds use 80% of the rows to fit B ℓ ( 1 ) B_\ell^{(1)} B ℓ ( 1 ) μ ℓ ( 1 ) , ω ℓ ( 1 ) \mu_\ell^{(1)}, \omega_\ell^{(1)} μ ℓ ( 1 ) , ω ℓ ( 1 )
Comparators. Resolved S4 ^4 4 4 ^4 4 ( F j , Q j ) (F_j,Q_j) ( F j , Q j ) ∣ Δ B ∣ \lvert\Delta B\rvert ∣ Δ B ∣ j j j ∥ B ^ ⋅ j ( 1 ) − B ^ ⋅ j ( 0 ) ∥ 2 \lVert\widehat{B}_{\cdot j}^{(1)}-\widehat{B}_{\cdot j}^{(0)}\rVert_2 ∥ B ⋅ j ( 1 ) − B ⋅ j ( 0 ) ∥ 2 F F F z z z j j j ∣ Z i j ( Z i W ) j ∣ \lvert Z_{ij}(Z_iW)_j\rvert ∣ Z ij ( Z i W ) j ∣ W W W [4 ] : it applies the fitted reference mechanism to both regimes and ranks D K L { N ( μ ^ j ( 1 ) , ω ^ j ( 1 ) ) ∥ N ( μ ^ j ( 0 ) , ω ^ j ( 0 ) ) } D_{\mathrm{KL}}\{\mathcal{N}(\widehat{\mu}_j^{(1)},\widehat{\omega}_j^{(1)})\mathbin{\|}\mathcal{N}(\widehat{\mu}_j^{(0)},\widehat{\omega}_j^{(0)})\} D KL { N ( μ j ( 1 ) , ω j ( 1 ) ) ∥ N ( μ j ( 0 ) , ω j ( 0 ) )}
The broader supplemental benchmark adds methods with adjacent but distinct outputs. Fixed-F F F [6 ] ). Linear-Gaussian CCA jointly switches fitted node conditionals using the incident-only endpoint, while MRCCA uses DoWhy 0.14’s multiply-robust score estimator with five folds, five Shapley permutations, a 2,000-row source cap, and propensity clipping at 0.02 ([4 , 25 ] ). We compute DCI’s estimated difference-DAG edge matrix Δ ^ \widehat{\Delta} Δ j ↦ ∥ Δ ^ ⋅ j ∥ 2 j \mapsto\lVert\widehat{\Delta}_{\cdot j}\rVert_2 j ↦ ∥ Δ ⋅ j ∥ 2 [33 ] ). We rank nodes by iSCAN’s returned variance-ratio shift statistic rather than thresholding its shifted-node set ([7 ] ). Our feature-shift SHAP baseline adapts SHAP by treating observed-coordinate distribution switches as players, without causal propagation ([16 ] ). FANS and do-Shapley remain conceptual comparisons in Table 1: the former does not return a scalar score attribution, and the latter targets variable interventions within one fixed SCM.
Protocol checks. Although the synthetic endpoints change only B B B ω \omega ω 10 − 10 10^{-10} 1 0 − 10 n 1 = 100 n_1=100 n 1 = 100
Results and interpretation. Main Figure 2 isolates the proposed method’s components: reference anchoring and score-aware propagation produce the clearest gains in function recovery, while the paired ablation shows that joint localization alone does not retain change type. Supplement Figure 1 places this result among broader alternatives. Fixed-F F F S 4 S^4 S 4
Figure 1 (PDF p. 30). Broad comparison on the same 30 scarce-incident SCMs as Main Figure 2 (mean ± SE). The main figure emphasizes transfer, resolution, and propagation ablations; this figure adds fixed-SCM event attribution, CCA/MRCCA, DCI, iSCAN, and feature-shift SHAP. Single-ranking outputs are copied to both F F F Q Q Q F F F n 1 = 100 n_1 = 100 n 1 = 100
Table 3 adds endpoint and full-vector error diagnostics. Recall improves with incident size even though full-vector error, especially for the law family, remains substantial. The evidence concerns leading effects, not uniform magnitude recovery or calibrated rank uncertainty.
Table 3. Calibrated d = 100 d = 100 d = 100 F F F Q Q Q F + Q F + Q F + Q ℓ 1 \ell_1 ℓ 1
Core-assumption stress tests
The guarantees assume a supplied common order, linear additive mechanisms, and a quadratic score. We vary these three conditions separately over ten independently generated d = 5 d = 5 d = 5 Oracle isolates the population attribution distortion caused by the tested departure, whereas E2E also estimates both SCM endpoints from samples. Relative ℓ 1 \ell_1 ℓ 1
All E2E runs use an 80/20 reference fit/moment split, reference Lasso and Delta-Lasso with five-fold internal CV, five-fold incident cross-fitting with held-out residual moments, and signed fold averaging. We set n 1 = 500 n_1 = 500 n 1 = 500 n 0 = 1000 n_0 = 1000 n 0 = 1000 n 0 n_0 n 0 Table 4 reports mean (SE) over SCMs.
Table 4. Core-assumption diagnostics over ten d = 5 d = 5 d = 5 ℓ 1 \ell_1 ℓ 1 ℓ 1 \ell_1 ℓ 1 R @ k R@k R @ k
Reference-only order estimation. We generate unique-order linear SCMs with random variable relabeling. The chain coefficients have independently random signs and magnitudes Unif ( 0.55 , 0.80 ) \operatorname{Unif}(0.55, 0.80) Unif ( 0.55 , 0.80 ) Unif ( 0.15 , 0.35 ) \operatorname{Unif}(0.15, 0.35) Unif ( 0.15 , 0.35 ) Unif ( 0.7 , 1.3 ) \operatorname{Unif}(0.7, 1.3) Unif ( 0.7 , 1.3 ) Unif ( 0.6 , 1.4 ) / d \operatorname{Unif}(0.6, 1.4)/d Unif ( 0.6 , 1.4 ) / d lingam 1.12.1 ([30 ] ) receives the 80% fit portion of nested total reference samples with n 0 ∈ { 100 , 200 , 300 , 400 } n_0 \in\{100, 200, 300, 400\} n 0 ∈ { 100 , 200 , 300 , 400 } Table 4 shows that improved order recovery removes the Oracle distortion, whereas E2E error remains governed largely by endpoint estimation. This is a favorable linear non-Gaussian identification check, not generic order-recovery evidence.
Nonlinear score approximation. Each zero-mean ordered linear-Gaussian SCM has positive chain coefficients B j − 1 , j ∼ Unif ( 0.20 , 0.38 ) B_{j-1,j} \sim\operatorname{Unif}(0.20, 0.38) B j − 1 , j ∼ Unif ( 0.20 , 0.38 ) Unif ( 0.08 , 0.20 ) \operatorname{Unif}(0.08, 0.20) Unif ( 0.08 , 0.20 ) Unif ( 0.75 , 1.25 ) \operatorname{Unif}(0.75, 1.25) Unif ( 0.75 , 1.25 ) d − 1 d - 1 d − 1 d d d B B B j ∈ [ d ] j \in[d] j ∈ [ d ]
s η ( x ) = x T W x + η ∑ j v j x j 4 . s_\eta(x) = x^\mathsf{T}Wx + \eta\sum_j v_j x_j^4. s η ( x ) = x T W x + η j ∑ v j x j 4 . Independently draw w ~ j , v ~ j ∼ Unif ( 0.7 , 1.3 ) \widetilde{w}_j,\widetilde{v}_j \sim\operatorname{Unif}(0.7, 1.3) w j , v j ∼ Unif ( 0.7 , 1.3 ) W = Diag ( w ~ ) / d W = \operatorname{Diag}(\widetilde{w})/d W = Diag ( w ) / d v j = v ~ j / d > 0 v_j = \widetilde{v}_j/d > 0 v j = v j / d > 0 Σ 0 \Sigma_0 Σ 0 q 0 = tr ( W Σ 0 ) q_0 = \operatorname{tr}(W\Sigma_0) q 0 = tr ( W Σ 0 ) f 0 = 3 ∑ j = 1 d v j Σ 0 , j j 2 f_0 = 3\sum_{j=1}^{d} v_j\Sigma_{0,jj}^2 f 0 = 3 ∑ j = 1 d v j Σ 0 , j j 2 η = r q 0 / f 0 \eta= rq_0/f_0 η = r q 0 / f 0 r ∈ { 0.2 , 0.4 , 0.6 } r \in\{0.2, 0.4, 0.6\} r ∈ { 0.2 , 0.4 , 0.6 } r r r x T W x x^\mathsf{T}Wx x T W x 2 2 d 2^{2d} 2 2 d
Nonlinear mechanism stress. We next use
x j = ∑ i < j B i j x i + ∑ i < j G i j tanh ( x i ) + e j x_j = \sum_{i<j} B_{ij}x_i + \sum_{i<j} G_{ij}\tanh(x_i) + e_j x j = i < j ∑ B ij x i + i < j ∑ G ij tanh ( x i ) + e j on a chain-plus-extra-edge graph: chain B B B Unif ( 0.22 , 0.38 ) \operatorname{Unif}(0.22, 0.38) Unif ( 0.22 , 0.38 ) G G G Unif ( 0.25 , 0.45 ) \operatorname{Unif}(0.25, 0.45) Unif ( 0.25 , 0.45 ) B i j ∼ Unif ( − 0.16 , 0.16 ) B_{ij} \sim\operatorname{Unif}(-0.16, 0.16) B ij ∼ Unif ( − 0.16 , 0.16 ) G i j ∼ Unif ( − 0.25 , 0.25 ) G_{ij} \sim\operatorname{Unif}(-0.25, 0.25) G ij ∼ Unif ( − 0.25 , 0.25 ) c ∈ { 0.3 , 0.6 , 0.9 } c \in\{0.3, 0.6, 0.9\} c ∈ { 0.3 , 0.6 , 0.9 } G ( 0 ) = c G u n i t G^{(0)} = cG^{\mathrm{unit}} G ( 0 ) = c G unit d − 1 , d d - 1,d d − 1 , d B B B G G G Unif ( 0.70 , 1.10 ) \operatorname{Unif}(0.70, 1.10) Unif ( 0.70 , 1.10 ) d d d
Zero-departure controls agree numerically, and the higher-order Gauss–Hermite audit leaves the nonlinear-mechanism conclusion unchanged. With ten SCMs, these tests expose qualitative boundaries rather than support a general robustness claim. Latent confounding, cycles, dependent disturbances, and general nonlinear or nonadditive identification remain outside scope.
Causal Chamber: measured Light-Tunnel interventions
Question, data, and score. We ask whether a resolved report can separate the manipulated target-law component of a fixed all-channel alarm from propagated fitted-response components. We use 16 Light-Tunnel set-point interventions whose manipulated variable is among the 22 measured channels [9 ] . Disjoint reference subsets fit a 12-component PCA T 2 + S P E T^2 + \mathrm{SPE} T 2 + SPE T 2 / s d c a l ( T 2 ) + S P E / s d c a l ( S P E ) T^2 / \mathrm{sd}_{\mathrm{cal}}(T^2) + \mathrm{SPE} / \mathrm{sd}_{\mathrm{cal}}(\mathrm{SPE}) T 2 / sd cal ( T 2 ) + SPE / sd cal ( SPE )
Methods, metrics, and checks. S⁴-RCA ranks the flat report by a j = ∣ Φ ‾ F j ∣ + ∣ Φ ‾ Q j ∣ a_j = \lvert\overline{\Phi}_{F_j}\rvert+ \lvert\overline{\Phi}_{Q_j}\rvert a j = ∣ Φ F j ∣ + ∣ Φ Q j ∣ F F F B B B Q Q Q
Localization and resolved allocation. Main Figure 3 reports target ranks and the harder polarization allocations. Across all 16 interventions, S⁴-RCA ranks the target node and its Q j Q_j Q j F F F
Optical interpretation and boundary. The supplied topology places the ir_3 infrared sensor downstream of both polarizer controls. Their measured relative set-point angle moves from near parallel to about ± 40 ∘ \pm40^\circ ± 4 0 ∘ ± 79 ∘ \pm79^\circ ± 7 9 ∘ cos 2 ( θ 1 − θ 2 ) \cos^2(\theta_1-\theta_2) cos 2 ( θ 1 − θ 2 ) ir_3, but does not establish that the optical law itself changed. The resolved allocation in Main Figure 3 and the projection audit below are consistent with this relinearization interpretation.
Projection audit. Table 5 replaces the dense graph display with the coefficients relevant to the optical interpretation. The supplied topology is fixed, and coefficients are all-row OLS projections in reference-standardized coordinates. Mid-regime signs mirror across the manipulated polarizer, while strong-regime R 2 R^2 R 2
Table 5. Reference-standardized polarization projections on the fixed official subgraph. P1/P2 are pol_1/pol_2; coefficients and R 2 R^2 R 2
External intervention stress tests: Wind Tunnel and Sachs
Purpose, methods, and metric. We use Wind Tunnel to test whether the all-channel operational-score analysis extends to a second physical system, and Sachs to examine the same pipeline under substantially weaker model assumptions. We compare S⁴-RCA, Fixed-F F F B B B k k k k k k
Table 6. All-channel external interventions. H@k k k k k k
Wind Tunnel. We retain all 26 nonempty 10,000-row single-target regimes in wt_walks_v1, covering 11 targets and 26 nonconstant channels after excluding six constants [9 ] . Disjoint reference subsets fit/calibrate an all-channel PCA T 2 + S P E T^2 + \mathrm{SPE} T 2 + SPE
Table 6 shows that direct endpoint and residual-KL detectors localize these single-target interventions at least as well as the score-attribution methods. S⁴-RCA nevertheless places every target in the top three and assigns little mass to F F F Q Q Q
Sachs. We log-transform 11 phosphoproteins and use condition 1 as reference; the seven tests are AKT inhibitor, G06976, Psitectorigenin, U0126, LY294002, PMA, and the dataset label b2cAMP [27 ] . Labels are documented measured targets, except that measured Akt serves as a proxy for unmeasured PI3K in LY294002, and, for U0126, either measured MEK, the direct intervention target, or downstream ERK counts as a hit. The 17-edge consensus DAG supplies support/order. Disjoint reference subsets of 200/150/250/126/127 cells serve PCA fitting, score calibration, SCM fitting, residual-moment estimation, and monitor evaluation, respectively; the monitor has eight components, and 400/200 incident cells are used for endpoint updates and evaluation. Sachs mean/minimum monitor AUROC is 0.844/0.630, and only four conditions reach 0.75. No method in Table 6 localizes more than three of the seven targets within the top three. We retain all seven in the common denominator; the weak monitor and localization make Sachs an assumption-boundary stress test, not positive cross-domain validation.
Interpretation boundary. Light and Wind allocate a fixed all-channel score under a supplied ordered linear reapproximation; endpoint shift, Residual local KL, and Observed PCA do not share the resolved function/law estimand. Sachs additionally involves proxy targets, latent confounding, and nonlinear biology. These studies do not establish hidden-target localization, causal discovery, physical mechanism recovery, universal superiority, or calibrated rank uncertainty.
Reproducibility
Software and traceability. Experiments use Python 3.12.4 with NumPy 2.3.5, SciPy 1.15.3, scikit-learn 1.8.0, pandas 2.3.3, and Matplotlib 3.10.8; the solver audit additionally uses DoWhy 0.14, the extended benchmark uses causaldag 0.1a163 and iscan-dag 1.0.0, and the order diagnostic uses lingam 1.12.1. Internal workflows retain raw and aggregated outputs and run manifests linking commands, settings, seeds, and paper artifacts. The principal synthetic seeds are 270730 for the action-ambiguity construction, 270713 for computation, 20270629 for scarce-incident recovery, and 20260728 for the assumption diagnostics. Light-Tunnel splits use seeds 20260715–20260719; Wind-Tunnel and Sachs splits use 20260720–20260724.
Penalty selection. Nodewise LassoCV RMS-scales predictors within the training fold, fits an intercept, and minimizes mean validation error over 100 log-spaced penalties from α max = ∥ X c T y c ∥ ∞ / n \alpha_{\max} = \lVert X_c^{\mathsf{T}} y_c \rVert_{\infty}/n α m a x = ∥ X c T y c ∥ ∞ / n 10 − 3 α max 10^{-3}\alpha_{\max} 1 0 − 3 α m a x max_iter = 100,000, and tolerance 10 − 4 10^{-4} 1 0 − 4
AI disclosure. Generative AI assisted manuscript editing, code review, and experiment orchestration. It was neither treated as an author nor used as a cited source; the human authors retain full responsibility.
