Dimension Reduction Speeds Up Causal Graph Discovery

Pillar 5, graph learning and probabilistic methods. Analysis by the aitrendblend editorial team. Reading time about 14 minutes.

Why finding cause and effect from a spreadsheet is so hard

Most of the data scientists and researchers ever collect is observational. Nobody randomized which patients got which gene variant, which households ended up with which income, or which neurons fired before which others. All you get is a snapshot of variables sitting next to each other, and from that snapshot you are supposed to reconstruct which variable actually influences which. The standard mathematical object for representing that structure is a directed acyclic graph, a DAG, where an arrow from variable i to variable j means i has a direct causal effect on j and where the graph never loops back on itself.

The trouble is that many different DAGs can produce the exact same pattern of statistical dependencies in the data. Two DAGs that cannot be told apart by any conditional independence test are called Markov equivalent, and the most anyone can honestly claim to recover from pure observational data is the equivalence class, not a single definitive graph. That equivalence class is usually drawn as a completed partially directed graph, a CPDAG, where some edges have a clear direction and others are left undirected because the data genuinely cannot settle the question.

The classic tool for estimating this structure is the PC algorithm, introduced by Spirtes, Glymour, and Scheines and later shown by Kalisch and Buhlmann to scale well with network size as long as the true graph is sparse. It works by repeatedly asking whether two variables are conditionally independent given various subsets of the other variables, and removing an edge whenever it finds a subset that breaks the dependence. That question, is variable i independent of variable j given some set S, is easy to answer under a Gaussian assumption using partial correlations. It gets much harder once you refuse to assume Gaussianity, because a fully general conditional independence test needs to condition on a potentially large, high dimensional set of variables, and nonparametric tests degrade badly as that conditioning set grows. This is the curse of dimensionality the paper is built to address.

What the paper actually proposes

The authors, Eftychia Solea, Bing Li, and Kyongwon Kim, split the estimation problem into two stages that they run for every candidate edge and every candidate conditioning set the PC algorithm considers.

Stage one, shrink the conditioning set

Instead of testing independence conditional on the raw, possibly high dimensional set of variables X_S, the method first runs a technique called generalized sliced inverse regression, GSIR, to find a small number of nonlinear functions of X_S that capture everything relevant to the relationship between the pair of variables being tested. This is nonlinear sufficient dimension reduction, an extension of the classical sliced inverse regression method by Ker-Chau Li into reproducing kernel Hilbert spaces, developed by Lee, Li, and Chiaromonte. The output is a low dimensional vector, often just one or two numbers, that the paper calls the sufficient predictor for that particular pair and conditioning set.

Stage two, test independence in the reduced space

Once the conditioning set has been compressed, the method tests conditional independence using an operator called the conjoined conditional cross covariance operator, CCCO, originally proposed by Fukumizu and coauthors for kernel based independence testing. The Hilbert Schmidt norm of this operator is zero if and only if the two variables are conditionally independent given the reduced predictor, under a set of regularity assumptions the paper lays out carefully. Thresholding that norm against a permutation based cutoff gives a practical independence test that only ever has to work with a handful of reduced dimensions rather than the full conditioning set.

The two stages together define what the authors call the DAG-PC algorithm, a modified version of the classical PC algorithm where each conditional independence test is preceded by this dimension reduction step. Crucially, the paper proves that the computational complexity of the whole procedure still depends on how sparse the true graph is rather than on the total number of variables, which is the same appealing property that made the original PC algorithm practical for large networks in the first place.

The math underneath

The central object the paper is trying to characterize is a sub sigma field, not just a fixed set of variables. Assumption 1 in the paper states that for any conditioning set X_S, there exists a smallest sub sigma field of X_S, denoted G_S, such that conditioning on G_S carries exactly the same information about the pair (X_i, X_j) as conditioning on the whole of X_S.

This is the theoretical basis that lets the method replace a high dimensional conditioning set with a low dimensional summary without losing any information relevant to conditional independence. The regression operator that the GSIR step estimates is defined as a composition of covariance operators on reproducing kernel Hilbert spaces.

At the sample level, this turns into a generalized eigenvalue problem solved with Tikhonov regularized matrix inverses, since the operator inverse would otherwise be unbounded.

Once the sufficient predictor has been extracted, the conjoined conditional cross covariance operator itself is defined as an ordinary conditional covariance, but between the joint pairs (X_i, X_S) and (X_j, X_S) rather than between X_i and X_j directly.

The key theoretical payoff, proved as Proposition 4 following Fukumizu and coauthors, is that this operator is exactly zero if and only if X_i and X_j are conditionally independent given X_S, under mild conditions on the kernels involved. Thresholding its Hilbert Schmidt norm against a data driven cutoff is what actually decides whether an edge stays or goes in the estimated skeleton.

What the simulations actually show

The authors test four methods across three data generating regimes, Gaussian, non Gaussian, and a more elaborate mixed non Gaussian setting with random parent interactions and randomly sampled error variances. The four methods are the Gaussian PC algorithm based on partial correlations, a rank correlation based PC algorithm for the nonparanormal case, a fully nonparametric kernel PC algorithm with no dimension reduction, and the paper’s own DAG-PC method. They measure two things, the structural Hamming distance between the estimated and true CPDAG, where lower is better, and the area under the ROC curve for skeleton recovery, where higher is better.

Table 5, reproduced above from the paper’s own results, is where the method’s real advantage shows up. Across non Gaussian data, DAG-PC generally beats the two linear baselines by a wide margin and holds a smaller but consistent edge over the fully nonparametric kernel method, especially as p grows toward 50. That pattern matches the paper’s core argument reasonably well.

The Gaussian setting tells a different story. The paper’s own Table 1 reports AUC values where, at n equal 100 and p equal 5 under E[N] equal 2, the Gaussian PC method scores 0.73, the rank based method scores 0.71, the kernel method scores 0.59, and DAG-PC scores 0.64. DAG-PC beats the fully nonparametric kernel method but sits behind both linear methods, which makes sense since the data really is Gaussian in that scenario and the linear methods are correctly specified for it. The authors are upfront about this in their own discussion of Figure 2, writing that Methods A, B, and C perform slightly better than the DAG method in the E[N] equal 4 case. That is a meaningful caveat sitting a few paragraphs after a comparison table that otherwise reads as a clean win.

One more detail worth flagging for anyone who goes back to the source. In the paper’s printed Table 1, two entries for the kernel PC and DAG-PC methods at n equal 100 read as 0.61(.0.3) and 0.64(0.3), where every other standard deviation in the same table is a two digit value like 0.06 or 0.02. Those two entries almost certainly contain a decimal point typo and were meant to read 0.61(0.03) and 0.64(0.03) to match the surrounding pattern. It is a minor production error rather than a substantive issue, but worth knowing if you are citing the exact numbers from that table.

The real data application, and why the win there is smaller than it looks

Beyond simulation, the authors apply all four methods to the Sachs et al flow cytometry dataset, a well known benchmark of eleven phosphorylated proteins and phospholipids measured in single cells, with a causal structure that has been independently validated by prior biological research. They subsample the data ten times to 90 observations, run all four methods, and report the mean and standard deviation of the structural Hamming distance against the accepted true graph.

DAG-PC does post the lowest mean structural Hamming distance, and the authors describe this as evidence of competitiveness for investigating causal relations among cells. Read the numbers themselves, though, and the picture softens. DAG-PC’s mean of 20.6 sits well inside one standard deviation of the kernel PC method’s 21.1, and not much further from the Gaussian PC method’s 21.2. With only ten subsampling repetitions and overlapping error bars this size, the honest description is that DAG-PC is roughly tied with the two nearest competitors on this particular dataset, with a small numerical edge that the paper does not test for statistical significance. That is a legitimate, useful result. It is a considerably smaller claim than what a bare comparison table implies at a glance.

A 20.6 average sitting inside one standard deviation of a 21.1 average is a tie dressed up as a win, and the paper never runs the significance test that would settle it either way. Editorial analysis of Table 7, Solea, Li, and Kim, JMLR 2025

What the paper does not compare against

The three baselines in this study, the Gaussian PC algorithm, a rank correlation PC algorithm, and a fully nonparametric kernel PC algorithm, are all constraint based methods from the same family as the paper’s own proposal. That is a fair and useful comparison for isolating the specific contribution of dimension reduction within that family. It leaves out an entire other branch of causal discovery research that has grown substantially in recent years, the continuous optimization and score based methods that formulate DAG learning as a differentiable constrained optimization problem, along with structural equation model approaches built around non Gaussian noise assumptions like LiNGAM, which the paper cites in its related work section but never runs as a numerical baseline. A reader trying to decide which tool to reach for in 2026 would want to know how DAG-PC stacks up against that broader landscape, not just against other constraint based methods, and the paper as published does not answer that question.

A demanding stack of assumptions

The theoretical guarantees in the paper, consistency of the estimated skeleton, uniform consistency of the estimated CPDAG, and the associated convergence rates, rest on eleven numbered assumptions spread across the methodology and theory sections. Individually each one is a standard and defensible condition in the sufficient dimension reduction and kernel operator literature, things like the completeness of the central class, finite rank regression operators, and bounded eigenvalues of certain kernel Hessians. Collectively, though, verifying that a real dataset satisfies all eleven is not something a practitioner can easily check, and the paper is honest that strong faithfulness itself, the condition underlying the uniform consistency result, is known from prior work by Uhler and coauthors to be a demanding assumption that fails for many DAGs. The gap between what the theory requires and what an analyst can confirm in practice is a real limitation, and it is one shared by most of the theoretical causal discovery literature, not a flaw unique to this paper.

Honest limitations

Beyond the points raised above, a few constraints are worth naming plainly. Choosing the dimension of the sufficient predictor, the number d_S^ij in the paper’s notation, has no principled data driven procedure attached to it. The authors acknowledge this directly and simply report that a value of 1 or 2 worked well in their simulations, which is a reasonable default but not a solved problem. The tuning of the various regularization constants, eta, epsilon, and delta, relies on a generalized cross validation scheme and a spectral cutoff at 99 percent of cumulative variation, both sensible choices that nonetheless add several more knobs an analyst has to trust are behaving well on their own data. Finally, every simulation and the real data example in this paper max out at p equal 50 variables. Modern genomic and neuroimaging applications routinely involve hundreds or thousands of variables, and while the theory promises complexity scaling with sparsity rather than network size, that promise has not yet been demonstrated numerically at the scale where it would matter most.

The implementation

The code below reconstructs the two core estimation steps described in Section 4 of the paper, the GSIR sufficient predictor via the regularized eigenvalue problem in Equation 14, and the conjoined conditional cross covariance operator via the coordinate representation in Section 4.3, along with a simplified version of the DAG-PC skeleton search from Algorithm 1. This is a nonparametric statistical estimator built from kernel matrix operations rather than a trained neural network, so there is no gradient based training loop in the usual deep learning sense. In its place, the code below has an estimation loop that plays the equivalent role, iterating the GSIR and CCCO steps over candidate edges and conditioning sets, exactly as the paper’s Algorithm 1 describes. A permutation test stands in for the paper’s data driven thresholding procedure, and a runnable smoke test at the end generates synthetic data from the same non Gaussian model the paper itself uses in Section 6 and confirms the full pipeline executes end to end.

import torch

torch.manual_seed(0)

# ---------------------------------------------------------------------
# Kernel utilities
# All Gram matrices use the Gaussian RBF kernel and the bandwidth rule
# in Equation 19 of the paper.
# ---------------------------------------------------------------------

def gaussian_rbf_gram(X, gamma):
    sq = torch.cdist(X, X, p=2) ** 2
    return torch.exp(-gamma * sq)

def center_gram(K):
    n = K.shape[0]
    Q = torch.eye(n) - torch.ones(n, n) / n
    return Q @ K @ Q

def median_heuristic_gamma(X):
    # Equation 19, the bandwidth rule the paper uses for every kernel.
    n = X.shape[0]
    dists = torch.cdist(X, X, p=2)
    iu = torch.triu_indices(n, n, offset=1)
    mean_dist = torch.clamp(dists[iu[0], iu[1]].mean(), min=1e-6)
    return 1.0 / (mean_dist ** 2)

def tik_inverse(G, reg):
    n = G.shape[0]
    return torch.linalg.inv(G + reg * torch.eye(n))


# ---------------------------------------------------------------------
# Step 1, GSIR
# Solves the eigenvalue problem in Equation 14 to obtain the sufficient
# predictor that summarizes a high dimensional conditioning set X_S.
# ---------------------------------------------------------------------

def gsir_sufficient_predictors(X_S, X_i, X_j, d_reduced, eta, eps):
    n = X_S.shape[0]
    Q = torch.eye(n) - torch.ones(n, n) / n
    gamma_S = median_heuristic_gamma(X_S)
    gamma_ij = median_heuristic_gamma(torch.cat([X_i, X_j], dim=1))
    K_S_raw = gaussian_rbf_gram(X_S, gamma_S)
    G_S = center_gram(K_S_raw)
    G_ij = center_gram(gaussian_rbf_gram(torch.cat([X_i, X_j], dim=1), gamma_ij))

    inv_S = tik_inverse(G_S, eta)
    inv_ij = tik_inverse(G_ij, eps)
    M = inv_S @ G_S @ G_ij @ inv_ij @ G_S @ inv_S
    M = 0.5 * (M + M.T)
    eigvals, eigvecs = torch.linalg.eigh(M)
    order = torch.argsort(eigvals, descending=True)
    top_v = eigvecs[:, order[:d_reduced]]
    U = (top_v.T @ inv_S @ (K_S_raw @ Q)).T
    return U


# ---------------------------------------------------------------------
# Step 2, conjoined conditional cross covariance operator
# Follows the coordinate representation in Section 4.3 and returns the
# Hilbert Schmidt norm used as the conditional independence statistic.
# ---------------------------------------------------------------------

def ccco_hs_norm(X_i, X_j, U, delta):
    n = U.shape[0]
    iu_pair = torch.cat([X_i, U], dim=1)
    ju_pair = torch.cat([X_j, U], dim=1)
    K_iu = center_gram(gaussian_rbf_gram(iu_pair, median_heuristic_gamma(iu_pair)))
    K_ju = center_gram(gaussian_rbf_gram(ju_pair, median_heuristic_gamma(ju_pair)))
    K_u = center_gram(gaussian_rbf_gram(U, median_heuristic_gamma(U)))
    inv_u = tik_inverse(K_u, delta)
    ccco = K_ju - K_u @ inv_u @ K_ju
    return (torch.linalg.norm(K_iu @ ccco, ord='fro') / n).item()

def unconditional_hs_norm(X_i, X_j):
    # A degenerate CCCO with no conditioning set, equivalent to the
    # Hilbert Schmidt norm of the plain cross covariance operator.
    K_i = center_gram(gaussian_rbf_gram(X_i, median_heuristic_gamma(X_i)))
    K_j = center_gram(gaussian_rbf_gram(X_j, median_heuristic_gamma(X_j)))
    return (torch.linalg.norm(K_i @ K_j, ord='fro') / X_i.shape[0]).item()

def permutation_threshold(stat_fn, n, n_perm=12, alpha=0.1):
    # The permutation test based thresholding procedure from Section 4.4.
    stats = [stat_fn(torch.randperm(n)) for _ in range(n_perm)]
    stats.sort()
    idx = max(0, int((1 - alpha) * len(stats)) - 1)
    return stats[idx]


# ---------------------------------------------------------------------
# The DAG-PC algorithm
# A simplified reproduction of Algorithm 1. For tractability in a smoke
# test we search conditioning sets up to a small maximum size rather
# than the full adaptive search the paper describes, but the
# conditional independence test itself, GSIR followed by the CCCO, is
# the real estimating equation from the paper.
# ---------------------------------------------------------------------

def dag_pc_skeleton(X, d_reduced=1, eta=0.05, eps=0.05, delta=0.05,
                     alpha=0.1, max_cond_size=1):
    n, p = X.shape
    adj = torch.ones(p, p) - torch.eye(p)

    for i in range(p):
        for j in range(i + 1, p):
            if adj[i, j] == 0:
                continue
            X_i, X_j = X[:, i:i + 1], X[:, j:j + 1]

            # Order 0, no conditioning set.
            hs0 = unconditional_hs_norm(X_i, X_j)
            rho0 = permutation_threshold(
                lambda perm: unconditional_hs_norm(X_i, X_j[perm]), n
            )
            if hs0 < rho0:
                adj[i, j] = 0
                adj[j, i] = 0
                continue

            # Higher order tests, conditioning on neighboring nodes.
            neighbors = [k for k in range(p) if k != i and k != j and adj[i, k] == 1]
            for cond_size in range(1, min(max_cond_size, len(neighbors)) + 1):
                if adj[i, j] == 0:
                    break
                for S in [neighbors[:cond_size]]:
                    X_S = X[:, S]
                    U = gsir_sufficient_predictors(X_S, X_i, X_j, d_reduced, eta, eps)
                    hs = ccco_hs_norm(X_i, X_j, U, delta)
                    rho = permutation_threshold(
                        lambda perm: ccco_hs_norm(X_i, X_j[perm], U[perm], delta), n
                    )
                    if hs < rho:
                        adj[i, j] = 0
                        adj[j, i] = 0
                        break
    return adj


# ---------------------------------------------------------------------
# Smoke test
# Generates synthetic data from the paper's own non Gaussian generative
# model, Step I through III in Section 6, then estimates the skeleton.
# ---------------------------------------------------------------------

if __name__ == "__main__":
    n, p = 150, 5
    D = torch.zeros(p, p)
    for i in range(p):
        for parent in range(i):
            if torch.rand(1).item() < 0.4:
                D[parent, i] = torch.rand(1).item() * 1.5 + 0.3

    noise = torch.randn(n, p)
    X = torch.zeros(n, p)
    X[:, 0] = noise[:, 0]
    for i in range(1, p):
        contribution = torch.zeros(n)
        for parent in range(i):
            if D[parent, i] != 0:
                contribution += D[parent, i] * torch.cos(X[:, parent])
        X[:, i] = contribution + noise[:, i]

    print("Running DAG-PC skeleton estimation smoke test on synthetic non Gaussian data")
    adj = dag_pc_skeleton(X, d_reduced=1, max_cond_size=1)
    print("True adjacency, from the generating DAG")
    print((D != 0).int())
    print("Estimated skeleton adjacency")
    print(adj.int())
    print("Smoke test complete, the GSIR, CCCO, and DAG-PC pipeline runs without errors")

What this means for anyone building on the paper

If the data you are working with is genuinely close to Gaussian, this paper’s own results suggest reaching for the classic Gaussian PC algorithm or the rank based variant first, since both outperformed DAG-PC in that regime in the paper’s own tables. If your data shows real nonlinear or non Gaussian structure, which is common in biological, financial, and behavioral measurements, the simulation evidence for DAG-PC over the fully nonparametric kernel alternative is more consistent and gets stronger as the number of variables grows, which is exactly where dimension reduction should help most. Before trusting a single real data comparison the way Table 7 might tempt you to, check whether the reported margin sits inside the reported standard deviation, the way DAG-PC’s edge over the kernel method does here. And if you plan to compare this family of constraint based methods against newer continuous optimization approaches to causal discovery, budget for running those baselines yourself, since this paper does not include them.

Conclusion

The core methodological idea in this paper holds up well under scrutiny. Compressing a high dimensional conditioning set into a low dimensional sufficient predictor before testing conditional independence is a genuinely elegant way to sidestep the curse of dimensionality that has limited fully nonparametric causal discovery methods for years, and the asymptotic theory backing it, consistency of the estimated skeleton and CPDAG with explicit convergence rates, is carefully worked out across eleven numbered assumptions and several dozen pages of proofs in the appendix.

The conceptual shift this represents, connecting nonlinear sufficient dimension reduction directly to the classical PC algorithm rather than relying on additive structure as prior work by Lee and coauthors did, opens a path toward nonparametric causal discovery that could plausibly extend well beyond the specific operator used here. The authors themselves point toward extensions using the normalized version of the CCCO, sparser penalty based thresholding through LASSO type methods, and eventually functional data settings relevant to electroencephalography and functional MRI, all of which sound like natural next steps built on a sound foundation.

Whether this specific implementation transfers well to other domains depends on how close real world data actually gets to the non Gaussian, moderately sparse regime where the simulations show DAG-PC’s clearest advantage. In settings closer to Gaussian, this paper’s own numbers say a simpler classical method may serve just as well or better. In settings with hundreds or thousands of variables, the numerical evidence for the promised sparsity based scaling has not yet been shown, only argued theoretically.

The honest remaining limitations are the ones already discussed. No principled method exists yet for choosing the dimension of the sufficient predictor. The real data application shows a margin of victory that overlaps with its nearest competitor’s uncertainty. And the comparison set leaves out an entire branch of modern causal discovery research built on continuous optimization, which any practitioner surveying the field in 2026 would want to see addressed.

None of that erases the value of the contribution. This is a mathematically careful paper that identifies a real gap, fully nonparametric conditional independence testing does not scale, and offers a principled, theoretically grounded fix. The next useful step is exactly what the discussion section gestures toward but does not yet deliver, a numerical demonstration at genuinely large p, a head to head comparison against continuous optimization based causal discovery methods, and a statistical significance test on the real data comparison that currently rests on ten subsamples and overlapping error bars. Until those experiments exist, the fair way to describe this paper is that it built a theoretically sound and computationally efficient nonparametric causal discovery method, and it has shown that method winning clearly in some regimes and only modestly, or not at all, in others.

Frequently asked questions

What problem is this paper actually solving

It estimates the causal structure behind observational data, represented as a directed acyclic graph, without assuming the data follows a Gaussian distribution and without the usual curse of dimensionality that slows down other nonparametric methods when the conditioning set gets large.

How does the method avoid the curse of dimensionality

It first applies nonlinear sufficient dimension reduction, called GSIR, to compress a large conditioning set down to a small number of summary values, then tests conditional independence in that reduced space using a kernel operator called the conjoined conditional cross covariance operator.

Does this method always outperform simpler causal discovery methods

No. In the paper’s own Gaussian data simulations, the classic Gaussian PC algorithm and a rank correlation based method sometimes outperform this new method, since those simpler methods are correctly specified for Gaussian data. The advantage of the new method shows up most clearly in non Gaussian and mixed non Gaussian settings.

How strong is the real data result on the flow cytometry dataset

The new method achieves the lowest average structural Hamming distance, 20.6 against 21.1 for the nearest competing method, but that margin sits inside one standard deviation of the competing method’s result, and the paper does not report a formal significance test on the difference.

Was this method compared against modern deep learning based causal discovery approaches

No. The paper compares only against other constraint based PC algorithm variants, a Gaussian version, a rank correlation version, and a fully nonparametric kernel version. It does not benchmark against continuous optimization or score based causal discovery methods that have become common in recent years.

What is the largest number of variables this method has been tested on

The simulations and the real data application in this paper go up to 50 variables. The theory argues that computational complexity should scale with the sparseness of the true graph rather than the number of variables, but that claim has not yet been demonstrated numerically at a much larger scale.

The full methodology, including the complete reference list and the appendix proofs, is available in the published paper. You can read it directly at the Journal of Machine Learning Research page for volume 26, paper 24-0048.

Related reading