1
Exploiting Non-Negativity in DAG Structure Learning
arXiv:2605.19947v1 [cs.LG] 19 May 2026
Samuel Rey, Member, IEEE, Madeline Navarro, Student Member, IEEE, and Gonzalo Mateos, Senior Member, IEEE
Abstract—This work addresses the problem of learning directed acyclic graphs (DAGs) from nodal observations generated by a linear structural equation model. DAG learning is a central task in signal processing, machine learning, and causal inference, but it remains challenging because acyclicity is a global combinatorial property. Continuous acyclicity constraints have led to important algorithmic advances by replacing the discrete DAG constraint with smooth equality constraints. However, existing formulations still involve difficult non-convex optimization landscapes and may suffer from degenerate first-order optimality conditions. Here, we restrict attention to DAGs with nonnegative edge weights and exploit this additional structure to obtain a simpler characterization of acyclicity. Building on this characterization, we formulate a regularized non-negative DAG learning problem and develop an algorithm based on the method of multipliers. We further analyze the benign optimization landscape induced by non-negativity. In the population regime, we show that the true DAG is the unique global minimizer of the proposed augmented-Lagrangian formulation; moreover, the landscape contains no spurious interior stationary points, and the true DAG is the only acyclic KKT point. Numerical experiments on synthetic and real-world data show that the proposed method improves over state-of-the-art continuous DAGlearning alternatives. Index Terms—DAG learning, network topology inference, causal discovery, graph signal processing, structural equation model
I. I NTRODUCTION Directed acyclic graphs (DAGs) provide a natural language for representing directed dependencies in complex systems, where edge directionality is an essential part of the model [2]. They are central to Bayesian networks and structural equation models [3]–[5], and have become standard tools in biology, genetics, machine learning, signal processing, and causal inference [6]–[13]. In many applications, however, the graph is not available a priori and must be inferred from nodal observations. From the broader graph signal processing viewpoint, this task is part of graph learning, where the goal is to infer network structure from the statistical or structural properties This work was supported by the NSF under Award ECCS 2231036, by the Spanish AEI Grants PID2022-136887NB-I00 and PID2023-149457OB-I00, and by the Community of Madrid (via grants CAM-URJC F1180 (CP2301), TEC-2024/COM-89, and Madrid ELLIS Unit). Part of the results in this paper appeared at the 2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2025) [1]. (Corresponding author: Samuel Rey.) Samuel Rey is with the Dept. of Signal Theory and Communications, Universidad Rey Juan Carlos, Madrid, Spain (e-mail: [email protected]). Madeline Navarro is with the Dept. of Electrical and Computer Engineering, Rice University, Houston, TX 77005, USA (e-mail: [email protected]), and Gonzalo Mateos is with the Dept. of Electrical and Computer Engineering, University of Rochester, Rochester, NY 14627, USA (e-mail: [email protected]).
of observed signals [14]. In the DAG setting, the problem becomes estimating a directed graph that explains the observed data while satisfying the global acyclicity requirement. The latter requirement is the main source of difficulty. Indeed, exact score-based DAG learning is NP-hard in general [15], [16], since acyclicity couples all edges through a combinatorial constraint. A large body of work has addressed this difficulty by designing score functions, greedy search strategies, and maximumlikelihood estimators with favorable statistical properties under suitable identifiability assumptions [4], [17]. These results clarify when the underlying DAG can be recovered in principle, but the corresponding optimization problems remain highly non-convex; as a result, attaining the global optimum is often not computationally straightforward. A major algorithmic step was introduced by NOTEARS [18], which replaced the discrete acyclicity constraint with a smooth equality constraint based on the matrix exponential. This idea opened the door to continuous optimization methods for DAG learning and inspired several refinements, including acyclicity functions based on matrix powers [19], [20], logdeterminant barriers [21], likelihood-based formulations such as GOLEM [22], and more recent continuous formulations such as CoLiDE [23]. Despite their empirical success, these approaches still inherit a challenging non-convex optimization landscape and may exhibit degenerate first-order optimality conditions, so the acyclicity constraint can provide limited local information at feasible DAGs. An alternative route bypasses explicit acyclicity constraints by first estimating a topological order, for instance through variance-based ordering ideas such as VARSORT and related Gaussian-DAG estimators [24], [25]. Order-based approaches can lead to simple and statistically sharp procedures, but their performance depends strongly on the validity and stability of the ordering criterion. Recent alternatives further exploit particular algebraic relations between the DAG adjacency matrix and the precision matrix to guarantee population recovery in linear models [26]. However, these procedures can be sensitive in low-sample regimes. This paper takes a complementary route. Rather than searching for a more sophisticated smooth relaxation over signed weighted graphs, we impose additional structure on the problem and focus on DAGs with non-negative edge weights. The non-negativity assumption removes the cancellations that occur in powers of signed matrices. As a consequence, acyclicity can be characterized by a simpler function applied directly to the adjacency matrix, yielding a smooth acyclicity constraint whose gradient does not vanish at DAGs, thereby avoiding the KKT degeneracy that
2
affects several existing continuous formulations. Building on this observation, we formulate a non-negative DAG learning estimator based on a regularized least-squares score and a log-determinant acyclicity constraint. We solve the resulting constrained problem using the method of multipliers, which provides a principled mechanism for promoting acyclicity. Although the optimization problem remains non-convex, the proposed formulation has a favorable population landscape. In particular, under standard identifiability assumptions, we show that the true DAG is the unique global minimizer of the population augmented Lagrangian, that there are no spurious interior stationary points, and that the only acyclic KKT point is the true DAG. Our numerical results further indicate that this favorable population behavior is reflected in finite-sample regimes, where the proposed method approaches the true graph as the number of observations increases and compares favorably with established continuous DAG-learning baselines. To summarize, our main contributions are as follows. • We show that, in the non-negative setting, acyclicity can be characterized by a simple smooth function that does not rely on Hadamard products, thus avoiding degeneracy of the KKT conditions. • We propose an augmented-Lagrangian algorithm based on the method of multipliers for non-negative DAG structure learning. • We characterize the population optimization landscape and prove that, under identifiability, the true DAG is the unique global minimizer of the population augmented Lagrangian, there are no spurious interior stationary points, and the true DAG is the only acyclic KKT point. • We provide a numerical evaluation showing that the proposed method exhibits consistent finite-sample behavior and competitive performance against continuous DAGlearning baselines. The remainder of the paper is organized as follows. Section II reviews DAG structure learning and continuous acyclicity constraints. Section III introduces the non-negative formulation, the proposed log-determinant acyclicity function, and the method-of-multipliers algorithm. Section IV studies the population optimization landscape. Section V reports numerical experiments, and the appendix contains the proofs of the main population results. II. F UNDAMENTALS OF DAG S TRUCTURE L EARNING This section introduces fundamental concepts related to DAGs and formally states the DAG structure learning problem. A. DAGs, structural equation models, and problem statement Let D = (V, E) be a directed graph with node set V = {1, . . . , d} and directed edge set E ⊆ V × V. An edge (i, j) ∈ E represents a directed link from node i to node j. The graph topology is encoded by a weighted adjacency matrix W ∈ Rd×d , whose entries satisfy Wij ̸= 0 if and only if (i, j) ∈ E. A graph is acyclic if it contains no directed cycles. Throughout, we let D denote the set of matrices in Rd×d whose support corresponds to a DAG.
DAG structure learning aims to infer the unknown connectivity encoded in W from nodal observations. Let X := [x1 , . . . , xn ] ∈ Rd×n collect n observed graph signals, with xi ∈ Rd denoting the i-th observation. We focus on the common linear structural equation model (SEM) X = W⊤ X + Z,
(1)
where Z ∈ Rd×n collects exogenous noise samples. Unless otherwise stated, the columns of Z are assumed independent and identically distributed (i.i.d.) random vectors with covariance matrix Σz = σ 2 I. Independence of the exogenous variables is a crucial requirement for causal interpretation and identifiability in SEMs [4], [17, pp. 83–84]. Given X, the DAG structure encoded in W can be inferred from the observed data by solving the optimization problem min F (W; X) W
s.t.
W ∈ D,
(2)
where F (W; X) is a data-dependent score function that captures the relation between W and the observations. The difficulty in (2) lies in the feasible set D. Directly enforcing acyclicity is a combinatorial endeavor and exact score-based DAG learning is NP-hard in general [15], [16]. B. Continuous acyclicity constraints Modern approaches to DAG structure learning advocate for replacing the discrete constraint W ∈ D with a smooth equality constraint h(W) = 0, where h : Rd×d 7→ R is a differentiable function whose zero level set coincides with D. As a consequence, the acyclicity function h satisfies h(W) = 0 if and only if W ∈ D. The discrete constraint in (2) can then be replaced by the continuous formulation min F (W; X) W
s.t.
h(W) = 0.
(3)
Under this reformulation, DAG learning becomes a continuous constrained optimization problem that can be addressed with first-order methods. Relying on continuous functions to enforce acyclicity was pioneered in [18] via the acyclicity function hnotears (W) = tr eW◦W − d, (4) where ◦ denotes the Hadamard (entry-wise) product. Subsequent works proposed alternative smooth characterizations based on matrix powers [19], [20] and log-determinants. In particular, [21] considers hsdagma (W) = d log(s) − log det (sI − W ◦ W) ,
(5)
with s ∈ R+ chosen so that the log-determinant is well defined. This barrier-type characterization has favorable numerical properties and has been shown to perform well in practice. Continuous acyclicity constraints have substantially changed the algorithmic landscape of DAG structure learning. Replacing the combinatorial constraint W ∈ D with the acyclicity condition h(W) = 0 enables the use of standard continuous optimization methods. This viewpoint has led to scalable algorithms with competitive empirical
3
performance in settings where exhaustive combinatorial search is infeasible. However, the resulting optimization problems remain fundamentally challenging. In particular, the term W ◦ W introduces a non-convex composition that complicates the optimization landscape. Moreover, this term makes the gradient of the acyclicity constraint vanish at every DAG, i.e., ∇h(W) = 0 for every W ∈ D. As observed in [19], this leads to degenerate Karush–Kuhn–Tucker (KKT) conditions ∇F (W) + λ∇h(W) = 0,
(6)
since a DAG W ∈ D can be a KKT point of (3) only if it is an unconstrained stationary point of the score function, i.e., ∇F (W) = 0. Consequently, when F is convex, this condition holds only at unconstrained global minimizers of F . To overcome these limitations, we henceforth assume that W has non-negative weights and propose a simpler acyclicity function that does not depend on the Hadamard product. This yields a more amenable optimization landscape, avoids the aforementioned degeneracy at feasible DAGs, and opens the door to stronger guarantees for the global behavior of the resulting estimator.
We now specialize the DAG structure learning problem to the case where the edge weights are non-negative. This restriction is natural in settings where directed interactions encode excitatory, additive, or non-inhibitory relations, and includes the important case of binary adjacency matrices; see e.g., [13]. From an optimization viewpoint, non-negativity removes the sign cancellations that arise in matrix powers, enabling acyclicity to be characterized through simpler smooth functions that do not require the Hadamard product. When the observations follow the linear SEM in (1), a standard choice is to estimate the weighted adjacency matrix by combining a least-squares data-fitting term with a sparsitypromoting penalty. Since W is constrained to be entrywise non-negative, the usual ℓ1 penalty reduces to a linear penalty in the entries of W. This yields the estimator
W
s.t.
X 1 ∥X − W⊤ X∥2F + αn Wij 2n
(7)
i̸=j
W ≥ 0,
h(W) = 0,
Smooth acyclicity constraints are central to modern continuous approaches for DAG learning. As discussed in [18], an effective acyclicity function h(W) should be smooth, have a gradient that can be evaluated efficiently, and satisfy h(W) = 0 if and only if W ∈ D. By exploiting the non-negativity of W and drawing inspiration from the logdeterminant characterization in (5), we introduce a zero-level acyclicity constraint that satisfies these requirements without relying on Hadamard products. Proposition 1. For any matrix W ∈ Rd×d whose spectral + radius satisfies ρ(W) < s for some s ∈ R+ , define h(W) := d log(s) − log det(sI − W).
(8)
The proposed acyclicity constraint is h(W) = 0, and its gradient is given by ∇h(W) = (sI − W)−⊤ .
(9)
Moreover, h(W) ≥ 0 for every W ∈ Rd×d such that ρ(W) < + s, and h(W) = 0
⇐⇒
W ∈ D.
(10)
Proof. The gradient expression in (9) follows directly from standard tools from differential calculus. We now prove the acyclicity characterization. First, rewrite h as
III. N ON -N EGATIVE DAG S TRUCTURE L EARNING
Ŵ = arg min
A. Acyclicity over non-negative matrices
ρ(W) < s,
where αn > 0 is a tunable parameter that controls the tradeoff between data fidelity and sparsity and is chosen so that αn → 0 as n → ∞. Moreover, h is a continuous acyclicity function, and the spectral-radius condition specifies the domain on which h, introduced below, is well defined. Note that the diagonal entries are excluded from the sparsity penalty since self-loops are incompatible with acyclicity. Equivalently, one may explicitly impose diag(W) = 0. By writing the estimator as a single global minimizer, (7) anticipates a uniqueness property. Although this property is not immediate from the formulation, it will be established in Section IV.
h(W) = d log(s) − log(sd ) − log det(I − s−1 W) = − log det(I − s−1 W).
(11)
Since ρ(s−1 W) < 1, the matrix logarithm is well defined, and applying the Mercator series gives log(I − s−1 W) = −
∞ X s−k Wk k=1
k
.
(12)
Then, using the identity log det(M) = tr(log M), yields ! ∞ X s−k Wk −1 h(W) = −tr log I − s W = tr k k=1
=
∞ X tr(Wk ) k=1
ksk
.
(13)
Because W is entrywise non-negative, every diagonal entry of Wk is non-negative, and hence tr(Wk ) ≥ 0 for all k ≥ 1. It follows from (13) that h(W) ≥ 0. It remains to characterize the zero level set. Since all terms in (13) are non-negative, h(W) = 0 holds if and only if tr(Wk ) = 0 for every k ≥ 1. The diagonal entry [Wk ]ii is the total weight of all directed closed walks of length k that start and end at node i. Thus, if the graph associated with W contained a directed cycle, the trace of a suitable power of W would be strictly positive. Therefore, h(W) = 0 implies that W represents a DAG. Conversely, if W ∈ D, then W is nilpotent up to a permutation, so tr(Wk ) = 0 for every k ≥ 1 and (13) gives h(W) = 0. The non-negativity assumption is essential in Proposition 1. Indeed, the series in (13) expresses h as a weighted sum
4
of closed-walk contributions. When edge weights are nonnegative, every directed cycle contributes a non-negative term and cannot be masked by cancellations from other closed walks. As a result, the zero level set of h coincides exactly with the set of DAGs. This cancellation-free interpretation is precisely what is lost when signed weights are allowed. Although h is not convex, it retains a favorable first-order structure that is absent from Hadamard-based acyclicity functions. In particular, (9) does not vanish at DAGs, so feasible acyclic matrices are not automatically stationary points of the acyclicity function, hence preventing the degeneracy of the KKT conditions. B. Algorithmic implementation via the method of multipliers We solve (7) using the method of multipliers, an augmented-Lagrangian approach for equality-constrained optimization [27, Ch. 4.2]. The appeal of this method is that, in settings where the required regularity and convexity assumptions hold, minimizing a sequence of augmented Lagrangians recovers a solution of the original constrained problem. In the present setting, the acyclicity constraint remains nonconvex, so these classical guarantees cannot be invoked directly. Nevertheless, the method provides a principled way to drive the acyclicity violation to zero while retaining a smooth objective at each iteration. Let us denote the score function in (7) by Fn (W) :=
X 1 ∥X − W⊤ X∥2F + αn Wij . 2n
(14)
i̸=j
For a multiplier λ ∈ R and a penalty parameter c > 0, the augmented Lagrangian associated with the acyclicity constraint is c Lc (W, λ) := Fn (W) + λh(W) + h(W)2 . 2
(15)
Only the acyclicity constraint is incorporated into (15). The non-negativity constraint is kept explicit, since it can be handled directly by projection onto the non-negative orthant. Similarly, the spectral-radius condition is treated as a domain restriction for the log-determinant term. Starting from λ(0) and c(0) > 0, the method generates a sequence of iterates by repeating the three following steps. Step 1. First, the weighted adjacency matrix is updated by solving W(k+1) =
arg min
Lc(k) (W, λ(k) ).
(16)
W≥0, ρ(W)<s
In practice, this inner problem can be addressed with projected first-order iterations. For the log-determinant acyclicity function in Proposition 1, the gradient of the augmented Lagrangian is ∇W Lc (W, λ) = ∇Fn (W) + (λ + ch(W)) ∇h(W). (17) The score-gradient term is given by ∇Fn (W) = Σ̂x (W − I) + αn (11⊤ − I),
(18)
where Σ̂x = n−1 XX⊤ is the sample covariance matrix. Starting from W0 = W(k) and using an inner stepsize η > 0, the projected-gradient update for (16) is h i Wℓ+1 = Wℓ − η∇W Lc(k) (Wℓ , λ(k) ) , (19) +
where [·]+ denotes the entrywise projection onto the nonnegative orthant. The stepsize is chosen so that the iterates remain in the domain ρ(Wℓ ) < s of the log-determinant term. Once the inner stopping criterion is met, the next outer iterate is set to the last inner iterate, i.e., W(k+1) = Wℓ+1 . Alternatively, accelerated methods such as FISTA can also be employed [28], [29]. Step 2. Second, the Lagrange multiplier is updated according to the current constraint violation, λ(k+1) = λ(k) + c(k) h(W(k+1) ).
(20)
This update can also be interpreted as a gradient ascent step since the constraint violation corresponds to the gradient of Lc(k) (W(k+1) , λ) with respect to λ. Step 3. Finally, the penalty parameter is increased only when the acyclicity violation is not reduced sufficiently. A standard update is ( βc(k) , if h(W(k+1) ) > γh(W(k) ), (k+1) c = (21) c(k) , otherwise, where β > 1 and 0 < γ < 1. This adaptive rule increases the penalty only when the current value of c(k) is not producing a sufficient decrease in the acyclicity residual. The resulting procedure follows the same logic as classical multiplier methods: the quadratic term penalizes violations of h(W) = 0, while the multiplier update corrects the linearization of the constraint across iterations. The computational complexity is determined by the solution of the inner problem in (16). Each projected-gradient step requires evaluating (17), with cost O(d3 ) due to the computation of (sI − W)−⊤ and the associated matrix products. In practice, we observe that only a small number of outer iterations is needed to satisfy the acyclicity constraint. The nonconvexity of the h constraint prevents us from invoking the classical convergence theory for augmented Lagrangian methods. Nevertheless, the numerical evaluation in Section V suggests that, in practice, the estimate provided by the proposed algorithm converges to the true DAG as the number of samples goes to infinity. To further understand this promising behavior, we next study the optimization landscape of the problem at hand. IV. P OPULATION L ANDSCAPE A NALYSIS This section analyzes how the non-negativity assumption and the simpler acyclicity constraint yield a benign optimization landscape. To that end, we now shift our attention to the augmented Lagrangian of the DAG learning problem in the population regime. Let W0 denote the weighted adjacency matrix of the true non-negative DAG, and let x be a generic observation generated according to the SEM in (1). At the population level,
5
the empirical least-squares term is replaced by its expectation, leading to the score function F̄ (W) := E ∥x−W⊤x∥22 = tr (I−W)⊤ Σx (I−W) , (22) where Σx = E[xx⊤ ] denotes the population covariance matrix of the observations. The identifiability of W0 is tied to the exogenous noise distribution. In the analysis below, we focus on the normalized case where the exogenous noise has identity covariance, meaning that Σz = E[zz⊤ ] = I. The extension to homoscedastic noise with covariance σ 2 I, or to heteroscedastic noise with known covariance up to a scalar factor, is immediate after the corresponding rescaling. What is essential is that the noise model ensures identifiability of the true DAG in the population regime [4], [17]. Considering the score function in (22) and the acyclicity constraint in Proposition 1, the population augmented Lagrangian for a multiplier λ ∈ R+ and a penalty parameter c > 0 is given by c (23) L̄c (W, λ) := F̄ (W) + λh(W) + h(W)2 . 2 The results below characterize the global and first-order behavior of (23) over the non-negative domain of the logdeterminant constraint. For the landscape analysis, we focus on the acyclicity constraint with s = 1, corresponding to the function h(W) = − log det(I − W). The feasible set over which the augmented Lagrangian will be optimized can be compactly written as W1 := W ∈ Rd×d : ρ(W) < 1 . (24) + Under this setting, the population DAG learning model is identifiable in the following sense [17]: W0 = arg min F̄ (W) W
s.t.
W ∈ W1 , h(W) = 0. (25)
We next characterize the population augmented Lagrangian through three complementary properties: its global minimizer, its classical stationary points, and its acyclic KKT points. The first result states that the true DAG W0 is the unique global minimizer of L̄c (·, λ) over the full domain W1 . Theorem 1 (Unique global minimizer). Assume that the data follow the SEM in (1), where W0 ∈ W1 is the weighted adjacency matrix of the true DAG and the covariance of the exogenous noise is E[zz⊤ ] = I. Let the identifiability condition in (25) hold. Then, for every c > 0 and every λ ≥ 2, W0 is the unique global minimizer of the population augmented Lagrangian over W1 , namely W0 = arg min L̄c (W, λ).
(26)
W∈W1
See Section A for the proof. Theorem 1 shows that the true DAG structure encoded in W0 is the only matrix in W1 that minimizes the population augmented Lagrangian. Thus, the identifiability of W0 over the constrained DAG set is lifted to the augmented-Lagrangian objective over the larger nonnegative domain. However, global uniqueness alone does not describe what first-order methods may encounter in a non-convex landscape.
We therefore turn to classical stationary points of the smooth augmented-Lagrangian objective. Theorem 2 (Classical stationary points). Assume that the data follow the SEM in (1), where W0 ∈ W1 is the weighted adjacency matrix of the true DAG and the covariance of the exogenous noise is E[zz⊤ ] = I. Let the identifiability condition in (25) hold, and let c > 0. Then the following statements hold: (i) If λ > 2, there is no W ∈ W1 such that ∇W L̄c (W, λ) = 0. (ii) If λ = 2, then ∇W L̄c (W, 2) = 0
⇐⇒
W = W0 .
(27)
See Section B for the proof. Theorem 2 rules out spurious interior stationary points of the population augmented Lagrangian. Thus, any additional first-order candidates relevant to projected methods must be understood through the boundary of the non-negative cone, which motivates the KKT analysis below. Since the spectral-radius condition in W1 is open, the only active constraint in the first-order analysis over W1 is W ≥ 0. Thus, a matrix W ∈ W1 is a KKT point of the problem of minimizing L̄c (·, λ) over W1 if and only if there exists a multiplier matrix Γ ∈ Rd×d such that + ∇W L̄c (W, λ) − Γ = 0,
Γ ◦ W = 0.
(28)
Equivalently, ∇W L̄c (W, λ) ≥ 0,
W ◦ ∇W L̄c (W, λ) = 0.
(29)
Theorem 3 (Acyclic KKT points). Assume that the data follow the SEM in (1), where W0 ∈ W1 is the weighted adjacency matrix of the true DAG and the covariance of the exogenous noise is E[zz⊤ ] = I. Let the identifiability condition in (25) hold, and let c > 0. If W ∈ W1 is a KKT point of L̄c (·, 2) over W1 and h(W) = 0, then W = W0 .
(30)
See Section C for the proof. In summary, Theorems 1 to 3 show that the non-negativity assumption, combined with the log-determinant acyclicity constraint, induces a favorable population landscape: W0 is the unique global minimizer of the augmented Lagrangian, there are no spurious interior stationary points, and every acyclic KKT point is the true DAG. These results do not by themselves guarantee convergence of a first-order method to W0 , but they indicate that the population augmented Lagrangian is more amenable to optimization than formulations based on Hadamard-product acyclicity constraints. They also help explain the numerical behavior in Section V, where the estimation error decreases as the number of samples grows. V. N UMERICAL EVALUATION We now provide a numerical evaluation of the proposed method in finite-sample regimes. We first evaluate performance in a controlled synthetic setting, where the groundtruth DAG and the data-generating process are known. We then assess the method on the Sachs protein-signaling benchmark,
6
which provides a real-data test case with a widely used reference DAG. We measure estimation accuracy through the normalized Frobenius error ∥Ŵ − W0 ∥2F nerr(Ŵ, W0 ) := , ∥W0 ∥2F
(31)
and structural accuracy through the structural Hamming distance (SHD), normalized by the number of nodes. The SHD counts the number of edge insertions, deletions, and reversals needed to match the support of the true DAG. We compare the proposed log-determinant method, denoted as “Logdet”, with three continuous DAG-learning baselines: NOTEARS [18], DAGMA [21], and CoLiDE [23]. For completeness, we also include the “Matexp” variant from the conference experiments, which replaces the log-determinant acyclicity function with its matrix-exponential counterpart while retaining the nonnegativity constraint. The synthetic study is designed to isolate three complementary effects: the number of available samples, the size and topology of the underlying DAG, and the variance of the exogenous noise. Unless otherwise stated, the data are generated from Erdős–Rényi DAGs with d = 100 nodes, average degree equal to 4, and n = 1000 samples from the linear SEM in (1) with standard Gaussian exogenous noise. Each curve reports the median together with the 25th and 75th percentiles over 100 independent realizations. Test case 1 - Sample complexity. Figure 1(a) shows the normalized Frobenius error as the number of samples increases. The proposed Logdet method consistently attains the lowest estimation error among the compared methods. More importantly, its error keeps decreasing as n grows, in contrast with the saturation observed for the non-convex baselines. This trend is consistent with the population result in Theorem 1: although the theorem does not itself provide a finitesample rate, the experiment suggests that the finite-sample estimator approaches the correct population solution as more observations become available. The comparison with Matexp also indicates that the log-determinant barrier is numerically preferable to the matrix-exponential alternative in this setting. Test case 2 - Graph size and topology. Figure 1(b) evaluates structural recovery as the number of nodes increases. The comparison focuses on Logdet and DAGMA, since both use a log-determinant acyclicity principle but lead to different optimization problems. For ER graphs, Logdet recovers the correct support with essentially zero normalized SHD across the considered graph sizes. For SF graphs, the task becomes more challenging for both methods, but the degradation is substantially milder for Logdet. This supports the practical value of exploiting non-negativity to obtain a better-conditioned logdeterminant formulation, rather than relying on Hadamardproduct acyclicity constraints. Test case 3 - Noise variance. Figure 1(c) studies the effect of increasing the exogenous noise variance in the homoscedastic setting Σz = σ 2 I. As expected, the estimation error of methods that do not explicitly account for the noise level increases with σ 2 . The variant “Logdet-σ”, which uses the
TABLE I: Results on the Sachs protein-signaling benchmark. Lower values are better for SHD, FDR, Err., and time; higher values are better for TPR and F1. Method SHD TPR FDR F1 Err. Time (s) NOMAD (ours) 10 0.412 0.000 0.583 0.898 34.12 CoLiDE-Fix 13 0.412 0.588 0.412 1.984 5.62 DAGMA 15 0.118 0.600 0.182 1.439 1.35
known noise covariance, remains stable across the tested range and performs similarly to CoLiDE. This behavior is aligned with the discussion in the population analysis: what matters for identifiability is not the normalization Σz = I itself, but having a noise model that makes the true DAG identifiable. When the covariance structure is known, it can be incorporated directly into the score function. A. Sachs protein-signaling benchmark We next evaluate the proposed method on the Sachs proteinsignaling dataset [6], a standard real-data benchmark for DAG structure learning and causal discovery. The dataset contains single-cell measurements of 11 phosphorylated proteins and phospholipids under experimental perturbations, together with a reference signaling network validated by biological knowledge. The benchmark DAG contains 11 nodes and 17 directed edges, and the observational subset used here contains 853 samples. Because it combines real biological measurements with a widely used reference DAG, this dataset is the preferred real-data benchmark for comparing DAG-learning methods in this work. For this experiment, we refer to the proposed method as NOMAD, short for Non-negative Optimization via Multipliers for Acyclic Digraphs. Table I reports structural and estimation metrics for NOMAD and two continuous DAG-learning baselines, DAGMA [21] and CoLiDE [23]. NOMAD attains the lowest SHD, FDR, and estimation error, while matching CoLiDE-Fix in TPR. Although CoLiDE-Fix is faster in this small-scale setting, its recovered graph contains a substantially larger fraction of false discoveries. The graph visualizations in Fig. 2 illustrate this behavior: NOMAD recovers a sparse estimate aligned with the reference network and, according to Table I, achieves an FDR of zero. VI. C ONCLUSIONS This paper studied DAG structure learning under a nonnegativity assumption on the edge weights. This additional structure provides more than a modeling prior: by ruling out sign cancellations, it enables a simple smooth characterization of acyclicity applied directly to the weighted adjacency matrix, leading to a constraint whose gradient remains informative at feasible DAGs. Building on this observation, we formulated a regularized non-negative DAG-learning problem and proposed NOMAD, an augmented-Lagrangian algorithm based on the method of multipliers. We further showed that the resulting formulation has a favorable population optimization landscape: under identifiability, the true DAG is the unique global minimizer, there are no spurious interior stationary points, and the true DAG is the only acyclic KKT point. Experiments on
7
Normalized SHD
nerr(Ŵ, W0 )
1.5 10−1
10−2
Logdet Matexp CoLiDE DAGMA NoTears
10−3 50 100 500 1000 5000 (a) Number of samples n
Logdet-ER DAGMA-ER
Logder-SF DAGMA-SF
10−1 nerr(Ŵ, W0 )
100
1
0.5
0 50
100 250 (b) Number of nodes d
500
Logdet CoLiDe
Logdet-σ DAGMA
10−2
2
4
6
8
10
(c) Exogenous noise variance σ 2
Fig. 1: Synthetic evaluation of the proposed DAG-learning method. (a) Normalized Frobenius error as the number of samples increases. (b) Normalized SHD as the graph size increases, for Erdős–Rényi (ER) and scale-free (SF) DAGs. (c) Normalized Frobenius error as the variance of the exogenous noise increases.
(a) Reference DAG
(b) DAGMA estimate
(c) NOMAD estimate
Fig. 2: Qualitative comparison on the Sachs protein-signaling benchmark. The proposed NOMAD method recovers a graph that is structurally closer to the reference DAG than DAGMA and introduces no false discoveries, in agreement with the SHD, F1, and FDR scores reported in Table I.
controlled synthetic data and on the Sachs protein-signaling benchmark support these theoretical findings, showing that exploiting non-negativity can improve both weighted estimation and structural recovery relative to established continuous DAG-learning baselines. Overall, the results suggest that incorporating meaningful sign information can turn DAG learning into a substantially better-conditioned optimization problem, providing a promising route for robust structure learning in applications where interactions are known to be non-inhibitory. A PPENDIX A P ROOF OF T HEOREM 1 We first state a technical lower bound that reduces the analysis of the augmented Lagrangian to a scalar function of the acyclicity value. Moreover, the following definitions will be used in the different proofs. M := I − W,
M0 := I − W0 ,
C := M−1 0 M. (32)
where c ϕλ,c (t) := de−2t/d + λt + t2 , t ≥ 0. (34) 2 Moreover, ϕλ,c (t) ≥ d for every t ≥ 0, with equality if and only if t = 0. Proof. The proof proceeds in two steps. First, we relate the population score to the acyclicity value by rewriting the score via the matrix C, which enables relating its determinant to h(W). Second, after substituting this bound into the augmented Lagrangian, it remains only to analyze a onedimensional function of h(W). We begin by rewriting the population score. Since the exogenous noise has identity covariance, the SEM implies that the covariance matrix of the data is given by −1 Σx = (I − W0 )−⊤ (I − W0 )−1 = M−⊤ 0 M0 .
(35)
Therefore, the score function can be expressed as F̄ (W) = tr(M⊤ Σx M) = tr(C⊤ C).
(36)
Lemma 1. Assume that the data follow the SEM in (1), where W0 ∈ W1 is the weighted adjacency matrix of the true DAG and the exogenous noise has identity covariance. Let c > 0 and λ ≥ 2. Then, for every W ∈ W1 ,
Moving our attention to the acyclicity function, first notice that W0 is a DAG. It follows that W0 is a nilpotent matrix, and hence det(M0 ) = det(I − W0 ) = 1. Thus,
L̄c (W, λ) ≥ ϕλ,c (h(W)),
det(C) = det(M−1 0 M) = det(M).
(33)
(37)
8
Since W ∈ W1 , the log-determinant is well defined and h(W) = − log det(I − W) = − log det(M) = − log det(C).
(38)
Since W is a global minimizer of L̄c , it follows that W also minimizes F̄ over the feasible DAG set. The identifiability condition in (25) then implies W = W0 . This proves that W0 is the unique global minimizer.
Letting σ1 , . . . , σd denote the singular values of C, and noticing that det(C) > 0 since det(M) > 0, we have d Y
σi = det(C) = e−h(W) .
(39)
i=1
Next, to relatePthe score and the acyclicity function, notice d 2 that F̄ (W) = i=1 σi . Applying the arithmetic-geometric mean inequality to the non-negative numbers σ12 , . . . , σd2 yields !1/d d d X Y F̄ (W) = σi2 ≥ d = de−2h(W)/d . (40) σi2 i=1
i=1
Substituting this bound into the definition of L̄c and using the function ϕλ,c yields c L̄c (W, λ) = F̄ (W) + λh(W) + h(W)2 2 c ≥ de−2h(W)/d + λh(W) + h(W)2 2 = ϕλ,c (h(W)), (41) which reduces the augmented-Lagrangian bound to the scalar function ϕλ,c evaluated at the acyclicity value. It remains to show that this scalar lower bound is never below its acyclic value. For t ≥ 0, the first and second order derivatives are ϕ′λ,c (t) = −2e−2t/d + λ + ct, (42) 4 −2t/d ′′ + c > 0. (43) ϕλ,c (t) = e d Thus ϕλ,c is strictly convex and ϕ′λ,c is strictly increasing. Since λ ≥ 2, we have ϕ′λ,c (0) = λ − 2 ≥ 0. Hence ϕλ,c is minimized at t = 0 over [0, ∞), and strict convexity gives ϕλ,c (t) ≥ ϕλ,c (0) = d,
(44)
with equality if and only if t = 0. This concludes the proof. We now leverage the result from Lemma 1 to first prove that W0 is a global minimizer of L̄c (·, λ) over W1 , and then prove that it is the unique global minimizer. Since W0 is a DAG, h(W0 ) = 0. Moreover, from (I − W0⊤ )x = z and E[zz⊤ ] = I, it follows that L̄c (W0 , λ) = F̄ (W0 ) = E[∥z∥22 ] = d.
(45)
Proof. The proof uses the scalar function introduced in Lemma 1 to analyze the first-order equation ∇W L̄c (W, λ) = 0. The main idea is to express the augmented Lagrangian in terms of the normalized matrix C = M−1 0 M [see (32)]. In these coordinates, stationarity forces all singular values of C to be equal. Combining this with the determinant identity encoded by h(W) reduces the matrix stationarity condition to a scalar equation involving ϕ′λ,c (h(W)). As in the proof of Lemma 1 [see (36) and (38)], leveraging the definitions in (32) we can rewrite the score and the acyclicity function as F̄ (W) = tr(C⊤ C),
h(W) = − log det(C).
(48)
Thus, with a slight abuse of notation, the population augmented Lagrangian can be written as a function of the matrix C, leading to 2 c L̄c (C, λ) = tr(C⊤ C)−λlog det(C)+ log det(C) . (49) 2 The affine map W 7→ C = M−1 0 (I−W) is invertible. Consequently, ∇W L̄c (W, λ) = 0 if and only if ∇C L̄c (C, λ) = 0. Differentiating with respect to C gives ∇C L̄c (C, λ) = 2C − λC−⊤ + c log det(C)C−⊤ = 2C − λ + ch(W) C−⊤ , (50) where the last equality uses the equivalence h(W) = − log det(I − W) = − log det(C).
(51)
Hence the stationarity condition ∇C L̄c (C, λ) = 0 amounts to 2C = λ + ch(W) C−⊤ . (52) We now turn this matrix equation into a scalar condition. Multiplying both sides of (52) on the right by C⊤ yields 2CC⊤ = λ + ch(W) I, (53) which implies that all singular values of C are equal. Letting this common singular value be σ, and using the relation in (51) yields det(C) = e−h(W) > 0, (54) which in turn results in
For any W ∈ W1 and λ ≥ 2, Lemma 1 gives L̄c (W, λ) ≥ ϕλ,c (h(W)) ≥ d = L̄c (W0 , λ).
A PPENDIX B P ROOF OF T HEOREM 2
σ = e−h(W)/d .
(46)
Therefore, W0 is a global minimizer of L̄c (·, λ) over W1 . To prove uniqueness, let W ∈ W1 be any global minimizer. Then L̄c (W, λ) = d. The inequalities above can be tight only if h(W) = 0. By Proposition 1, this means that W is a DAG. f ∈ W1 satisfies h(W) f = Furthermore, every feasible DAG W 0, so f λ) = F̄ (W). f L̄c (W, (47)
(55)
Since all singular values of C are equal to σ, the left-hand side of (53) is 2CC⊤ = 2σ 2 I. Thus, comparing the scalar multiplying I on both sides of (53) gives 2e−2h(W)/d = λ + ch(W).
(56)
ϕ′λ,c h(W) = 0,
(57)
Equivalently,
9
where ϕλ,c is the scalar function in Lemma 1, whose derivative is given in (42). It remains to analyze this scalar equation. Since W ∈ W1 , Proposition 1 gives h(W) ≥ 0. For every t ≥ 0, ϕ′′λ,c (t) =
4 −2t/d e + c > 0, d
Therefore, (60)
Since h(W) = 0, Proposition 1 implies that W is a DAG. Thus W is feasible for the identifiable population DAG problem in (25). Moreover, F̄ (W) = d = F̄ (W0 ), where the last equality follows from the SEM relation at the true DAG and identity noise covariance. By the identifiability condition in (25), we conclude that W = W0 . Conversely, at W = W0 we have h(W0 ) = 0. Using ∇W F̄ (W0 ) = −2M−⊤ 0 ,
∇h(W0 ) = M−⊤ 0 ,
(61)
∇W L̄c (W0 , 2) = ∇W F̄ (W0 ) + 2∇h(W0 ) = 0.
(62)
we obtain
This proves that, when λ = 2, the only classical stationary point is W0 . A PPENDIX C P ROOF OF T HEOREM 3 Proof. The proof has two steps. First, we use the KKT conditions to derive a necessary inequality relating the population score to the acyclicity value. Second, when the KKT point is acyclic, this inequality forces the score to match the value attained by the true DAG, and identifiability then yields W = W0 . Let W ∈ W1 be a KKT point of L̄c (·, 2) over W1 , and let M = I − W. Define A(W) :=
1 ∇W L̄c (W, 2). 2
(63)
Using ∇W F̄ (W) = −2Σx M,
∇h(W) = M−⊤ ,
0 = ⟨W, A(W)⟩ = tr(A(W)) − ⟨M, A(W)⟩.
(67)
Hence, (58)
so ϕ′λ,c is strictly increasing on [0, ∞). If λ > 2, then ϕ′λ,c (0) = λ − 2 > 0. Therefore ϕ′λ,c (t) > 0 for every t ≥ 0, and (57) has no solution. This proves that no classical stationary point exists in W1 when λ > 2. Consider now the threshold case λ = 2. Then ϕ′2,c (0) = 0, and strict monotonicity implies that the only solution of (57) on [0, ∞) is h(W) = 0. Substituting h(W) = 0 into (53) gives CC⊤ = I. (59) F̄ (W) = tr(C⊤ C) = d.
We next convert complementarity into a scalar inequality. From W = I − M and (66),
(64)
we can write c A(W) = −Σx M + 1 + h(W) M−⊤ . (65) 2 Since W is KKT, it must satisfy the complementarity conditions A(W) ≥ 0, W ◦ A(W) = 0. (66)
tr(A(W)) = ⟨M, A(W)⟩.
(68)
Substituting (65) into the right-hand side gives c ⟨M, A(W)⟩ = −⟨M, Σx M⟩ + 1 + h(W) ⟨M, M−⊤ ⟩ 2 c = −F̄ (W) + d 1 + h(W) . (69) 2 Since A(W) ≥ 0, its trace is non-negative, and (68) yields c (70) F̄ (W) ≤ d 1 + h(W) . 2 We now specialize to acyclic KKT points. If h(W) = 0, then (70) gives F̄ (W) ≤ d. On the other hand, Lemma 1 with λ = 2 gives L̄c (W, 2) = F̄ (W) ≥ d, (71) where we used h(W) = 0 in the equality. Therefore F̄ (W) = d. Since h(W) = 0, Proposition 1 implies that W is a DAG. Moreover, the true DAG satisfies F̄ (W0 ) = d by the SEM relation and identity noise covariance. Thus W is a feasible DAG attaining the same population score as W0 . The identifiability condition in (25) then implies W = W0 . R EFERENCES [1] S. Rey, S. S. Saboksayr, and G. Mateos, “Non-negative weighted dag structure learning,” in IEEE Intl. Conf. Acoustics, Speech and Signal Process. (ICASSP). IEEE, 2025, pp. 1–5. [2] A. G. Marques, S. Segarra, and G. Mateos, “Signal processing on directed graphs: The role of edge directionality when processing and learning from network data,” IEEE Signal Process. Mag., vol. 37, no. 6, pp. 99–116, 2020. [3] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles and Techniques. MIT Press, 2009. [4] J. Peters, D. Janzing, and B. Schölkopf, Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press, 2017. [5] P. Spirtes, C. Glymour, and R. Scheines, Causation, Prediction, and Search. MIT Press, 2001. [6] K. Sachs, O. Perez, D. Pe’er, D. A. Lauffenburger, and G. P. Nolan, “Causal protein-signaling networks derived from multiparameter singlecell data,” Science, vol. 308, no. 5721, pp. 523–529, 2005. [7] P. J. F. Lucas, L. C. Van der Gaag, and A. Abu-Hanna, “Bayesian networks in biomedicine and health-care,” Artif. Intell. Med., vol. 30, no. 3, pp. 201–214, 2004. [8] B. Zhang, C. Gaiteri, L.-G. Bodea, Z. Wang, J. McElwee, A. A. Podtelezhnikov, C. Zhang, T. Xie, L. Tran, R. Dobrin et al., “Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer’s disease,” Cell, vol. 153, no. 3, pp. 707–720, 2013. [9] Y. Yu, J. Chen, T. Gao, and M. Yu, “DAG-GNN: DAG structure learning with graph neural networks,” in Proc. Intl. Conf. Mach. Learn. (ICML), 2019, pp. 7154–7163. [10] S. Rey, H. Ajorlou, and G. Mateos, “Directed acyclic graph convolutional networks,” IEEE Transactions on Signal Processing, 2026. [11] P. Misiakos, V. Mihal, and M. Püschel, “Learning signals and graphs from time-series graph data with few causes,” in Proc. IEEE Intl. Conf. Acoust., Speech Signal Process. (ICASSP), 2024, pp. 9681–9685. [12] L. Yao, Z. Chu, S. Li, Y. Li, J. Gao, and A. Zhang, “A survey on causal inference,” ACM Trans. Knowl. Discovery Data, vol. 15, no. 5, pp. 1–46, 2021. [13] B. Seifert, C. Wendler, and M. Püschel, “Causal Fourier analysis on directed acyclic graphs and posets,” IEEE Trans. Signal Process., vol. 71, pp. 3805–3820, 2023. [14] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro, “Connecting the dots: Identifying network structure via graph signal processing,” IEEE Signal Process. Mag., vol. 36, no. 3, pp. 16–43, 2019.
10
[15] D. M. Chickering and D. Heckerman, “Efficient approximations for the marginal likelihood of Bayesian networks with hidden variables,” Mach. Learn., vol. 29, pp. 181–212, 1997. [16] D. M. Chickering, D. Heckerman, and C. Meek, “Large-sample learning of Bayesian networks is NP-hard,” J. Mach. Learning Res., vol. 5, pp. 1287–1330, 2004. [17] P.-L. Loh and P. Bühlmann, “High-dimensional learning of linear causal networks via inverse covariance estimation,” J. Mach. Learning Res., vol. 15, no. 140, pp. 3065–3105, 2014. [18] X. Zheng, B. Aragam, P. K. Ravikumar, and E. P. Xing, “DAGs with NO TEARS: Continuous optimization for structure learning,” Conf. Neural Inform. Process. Syst., vol. 31, 2018. [19] D. Wei, T. Gao, and Y. Yu, “DAGs with No Fears: A closer look at continuous optimization for learning Bayesian networks,” Conf. Neural Inform. Process. Syst., vol. 33, pp. 3895–3906, 2020. [20] R. Pamfil, N. Sriwattanaworachai, S. Desai, P. Pilgerstorfer, K. Georgatzis, P. Beaumont, and B. Aragam, “DYNOTEARS: Structure learning from time-series data,” in Intl. Conf. Artif. Intel. Statist. (AISTATS), 2020, pp. 1595–1605. [21] K. Bello, B. Aragam, and P. Ravikumar, “DAGMA: Learning DAGs via M-matrices and a log-determinant acyclicity characterization,” Conf. Neural Inform. Process. Syst., vol. 35, pp. 8226–8239, 2022. [22] I. Ng, A. Ghassami, and K. Zhang, “On the role of sparsity and DAG constraints for learning linear DAGs,” Conf. Neural Inform. Process. Syst., vol. 33, pp. 17 943–17 954, 2020. [23] S. S. Saboksayr, G. Mateos, and M. Tepper, “CoLiDE: Concomitant linear DAG estimation,” Intl. Conf. Learn. Repr. (ICLR), 2024. [24] A. G. Reisach, C. Seiler, and S. Weichwald, “Beware of the Simulated DAG! Causal Discovery Benchmarks May Be Easy to Game,” Conf. Neural Inform. Process. Syst., vol. 34, pp. 27 772–27 784, 2021. [25] M. Gao, W. M. Tai, and B. Aragam, “Optimal estimation of Gaussian DAG models,” arXiv preprint arXiv:2201.10548, 2022. [26] H. Ajorlou, S. Rey, G. Mateos, G. Leus, and A. G. Marques, “Build With Precision: Bottom-Up Inference of Linear DAGs,” in IEEE Intl. Conf. Acoustics, Speech and Signal Process. (ICASSP). IEEE, 2026, pp. 446–450. [27] D. P. Bertsekas, Nonlinear Programming, 3rd ed. Belmont, MA, USA: Athena Scientific, 2016. [28] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci., vol. 2, no. 1, pp. 183–202, 2009. [29] A. Beck, First-Order Methods in Optimization. Philadelphia, PA, USA: SIAM, 2017.