Discrete Score Matching Enables Causal Discovery from Count Data

Euijong Song; Hyewon Park; Gunwoong Park

stat.MLSep 30, 2026

Discrete Score Matching Enables Causal Discovery from Count Data

Authors: Euijong Song, Hyewon Park, Gunwoong Park

Organizations: Department of Statistics Seoul National University, Korea · Interdisciplinary Program in Artificial Intelligence; Institute for Data Innovation in Science, Seoul National University, Korea

Abstract

Count data pose a challenge for score-matching-based causal discovery: derivatives are unavailable, and simply replacing them with finite differences does not generally suffice for causal discovery. We generalize SCORE's constant-curvature criterion (Rolland et al., 2022) by conditioning on the node's value, yielding the conditional curvature score (CCS) for ordering. We also extend curvature-based parent recovery through the off-diagonal curvature score (OCS), enabling directed acyclic graph (DAG) recovery with both scores constructed from score functions for continuous data and concrete scores for counts. In the bivariate setting, zero CCS exactly characterizes a semiparametric generalized linear model (GLM) conditional form in which the conditional family need not be specified in advance, unlike in classical GLMs. For bivariate semiparametric GLM DAGs under our regularity condition, canonical-parameter nonlinearity is necessary and sufficient for identifiability. In multivariate DAGs, this nonlinearity enables DAG recovery through CCS and OCS. Our framework identifies a new class of semiparametric GLM DAGs that strictly contains the nonlinear Gaussian ANM class identified by SCORE. We introduce DISCO (DIscrete SCOre), a count-DAG recovery algorithm that estimates CCS and OCS using discrete diffusion. Experiments demonstrate accurate DAG recovery across Poisson, negative binomial, binomial, and mixed-family settings, as well as scalability to 1,000-node DAGs on a single GPU.

Figures & tables

Conditional $X_{j}\mid x_{\operatorname{Pa}_{G}(j)}$	Conditional family	Support	Link $g$
$\mathrm{Poisson}(\exp f)$	Poisson	$\mathbb{N}_{0}$	$\log$
$\mathrm{Poisson}(\operatorname{softplus}f)$	Poisson	$\mathbb{N}_{0}$	$\operatorname{softplus}^{-1}$
$\mathrm{NB}(r,\exp f)$	NB $(r)$	$\mathbb{N}_{0}$	$\log$
$\mathrm{Binomial}(M,\operatorname{sigmoid}(f))$	Binomial $(M)$	$\{0,\dots,M\}$	$\mu\mapsto\operatorname{logit}(\mu/M)$
$\mathcal{N}(f,\tau^{2})$	Gaussian	$\mathbb{R}$	identity
$\mathrm{Gamma}(\alpha,\exp f)$	Gamma $(\alpha)$	$(0,\infty)$	$\log$

Table 1: Example conditional distributions, with index

f=f_{j}(x_{\operatorname{Pa}_{G}(j)})

and logistic function

\operatorname{sigmoid}

. Noncanonical links can yield nonlinear

\eta

even for affine

f

. NB and Gamma use the mean as their second parameter, with fixed size

r

and shape

\alpha

, respectively.

Figure 1: End-to-end

F_{1}

versus sample size on DAGs with

d=100

: (a) NB with

r=1

; (b) NB with

r=6

; and (c) the mixed NB

\to

Bin

\to

Poi configuration. Points and error bars show the median and interquartile range. ODS (Oracle) uses the true nodewise conditional families and their fixed parameters.

Family	Constant	CCS	Gain
Poi	68.3	81.7	$+13.3$
NB	73.3	93.3	$+20.0$
Bin	80.0	83.3	$+3.3$

Table 2: Ordering and parent-selection comparisons (

d=50

n=5000

). (a) Sink-selection accuracy: Constant denotes the constant-curvature criterion; gain is in percentage points. (b) Mean

F_{1}

given the same estimated DISCO order.

Method	$F_{1}$	Total time (min)	$F_{1}$	Total time (min)	$F_{1}$	Total time (min)	$F_{1}$	Total time (min)
	$d=100$		$d=200$		$d=500$		$d=1000$
DISCO	0.915	2.7	0.835	3.5	0.769	10.8	0.769	57.4
ODS (Oracle)	0.725	12.5	0.713	37.0	0.729	204.7	–	Timeout
MRS	0.786	3.5	0.763	10.8	0.742	96.3	–	Timeout
NOTEARS-MLP	0.467	174.1	–	Timeout	–	Timeout	–	Timeout
DiffAN	–	Timeout	–	Timeout	–	Timeout	–	Timeout

Table 3: Scalability on Poisson DAGs (

n=5000

; numerical entries are medians). Timeout : runtime exceeds 240 minutes.

Appendix figures & tables18 assets

Supplementary material from the paper’s appendix.

Appendix

Symbol	Meaning (defined in)
$G=(V,E)$ , $\|V\|=d$	DAG and number of nodes (Sec. 2.2 )
$\widehat{G}$	estimated DAG (Sec. 6.3 )
$\operatorname{Pa}_{G}(j),\ \operatorname{Ch}_{G}(j)$	parents and children of $j$ (a sink has $\operatorname{Ch}_{G}(j)=\emptyset$ )
source node	a node $j$ with $\operatorname{Pa}_{G}(j)=\emptyset$
$R\subseteq V$	remaining node set, ancestral when obtained by recursive sink removal
$G_{R},\ P_{R},\ p_{R}$	induced subgraph, marginal distribution of $X_{R}$ , and its density

Appendix

Table 4: Recurring symbols.

Conditional family	Conditional response	$\phi(x)$	$b_{0}$	Coefficients
Poi	$\mu_{j}=0.2+\operatorname{softplus}(z_{j})$	$x$	$0.5$	$[0.30,0.45]$
Poi	$\mu_{j}=\exp(\widetilde{z}_{j})$	$\log(1+x)$	$0.4$	$[0.15,0.30]$
NB	$\mu_{j}=0.2+\operatorname{softplus}(z_{j})$	$x$	$0.5$	$[0.30,0.45]$
NB	$\mu_{j}=\exp(\widetilde{z}_{j})$	$\log(1+x)$	$0.4$	$[0.15,0.30]$
Bin	$q_{j}=\Phi(z_{j})$	$x/50$	$-1.28$	$[0.50,1.00]$
Bin	$q_{j}=\operatorname{sigmoid}(z_{j})$	$\sqrt{x/50}$	$-2.197$	$[0.60,0.90]$

Appendix

Table 5: The six single-family benchmark mechanisms. Coefficients are drawn independently and uniformly from the stated intervals.

Conditional family	Criterion	$A_{\mathrm{top}}$	$F_{1}$	SHD
Poi	Constant	$0.915\pm 0.022$	$0.883\pm 0.037$	$31.5\pm 9.4$
Poi	CCS	$0.931\pm 0.017$	$0.892\pm 0.034$	$28.9\pm 8.7$
NB	Constant	$0.907\pm 0.019$	$0.763\pm 0.070$	$57.9\pm 15.5$
NB	CCS	$0.931\pm 0.017$	$0.792\pm 0.066$	$51.2\pm 15.4$
Bin	Constant	$0.899\pm 0.028$	$0.884\pm 0.027$	$30.8\pm 7.7$
Bin	CCS	$0.901\pm 0.027$	$0.877\pm 0.031$	$32.3\pm 8.4$

Appendix

Table 6: DAG recovery with CCS and the constant-curvature criterion (

d=50

n=5000

). Entries are means

\pm

sample standard deviations over 10 seeds. Constant denotes the constant-curvature criterion.

	Raw counts		$\log(1+x)$		Anscombe		$x^{2}+x$
Method	$A_{\mathrm{top}}$	$F_{1}$	$A_{\mathrm{top}}$	$F_{1}$	$A_{\mathrm{top}}$	$F_{1}$	$A_{\mathrm{top}}$	$F_{1}$
DISCO	0.922	0.894	0.922	0.894	0.922	0.894	0.922	0.894
DISCO-NoRank	0.937	0.899	0.897	0.726	0.890	0.860	0.834	0.301

Appendix

Table 7: Invariance to strictly increasing transformations on Poisson DAGs (

d=50

n=5000

). Entries are medians. DISCO-NoRank omits the rank transform.

Figure 2: Extended single-family comparison across all six mechanisms and six sample sizes. Columns correspond to the Poi, NB, and Bin conditional families; the top and bottom rows use the first and second mechanisms for each family in Table 5 , respectively. Points and error bars show the median and interquartile range.

Figure 3: Mixed-family DAG recovery across six conditional-family configurations. Each panel compares DISCO with ODS (Oracle), ODS (Poi), ODS (NB), and ODS (Bin). ODS (Oracle) uses the true nodewise conditional families and their fixed parameters. Points and error bars show the median and interquartile range.

Method	$d$	$A_{\mathrm{top}}$
DISCO	100	0.937
	200	0.902
	500	0.914
	1000	0.908
ODS (Oracle)	100	0.819
	200	0.816

Appendix

Table 8: Ordering accuracy in the scalability experiments at

n=5000

. Entries are medians; configurations marked Timeout in Table 3 are omitted.

Setting	Conditional family	Index/configuration	$A_{\mathrm{top}}$	$F_{1}$	SHD
Single- family	Poi	Affine	0.816	0.750	49.0
	Poi	Nonlinear	0.897	0.873	25.0
	NB	Affine	0.841	0.788	42.5
	NB	Nonlinear	0.907	0.884	22.0
	Bin	Affine	0.886	0.860	28.5
	Bin	Nonlinear	0.904	0.886	23.0

Appendix

Table 9: Single-family and mixed-family DAG recovery on ER1 graphs (

d=100

n=5000

). Entries are medians for DISCO .

Graph	Conditional family	$A_{\mathrm{top}}$	$F_{1}$	SHD
SF1	Poi	$0.784\pm 0.065$	$0.705\pm 0.091$	$28.6\pm 8.5$
	NB	$0.822\pm 0.054$	$0.741\pm 0.094$	$24.8\pm 8.9$
	Bin	$0.778\pm 0.060$	$0.782\pm 0.043$	$20.1\pm 4.0$
SF3	Poi	$0.918\pm 0.019$	$0.719\pm 0.074$	$65.3\pm 14.4$
	NB	$0.904\pm 0.035$	$0.620\pm 0.114$	$82.4\pm 19.7$
	Bin	$0.769\pm 0.050$	$0.693\pm 0.077$	$71.2\pm 14.8$

Appendix

Table 10: Scale-free count-DAG recovery (

d=50

n=5000

). Entries are means

\pm

sample standard deviations.

Conditional family	$A_{\mathrm{top}}$	$F_{1}$	SHD
COM–Poisson ( $\nu=0.5$ )	$0.934\pm 0.018$	$0.853\pm 0.033$	$37.8\pm 8.5$
COM–Poisson ( $\nu=1$ )	$0.928\pm 0.020$	$0.877\pm 0.043$	$32.4\pm 10.2$
COM–Poisson ( $\nu=2$ )	$0.920\pm 0.019$	$0.873\pm 0.026$	$33.6\pm 6.5$

Appendix

Table 11: Recovery beyond the three benchmark conditional families (

d=50

n=5000

). Entries are means

\pm

sample standard deviations. The

\nu=1

row uses results from the corresponding Poisson experiment; graph metrics include DISCO parent selection.

Graph	Method	Gaussian	Gamma
ER1	DISCO-Cont	$0.877\pm 0.063$	$0.640\pm 0.126$
ER1	DiffAN	$0.850\pm 0.072$	$0.540\pm 0.137$
ER3	DISCO-Cont	$0.866\pm 0.041$	$0.718\pm 0.021$
ER3	DiffAN	$0.882\pm 0.049$	$0.657\pm 0.057$

Appendix

Table 12: Continuous-data ordering recovery (

d=30

n=5000

). Entries report

A_{\mathrm{top}}

(higher is better) as means

\pm

sample standard deviations.

Panel of Figure 4 — Figure 4: Ordering and parent-edge diagnostics on the Lahman 2019–2025 q4 representation. The top row shows ODS with Poisson and negative-binomial specifications; the bottom row shows ODS with a binomial specification and DISCO . Blue arrows agree with the reference direction and orange arrows reverse it. Arrow labels report ordering support across the three runs. Solid arrows denote parent edges selected in at least two runs, whereas dashed arrows denote majority ordering without majority parent-edge selection. The reference pairs are used only for post-fit evaluation.

ODS pair	Order agreement	Edge Jaccard	Pairwise SHD
Poi–NB	$0.882\pm 0.007$	$0.747\pm 0.023$	$16.7\pm 1.5$
Poi–Bin	$0.988\pm 0.008$	$0.924\pm 0.021$	$4.3\pm 1.2$
NB–Bin	$0.880\pm 0.004$	$0.713\pm 0.019$	$19.0\pm 1.0$

Appendix

Table 13: Pairwise sensitivity of ODS to the supplied conditional family on the Lahman q4 representation. Entries are mean

\pm

standard deviation across three paired cross-validation assignments. Pairwise SHD compares two fitted directed graphs, not an estimated graph against a causal ground truth.

Representation	Boundary	$A_{\mathrm{top}}^{\mathcal{R}_{7}}$	$A_{\mathrm{top}}^{\mathcal{R}_{17}}$
q4	clipped	0.905	0.843
q4	admissible	0.762	0.804
q6	admissible	0.762	0.765
q8	admissible	0.762	0.784
raw	admissible	0.000	0.137

Appendix

Table 14: Boundary and support-resolution sensitivity of DISCO on Lahman. The first row reproduces the main q4 analysis; the remaining rows use admissible finite differences with the same DISCO training configuration.

Supplied family	$A_{\mathrm{top}}^{\mathcal{R}_{7}}$	$A_{\mathrm{top}}^{\mathcal{R}_{17}}$	Mean $\|\widehat{E}\|$
Poi	0.333	0.412	73.7
NB	0.857	0.843	69.3
Bin	0.286	0.392	79.7

Appendix

Table 15: ODS sensitivity on the raw Lahman counts. The binomial row uses empirical nodewise support maxima and is not an oracle binomial specification.

Conditional family	Projection	$A_{\mathrm{top}}$	$F_{1}$	SHD	Projection time (s)
Poi	Amortized	$0.931\pm 0.017$	$0.892\pm 0.034$	$28.9\pm 8.7$	92.6
	Stagewise direct	$0.917\pm 0.016$	$0.874\pm 0.034$	$34.0\pm 9.2$	511.9
	Stagewise recursive	$0.916\pm 0.013$	$0.869\pm 0.032$	$35.1\pm 8.2$	544.5
NB	Amortized	$0.931\pm 0.017$	$0.792\pm 0.066$	$51.2\pm 15.4$	94.0
	Stagewise direct	$0.923\pm 0.019$	$0.785\pm 0.055$	$53.3\pm 13.1$	515.0
	Stagewise recursive	$0.924\pm 0.021$	$0.788\pm 0.053$	$52.6\pm 12.9$	545.2

Appendix

Table 16: Comparison of marginal concrete-score projection methods (

d=50

n=5000

). Recovery metrics are means

\pm

sample standard deviations, with seeds paired across methods within each conditional family. Projection time is the mean wall-clock time for the projection stage, in seconds.

	Poi			NB			Bin
Rule	$F_{1}$	Prec.	Rec.	$F_{1}$	Prec.	Rec.	$F_{1}$	Prec.	Rec.
$\tau_{Z}=2,\ \pi_{\min}=1$	0.877	0.934	0.827	0.797	0.917	0.708	0.878	0.943	0.821
$\tau_{Z}=2,\ \pi_{\min}=2/3$	0.891	0.898	0.884	0.870	0.898	0.845	0.890	0.895	0.885
$\tau_{Z}=1,\ \pi_{\min}=1$	0.900	0.898	0.901	0.866	0.892	0.843	0.896	0.906	0.885
$\tau_{Z}=1,\ \pi_{\min}=2/3$	0.800	0.709	0.917	0.813	0.735	0.910	0.809	0.732	0.907

Appendix

Table 17: Parent-selection sensitivity for the threshold settings shown, with fixed learned orders and OCS estimates (

d=50

n=5000

). Entries are arithmetic means over the same seeds. Prec. and Rec. denote precision and recall against the full true DAG, including true edges excluded by the estimated order.

Conditional family	Order	Parent selector	$F_{1}$	SHD
Poi	True	OCS	0.878	30.8
		CAM	0.989	3.2
		PCM-GAM	0.839	54.1
	Estimated	OCS	0.858	37.2
		CAM	0.850	46.3
		PCM-GAM	0.712	105.6

Appendix

Table 18: Parent recovery under fixed orders (

d=50

n=5000

). Entries are arithmetic means over seeds shared by all methods and both order conditions. Metrics use the full true DAG, including edges excluded by an estimated order. Estimated denotes the shared DISCO order; OCS denotes the default DISCO parent selector.

Explore similar work

Sep 17, 2026cs.LG

Epidemiological Causal Graph Identification: Challenges, Identifiability and Algorithms

Causal discovery from observational data is fundamental to statistics and machine learning, yet determining causal direction without interventions necessitates structural assumptions. Existing identifiability research primarily focuses on continuous variables under additive noise models, often neglecting mixed datasets containing ordinal scales, counts, and continuous measurements. This paper investigates causal discovery in Directed Acyclic Graphs (DAGs) where nodes follow either an ordinal distribution (via an ordered logit model) or a regular one-parameter exponential family distribution. We prove that the edge direction between an ordinal and an exponential family node is distributionally identifiable for generic parameter values. Our findings generalize previous Ordinal-Poisson results to the broader exponential family. Computationally, we introduce a score-based exhaustive search and a masked continuous optimization framework using DAGMA for larger graphs. Numerical results validate the theory, recovering edge orientations within a Markov equivalence class that are unidentifiable under classical structural equation models.

Sambit Mishra, Yingying Wang, Christine K. Johnson +1

University of Southern California · University of California, Davis

Causal Discovery Methods Causal Graph Linear Structural Equation Models Directed Acyclic Graph Identifiability Distributions+1

Ask PDF

Apr 28, 2026cs.LG

Optimization-Free Topological Sort for Causal Discovery via the Schur Complement of Score Jacobians

Continuous causal discovery typically couples representation learning with structural optimization via non-convex acyclicity penalties, which subjects solvers to local optima and restricts scalability in high-dimensional regimes. We propose a decoupled paradigm that shifts the causal discovery bottleneck from non-convex optimization to statistical score estimation. We introduce the Score-Schur Topological Sort (SSTS), an algorithm that extracts topological order directly from unconstrained generative models, bypassing constrained structure optimization. We establish that the causal hierarchy leaves a geometric signature within the score function: iterative graph marginalization is mathematically equivalent to computing the Schur complement of the Score-Jacobian Information Matrix (SJIM) under linear conditions. This translates the acyclicity constraint into an algebraic procedure with a dominant cost of O(d^3) operations. For non-linear systems, we formulate the expectation gap of Schur marginalization and introduce Block-SSTS to compress extraction depth, bounding structural error. Empirically, SSTS allows causal structural analysis on non-linear graphs up to d=1000. At this scale, our framework indicates that once the non-convex optimization bottleneck is mathematically bypassed, the structural fidelity of continuous causal discovery is bounded by the finite-sample estimation variance of the global score geometry. By reducing graph extraction to matrix operations, this work reframes scalable causal discovery from a constrained optimization problem to a statistical estimation challenge.

Rui Wu, Hong Xie

School of Computer Science and Engineering University of Science and Technology of China

Causal Discovery Methods Jacobian Fidelity

Ask PDF

Sep 20, 2026cs.LG

Decoupled Causal Discovery

Causal discovery from observational data is a fundamental yet challenging task in scientific research. While existing approaches are primarily based on conditional independence tests, structure scores, or restrictive functional assumptions, we propose Decoupled Causal Discovery (DCD), a novel decoupling-based perspective that does not rely on these methodologies. DCD directly identifies the Markov boundary (MB) by decoupling non-target variables via weighting functions, such that only variables within the MB preserve dependence with the target under the decoupled distribution. Building on this, DCD iteratively constructs the Completed Partially Directed Acyclic Graph (CPDAG) by exploiting structural asymmetries within the MBs. We establish the theoretical identifiability, soundness, and completeness of DCD. Empirical evaluations demonstrate that DCD achieves strong performance, particularly excelling in challenging noise regimes.

Zhengkang Guan, Fei Wu, Kun Kuang

College of Computer Science and Technology Zhejiang University

Causal Discovery Methods Decoupling Markov

Ask PDF