Paper deep dive
Structure-agnostic Causal Representation Learning
Arman Behnam, Binghui Wang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 10/3/2026, 1:50:28 AM
Summary
The paper introduces SaCRL (Structure-agnostic Causal Representation Learning), a framework that jointly identifies the underlying causal structure of data and learns invariant representations without prior structural knowledge. It uses HSIC-based violation metrics and a SoftMin optimization to automatically select the correct invariance constraint among three canonical causal structures (anti-causal, confounded-descendant, confounded-outcome). The method provides theoretical guarantees for structure identification and out-of-distribution generalization, demonstrating state-of-the-art performance on benchmarks like PACS, VLCS, and OfficeHome.
Entities (11)
Relation Signals (8)
SaCRL → achievesstateofthearton → VLCS
confidence 95% · achieves state-of-the-art accuracy on three DomainBed benchmarks (PACS, VLCS, OfficeHome)
SaCRL → achievesstateofthearton → OfficeHome
confidence 95% · achieves state-of-the-art accuracy on three DomainBed benchmarks (PACS, VLCS, OfficeHome)
SaCRL → achievesstateofthearton → PACS
confidence 95% · achieves state-of-the-art accuracy on three DomainBed benchmarks (PACS, VLCS, OfficeHome)
SaCRL → uses → HSIC
confidence 95% · Our approach formulates structure selection as a soft optimization over candidate invariances using HSIC-based violation metrics
SaCRL → employs → SoftMin
confidence 90% · This asymmetry enables structure identification via a SoftMin objective
SaCRL → generalizes → ACIA
confidence 90% · SaCRL generalizes Anti-Causal Invariant Abstractions (ACIA) [7], which learns invariant representations for the anti-causal structure
SaCRL → outperforms → Colored MNIST
confidence 90% · outperforms fixed-invariance baselines on Colored MNIST
Anti-causal → istargetedby → ACIA
confidence 85% · ACIA [7], which learns invariant representations for the anti-causal structure
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Causal representation learning aims to discover robust features by exploiting the causal structure underlying data generation. Existing methods require specifying the causal structure a priori, yet different structures demand fundamentally incompatible invariance constraints, and misspecification leads to representations that discard predictive information. We introduce SaCRL, a framework that jointly identifies the causal structure and learns the corresponding invariant representation without prior structural knowledge. Our approach formulates structure selection as a soft optimization over candidate invariances using HSIC-based violation metrics, with adaptive weights that automatically concentrate on the achievable structure. We provide theoretical guarantees for structure identification, including under random-feature approximation, invariance satisfaction, and out-of-distribution generalization. Empirically, SaCRL recovers the true structure on synthetic and semi-synthetic Bayesian-network benchmarks, outperforms fixed-invariance baselines on Colored MNIST, achieves state-of-the-art accuracy on three DomainBed benchmarks (PACS, VLCS, OfficeHome), and degrades gracefully under structural misspecification and limited environment diversity. Code is available at: this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2610.00968v1
- Canonical: https://arxiv.org/abs/2610.00968v1
Trouble viewing inline? Open PDF directly →
Full Text
178,511 characters extracted from source content.
Expand or collapse full text
Structure-agnostic Causal Representation Learning Arman Behnam Affiliation: Department of Computer Science Affiliation: Illinois Institute of Technology Affiliation: Chicago, IL, USA Email: abehnam@hawk.illinoistech.edu Binghui Wang Affiliation: Department of Computer Science Affiliation: Illinois Institute of Technology Affiliation: Chicago, IL, USA Email: bwang70@illinoistech.edu Abstract Causal representation learning aims to discover robust features by exploiting the causal structure underlying data generation. Existing methods require specifying the causal structure a priori, yet different structures demand fundamentally incompatible invariance constraints, and misspecification leads to representations that discard predictive information. We introduce SaCRL, a framework that jointly identifies the causal structure and learns the corresponding invariant representation without prior structural knowledge. Our approach formulates structure selection as a soft optimization over candidate invariances using HSIC-based violation metrics, with adaptive weights that automatically concentrate on the achievable structure. We provide theoretical guarantees for structure identification, including under random-feature approximation, invariance satisfaction, and out-of-distribution generalization. Empirically, SaCRL recovers the true structure on synthetic and semi-synthetic Bayesian-network benchmarks, outperforms fixed-invariance baselines on Colored MNIST, achieves state-of-the-art accuracy on three DomainBed benchmarks (PACS, VLCS, OfficeHome), and degrades gracefully under structural misspecification and limited environment diversity. Code is available at: https://github.com/ArmanBehnam/sacrl. 1 Introduction Causal representation learning [55] has emerged as a principled approach to building robust models by uncovering the underlying mechanisms that generate data, enabling generalization beyond training distributions [55, 45]. Current approaches are broadly categorized into two paradigms, each with distinct structural assumptions. Structure-based methods [73, 61, 6] require explicit knowledge of Directed Acyclic Graphs (DAGs) and complete Structural Causal Models (SCMs), imposing significant barriers when causal structures are unknown. Distribution-invariant methods [4, 30] avoid an explicit graph but enforce a fixed invariance across environments, committing to one causal structure. When that commitment is wrong, the enforced invariance discards predictive signal. Multi-HospitalClinical Data=Deii=13D=\D_e_i\_i=1^3e1e_1: Generale2e_2: ICUe3e_3: RuralUnknown Structurefor Each VariableStructure-AgnosticCausal Representation LearningDiscoveredStructuresSymptoms: 1G_1Biomarker: 2G_2Vitals: 3G_3OptimalRepresentation ϕ∗φ^*Prediction on New Hospital eteste_testLow Test Error Figure 1: Clinical toy example. Consider predicting patient outcomes from clinical data collected at several hospitals (Figure 1). Different variables follow different causal patterns: symptoms arise anti-causally from the disease (1G_1); biomarkers are confounded descendants that share a genetic cause with the diagnosis (2G_2); and vital signs are confounded outcomes, responding to the disease and to unobserved severity (3G_3). The structure for each variable is unknown to the learner. Existing methods require manually partition variables by (a priori) causal type, a process that is error-prone (misspecification yields suboptimal representations [78, 32]), computationally prohibitive (O(mk)O(m^k) evaluations for m structures over k variable subsets [88, 66]), and theoretically unsound (no guarantees exist when structure identification introduces uncertainty [81, 14]). Struc.-Specific (Correct)UUXXYYTrue: Conf-DescApply R2R_G_2:P(Y|ϕ(X))P(Y|φ(X)) invariantGood OODTest Error: LowStruc.-Specific (Wrong)UUXXYYTrue: Conf-DescApply R1R_G_1:ϕ(X)⟂E|Yφ(X) E|YPoor OODTest Error: HighStruc.-Agnostic (Ours)XXYYUnknown StructureLearn α:Auto-select RkR_G_kGood OODRobust to Structure Figure 2: Left: When true structure (confounded-descendant) is known, applying the correct regularizer yields good Out-of-Distribution performance. Middle: Misspecifying structure leads to wrong invariance constraints and poor generalization. Right: Our structure-agnostic approach automatically discovers the structure and applies appropriate regularization. XzX_zYYZZXz⟂X_z UUEE X (a) 1G_1 XzX_zYYZZXz⟂X_z UUEE X (c) 2G_2 XzX_zYYZZXz⟂X_z UUEE X (b) 3G_3 Figure 3: Three canonical causal structures: (a) anti-causal 1G_1, (b) confounded-descendant 2G_2 (c) confounded-outcome 3G_3. Latent variables U,ZU,Z (shaded dashed nodes) influence the relationship between observed features X (dotted box) and target Y, with environment E capturing distributional shifts. Magenta arrows denote causal mechanisms that distinguish each structure: Y→XY→ X for anti-causal, X→YX→ Y for both confounded types. Dashed arrows represent confounding paths through U and Z. Different causal structures demand incompatible invariance properties [40, 74]. A representation optimized for one structure performs poorly under another, as shown in Figure 3. The central question we address is can we learn optimal invariant representations without knowing the causal structure? This requires a paradigm shift from structure-specific to structure-agnostic learning11 1 It refers to the learner’s prior knowledge, not the space of all possible structures. Prior methods require specifying the true structure that governs the data in advance, whereas we eliminate this requirement by learning the structure directly from data.: jointly discover causal mechanisms generating the data and learn representations tailored to those mechanisms. The challenge is both computational (avoiding exponential structure search) and theoretical (providing guarantees when structure itself is uncertain). We introduce a kernel-based framework, Structure-agnostic Causal Representation Learning (SaCRL), for joint structure discovery and representation learning. SaCRL generalizes Anti-Causal Invariant Abstractions (ACIA) [7], which learns invariant representations for the anti-causal structure 1G_1, to settings where the structure itself is unknown. The key theoretical insight is that under the true causal structure, the corresponding invariance is achievable by predictive representations, whereas wrong invariances incur irreducible violation. This asymmetry enables structure identification via a SoftMin objective that automatically concentrates on the correct structure. We prove: (i) consistent structure identification, which also holds under random-feature approximation of HSIC, (i) bounds on the remaining invariance violation, (i) out-of-distribution generalization bounds, and (iv) sample complexity for joint identification and generalization. Empirically, SaCRL also recovers the true structure on semi-synthetic Bayesian-network data with known ground truth, and its adaptive weights provide a calibrated confidence for the selected structure. 2 Background Random variables are uppercase (X,Y,EX,Y,E), realizations lowercase (x,y,ex,y,e), and spaces calligraphic (,,ℰX,Y,E). We write |ℰ||E| for set cardinality, P[⋅]E_P[·] for expectation under P, and AI_A for the indicator of set A. Statistical independence is denoted X⟂YX \!\!\! Y (unconditional) and X⟂Y|ZX \!\!\! Y Z (conditional). Complete notations appear in Appendix B. 2.1 Measure Theory Following ACIA [7], this work uses measure-theoretic probability. See Appendix C.1 for formal definitions of measurable space, probability measure, product space, σ-algebra, and Markov kernel. 2.2 Canonical Causal Structures We formalize causal relationships using DAGs in Appendix C.2, which encode the qualitative structure of causal mechanisms. There exist three canonical causal structures [74], with definitions below: Definition 2.1. Three Canonical Causal Structures. (a) Anti-causal (1G_1): DAG contains edges Y→X←EY→ X← E, capturing settings where labels generate observations (e.g., disease causes symptoms). (b) Confounded-descendant (2G_2): Observed features decompose as X=(Xz,Xz⟂)X=(X_z,X_z ) where Xz⟂→YX_z → Y, and a latent confounder U influences both Xz⟂X_z and other components. DAG includes Xz⟂→YX_z → Y with U→Xz⟂U→ X_z . (c) Confounded-outcome 3G_3: Observed features decompose as X=(Xz,Xz⟂)X=(X_z,X_z ) where Xz⟂→Y→XzX_z → Y→ X_z, and a latent confounder U influences both Y and the latent variable Z (which in turn affects XzX_z). DAG includes Xz⟂→Y→XzX_z → Y→ X_z with U→YU→ Y and U→Z→XzU→ Z→ X_z. Here, XZ⟂X_Z denotes the component of X that is causally orthogonal to the Z-mediated pathway. Specifically, XZX_Z captures variation in X that flows through the latent Z (e.g., features influenced by Y via Z in the anti-causal structure), and XZ⟂X_Z captures variation from other sources, including direct effects from Y, confounders U, or environment E. Figure 3 illustrates the three canonical structures. The three canonical structures form an exhaustive partition of valid DAGs under two mild conditions: (i) partial observability: the learner observes only (X,Y,E)(X,Y,E) with latent confounders U unobserved; (i) environment mechanism: E influences the system through latent confounders rather than directly rewiring the X–Y relationship. Under (i)–(i), the causal relationship between X and Y must take one of three forms: (1) the X–Y edge is either Y→XY→ X (anti-causal 1G_1) or X→YX→ Y (causal). In the causal direction, confounding either operates through features (2G_2) or bypasses the label via a latent mediator (3G_3). Bidirectional X↔YX Y violates acyclicity; unconfounded X→YX→ Y is recovered by 2G_2 with zero confounder influence. These conditions are standard in causal representation learning [74]. Our theoretical analysis relies on the distinction between observational conditioning P(Y∣X=x)P(Y X=x) and interventional conditioning P(Y∣do(X=x))P(Y (X=x)), where the latter severs incoming edges to X in the causal graph. While this distinction motivates why different causal structures require different invariance properties, our method circumvents explicit do-calculus. Formal definitions of interventions and causal independence appear in Appendix C.3. Definition 2.2 (Environment Space). An environment space consists of a finite collection ℰ=e1,…,emE=\e_1,…,e_m\ of environments, where each environment ei∈ℰe_i has a measurable input space (ei,ℱei)(X_e_i, F_X_e_i), a measurable output space (ei,ℱei)(Y_e_i, F_Y_e_i), a probability measure PeiP_e_i on the product space (ei×ei,ℱei⊗ℱei)(X_e_i×Y_e_i, F_X_e_i F_Y_e_i), and a dataset Dei=(xjei,yjei)j=1neiD_e_i=\(x_j^e_i,y_j^e_i)\_j=1^n_e_i sampled i.i.d. from PeiP_e_i. Each environment e induces a joint distribution Pe(X,Y)P_e(X,Y) that factorizes according to an underlying causal DAG eG_e. We formalize this via causal spaces that pair probability measures with their generating mechanisms. The key quantity is the interventional kernel distance d(e,e′)d_K(e,e ) measuring how causal mechanisms differ across environments. Full measure-theoretic causality definitions are in Appendix C.4. One should understand the meaning of representation for this purpose. Definition 2.3 (Representation). It is a measurable function ϕ:→φ:X mapping input space X to a latent representation space Z, where (,ℱ)(Z, F_Z) is a measurable space. 3 Problem Formulation 3.1 Motivation Different causal structures demand fundamentally incompatible invariance properties. An anti-causal setting (Y→XY→ X) requires representations conditionally independent of environments given the target. A confounded-descendant setting (X→YX→ Y with confounding) requires the predictive distribution P(Y|ϕ(X))P(Y|φ(X)) to be invariant. A confounded-outcome setting requires the marginal P(ϕ(X))P(φ(X)) to be invariant. A representation optimized for one structure typically performs poorly under another. Suppose the true structure is confounded-descendant (2G_2), but we enforce anti-causal invariance (ℐ1:ϕ(X)⟂E|YI_1:φ(X) E|Y). This constraint forces the representation to discard information that varies with E conditional on Y, but in 2G_2, this includes precisely the causal signal X→YX→ Y. The learned representation becomes predictively useless. We observe that the correct invariance is achievable (the representation can satisfy it while maintaining prediction accuracy), whereas incorrect invariances are unachievable (satisfying them requires discarding predictive information). 3.2 Problem Setup Let ℰ=e1,…,emE=\e_1,…,e_m\ with m≥2m≥ 2 denote the set of environments. For each environment e∈ℰe , we observe a dataset De=(xje,yje)j=1neD_e=\(x_j^e,y_j^e)\_j=1^n_e drawn i.i.d. from an unknown joint distribution Pe(X,Y)P_e(X,Y) on the product space (×,ℱ⊗ℱ)(X×Y, F_X F_Y). The data-generating process follows one of three canonical causal structures ∗∈1,2,3G^*∈\G_1,G_2,G_3\, but the true structure ∗G^* is unknown to the learner. Definition 3.1 (Structure-Specific Invariance). For each causal structure kG_k, the corresponding invariance condition ℐkI_k on representation ϕ:→φ:X is: ℐ1(ϕ):ϕ(X)⟂E|YI_1(φ):φ(X) \!\!\! E Y for anti-causal, ℐ2(ϕ):Y⟂E|ϕ(X)I_2(φ):Y \!\!\! E φ(X) for confounded-descendant, and ℐ3(ϕ):ϕ(X)⟂EI_3(φ):φ(X) \!\!\! E for confounded-outcome case, where ⟂⟂ \!\!\! denotes statistical independence and E is the environment indicator. We formalize the structure-agnostic setting with the product causal space framework. The learner’s task is to identify which kernel family governs the observed data. Given data =Deii=1mD=\D_e_i\_i=1^m from m environments generated by an unknown causal structure ∗∈1,2,3G^*∈\G_1,G_2,G_3\, we want to identify the true causal structure ∗G^* (structure discovery), find ϕ∗:→φ^*:X satisfying the invariance ℐk∗I_k^* corresponding to ∗G^*, and learn predictor w∗:→w^*:Z achieving optimal out-of-distribution generalization simultaneously. Formally, we solve the below objective function: minϕ,w,kmaxe∈ℰℛe(w∘ϕ)+λ⋅[¬ℐk(ϕ)] _φ,w,k \ _e R_e(w φ)+λ·I[ _k(φ)] \ (1) where ℛe(f)=(X,Y)∼Pe[ℓ(f(X),Y)]R_e(f)=E_(X,Y) P_e[ (f(X),Y)] is the risk in environment e. The objective jointly selects the causal structure kG_k, learns a representation ϕφ satisfying the invariance ℐkI_k, and optimizes worst-case risk across environments. Two challenges arise: (i) the structure k is discrete, requiring search over candidates; (i) the indicator [⋅]I[·] is non-differentiable. We address both by observing that each ℐkI_k is an independence statement, which can be quantified continuously via kernel independence testing. 4 Kernel Independence Framework We develop a kernel-based framework offering rigorous independence characterization and differentiable optimization. Reproducing kernel Hilbert spaces (RKHS) (See Appendix C.5) embed probability measures as unique elements, enabling independence testing via norm computations. Definition 4.1 (Characteristic Kernel). A kernel k:×→ℝk:X×X is characteristic if the kernel mean embedding μ:()→ℋμ: P(X) is injective on the space of probability measures () P(X) defined by μP:=X∼P[k(⋅,X)]=∫k(⋅,x)P(x) _P:=E_X P[k(·,X)]= _Xk(·,x)\,dP(x) (2) Common characteristic kernels include the Gaussian RBF kernel k(x,x′)=exp(−∥x−x′∥2/2σ2)k(x,x )= (-\|x-x \|^2/2σ^2). For characteristic kernels, μP=μQ _P= _Q if and only if P=QP=Q. We quantify dependence between random variables via the cross-covariance operator XYC_XY, whose Hilbert-Schmidt norm measures dependence strength. Definition 4.2 (Hilbert-Schmidt Independence Criterion (HSIC)). Let ∥⋅∥HS\|·\|_HS be the Hilbert-Schmidt norm. The HSIC between random variables X and Y w.r.t. characteristic kernels k_X and k_Y is: HSIC(X,Y):=‖XY‖HS2HSIC(X,Y):=\|C_XY\|_HS^2 (3) Key property of HSIC is it completely characterizes statistical independence for characteristic kernels. Lemma 1 (HSIC Characterizes Independence). Let k⊗k_X k_Y be characteristic (e.g., Gaussian or delta kernels [64]). Then: HSIC(X,Y)=0⟺X⟂YHSIC(X,Y)=0 X \!\!\! Y (4) In practice, we estimate HSIC from finite samples using kernel matrices. Definition 4.3 (Empirical HSIC). Given samples (xi,yi)i=1n\(x_i,y_i)\_i=1^n, the empirical HSIC estimator is: HSIC^(X,Y)=1n2tr(KXHKYH) HSIC(X,Y)= 1n^2tr(K_XHK_YH) (5) where [KX]ij=k(xi,xj)[K_X]_ij=k_X(x_i,x_j), [KY]ij=k(yi,yj)[K_Y]_ij=k_Y(y_i,y_j), and H=In−1nnn⊤H=I_n- 1n1_n1_n is the centering matrix. To test conditional independence, we employ a residualization approach based on kernel regression. Definition 4.4 (Kernel Ridge Regression). Let (xi,zi)i=1n\(x_i,z_i)\_i=1^n be observations. The kernel ridge regression estimate is: m^X(z):=z⊤(KZ+γIn)−1 m_X(z):=k_z (K_Z+γ I_n)^-1x (6) where γ>0γ>0 is the regularization parameter, chosen via cross-validation or set to O(n−1/(2+d))O(n^-1/(2+d)) for d-dimensional inputs. The residual at observation i is x~i:=xi−m^X(zi) x_i:=x_i- m_X(z_i). These residuals remove the dependence on Z, allowing us to test conditional independence by applying standard HSIC to the residualized variables. Definition 4.5 (Conditional HSIC). The conditional HSIC for testing X⟂Y|ZX \!\!\! Y Z is: HSIC(X,Y∣Z):=HSIC(X~,Y~)HSIC(X,Y Z):=HSIC( X, Y) (7) where X~=X−m^X(Z) X=X- m_X(Z) and Y~=Y−m^Y(Z) Y=Y- m_Y(Z) are kernel regression residuals. Below lemma states kernel regression residuals yield a consistent test for conditional independence. Lemma 2 (Conditional HSIC Consistency). Assume additive noise, X=mX(Z)+εX=m_X(Z)+ _X and Y=mY(Z)+εY=m_Y(Z)+ _Y with (εX,εY)⟂Z( _X, _Y) \!\!\! Z, bounded Lipschitz kernels with k⊗k_X k_Y characteristic, and L2L^2-consistent regression estimates m^X,m^Y m_X, m_Y. Then X⟂Y|Z⇔εX⟂εYX \!\!\! Y Z _X \!\!\! _Y, and: 1. If X⟂Y|ZX \!\!\! Y Z, then HSIC^(X~,Y~)→0 HSIC( X, Y) p0. 2. If X⟂̸⟂Y|ZX \!\!\! Y Z, then HSIC^(X~,Y~)→HSIC(εX,εY)>0 HSIC( X, Y) pHSIC( _X, _Y)>0. Under this model, the violation metrics V1V_1 and V2V_2, which rely on conditional HSIC, correctly detect whether the corresponding invariances hold. 5 Structure-agnostic Causal Representation Learning (SaCRL) Our theoretical findings rely on standard regularity conditions for kernel methods and the representation class (Appendix C.6). The key structural requirement is: Assumption 1 (Sufficient Environment Diversity). The training environments exhibit variation along three axes: marginal distribution P(X)P(X), conditional P(X|Y)P(X|Y), and outcome P(Y|X)P(Y|X). We assume the three causal structures are distinguishable: testing ℐ1I_1 requires variation in P(X|Y)P(X|Y), testing ℐ2I_2 requires variation in P(Y|X)P(Y|X), and testing ℐ3I_3 requires variation in P(X)P(X). 5.1 Structure Violation Metrics We quantify the degree to which a representation ϕφ violates structure-specific invariance using HSIC. Definition 5.1 (Population Violation Functions). Let ϕ:→φ:X be a representation. The population violations are V1∗(ϕ):=HSIC(ϕ(X),E∣Y)V_1^*(φ):=HSIC(φ(X),E Y) (anti-causal), V2∗(ϕ):=HSIC(Y,E∣ϕ(X))V_2^*(φ):=HSIC(Y,E φ(X)) (confounded-descendant), and V3∗(ϕ):=HSIC(ϕ(X),E)V_3^*(φ):=HSIC(φ(X),E) (confounded-outcome). For empirical violation functions, given pooled observations (xi,yi,ei)i=1N\(x_i,y_i,e_i)\_i=1^N from all environments, the empirical violations are V1(ϕ):=HSIC^(ϕ~,E~)V_1(φ):= HSIC( φ, E), V2(ϕ):=HSIC^(Y~,E~′)V_2(φ):= HSIC( Y, E ), and V3(ϕ):=HSIC^(ϕ(X),E)V_3(φ):= HSIC(φ(X),E). Here, ϕ~i:=ϕ(xi)−m^ϕ(yi) φ_i:=φ(x_i)- m_φ(y_i) are residuals from regressing ϕ(X)φ(X) on Y, E~i:=ei−m^E(yi) E_i:=e_i- m_E(y_i) are residuals from regressing E on Y, Y~i:=yi−m^Y(ϕ(xi)) Y_i:=y_i- m_Y(φ(x_i)) are residuals from regressing Y on ϕ(X)φ(X) E~i′:=ei−m^E′(ϕ(xi)) E _i:=e_i- m_E (φ(x_i)) are residuals from regressing E on ϕ(X)φ(X), and m^(⋅) m(·) kernel ridge regression estimates. The distinction between population (Vk∗V_k^*) and empirical (VkV_k) violations is essential. Theoretical results concern Vk∗V_k^*, while the algorithm optimizes VkV_k. Consistency ensures Vk(ϕ)→Vk∗(ϕ)V_k(φ) pV_k^*(φ) as n→∞n→∞. 5.2 Soft Structure Selection SoftMin offers a differentiable relaxation for joint structure discovery and representation learning. Definition 5.2 (SoftMin Aggregation). SoftMin operator on violation metrics with temperature β is: SoftMinβ(Vk):=−1βlog∑kexp(−βVk)SoftMin_β(\V_k\):=- 1β _k (-β V_k) (8) Definition 5.3 (Adaptive Structure Weights). The adaptive weight for structure k given representation ϕφ is: αk(ϕ):=exp(−βVk(ϕ))∑jexp(−βVj(ϕ)) _k(φ):= (-β V_k(φ)) _j (-β V_j(φ)) (9) These weights satisfy αk≥0 _k≥ 0 and ∑kαk=1 _k _k=1. Lemma 3 (SoftMin Properties). SoftMin operator satisfies: 1. Bounds: minkVk−log3β≤SoftMinβ(Vk)≤minkVk _kV_k- 3β _β(\V_k\)≤ _kV_k 2. Limiting behavior: limβ→∞SoftMinβ(Vk)=minkVk _β→∞SoftMin_β(\V_k\)= _kV_k 3. Gradient: ∇ϕSoftMinβ=∑kαk(ϕ)∇ϕVk _φSoftMin_β= _k _k(φ) _φV_k 5.3 Unified Objective Function Definition 5.4 (SaCRL Objective). The structure-agnostic representation learning objective is: ℒ(ϕ,w):=maxe∈ℰℛ^e(w∘ϕ)⏟worst-case prediction loss+λSoftMinβ(Vk(ϕ))⏟adaptive invariance penaltyL(φ,w):= _e R_e(w φ)_worst-case prediction loss+λ SoftMin_β (\V_k(φ)\ )_adaptive invariance penalty (10) where ℛ^e(f)=1ne∑j=1neℓ(f(xje),yje) R_e(f)= 1n_e _j=1^n_e (f(x_j^e),y_j^e) is the empirical risk in the environment e. The objective can be written with explicit structure weights: ℒ(ϕ,w)=maxe∈ℰℛ^e(w∘ϕ)+λ∑kαk(ϕ)Vk(ϕ)L(φ,w)= _e R_e(w φ)+λ _k _k(φ)V_k(φ) (11) where αk(ϕ) _k(φ) are the adaptive weights (Definition 5.3). Definition 5.5 (Distributionally Robust Risk). For risk level ρ∈(0,1]ρ∈(0,1], the conditional value-at-risk (CVaR) over environments is defined as: CVaRρ(ℛ^e):=inft∈ℝ(t+1ρ|ℰ|∑e∈ℰmax(ℛ^e(w∘ϕ)−t,0))CVaR_ρ(\ R_e\):= _t (t+ 1ρ|E| _e ( R_e(w φ)-t,0) ) (12) CVaR interpolates between the average risk (ρ=1ρ=1) and the worst-case risk maxeℛ^e _e R_e (ρ=1/|ℰ|ρ=1/|E|); ρ=0.5ρ=0.5 (our default) averages the worst half of the environments. Definition 5.6 (Distributionally Robust Objective). The distributionally robust structure-agnostic objective is: ℒρ(ϕ,w):=CVaRρ(ℛ^e(w∘ϕ)e∈ℰ)+λSoftMinβ(Vk(ϕ)k)L_ρ(φ,w):=CVaR_ρ (\ R_e(w φ)\_e )+λ\,SoftMin_β (\V_k(φ)\_k ) (13) ρ=1/|ℰ|ρ=1/|E| recovers Definition 5.4. We use ρ=1/|ℰ|ρ=1/|E| for theoretical analysis and ρ=0.5ρ=0.5 in practice. CVaRρCVaR_ρ is convex with a closed-form subgradient (Proposition 1, Appendix F). Optimization: The problem admits a bilevel decomposition, with the structure weights α in the outer problem and (ϕ,w)(φ,w) in the inner problem (Appendix F.1).The optimization problem is challenging due to (i) the non-convex loss landscape with multiple local minima, and (i) the difficulty of jointly optimizing all variables. To address these challenges, we adopt the following observations and optimization strategies: 1. Single-loop relaxation: Instead of solving the bilevel optimization problem explicitly, we employ a single-loop algorithm in which α is implicitly determined via a SoftMin operator. 2. Multi-start initialization: We run the optimization from M≥10M≥ 10 random initializations and select the solution with the lowest objective value ℒL. 3. Temperature annealing: We initialize the SoftMin temperature β with a small value to allow soft selection, and gradually increase it to βmax _ to encourage hard selection. 4. Validation-based selection: Among the resulting stationary points, we select the model that achieves the best performance on held-out environments. Because V1V_1 and V2V_2 condition on the current ϕφ, the conditioning target moves during training. Our guarantees (Section 6) are stated at the converged representation, where conditioning is on a fixed ϕφ. In practice, the ERM warm start, temperature annealing, and cross-fitted residuals (Appendix H.7) stabilize training, and conditioning on an exponential-moving-average (EMA) copy of the encoder gives nearly identical results (Appendix H.8). The soft weights αk _k only guide optimization: the identified structure is the hard decision k^=argminkVk(ϕT) k= _kV_k( _T) made after training (Algorithm 1). Computation Speedup: Computing empirical HSIC (Definition 4.3) requires O(n2)O(n^2) time and storage, which is prohibitive for large-scale datasets. To improve scalability, we employ efficient approximation techniques using the random Fourier features, with Algorithm details in Appendix F. With D random Fourier features ψ [50] (Definition F.1, Appendix F), the approximate HSIC is HSIC^RFF(X,Y)=‖Σ^ψXψY−μ^ψXμ^ψY⊤‖F2 HSIC_RFF(X,Y)=\| _ _X _Y- μ_ _X μ_ _Y \|_F^2, which costs O(nD2)O(nD^2) time and O(nD)O(nD) memory, versus O(n2)O(n^2) for both with exact HSIC; memory saving makes kernel penalties feasible on image data (Appendix H.11). The identification guarantee survives this approximation: with D=Ω(κ4η−2dlog(n/δ′))D= (κ^4η^-2d (n/δ )) random features, where η is the faithfulness gap, k^=k∗ k=k^* with probability at least 1−6exp(−nη2/(8192κ4))−δ′1-6 (-nη^2/(8192κ^4))-δ (Corollary 1). Algorithm details are deferred to Appendix F. Selection Confidence: SaCRL reports two confidence measures for k k: the weight margin α(1)−α(2) _(1)- _(2) between the two largest adaptive weights, and the permutation p-value of minkVk _kV_k against the null threshold τα _α (Definition G.2). A large margin means one structure clearly fits best. If minkVk>τα _kV_k> _α, no candidate invariance holds exactly, and k k should be read as the closest canonical structure. We calibrate the margin against identification accuracy on data with known structure (Section 7). Heterogeneous Structures: When different feature subsets follow different causal patterns, we extend the framework to per-feature structure assignment σ∗:1,…,d→1,2,3,⊥σ^*:\1,…,d\→\1,2,3, \ via the objective ℒhet(ϕ,w)=maxeℛ^e(w∘ϕ)+λ∑j=1dSoftMinβ(Vk(j)(ϕ)k)L_het(φ,w)= _e R_e(w φ)+λ _j=1^dSoftMin_β(\V_k^(j)(φ)\_k), applying SoftMin independently per feature. A calibrated permutation threshold τα(j) _α^(j) controls the false-positive rate for features with σ∗(j)=⊥σ^*(j)= . Per-feature identification guarantees mirror Theorem D.2 with sample cost scaling as O(d)O(d). Full treatment is in Appendix G. 6 Theoretical Guarantees We establish four theoretical results: (1) violation characterization showing the correct structure is identifiable, (2) finite-sample structure identification holding under random-feature approximation (Corollary 1), (3) representation optimality, and (4) out-of-distribution generalization bounds. The details of the theoretical results and proofs are discussed in Appendix D and Appendix E, respectively. 7 Experiments We evaluate SaCRL across synthetic and real-world settings, addressing four questions: Q1: Can the method identify the true causal structure from data? Q2: Does correct structure identification improve OOD generalization? Q3: Does SaCRL outperform fixed-invariance baselines on standard domain generalization benchmarks? Q4: How does it behave under structural misspecification and limited environment diversity? Additional experiments are deferred to Appendix H. 7.1 Experimental Setup Baselines. We compare against: ERM (no invariance), IRM [4] (enforces ℐ2I_2, i.e., assumes 2G_2), VREx [30] (variance-based regularization), and CIRCE [47] (enforces ℐ1I_1, i.e., assumes 1G_1). Metrics. For synthetic data: we use structure identification accuracy (percentage of runs identifying the correct kG_k) and OOD test accuracy on a held-out environment with larger distribution shift. For real-world datasets, we use OOD accuracy on the anti-correlated test environment. Results are means over 5 random seeds, and ± denotes the standard deviation across seeds. Implementation Details. For synthetic data, we use 3-layer MLPs with hidden dimensions [256, 128] for both encoder and predictor. For Colored MNIST, we use a 4-layer CNN encoder with base channels 32. For PACS, VLCS, and OfficeHome, we use a frozen ImageNet-pretrained ResNet18 followed by an MLP encoder with hidden dimensions [256, 128] and a linear predictor; the same model produces both the identified structures (Table 4) and the accuracies (Table 5). We select the invariance penalty weight λ∈0.5,1.0,2.0,5.0λ∈\0.5,1.0,2.0,5.0\ per experiment by OOD accuracy, i.e., test-domain (oracle) selection in the terminology of Gulrajani and Lopez-Paz [20] (λ=5.0λ=5.0 for Section 7.6; sensitivity in Appendix H.3), and use CVaR with ρ=0.5ρ=0.5. We train with Adam optimizer (learning rate 10−410^-4 for synthetic, 5×10−55× 10^-5 for images) using balanced environment sampling. Random Fourier Features (RFF) with D=1000D=1000 dimensions are used to approximate HSIC, reducing memory from O(n2)O(n^2) to O(nD)O(nD) with time linear in n. Full hyperparameters are provided in Appendix H.1. 7.2 Experiment 1: Structure Identification Setup. We generate synthetic data from each of the three canonical causal structures 1,2,3G_1,G_2,G_3 (Definition 3.1). For each structure, we create m=3m=3 environments with n total samples (divided equally across environments) and d=10d=10 features. We vary the sample size n∈500,1000,2000,5000n∈\500,1000,2000,5000\ to study the effect of data availability on structure identification accuracy. For each kG_k, environments differ in noise scale (1G_1, 3G_3) or confounder mean (2G_2); weight matrices are sampled from (0,0.52)N(0,0.5^2) and fixed across runs. Full generative equations are in Appendix H.2. Results. Table 1 reports structure identification accuracy. Identification accuracy reaches at least 95% at n=5000n=5000 and at least 89% at n=2000n=2000, demonstrating that the violation-based approach successfully distinguishes the three canonical structures. Table 1: Structure identification accuracy (%) across sample sizes and true structures. Results averaged over 20 random seeds. Standard deviations in parentheses. True Structure n=500n=500 n=1000n=1000 n=2000n=2000 n=5000n=5000 1G_1 (anti-causal) 73.5 (2.8) 85.0 (2.1) 92.5 (1.4) 96.5 (0.9) 2G_2 (conf-desc) 68.0 (3.2) 81.5 (2.5) 89.0 (1.8) 95.0 (1.1) 3G_3 (conf-out) 78.5 (2.4) 88.5 (1.9) 94.0 (1.2) 97.5 (0.7) Performance improves monotonically with sample size, consistent with the consistency result in Theorem D.2. Structure 3G_3 is identified most reliably because its invariance ϕ(X)⟂Eφ(X) E is unconditional and thus has higher statistical power. Structure 2G_2 is most challenging in small-sample regimes due to the need to estimate conditional independence Y⟂E|ϕ(X)Y E φ(X). Semi-synthetic Bayesian networks. To test identification where the ground truth is fixed independently of our method, we build five benchmarks from three published Bayesian networks (ASIA, SACHS, ALARM). For each, we choose a target Y, take X as its Markov blanket, induce environments by soft interventions on a latent parent outside the blanket, and read the canonical structure off the published DAG by the rules of Section 2; each structure is verified in simulation before training (Appendix H.6). These topologies admit 2G_2 and 3G_3 but not 1G_1, which would require the Markov blanket of Y to consist only of its descendants. SaCRL selects the correct structure in the majority of runs (75–95%, Table 2), whereas returning any single structure by default is correct on at most 3/5. Accuracy and margin are highest on the most separated benchmark (asia_g3) and lowest on the least separated one (alarm_g2), as predicted by the dependence of Theorem D.2 on the faithfulness gap η. Table 2: Semi-synthetic Bayesian-network benchmarks. Violation: the competing invariance’s violation before training, measured in standard deviations of the permutation null. Margin: α(1)−α(2) _(1)- _(2). Benchmark Network (nodes) Y True Violation ID acc. (%) Margin asia_g3 ASIA (8) either 3G_3 56.3 95.0 (19/20) 0.45 asia_g2 ASIA (8) either 2G_2 24.6 90.0 (18/20) 0.42 alarm_co ALARM (37) CO 3G_3 28.3 85.0 (17/20) 0.35 sachs_g2 SACHS (11) Akt 2G_2 6.4 85.0 (17/20) 0.25 alarm_g2 ALARM (37) EXPCO2 2G_2 4.2 75.0 (15/20) 0.18 7.3 Experiment 2: OOD Generalization Setup. Following [4], we construct a binary classification task on MNIST digits (0–4 vs. 5–9) where color is spuriously correlated with the label. Specifically, for each digit, we assign color (red or green) based on its label with probability ρe _e in environment e. Training uses two environments with ρ1=0.9 _1=0.9 and ρ2=0.8 _2=0.8, while the test environment has ρtest=0.1 _test=0.1 (anti-correlation). The true causal structure is anti-causal (1G_1): the digit label causes the color assignment, not vice versa. Results. We subsample 50005000 examples per training environment and 50005000 test examples to ensure balanced evaluation. Table 3 reports the results. Table 3: Colored MNIST OOD Acc. on anti-correlated test environment. Method OOD Acc (%) ID. Stru. ERM 27.2 ± 11.3 — IRM 21.3 ± 13.2 2G_2 (fixed) VREx 24.3 ± 12.8 — CIRCE 10.3 ± 0.5 1G_1 (fixed) SaCRL 45.2 ± 10.0 1G_1 (92%) SaCRL identifies 1G_1 in 92% of runs and attains the highest OOD accuracy (45.2%, +18p over ERM). ERM relies on the spurious color. IRM enforces ℐ2I_2, the invariance of the wrong structure, and falls below ERM, which is the misspecification cost illustrated in Figure 3. CIRCE targets the correct invariance ℐ1I_1, yet its fixed penalty does not overcome the color shortcut under our shared training budget. SaCRL selects the invariance from data, without prior knowledge of the structure. 7.4 Experiment 3: Structure Analysis on PACS The PACS dataset [31] is an important benchmark for image classification under domain shift. It contains images from four distinct visual domains (Photo, Art Painting, Cartoon, and Sketch), each depicting the seven object categories (dog, elephant, giraffe, guitar, horse, house, and person) under different styles and appearances. Overall, the dataset includes 9,991 images, with domain-specific subsets ranging from natural photographs to artistic and hand-drawn sketches, making it a challenging testbed for methods to generalize across unseen domains. Table 4: Structure selected on PACS; each row trains on the other three domains (5 seeds). Margin: α(1)−α(2) _(1)- _(2); p: permutation p-value of minkVk _kV_k; Seeds: agreement of k k across seeds; Cal. acc.: identification accuracy of synthetic runs in the same margin bin (Appendix H.5). Held-out α1 _1 α2 _2 α3 _3 k k Margin p Seeds Cal. acc. (%) Photo 0.21 0.68 0.11 2G_2 0.47 0.008 5/5 86.1 Art Painting 0.59 0.31 0.10 1G_1 0.28 0.021 5/5 65.0 Cartoon 0.35 0.52 0.13 2G_2 0.17 0.037 4/5 65.0 Sketch 0.63 0.25 0.12 1G_1 0.38 0.014 5/5 86.1 Table 4 reports the adaptive weights and the confidence measures of Section 5 for each leave-one-domain-out split. Since PACS has no ground-truth DAG, we make three verifiable claims. (i) The selection responds to the data: two different structures are selected across splits, so SaCRL does not return a default. (i) Confidence tracks the margin: seeds agree 5/5 on the three splits with margin at least 0.280.28 and 4/5 on Cartoon, the split with the smallest margin; on synthetic data with known structure, the margins of Photo and Sketch fall in a bin with 86.1% identification accuracy, and those of Art and Cartoon in a bin with 65.0%. (i) No canonical invariance holds exactly: p<0.05p<0.05 on every split, so k k is the closest canonical structure for each training set rather than an identified one. 7.5 Experiment 4: Domain Generalization on DomainBed Setup. We further evaluate SaCRL on DomainBed [20], a benchmark suite to test how well models generalize to unseen domains. Besides PACS [31] (4 domains, 7 classes), we choose two additional benchmarks: VLCS [67] (4 domains: Caltech101, LabelMe, SUN09, VOC2007; 5 classes), and OfficeHome [70] (4 domains: Art, Clipart, Product, Real-World; 65 classes). Each spans diverse domain shifts in style, viewpoint, and resolution. We follow the leave-one-domain-out protocol with the model of Section 7.1. Freezing the backbone holds representation capacity fixed, so differences across methods trace to the objective rather than to encoder optimization. Table 5: Results on DomainBed benchmarks: average leave-one-domain-out OOD accuracy (%). PACS VLCS OfficeHome Method Photo Art Cartoon Sketch Avg Caltech LabelMe SUN09 VOC Avg Art Clipart Product Real Avg ERM 93.2 66.0 51.8 37.6 62.2 94.0 60.2 69.9 73.0 74.3 49.9 41.7 66.7 67.7 56.5 IRM 93.2 66.2 51.7 37.6 62.2 95.0 60.1 70.1 72.6 74.5 49.0 41.2 66.5 67.7 56.1 VREx 93.0 65.5 51.6 37.9 62.0 93.3 60.2 69.8 71.0 73.6 48.6 41.2 66.2 67.6 55.9 CIRCE 92.5 64.5 51.4 37.3 61.4 93.8 60.0 67.7 69.7 72.8 49.0 41.9 66.0 67.4 56.1 SaCRL 94.0 68.5 54.0 40.5 64.3 97.7 64.4 75.7 77.1 78.7 52.9 46.5 69.1 69.8 59.6 Table 5 shows the comparison results with recent baselines. SaCRL achieves the best per-domain and average accuracy on all benchmarks, with average improvements over +2.1p (PACS), +4.2p (VLCS), and +3.1p (OfficeHome). Fixed-invariance baselines (IRM, VREx, CIRCE) show no consistent improvement, confirming that adaptive structure is essential when the governing causal structure varies across domains. Fine-tuned backbone. Fine-tuning ResNet18 end-to-end (PACS, Photo held out, 3 seeds) raises all methods by 1.0–2.0p; SaCRL (95.0) is on par with ERM (95.1) and IRM (95.2) (Appendix H.12). Photo is ImageNet-like, so the pretrained backbone is near saturation and the choice of objective matters little; SaCRL’s gains concentrate on the stylized domains and on VLCS (Table 5). 7.6 Experiment 5: Robustness to Misspecification and Limited Diversity Table 6: Robustness under misspecification (ν) and limited diversity (δmin _ ). OOD accuracy (%) over 5 seeds. ID. Stru.: majority structure (seeds agreeing). Setting ERM IRM VREx CIRCE SaCRL ID. Stru. Misspecification (δmin=1.0 _ =1.0): ν=0.25ν=0.25 (near-causal) 34.4 ± 2.7 — — — 50.2 ± 3.7 2G_2 (4/5) ν=0.50ν=0.50 (bidirectional) 61.2 ± 3.8 61.3 ± 4.2 61.4 ± 4.4 66.5 ± 2.6 73.9 ± 4.1 2G_2 (5/5) ν=0.75ν=0.75 (near-anticausal) 98.5 ± 0.6 — — — 99.6 ± 0.3 1G_1 (4/5) Limited diversity (ν=0.50ν=0.50): δmin=1.0 _ =1.0 (full) 61.2 ± 3.8 61.3 ± 4.2 61.4 ± 4.4 66.5 ± 2.6 73.9 ± 4.1 2G_2 (5/5) δmin=0.5 _ =0.5 (moderate) 56.9 ± 2.5 — — — 67.2 ± 1.8 2G_2 (5/5) δmin=0.2 _ =0.2 (low) 53.1 ± 1.7 — — — 59.7 ± 3.7 2G_2 (4/5) δmin=0.05 _ =0.05 (minimal) 50.2 ± 2.0 — — — 57.2 ± 3.7 2G_2 (5/5) We stress-test SaCRL when (a) the true DGP falls outside the three canonical structures and (b) environment diversity is reduced. We use a bidirectional feedback DGP X↔YX Y with mixing coefficient ν∈[0,1]ν∈[0,1]: ν=0ν=0 is causal X→YX→Y (2G_2), ν=1ν=1 is anti-causal Y→XY→X (1G_1), and intermediate values fit none of 1G_1–3G_3. Diversity is controlled by δmin _ (Definition C.19). SaCRL outperforms every evaluated baseline in each setting (Table 6). Near the endpoints it selects the nearer structure (2G_2 at ν=0.25ν=0.25, 1G_1 at ν=0.75ν=0.75). At ν=0.50ν=0.50, where no candidate holds, the weights stay spread (α1≈α2≈0.42 _1≈ _2≈ 0.42, α3≈0.16 _3≈ 0.16): SaCRL acts as a soft model average instead of committing to a wrong structure, and still exceeds ERM by 12.7p. A finer sweep (Appendix H.4) shows α1 _1 increasing and α2 _2 decreasing monotonically in ν, crossing at ν=0.50ν=0.50, with SoftMin outperforming hard argmin selection throughout and by the widest margin at the crossover. Under reduced diversity, identification remains stable and the gain over ERM persists at δmin=0.05 _ =0.05, consistent with Assumption 1 governing the informativeness of the penalty rather than its feasibility. 7.7 Ablation Studies Table 7: Ablation on 2G_2 (n=2000n=2000). Default: SoftMin, β annealing 1→501→50, λ=1.0λ=1.0, ρ=0.5ρ=0.5. Variant Stru. ID (%) OOD Acc (%) Full method 89.0 84.2 Hard argmink _k 79.5 79.8 Fixed β=10β=10 83.0 81.5 Fixed λ=0.1λ=0.1 72.5 78.3 Fixed λ=5.0λ=5.0 88.5 81.0 Uniform (ρ=1ρ=1) 88.0 82.1 Worst-case (ρ=1/|ℰ|ρ=1/|E|) 87.5 81.4 We conduct ablations on synthetic 2G_2 data with n=2000n=2000 samples to isolate the contribution of each component. Results are shown in Table 7. Replacing SoftMin with hard argmink _k selection during training lowers identification accuracy from 89.0% to 79.5%: early in training, violation estimates are noisy, and hard selection commits prematurely to whichever structure has the smallest estimate. SoftMin delays commitment until the representation is informative; the reported structure is still the hard decision k k made after training. The effect of varying the number of training environments m is reported in Appendix H.10, where we observe that m≥3m≥ 3 yields reliable structure identification while m=2m=2 is unstable. Hyperparameter sensitivity sweeps and a cross-fitting variant of conditional HSIC are reported in Appendix H.3 and Appendix H.7, respectively. Computational efficiency analysis (training time, memory, runtime breakdown) is in Appendix H.11. 8 Conclusion We proposed a structure-agnostic approach to causal representation learning that formulates causal structure selection as a soft optimization problem over HSIC-based invariance violation metrics. The key insight is that correct invariances can be achieved with zero violation while preserving predictive accuracy, whereas incorrect invariances cannot, enabling automatic structure identification through optimization. We developed (i) a unified measure-theoretic framework that extends prior work on anti-causal learning to confounded-descendant and confounded-outcome settings; (i) an efficient learning algorithm whose complexity scales linearly with the number of candidate structures, leveraging SoftMin aggregation and RFF-based HSIC approximations; and (i) theoretical guarantees for structure identification, under random-feature approximation, invariance satisfaction, and OOD generalization. Empirically, SaCRL recovers the true structure on synthetic and semi-synthetic Bayesian-network benchmarks with a calibrated confidence, surpasses fixed-invariance baselines on Colored MNIST and three DomainBed benchmarks, and remains robust under bidirectional dynamics and limited environment diversity. Acknowledgments and Disclosure of Funding We thank the anonymous reviewers for their constructive comments. This work was supported in part by the National Science Foundation under Grant Nos. ECCS-2216926, CCF-2331302, CNS-2241713, and CNS-2339686. References [1] K. Ahuja, E. Caballero, D. Zhang, J. Gagnon-Audet, Y. Bengio, I. Mitliagkas, and I. Rish (2021) Invariance principle meets information bottleneck for out-of-distribution generalization. Advances in Neural Information Processing Systems. Cited by: Appendix A. [2] K. Ahuja, D. Mahajan, Y. Wang, and Y. Bengio (2023) Interventional causal representation learning. In International conference on machine learning, Cited by: Appendix A, Appendix J. [3] K. Ahuja, A. Mansouri, and Y. Wang (2024) Multi-domain causal representation learning via weak distributional invariances. In Artificial Intelligence and Statistics, Cited by: Appendix A. [4] M. Arjovsky, L. Bottou, I. Gulrajani, and D. Lopez-Paz (2019) Invariant risk minimization. arXiv preprint arXiv:1907.02893. Cited by: Appendix A, Appendix J, §1, §7.1, §7.3. [5] S. Barocas, M. Hardt, and A. Narayanan (2023) Fairness and machine learning: limitations and opportunities. MIT press. Cited by: Appendix J. [6] A. Behnam and B. Wang (2024) Graph neural network causal explanation via neural causal models. In European conference on computer vision, p. 410–427. Cited by: §1. [7] A. Behnam and B. Wang (2025) Measure-theoretic anti-causal representation learning. Advances in Neural Information Processing Systems 38, p. 61375–61431. Cited by: Appendix A, §C.4, §1, §2.1. [8] J. Brehmer, P. De Haan, P. Lippe, and T. S. Cohen (2022) Weakly supervised causal representation learning. Advances in Neural Information Processing Systems. Cited by: Appendix A. [9] S. Buchholz, G. Rajendran, E. Rosenfeld, B. Aragam, B. Schölkopf, and P. Ravikumar (2024) Learning linear causal representations from interventions under general nonlinear mixing. Advances in Neural Information Processing Systems 36. Cited by: Appendix A. [10] T. Cai, H. Namkoong, and S. Yadlowsky (2026) Diagnosing model performance under distribution shift. Operations Research 74 (2), p. 898–916. Cited by: Appendix J. [11] D. C. Castro, I. Walker, and B. Glocker (2020) Causality matters in medical imaging. Nature Communications 11 (1), p. 3673. Cited by: Appendix J. [12] J. Cha, S. Chun, K. Lee, H. Cho, S. Park, Y. Lee, and S. Park (2021) Swad: domain generalization by seeking flat minima. Advances in Neural Information Processing Systems 34, p. 22405–22418. Cited by: Appendix A. [13] Y. Chen, W. Huang, K. Zhou, Y. Bian, B. Han, and J. Cheng (2024) Understanding and improving feature learning for out-of-distribution generalization. Advances in Neural Information Processing Systems 36. Cited by: Appendix A. [14] Y. Chen, Y. Zhang, Y. Bian, H. Yang, M. Kaili, B. Xie, T. Liu, B. Han, and J. Cheng (2022) Learning causally invariant representations for out-of-distribution generalization on graphs. In Advances in Neural Information Processing Systems, Cited by: §1. [15] A. Chouldechova (2017) Fair prediction with disparate impact: a study of bias in recidivism prediction instruments. Big data 5 (2), p. 153–163. Cited by: Appendix J. [16] C. M. Fang, A. R. Liu, V. Danry, E. Lee, S. W. Chan, P. Pataranutaporn, P. Maes, J. Phang, M. Lampe, L. Ahmad, et al. (2025) How ai and human behaviors shape psychosocial effects of extended chatbot use: a longitudinal randomized controlled study. arXiv preprint arXiv:2503.17473. Cited by: Appendix J. [17] S. G. Finlayson, A. Subbaswamy, K. Singh, J. Bowers, A. Kupke, J. Zittrain, I. S. Kohane, and S. Saria (2021) The clinician and dataset shift in artificial intelligence. New England Journal of Medicine 385 (3), p. 283–286. Cited by: Appendix J. [18] Y. Ganin, E. Ustinova, H. Ajakan, P. Germain, H. Larochelle, F. Laviolette, M. March, and V. Lempitsky (2016) Domain-adversarial training of neural networks. Journal of machine learning research 17 (59), p. 1–35. Cited by: Appendix J. [19] S. Gui, M. Liu, X. Li, Y. Luo, and S. Ji (2024) Joint learning of label and environment causal independence for graph out-of-distribution generalization. Advances in Neural Information Processing Systems 36. Cited by: Appendix A. [20] I. Gulrajani and D. Lopez-Paz (2021) In search of lost domain generalization. In International Conference on Learning Representations, Cited by: Appendix A, Appendix J, §7.1, §7.5. [21] C. Heinze-Deml, J. Peters, and N. Meinshausen (2018) Invariant causal prediction for nonlinear models. Journal of Causal Inference 6 (2), p. 20170016. Cited by: Appendix A, Appendix J. [22] A. Hyvärinen and P. Pajunen (1999) Nonlinear independent component analysis: existence and uniqueness results. Neural networks 12 (3), p. 429–439. Cited by: Appendix J. [23] R. Jiang, S. Chiappa, T. Lattimore, A. György, and P. Kohli (2019) Degenerate feedback loops in recommender systems. In Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society, p. 383–390. Cited by: Appendix J. [24] Y. Jiang and V. Veitch (2022) Invariant and transportable representations for anti-causal domain shifts. Advances in Neural Information Processing Systems 35, p. 20782–20794. Cited by: Appendix A. [25] C. Jones, D. C. Castro, F. De Sousa Ribeiro, O. Oktay, M. McCradden, and B. Glocker (2024) A causal perspective on dataset bias in machine learning for medical imaging. Nature Machine Intelligence 6 (2), p. 138–146. Cited by: Appendix J. [26] J. N. Kaur, E. Kiciman, and A. Sharma (2022) Modeling the data-generating process is necessary for out-of-distribution generalization. arXiv preprint arXiv:2206.07837. Cited by: Appendix J. [27] I. Khemakhem, D. Kingma, R. Monti, and A. Hyvarinen (2020) Variational autoencoders and nonlinear ica: a unifying framework. In International conference on artificial intelligence and statistics, p. 2207–2217. Cited by: Appendix J. [28] H. R. Kirk, A. Whitefield, P. Röttger, A. Bean, K. Margatina, J. Ciro, R. Mosquera, M. Bartolo, A. Williams, H. He, et al. (2024) The prism alignment dataset: what participatory, representative and individualised human feedback reveals about the subjective and multicultural alignment of large language models. Advances in Neural Information Processing Systems 37, p. 105236–105344. Cited by: Appendix J. [29] J. Kleinberg, S. Mullainathan, and M. Raghavan (2016) Inherent trade-offs in the fair determination of risk scores. arXiv preprint arXiv:1609.05807. Cited by: Appendix J. [30] D. Krueger, E. Caballero, J. Jacobsen, A. Zhang, J. Binas, D. Zhang, R. Le Priol, and A. Courville (2021) Out-of-distribution generalization via risk extrapolation (rex). In International Conference on Machine Learning, Cited by: Appendix A, §1, §7.1. [31] D. Li, Y. Yang, Y. Song, and T. M. Hospedales (2017) Deeper, broader and artier domain generalization. In IEEE international conference on computer vision, Cited by: §7.4, §7.5. [32] H. Li, Z. Zhang, X. Wang, and W. Zhu (2022) Learning invariant graph representations for out-of-distribution generalization. Advances in Neural Information Processing Systems. Cited by: §1. [33] Z. Lipton, Y. Wang, and A. Smola (2018) Detecting and correcting for label shift with black box predictors. In International conference on machine learning, p. 3122–3130. Cited by: Appendix J. [34] C. Liu, X. Sun, J. Wang, H. Tang, T. Li, T. Qin, W. Chen, and T. Liu (2021) Learning causal semantic representation for out-of-distribution prediction. Advances in Neural Information Processing Systems 34, p. 6155–6170. Cited by: Appendix A. [35] J. Liu, T. Wang, P. Cui, and H. Namkoong (2023) On the need for a language describing distribution shifts: illustrations on tabular datasets. Advances in Neural Information Processing Systems 36, p. 51371–51408. Cited by: Appendix J. [36] F. Locatello, S. Bauer, M. Lucic, G. Raetsch, S. Gelly, B. Schölkopf, and O. Bachem (2019) Challenging common assumptions in the unsupervised learning of disentangled representations. In international conference on machine learning, p. 4114–4124. Cited by: Appendix J. [37] M. Long, Z. Cao, J. Wang, and M. I. Jordan (2018) Conditional adversarial domain adaptation. Advances in neural information processing systems 31. Cited by: Appendix A. [38] C. Lu, Y. Wu, J. M. Hernández-Lobato, and B. Schölkopf (2021) Invariant causal representation learning for out-of-distribution generalization. In International Conference on Learning Representations, Cited by: Appendix A. [39] F. Lv, J. Liang, S. Li, B. Zang, C. H. Liu, Z. Wang, and D. Liu (2022) Causality inspired representation learning for domain generalization. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, Cited by: Appendix A. [40] M. Makar and A. D’Amour (2022) Fairness and robustness in anti-causal prediction. arXiv preprint arXiv:2209.09423. Cited by: Appendix A, Appendix J, §1. [41] J. Mitrovic, D. Sejdinovic, and Y. W. Teh (2018) Causal inference via kernel deviance measures. Advances in neural information processing systems 31. Cited by: Appendix A. [42] H. Nam, H. Lee, J. Park, W. Yoon, and D. Yoo (2021) Reducing domain gap by reducing style bias. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 8690–8699. Cited by: Appendix A. [43] J. Perdomo, T. Zrnic, C. Mendler-Dünner, and M. Hardt (2020) Performative prediction. In International conference on machine learning, p. 7599–7609. Cited by: Appendix J. [44] J. Peters, P. Bühlmann, and N. Meinshausen (2016) Causal inference by using invariant prediction: identification and confidence intervals. Journal of the Royal Statistical Society Series B: Statistical Methodology 78 (5), p. 947–1012. Cited by: Appendix A, Appendix J. [45] J. Peters, D. Janzing, and B. Schölkopf (2017) Elements of causal inference: foundations and learning algorithms. MIT Press. Cited by: §1. [46] M. Pezeshki, D. Bouchacourt, M. Ibrahim, N. Ballas, P. Vincent, and D. Lopez-Paz (2023) Discovering environments with xrm. arXiv preprint arXiv:2309.16748. Cited by: Appendix A. [47] R. Pogodin, N. Deka, Y. Li, D. J. Sutherland, V. Veitch, and A. Gretton (2022) Efficient conditionally invariant representation learning. arXiv. Cited by: §7.1. [48] R. Pogodin, A. Schrab, Y. Li, D. J. Sutherland, and A. Gretton (2024) Practical kernel tests of conditional independence. arXiv preprint arXiv:2402.13196. Cited by: Appendix J, §H.7. [49] S. Rabanser, S. Günnemann, and Z. Lipton (2019) Failing loudly: an empirical study of methods for detecting dataset shift. Advances in Neural Information Processing Systems 32. Cited by: Appendix J. [50] A. Rahimi and B. Recht (2007) Random features for large-scale kernel machines. Advances in neural information processing systems 20. Cited by: §E.6, §5.3, Corollary 1. [51] A. Ramé, K. Ahuja, J. Zhang, M. Cord, L. Bottou, and D. Lopez-Paz (2023) Model ratatouille: recycling diverse models for out-of-distribution generalization. In International Conference on Machine Learning, p. 28656–28679. Cited by: Appendix A. [52] K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan (2005) Causal protein-signaling networks derived from multiparameter single-cell data. Science 308 (5721), p. 523–529. Cited by: Appendix J. [53] N. Scherrer, O. Bilaniuk, Y. Annadani, A. Goyal, P. Schwab, B. Schölkopf, M. C. Mozer, Y. Bengio, S. Bauer, and N. R. Ke (2021) Learning neural causal models with active interventions. arXiv preprint arXiv:2109.02429. Cited by: Appendix A. [54] B. Schölkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and J. Mooij (2012) On causal and anticausal learning. arXiv preprint arXiv:1206.6471. Cited by: Appendix J. [55] B. Schölkopf, F. Locatello, S. Bauer, N. R. Ke, N. Kalchbrenner, A. Goyal, and Y. Bengio (2021) Toward causal representation learning. Proceedings of the IEEE 109 (5), p. 612–634. Cited by: Appendix J, §1. [56] J. Schrouff, A. Bellot, A. Rannen-Triki, A. Malek, I. Albuquerque, A. Gretton, A. D’Amour, and S. Chiappa (2024) Mind the graph when balancing data for fairness or robustness. Advances in Neural Information Processing Systems 37, p. 29913–29947. Cited by: Appendix J. [57] R. D. Shah and J. Peters (2020) The hardness of conditional independence testing and the generalised covariance measure. Cited by: Appendix J. [58] M. Sharma, M. Tong, T. Korbak, D. Duvenaud, A. Askell, S. Bowman, E. Durmus, Z. Hatfield-Dodds, S. Johnston, S. Kravec, et al. (2024) Towards understanding sycophancy in language models. In International Conference on Learning Representations, Vol. 2024, p. 110–144. Cited by: Appendix J. [59] J. Shi, G. Gare, J. Tian, S. Chai, Z. Lin, A. Vasudevan, D. Feng, F. Ferroni, and S. Kong (2024) LCA-on-the-line: benchmarking out-of-distribution generalization with class taxonomies. arXiv preprint arXiv:2407.16067. Cited by: Appendix A. [60] C. Squires, A. Seigal, S. S. Bhate, and C. Uhler (2023) Linear causal disentanglement via interventions. In International Conference on Machine Learning, Cited by: Appendix A. [61] C. Squires and C. Uhler (2023) Causal structure learning: a combinatorial perspective. Foundations of Computational Mathematics 23 (5), p. 1781–1815. Cited by: §1. [62] A. Subbaswamy and S. Saria (2020) From development to deployment: dataset shift, causality, and shift-stable models in health ai. Biostatistics 21 (2), p. 345–352. Cited by: Appendix J. [63] Y. Sui, Q. Wu, J. Wu, Q. Cui, L. Li, J. Zhou, X. Wang, and X. He (2024) Unleashing the power of graph data augmentation on covariate distribution shift. Advances in Neural Information Processing Systems 36. Cited by: Appendix A. [64] Z. Szabó and B. K. Sriperumbudur (2018) Characteristic and universal tensor product kernels. Journal of Machine Learning Research 18 (233), p. 1–29. Cited by: §C.6, §E.1, Lemma 1. [65] R. Tachet des Combes, H. Zhao, Y. Wang, and G. J. Gordon (2020) Domain adaptation with conditional distribution matching and generalized label shift. Advances in Neural Information Processing Systems 33, p. 19276–19289. Cited by: Appendix J. [66] E. J. T. Tchetgen, A. Ying, Y. Cui, X. Shi, and W. Miao (2020) An introduction to proximal causal learning. arXiv preprint arXiv:2009.10982. Cited by: §1. [67] A. Torralba and A. A. Efros (2011) Unbiased look at dataset bias. In CVPR, Cited by: §7.5. [68] R. Vedantam, D. Lopez-Paz, and D. J. Schwab (2021) An empirical investigation of domain generalization with empirical risk minimizers. In Advances in Neural Information Processing Systems, Cited by: Appendix A. [69] V. Veitch, A. D’Amour, S. Yadlowsky, and J. Eisenstein (2021) Counterfactual invariance to spurious correlations: why and how to pass stress tests. arXiv preprint arXiv:2106.00545. Cited by: Appendix A. [70] H. Venkateswara, J. Eusebio, S. Chakraborty, and S. Panchanathan (2017) Deep hashing network for unsupervised domain adaptation. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 5018–5027. Cited by: §7.5. [71] J. von Kügelgen, M. Besserve, L. Wendong, L. Gresele, A. Kekić, E. Bareinboim, D. Blei, and B. Schölkopf (2023) Nonparametric identifiability of causal representations from unknown interventions. In Advances in Neural Information Processing Systems, Cited by: Appendix A. [72] J. Von Kügelgen, Y. Sharma, L. Gresele, W. Brendel, B. Schölkopf, M. Besserve, and F. Locatello (2021) Self-supervised learning with data augmentations provably isolates content from style. Advances in neural information processing systems 34, p. 16451–16467. Cited by: Appendix J. [73] X. Wang, H. Chen, Z. Wu, W. Zhu, et al. (2024) Disentangled representation learning. IEEE Transactions on Pattern Analysis and Machine Intelligence. Cited by: §1. [74] Z. Wang and V. Veitch (2023) The causal structure of domain invariant supervised representation learning. arXiv preprint arXiv:2208.06987. Cited by: Appendix A, Appendix A, Appendix J, §1, §2.2, §2.2. [75] L. Wendong, A. Kekić, J. von Kügelgen, S. Buchholz, M. Besserve, L. Gresele, and B. Schölkopf (2023) Causal component analysis. Advances in Neural Information Processing Systems 36, p. 32481–32520. Cited by: Appendix A. [76] A. Wong, E. Otles, J. P. Donnelly, A. Krumm, J. McCullough, O. DeTroyer-Cooley, J. Pestrue, M. Phillips, J. Konye, C. Penoza, et al. (2021) External validation of a widely implemented proprietary sepsis prediction model in hospitalized patients. JAMA internal medicine 181 (8), p. 1065–1070. Cited by: Appendix J. [77] Q. Xu, R. Zhang, Y. Zhang, Y. Wang, and Q. Tian (2021) A fourier-based framework for domain generalization. In IEEE/CVF conference on computer vision and pattern recognition, Cited by: Appendix A. [78] N. Yang, K. Zeng, Q. Wu, X. Jia, and J. Yan (2022) Learning substructure invariance for out-of-distribution molecular representations. Advances in Neural Information Processing Systems 35, p. 12964–12978. Cited by: §1. [79] J. R. Zech, M. A. Badgeley, M. Liu, A. B. Costa, J. J. Titano, and E. K. Oermann (2018) Variable generalization performance of a deep learning model to detect pneumonia in chest radiographs: a cross-sectional study. PLoS medicine 15 (11), p. e1002683. Cited by: Appendix J. [80] H. Zhang, H. Singh, M. Ghassemi, and S. Joshi (2023) " Why did the model fail?": attributing model performance changes to distribution shifts. Cited by: Appendix J. [81] J. Zhang, K. Greenewald, C. Squires, A. Srivastava, K. Shanmugam, and C. Uhler (2023) Identifiability guarantees for causal disentanglement from soft interventions. In Advances in Neural Information Processing Systems, Cited by: §1. [82] K. Zhang, J. Peters, D. Janzing, and B. Schölkopf (2012) Kernel-based conditional independence test and application in causal discovery. arXiv preprint arXiv:1202.3775. Cited by: Appendix J. [83] K. Zhang, B. Schölkopf, K. Muandet, and Z. Wang (2013) Domain adaptation under target and conditional shift. In International conference on machine learning, p. 819–827. Cited by: Appendix J. [84] X. Zhang, P. Cui, R. Xu, L. Zhou, Y. He, and Z. Shen (2021) Deep stable learning for out-of-distribution generalization. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, p. 5372–5382. Cited by: Appendix A. [85] H. Zhao, R. T. Des Combes, K. Zhang, and G. Gordon (2019) On learning invariant representations for domain adaptation. In International conference on machine learning, p. 7523–7532. Cited by: Appendix J. [86] K. Zhou, Y. Yang, T. Hospedales, and T. Xiang (2020) Deep domain-adversarial image generation for domain generalisation. In Proceedings of the AAAI conference on artificial intelligence, Vol. 34, p. 13025–13032. Cited by: Appendix A. [87] K. Zhou, Y. Yang, T. Hospedales, and T. Xiang (2020) Learning to generate novel domains for domain generalization. In ECCV, p. 561–578. Cited by: Appendix A. [88] Y. Zhu, Y. He, J. Ma, M. Hu, S. Li, and J. Li (2024) Causal inference with latent variables: recent advances and future prospectives. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, p. 6677–6687. Cited by: §1. Appendix Appendix A Related Work Distribution-Invariant Learning. A large body of work in domain adaptation and domain generalization aims to learn representations that are invariant across environments or distributions. Early non-causal approaches such as (C-)ADA [37] and DDAIG [86] focus on aligning feature distributions across domains. More recent methods, including FeAT [13] and AIA [63], design more expressive objectives to improve robustness under distribution shift. Representative domain generalization methods include ERM [68], DOMAINBED [20], L2A-OT [87], StableNet [84], SagNets [42], SWAD [12], FACT [77], Ratatouille [51], XRM [46], and IRM [4], along with extensions such as REx [30], counterfactual invariance [69], information bottleneck approaches [1], causal domain adaptation [74], and transportable representations [24]. While effective in practice, these methods generally enforce distributional invariance without explicitly reasoning about the underlying causal structure or validating whether the assumed invariance is causally appropriate. Weight-averaging methods such as SWAD [12] modify the optimization trajectory rather than the objective, and are complementary to SaCRL, which determines which invariance is enforced. Intervention-based Causal Learning. Intervention-based causal representation learning leverages environment changes or interventions to identify causal factors. Foundational works include nonparametric ICR [71], which jointly learns representations and intervention targets under SCM assumptions, and general nonlinear mixing models [9], which address nonlinear latent interactions. Several weakly supervised approaches [8, 38, 39, 59] exploit partial intervention information to learn invariant representations, while ICRL [2] extends these ideas to imperfect interventions. Weak distributional invariance [3] studies identifiability under perfect and imperfect interventions by identifying latent variables with stable distributional properties. Unsupervised approaches such as ICA-based causal learning [75] aim to recover latent causal variables through component analysis and statistical independence criteria. AIT [53] assumes explicit SCMs with known DAG structures and focuses on identifying standard causal directions. These methods typically rely on either explicit intervention signals or strong assumptions about the form of the causal graph. Structure-based Causal Representation Learning. Structure-based causal representation learning methods explicitly assume or model a particular causal structure. DISRL [74] and Causal Disentanglement [60] learn representations aligned with latent causal graphs, while classical approaches such as ICP [44] and its nonlinear extension [21] exploit invariance across environments to identify causal mechanisms. Related works such as CSG [34] and LECI [19] focus on identifying causal subgraphs and removing spurious correlations, often under DAG assumptions. KCDC [41] employs kernel-based independence tests for causal discovery and orientation. Anti-causal modeling approaches [40] study settings where inputs are generated from labels and protected attributes using conventional causal models and do-calculus, while [7] develops measure-theoretic formulations for anti-causal representation learning. Research Gap. Existing approaches face a trilemma. Structure-based methods achieve strong guarantees but require knowing the causal graph. Invariant learning methods avoid explicit graphs but hardcode implicit structural assumptions that cause failures under mismatch. Causal discovery methods learn structure but decouple it from representation learning, preventing joint optimization. No existing framework simultaneously discovers causal structure and learns representations adapted to that structure in a unified, end-to-end differentiable objective. Our work fills this gap: we introduce a structure-agnostic objective that evaluates representations against all candidate causal structures, automatically concentrating on the correct one through adaptive weighting, while jointly optimizing predictive performance. Appendix B Notation Table 8 summarizes notation used throughout. We organize symbols into six categories: spaces and sets defining the mathematical objects; random variables and observed data; causal graph components encoding structural assumptions; kernels and independence testing machinery; learning components including representations, predictors, and losses; and theoretical quantities appearing in our analysis. Lowercase denotes realizations (x,yx,y), uppercase denotes random variables (X,YX,Y), and calligraphic denotes spaces (,X,Y). Subscripts index environments (PeP_e) or structures (kG_k); superscripts denote components (XZ,XZ⟂X_Z,X_Z ) or iterates (ϕt _t). Table 8: Notation summary. Spaces & Sets Random Variables & Data ,,X,Y,Z Input, output, representation spaces X,Y,EX,Y,E Features, target, environment ℰE Environment set e1,…,em\e_1,…,e_m\ U,ZU,Z Latent confounder, mediator (Ω,ℱ,P)( , F,P) Probability space XZ,XZ⟂X_Z,X_Z Mediated, orthogonal components ℋ,ℋH,H_X RKHS, RKHS on X DeD_e Dataset (xje,yje)j=1ne\(x_j^e,y_j^e)\_j=1^n_e Δ2 ^2 Probability simplex n,m,dn,m,d Samples, environments, dimension Causal Structure Kernels & Independence =(V,E)G=(V,E) DAG with vertices, edges k,k,k_X Reproducing kernel 1,2,3G_1,G_2,G_3 Anti-causal, conf-desc, conf-out κ Kernel bound: k(x,x)≤κ2k(x,x)≤κ^2 ∗=k∗G^*=G_k^* True structure μP _P Mean embedding of P Pa(v)Pa(v) Parents of v XYC_XY Cross-covariance operator KS,KSℚK_S,K_S^Q Causal, interventional kernel HSIC(X,Y)HSIC(X,Y) Independence criterion ℐkI_k Invariance for structure k X~,Y~ X, Y Regression residuals Learning Components Theoretical Quantities ϕ,wφ,w Representation, predictor D Diversity matrix (m×3m× 3) Φ, ,W Function classes δmin _ Min diversity parameter ℓ,ℛe,ℛ^e ,R_e, R_e Loss, population/empirical risk η Faithfulness gap Vk,Vk∗V_k,V_k^* Empirical/population violation d(e,e′)d_K(e,e ) Interventional distance αk _k Adaptive weight ∝e−βVk e^-β V_k ξn _n Estimation error λ,β,ρ,γλ,β,ρ,γ Regularization, temp, CVaR, ridge LwL_w Lipschitz constant Appendix C Mathematical Preliminaries C.1 Measure-Theoretic Foundations Definition C.1 (Measurable Space). A pair (Ω,ℱ)( , F) is a measurable space where Ω is a non-empty set and ℱ F is a σ-algebra on Ω , i.e., a collection of subsets satisfying: (i) ∅,Ω∈ℱ , ∈ F; (i) A∈ℱ⇒Ac∈ℱA∈ F A^c∈ F; (i) Aii=1∞⊂ℱ⇒⋃i=1∞Ai∈ℱ\A_i\_i=1^∞⊂ F _i=1^∞A_i∈ F. Definition C.2 (Probability Measure). A probability measure on a measurable space (Ω,ℱ)( , F) is a function μ:ℱ→[0,1]μ: F→[0,1] satisfying μ(Ω)=1μ( )=1. For any countable collection Aii=1∞\A_i\_i=1^∞ of pairwise disjoint sets in ℱ F, μ(⋃i=1∞Ai)=∑i=1∞μ(Ai)μ ( _i=1^∞A_i )= _i=1^∞μ(A_i) (14) The triple (Ω,ℱ,μ)( , F,μ) is a probability space. Definition C.3 (Product Space). Given measurable spaces (Ω1,ℱ1)( _1, F_1) and (Ω2,ℱ2)( _2, F_2), the product σ-algebra ℱ1⊗ℱ2 F_1 F_2 on Ω1×Ω2 _1× _2 is the σ-algebra generated by measurable rectangles A1×A2A_1× A_2 where A1∈ℱ1A_1∈ F_1 and A2∈ℱ2A_2∈ F_2. For probability measures μ1 _1 on (Ω1,ℱ1)( _1, F_1) and μ2 _2 on (Ω2,ℱ2)( _2, F_2), the product measure μ1⊗μ2 _1 _2 on (Ω1×Ω2,ℱ1⊗ℱ2)( _1× _2, F_1 F_2) is uniquely determined by (μ1⊗μ2)(A1×A2)=μ1(A1)⋅μ2(A2)( _1 _2)(A_1× A_2)= _1(A_1)· _2(A_2) for all A1∈ℱ1A_1∈ F_1, A2∈ℱ2A_2∈ F_2. Definition C.4 (Markov Kernel). A (Markov) kernel from a measurable space (Ω1,ℱ1)( _1, F_1) to a measurable space (Ω2,ℱ2)( _2, F_2) is a mapping K:Ω1×ℱ2→[0,1]K: _1× F_2→[0,1]. For each fixed ω∈Ω1ω∈ _1, the function A↦K(ω,A)A K(ω,A) is a probability measure on (Ω2,ℱ2)( _2, F_2). For each fixed A∈ℱ2A∈ F_2, the function ω↦K(ω,A)ω K(ω,A) is ℱ1 F_1-measurable. Intuitively, K(ω,A)K(ω,A) represents the conditional probability of event A∈ℱ2A∈ F_2 given state ω∈Ω1ω∈ _1. C.2 Causal Structures and Directed Acyclic Graphs Definition C.5 (Directed Acyclic Graph). A directed acyclic graph (DAG) is a pair =(V,E)G=(V,E) where V is a finite set of vertices (nodes) and E⊆V×VE V× V is a set of directed edges such that there exists no sequence of vertices v1,…,vkv_1,…,v_k with k≥2k≥ 2, v1=vkv_1=v_k, and (vi,vi+1)∈E(v_i,v_i+1)∈ E for all i=1,…,k−1i=1,…,k-1. Definition C.6 (Causal Structure). A causal structure is a tuple =(V,E,ℱvv∈V,Kvv∈V)S=(V,E,\ F_v\_v∈ V,\K_v\_v∈ V) where (V,E)(V,E) is a DAG with vertices representing random variables. For each v∈Vv∈ V, (Ωv,ℱv)( _v, F_v) is a measurable space. For each v∈Vv∈ V, KvK_v is a kernel from (∏u∈Pa(v)Ωu,⨂u∈Pa(v)ℱu) ( _u (v) _u, _u (v) F_u ) to (Ωv,ℱv)( _v, F_v) representing the causal mechanism by which parents generate the variable v. For root nodes with Pa(v)=∅Pa(v)= , the kernel KvK_v reduces to a probability measure. We now formalize how observed features arise from multiple generating mechanisms. Definition C.7 (Causal Feature Decomposition). Let Z be a latent variable in the causal graph G. The observed features X admit a causal decomposition with respect to Z if there exist measurable functions πZ:→Z _Z:X _Z and πZ⟂:→Z⟂ _Z :X _Z such that X=(πZ(X),πZ⟂(X))=:(XZ,XZ⟂)X=( _Z(X), _Z (X))=:(X_Z,X_Z ). XZX_Z contains all and only the information in X that is causally downstream of Z, and XZ⟂X_Z contains the remaining information in X that is not mediated by Z. We write X=(XZ,XZ⟂)X=(X_Z,X_Z ) to denote this decomposition, where ≅Z×Z⟂X _Z×X_Z . C.3 Interventions and the Do-Calculus Interventions distinguish causal relationships from statistical associations. We formalize both perfect (hard) and imperfect (soft) interventions within the generalized framework. Definition C.8 (Perfect and Imperfect Interventions). An intervention on a causal structure S is specified by a measurable mapping ℚ:ℱ×Ω→[0,1]Q: F× →[0,1] where ℱ F is a suitable σ-algebra. (i) A perfect (hard) intervention on variable V replaces its causal mechanism KVK_V with a fixed distribution ℚQ independent of parent values: ℚ(A∣ω)=ℚ(A)Q(A ω)=Q(A) for all ω∈Ωω∈ , A∈ℱA∈ F (i) An imperfect (soft) intervention on variable V modifies its causal mechanism to depend on both parents and intervention parameters: ℚ(A∣ω)Q(A ω) varies with ω∈Ωω∈ The do-operator do(X=x)do(X=x) corresponds to a perfect intervention that sets X to value x by replacing the kernel KXK_X with a point mass at x, formally: ℚ(A)=x(A)Q(A)=I_\x\(A) for measurable sets A. Definition C.9 (Causal Independence). Random variables X and Y are causally independent given Z, denoted X⟂cY∣ZX \!\!\! _cY Z, if P(Y∣do(X=x),Z=z)=P(Y∣Z=z)P(Y (X=x),Z=z)=P(Y Z=z) (15) for all x in the support of X and all z in the support of Z, where do(X=x)do(X=x) denotes an intervention setting X to value x. The do-operator do(X=x)do(X=x) differs fundamentally from P(Y∣X=x)P(Y X=x) which observes X=xX=x while preserving causal relationships, whereas P(Y∣do(X=x))P(Y (X=x)) forces X=xX=x by intervention. C.4 Measure-Theoretic Causality The causal-space, product-causal-space, and interventional-kernel constructions below follow Behnam and Wang [7]; we extend them with structure-indexed kernel families (Definition C.14). Definition C.10 (Sub-σ-algebra). Given a product causal space (Ω,ℱ,ℙ,,)( , F,P,K,G) with index set T, for any subset S⊆TS T, the sub-σ-algebra ℱS F_S is the σ-algebra generated by measurable rectangles of the form ∏t∈SAt×∏t∈T∖SΩt _t∈ SA_t× _t∈ T S _t where At∈tA_t∈ A_t for t∈St∈ S. Intuitively, ℱS F_S captures events that depend only on the components indexed by S. Definition C.11 (Causal Kernel). A causal kernel for index set S⊆TS T is a function KS:Ω×ℱ→[0,1]K_S: × F→[0,1]. For each fixed ω∈Ωω∈ , the mapping A↦KS(ω,A)A K_S(ω,A) is a probability measure on (Ω,ℱ)( , F). For each fixed A∈ℱA∈ F, the mapping ω↦KS(ω,A)ω K_S(ω,A) is ℱS F_S-measurable. Definition C.12 (Causal Space). For an environment e∈ℰe , a causal space is a tuple (Ωe,ℱe,Pe,e,e)( _e, F_e,P_e,K_e,G_e) where Ωe=∏v∈VeΩv _e= _v∈ V_e _v is the sample space with product structure, and VeV_e is the vertex set of the causal DAG and each Ωv _v is the state space of variable v. Here, ℱe=⨂v∈Veℱv F_e= _v∈ V_e F_v is the product σ-algebra, where ℱv F_v is the σ-algebra on Ωv _v, PeP_e is a probability measure on (Ωe,ℱe)( _e, F_e) compatible with eG_e, e=Kv(e)v∈VeK_e=\K_v^(e)\_v∈ V_e is the family of causal mechanisms. Each Kv(e)K_v^(e) is a Markov kernel (Definition C.4) from the parent space (∏u∈Pa(v)Ωu,⨂u∈Pa(v)ℱu) ( _u (v) _u, _u (v) F_u ) to (Ωv,ℱv)( _v, F_v), encoding the conditional distribution Pe(V=⋅∣Pa(v)=⋅)P_e(V=· (v)=·), and e=(Ve,Ee)G_e=(V_e,E_e) is the causal DAG specifying the qualitative causal structure. The joint distribution factorizes according to the DAG: Pe(ω)=∏v∈VeKv(e)(ωPa(v),ωv)P_e(ω)= _v∈ V_eK_v^(e) ( _Pa(v),\ _v\ ) (16) where ω=(ωv)v∈Ve∈Ωeω=( _v)_v∈ V_e∈ _e and ωPa(v) _Pa(v) denotes the components corresponding to parents of v. Definition C.13 (Product Causal Space). Given causal spaces (Ωei,ℱei,Pei,Kei,ei)i=1m\( _e_i, F_e_i,P_e_i,K_e_i,G_e_i)\_i=1^m for environments eii=1m\e_i\_i=1^m, the product causal space is the tuple (Ω,ℱ,ℙ,,)( , F,P,K,G) where Ω=∏i=1mΩei = _i=1^m _e_i is the product sample space, ℱ=⨂i=1mℱei F= _i=1^m F_e_i is the product σ-algebra, ℙ=⨂i=1mPeiP= _i=1^mP_e_i is the product measure, =KS:S∈(T)K=\K_S:S∈ P(T)\ is a family of causal kernels indexed by subsets S of the combined index set T=⋃i=1mTeiT= _i=1^mT_e_i, and =(V,E)G=(V,E) is a causal DAG compatible with the individual environment graphs. The product causal space admits sub-σ-algebras and causal kernels indexed by variable subsets; see Appendix C for formal definitions. Definition C.14 (Structure-Agnostic Causal Space). A structure-agnostic causal space extends the product causal space (Definition C.13) with structure-indexed components: (Ω,ℱ,ℙ,kk=13,kk=13)( , F,P,\K_k\_k=1^3,\G_k\_k=1^3) where (Ω,ℱ,ℙ)( , F,P) is the product space over environments, kk=13\G_k\_k=1^3 are the three candidate causal DAGs (Definition 2.1), and k=KS(k):S⊆TK_k=\K_S^(k):S T\ is the causal kernel family encoding how interventions propagate under structure kG_k. The ℱS F_S-measurability condition ensures KS(ω,A)K_S(ω,A) depends only on components indexed by S, encoding graph-specified causal dependencies. Definition C.15 (Interventional Kernel). Let (Ω,ℱ,ℙ,,)( , F,P,K,G) be a product causal space with causal kernel KSK_S for index set S⊆TS T. Let ℚ:Ω×ℱ→[0,1]Q: × F→[0,1] be an intervention kernel, where ℚ(ω′,A)Q(ω ,A) specifies the probability of event A when intervening at state ω′ω . The interventional kernel KSℚ:Ω×ℱ→[0,1]K_S^Q: × F→[0,1] is defined by: KSℚ(ω,A):=∫Ωℚ(ω′,A)KS(ω,dω′)K_S^Q(ω,A):= _ Q(ω ,A)\,K_S(ω,dω ) (17) where ω∈Ωω∈ is the conditioning state (values of variables in S), KS(ω,dω′)K_S(ω,dω ) is the original causal kernel giving the distribution over ω′ω given ω, and ℚ(ω′,A)Q(ω ,A) is the intervention’s effect on event A at state ω′ω . This construction unifies perfect and imperfect interventions: when ℚ(A∣ω′)=ℚ(A)Q(A ω )=Q(A) is constant, the interventional kernel corresponds to a hard intervention; when ℚ(A∣ω′)Q(A ω ) varies with ω′ω , it represents a soft intervention. For multi-level representations, we consider compositions ϕ=ϕH∘ϕLφ= _H _L where ϕL:→L _L:X _L is the low-level representation mapping inputs to an intermediate feature space, and ϕH:L→H _H:Z_L _H is the high-level representation mapping intermediate features to a refined representation space. Definition C.16 (Sufficiency and Minimality). A representation ϕφ is sufficient for predicting Y if there exists a predictor f:→f:Z such that (f∘ϕ)(X)(f φ)(X) achieves the Bayes-optimal prediction error for Y. It is minimal sufficient if it is sufficient and for any other sufficient representation ϕ′:→′φ :X , there exists a measurable function g:′→g:Z such that ϕ=g∘ϕ′φ=g φ . C.5 Kernel Independence Framework Definitions Definition C.17 (Reproducing Kernel Hilbert Space). Let X be a non-empty set. A Hilbert space ℋH of functions f:→ℝf:X is a reproducing kernel Hilbert space (RKHS) if there exists a function k:×→ℝk:X×X such that: 1. k(⋅,x)∈ℋk(·,x) for all x∈x 2. f(x)=⟨f,k(⋅,x)⟩ℋf(x)= f,k(·,x) _H for all f∈ℋf , x∈x (reproducing property) Definition C.18 (Cross-Covariance Operator). Let (,ℋ,k)(X,H_X,k_X) and (,ℋ,k)(Y,H_Y,k_Y) be RKHS pairs. The cross-covariance operator XY:ℋ→ℋC_XY:H_Y _X is defined by XY:=XY[(k(⋅,X)−μX)⊗(k(⋅,Y)−μY)]C_XY:=E_XY[(k_X(·,X)- _X) (k_Y(·,Y)- _Y)] (18) where ⊗ denotes the tensor product. C.6 Assumptions Here we state assumptions required for theoretical findings. Assumption 2 (Characteristic Kernels). The kernels k:×→ℝk_X:X×X , k:×→ℝk_Y:Y×Y , k:×→ℝk_Z:Z×Z , and kℰ:ℰ×ℰ→ℝk_E:E×E are characteristic and bounded: supxk(x,x)≤κ2<∞ _xk(x,x)≤κ^2<∞. Characteristic kernels ensure HSIC fully characterizes independence (Lemma 1). We further require the tensor products of these kernels used in HSIC to be characteristic, which holds for Gaussian and delta kernels [64]. For the conditional HSIC tests underlying V1V_1 and V2V_2, we additionally require the regression operator to be well-behaved. Assumption 3 (Kernel Regularity). The kernel regression operator satisfies: 1. The conditional expectation operator [⋅∣Z]E[· Z] lies in the RKHS ℋH_Z. 2. The eigenvalues λj\ _j\ of the kernel integral operator satisfy ∑jλj−1/2<∞ _j _j^-1/2<∞. Beyond kernel conditions, structure identification requires environments that vary in distinguishable ways. We formalize this via a diversity matrix measuring distributional shifts across environments. Definition C.19 (Environment Diversity Matrix). Let ℰ=e1,…,emE=\e_1,…,e_m\ be the training environments and P¯ P the pooled distribution. Center each kernel mean embedding at its average over environments: μ¯X:=1m∑eμPe(X) μ_X:= 1m _e _P_e(X), μ¯X|Y=y:=1m∑eμPe(X|Y=y) μ_X|Y=y:= 1m _e _P_e(X|Y=y), and μ¯Y|X=x:=1m∑eμPe(Y|X=x) μ_Y|X=x:= 1m _e _P_e(Y|X=x). The m×3m× 3 diversity matrix D has entries De1 D_e1 :=‖μPe(X)−μ¯X‖ℋ :=\| _P_e(X)- μ_X\|_H (marginal shift), (marginal shift), (19) De2 D_e2 :=(y∼P¯(Y)‖μPe(X|Y=y)−μ¯X|Y=y‖ℋ2)1/2 := (E_y P(Y)\| _P_e(X|Y=y)- μ_X|Y=y\|_H^2 )^1/2 (conditional shift), (conditional shift), De3 D_e3 :=(x∼P¯(X)‖μPe(Y|X=x)−μ¯Y|X=x‖ℋ2)1/2 := (E_x P(X)\| _P_e(Y|X=x)- μ_Y|X=x\|_H^2 )^1/2 (outcome shift). (outcome shift). For characteristic kernels, column j of D is zero if and only if the corresponding distribution is identical across environments (P¯ P-almost everywhere for j=2,3j=2,3). The matching invariance is then untestable: a zero first column makes ℐ3I_3 untestable, a zero second column ℐ1I_1, and a zero third column ℐ2I_2. Assumption 1 requires every column to be nonzero, quantified by the minimum diversity parameter δmin:=minj∈1,2,3maxe∈ℰDej>0. _ := _j∈\1,2,3\ _e D_ej>0. (20) Full column rank, rank()=3rank(D)=3, is sufficient for δmin>0 _ >0 but not necessary. Assumption 4 (Faithfulness). The joint distribution P(X,Y,E)P(X,Y,E) is faithful to the true causal structure ∗G^*: all and only the conditional independencies implied by ∗G^* hold in P. Faithfulness ensures that violations of wrong invariances are detectable. Without it, a distribution could satisfy independencies not implied by the graph, making structures indistinguishable. The remaining assumptions concern the function classes used for prediction and representation. Assumption 5 (Predictor Lipschitz Continuity). The predictor class W consists of LwL_w-Lipschitz functions: for all w∈w and z,z′∈z,z , ‖w(z)−w(z′)‖2≤Lw‖z−z′‖2\|w(z)-w(z )\|_2≤ L_w\|z-z \|_2 (21) with Lw≤Lmax<∞L_w≤ L_ <∞ uniformly over W. This Lipschitz condition controls how prediction errors propagate through distribution shift, appearing directly in our generalization bound (Theorem D.5). Assumption 6 (Representation Class Regularity). The representation class Φ is parameterized by θ∈Θ⊆ℝpθ∈ ^p where: 1. Θ is compact and convex, 2. ϕθ(x) _θ(x) is continuously differentiable in θ, 3. The violation functions Vk(ϕθ)V_k( _θ) are LVL_V-smooth in θ. These regularity conditions are standard for neural network function classes and ensure convergence of gradient-based optimization. Assumption 1 can be checked before training: (1) estimate D with empirical (conditional) mean embeddings; (2) for each column j, compare maxeD^ej _e D_ej with its permutation null obtained by shuffling environment labels; (3) if column j is not significantly nonzero, remove from the candidate set the structure whose invariance it tests (j=1j=1: 3G_3, j=2j=2: 1G_1, j=3j=3: 2G_2); (4) if no column is significant, the environments carry no information about the structure, and ERM is the appropriate fallback. Appendix D Theoretical Results D.1 Violation Characterization, Structure Identification & Representation Optimality The key insight is that under the true structure, the corresponding invariance is achievable without sacrificing predictive accuracy, while wrong invariances are not. We introduce between definitions and theorems. Definition D.1 (ϵε-Predictive Representation). A representation ϕ:→φ:X is ϵε-predictive if: infw∈maxe∈ℰℛe(w∘ϕ)≤ℛ∗+ϵ _w _e R_e(w φ) ^*+ε (22) where ℛ∗:=inff:→maxeℛe(f)R^*:= _f:X _eR_e(f) is the Bayes-optimal worst-case risk. Denote the class of ϵε-predictive representations as Φϵ _ε. We define the base predictive class Φ0:=ϕ:ϕ is ϵ0-predictive _0:=\φ:φ is _0-predictive\ with ϵ0:=δmin/8 _0:= _ /8 (a fixed constant determined by environment diversity). One needs to measure the minimum violation incurred by any predictive representation when the wrong invariance is imposed. Definition D.2 (Faithfulness Gap). Let ∗=k∗G^*=G_k^* be the true causal structure. The faithfulness gap is: η:=mininfϕ∈Φ0k≠k∗Vk∗(ϕ)η:= _k≠ k^* _φ∈ _0V_k^*(φ) (23) Under faithfulness (Assumption 4), η>0η>0. Larger η implies easier structure identification; η→0η→ 0 indicates near-Markov equivalence between structures. Theorem D.1 (Violation Characterization). Under Assumptions 2–4, for the class of ϵε-predictive representations Φϵ _ε: 1. Achievability: If ∗=kG^*=G_k, then infϕ∈ΦϵVk∗(ϕ)=0 _φ∈ _εV_k^*(φ)=0. 2. Non-achievability: If ∗≠kG^* _k, then infϕ∈ΦϵVk∗(ϕ)≥η/2>0 _φ∈ _εV_k^*(φ)≥η/2>0. The constraint ϕ∈Φϵφ∈ _ε prevents trivial solutions: a constant representation ϕ(x)=cφ(x)=c achieves Vk=0V_k=0 for all k, but has ℛe(w∘ϕ)=ℛtrivial≫ℛ∗+ϵR_e(w φ)=R_trivial ^*+ε, hence ϕ∉Φϵφ∉ _ε. The joint objective ℒL enforces this implicitly via the prediction loss. The asymmetry established in Theorem D.1 (achievable invariance for the true structure versus irreducible violations for wrong structures) directly enables consistent identification from finite samples. Theorem D.2 (Structure Identification). Under Assumptions 2–6, let ϕn _n minimize ℒ(ϕ,w)L(φ,w) with n samples per environment. Define the estimated structure: k^n:=argminkVk(ϕn) k_n:= _kV_k( _n) (24) Then for the true structure ∗=k∗G^*=G_k^*: ℙ(k^n=k∗)≥1−6exp(−nη22048κ4)P( k_n=k^*)≥ 1-6 (- nη^22048κ^4 ) (25) where η is the faithfulness gap (Definition D.2) and κ is the kernel bound (Assumption 2). Corollary 1 (Identification with Random Features). Suppose D random features (Definition F.1) achieve uniform kernel error εD _D with probability at least 1−δ′1-δ , so that |HSIC^RFF−HSIC^|≤Cκ2εD| HSIC_RFF- HSIC|≤ Cκ^2 _D. Under the conditions of Theorem D.2, if Cκ2εD≤η/16Cκ^2 _D≤η/16, then ℙ(k^nRFF=k∗)≥1−6exp(−nη28192κ4)−δ′,P( k_n^RFF=k^*)≥ 1-6 (- nη^28192κ^4 )-δ , (26) which holds for D=Ω(κ4η−2dlog(n/δ′))D= (κ^4η^-2d (n/δ )) [50]. The required D grows as η−2η^-2: near-equivalent structures need more random features, not only more samples. Corollary 2 (Adaptive Weight Concentration). Under the conditions of Theorem D.2, the adaptive weights satisfy: αk∗(ϕn)≥1−2exp(−βη/2) _k^*( _n)≥ 1-2 (-βη/2) (27) for n sufficiently large, where η is the faithfulness gap. Via established structure identification, we now quantify how well the learned representation achieves the corresponding invariance and characterize properties of global optima. Theorem D.3 (Invariance Achievement). Let ϕn _n minimize ℒ(ϕ,w)L(φ,w) with regularization λ>0λ>0. Under Assumptions 2–6, for the true structure ∗=k∗G^*=G_k^*: Vk∗(ϕn)≤C1λ⋅maxeℛ^e(w∗∘ϕ∗)+C2n−1/2V_k^*( _n)≤ C_1λ· _e R_e(w^* φ^*)+C_2n^-1/2 (28) where ϕ∗∈Φϵφ^*∈ _ε achieves infϕ∈ΦϵVk∗(ϕ)=0 _φ∈ _εV_k^*(φ)=0 and C1,C2>0C_1,C_2>0 are constants depending on κ and the function class complexity. The violation decreases at rate O(n−1/2)O(n^-1/2) and can be made arbitrarily small by increasing λ, though at the cost of prediction accuracy. When the function classes are sufficiently expressive, the global minimum achieves exact invariance. Theorem D.4 (Global Optimality). If the representation class Φ includes the true invariant representation ϕ∗φ^* satisfying ℐk∗I_k^*, and the predictor class W includes the Bayes-optimal predictor, then the global minimum of ℒL achieves: 1. Vk∗(ϕ)=0V_k^*(φ)=0 (exact invariance) 2. maxeℛe(w∘ϕ)=ℛ∗ _eR_e(w φ)=R^* (Bayes risk) We know the objective learns the optimal invariant representation when the function classes are well-specified, with finite-sample guarantees that degrade under approximation. D.2 Generalization Bounds We now bound out-of-distribution (OOD) generalization by quantifying how distribution shift affects test performance. To measure the similarity between environments, we introduce a causal distance metric below. Definition D.3 (Interventional Kernel Distance). For environments e,e′∈ℰe,e , the interventional kernel distance is: d(e,e′):=supS⊆T‖KS(e)−KS(e′)‖opd_K(e,e ):= _S T\|K_S^(e)-K_S^(e )\|_op (29) where KS(e)K_S^(e) is the causal kernel for environment e and ∥⋅∥op\|·\|_op is the operator norm. For a test environment, define d(etest,ℰtrain):=mine∈ℰtraind(etest,e)d_K(e_test,E_train):= _e _traind_K(e_test,e). This distance captures the maximal change in causal mechanisms across environments. Using this metric, we decompose test error into training error, distribution shift, and statistical estimation error. Theorem D.5 (Generalization Bound). Let (ϕn,wn)( _n,w_n) minimize ℒL over training environments ℰtrainE_train. For a test environment eteste_test, under Assumptions 2–6: ℛetest(wn∘ϕn)≤maxe∈ℰtrainℛ^e(wn∘ϕn)⏟training error+Lw⋅d(etest,ℰtrain)⏟shift penalty+ξn⏟estimation error _e_test(w_n _n)≤ _e _train R_e(w_n _n)_training error+ L_w· d_K(e_test,E_train)_shift penalty+ _n_estimation error (30) where the estimation error satisfies: ξn=O((Φ,)+log(1/δ)n)+Cη⋅Vk∗(ϕn) _n=O ( C( ,W)+ (1/δ)n )+C_η· V_k^*( _n) (31) with (Φ,)C( ,W) the function class complexity and Cη>0C_η>0 a constant depending on the faithfulness gap. The bound holds with probability at least 1−δ1-δ. The invariance violation Vk∗(ϕn)V_k^*( _n) directly contributes to generalization error, motivating the joint optimization of prediction and invariance. Combining all preceding results yields asymptotic optimality. Corollary 3 (Asymptotic Optimality). Under the conditions of Theorems D.2–D.5, as n→∞n→∞: ℛetest(wn∘ϕn)→infϕ∈Φk∗ℛetest(w∗(ϕ)∘ϕ)R_e_test(w_n _n) p _φ∈ _k^*R_e_test(w^*(φ) φ) (32) where Φk∗=ϕ:ℐk∗(ϕ) holds _k^*=\φ:I_k^*(φ) holds\ is the class of representations satisfying the true invariance. D.3 Sample Complexity Theorem D.6 (Sample Complexity). To achieve structure identification with probability ≥1−δ≥ 1-δ and generalization error ≤ϵ≤ε, the required sample size per environment is: n=O(κ4η2log12δ+κ4ϵ2log|Φ|δ)n=O ( κ^4η^2 12δ+ κ^4ε^2 | |δ ) (33) where |Φ|| | is the covering number of the representation class. Appendix E Proofs This section contains complete proofs of all theoretical results. We organize proofs by building from foundational kernel independence results to the main structure identification and generalization theorems. We first establish auxiliary concentration results that underpin our main theorems. Lemma 4 (HSIC Concentration). Let (Xi,Yi)i=1n(X_i,Y_i)_i=1^n be i.i.d. samples from PXYP_XY. Under Assumption 2 with kernel bound κ, the empirical HSIC satisfies: ℙ(|HSIC^(X,Y)−HSIC(X,Y)|>ϵ)≤2exp(−nϵ232κ4)P ( | HSIC(X,Y)-HSIC(X,Y) |>ε )≤ 2 (- nε^232κ^4 ) (34) for all ϵ>0ε>0. Proof. The empirical HSIC can be written as a V-statistic: HSIC^(X,Y)=1n2∑i,j=1nh((Xi,Yi),(Xj,Yj)) HSIC(X,Y)= 1n^2 _i,j=1^nh((X_i,Y_i),(X_j,Y_j)) (35) where the kernel h is defined by: h((xi,yi),(xj,yj)) h((x_i,y_i),(x_j,y_j)) =k(xi,xj)k(yi,yj)+[k(X,X′)][k(Y,Y′)] =k_X(x_i,x_j)k_Y(y_i,y_j)+E[k_X(X,X )]E[k_Y(Y,Y )] (36) −k(xi,xj)[k(Y,Y′)]−[k(X,X′)]k(yi,yj) -k_X(x_i,x_j)E[k_Y(Y,Y )]-E[k_X(X,X )]k_Y(y_i,y_j) (37) Under Assumption 2, we have |k(x,x′)|≤κ2|k_X(x,x )|≤κ^2 and |k(y,y′)|≤κ2|k_Y(y,y )|≤κ^2. Therefore: |h((xi,yi),(xj,yj))|≤4κ4|h((x_i,y_i),(x_j,y_j))|≤ 4κ^4 (38) We apply McDiarmid’s inequality. Changing a single observation (Xi,Yi)(X_i,Y_i) affects at most 2n2n terms in the double sum, each bounded by 4κ44κ^4, giving a bounded difference of: ci=2n⋅4κ4n2=8κ4nc_i= 2n· 4κ^4n^2= 8κ^4n (39) By McDiarmid’s inequality: ℙ(|HSIC^−[HSIC^]|>ϵ)≤2exp(−2ϵ2∑i=1nci2)=2exp(−2ϵ2n⋅64κ8/n2)P ( | HSIC-E[ HSIC] |>ε )≤ 2 (- 2ε^2 _i=1^nc_i^2 )=2 (- 2ε^2n· 64κ^8/n^2 ) (40) Since the empirical HSIC is an asymptotically unbiased estimator with |[HSIC^]−HSIC|=O(n−1)|E[ HSIC]-HSIC|=O(n^-1), absorbing this bias into the concentration bound yields: ℙ(|HSIC^−HSIC|>ϵ)≤2exp(−nϵ232κ4)P ( | HSIC-HSIC |>ε )≤ 2 (- nε^232κ^4 ) (41) ∎ Lemma 5 (Kernel Regression Consistency). Under Assumptions 2 and 3, let m^(z)=z⊤(K+γI)−1 m(z)=k_z (K+γ I)^-1y be the kernel ridge regression estimate. For regularization γ=n−1/(2+d)γ=n^-1/(2+d) where d is the effective dimension: [∥m^−m∗∥L2(PZ)2]=O(n−2/(2+d))E [\| m-m^*\|_L^2(P_Z)^2 ]=O (n^-2/(2+d) ) (42) where m∗(z)=[Y|Z=z]m^*(z)=E[Y|Z=z] is the true conditional expectation. Proof. The kernel ridge regression estimator admits the bias-variance decomposition: [‖m^−m∗‖2]=‖(I+γ−1)−1m∗−m∗‖2⏟bias2+σ2ntr((K+γI)−1K(K+γI)−1)⏟varianceE [\| m-m^*\|^2 ]= \|(I+ ^-1)^-1m^*-m^*\|^2_bias^2+ σ^2ntr ((K+γ I)^-1K(K+γ I)^-1 )_variance (43) where T is the integral operator with kernel k. Under Assumption 3, the eigenvalues λj\ _j\ of T satisfy the decay condition. The bias term satisfies: bias2≤γ2‖m∗‖ℋ2bias^2≤γ^2\|m^*\|_H^2 (44) For the variance, using the eigenvalue decay: variance≤σ2n∑j=1∞λj(λj+γ)2≤σ2nγvariance≤ σ^2n _j=1^∞ _j( _j+γ)^2≤ σ^2nγ (45) Balancing bias and variance by setting γ2=(nγ)−1γ^2=(nγ)^-1, i.e., γ=n−1/3γ=n^-1/3 for the general case, or γ=n−1/(2+d)γ=n^-1/(2+d) accounting for effective dimension d: [∥m^−m∗∥2]=O(γ2)+O((nγ)−1)=O(n−2/(2+d))E [\| m-m^*\|^2 ]=O(γ^2)+O((nγ)^-1)=O (n^-2/(2+d) ) (46) ∎ E.1 Proof of Lemma 1 (HSIC Characterizes Independence) Proof. We prove both directions. (⇐ ) Independence implies zero HSIC: Assume X⟂YX \!\!\! Y, so PXY=PX⊗PYP_XY=P_X P_Y. The cross-covariance operator is: XY _XY =XY[ϕ(X)⊗ψ(Y)]−μX⊗μY =E_XY[φ(X) ψ(Y)]- _X _Y (47) where ϕ(x)=k(⋅,x)φ(x)=k_X(·,x) and ψ(y)=k(⋅,y)ψ(y)=k_Y(·,y) are the canonical feature maps. Under independence: XY[ϕ(X)⊗ψ(Y)] _XY[φ(X) ψ(Y)] =X[ϕ(X)]⊗Y[ψ(Y)]=μX⊗μY =E_X[φ(X)] _Y[ψ(Y)]= _X _Y (48) Therefore XY=0C_XY=0, and HSIC(X,Y)=‖XY‖HS2=0HSIC(X,Y)=\|C_XY\|_HS^2=0. (⇒ ) Zero HSIC implies independence: Assume HSIC(X,Y)=0HSIC(X,Y)=0. By definition: ‖XY‖HS2=XX′YY′[k(X,X′)k(Y,Y′)] \|C_XY\|_HS^2=E_X Y [k_X(X,X )k_Y(Y,Y )] −2XYX′Y′[k(X,X′)k(Y,Y′)] -2E_XYX Y [k_X(X,X )k_Y(Y,Y )] (49) +XX′[k(X,X′)]YY′[k(Y,Y′)] +E_X [k_X(X,X )]E_Y [k_Y(Y,Y )] where (X,Y)(X,Y) and (X′,Y′)(X ,Y ) are independent copies. This equals the squared MMD between PXYP_XY and PX⊗PYP_X P_Y in the product RKHS ℋ⊗ℋH_X _Y: HSIC(X,Y)=MMD2(PXY,PX⊗PY)HSIC(X,Y)=MMD^2(P_XY,P_X P_Y) (50) Since k⊗k_X k_Y is characteristic (by assumption; e.g., for Gaussian and delta kernels [64]), MMD2=0MMD^2=0 implies PXY=PX⊗PYP_XY=P_X P_Y, hence X⟂YX \!\!\! Y. ∎ E.2 Proof of Lemma 2 (Conditional HSIC Consistency) Proof. Equivalence. Given Z=zZ=z, X and Y are εX _X and εY _Y shifted by the constants mX(z)m_X(z) and mY(z)m_Y(z), so X⟂Y|ZX \!\!\! Y Z if and only if εX⟂εY|Z _X \!\!\! _Y Z. Since (εX,εY)⟂Z( _X, _Y) \!\!\! Z, the conditional law of (εX,εY)( _X, _Y) given Z equals its marginal law, so this is equivalent to εX⟂εY _X \!\!\! _Y. Consistency. The population residuals are εX _X and εY _Y. For bounded L-Lipschitz kernels, replacing εX,i _X,i by x~i=xi−m^X(zi) x_i=x_i- m_X(z_i) changes each Gram entry by at most L(‖x~i−εX,i‖+‖x~j−εX,j‖)L(\| x_i- _X,i\|+\| x_j- _X,j\|). Since centering satisfies n−2∑i,j|(HAH)ij|≤4n−2∑i,j|Aij|n^-2 _i,j|(HAH)_ij|≤ 4n^-2 _i,j|A_ij| and the kernels are bounded, |HSIC^(X~,Y~)−HSIC^(εX,εY)|≤Cn∑i=1n(‖m^X(zi)−mX(zi)‖+‖m^Y(zi)−mY(zi)‖) | HSIC( X, Y)- HSIC( _X, _Y) |≤ Cn _i=1^n (\| m_X(z_i)-m_X(z_i)\|+\| m_Y(z_i)-m_Y(z_i)\| ) for a constant C depending on L and the kernel bounds. By Cauchy–Schwarz and L2L^2-consistency, the right-hand side converges to 00 in probability. By Lemma 4, HSIC^(εX,εY)→HSIC(εX,εY) HSIC( _X, _Y) pHSIC( _X, _Y). Hence HSIC^(X~,Y~)→HSIC(εX,εY) HSIC( X, Y) pHSIC( _X, _Y), which by Lemma 1 and the equivalence above is 00 if X⟂Y|ZX \!\!\! Y Z and positive otherwise. ∎ E.3 Proof of Lemma 3 (SoftMin Properties) Proof. We establish each property. Part 1: Bounds. For the lower bound, note that for all k: ∑j=13exp(−βVj)≥exp(−βVk) _j=1^3 (-β V_j)≥ (-β V_k) (51) Taking logarithms and the infimum over k: −1βlog∑j=13exp(−βVj)≤−1βlogexp(−βminkVk)=minkVk- 1β _j=1^3 (-β V_j)≤- 1β (-β _kV_k)= _kV_k (52) For the upper bound: ∑j=13exp(−βVj)≤3exp(−βminkVk) _j=1^3 (-β V_j)≤ 3 (-β _kV_k) (53) Therefore: −1βlog∑j=13exp(−βVj)≥−1βlog(3exp(−βminkVk))=minkVk−log3β- 1β _j=1^3 (-β V_j)≥- 1β (3 (-β _kV_k) )= _kV_k- 3β (54) Combining: minkVk−log3β≤SoftMinβ≤minkVk _kV_k- 3β _β≤ _kV_k. Part 2: Limiting behavior. Let k∗=argminkVkk^*= _kV_k and assume uniqueness (the general case follows by continuity). Then: ∑j=13exp(−βVj)=exp(−βVk∗)(1+∑j≠k∗exp(−β(Vj−Vk∗))) _j=1^3 (-β V_j)= (-β V_k^*) (1+ _j≠ k^* (-β(V_j-V_k^*)) ) (55) Since Vj−Vk∗>0V_j-V_k^*>0 for j≠k∗j≠ k^*: ∑j≠k∗exp(−β(Vj−Vk∗))→0as β→∞ _j≠ k^* (-β(V_j-V_k^*))→ 0 β→∞ (56) Therefore: SoftMinβ=−1βlogexp(−βVk∗)−1βlog(1+o(1))=Vk∗+o(1)→minkVkSoftMin_β=- 1β (-β V_k^*)- 1β (1+o(1) )=V_k^*+o(1)→ _kV_k (57) Part 3: Gradient. By direct computation: ∇ϕSoftMinβ _φSoftMin_β =−1β⋅∑k=13exp(−βVk)⋅(−β∇ϕVk)∑j=13exp(−βVj) =- 1β· _k=1^3 (-β V_k)·(-β _φV_k) _j=1^3 (-β V_j) (58) =∑k=13exp(−βVk)∑j=13exp(−βVj)∇ϕVk = _k=1^3 (-β V_k) _j=1^3 (-β V_j) _φV_k (59) =∑k=13αk(ϕ)∇ϕVk = _k=1^3 _k(φ) _φV_k (60) where αk(ϕ)=exp(−βVk)/∑jexp(−βVj) _k(φ)= (-β V_k)/ _j (-β V_j) are the adaptive weights. ∎ E.4 Proof of Theorem D.1 (Violation Characterization) Proof. By the Markov property, d-separation in k∗G_k^* implies conditional independence; by faithfulness, d-connection implies dependence. In all three graphs, E enters only through E→UE→ U. 1G_1. The only parent of Xz⟂X_z is Y. Every path from E to Xz⟂X_z either passes through Y as a non-collider or contains an unconditioned collider (Z or XzX_z); hence Xz⟂E|YX_z \!\!\! E Y (ℐ1I_1 holds). The path E→U→Y→Xz⟂E→ U→ Y→ X_z is open without conditioning, so Xz⟂̸⟂EX_z \!\!\! E (ℐ3I_3 fails); E→U→YE→ U→ Y is open given Xz⟂X_z , so Y⟂̸⟂E|Xz⟂Y \!\!\! E X_z (ℐ2I_2 fails). 2G_2. Given Xz⟂X_z , the path E→U→Xz⟂→YE→ U→ X_z → Y is blocked at Xz⟂X_z , and every other path from E to Y contains an unconditioned collider (XzX_z or Z); hence Y⟂E|Xz⟂Y \!\!\! E X_z (ℐ2I_2 holds). The path E→U→Xz⟂E→ U→ X_z is open both unconditionally and given Y, so ℐ3I_3 and ℐ1I_1 fail. 3G_3. Xz⟂X_z has no parents, and every path from E to Xz⟂X_z contains a collider (Y, Z, or XzX_z); hence Xz⟂EX_z \!\!\! E (ℐ3I_3 holds). Conditioning on the collider Y opens E→U→Y←Xz⟂E→ U→ Y← X_z , so ℐ1I_1 fails; E→U→YE→ U→ Y is open given Xz⟂X_z , so ℐ2I_2 fails. Lemma 1 (for V3V_3) and Lemma 2 (for V1,V2V_1,V_2) convert independence into zero violation and dependence into positive violation. Achievability follows because Xz⟂∈ΦϵX_z ∈ _ε attains Vk∗(Xz⟂)=0V_k^*^*(X_z )=0 and Vk∗≥0V_k^*≥ 0. Non-achievability follows from Definition D.2, since ϵ≤ϵ0ε≤ _0 implies Φϵ⊆Φ0 _ε _0. ∎ E.5 Proof of Theorem D.2 See D.2 Proof. Let ϕn∈Φϵ _n∈ _ε with ϵ≤ϵ0ε≤ _0 attain the infimum of Theorem D.1, as the global minimizer does (Theorem D.4), so that Vk∗(ϕn)=0V_k^*^*( _n)=0; non-achievability gives Vk∗(ϕn)≥η/2V_k^*( _n)≥η/2 for k≠k∗k≠ k^*. Apply Lemma 4 with ϵ=η/8ε=η/8 to each k and take a union bound over k∈1,2,3k∈\1,2,3\: with probability at least 1−3⋅2exp(−nη2/(64⋅32κ4))=1−6exp(−nη2/(2048κ4))1-3· 2 (-nη^2/(64· 32κ^4))=1-6 (-nη^2/(2048κ^4)), we have |Vk(ϕn)−Vk∗(ϕn)|≤η/8|V_k( _n)-V_k^*( _n)|≤η/8 for all k. On this event, Vk∗(ϕn)≤η/8V_k^*( _n)≤η/8 while Vk(ϕn)≥η/2−η/8=3η/8V_k( _n)≥η/2-η/8=3η/8 for k≠k∗k≠ k^*, so k^n=k∗ k_n=k^*. Since ∑nexp(−cn)<∞ _n (-cn)<∞ for c>0c>0, the Borel–Cantelli lemma gives k^n→a.s.k∗ k_n a.s.k^*. ∎ E.6 Proof of Corollary 1 Proof. Split the tolerance η/8η/8 of the proof of Theorem D.2 into two halves. Sampling error: Lemma 4 with ϵ=η/16ε=η/16 and a union bound over k give |Vk(ϕn)−Vk∗(ϕn)|≤η/16|V_k( _n)-V_k^*( _n)|≤η/16 for all k with probability at least 1−6exp(−nη2/(256⋅32κ4))=1−6exp(−nη2/(8192κ4))1-6 (-nη^2/(256· 32κ^4))=1-6 (-nη^2/(8192κ^4)). Approximation error: with probability at least 1−δ′1-δ , |VkRFF(ϕn)−Vk(ϕn)|≤Cκ2εD≤η/16|V_k^RFF( _n)-V_k( _n)|≤ Cκ^2 _D≤η/16. On both events, Vk∗RFF(ϕn)≤η/8V_k^*^RFF( _n)≤η/8 and VkRFF(ϕn)≥η/2−η/8=3η/8V_k^RFF( _n)≥η/2-η/8=3η/8 for k≠k∗k≠ k^*, so k^nRFF=k∗ k_n^RFF=k^*. By the uniform convergence of random Fourier features [50], D=Ω(εD−2dlog(n/δ′))D= ( _D^-2d (n/δ )) features achieve error εD _D with probability 1−δ′1-δ ; setting εD=η/(16Cκ2) _D=η/(16Cκ^2) gives D=Ω(κ4η−2dlog(n/δ′))D= (κ^4η^-2d (n/δ )). ∎ E.7 Proof of Theorem D.3 See D.3 Proof. Let (ϕn,wn)( _n,w_n) minimize ℒ(ϕ,w)=maxeℛ^e(w∘ϕ)+λ⋅SoftMinβ(Vk(ϕ))L(φ,w)= _e R_e(w φ)+λ·SoftMin_β(\V_k(φ)\). Step 1: First-order optimality. At a stationary point, first-order optimality gives:: λ∑k=13αk(ϕn)∇ϕVk(ϕn)=−∇ϕmaxeℛ^e(wn∘ϕn)λ _k=1^3 _k( _n) _φV_k( _n)=- _φ _e R_e(w_n _n) (61) Step 2: Bounding the violation gradient. Under Assumption 6, the violation functions are LVL_V-smooth: ‖∇ϕVk(ϕ)‖≤LV\| _φV_k(φ)\|≤ L_V (62) From the optimality condition: ‖∑k=13αk(ϕn)∇ϕVk(ϕn)‖≤1λ‖∇ϕmaxeℛ^e‖ \| _k=1^3 _k( _n) _φV_k( _n) \|≤ 1λ\| _φ _e R_e\| (63) Step 3: Relating violation to gradient. For the true structure k∗k^*, by Corollary 2, αk∗(ϕn)→1 _k^*( _n)→ 1 as n→∞n→∞. The SoftMin gradient is approximately: ∇ϕSoftMinβ≈∇ϕVk∗(ϕn) _φSoftMin_β≈ _φV_k^*( _n) (64) By the LVL_V-smoothness and the connection between gradient norms and function values for convex-like functions: Vk∗(ϕn)≤Vk∗(ϕ∗)+⟨∇Vk∗(ϕ∗),ϕn−ϕ∗⟩+LV2‖ϕn−ϕ∗‖2V_k^*( _n)≤ V_k^*(φ^*)+ ∇ V_k^*(φ^*), _n-φ^* + L_V2\| _n-φ^*\|^2 (65) Since Vk∗(ϕ∗)=0V_k^*(φ^*)=0 and ∇Vk∗(ϕ∗)=0∇ V_k^*(φ^*)=0 at the optimal invariant representation: Vk∗(ϕn)≤LV2‖ϕn−ϕ∗‖2V_k^*( _n)≤ L_V2\| _n-φ^*\|^2 (66) Step 4: Bounding representation error. The optimality of (ϕn,wn)( _n,w_n) for ℒL implies: ℒ(ϕn,wn)≤ℒ(ϕ∗,w∗)L( _n,w_n) (φ^*,w^*) (67) This gives: maxeℛ^e(wn∘ϕn)+λ⋅SoftMinβ(Vk(ϕn))≤maxeℛ^e(w∗∘ϕ∗)+λ⋅0 _e R_e(w_n _n)+λ·SoftMin_β(\V_k( _n)\)≤ _e R_e(w^* φ^*)+λ· 0 (68) Therefore: λ⋅SoftMinβ(Vk(ϕn))≤maxeℛ^e(w∗∘ϕ∗)λ·SoftMin_β(\V_k( _n)\)≤ _e R_e(w^* φ^*) (69) By Lemma 3 Part 1: λ(Vk∗(ϕn)−log3β)≤λ⋅SoftMinβ≤maxeℛ^e(w∗∘ϕ∗)λ (V_k^*( _n)- 3β )≤λ·SoftMin_β≤ _e R_e(w^* φ^*) (70) Rearranging: Vk∗(ϕn)≤1λmaxeℛ^e(w∗∘ϕ∗)+log3βV_k^*( _n)≤ 1λ _e R_e(w^* φ^*)+ 3β (71) Step 5: Adding estimation error. The empirical violation Vk∗(ϕn)V_k^*( _n) differs from the population violation Vk∗(ϕn)V_k^*^*( _n) by O(n−1/2)O(n^-1/2) (Lemma 4). Therefore: Vk∗(ϕn)≤C1λmaxeℛ^e(w∗∘ϕ∗)+C2n−1/2V_k^*( _n)≤ C_1λ _e R_e(w^* φ^*)+C_2n^-1/2 (72) where C1=1C_1=1 and C2C_2 absorbs the log3/β 3/β term and concentration error. ∎ E.8 Proof of Theorem D.4 (Global Optimality) Proof. Assume the representation class Φ contains ϕ∗φ^* satisfying ℐk∗I_k^*, and W contains the Bayes-optimal predictor w∗w^*. Part 1: Exact invariance. By Theorem D.1 Part 1, Vk∗(ϕ∗)=0V_k^*^*(φ^*)=0. Since ϕ∗∈Φφ^*∈ : infϕ∈ΦVk∗(ϕ)≤Vk∗(ϕ∗)=0 _φ∈ V_k^*^*(φ)≤ V_k^*^*(φ^*)=0 (73) Since Vk∗≥0V_k^*^*≥ 0 by definition (HSIC is non-negative), equality holds: infϕ∈ΦVk∗(ϕ)=0 _φ∈ V_k^*^*(φ)=0 (74) At the global minimum (ϕ†,w†)(φ ,w ) of ℒL, if Vk∗(ϕ†)>0V_k^*(φ )>0, then replacing ϕ†φ with ϕ∗φ^* would strictly decrease the SoftMin term while achieving no worse prediction loss (since ϕ∗φ^* is sufficient). This contradicts optimality of ϕ†φ . Therefore Vk∗(ϕ†)=0V_k^*(φ )=0, which implies ℐk∗(ϕ†)I_k^*(φ ) holds. Part 2: Bayes risk. At the global minimum, the prediction loss component satisfies: maxeℛe(w†∘ϕ†)≤maxeℛe(w∗∘ϕ∗) _eR_e(w φ )≤ _eR_e(w^* φ^*) (75) Since ϕ∗φ^* achieves the invariance ℐk∗I_k^* and w∗w^* is Bayes-optimal: maxeℛe(w∗∘ϕ∗)=ℛ∗ _eR_e(w^* φ^*)=R^* (76) where ℛ∗R^* is the Bayes risk (the irreducible error). For the reverse inequality: since ϕ†φ satisfies ℐk∗I_k^* (from Part 1) and invariant representations cannot improve upon the Bayes risk: maxeℛe(w†∘ϕ†)≥ℛ∗ _eR_e(w φ ) ^* (77) Combining: maxeℛe(w†∘ϕ†)=ℛ∗ _eR_e(w φ )=R^*. ∎ E.9 Proof of Theorem D.5 (Generalization Bound) Proof. Let (ϕn,wn)( _n,w_n) minimize ℒL over training environments ℰtrainE_train. We bound the risk on a test environment eteste_test. Step 1: Decomposition. ℛetest(wn∘ϕn) _e_test(w_n _n) =ℛetest(wn∘ϕn)−ℛe∗(wn∘ϕn)⏟(I) distribution shift = R_e_test(w_n _n)-R_e^*(w_n _n)_(I) distribution shift (78) +ℛe∗(wn∘ϕn)−ℛ^e∗(wn∘ϕn)⏟(I) generalization gap + R_e^*(w_n _n)- R_e^*(w_n _n)_(I) generalization gap (79) +ℛ^e∗(wn∘ϕn)⏟(I) training error + R_e^*(w_n _n)_(I) training error (80) where e∗=argmine∈ℰtraind(etest,e)e^*= _e _traind_K(e_test,e) is the closest training environment. Step 2: Bounding (I) - Distribution shift. By the definition of interventional kernel distance (Definition D.3) and Lipschitz continuity of wnw_n (Assumption 5): The change in risk across environments depends on how the representation changes. Under the invariance condition ℐk∗I_k^*: |ℛetest(wn∘ϕn)−ℛe∗(wn∘ϕn)|≤Lw⋅[‖ϕn(Xetest)−ϕn(Xe∗)‖]|R_e_test(w_n _n)-R_e^*(w_n _n)|≤ L_w·E[\| _n(X_e_test)- _n(X_e^*)\|] (81) The representation difference is bounded by the kernel distance: [‖ϕn(Xetest)−ϕn(Xe∗)‖]≤Cϕ⋅d(etest,e∗)E[\| _n(X_e_test)- _n(X_e^*)\|]≤ C_φ· d_K(e_test,e^*) (82) where CϕC_φ depends on the representation class complexity. Therefore: (I)≤Lw⋅d(etest,ℰtrain)(I)≤ L_w· d_K(e_test,E_train) (83) Step 3: Bounding (I) - Generalization gap. By standard Rademacher complexity bounds, for the function class (x,y)↦ℓ(w(ϕ(x)),y):ϕ∈Φ,w∈\(x,y) (w(φ(x)),y):φ∈ ,w \: ℛe∗(wn∘ϕn)−ℛ^e∗(wn∘ϕn)≤2ℜn(ℓ∘Φ)+log(2/δ)2nR_e^*(w_n _n)- R_e^*(w_n _n)≤ 2 R_n( )+ (2/δ)2n (84) with probability ≥1−δ/2≥ 1-δ/2. The Rademacher complexity satisfies: ℜn(ℓ∘Φ)≤LℓLwℜn(Φ)=O((Φ,)n) R_n( )≤ L_ L_w R_n( )=O ( C( ,W)n ) (85) where (Φ,)C( ,W) is the combined complexity (e.g., VC dimension, covering number). Step 4: Bounding (I) - Training error. By definition: ℛ^e∗(wn∘ϕn)≤maxe∈ℰtrainℛ^e(wn∘ϕn) R_e^*(w_n _n)≤ _e _train R_e(w_n _n) (86) Step 5: Invariance penalty. The violation Vk∗(ϕn)V_k^*( _n) contributes to generalization error when the invariance is not exactly satisfied. Under approximate invariance: ℛetest(wn∘ϕn)−ℛe∗(wn∘ϕn)≤Lw⋅d(etest,e∗)+Cη⋅Vk∗(ϕn)R_e_test(w_n _n)-R_e^*(w_n _n)≤ L_w· d_K(e_test,e^*)+C_η· V_k^*( _n) (87) where CηC_η depends on the faithfulness gap (larger η means violations more directly impact generalization). Step 6: Combining bounds. Adding the three terms: ℛetest(wn∘ϕn) _e_test(w_n _n) ≤maxe∈ℰtrainℛ^e(wn∘ϕn)+Lw⋅d(etest,ℰtrain) ≤ _e _train R_e(w_n _n)+L_w· d_K(e_test,E_train) (88) +O((Φ,)+log(1/δ)n)+Cη⋅Vk∗(ϕn) +O ( C( ,W)+ (1/δ)n )+C_η· V_k^*( _n) (89) Defining ξn _n as the estimation error term: ξn=O((Φ,)+log(1/δ)n)+Cη⋅Vk∗(ϕn) _n=O ( C( ,W)+ (1/δ)n )+C_η· V_k^*( _n) (90) completes the proof. ∎ E.10 Proof of Corollary 3 (Asymptotic Optimality) Proof. We show convergence in probability to the oracle risk. Step 1: Vanishing estimation error. From Theorem D.5, the estimation error satisfies: ξn=O((Φ,)n)+Cη⋅Vk∗(ϕn)→0 _n=O ( C( ,W)n )+C_η· V_k^*( _n) p0 (91) as n→∞n→∞, since: • The first term is O(n−1/2)→0O(n^-1/2)→ 0 • Vk∗(ϕn)→0V_k^*( _n)→ 0 by Theorem D.3 Step 2: Training error convergence. By uniform convergence of empirical risk to population risk: maxe∈ℰtrainℛ^e(wn∘ϕn)→maxe∈ℰtrainℛe(w∗∘ϕ∗) _e _train R_e(w_n _n) p _e _trainR_e(w^* φ^*) (92) where (ϕ∗,w∗)(φ^*,w^*) is the optimal invariant solution. Step 3: Invariance and optimality. As n→∞n→∞: • ϕn→ϕ∗∈Φk∗ _n→φ^*∈ _k^* (the class of representations satisfying ℐk∗I_k^*) • wn→w∗(ϕ∗)=argminwℛ(w∘ϕ∗)w_n→ w^*(φ^*)= _wR(w φ^*) Step 4: Limit. Taking n→∞n→∞ in Theorem D.5: limn→∞ℛetest(wn∘ϕn) _n→∞R_e_test(w_n _n) =limn→∞[maxeℛ^e(wn∘ϕn)+Lw⋅d(etest,ℰtrain)+ξn] = _n→∞ [ _e R_e(w_n _n)+L_w· d_K(e_test,E_train)+ _n ] (93) =maxeℛe(w∗∘ϕ∗)+Lw⋅d(etest,ℰtrain) = _eR_e(w^* φ^*)+L_w· d_K(e_test,E_train) (94) For invariant representations, the shift penalty captures the irreducible distribution mismatch. The oracle risk over Φk∗ _k^* is: infϕ∈Φk∗ℛetest(w∗(ϕ)∘ϕ)=maxeℛe(w∗∘ϕ∗)+Lw⋅d(etest,ℰtrain) _φ∈ _k^*R_e_test(w^*(φ) φ)= _eR_e(w^* φ^*)+L_w· d_K(e_test,E_train) (95) Therefore: ℛetest(wn∘ϕn)→infϕ∈Φk∗ℛetest(w∗(ϕ)∘ϕ)R_e_test(w_n _n) p _φ∈ _k^*R_e_test(w^*(φ) φ) (96) ∎ E.11 Proof of Theorem D.6 (Sample Complexity) Proof. We derive the sample complexity for joint structure identification and generalization. Requirement 1: Structure identification. By Theorem D.2, ℙ(k^n≠k∗)≤6exp(−nη2/(2048κ4))≤δ/2P( k_n≠ k^*)≤ 6 (-nη^2/(2048κ^4))≤δ/2 whenever n≥2048κ4η2log12δ.n≥ 2048κ^4η^2 12δ. (97) For this to be ≤δ/2≤δ/2: 3exp(−nδmin232κ4)≤δ23 (- n _ ^232κ^4 )≤ δ2 (98) Solving for n: n≥32κ4δmin2log6δn≥ 32κ^4 _ ^2 6δ (99) Requirement 2: Generalization error. From Theorem D.5, the estimation error is: ξn=O((Φ,)+log(2/δ)n)+Cη⋅Vk∗(ϕn) _n=O ( C( ,W)+ (2/δ)n )+C_η· V_k^*( _n) (100) For this to be ≤ϵ≤ε with probability ≥1−δ/2≥ 1-δ/2, we need: (Φ,)+log(2/δ)n≤ϵ2 C( ,W)+ (2/δ)n≤ ε2 (101) and Cη⋅Vk∗(ϕn)≤ϵ2C_η· V_k^*( _n)≤ ε2 (102) The first condition requires: n≥4((Φ,)+log(2/δ))ϵ2n≥ 4(C( ,W)+ (2/δ))ε^2 (103) For the second condition, using the violation bound from Theorem D.3: Vk∗(ϕn)≤C1λℛ∗+C2n−1/2V_k^*( _n)≤ C_1λR^*+C_2n^-1/2 (104) For fixed λ and ℛ∗R^*, requiring Vk∗(ϕn)≤ϵ/(2Cη)V_k^*( _n)≤ε/(2C_η) gives another O(n−1/2)O(n^-1/2) requirement, absorbed into the first bound. Combining requirements. Taking the maximum: n=O(κ4η2log12δ+(Φ,)+log(2/δ)ϵ2)n=O ( κ^4η^2 12δ+ C( ,W)+ (2/δ)ε^2 ) (105) Simplifying n=O(κ4η2log12δ+κ4ϵ2log|Φ|δ)n=O ( κ^4η^2 12δ+ κ^4ε^2 | |δ ) (106) where |Φ|| | denotes the covering number (related to (Φ,)C( ,W) via standard bounds). Union bound. The total failure probability is: ℙ(k^n≠k∗ or ξn>ϵ)≤δ2+δ2=δP( k_n≠ k^* or _n>ε)≤ δ2+ δ2=δ (107) Therefore, with probability ≥1−δ≥ 1-δ, both structure identification succeeds and generalization error is ≤ϵ≤ε. ∎ E.12 Proof of Proposition 1 (Differentiable CVaR) Proof. Recall the CVaR definition: CVaRρ(re)=inft∈ℝt+1ρ|ℰ|∑e∈ℰmax(re−t,0)CVaR_ρ(\r_e\)= _t \t+ 1ρ|E| _e (r_e-t,0) \ (108) Convexity. The function g(t)=t+1ρ|ℰ|∑emax(re−t,0)g(t)=t+ 1ρ|E| _e (r_e-t,0) is convex in t (sum of affine and piecewise-linear convex functions). The infimum of convex functions preserves convexity in the parameters re\r_e\. Optimal threshold. Taking the derivative w.r.t. t: ∂g∂t=1−1ρ|ℰ|∑e:re>t1=1−|e:re>t|ρ|ℰ| ∂ g∂ t=1- 1ρ|E| _e:r_e>t1=1- |\e:r_e>t\|ρ|E| (109) Setting to zero: |e:re>t∗|=ρ|ℰ||\e:r_e>t^*\|=ρ|E|. Thus t∗t^* is the (1−ρ)(1-ρ)-quantile of re\r_e\. Subgradient. By Danskin’s theorem, the subgradient of CVaRρCVaR_ρ w.r.t. rer_e at the optimal t∗t^* is: ∂CVaRρ∂re=1ρ|ℰ|[re≥t∗] _ρ∂ r_e= 1ρ|E|I[r_e≥ t^*] (110) By the chain rule, the subgradient w.r.t. ϕφ is: ∂ϕCVaRρ=∑e∈ℰ∂CVaRρ∂re∇ϕre=∑e∈ℰpe∗∇ϕℛ^e _φCVaR_ρ= _e _ρ∂ r_e _φr_e= _e p_e^* _φ R_e (111) where pe∗=1ρ|ℰ|[re≥t∗]p_e^*= 1ρ|E|I[r_e≥ t^*]. Note that ∑epe∗=1 _ep_e^*=1 (the weights form a probability distribution concentrated on the worst ρ-fraction of environments). ∎ E.13 Proof of Theorem F.1 (Algorithm Convergence) Proof. We analyze convergence of Algorithm 1 using a Lyapunov argument. Step 1: Lyapunov function. Define: Ψt:=ℒ(ϕt,wt)+γDKL((t)∥k∗) _t:=L( _t,w_t)+γ D_KL( α^(t)\|e_k^*) (112) where k∗e_k^* is the indicator vector for the true structure and γ>0γ>0 is a coupling constant. The KL divergence term penalizes deviation from the correct structure: DKL((t)∥k∗)=−logαk∗(t)≥0D_KL( α^(t)\|e_k^*)=- _k^*^(t)≥ 0 (113) with equality iff αk∗(t)=1 _k^*^(t)=1. Step 2: Per-step descent. Under smoothness (Assumption 6), for step size ηt _t: ℒ(ϕt+1,wt+1) ( _t+1,w_t+1) ≤ℒ(ϕt,wt)−ηt‖∇ℒ(ϕt,wt)‖2+Lηt22‖∇ℒ‖2+noise ( _t,w_t)- _t\| ( _t,w_t)\|^2+ L _t^22\| \|^2+noise (114) where L is the smoothness constant. Taking expectations over minibatch randomness: [ℒ(ϕt+1,wt+1)∣ℱt]≤ℒ(ϕt,wt)−ηt(1−Lηt/2)‖∇ℒ‖2+ηt2σ2E[L( _t+1,w_t+1) _t] ( _t,w_t)- _t(1-L _t/2)\| \|^2+ _t^2σ^2 (115) where σ2σ^2 is the gradient variance bound. Step 3: Weight evolution. The adaptive weights evolve as: αk(t+1)=exp(−βt+1Vk(ϕt+1))∑jexp(−βt+1Vj(ϕt+1)) _k^(t+1)= (- _t+1V_k( _t+1)) _j (- _t+1V_j( _t+1)) (116) As Vk∗(ϕt)→0V_k^*( _t)→ 0 (by the invariance optimization) and Vk(ϕt)≥η/2V_k( _t)≥η/2 for k≠k∗k≠ k^* (by Theorem D.1): αk∗(t)→1as t→∞ _k^*^(t)→ 1 t→∞ (117) The temperature annealing βt→βmax _t→ _ accelerates this concentration. Step 4: Lyapunov descent. Combining the objective descent and weight concentration: [Ψt+1−Ψt∣ℱt] [ _t+1- _t _t] ≤−cηt‖∇ℒ‖2+O(ηt2)−γ⋅(KL decrease) ≤-c _t\| \|^2+O( _t^2)-γ·(KL decrease) (118) The KL term decreases as weights concentrate, providing additional descent. Step 5: Summing over iterations. For ηt=c0/t _t=c_0/ t: ∑t=1Tηt(1−Lηt/2)=Ω(T) _t=1^T _t(1-L _t/2)= ( T) (119) ∑t=1Tηt2=O(logT) _t=1^T _t^2=O( T) (120) Summing the descent inequality: ∑t=1Tηt[‖∇ℒ(ϕt,wt)‖2]≤Ψ1−[ΨT+1]+O(logT)≤Ψ1+O(logT) _t=1^T _tE[\| ( _t,w_t)\|^2]≤ _1-E[ _T+1]+O( T)≤ _1+O( T) (121) Therefore: mint≤T[∥∇ℒ(ϕt,wt)∥2]≤Ψ1+O(logT)∑t=1Tηt=O(T−1/2) _t≤ TE[\| ( _t,w_t)\|^2]≤ _1+O( T) _t=1^T _t=O(T^-1/2) (122) This establishes Part 1 (stationarity). Parts 2 and 3. Part 2 (weight concentration) follows from Corollary 2 applied at iteration T. Part 3 (violation decrease) follows from the objective structure: at a stationary point, the SoftMin term (dominated by Vk∗V_k^* for large β) must be small for the total gradient to vanish. Specifically, from the optimality condition: λ∇ϕVk∗(ϕT)≈−∇ϕmaxeℛ^eλ _φV_k^*( _T)≈- _φ _e R_e (123) By smoothness of Vk∗V_k^*: Vk∗(ϕT)≤12LV∥∇Vk∗(ϕT)∥2=O(λ−2)⋅O(T−1/2)=O(T−1/2)V_k^*( _T)≤ 12L_V\|∇ V_k^*( _T)\|^2=O(λ^-2)· O(T^-1/2)=O(T^-1/2) (124) for fixed λ. ∎ E.14 Proof of Theorem G.1 (Heterogeneous Structure Identification) Proof. We extend the homogeneous case to per-feature structure identification. Step 1: Per-feature violation separation. For each feature j with true structure σ∗(j)=kj∗≠⊥σ^*(j)=k_j^*≠ , by Theorem D.1 applied to the j-th component: Vkj∗(j)(ϕ∗) V_k_j^*^(j)(φ^*) =0 =0 (125) Vk(j)(ϕ∗) V_k^(j)(φ^*) ≥η(j)/2for k≠kj∗ ≥η^(j)/2 k≠ k_j^* (126) where η(j)η^(j) is the per-feature faithfulness gap. Step 2: Concentration for each feature. By Lemma 4 applied to feature j: ℙ(|Vk(j)(ϕn)−Vk(j)∗(ϕn)|>η(j)8)≤2exp(−n(η(j))22048κ4)P (|V_k^(j)( _n)-V_k^(j)*( _n)|> η^(j)8 )≤ 2 (- n(η^(j))^22048κ^4 ) (127) Step 3: Threshold calibration. Under H0:σ∗(j)=⊥H_0:σ^*(j)= (no structure applies), all three violations are positive. The permutation null distribution of HSIC HSIC concentrates around zero, and the (1−α)(1-α)-quantile τα(j) _α^(j) provides a valid threshold. By construction of permutation tests: ℙ(minkVk(j)(ϕn)≤τα(j)∣σ∗(j)=⊥)≤αP ( _kV_k^(j)( _n)≤ _α^(j) σ^*(j)= )≤α (128) Step 4: Per-feature identification. For feature j with σ∗(j)=kj∗≠⊥σ^*(j)=k_j^*≠ : On the high-probability event (concentration holds): • Vkj∗(j)(ϕn)≤η(j)/8<τα(j)V_k_j^*^(j)( _n)≤η^(j)/8< _α^(j) for appropriate calibration • Vk(j)(ϕn)≥η(j)/2−η(j)/8=3η(j)/8>τα(j)V_k^(j)( _n)≥η^(j)/2-η^(j)/8=3η^(j)/8> _α^(j) for k≠kj∗k≠ k_j^* Therefore σ^n(j)=kj∗=σ∗(j) σ_n(j)=k_j^*=σ^*(j). Step 5: Probability bound. For each feature j with σ∗(j)≠⊥σ^*(j)≠ : ℙ(σ^n(j)=σ∗(j)) ( σ_n(j)=σ^*(j)) ≥1−ℙ(concentration fails)−ℙ(threshold miscalibrated) ≥ 1-P(concentration fails)-P(threshold miscalibrated) (129) ≥1−6exp(−n(η(j))22048κ4)−α ≥ 1-6 (- n(η^(j))^22048κ^4 )-α (130) For features with σ∗(j)=⊥σ^*(j)= : ℙ(σ^n(j)≠⊥)=ℙ(minkVk(j)≤τα(j))≤αP( σ_n(j)≠ )=P ( _kV_k^(j)≤ _α^(j) )≤α (131) by the threshold calibration. ∎ Appendix F Algorithm Details This section presents the complete optimization procedure for structure-agnostic representation learning, including the main algorithm, subroutines, complexity analysis, and convergence guarantees. F.1 Main Algorithm Algorithm 1 shows a single run; in practice we run it from M≥10M≥ 10 initializations and select the final model as described in Section 5.3. Proposition 1 (Differentiable CVaR). CVaRρCVaR_ρ is convex and subdifferentiable. The sub-gradient w.r.t. ϕφ is ∂ϕCVaRρ=∑e∈ℰpe∗∇ϕℛ^e _φCVaR_ρ= _e p_e^* _φ R_e, where pe∗=1ρ|ℰ|[ℛ^e≥t∗]p_e^*= 1ρ|E|I[ R_e≥ t^*] and t∗t^* is the optimal threshold. The structure-agnostic learning problem admits a bi-level decomposition: Outer:∗=argmin∈Δ2ℒ(ϕ∗(),w∗();) : 9.24994pt α^*= _ α∈ ^2L(φ^*( α),w^*( α); α) (132) Inner:(ϕ∗(),w∗())=argminϕ,w(CVaRρ(ℛ^e)+λ∑kαkVk(ϕ)) : 9.24994pt(φ^*( α),w^*( α))= _φ,w (CVaR_ρ(\ R_e\)+λ _k _kV_k(φ) ) (133) where Δ2=∈ℝ+3:∑kαk=1 ^2=\ α ^3_+: _k _k=1\ is the probability simplex. Algorithm 1 Structure-agnostic Causal Representation Learning (SaCRL) 0: Multi-environment data =Dee∈ℰD=\D_e\_e , regularization λ, initial temperature β0 _0, final temperature βmax _ , learning rate η, iterations T, RFF dimension D, CVaR level ρ 0: Learned representation ϕ∗φ^*, predictor w∗w^*, identified structure k k 1: // Warm-start with ERM 2: (ϕ0,w0)←argminϕ,w1|ℰ|∑e∈ℰℛ^e(w∘ϕ)( _0,w_0)← _φ,w 1|E| _e R_e(w φ) 3: // Initialize RFF maps 4: Sample ωd,bdd=1D\ _d,b_d\_d=1^D for kernels k,k,kℰk_Z,k_Y,k_E 5: for t=1,…,Tt=1,…,T do 6: // Temperature annealing 7: βt←β0+(βmax−β0)⋅min(2t/T,1) _t← _0+( _ - _0)· (2t/T,1) 8: // Compute representations across all environments 9: for each environment e∈ℰe do 10: Sample minibatch (xje,yje)j=1B\(x_j^e,y_j^e)\_j=1^B from DeD_e 11: zje←ϕt−1(xje)z_j^e← _t-1(x_j^e) for all j 12: end for 13: Pool data: (zi,yi,ei)i=1N←\(z_i,y_i,e_i)\_i=1^N← all minibatch samples with environment labels 14: // Compute structure violations (Algorithm 2) 15: (V1,V2,V3)←ComputeViolations((zi,yi,ei)i=1N)(V_1,V_2,V_3)← ComputeViolations(\(z_i,y_i,e_i)\_i=1^N) 16: // Compute adaptive weights 17: αk←exp(−βtVk)/∑j=13exp(−βtVj) _k← (- _tV_k)/ _j=1^3 (- _tV_j) for k=1,2,3k=1,2,3 18: // Compute distributionally robust prediction loss 19: re←1B∑j=1Bℓ(wt−1(zje),yje)r_e← 1B _j=1^B (w_t-1(z_j^e),y_j^e) for each e∈ℰe 20: ℒpred←CVaRρ(ree∈ℰ)L_pred _ρ(\r_e\_e ) 21: // Total loss 22: ℒ←ℒpred+λ⋅SoftMinβt(Vkk=13)L _pred+λ·SoftMin_ _t(\V_k\_k=1^3) 23: // Gradient update 24: ϕt←ϕt−1−η∇ϕℒ _t← _t-1-η _φL 25: wt←wt−1−η∇wℒw_t← w_t-1-η _wL 26: end for 27: // Final structure identification 28: Compute final violations Vk(ϕT)V_k( _T) for k=1,2,3k=1,2,3 29: k^←argmink∈1,2,3Vk(ϕT) k← _k∈\1,2,3\V_k( _T) 30: return ϕT _T, wTw_T, k k F.2 Subroutines See Algorithm 2-Algorithm 4. Definition F.1 (Random Fourier Features). For shift-invariant kernels k(x,x′)=k(x−x′)k(x,x )=k(x-x ), the random feature map is ψ:→ℝDψ:X ^D: ψ(x)=2D[cos(ω1⊤x+b1),…,cos(ωD⊤x+bD)]⊤ψ(x)= 2D [ ( _1 x+b_1),…, ( _D x+b_D) ] , where ωd∼p(ω) _d p(ω) is the spectral density and bd∼Uniform[0,2π]b_d [0,2π]. Algorithm 2 Compute Structure Violations 0: Pooled data (zi,yi,ei)i=1N\(z_i,y_i,e_i)\_i=1^N, regularization γ 0: Violations (V1,V2,V3)(V_1,V_2,V_3) 1: // Structure 1 (anti-causal): Test ϕ(X)⟂E|Yφ(X) \!\!\! E Y 2: z~←KernelResiduals(zi,yi,γ) z← KernelResiduals(\z_i\,\y_i\,γ) ⊳ Regress Z on Y 3: e~←KernelResiduals(ei,yi,γ) e← KernelResiduals(\e_i\,\y_i\,γ) ⊳ Regress E on Y 4: V1←HSIC^RFF(z~,e~)V_1← HSIC_RFF( z, e) 5: // Structure 2 (confounded-descendant): Test Y⟂E|ϕ(X)Y \!\!\! E φ(X) 6: y~←KernelResiduals(yi,zi,γ) y← KernelResiduals(\y_i\,\z_i\,γ) ⊳ Regress Y on Z 7: e~′←KernelResiduals(ei,zi,γ) e ← KernelResiduals(\e_i\,\z_i\,γ) ⊳ Regress E on Z 8: V2←HSIC^RFF(y~,e~′)V_2← HSIC_RFF( y, e ) 9: // Structure 3 (confounded-outcome): Test ϕ(X)⟂Eφ(X) \!\!\! E 10: V3←HSIC^RFF(z,e)V_3← HSIC_RFF(z,e) ⊳ No residualization needed 11: return (V1,V2,V3)(V_1,V_2,V_3) Algorithm 3 Kernel Regression Residuals 0: Target values aii=1n\a_i\_i=1^n, conditioning values cii=1n\c_i\_i=1^n, regularization γ 0: Residuals a~ii=1n\ a_i\_i=1^n 1: Compute kernel matrix: [K]ij←k(ci,cj)[K]_ij← k(c_i,c_j) 2: Solve linear system: ←(K+γIn)−1 α←(K+γ I_n)^-1a where =[a1,…,an]⊤a=[a_1,…,a_n] 3: Compute predictions: a^i←∑j=1nαjk(ci,cj) a_i← _j=1^n _jk(c_i,c_j) for i=1,…,ni=1,…,n 4: Compute residuals: a~i←ai−a^i a_i← a_i- a_i for i=1,…,ni=1,…,n 5: return a~ii=1n\ a_i\_i=1^n Algorithm 4 RFF-HSIC Computation 0: Samples uii=1n\u_i\_i=1^n, vii=1n\v_i\_i=1^n, RFF parameters (ωdu,bdu)d=1D\( _d^u,b_d^u)\_d=1^D, (ωdv,bdv)d=1D\( _d^v,b_d^v)\_d=1^D 0: HSIC^RFF(U,V) HSIC_RFF(U,V) 1: // Compute random feature embeddings 2: for i=1,…,ni=1,…,n do 3: ψiu←2/D⋅[cos(ω1u⊤ui+b1u),…,cos(ωDu⊤ui+bDu)]⊤ψ^u_i← 2/D·[ ( _1^u u_i+b_1^u),…, ( _D^u u_i+b_D^u)] 4: ψiv←2/D⋅[cos(ω1v⊤vi+b1v),…,cos(ωDv⊤vi+bDv)]⊤ψ^v_i← 2/D·[ ( _1^v v_i+b_1^v),…, ( _D^v v_i+b_D^v)] 5: end for 6: // Compute mean embeddings 7: μ^u←1n∑i=1nψiu μ_u← 1n _i=1^nψ^u_i, μ^v←1n∑i=1nψiv μ_v← 1n _i=1^nψ^v_i 8: // Compute cross-covariance 9: Σ^uv←1n∑i=1nψiu(ψiv)⊤ _uv← 1n _i=1^nψ^u_i(ψ^v_i) 10: // HSIC estimate 11: HSIC^←‖Σ^uv−μ^uμ^v⊤‖F2 HSIC←\| _uv- μ_u μ_v \|_F^2 12: return HSIC HSIC F.3 Complexity Analysis Each iteration of Algorithm 1 has time complexity O(|ℰ|⋅B⋅(dx+dz)+ND2+N3),O (|E|· B·(d_x+d_z)+ND^2+N^3 ), (134) where |ℰ||E| is the number of environments, B the batch size per environment, N=|ℰ|⋅BN=|E|· B the pooled batch size, D the RFF dimension, and dx,dzd_x,d_z the input and representation dimensions. The ND2ND^2 term is RFF-HSIC (Algorithm 4); the N3N^3 term is the kernel ridge solve (Algorithm 3), where one Cholesky factorization per conditioning variable is shared across all regression targets. Both terms depend on the minibatch size N rather than the dataset size, and the total cost is linear in the number of candidate structures. F.4 Convergence Guarantees Proposition 2 (Warm-Start Guarantee). The ERM initialization (ϕ0,w0)( _0,w_0) satisfies ϕ0∈Φϵ _0∈ _ε (Definition D.1) with ϵ=O(n−1/2)ε=O(n^-1/2). This ensures optimization begins in the predictive region, avoiding degenerate solutions where Vk=0V_k=0 for all k but prediction accuracy is trivial. Theorem F.1 (Algorithm Convergence). Under Assumptions 2–6, with step size ηt=O(t−1/2) _t=O(t^-1/2) and linear annealing βt=β0+(βmax−β0)⋅min(2t/T,1) _t= _0+( _ - _0)· (2t/T,1), Algorithm 1 satisfies: 1. Stationarity: [∥∇ℒ(ϕt,wt)∥2]≤O(T−1/2)E[\| ( _t,w_t)\|^2]≤ O(T^-1/2) 2. Weight concentration: αk∗(T)≥1−2exp(−βTη/2) _k^*^(T)≥ 1-2 (- _Tη/2) where k∗k^* is the true structure 3. Violation decrease: [Vk∗(ϕT)]≤O(T−1/2)E[V_k^*( _T)]≤ O(T^-1/2) where expectations are over stochastic minibatches. Proof sketch. Define the Lyapunov function: Ψt:=ℒ(ϕt,wt)+γDKL((t)∥∗) _t:=L( _t,w_t)+γ D_KL( α^(t)\| α^*) (135) where ∗=ek∗ α^*=e_k^* is the indicator for the true structure and γ>0γ>0 is a coupling constant. Under smoothness (Assumption 6), the expected per-step decrease satisfies: [Ψt+1−Ψt∣ℱt]≤−cηt‖∇ℒ(ϕt,wt)‖2+O(ηt2)E[ _t+1- _t _t]≤-c _t\| ( _t,w_t)\|^2+O( _t^2) (136) where ℱtF_t is the filtration up to iteration t. Summing over t=1,…,Tt=1,…,T and using ∑tηt=O(T) _t _t=O( T), ∑tηt2=O(logT) _t _t^2=O( T) yields Part 1. Part 2 follows from Theorem D.1: as Vk∗(ϕt)→0V_k^*( _t)→ 0 and Vk(ϕt)≥η/2V_k( _t)≥η/2 for k≠k∗k≠ k^*, the SoftMin weights concentrate on k∗k^*. Part 3 follows from the objective structure: reducing ℒL requires reducing SoftMinβ(Vk)SoftMin_β(\V_k\), which for large β is dominated by minkVk=Vk∗ _kV_k=V_k^*. ∎ Appendix G Extension to Heterogeneous Structures When different feature subsets follow different causal structures, we extend the framework to per-variable structure assignment. Definition G.1 (Variable-Wise Structure Assignment). Let X=(X(1),…,X(d))X=(X^(1),…,X^(d)) be d-dimensional features. A heterogeneous causal structure is a mapping σ∗:1,…,d→1,2,3,⊥σ^*:\1,…,d\→\1,2,3, \ assigning each feature to a canonical structure or ⊥ (none applies). The variable-wise violation for feature j is Vk(j)(ϕ):=Vk(ϕ(j))V_k^(j)(φ):=V_k(φ^(j)). Definition G.2 (Calibrated Threshold). For significance level α∈(0,1)α∈(0,1), the calibrated threshold for feature j is: τα(j):=quantile1−α(HSIC^null(j)) _α^(j):=quantile_1-α ( HSIC_null^(j) ) (137) where HSIC^null(j) HSIC_null^(j) is the permutation distribution under the null hypothesis of independence. Definition G.3 (Heterogeneous Objective). The structure-agnostic objective with heterogeneous assignment is: ℒhet(ϕ,w):=maxe∈ℰℛ^e(w∘ϕ)+λ∑j=1dSoftMinβ(Vk(j)(ϕ)k=13)L_het(φ,w):= _e R_e(w φ)+λ _j=1^dSoftMin_β (\V_k^(j)(φ)\_k=1^3 ) (138) Theorem G.1 (Heterogeneous Structure Identification). Under Assumptions 2–6 applied to each feature subset, let ϕn _n minimize ℒhetL_het with n samples per environment. Define the estimated per-variable structure: σ^n(j):=argmink∈1,2,3Vk(j)(ϕn)if minkVk(j)(ϕn)≤τα(j)⊥otherwise σ_n(j):= cases _k∈\1,2,3\V_k^(j)( _n)&if _kV_k^(j)( _n)≤ _α^(j)\\ &otherwise cases (139) Then for each feature j with true structure σ∗(j)≠⊥σ^*(j)≠ : ℙ(σ^n(j)=σ∗(j))≥1−6exp(−n(η(j))22048κ4)−αP ( σ_n(j)=σ^*(j) )≥ 1-6 (- n(η^(j))^22048κ^4 )-α (140) where η(j)η^(j) is the faithfulness gap of feature j. For features with σ∗(j)=⊥σ^*(j)= , the false-positive rate satisfies ℙ(σ^n(j)≠⊥)≤αP( σ_n(j)≠ )≤α. Corollary 4 (Heterogeneous Generalization). Under the conditions of Theorem G.1, the generalization bound (Theorem D.5) holds with estimation error: ξnhet=O(d⋅(Φ,)+log(d/δ)n)+Cη∑j:σ∗(j)≠⊥Vσ∗(j)(j)(ϕn) _n^het=O ( d·C( ,W)+ (d/δ)n )+C_η _j:σ^*(j)≠ V_σ^*(j)^(j)( _n) (141) The factor d reflects the cost of per-feature structure identification. Appendix H Additional Experiments H.1 Implementation H.1.1 Hyperparameter Selection The objective ℒρ(ϕ,w)L_ρ(φ,w) involves hyperparameters: λ (invariance weight), β (SoftMin temperature), γ (kernel regression regularization), D (RFF dimension), and ρ (CVaR level). Table 9 provides defaults and selection guidance. Table 9: Recommended hyperparameter settings for Algorithm 1. Parameter Default Range Selection Strategy λ 1.0 [0.5,5][0.5,5] Per experiment, by OOD accuracy (oracle selection) β0 _0 1.0 [0.5,2][0.5,2] Fixed (soft selection initially) βmax _ 50 [20,100][20,100] Higher for sharper selection γ n−1/3n^-1/3 [n−1/2,n−1/5][n^-1/2,n^-1/5] Cross-validation for kernel regression ρ 0.5 [1/|ℰ|,1][1/|E|,1] Smaller ρ when all environments are trusted D 1000 [500,2000][500,2000] Increase if approximation error high H.1.2 Numerical Stability • Add small ϵ=10−6ε=10^-6 to diagonal of kernel matrices before inversion • Clip βVkβ V_k to [−20,20][-20,20] before exp to prevent overflow • Normalize HSIC values by their empirical standard deviation for comparable scales H.1.3 Kernel Selection We use the following kernels: • Representations k_Z: Gaussian RBF with bandwidth σ=median(‖zi−zj‖)σ=median(\|z_i-z_j\|) • Labels k_Y: Delta kernel for classification; Gaussian RBF for regression • Environments kℰk_E: Delta kernel (discrete environments) H.2 Synthetic DGP Equations For each canonical structure, we generate m=3m=3 training environments with d=10d=10 features: • 1G_1 (Anti-causal): Y∼Bernoulli(0.5)Y (0.5), Z∼(0,I5)Z (0,I_5), ϵ∼(0,σe2Id)ε (0, _e^2I_d) where σe∈0.5,1.0,1.5 _e∈\0.5,1.0,1.5\ varies across environments. X=W1⋅onehot(Y)+W2⋅Z+ϵX=W_1·onehot(Y)+W_2· Z+ε where W1∈ℝd×2W_1 ^d× 2 and W2∈ℝd×5W_2 ^d× 5 are fixed random matrices. • 2G_2 (Confounded-descendant): U∼(μe,I3)U ( _e,I_3) with environment-specific means μe∈(−1,0,0),(0,0,0),(1,0,0) _e∈\(-1,0,0),(0,0,0),(1,0,0)\, X=W1⋅U+ϵX=W_1· U+ _X where ϵX∼(0,0.1Id) _X (0,0.1I_d), and Y=σ(wxy⊤X+wuy⊤U+ϵY)Y=σ(w_xy X+w_uy U+ _Y) where ϵY∼(0,0.1) _Y (0,0.1) and σ is the sigmoid function, with the outcome binarized via Bernoulli sampling. • 3G_3 (Confounded-outcome): U∼(0,I3)U (0,I_3), XZ⟂=W1⋅U+ϵ1X_Z =W_1· U+ _1 where ϵ1∼(0,σe2Id/2) _1 (0, _e^2I_d/2) with σe∈0.5,1.0,1.5 _e∈\0.5,1.0,1.5\ varying across environments, Y=σ(w⊤XZ⟂)Y=σ(w X_Z ) binarized, XZ=W2⋅onehot(Y)+W3⋅U+ϵ2X_Z=W_2·onehot(Y)+W_3· U+ _2, and X=(XZ,XZ⟂)X=(X_Z,X_Z ). All weight matrices are sampled from (0,0.52)N(0,0.5^2) and fixed across runs. H.3 Hyperparameter Sensitivity Results We sweep the invariance penalty weight λ and SoftMin temperature βmax _ across varying degrees of structural misspecification using the bidirectional DGP from Section 7.6 with ν∈0.25,0.50,0.75ν∈\0.25,0.50,0.75\. Results are shown in Table 10 to Table 11. Table 10: Sensitivity to invariance penalty weight λ (βmax=50 _ =50 fixed). OOD accuracy (%). Setting λ=0.5λ=0.5 λ=1.0λ=1.0 λ=2.0λ=2.0 λ=5.0λ=5.0 ν=0.25ν=0.25 (near-causal) 46.0 46.8 47.8 50.2 ν=0.50ν=0.50 (bidirectional) 69.0 70.3 70.4 73.9 ν=0.75ν=0.75 (near-anticausal) 99.7 99.5 99.5 99.6 Table 11: Sensitivity to SoftMin temperature βmax _ (λ=5.0λ=5.0 fixed). OOD accuracy (%). Setting βmax=20 _ =20 βmax=50 _ =50 βmax=100 _ =100 ν=0.25ν=0.25 (near-causal) 50.3 50.2 51.3 ν=0.50ν=0.50 (bidirectional) 71.4 73.9 72.9 ν=0.75ν=0.75 (near-anticausal) 99.6 99.6 99.6 λ sensitivity is mild and monotonic: total variation from λ=0.5λ=0.5 to λ=5.0λ=5.0 is at most ∼5 5p in the hardest settings, with no catastrophic failure. βmax _ is even more stable (variation ≤2.5≤ 2.5p). When the anti-causal component dominates (ν=0.75ν=0.75), all configurations reach ∼99.5 99.5%. We recommend λ∈[1.0,5.0]λ∈[1.0,5.0] and βmax=50 _ =50 as robust defaults. H.4 Crossover Between Candidate Structures Table 12 sweeps ν across the crossover between 2G_2 (ν=0ν=0) and 1G_1 (ν=1ν=1) in the bidirectional DGP of Section 7.6, comparing SoftMin with hard argminkVk _kV_k selection. The sweep uses independent seeds; at ν=0.50ν=0.50 it agrees with Table 6 within one standard error. The last column is the permutation p-value of minkVk _kV_k against the null threshold τα _α (Definition G.2). Table 12: Adaptive weights and OOD accuracy (%) across the 2G_2–1G_1 crossover. ν α1 _1 α2 _2 α3 _3 SoftMin Hard argmin Perm. p 0.40 0.407± 0.041 0.489± 0.080 0.104± 0.058 59.5± 2.0 56.5± 2.7 0.013± 0.018 0.45 0.415± 0.073 0.478± 0.139 0.107± 0.075 66.2± 1.5 62.5± 2.9 0.011± 0.012 0.50 0.421± 0.089 0.422± 0.143 0.157± 0.071 72.4± 1.5 68.6± 3.6 0.181± 0.011 0.55 0.484± 0.079 0.417± 0.143 0.099± 0.077 79.8± 2.0 76.5± 2.2 0.033± 0.039 0.60 0.549± 0.069 0.346± 0.126 0.105± 0.077 87.4± 1.9 84.7± 2.2 0.016± 0.033 H.5 Margin Calibration Table 13 bins synthetic runs with known structure by the weight margin α(1)−α(2) _(1)- _(2). Identification accuracy increases monotonically with the margin, so the margin measures how strongly the data support the selected structure. Table 13: Identification accuracy by weight margin on synthetic data. Margin Runs Correct ID accuracy (%) <0.1<0.1 48 20 41.7 0.10.1–0.30.3 60 39 65.0 0.30.3–0.50.5 72 62 86.1 >0.5>0.5 60 57 95.0 H.6 Semi-Synthetic Bayesian-Network Benchmarks Construction. For each benchmark, we select a target Y in a published network, take X as its Markov blanket, and induce environments by soft interventions on a latent parent outside the blanket. The canonical structure is read off the published DAG by the rules of Section 2. 3G_3 requires two independent root components reconverging at Y, which ALARM (12 roots) supplies twice and ASIA (2 roots) once; 1G_1 is not derivable under this recipe, since it requires the Markov blanket of Y to consist only of descendants of Y. Before training, each derived structure is verified in simulation: its own invariance holds while the competing one fails. We keep a benchmark only if the competing-invariance violation exceeds 5.0 standard deviations of the permutation null. Departures from the published networks. The SACHS consensus DAG fits its own observational data weakly, so SACHS requires larger interventions to produce contrast. For ASIA, we keep the published structure but re-parameterize two conditional probability tables whose published marginals leave the exogenous branch nearly constant. ALARM requires no adjustment. Discussion. Five benchmarks span two derivable structures across three networks (three 2G_2, two 3G_3), so a method that returns a single structure by default is correct on at most 3/5. alarm_co separates more cleanly than asia_g2 (28.3 vs. 24.6 null standard deviations) but has lower accuracy, because its target has three classes over 12 one-hot columns rather than two. alarm_g2 barely clears the 5.0 threshold and is correspondingly the least reliable. The two 3G_3 benchmarks come from different networks, so their agreement is not an artifact of a single graph. H.7 Cross-Fitting Variant Conditional HSIC via kernel ridge regression residualization can be sensitive to estimation bias. We implement a cross-fitting variant following the SplitKCI protocol [48]: fit kernel regression on fold 1, compute residuals on fold 2, swap and average. Results are displayed in Table 14. Table 14: Cross-fitting comparison on the bidirectional DGP (ν=0.5ν=0.5, n=1000n=1000 per environment). Independent run of the default configuration. Variant OOD Acc (%) No splitting (default) 73.9 ± 2.7 Cross-fitted (2-fold) 74.0 ± 2.7 ERM (reference) 61.2 ± 4.5 OOD accuracy is virtually identical (73.9 vs. 74.0), confirming that estimation-induced bias from joint residualization does not meaningfully affect downstream performance at the sample sizes used in our experiments. Both variants outperform ERM by ∼13 13p. H.8 EMA Conditioning Target Table 15 conditions V1V_1 and V2V_2 on an exponential moving average (EMA) of the encoder weights with momentum m, while gradients flow only through the online encoder. Results are nearly unchanged, so SaCRL is robust to this choice. Table 15: EMA conditioning target on synthetic 2G_2 (n=2000n=2000). m=0m=0 is the default online encoder. Momentum m Stru. ID (%) OOD Acc (%) 0 (default) 89.0 84.2 0.99 89.4 84.3 0.999 89.8 84.5 H.9 RFF Approximation Error Table 16 reports the mean absolute difference between RFF-based and exact HSIC on synthetic data (n=2000n=2000). The error decreases steadily with D, from 1.42×10−31.42× 10^-3 at D=100D=100 to 4.4×10−44.4× 10^-4 at D=2000D=2000, consistent with the η−2η^-2 scaling of the required D in Corollary 1. Table 16: RFF approximation error versus exact HSIC (synthetic, n=2000n=2000). D 100 250 500 1000 2000 Mean |HSIC^RFF−HSIC^|| HSIC_RFF- HSIC| (×10−4× 10^-4) 14.2 12.9 8.4 7.7 4.4 H.10 Effect of the Number of Training Environments Table 17 reports the effect of varying the number of training environments m on synthetic 2G_2 data (n=2000n=2000). With only two environments, structure identification is unstable due to insufficient contrast for testing the three invariances (consistent with Assumption 1 requiring rank()=3rank(D)=3). Using three or more environments yields reliable identification (≥89%≥ 89\%) and OOD accuracy (≥84%≥ 84\%), with diminishing returns beyond m=4m=4. Table 17: Effect of number of training environments m on synthetic 2G_2 data (n=2000n=2000). m Stru. ID. Acc (%) OOD Acc (%) 2 71.5 ± 4.2 74.3 ± 2.8 3 89.0 ± 1.6 84.2 ± 1.3 4 96.5 ± 1.2 85.1 ± 1.1 5 97.0 ± 0.9 85.4 ± 0.9 H.11 Computational Efficiency and Runtime Analysis H.11.1 Training Time and Memory Table 18 compares training time and memory. With D=1000D=1000 RFF features, our method adds modest overhead (∼30% 30\%) over ERM while avoiding the O(n2)O(n^2) complexity of exact HSIC. For Colored MNIST (10K images), exact kernel matrices would require >10>10 GB GPU memory and 4×4× training time; RFF introduces negligible accuracy loss while enabling scalability. Table 18: Training time and memory for 100 epochs on Colored MNIST (10K samples). Method Time (min) GPU Memory Relative ERM 8.3 2.1 GB 1.0×1.0× IRM 9.8 2.4 GB 1.18×1.18× SaCRL (RFF, D=1000D=1000) 10.9 2.7 GB 1.31×1.31× SaCRL (exact HSIC) 33.2 11.8 GB 4.0×4.0× H.11.2 Per-Epoch Runtime Breakdown Table 19 breaks down the training time into different components. We observe that computing HSIC via random Fourier features (RFF) significantly reduces the computational cost. Table 19: Runtime breakdown per epoch on PACS. Component Time (s) Fraction Forward pass 4.59 42% HSIC computation (RFF) 3.02 28% Kernel regression 1.75 16% Backward pass 1.19 11% Other 0.34 3% Total 10.9 100% H.12 Fine-Tuned Backbone Table 20: PACS (Photo held out): frozen vs. fine-tuned ResNet18. OOD accuracy (%), 3 seeds. Method Frozen Fine-tuned Δ ERM 93.2 95.1 +1.9 IRM 93.2 95.2 +2.0 VREx 93.0 94.9 +1.9 CIRCE 92.5 94.1 +1.6 SaCRL 94.0 95.0 +1.0 Appendix I Limitation Our framework targets three canonical structures arising under partial observability and confounding-mediated environment effects; cyclic systems, multi-confounder dynamics, or settings where E directly rewires X–Y fall outside this scope, though the heterogeneous extension and the empirical robustness study (Section 7.6) show graceful degradation under misspecification. In such cases, SaCRL signals the mismatch: if minkVk _kV_k exceeds the permutation threshold τα _α (Definition G.2), the data do not support a canonical structure, and the selected structure should not be interpreted as identified. Structure identification further requires sufficient environment diversity (Assumption 1) and at least m≥3m≥ 3 training environments (Appendix H.10); under near-identical environments, the invariance penalty loses informativeness even though identification remains stable. Conditional HSIC estimation depends on kernel regression quality and may benefit from cross-fitting in low-sample regimes, per results in Table 14. Appendix J Applications and Broader Connections SaCRL addresses a question that recurs across several fields: which invariance should a model enforce when the causal structure of the data is unknown? We discuss five areas where this question arises, review recent work in each, and state what SaCRL adds. Identifiability beyond an assumed structure. A central open problem in representation learning is identifiability: whether the data determine the learned representation uniquely, up to simple transformations such as permutation and element-wise reparameterization, so that it recovers the factors that generated the data rather than one of many equally predictive alternatives. Without identifiability, a representation carries no guarantee of meaning, and the guarantees built on it, including robustness under distribution shift, do not follow. The foundational results are negative: nonlinear independent component analysis is not identifiable in general [22], and unsupervised disentanglement is impossible without inductive biases on both the model and the data [36]. Progress since then has come from adding structure. Auxiliary variables such as environment labels [27], paired views from data augmentation [72], and interventions [2] each restore identifiability, and together they define the program of causal representation learning [55]. Each of these results, however, relies on a generative structure supplied by the analyst: which variables cause which, and where environments intervene. The structure itself is never identified; it enters as an assumption, and when the assumption is wrong the guarantees no longer apply. SaCRL moves identification one level up. Within the three canonical structures, it identifies from multi-environment data which structure generated the data, and therefore which invariance a representation should satisfy. Human–AI interaction and longitudinal behavioral modeling. Conversational AI systems now interact with the same people over weeks and months, which makes human–AI interaction a longitudinal, two-way process: users adapt to the system, and the system adapts to users through personalization, memory, and training on their feedback. Recent evidence shows why this process must be modeled causally. In a four-week randomized study with about 1,000 participants and more than 300,000 messages, the assigned chatbot conditions had no significant effect on psychosocial outcomes, whereas participants who chose to use the chatbot more reported worse outcomes [16], which leaves open whether heavier use drives these outcomes or reflects them. Language models tailor their answers to the beliefs users express [58], so the user’s state shapes the model’s behavior; conversely, deployed predictions and recommendations change the very behavior they later learn from [43, 23]. Preferences also differ systematically across people and cultures, as PRISM documents by linking 1,500 participants from 75 countries to 8,011 live conversations with 21 language models [28]. The gap is that any model of user or assistant behavior over time must decide what to hold invariant across users, cultures, and time periods, and that decision depends on the direction of influence. If the modeled signal is generated by the target, as when a user’s latent state generates their messages, the anti-causal invariance ℐ1I_1 applies; if the signal drives the target, as when the assistant’s responses drive the user’s next action, the causal invariances ℐ2I_2 or ℐ3I_3 apply; and if influence runs both ways, the system is in a feedback loop that no single invariance describes. SaCRL gives this decision a statistical footing. Treating users, cultural groups, or time windows as environments, it selects the invariance supported by the data, assigns different structures to different behavioral signals through its heterogeneous extension (Appendix G), and, when influence is bidirectional, keeps its structure weights spread instead of committing to a wrong structure (Section 7.6, Appendix H.4). The feedback regime thus becomes observable rather than silently misspecified. Cyclic dynamics lie outside the three canonical structures, and extending the candidate set to dynamic, feedback-coupled structures is a natural next step toward causally grounded models of long-term human–AI relationships. Multi-site clinical prediction and the direction of causation. Clinical models are trained at some hospitals and deployed at others, and dataset shift between sites is a leading cause of failure. A pneumonia detector trained on chest radiographs degraded at external hospitals, in part because it exploited site-specific signals [79]; a widely deployed proprietary sepsis model performed markedly worse under external validation than its developer reported [76]; and clinicians are now urged to treat dataset shift as a routine deployment risk [17]. A causal view explains why remedies differ. In anti-causal tasks the disease generates the image, so shifts appear as changes in prevalence or manifestation, whereas in causal tasks the inputs generate the outcome [11]. Dataset biases that look identical in the data can therefore demand different mitigation strategies [25], and shift-stable modeling requires the analyst to specify the causal graph [62]. The gap is that the causal direction must be declared by hand, per task and often per variable: symptoms are anti-causal, biomarkers are not (Figure 1). SaCRL treats each hospital as an environment and infers the direction from data. Theorem D.2 bounds the probability of selecting a wrong structure by 6exp(−nη2/(2048κ4))6 (-nη^2/(2048κ^4)), the heterogeneous extension assigns a structure to each variable (Appendix G), and the permutation test of Section 5 flags cases where no canonical structure fits, so a clinical team knows when the selected structure should not be trusted. Domain adaptation and generalization: which distribution to align. Methods for learning under distribution shift differ in which part of the joint distribution they hold fixed across domains. Marginal alignment makes P(ϕ(X))P(φ(X)) invariant, as in domain-adversarial training [18]; conditional alignment makes P(ϕ(X)∣Y)P(φ(X) Y) invariant, as under target and conditional shift [83] and generalized label shift [65]; and covariate-shift methods such as IRM make P(Y∣ϕ(X))P(Y φ(X)) invariant [4]. The appropriate assumption depends on the causal direction [54]: label shift is natural in anti-causal tasks and can be corrected with black-box predictors [33], while marginal alignment provably fails when label distributions differ across domains [85]. Recent work makes this dependence explicit: each causal graph entails different independence constraints, so any algorithm built on a single fixed constraint fails on some shifts [26, 74], consistent with the finding that no domain generalization method reliably outperforms ERM [20]. Adaptive methods such as CACM [26] still take the graph as input. SaCRL removes this requirement. Its three invariances are exactly these three alignment targets, ℐ3I_3 (marginal), ℐ1I_1 (conditional), and ℐ2I_2 (covariate shift), measured on a common HSIC scale and selected by SoftMin without the graph. Theorem D.1 shows that a predictive representation can satisfy only the invariance of the true structure, and the experiments show the cost of a wrong choice: on anti-causal Colored MNIST, IRM, which enforces ℐ2I_2, falls below ERM (21.3% vs. 27.2%), while SaCRL reaches 45.2% (Table 3). Algorithmic fairness: independence, separation, and sufficiency. The three main group-fairness criteria constrain a score or representation R with respect to a sensitive attribute A: independence (demographic parity) requires R⟂AR \!\!\! A, separation (equalized odds) requires R⟂A|YR \!\!\! A Y, and sufficiency (calibration within groups) requires Y⟂A|RY \!\!\! A R [5]. Impossibility results show that these criteria cannot hold simultaneously outside degenerate cases [29, 15]. Recent causal analyses show that the appropriate criterion depends on the data-generating graph: in anti-causal prediction, separation coincides with a notion of risk invariance [40], and whether data balancing achieves a given criterion depends on the graph [56]. The gap is that practitioners choose a criterion without knowing which one the graph permits at no cost in accuracy. With A in the role of E, SaCRL’s invariances are precisely these criteria: ℐ3I_3 is independence, ℐ1I_1 is separation, and ℐ2I_2 is sufficiency, with V3V_3, V1V_1, and V2V_2 measuring their violations. When A enters the system as E does in Figure 3, Theorem D.1 gives a causal counterpart of the impossibility results: a predictive representation can satisfy at most the criterion matching the true structure, that criterion is achievable when Xz⟂X_z is predictive, and the adaptive weights reveal which one it is. Since the predictor is a function of ϕ(X)φ(X), independence and separation of ϕ(X)φ(X) carry over to the score. The choice of a fairness criterion remains a normative decision; SaCRL informs it by showing which criterion the data can support. Diagnosing and explaining distribution shift. A related line of work asks not how to be robust, but what changed. Shift detection tests whether training and deployment distributions differ [49]. Shift attribution explains a performance drop by the mechanisms that changed, either using a known causal graph [80] or by decomposing the drop into a change in the covariate distribution (X-shift) and in the outcome given the covariates (Y|XY X-shift) [10]. A benchmark of natural tabular shifts found Y|XY X-shifts to be the most prevalent and showed that conclusions about algorithms depend on the type of shift [35]. These tools explain the failure of a fixed model after deployment; they do not decide what a model should be invariant to during training. SaCRL connects the two. Its diversity matrix D (Definition C.19) measures, across training environments, the marginal shift in P(X)P(X), the conditional shift in P(X∣Y)P(X Y), and the outcome shift in P(Y∣X)P(Y X) with kernel mean embeddings; the first and third are the X- and Y|XY X-shifts of the attribution literature, computed between training environments rather than between training and deployment. A permutation check on D before training determines which invariances are testable (Appendix C.6), and training then enforces the one the data support, so diagnosis and remedy become parts of a single objective. Invariance-based causal discovery and conditional independence testing. Heterogeneity across environments is also a source of causal knowledge. Invariant causal prediction identifies the causal parents of a target as the variables whose conditional distribution is stable across environments [44], and its nonlinear extension replaces linear models with nonparametric conditional independence tests [21]; interventional single-cell data, such as the protein-signaling measurements of Sachs et al. [52], remain standard testbeds. These methods rest on conditional independence testing, where kernel tests are widely used [82] and recent split-sample variants improve calibration [48]. A fundamental limit applies: no conditional independence test can control its level over all null distributions while retaining power against all alternatives [57], so useful guarantees require explicit assumptions. Existing methods test invariances of observed variables and return a set of parents; they neither handle high-dimensional inputs such as images nor learn a representation. SaCRL tests invariances of a learned representation ϕ(X)φ(X) and learns ϕφ jointly with the test. Its guarantees are stated under explicit assumptions consistent with this limit, namely an additive-noise model for conditional HSIC (Lemma 2) and a positive faithfulness gap η (Definition D.2). On semi-synthetic benchmarks built from the ASIA, ALARM, and SACHS networks, SaCRL selects the structure fixed by the published graph in 75–95% of runs (Table 2).