Causal Discovery in Equal Variance Linear Gaussian DAGs via SURE-Tuned Ridge Regression

# Causal Discovery in Equal Variance Linear Gaussian DAGs via SURE-Tuned Ridge Regression
Source: [https://arxiv.org/html/2608.17132](https://arxiv.org/html/2608.17132)
## Causal Discovery in Equal Variance Linear Gaussian DAGs via SURE\-Tuned Ridge RegressionThis work has been funded by one or all of the following grants: ARO W911NF1910269, ARO W911NF2410094, ONR N00014\-22\-1\-2363, NSF CIF\-2311653, NSF CIF\-2148313, NSF RINGS\-2148313, NSF DBI\-2412522, and is also supported in part by funds from federal agency and industry partners as specified in the RINGS program\.

Sambit Mishra & Urbashi MitraAffiliation:Department of Electrical and Computer Engineering University of Southern California \{sambitmi, ubli\}@usc\.edu

###### Abstract

Recovering the directed acyclic graph \(DAG\) of a structural equation model \(SEM\) from observational data is a central problem in causal discovery\. The iterative gradient descent and per\-problem hyperparameter tuning of continuous\-optimization methods are poorly suited to two practically important regimes: the sample\-limited regime, where the number of samples is comparable to or smaller than the number of nodes in the DAG, and the compute\-limited regime\. This work proposes SURE\-Ridge, a non\-iterative, closed\-form estimator for equal variance linear Gaussian SEM\. The method performs parallel node\-wise regressions with regularization parameters chosen adaptively by Stein’s unbiased risk estimate \(SURE\), and applies an adaptive thresholding procedure to extract a DAG from the resulting soft adjacency matrix\. Numerical results show that SURE\-Ridge achieves the lowest structural Hamming distance in the small\-sample regime and the lowest run time across all sample sizes tested, compared with NOTEARS, DAGMA, and GBNSL baselines\.

###### Index Terms:

Causal Discovery, Structural Equation Models, Directed Acyclic Graphs, Stein’s unbiased risk estimate, ridge regression

## IIntroduction

Linear structural equation models \(SEMs\) provide a flexible framework for representing causal relationships among random variables, and recovering the underlying directed acyclic graph \(DAG\) from the observational data is a central problem in causal discovery, with applications spanning epidemiology\[[12](https://arxiv.org/html/2608.17132#bib.bib7)\], wireless networks\[[9](https://arxiv.org/html/2608.17132#bib.bib8)\], and financial fraud detection\[[16](https://arxiv.org/html/2608.17132#bib.bib9)\]\. For linear Gaussian SEMs with equal noise variance across all nodes, the true DAG is identifiable from the joint distribution alone\[[11](https://arxiv.org/html/2608.17132#bib.bib1)\]\. We consider two key regimes herein: \(a\) sample\-limited regimes \(as considered in\[[13](https://arxiv.org/html/2608.17132#bib.bib3),[8](https://arxiv.org/html/2608.17132#bib.bib4)\]\), where the number of samples is often comparable to or even smaller than the number of nodes, as in\[[10](https://arxiv.org/html/2608.17132#bib.bib13),[19](https://arxiv.org/html/2608.17132#bib.bib14)\], and \(b\) compute\-limited regimes where a DAG must be continuously estimated as in online and streaming applications\[[17](https://arxiv.org/html/2608.17132#bib.bib15)\]\.

Classical approaches such as the PC\[[14](https://arxiv.org/html/2608.17132#bib.bib10)\]and GES\[[5](https://arxiv.org/html/2608.17132#bib.bib11)\]have computational costs that scale poorly with the number of nodes in the DAG\. Continuous optimization formulations introduced by NOTEARS\[[20](https://arxiv.org/html/2608.17132#bib.bib6)\]and DAGMA\[[2](https://arxiv.org/html/2608.17132#bib.bib12)\]recast structure learning as a smooth, constrained optimization problem via a differentiable acyclicity functional, but they require thousands of gradient steps and per\-problem hyperparameter tuning, and are thus not well\-matched to the two regimes under consideration\.

Work that directly exploits the equal\-variance assumption directly considers two principles\. The first is iterative peeling from the inverse covariance matrix: GBNSL in\[[6](https://arxiv.org/html/2608.17132#bib.bib16)\]identifies terminal nodes one at a time, and BUILD\[[1](https://arxiv.org/html/2608.17132#bib.bib17)\]prunes leaves bottom\-up\. The second approach is staged regression, in which ordering and parent selection are separated into sequential phases: first ranking nodes by conditional variances and then fitting parents on the recovered order\[[4](https://arxiv.org/html/2608.17132#bib.bib18)\]\. Both strategies treat ordering recovery as an explicit intermediate step, in contrast to SURE\-Ridge, and both compound finite\-sample errors across stages when the sample size is comparable to the number of nodes in the DAG\.

We proposeSURE\-Ridge, a non\-iterative, closed\-form estimator for equal\-variance linear Gaussian DAG recovery\. The method performs parallel node\-wise ridge regressions, with theL2L\_\{2\}regularization parameter chosen by minimizing Stein’s unbiased risk estimate \(SURE\)\[[15](https://arxiv.org/html/2608.17132#bib.bib2)\]\. SURE provides a tractable estimate of the prediction error yielding a tuning rule that requires neither held\-out data nor iterative search, thus providing a decisive advantage over cross\-validation when samples are scarce and over gradient\-based tuning when many DAGs must be estimated\. Pairing SURE with ridge regression is natural under the equal\-variance assumption, where each node’s structural equation is itself a linear Gaussian regression on its parents, and subproblems decouple, enabling parallel discovery\. The contributions of this paper are as follows:

1. 1\.A closed\-form formulation of equal variance linear Gaussian DAG recovery via SURE\-tuned node\-wise regression with total computational cost𝒪⁡\(d4\+d3​n\+d​M​\(d−1\)\+d3​log⁡d\)\\mathcal\{O\}\\left\(d^\{4\}\+d^\{3\}n\+dM\\left\(d\-1\\right\)\+d^\{3\}\\log d\\right\)per DAG;
2. 2\.An adaptive thresholding procedure that requires only two scale\-invariant constants and produces a valid DAG in the typical case;
3. 3\.Numerical results demonstrating that SURE\-Ridge achieves the lowest structural Hamming distance \(SHD\) in the small\-sample regime and the lowest run time across all sample sizes tested, against NOTEARS, DAGMA, and GBNSL baselines\.

This paper is organized as follows\. Section[II](https://arxiv.org/html/2608.17132#S2)describes the system model and notation\. Section[III](https://arxiv.org/html/2608.17132#S3)derives the SURE\-tuned ridge estimator\. Section[IV](https://arxiv.org/html/2608.17132#S4)presents the adaptive thresholding procedure\. Section[V](https://arxiv.org/html/2608.17132#S5)reports numerical results and Section[VI](https://arxiv.org/html/2608.17132#S6)concludes the work\.

## IISystem Model

We consider ad−d\-node causal DAG𝒢=\(𝒱,ℰ\)\\mathcal\{G\}=\\left\(\\mathcal\{V\},\\mathcal\{E\}\\right\), where𝒱=\{1,…,d\}\\mathcal\{V\}=\\left\\\{1,\\dots,d\\right\\\}represents the node set andℰ=\{\(i,j\):edge​i→j​exists\}\\mathcal\{E\}=\\left\\\{\\left\(i,j\\right\):\\text\{edge \}i\\rightarrow j\\text\{ exists\}\\right\\\}represents the edge set\. We represent the graph using a weighted adjacency matrix𝐖∈ℝd×d\\mathbf\{W\}\\in\\mathbb\{R\}^\{d\\times d\}under the convention

𝐖i,j≠0⇔\(j,i\)∈ℰ\.\\mathbf\{W\}\_\{i,j\}\\neq 0\\iff\\left\(j,i\\right\)\\in\\mathcal\{E\}\.\(1\)For each nodei∈𝒱i\\in\\mathcal\{V\}, we define the row\-vector of causal coefficients,𝐰i∈ℝ1×\(d−1\)\\mathbf\{w\}\_\{i\}\\in\\mathbb\{R\}^\{1\\times\\left\(d\-1\\right\)\}, \(it​hi^\{th\}row of𝐖\\mathbf\{W\}with𝐖i,i\\mathbf\{W\}\_\{i,i\}removed\),

𝐰i≜𝐖i,−i\.\\mathbf\{w\}\_\{i\}\\triangleq\\mathbf\{W\}\_\{i,\-i\}\.\(2\)LetXiX\_\{i\}be the random variable associated with nodeiiand𝐱i∈ℝ1×n\\mathbf\{x\}\_\{i\}\\in\\mathbb\{R\}^\{1\\times n\}be the row\-vector ofnni\.i\.d\. observations ofXiX\_\{i\}\. The concatenation of observations is denoted𝐗∈ℝd×n\\mathbf\{X\}\\in\\mathbb\{R\}^\{d\\times n\}, where𝐗i,:=𝐱i\\mathbf\{X\}\_\{i,:\}=\\mathbf\{x\}\_\{i\}\. Under a linear Gaussian SEM with equal noise variance,

𝐱i=𝐰i𝐗−i,:\+𝐧i,\\mathbf\{x\}\_\{i\}=\\mathbf\{w\}\_\{i\}\\mathbf\{X\}\_\{\-i,:\}\+\\mathbf\{n\}\_\{i\},\(3\)where𝐗−i,:\\mathbf\{X\}\_\{\-i,:\}represents the data matrix𝐗\\mathbf\{X\}with theit​hi^\{th\}row removed, and𝐧i∼𝒩⁡\(𝟎,σ2​𝐈n\)\\mathbf\{n\}\_\{i\}\\sim\\mathcal\{N\}\\left\(\\mathbf\{0\},\\sigma^\{2\}\\mathbf\{I\}\_\{n\}\\right\)and i\.i\.d\. with respect toii\. This case is identifiable as per\[[11](https://arxiv.org/html/2608.17132#bib.bib1)\]\.

We see that \([3](https://arxiv.org/html/2608.17132#S2.E3)\) can be interpreted as a regression to determineXiX\_\{i\}from noisy observations governed by the causal coefficient vector𝐰i\\mathbf\{w\}\_\{i\}\. Thus, we can estimate𝐰i\\mathbf\{w\}\_\{i\}using ridge regression, with the tunable regularization parameterλi\\lambda\_\{i\}\. Let the Gram matrix corresponding to nodeiibe denoted as𝐆i=𝐗−i,:𝐗−i,:T\\mathbf\{G\}\_\{i\}=\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{X\}\_\{\-i,:\}^\{T\}, then we have the ridge regression\-based estimator of𝐰i\\mathbf\{w\}\_\{i\}given by

𝐰^i\(λi\)=\(𝐆i\+λi𝐈d−1\)−1𝐗−i,:𝐱iT\.\\hat\{\\mathbf\{w\}\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)=\\left\(\\mathbf\{G\}\_\{i\}\+\\lambda\_\{i\}\\mathbf\{I\}\_\{d\-1\}\\right\)^\{\-1\}\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{x\}\_\{i\}^\{T\}\.\(4\)A principled choice of the tunable parameterλi\\lambda\_\{i\}is critical, as it directly controls the bias\-variance trade\-off of the resulting estimator𝐰^i​\(λi\)\\hat\{\\mathbf\{w\}\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)\. We adopt SURE as a data\-driven criterion for selectingλi\\lambda\_\{i\}\. The following proposition formalizes the per\-node SURE objective\.

###### Proposition 1\(Per\-node SURE Objective\)\.

Let the noise vector satisfy𝐧iT∼𝒩⁡\(𝟎,σ2​𝐈n\)\\mathbf\{n\}\_\{i\}^\{T\}\\sim\\mathcal\{N\}\\left\(\\mathbf\{0\},\\sigma^\{2\}\\mathbf\{I\}\_\{n\}\\right\), and let𝐰^i​\(λi\)\\hat\{\\mathbf\{w\}\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)be the ridge estimator defined in \([4](https://arxiv.org/html/2608.17132#S3.E4)\)\. Then an unbiased estimator of the expected squared\-error risk𝔼\[∥𝐱iT−𝐗−i,:T𝐰^iT\(λi\)∥22\]−nσ2\\mathbb\{E\}\\left\[\\\|\\mathbf\{x\}\_\{i\}^\{T\}\-\\mathbf\{X\}\_\{\-i,:\}^\{T\}\\hat\{\\mathbf\{w\}\}\_\{i\}^\{T\}\\left\(\\lambda\_\{i\}\\right\)\\\|\_\{2\}^\{2\}\\right\]\-n\\sigma^\{2\}is given by

SUREi⁡\(λi\)\\displaystyle\\operatorname\{SURE\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)=\\displaystyle=−n​σ2\\displaystyle\-n\\sigma^\{2\}\(5\)\+∥𝐱iT−𝐗−i,:T\(𝐆i\+λi𝐈d−1\)−1𝐗−i,:𝐱iT∥22\\displaystyle\\\!\\\!\\\!\\\!\\\!\\\!\+\\\|\\mathbf\{x\}\_\{i\}^\{T\}\-\\mathbf\{X\}\_\{\-i,:\}^\{T\}\\left\(\\mathbf\{G\}\_\{i\}\+\\lambda\_\{i\}\\mathbf\{I\}\_\{d\-1\}\\right\)^\{\-1\}\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{x\}\_\{i\}^\{T\}\\\|\_\{2\}^\{2\}\+2​σ2​tr⁡\[\(𝐆i\+λi​𝐈d−1\)−1​𝐆i\]\.\\displaystyle\\\!\\\!\\\!\\\!\\\!\\\!\+2\\sigma^\{2\}\\operatorname\{tr\}\\left\[\\left\(\\mathbf\{G\}\_\{i\}\+\\lambda\_\{i\}\\mathbf\{I\}\_\{d\-1\}\\right\)^\{\-1\}\\mathbf\{G\}\_\{i\}\\right\]\.

The proposition follows from applying the SURE construction in\[[15](https://arxiv.org/html/2608.17132#bib.bib2)\]to the regression model in \([3](https://arxiv.org/html/2608.17132#S2.E3)\)\. The objective in \([5](https://arxiv.org/html/2608.17132#S3.E5)\) thus gives us a principled risk metric for choosing the optimal parameterλi\\lambda\_\{i\}for each nodei∈𝒱i\\in\\mathcal\{V\}, which we obtain by solving the following optimization problem

λi∗=arg⁡minλi\>0​SUREi⁡\(λi\)\.\\lambda\_\{i\}^\{\*\}=\\arg\\min\_\{\\lambda\_\{i\}\>0\}\\operatorname\{SURE\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)\.\(6\)
We further massage the objective function via the eigendecomposition𝐆i=𝐔i​𝚪i​𝐔iT\\mathbf\{G\}\_\{i\}=\\mathbf\{U\}\_\{i\}\\boldsymbol\{\\Gamma\}\_\{i\}\\mathbf\{U\}\_\{i\}^\{T\}, where𝐔i\\mathbf\{U\}\_\{i\}is the orthogonal matrix of eigenvectors of𝐆i\\mathbf\{G\}\_\{i\}and𝚪i=diag⁡\(γi,1,…,γi,d−1\)\\boldsymbol\{\\Gamma\}\_\{i\}=\\operatorname\{diag\}\\left\(\\gamma\_\{i,1\},\\dots,\\gamma\_\{i,d\-1\}\\right\)collects the eigenvalues\. Defining the projection𝐳i=𝐔iT𝐗−i,:𝐱iT\\mathbf\{z\}\_\{i\}=\\mathbf\{U\}\_\{i\}^\{T\}\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{x\}\_\{i\}^\{T\}translates the SURE objective to a sum of scalar functions ofλi\\lambda\_\{i\}that can be evaluated at𝒪⁡\(d−1\)\\mathcal\{O\}\\left\(d\-1\\right\)cost per grid point\[[7](https://arxiv.org/html/2608.17132#bib.bib19)\]\. The SURE objective in \([5](https://arxiv.org/html/2608.17132#S3.E5)\) admits the following closed\-form expression

SUREi⁡\(λi\)\\displaystyle\\operatorname\{SURE\}\_\{i\}\\left\(\\lambda\_\{i\}\\right\)\\\!\\\!\\\!=\\displaystyle=−n​σ2\+‖𝐱i‖22−2​∑k=1d−1zi,k2γi,k\+λi\\displaystyle\\\!\\\!\\\!\-n\\sigma^\{2\}\+\\\|\\mathbf\{x\}\_\{i\}\\\|\_\{2\}^\{2\}\-2\\sum\_\{k=1\}^\{d\-1\}\\frac\{z\_\{i,k\}^\{2\}\}\{\\gamma\_\{i,k\}\+\\lambda\_\{i\}\}\(7\)\+∑k=1d−1γi,k​zi,k2\(γi,k\+λi\)2\+2σ2∑k=1d−1γi,kγi,k\+λi\.\\displaystyle\\\!\\\!\\\!\\\!\\\!\\\!\\\!\\\!\\\!\\\!\\\!\\\!\+\\sum\_\{k=1\}^\{d\-1\}\\frac\{\\gamma\_\{i,k\}\\,z\_\{i,k\}^\{2\}\}\{\\left\(\\gamma\_\{i,k\}\+\\lambda\_\{i\}\\right\)^\{2\}\}\+2\\sigma^\{2\}\\sum\_\{k=1\}^\{d\-1\}\\frac\{\\gamma\_\{i,k\}\}\{\\gamma\_\{i,k\}\+\\lambda\_\{i\}\}\.Given the optimal parameterλi∗\\lambda\_\{i\}^\{\*\}, the corresponding closed form ridge estimate is

𝐰^i∗=𝐳iT​diag⁡\(1γi,k\+λi∗\)k=1d−1​𝐔iT\.\\displaystyle\\hat\{\\mathbf\{w\}\}\_\{i\}^\{\*\}=\\mathbf\{z\}\_\{i\}^\{T\}\\operatorname\{diag\}\\left\(\\frac\{1\}\{\\gamma\_\{i,k\}\+\\lambda\_\{i\}^\{\*\}\}\\right\)\_\{k=1\}^\{d\-1\}\\mathbf\{U\}\_\{i\}^\{T\}\.\(8\)Finally, concatenating all the𝐰^i∗\\hat\{\\mathbf\{w\}\}\_\{i\}^\{\*\}row\-wise with additional diagonal element set to00leads us to the final soft estimate of the adjacency matrix𝐖~est\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\.

The SURE\-tuned ridge estimator captures the population\-level regression structure that ensures identifiability under the noise\-equal\-variance assumption\[[11](https://arxiv.org/html/2608.17132#bib.bib1)\], and the adaptive thresholding stage of Section[IV](https://arxiv.org/html/2608.17132#S4)extracts a DAG estimate from the resulting soft adjacency matrix\.

## IVAdaptive Thresholding

The soft adjacency matrix𝐖~est\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}is typically not a DAG\. To extract a valid DAG structure, we propose a three\-stage procedure\. We start with an absolute floor threshold that scales with the signal magnitude, then perform a discrete bisection search over the above\-floor magnitudes, and lastly consider a fallback that preserves informative edges when the bisection becomes too aggressive\.

### IV\-AAcyclicity Test

We use the polynomial acyclicity functional of\[[18](https://arxiv.org/html/2608.17132#bib.bib5)\],

h⁡\(𝐖\)=tr⁡\[\(𝐈d\+𝐖⊙𝐖d\)d\]−d,\\displaystyle h\\left\(\\mathbf\{W\}\\right\)=\\operatorname\{tr\}\\left\[\\left\(\\mathbf\{I\}\_\{d\}\+\\frac\{\\mathbf\{W\}\\odot\\mathbf\{W\}\}\{d\}\\right\)^\{d\}\\right\]\-d,\(9\)which is a numerically stable variant of the matrix exponential form originally proposed in\[[20](https://arxiv.org/html/2608.17132#bib.bib6)\]\. It can be verified thath⁡\(𝐖\)=0h\\left\(\\mathbf\{W\}\\right\)=0if and only if𝐖\\mathbf\{W\}represents a DAG\. Evaluatinghhat a candidate thresholded graph is𝒪⁡\(d3\)\\mathcal\{O\}\\left\(d^\{3\}\\right\)in computational complexity and vectorizes trivially across a batch of candidate matrices\.

### IV\-BSignal\-Scaled Floor Filter

A simple threshold of zero retains all the ridge\-shrinkage noise, while a bisection starting from zero can, in dense regimes, push the threshold above the true edge magnitudes and return an empty graph\. We therefore introduce a scale\-adaptive floor

ηmin=max⁡\(η0,β​maxi≠j​\|\(𝐖~est\)i,j\|\),\\displaystyle\\eta\_\{\\operatorname\{min\}\}=\\max\\left\(\\eta\_\{0\},\\beta\\max\_\{i\\neq j\}\\left\|\\left\(\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\\right\)\_\{i,j\}\\right\|\\right\),\(10\)whereη0\\eta\_\{0\}andβ\\betaare tunable hyperparameters governed by the number of nodesddand the expected signal\-to\-noise ratio \(SNR\)\. Since the ridge shrinkage is bounded above by OLS,max⁡\|𝐖~est\|\\max\\left\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\\right\|tracks the largest signal magnitude\. In small\-to\-moderate\-ddor high\-SNR settings, theβ\\beta\-scaled term is highly significant and acts as a proper scale\-adaptive floor, and therefore it is set high in such settings to aggressively filter out spurious correlations\. However, in large\-ddor low\-SNR regimes, heavy ridge regularization severely suppresses the maximum edge weights toward the noise floor, rendering theβ\\betaterm insufficient\. Consequently,η0\\eta\_\{0\}acts as an absolute safety bound and is tuned as per the expected noise parameters to separate ambient noise artifacts from the heavily penalized true values\.

The first stage then zeroes out all entries belowηmin\\eta\_\{\\operatorname\{min\}\}\. Formally, we have

𝐖floor=𝕀⁡\(\|𝐖~est\|\>ηmin\)⊙𝐖~est\.\\displaystyle\\mathbf\{W\}\_\{\\operatorname\{floor\}\}=\\mathbb\{I\}\\left\(\\left\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\\right\|\>\\eta\_\{\\operatorname\{min\}\}\\right\)\\odot\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\.\(11\)If we obtainh⁡\(𝐖floor\)=0h\\left\(\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\\right\)=0, in this thresholding step, the procedure stops with𝐖floor\\mathbf\{W\}\_\{\\operatorname\{floor\}\}as the final estimate\.

### IV\-CBisection Over Above\-Floor Magnitudes

When𝐖floor\\mathbf\{W\}\_\{\\operatorname\{floor\}\}is still cyclic, we search for a larger threshold that restores acyclicity\. The candidate set is\{ηmin\}∪\{\|𝐖~est,i​j\|:\|𝐖~est,i​j\|\>ηmin\}\\left\\\{\\eta\_\{\\min\}\\right\\\}\\cup\\left\\\{\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\},ij\}\|:\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\},ij\}\|\>\\eta\_\{\\min\}\\right\\\}, sorted in ascending order\. Because the mapη↦𝕀⁡\[𝕀⁡\(\|𝐖~est\|\>η\)⊙𝐖~est​is a DAG\]\\eta\\mapsto\\mathbb\{I\}\\left\[\\mathbb\{I\}\\left\(\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\|\>\\eta\\right\)\\odot\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\\text\{ is a DAG\}\\right\]is upward\-closed on this sorted list, the smallest acyclic candidateη∗\\eta^\{\*\}can be located by discrete binary search in at most⌈2​log2​d⌉\\lceil 2\\log\_\{2\}d\\rceilacyclicity checks\. We denote the resulting estimate𝐖bisect\\mathbf\{W\}\_\{\\operatorname\{bisect\}\}\.

### IV\-DFallback for Aggressive Cuts

Aggressive bisection can remove most floor\-thresholded edges to break a small number of cycles, yielding a near\-empty DAG that is less informative than a mildly cyclic matrix\. To guard against this, we return

𝐖est=\{𝐖floorif​‖𝐖bisect‖0<12​‖𝐖floor‖0,𝐖bisectotherwise,\\mathbf\{W\}\_\{\\operatorname\{est\}\}=\\begin\{cases\}\\mathbf\{W\}\_\{\\operatorname\{floor\}\}&\\text\{if \}\\\|\\mathbf\{W\}\_\{\\operatorname\{bisect\}\}\\\|\_\{0\}<\\tfrac\{1\}\{2\}\\\|\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\\\|\_\{0\},\\\\ \\mathbf\{W\}\_\{\\operatorname\{bisect\}\}&\\text\{otherwise,\}\\end\{cases\}\(12\)where∥⋅∥0\\\|\\cdot\\\|\_\{0\}counts the number of nonzero entries\. The fallback prioritizes edge preservation over strict acyclicity\.𝐖est\\mathbf\{W\}\_\{\\operatorname\{est\}\}is a valid DAG whenever the floor or bisection stage succeeds which commonly occurred in our numerical results\. In the rare fallback case𝐖est\\mathbf\{W\}\_\{\\operatorname\{est\}\}may contain a small number of cycles\.

The entire algorithm is summarized in Algorithm[1](https://arxiv.org/html/2608.17132#alg1)\. The dominant per\-node cost in the SURE\-tuned ridge stage is the eigendecomposition of𝐆i\\mathbf\{G\}\_\{i\}, which is𝒪⁡\(\(d−1\)3\)\\mathcal\{O\}\\left\(\\left\(d\-1\\right\)^\{3\}\\right\), followed by the projection𝐳i=𝐔iT𝐗−i,:𝐱iT\\mathbf\{z\}\_\{i\}=\\mathbf\{U\}\_\{i\}^\{T\}\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{x\}\_\{i\}^\{T\}at𝒪⁡\(\(d−1\)2​n\)\\mathcal\{O\}\\left\(\\left\(d\-1\\right\)^\{2\}n\\right\)and the SURE grid evaluation at𝒪⁡\(M⁡\(d−1\)\)\\mathcal\{O\}\\left\(M\\left\(d\-1\\right\)\\right\)\. Since these steps are repeated for each of theddnodes, the total cost of the SURE\-tuned ridge stage is𝒪⁡\(d​\(d−1\)3\+d​\(d−1\)2​n\+d​M​\(d−1\)\)\\mathcal\{O\}\\left\(d\\left\(d\-1\\right\)^\{3\}\+d\\left\(d\-1\\right\)^\{2\}n\+dM\\left\(d\-1\\right\)\\right\)\. The thresholding stage adds a single𝒪⁡\(d3\)\\mathcal\{O\}\\left\(d^\{3\}\\right\)acyclicity evaluation in the floor filter, and at most⌈2​log2​d⌉\\lceil 2\\log\_\{2\}d\\rceilfurther𝒪⁡\(d3\)\\mathcal\{O\}\\left\(d^\{3\}\\right\)evaluations during the bisection\. Therefore, the worst\-case overall complexity of Algorithm[1](https://arxiv.org/html/2608.17132#alg1)is𝒪⁡\(d4\+d3​n\+d​M​\(d−1\)\+d3​log⁡d\)\\mathcal\{O\}\\left\(d^\{4\}\+d^\{3\}n\+dM\\left\(d\-1\\right\)\+d^\{3\}\\log d\\right\)per DAG\.

Algorithm 1SURE\-Ridge DAG Recovery1:Data matrix

𝐗∈ℝd×n\\mathbf\{X\}\\in\\mathbb\{R\}^\{d\\times n\}, noise variance

σ2\\sigma^\{2\}, grid

Λ\\Lambda, floor

η0\\eta\_\{0\}, scale

β\\beta\.

2:Estimated adjacency matrix

𝐖est∈ℝd×d\\mathbf\{W\}\_\{\\operatorname\{est\}\}\\in\\mathbb\{R\}^\{d\\times d\}\.

3:for

i=1,…,di=1,\\ldots,ddo

4:Form

𝐗−i,:\\mathbf\{X\}\_\{\-i,:\}and

𝐆i=𝐗−i,:𝐗−i,:T\\mathbf\{G\}\_\{i\}=\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{X\}\_\{\-i,:\}^\{T\}; eigendecompose

𝐆i=𝐔i​𝚪i​𝐔iT\\mathbf\{G\}\_\{i\}=\\mathbf\{U\}\_\{i\}\\boldsymbol\{\\Gamma\}\_\{i\}\\mathbf\{U\}\_\{i\}^\{T\}\.

5:Compute

𝐳i=𝐔iT𝐗−i,:𝐱iT\\mathbf\{z\}\_\{i\}=\\mathbf\{U\}\_\{i\}^\{T\}\\mathbf\{X\}\_\{\-i,:\}\\mathbf\{x\}\_\{i\}^\{T\}and

SUREi⁡\(λ\)\\operatorname\{SURE\}\_\{i\}\\left\(\\lambda\\right\)for each

λ∈Λ\\lambda\\in\\Lambda\.

6:

λi∗←arg⁡minλ∈Λ​SUREi⁡\(λ\)\\lambda\_\{i\}^\{\*\}\\leftarrow\\arg\\min\_\{\\lambda\\in\\Lambda\}\\operatorname\{SURE\}\_\{i\}\\left\(\\lambda\\right\); recover

\(𝐰^iT\)∗\\left\(\\hat\{\\mathbf\{w\}\}\_\{i\}^\{T\}\\right\)^\{\*\}using the eigendecomposition\.

7:endfor

8:Assemble

𝐖~est\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}from

\{\(𝐰^iT\)∗\}\\left\\\{\\left\(\\hat\{\\mathbf\{w\}\}\_\{i\}^\{T\}\\right\)^\{\*\}\\right\\\}with zero diagonal\.

9:Compute

ηmin=max⁡\{η0,β​maxi≠j​\|𝐖~est,i​j\|\}\\eta\_\{\\min\}=\\max\\left\\\{\\eta\_\{0\},\\ \\beta\\max\_\{i\\neq j\}\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\},ij\}\|\\right\\\}and set

𝐖floor=𝕀⁡\(\|𝐖~est\|\>ηmin\)⊙𝐖~est\\mathbf\{W\}\_\{\\operatorname\{floor\}\}=\\mathbb\{I\}\\left\(\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\|\>\\eta\_\{\\min\}\\right\)\\odot\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\.

10:if

h⁡\(𝐖floor\)=0h\\left\(\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\\right\)=0then

11:return

𝐖floor\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\.

12:endif

13:Binary\-search the sorted above\-floor magnitudes for the smallest

η∗\\eta^\{\*\}giving an acyclic

𝐖bisect=𝕀⁡\(\|𝐖~est\|\>η∗\)⊙𝐖~est\\mathbf\{W\}\_\{\\operatorname\{bisect\}\}=\\mathbb\{I\}\\left\(\|\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\|\>\\eta^\{\*\}\\right\)\\odot\\tilde\{\\mathbf\{W\}\}\_\{\\operatorname\{est\}\}\.

14:if

‖𝐖bisect‖0<12​‖𝐖floor‖0\\\|\\mathbf\{W\}\_\{\\operatorname\{bisect\}\}\\\|\_\{0\}<\\tfrac\{1\}\{2\}\\\|\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\\\|\_\{0\}then

15:return

𝐖floor\\mathbf\{W\}\_\{\\operatorname\{floor\}\}\.

16:else

17:return

𝐖bisect\\mathbf\{W\}\_\{\\operatorname\{bisect\}\}\.

18:endif

## VNumerical Results and Discussion

### V\-AExperimental Setup

We simulate Gaussian SEMs withd∈\{20,50\}d\\in\\left\\\{20,50\\right\\\}nodes and⌈2​d/5⌉\\lceil 2d/5\\rceiledges per graph, drawn from a sparse random\-DAG distribution\. This corresponds to a sub\-Erdős–Rényi\-1 average in\-degree, matching the sparse regime targeted by equal\-variance methods\. Non\-zero entries of𝐖\\mathbf\{W\}are sampled uniformly from\[−2\.0,−1\.0\]∪\[1\.0,2\.0\]\\left\[\-2\.0,\-1\.0\\right\]\\cup\\left\[1\.0,2\.0\\right\]\. The noise variance is fixed atσ2=1\\sigma^\{2\}=1\. We have100100Monte Carlo trials per setting \(values ofnnanddd\)\. Each method observes the same\(𝐖true,𝐗\)\\left\(\\mathbf\{W\}\_\{\\operatorname\{true\}\},\\mathbf\{X\}\\right\)per trial to enable fair comparisons, where𝐖true\\mathbf\{W\}\_\{\\operatorname\{true\}\}is the ground truth DAG adjacency matrix\.

We compare SURE\-Ridge against three baselines: NOTEARS, DAGMA, and GBNSL\. NOTEARS\[[20](https://arxiv.org/html/2608.17132#bib.bib6)\]solves

min𝐖∈ℝd×d⁡12​n​‖𝐗−𝐖𝐗‖F2\+λ​‖𝐖‖1s\.t\.h⁡\(𝐖\)=0,\\displaystyle\\min\_\{\\mathbf\{W\}\\in\\mathbb\{R\}^\{d\\times d\}\}\\;\\tfrac\{1\}\{2n\}\\\|\\mathbf\{X\}\-\\mathbf\{WX\}\\\|\_\{F\}^\{2\}\+\\lambda\\\|\\mathbf\{W\}\\\|\_\{1\}\\quad\\text\{s\.t\.\}\\quad h\\left\(\\mathbf\{W\}\\right\)=0,\(13\)with the trace\-exponential acyclicity functionalh⁡\(𝐖\)=tr⁡\(e𝐖⊙𝐖\)−dh\\left\(\\mathbf\{W\}\\right\)=\\operatorname\{tr\}\\left\(e^\{\\mathbf\{W\}\\odot\\mathbf\{W\}\}\\right\)\-d, via an augmented Lagrangian\. We substitute the polynomial form in \([9](https://arxiv.org/html/2608.17132#S4.E9)\) for numerical stability\. DAGMA\[[2](https://arxiv.org/html/2608.17132#bib.bib12)\]swaps this constraint for a log\-determinant form,

hldet\(𝐖;s\)=−logdet\(s𝐈d−𝐖⊙𝐖\)\+dlogs,\\displaystyle h\_\{\\text\{ldet\}\}\\left\(\\mathbf\{W\};s\\right\)=\-\\log\\det\\left\(s\\mathbf\{I\}\_\{d\}\-\\mathbf\{W\}\\odot\\mathbf\{W\}\\right\)\+d\\log s,\(14\)which has better\-behaved gradients, and replaces the augmented Lagrangian with a central\-path scheme that solves a sequence of unconstrained subproblems\.

GBNSL\[[6](https://arxiv.org/html/2608.17132#bib.bib16)\]uses successive node pruning to determine the graph and does not require gradient\-based optimization\. Given a sparseℓ1\\ell\_\{1\}\-regularized estimate𝛀^\\hat\{\\boldsymbol\{\\Omega\}\}of the precision matrix and least\-squares estimates𝜽^i\\hat\{\\boldsymbol\{\\theta\}\}^\{i\}of the regression of eachXiX\_\{i\}on its Markov blanket \(parents and children\), GBNSL uses the identity that under equal noise varianceσ2\\sigma^\{2\}, nodeiiis a leaf node if and only if

θji=σ2​Ωi,j,∀j≠i\.\\displaystyle\\theta^\{i\}\_\{j\}=\\sigma^\{2\}\\,\\Omega\_\{i,j\},\\quad\\forall j\\neq i\.\(15\)At each iteration, leaf nodes are determined as those minimizing\|Ω^i,j/θ^ji\|\|\\hat\{\\Omega\}\_\{i,j\}/\\hat\{\\theta\}^\{i\}\_\{j\}\|over its Markov blanket, after which𝛀^\\hat\{\\boldsymbol\{\\Omega\}\}is updated by a rank\-1 Schur complement\. Conditioned on the found leaf nodes enables the next iteration which continues until the causal ordering is determined\. A final least\-squares estimation yields the graph edge weights\. NOTEARS and DAGMA are initialized to zero with a matched budget of3200032000Adam steps \(learning rate3×10−43\\times 10^\{\-4\},L1L\_\{1\}penalty5×10−25\\times 10^\{\-2\}\): NOTEARS uses3232outer iterations of10001000inner steps, DAGMA usesT=4T=4central\-path stages of80008000inner steps withμ0=1\.0\\mu\_\{0\}=1\.0and decayα=0\.1\\alpha=0\.1\. GBNSL’sℓ1\\ell\_\{1\}parameter and Markov\-blanket threshold are tuned once and held fixed\. All three baselines apply a post\-hoc threshold of0\.10\.1before scoring\. We measure accuracy by normalized SHD,

Normalized SHD=2​SHDd⁡\(d−1\),\\displaystyle\\text\{Normalized SHD\}=\\frac\{2\\,\\operatorname\{SHD\}\}\{d\\left\(d\-1\\right\)\},\(16\)where SHD is the minimum number of edge insertions, deletions, or reversals needed to transform one graph into another\. The code for the implementation can be found at[https://github\.com/SamMathelete/SURE\_Ridge\_Final](https://github.com/SamMathelete/SURE_Ridge_Final)\.

### V\-BSample Efficiency

Figure[1](https://arxiv.org/html/2608.17132#S5.F1)shows the normalized SHD as a function ofnnatd=20d=20andd=50d=50\. The four methods exhibit qualitatively distinct sample\-efficiency profiles, and the ranking is regime\-dependent\. In the sample\-limited regime \(nnon the order ofddor smaller\), SURE\-Ridge attains the lowest normalized SHD by a substantial margin, with the gap widening asddgrows\. This indicates that the advantage is not a smallddartifact\. The two iterative methods, NOTEARS and DAGMA, perform considerably worse here, and GBNSL performs the worst of the four at very smallnn, with normalized SHD close to11ford=50d=50\. This is consistent with the high\-dimensional regime, where the guarantees of GBNSL’s CLIME\-based\[[3](https://arxiv.org/html/2608.17132#bib.bib20)\]precision estimate do not yet apply\. Asnngrows, the ordering changes: GBNSL improves the fastest of any method and overtakes the iterative baselines nearn≈2​d,3​dn\\approx 2d,3d, eventually reaching parity with and slightly surpassing SURE\-Ridge near the upper end of the sweep, while SURE\-Ridge stays approximately flat at a small, but non\-zero value across the entire range, which is consistent with the small residual bias introduced by its post\-hoc thresholding step\. NOTEARS and DAGMA decrease monotonically, but neither reaches SURE\-Ridge or GBNSL\. Thus, SURE\-Ridge dominates when samples are scarce, GBNSL dominates when samples are plentiful, and the iterative methods are outperformed throughout the small\-sample regime\.

Fig\. 1:SHD vs\. sample sizennacross two problem sizes \(d=20,50d=20,50, top to bottom\)\. Lines and shaded bands show the mean and one standard deviation across 100 Monte Carlo trials per setting\.
### V\-CRun Time Comparison

Fig\. 2:Per\-DAG run time on an Intel i9\-13900F CPU atd=50d=50, on a logarithmic scale\.Figure[2](https://arxiv.org/html/2608.17132#S5.F2)reports per\-DAG run time atd=50d=50on a logarithmic scale, and three distinct tiers are visible\. NOTEARS has slightly improved run time over DAGMA, as the polynomial acyclicity functional is cheaper per gradient step than the log\-determinant constraints in DAGMA\. GBNSL sits roughly an order of magnitude below the iterative methods, dominated by the cost of solving linear programs and regression problems\. SURE\-Ridge is another magnitude below GBNSL, putting it roughly two orders of magnitude faster than the iterative baselines at this problem size\. All four runtimes are essentially flat innn\.

### V\-DDiscussion

SURE\-Ridge’s small\-sample dominance and speed advantage stem from two structural features\. First, sparse DAGs set most node\-wise regression coefficients to zero, the regime where shrinkage is most informative, and SURE adapts its strength without held\-out data\. Second, equal variance decouples the subproblems intoddindependent closed\-form regressions, eliminating the need for gradient descent altogether\. The contrast with GBNSL is a bias\-variance tradeoff\. GBNSL’s OLS estimator is asymptotically unbiased but inherits the variance of the sample covariance and precision estimates, which need enough samples relative to the dimension to stabilize, whereas SURE\-Ridge’s adaptive ridge trades a small bias for a large variance reduction\. The contrast with NOTEARS and DAGMA is computational\. Both pay a fixed gradient\-step budget regardless of problem difficulty, optimize a non\-convex objective vulnerable to finite\-sample noise in the small\-nnregime, and depend on a learning rate, a sparsity weight, and a continuation or augmented\-Lagrangian schedule, none of which can be tuned against ground truth in the causal\-discovery setting\. SURE\-Ridge therefore offers the best accuracy\-compute\-tunability trade\-off in the small\-sample regime targeted by this work\.

## VIConclusions

This work proposes SURE\-Ridge, a non\-iterative, closed\-form estimator for causal discovery in equal\-variance linear Gaussian SEMs that combinesddnode\-wise ridge regressions with SURE\-tuned regularization and an adaptive thresholding procedure for DAG extraction, requiring no per\-problem hyperparameter tuning\. On synthetic Gaussian SEMs, SURE\-Ridge achieves the lowest SHD in the sparse, sample\-limited regime against NOTEARS, DAGMA, and GBNSL baselines, and has a substantially faster run time than all three across every sample size tested, offering an attractive accuracy–compute tradeoff for constrained applications\. Future work will assess the sample\-complexity and convergence of SURE\-Ridge\.

## References

- \[1\]\(2025\)BUILD with precision: bottom\-up inference of linear DAGs\.arXiv preprint arXiv:2512\.16111\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p3.1)\.
- \[2\]K\. Bello, B\. Aragam, and P\. Ravikumar\(2022\)DAGMA: learning dags via m\-matrices and a log\-determinant acyclicity characterization\.Advances in Neural Information Processing Systems35,pp\. 8226–8239\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p2.1),[§V\-A](https://arxiv.org/html/2608.17132#S5.SS1.p2.2)\.
- \[3\]T\. Cai, W\. Liu, and X\. Luo\(2011\)A constrainedL1L\_\{1\}minimization approach to sparse precision matrix estimation\.Journal of the American Statistical Association106\(494\),pp\. 594–607\.Cited by:[§V\-B](https://arxiv.org/html/2608.17132#S5.SS2.p1.1)\.
- \[4\]W\. Chen, M\. Drton, and Y\. S\. Wang\(2019\)On causal discovery with an equal\-variance assumption\.Biometrika106\(4\),pp\. 973–980\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p3.1)\.
- \[5\]D\. M\. Chickering\(2002\)Optimal structure identification with greedy search\.Journal of machine learning research3\(Nov\),pp\. 507–554\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p2.1)\.
- \[6\]A\. Ghoshal and J\. Honorio\(2017\)Learning identifiable gaussian bayesian networks in polynomial time and sample complexity\.Advances in Neural Information Processing Systems30\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p3.1),[§V\-A](https://arxiv.org/html/2608.17132#S5.SS1.p3.1)\.
- \[7\]T\. Hastie, R\. Tibshirani, J\. Friedman,et al\.\(2009\)The elements of statistical learning\.Springer series in statistics New\-York\.Cited by:[§III](https://arxiv.org/html/2608.17132#S3.p3.1)\.
- \[8\]V\. Lungu, J\. Shaska, I\. Kontoyiannis, and U\. Mitra\(2026\)Bayesian structure learning and detection in the linear causal model\.InInternational Symposium on Information Theory,Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[9\]L\. Mata, M\. Sousa, P\. Vieira, M\. P\. Queluz, and A\. Rodrigues\(2025\)Integrating machine learning and observational causal inference for enhanced spectral and energy efficiency in wireless networks\.IEEE Access\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[10\]H\. Nguyen, D\. Tran, B\. Tran, B\. Pehlivan, and T\. Nguyen\(2021\)A comprehensive survey of regulatory network inference methods using single cell RNA sequencing data\.Briefings in bioinformatics22\(3\),pp\. bbaa190\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[11\]J\. Peters and P\. Bühlmann\(2014\)Identifiability of gaussian structural equation models with equal error variances\.Biometrika101\(1\),pp\. 219–228\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1),[§II](https://arxiv.org/html/2608.17132#S2.p1.4),[§III](https://arxiv.org/html/2608.17132#S3.p4.1)\.
- \[12\]K\. J\. Rothman and S\. Greenland\(2005\)Causation and causal inference in epidemiology\.American journal of public health95\(S1\),pp\. S144–S150\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[13\]J\. Shaska and U\. Mitra\(2025\)Causal link discovery with unequal edge error tolerance\.IEEE Transactions on Signal Processing\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[14\]P\. Spirtes and C\. Glymour\(1991\)An algorithm for fast recovery of sparse causal graphs\.Social science computer review9\(1\),pp\. 62–72\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p2.1)\.
- \[15\]C\. M\. Stein\(1981\)Estimation of the mean of a multivariate normal distribution\.The annals of Statistics,pp\. 1135–1151\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p4.1),[§III](https://arxiv.org/html/2608.17132#S3.p2.1)\.
- \[16\]Y\. Vivek, V\. Ravi, A\. Mane, and L\. R\. Naidu\(2024\)Explainable artificial intelligence and causal inference based atm fraud detection\.In2024 IEEE Symposium on Computational Intelligence for Financial Engineering and Economics \(CIFEr\),pp\. 1–7\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[17\]K\. Yu, X\. Wu, W\. Ding, and H\. Wang\(2012\)Exploring causal relationships with streaming features\.The Computer Journal55\(9\),pp\. 1103–1117\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[18\]Y\. Yu, J\. Chen, T\. Gao, and M\. Yu\(2019\)DAG\-GNN: DAG structure learning with graph neural networks\.InInternational conference on machine learning,pp\. 7154–7163\.Cited by:[§IV\-A](https://arxiv.org/html/2608.17132#S4.SS1.p1.1)\.
- \[19\]A\. Zhang, G\. Zhang, B\. Cai, W\. Hu, L\. Xiao, T\. W\. Wilson, J\. M\. Stephen, V\. D\. Calhoun, and Y\. Wang\(2020\)Causal inference of brain connectivity from fMRI withψ\\psi\-learning incorporated linear non\-gaussian acyclic model \(ψ\\psi\-lingam\)\.arXiv preprint arXiv:2006\.09536\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p1.1)\.
- \[20\]X\. Zheng, B\. Aragam, P\. K\. Ravikumar, and E\. P\. Xing\(2018\)DAGs with no tears: continuous optimization for structure learning\.Advances in neural information processing systems31\.Cited by:[§I](https://arxiv.org/html/2608.17132#S1.p2.1),[§IV\-A](https://arxiv.org/html/2608.17132#S4.SS1.p1.2),[§V\-A](https://arxiv.org/html/2608.17132#S5.SS1.p2.1)\.