While constraint-based algorithms leverage conditional independence tests to prune the space of possible causal structures, score-based discovery algorithms take a different route. They formulate causal discovery as an optimization problem: define a score function that measures how well a candidate graph structure fits the observed data, and then search the space of possible graphs for the one that maximizes (or minimizes) this score.

The primary challenge is the vast search space. The number of possible Directed Acyclic Graphs (DAGs) grows super-exponentially with the number of variables, making exhaustive search computationally infeasible for all but the smallest problems. Therefore, score-based methods rely on heuristic or greedy search strategies. We will discuss two prominent approaches: Greedy Equivalence Search (GES) and Linear Non-Gaussian Acyclic Model (LiNGAM).

Greedy Equivalence Search (GES)

Instead of searching the space of individual DAGs, the Greedy Equivalence Search (GES) algorithm cleverly operates on the space of Markov Equivalence Classes (MECs). An MEC represents a set of DAGs that encode the same set of conditional independence relations and, therefore, cannot be distinguished using observational data alone under standard assumptions. Searching over MECs significantly reduces the complexity of the search space.

GES typically proceeds in two phases:

Forward Equivalence Search (FES): Starting with an empty graph (no edges), the algorithm iteratively adds the single edge ( $X \rightarrow Y$ or $Y \rightarrow X$ ) that results in the largest increase in the score function, provided the addition moves the current graph to a different, higher-scoring MEC. This continues until no single edge addition improves the score.
Backward Equivalence Search (BES): Starting from the graph obtained at the end of the FES phase, the algorithm iteratively removes the single edge whose removal results in the largest increase (or smallest decrease) in the score function, ensuring the resulting graph belongs to a valid MEC. This phase continues until no single edge removal improves the score.

The final output of GES is typically represented as a Completed Partially Directed Acyclic Graph (CPDAG), which graphically characterizes the MEC identified.

Score Functions

The choice of score function is fundamental to GES. Common scores include:

Bayesian Information Criterion (BIC) / Schwarz Information Criterion (SIC): $Score_{BIC}(G; D) = \log L(\hat{\theta}_G | D) - \frac{\log n}{2} |G|$ Here, $L(\hat{\theta}_G | D)$ is the maximized likelihood of the data $D$ given the graph $G$ with estimated parameters $\hat{\theta}_G$ , $n$ is the sample size, and $|G|$ is the number of independent parameters in the model (related to graph complexity, often the number of edges plus parameters for each node's conditional distribution). BIC penalizes model complexity more strongly than AIC, often favoring sparser graphs, especially with large datasets.
Akaike Information Criterion (AIC): $Score_{AIC}(G; D) = \log L(\hat{\theta}_G | D) - |G|$ AIC has a smaller penalty term than BIC.
Bayesian Dirichlet equivalent uniform (BDeu) score: Specifically designed for discrete variables, this score calculates the marginal likelihood $P(D|G)$ under assumptions of parameter independence and Dirichlet priors. It possesses the desirable property of score equivalence, meaning all DAGs within the same MEC receive the same score.

Assumptions and Limitations

GES relies on several assumptions:

Acyclicity: The true causal structure is a DAG.
Causal Sufficiency: There are no unobserved common causes (latent confounders) of the observed variables.
Faithfulness: All conditional independencies present in the data correspond to d-separations in the true causal graph.
Distributional Assumptions: Specific scores often imply distributional assumptions (e.g., BIC/AIC are often calculated assuming Gaussian distributions for continuous data or multinomial for discrete data).

GES is statistically consistent under these assumptions, meaning it converges to the true MEC as the sample size increases. However, it performs a greedy search, which means it can get stuck in local optima. Its performance is sensitive to the choice of score function and potential violations of the underlying assumptions, particularly causal sufficiency.

Linear Non-Gaussian Acyclic Model (LiNGAM)

LiNGAM offers a distinct approach that, under specific assumptions, can identify the exact DAG structure, not just the MEC. It achieves this by leveraging information contained in the shape of the data distribution, specifically requiring non-Gaussianity.

The LiNGAM Model

The core assumption of LiNGAM is that the data generating process is linear and acyclic, with non-Gaussian error terms (innovations):

X_i = \sum_{X_j \in PA_i} \beta_{ij} X_j + \epsilon_i

Here, $X_i$ is the $i$ -th variable, $PA_i$ represents the set of its direct causal parents in the DAG, $\beta_{ij}$ are the linear causal coefficients representing the strength of the direct effect from $X_j$ to $X_i$ , and $\epsilon_i$ are independent, non-Gaussian error terms with zero mean and non-zero variance. The acyclicity means that if we order the variables according to the causal structure (e.g., $X_{k_1}, X_{k_2}, ..., X_{k_d}$ where $k_1, ..., k_d$ is a permutation of $1, ..., d$ ), the adjacency matrix $B$ (where $B_{ij} = \beta_{ij}$ if $X_j \in PA_i$ and 0 otherwise) can be made strictly lower triangular by applying the permutation.

This model can be written in matrix form:

X = BX + \epsilon

Which implies:

X = (I - B)^{-1} \epsilon

where $I$ is the identity matrix.

Identification via Non-Gaussianity and ICA

The insight behind LiNGAM comes from Independent Component Analysis (ICA). ICA is a technique used to separate multivariate signals into additive, independent, non-Gaussian source signals. In the LiNGAM context, the observed variables $X$ are linear mixtures of the independent, non-Gaussian error terms $\epsilon$ . The matrix $A = (I-B)^{-1}$ acts as the mixing matrix.

ICA algorithms aim to find an unmixing matrix $W$ such that $W X$ recovers the original independent sources (up to permutation and scaling). That is, $W \approx A^{-1} = (I-B)$ .

LiNGAM leverages the fact that if the correct unmixing matrix $W$ is found, it holds a specific relationship to the causal structure matrix $B$ : $W = P D (I - B)$ , where $P$ is a permutation matrix and $D$ is a diagonal scaling matrix. Critically, one can find a permutation of the rows of $W$ such that the permuted matrix $W'$ has zeros corresponding to the zero entries in $(I-B)$ . Since $B$ is lower triangular in the correct causal order, $(I-B)$ is also lower triangular (with ones on the diagonal). Thus, $W'$ can be permuted to become lower triangular.

The LiNGAM algorithm (in its basic form) typically involves these steps:

Center the data: Subtract the mean from each variable.
Perform ICA: Apply an ICA algorithm (like FastICA) to the centered data $X$ to obtain an estimated unmixing matrix $W$ .
Permutation Search: Find the permutation matrix $P$ such that $W' = PW$ is as close as possible to being lower triangular (e.g., minimizing the sum of squares of the elements above the diagonal after normalizing rows).
Estimate B: Reorder the variables according to the permutation $P$ . Estimate the connection strengths $B$ from the permuted, scaled unmixing matrix. $B$ can be derived from $W'$ using $B = I - (W')^{-1}$ (after appropriate scaling). Alternatively, estimate connections via regression based on the identified causal order.
Pruning (Optional): Perform statistical tests or use bootstrapping to prune weak or insignificant edges in the estimated $B$ .

Assumptions and Strengths

LiNGAM's main strengths are its ability to identify the full DAG (including edge directions) and estimate causal strengths under its assumptions:

Linearity: Relationships between variables are linear.
Acyclicity: No feedback loops.
Non-Gaussianity: All independent error terms $\epsilon_i$ must be non-Gaussian. At most one error term can be Gaussian.
Causal Sufficiency: No latent confounders.

If these assumptions hold, LiNGAM provides a powerful tool for discovery. Various extensions exist to relax some assumptions, such as handling time series data (VAR-LiNGAM) or improving computational efficiency (DirectLiNGAM).

Limitations

The strict linearity and non-Gaussianity requirements are LiNGAM's primary limitations. It fundamentally fails if the underlying error terms are Gaussian, as ICA cannot uniquely determine the rotation in this case. Performance also depends on the quality of the ICA estimation and the permutation search. Like GES, it assumes causal sufficiency.

A conceptual illustration of the LiNGAM data generating process. Independent non-Gaussian sources ( $\epsilon_i$ ) influence observed variables ( $X_i$ ) through a linear acyclic structure defined by coefficients ( $\beta_{ij}$ ). LiNGAM aims to recover this structure from the observed $X_i$ .

Choosing Between Score-Based Methods

GES is more general regarding functional forms (depending on the score) but typically requires Gaussian assumptions for standard scores like BIC/AIC calculated via linear models. It identifies the MEC, not necessarily the full DAG, and requires causal sufficiency.
LiNGAM identifies the full DAG but imposes stricter linearity and non-Gaussianity assumptions. It also requires causal sufficiency.

If your data is strongly believed to be generated by a linear process and exhibits non-Gaussian characteristics, LiNGAM is a compelling option for potentially recovering the full DAG structure. If linearity is questionable or data appears Gaussian, GES (with an appropriate score) might be more suitable, acknowledging that it will only identify the MEC.

Implementation Notes

Libraries like causal-learn provide implementations for GES and related algorithms, often allowing specification of different score functions (e.g., BIC for Gaussian data). The lingam Python package offers various implementations of LiNGAM and its extensions (DirectLiNGAM, VAR-LiNGAM, etc.).

When applying these methods, remember their sensitivity to assumptions. Assess data distributions (e.g., using normality tests) before choosing LiNGAM. For GES, consider the implications of the chosen score function and potential sensitivity to local optima. Validating discovered structures using domain knowledge or interventional data (if available) is always advisable. Score-based methods, while optimizing a global criterion, can still be influenced by data quality and the specific heuristic search strategy employed.