TSLiNGAM: Robust Causal Discovery

• TSLiNGAM is a robust causal discovery method that extends the linear non-Gaussian acyclic model (LiNGAM) by replacing OLS with robust slope estimators.
• It employs Theil–Sen and Repeated-Median estimators to achieve high breakdown points and resistance to outliers, reliably identifying exogenous variables and causal order.
• Empirical simulations demonstrate that TSLiNGAM outperforms DirectLiNGAM, with superior accuracy and efficiency in recovering true causal structures under heavy-tailed noise.

TSLiNGAM (Theil–Sen LiNGAM) is a causal discovery algorithm that extends the DirectLiNGAM procedure to settings characterized by heavy-tailed, skewed, or contaminated noise distributions. It achieves enhanced robustness and efficiency by replacing the ordinary least squares (OLS) regression step in DirectLiNGAM with robust regression slope estimators, specifically the Theil–Sen and Repeated-Median slopes. TSLiNGAM operates within the linear, non-Gaussian, acyclic model (LiNGAM) framework and is designed to recover the causal directed acyclic graph (DAG) underlying observed multivariate data, with strong identifiability and statistical guarantees in the presence of non-Gaussian and heavy-tailed noise (Leyder et al., 2023).

1. The LiNGAM Framework

LiNGAM models a random vector $X = (X_1,\dots,X_p)^\top\in\mathbb R^p$ via a linear structural equation model: $$X_i = \sum_{j: k(j)<k(i)} b_{ij} X_j + e_i, \quad i = 1, \dots, p$$ where $k: \{1,\dots,p\} \to \{1,\dots,p\}$ encodes the true causal order (with $k(j) < k(i)$ indicating that $X_j$ can be a direct cause of $X_i$), $B = (b_{ij})$ is a connection matrix that is strictly lower triangular in causal order, and $e = (e_1, \ldots, e_p)^\top$ is a vector of mutually independent, non-Gaussian, zero-mean noise terms. In matrix notation: $$X = B X + e, \qquad B\ \text{permuted strictly lower triangular},$$ and equivalently: $X = (I - B)^{-1} e = A e,$ with $A = (I-B)^{-1}$—the ICA (Independent Component Analysis) formulation.

Identifiability is achieved under two conditions: acyclicity of the causal graph (ensuring strict lower-triangularity of $B$ in the appropriate order) and independence and non-Gaussianity of the noise components. Under these, both the connection matrix $B$ and the error vector $e$ are uniquely determined up to permutation and scaling, which are fixed by the acyclic and zero-mean conventions.

2. Shortcomings of DirectLiNGAM Under Heavy-Tailed Distributions

DirectLiNGAM recovers the causal ordering by iteratively identifying exogenous variables via independence scoring between candidates $X_j$ and the residuals $r_i^{(j)}$ from OLS regressions of the remaining variables on $X_j$: $$r_i^{(j)} = X_i - \frac{\operatorname{Cov}(X_i, X_j)}{\operatorname{Var}(X_j)} X_j = X_i - \widehat\beta_{ij}^{OLS} X_j,$$ where

$$\widehat\beta_{ij}^{OLS} = \frac{\sum_{t=1}^n (X_i^{(t)} - \bar X_i)(X_j^{(t)} - \bar X_j)}{\sum_{t=1}^n (X_j^{(t)} - \bar X_j)^2}.$$

OLS is optimal under Gaussian noise but suffers from unbounded influence and zero breakdown against outliers or heavy-tailed errors. Such sensitivity can corrupt both the estimated regression coefficients and variance estimates, leading to incorrect identification of the causal order in scenarios with heavy-tailed or contaminated disturbances.

3. The TSLiNGAM Algorithm: Robust Causal Discovery

TSLiNGAM replaces the OLS slope estimator with a robust slope functional $T(X_j, X_i)$ that satisfies regression equivariance, Fisher consistency under independence ($T(X, Y) = 0$ if $X \perp Y$), and sign preservation under argument swapping.

Key robust estimators include:

• Theil–Sen Slope

$$T_{TS}(X, Y) = \operatorname{median}_{1 \leq i < j \leq n} \frac{Y_j - Y_i}{X_j - X_i}, \qquad X_j \neq X_i.$$

This estimator offers a breakdown point of 29.3%, bounded influence, and high efficiency for heavy-tailed or skewed distributions.

• Repeated-Median Slope

$$T_{RM}(X, Y) = \operatorname{median}_{i} \left[ \operatorname{median}_{j \neq i} \frac{Y_j - Y_i}{X_j - X_i} \right].$$

This provides a breakdown point of 50%, ensuring better resistance to corruption at a slight cost to efficiency under near-Gaussian errors.

At each iteration, for each candidate $j$, robust regression slopes $T(X_j, X_i)$ are computed and residuals formed: $r_i^{(j)} = X_i - T(X_j, X_i) X_j.$ An independence score, typically based on kernel mutual information or distance correlation, is calculated between $X_j$ and each $r_i^{(j)}$. The candidate with the lowest aggregate dependence is deemed exogenous and removed, with the process repeating recursively on the residualized data. The procedure continues until all variables are ordered, at which point each variable is regressed on its ordered predecessors to recover $B$.

Algorithmic Workflow Summary

Step	Operation	Details/Methods Used
1	Robust regression	Theil–Sen / Repeated-Median slope $T(X_j,X_i)$
2	Independence scoring	Kernel mutual information or distance correlation
3	Exogenous identification	Minimum aggregate dependence
4	Residualization	Use robust slopes for new data matrix
5	Final regression	OLS or adaptive Lasso for $B$ on predecessors

4. Theoretical Guarantees and Consistency

TSLiNGAM offers identifiability under the same structural conditions as LiNGAM: data generated by a linear, non-Gaussian, acyclic SCM with independent error terms, and correlation-faithful joint distributions. For any robust slope estimator adhering to regression equivariance, Fisher consistency, and sign preservation (including Theil–Sen and Repeated-Median), TSLiNGAM recovers the true causal order and the true support of the connection matrix $B$, with the probability of correct recovery tending to one as $n \to \infty$.

The statistical convergence rate for the Theil–Sen slope is $n^{-1/2}$ under finite-variance noise, matching the best case for OLS-based DirectLiNGAM in the Gaussian setting. For some discrete-type heavy tails, the rate is even faster. Consequently, the procedure achieves optimal efficiency in well-behaved settings and strictly better performance under heavy tails.

5. Empirical Evidence and Performance Analysis

The empirical evaluation of TSLiNGAM follows simulation protocols akin to Gnecco et al. (2021), using random DAGs with edge weights drawn from $\mathrm{Uniform}([-0.9,-0.1]\cup[0.1,0.9])$, error distributions including Student-$t$, centered log-Normal, centered Pareto, and centered Exponential, and intentional contamination scenarios.

Key metrics include adjacency recovery rate (fraction of simulations with exactly correct causal order) and structural Hamming distance (SHD) between estimated and ground-truth adjacency matrices. For $p = 10$ variables and Student-$t_1$ noise, the following results were recorded (number of correct orderings out of 1000):

sample size	50	100	200	300
DirectLiNGAM	286	432	555	640
TSLiNGAM (TS+KBI)	477	806	942	984
TSLiNGAM (RM+KBI)	379	733	915	969
TSLiNGAM (TS+dcorr)	232	530	747	841
EASE	21	124	307	450

TSLiNGAM (with Theil–Sen slope) outperforms both DirectLiNGAM and the alternative EASE method as the noise becomes heavier-tailed, with the accuracy advantage widening as sample size increases or tail-heaviness intensifies. Under near-Gaussian noises (e.g., Student-$t_5$), performance converges, but TSLiNGAM never underperforms DirectLiNGAM.

Robustness is further validated by single-outlier and small-sample experiments. In a bivariate chain contaminated with a single large outlier, DirectLiNGAM fails at sufficiently high contamination magnitudes, while TSLiNGAM (Theil–Sen or Repeated-Median) remains consistent. For real-world gene expression data subsampled to low $n$, TSLiNGAM exhibits approximately 10% higher causal order recovery than DirectLiNGAM over multiple repeats.

6. Computational and Practical Aspects

The main computational steps include calculation of robust slopes ($O(n \log n)$ per pair using Theil–Sen or Repeated-Median algorithms), and independence estimation using either kernel-based mutual information ($O(n M^2 + M^3)$, with $M \ll n$ for low-rank approximations) or distance correlation ($O(n \log n)$). The overall complexity per iteration is $O(p n \log n)$ for distance correlation, leading to $O(p^2 n \log n)$ total for the full algorithm.

In practice, TSLiNGAM (with Theil–Sen and distance correlation) exhibits runtime comparable to DirectLiNGAM and substantially lower than kernel mutual information-based scoring for large $n$.

Tuning parameters involve:

• Choice of slope functional: Theil–Sen is default for balancing robustness and efficiency. Repeated-Median is preferred under extreme contamination ($\gtrsim 25\%$).
• Independence measure: Distance correlation is preferred for large $n$ ($\gtrsim 500$), while kernel-MI offers maximal sensitivity for moderate $n$.
• Post-hoc sparsification: Adaptive Lasso or thresholding small coefficients in $\widehat B$ (using BIC or cross-validation for percentile selection).

A plausible implication is that the adaptation to robust regression slopes fundamentally broadens the practical scope of LiNGAM-based causal discovery by mitigating the vulnerabilities of OLS to heavy tails and outliers, thus supporting more reliable inference in realistic (non-Gaussian, contaminated) data environments (Leyder et al., 2023).