Technical Summary: Estimating Bidirectional Causal Effects with Large Scale Online Kernel Learning

Problem Statement

The paper addresses the challenge of estimating bidirectional causal effects in systems characterized by mutual dependence and heteroskedasticity. While traditional causal inference often focuses on unidirectional relationships, many real-world systems (e.g., predator-prey dynamics, policy interventions, and brain-behavior interactions) exhibit reciprocal causality. A significant methodological gap exists in applying machine learning to these scenarios, particularly when handling large-scale, high-dimensional, and streaming data.

Existing approaches face two primary limitations:

1. Instrumental Variable (IV) Methods: Conventional IV methods rely on exogenous shifts in the conditional mean of treatment variables. They often struggle with bidirectional systems unless valid instruments are available that shift one curve while leaving the other fixed.
2. Heteroskedasticity-Based Identification: While this approach (relying on exogenous variations in conditional variance) allows for the identification of simultaneous equation models (SEMs) without external instruments, traditional implementations often require dividing samples into subsamples or using complex polynomial approximations. These methods can lack flexibility in modeling nonlinear conditional variances and means, or they suffer from computational inefficiency in high-dimensional settings.

Methodology

The proposed framework integrates heteroskedasticity-based identification with large-scale online kernel learning. The core components are:

1. Model Specification

The authors define a Simultaneous Equation Model (SEM) with two endogenous variables ($y_1, y_2$) and exogenous covariates ($x$):
$$\begin{aligned} y_{1,i} &= \gamma_1 y_{2,i} + h_1(x_i, \beta_1) + \varepsilon_{1,i} \\ y_{2,i} &= \gamma_2 y_{1,i} + h_2(x_i, \beta_2) + \varepsilon_{2,i} \end{aligned}$$
Crucially, the conditional variances of the error terms are modeled as functions of the covariates:
$$g_{j,i} = V(\varepsilon_{j,i}|x_i) = \exp(f_j(x_i, \alpha_j))$$
The goal is to estimate the causal parameters $\gamma_1$ and $\gamma_2$ alongside the unknown functions $h_j$ and $f_j$.

2. Identification Strategy

The method relies on quasi-maximum likelihood (QML) estimation. Identification is achieved if the conditional variances of the two equations are not proportional to each other ($g_1 \not\propto g_2$). This allows the slopes of the structural equations to be identified through the variation in error variances, a variant of the instrumental variable approach that exploits heteroskedasticity rather than mean shifts.

3. Function Approximation via Random Fourier Features

To handle nonlinearity and high dimensionality, the unknown functions $h_j(\cdot)$ and $f_j(\cdot)$ are approximated using kernel-based functional representations.

• Instead of storing the full kernel matrix (which is computationally prohibitive for large $N$), the authors employ Random Fourier Features (RFF).
• Based on Bochner's theorem, shift-invariant kernels (e.g., Gaussian) are approximated by mapping inputs $x$ to a finite-dimensional feature space $z(x)$ using random sinusoidal projections:
  $$z(x) = (\sin(u_1^\top x), \cos(u_1^\top x), \dots, \sin(u_m^\top x), \cos(u_m^\top x))^\top$$
• This transforms the kernel method into a linear model in the feature space, enabling efficient computation.

4. Online Optimization

The parameters are estimated using an adaptive online gradient descent algorithm:

• Loss Function: Derived from the log-Gaussian quasi-likelihood.
• Optimization: The algorithm processes data in a streaming fashion (one observation at a time or in small batches).
• Adaptivity: It utilizes Adam optimization for step-size tuning and adaptive gradient clipping to ensure stability. The clipping threshold is updated as a bias-corrected exponential moving average of past gradient norms.

Key Contributions

1. Scalable Bidirectional Estimation: The paper proposes the first scalable online learning framework specifically designed to jointly estimate bidirectional causal effects in SEMs using heteroskedasticity-based identification.
2. Flexible Nonlinear Modeling: By integrating RFF with QML, the method flexibly models nonlinear conditional means and variances without the rigidity of polynomial approximations or the computational cost of full kernel matrices.
3. Computational Efficiency: The approach achieves near-linear computational scaling with respect to the number of observations and covariates, making it suitable for high-dimensional ($d=1000$) and large-sample ($N=20,000$) streaming data.
4. Theoretical Grounding: The method bridges econometric identification theory (specifically local identification conditions for SEMs under heteroskedasticity) with modern machine learning techniques.

Experimental Results

The authors evaluated the method against two baselines:

• Single-Kernel: Estimates equations independently (ignoring simultaneity).
• SEM-PAB: Uses polynomial approximations with beta function weights and Box–Cox transformations (the most flexible specification in prior work [11]).

Findings:

• Accuracy: The proposed SEM-Kernel method consistently achieved lower bias and Root Mean Squared Error (RMSE) across various Data Generating Processes (DGPs), including complex nonlinear and heteroskedastic scenarios (DGP-3).
• Bias in Baselines: The Single-Kernel method exhibited systematic bias, confirming that ignoring endogeneity leads to inconsistent estimates. The SEM-PAB method showed numerical instability and significantly higher errors.
• Scalability:
  • Accuracy vs. Features: Increasing the number of random Fourier components ($m$) improved accuracy up to $m \approx 500$, after which gains were negligible.
  • Computation Time: The SEM-Kernel method required only slightly more time than the Single-Kernel baseline (e.g., ~4.5s vs. ~3.0s for $N=20,000, d=1000$) but was orders of magnitude faster than SEM-PAB (which took >100s).
• Robustness: The method maintained stability across different sample sizes and feature dimensions without requiring recalibration.

Significance and Claims

The paper claims that the proposed framework offers a practical, scalable, and theoretically grounded solution for large-scale causal inference in natural/social sciences, policy making, and business applications where mutual dependence is present.

By combining econometric identification logic with modern machine learning, the method effectively captures complex bidirectional causal effects with near-linear computational scaling. The authors emphasize that this approach overcomes the limitations of traditional unidirectional causal inference and the computational bottlenecks of previous heteroskedasticity-based methods.

Limitations and Future Work:
The authors modestly note that the current framework estimates only linear causal effects (in terms of the structural parameters $\gamma$) and assumes symmetric error terms. They identify extending the framework to accommodate nonlinear causal relationships and strongly skewed data as important directions for future research.