Dec 2024, revised Jan 2026 · with Yifan Hu (equal contribution), Peter Bühlmann, Zijian Guo · arXiv:2412.11850
The short version
A company wants to know which ad channel causally drives sales \(Y\): clicks from search ads \(X_1\), or engagement from social ads \(X_2\)? A model fit purely to predict \(Y\) will happily use both, because both correlate with sales. But correlation cannot separate the mechanisms below, and spending against a non-causal metric wastes money.
A powerful principle rescues us: if a set of variables truly causes \(Y\), the way \(Y\) depends on them stays invariant across environments — time slots, regions, regimes — even as everything else shifts. Environments are modelled as additive interventions on the covariates: a linear structural equation model whose noise gains an environment-specific shift \(\delta^e\) on the \(X\) coordinates, never on \(Y\)'s own noise.
This paper uses the weakest useful form, risk invariance: only the residual second moment must match across environments. Writing \(R_e(b) = \mathbb E[(Y^e - b^\top X^e)^2]\), the invariant set is
\[ \mathcal B_{\mathrm{inv}} \;=\; \big\{\, b \in \mathbb R^p \ :\ R_e(b) = R_f(b) \quad \text{for all } e, f \in \mathcal E \,\big\}. \]This matters more than it looks. The conditional-mean version of invariance used by EILLS, and the stronger independence-based version underlying ICP — which asks that \(\varepsilon_Y\) be independent of \(X_{S^\star}\) with an identical distribution across environments — both break under hidden confounding. Risk invariance does not: the systematic noise \(\eta\) may correlate \(Y\) and \(X\) freely, and \(\beta^\star\) remains the unique risk-invariant model. Theorem 1 shows \(\mathcal B_{\mathrm{inv}} = \{\beta^\star\}\) under an identification condition, and Theorem 2 shows that condition is nearly necessary.
The identification condition (Condition 1b) is refreshingly concrete — it asks only that the environments differ enough in a spectral sense. With \(\boldsymbol A(w) = \sum_e \big(w_e - \tfrac{1}{|\mathcal E|}\big)\,\mathbb E[X^e X^{e\top}]\), it requires some simplex weight \(w^0\) with
\[ \lambda \;:=\; \lambda_{\min}\big(\boldsymbol A(w^0)\big) \;>\; 0 . \]The scalar \(\lambda\) is a measure of environmental heterogeneity, and it governs every rate below. Crucially, this can hold with as few as two environments — where IRM's guarantee wants \(|\mathcal E| \ge p\), and earlier invariance results wanted infinitely many.
Invariant causal prediction turns the principle into a procedure: for every subset \(S \subseteq \{1,\dots,p\}\), test whether the relationship is invariant, and intersect the sets that pass. It is correct, and it is a search over \(2^p\) subsets. EILLS replaces the tests with a regularized least-squares objective, but that objective is discontinuous and nonconvex, and solving it still amounts to enumeration — essentially a mixed-integer program. Its continuous Gumbel-softmax relaxation buys speed by giving up the identification guarantee. So the field had a choice between a method with theory and no runtime, and a method with runtime and no theory.
Start from the constrained program "minimize risk subject to risk invariance", and penalize the (nonconvex) equality constraint. What comes out is a minimax problem that looks like distributionally robust optimization, except that the weights on environments may be negative:
\[ b^\gamma_{\mathrm{Neg}} \;\in\; \arg\min_{b\in\mathbb R^p}\ \max_{w \in \mathcal U(\gamma)}\ \sum_{e\in\mathcal E} w_e\, R_e(b), \qquad \mathcal U(\gamma) = \Big\{ w : \textstyle\sum_e w_e = 1,\ \min_e w_e \ge -\gamma \Big\}. \]Set \(\gamma = 0\) and \(\mathcal U(0)\) is the simplex: this is exactly Group DRO, which merely protects the worst environment. Let \(\gamma > 0\) and something different happens. The inner maximum is attained at a vertex that loads \(1 + \gamma(|\mathcal E|-1)\) on the worst environment and \(-\gamma\) on every other, so the objective has the closed form
\[ \Phi_\gamma(b) \;=\; \big(1 + \gamma|\mathcal E|\big)\,\max_{e}\, R_e(b) \;-\; \gamma \sum_{e\in\mathcal E} R_e(b). \]Read it as: minimize the worst risk, while being rewarded for total risk. The only way to win on both counts is to pull the maximum risk down toward the average — that is, to equalize risks. As \(\gamma \to \infty\) the penalty enforces exact risk invariance and \(b^\gamma_{\mathrm{Neg}} \to \beta^\star\); at finite \(\gamma\), Proposition 2 gives \(\|b^\gamma_{\mathrm{Neg}} - \beta^\star\|_2 \lesssim \big[\lambda(1+\gamma|\mathcal E|)\big]^{-1}\). Robustness has become identification.
Nothing is free: negative weights make \(\Phi_\gamma\) a difference of convex functions, hence nonconvex. But unlike ICP and EILLS, it is a continuous problem in \(b\) — no subsets, no integers. In practice one adds a ridge term \(-\mu\|w\|_2^2\) inside the maximum, which makes the inner problem strongly concave and \(\Phi_\mu\) differentiable by Danskin's theorem, then alternates a closed-form simplex update for \(w\) with a gradient step on \(b\).
Solving a nonconvex problem to global optimality is generally hopeless, so a guarantee about the global minimizer would be of little use. The paper proves something stronger and stranger — a benign landscape. Under the same heterogeneity condition, for any point \(b\) whatsoever,
\[ \|b - \beta^\star\|_2 \;\lesssim\; \frac{1}{\lambda}\Big( \frac{1}{1+\gamma|\mathcal E|} \;+\; \|\nabla \Phi_\mu(b)\|_2 \Big) \;+\; \sqrt{\mu/\lambda}. \]Read the right-hand side: it is small whenever the gradient is small. Every stationary point — not just the global minimum, not just local minima — sits near \(\beta^\star\). Saddle points and spurious local minima, if they exist, are all in the right neighbourhood. This is what licenses a plain first-order method, which is only ever guaranteed to find stationary points. Note the bound is on the distance to \(\beta^\star\), not to the optimum: it is not a classical error bound, and it is why global optimization is unnecessary rather than merely unattractive.
Feeding in sampling error and a finite iteration budget, Theorem 5 gives the non-asymptotic rate. With \(T\) iterations, \(n\) samples per environment, and the optimal \(\mu \asymp T^{-1/2}\),
\[ \|\widehat b^\gamma - \beta^\star\|_2 \;\lesssim\; \underbrace{\frac{1}{1+\gamma|\mathcal E|}}_{\text{penalty}} \;+\; \underbrace{T^{-1/4}}_{\text{optimization}} \;+\; \underbrace{n^{-1/4}}_{\text{statistics}} . \]So \(\epsilon\)-accuracy needs \(\gamma = \Omega(\epsilon^{-1})\), \(T = \Omega(\epsilon^{-4})\) and \(n = \Omega(\epsilon^{-4})\) — all polynomial, with each iteration costing one gradient evaluation over \(p\) coordinates and \(|\mathcal E|\) risks. Whether the \(n^{-1/4}\) and \(T^{-1/4}\) can be improved to \(n^{-1/2}\), \(T^{-1/2}\) is stated as open.
The evidence is simulation — the paper has no real-data application. The headline experiment sweeps \(p\) from 5 to 100 with \(n = 20{,}000\), four environments, and a 30-minute cap per run. NegDRO and EILLS both use their recommended \(\gamma = 20\); the remaining causal-invariance baselines are tuned in an oracle manner, choosing whichever hyperparameter minimizes \(\|\widehat b - \beta^\star\|_2\) — a handicap in NegDRO's disfavour.
| Method | Confounding? | Environments | Search | Largest \(p\) in 30 min |
|---|---|---|---|---|
| NegDRO (ours) | yes | as few as 2 | continuous | 100+ |
| ICP | no | ≥ 2 | all \(2^p\) subsets | 10 |
| EILLS | no | ≥ 2 | mixed-integer | 20 |
| IRM | — | \(|\mathcal E| \ge p\) | continuous | — |
| ERM / Anchor | — | — | convex | 100+ |
"Confounding?" asks whether \(\beta^\star\) is still identified when the systematic noise correlates \(Y\) with \(X\); ICP's independence-based invariance and EILLS's conditional-mean invariance both fail there. In the paper's confounded simulation, NegDRO's error stays at 0.03–0.10 across \(p\), while ICP and EILLS both sit near 0.7. ERM and Anchor regression scale to any \(p\) but never identify \(\beta^\star\) at all — they are fast at solving the wrong problem.
A second experiment probes the regime the theory cares about: how strong must the interventions be? With two environments, four covariates, and interventions of varying strength, the competitors that need a dominated reference environment (DRIG) or an invertible difference of second moments (CausalDantzig) degrade badly, while NegDRO is nearly unaffected.
| Estimation error, 2 environments | Limited interventions | Weak interventions | Strong interventions |
|---|---|---|---|
| NegDRO (ours) | ≈ 0.05 | ≈ 0.05 | ≈ 0.05 |
| CausalDantzig | ≈ 0.75 | ≈ 0.75 | ≈ 0.05 |
| DRIG | grows to ≈ 2.5 | grows to ≈ 2.7 | ≈ 0.05 |
Values read from the paper's Figure 5 at large \(\gamma\), \(n_e = 10{,}000\). "Limited" means some covariates are never intervened on; "weak" means the intervention variance is 0.01. DRIG's error increases with \(\gamma\) in these regimes — its guarantee needs a strictly dominated environment that does not exist here. Section 5 extends NegDRO's identification to the limited regime under a weaker condition: it suffices to intervene on the children of \(Y\).
This closes a gap between an elegant statistical principle and usable methodology. Invariance-based discovery had the right theory but an exponential price tag; recasting it as a continuous nonconvex program — and then proving the landscape is benign — makes it feasible at realistic dimensions, with the multi-environment data modern applications already generate: multiple hospitals, time periods, regional markets, experimental regimes. The conceptual takeaway is that a negative weight is not a bug in the DRO formalism but the precise device that converts worst-case robustness into causal identification. To our knowledge this is the first causal-invariance method that provably reaches approximate global optimality for its nonconvex objective in polynomial time.
Honest caveats. The theory is for linear structural equation models, and the authors position it as a stepping stone toward nonlinear invariance learning. The \(n^{-1/4}\) and \(T^{-1/4}\) rates may well be suboptimal; closing them is open. Risk invariance hinges on the outcome's own noise never being intervened on, so all heterogeneity must come from additive interventions on the covariates. Identification conditions are sufficient and only nearly necessary — necessity is proven for single-coordinate interventions. And all evidence is simulated: there is no real-data application in the paper.
© 2026 Zhenyu (Zach) Wang 王振宇