MARCEDES: Score-based causal discovery under non-Gaussianity with continuous optimization

1 Introduction

• DAG learning remains challenging because the DAG space is super-exponential, acyclicity is combinatorial, and observational data generally identify only Markov equivalence classes.For non-Gaussian models, existing score-based methods must also choose a tractable working error model despite unknown true errors.Source: §1 · pp. 1–2
• Under scale-mixture-of-Gaussians errors, the true causal weighted matrix has almost surely lower asymptotic mean absolute residual risk than any other weighted causal matrix, cyclic or acyclic.The risk is motivated by an intentionally misspecified Laplace error-SEM.Source: §1 · p. 2
• MARCEDES combines the risk with equation-specific sparsity penalties and a soft acyclicity penalty in an unconstrained score-minimization framework.The equation-specific penalties allow more flexible regularization than a single sparsity parameter.Source: §1 · p. 2
• A smooth surrogate and generalized-Bayes MAP interpretation support an alternating minimization algorithm and cross-validation for tuning prior hyperparameters.The approach addresses the objective's nonconvexity and non-smoothness.Source: §1 · p. 2
• Numerical experiments show that MARCEDES improves both structure learning and parameter estimation over existing methods across various evaluation metrics.Source: §1 · p. 2

2.1 Structural equation model

• For each variable X_j, the linear recursive SEM uses a coefficient vector B_j and independent noise ε_j with distribution P_j.Source: §2.1 · p. 2
• The model is written X = BX + ε, where B collects the equation coefficients and ε collects the noise variables.Source: §2.1 · p. 2
• The dataset consists of n independent and identically distributed observations of X generated by the SEM.Source: §2.1 · p. 3
• A nonzero coefficient B_jk encodes a direct linear causal effect of X_k on X_j, and the resulting adjacency matrix represents the causal DAG.Source: §2.1 · p. 3
• Learning the true DAG is the target because it captures conditional-independence relationships in the distribution of X, while estimating B is a more specific objective.Source: §2.1 · p. 3
• Under the work's assumption that all errors are non-Gaussian, LiNGAM establishes unique identifiability of the causal DAG; Gaussian errors can leave Markov-equivalent DAGs distributionally equivalent.Source: §2.1 · p. 3

2.2 Score-based Causal DAG learning

• Score-based learning minimizes a data-dependent score over real matrices subject to the induced adjacency matrix being a DAG.Source: §2.2 · p. 3
• The discrete acyclicity constraint is computationally difficult, while continuous relaxation uses a smooth function that equals zero exactly when the induced graph is a DAG.The relaxed formulation supports gradient-based continuous optimization.Source: §2.2 · p. 3

Mean absolute residual risk

• The mean absolute residual risk evaluates mean absolute residuals when each variable is predicted from the others using a coefficient matrix B.It is motivated by fitting a misspecified Laplace-error SEM.Source: §3 · pp. 3–4
• The risk is connected to hierarchical Bayesian modeling through the likelihood under the misspecified Laplace-error working model.The proposition introduces this connection alongside the risk definition.Source: §3 · p. 4
• For scale-mixture-of-Gaussians errors with finite first moments, the true weighted causal matrix has a strictly smaller asymptotic risk than any distinct matrix.The pairwise risk gap is positive almost surely.Source: §3 · p. 4

Proposed method: MARCEDES

• MARCEDES combines this risk with row-specific sparsity penalties and a soft acyclicity penalty in an unconstrained score-minimization framework.The soft DAG penalty is retained for general non-Gaussian errors and finite-sample performance.Source: §3 · pp. 4–5
• Unlike GOLEM's single global sparsity parameter, MARCEDES assigns a separate sparsity penalty to each structural equation.The paper states that this adaptive regularization is particularly useful in high-dimensional settings and improves estimation accuracy in numerical experiments.Source: §3 · p. 5

4 Optimization

• The algorithm addresses a highly nonconvex objective that remains nonsmooth because of absolute residuals, while also requiring several penalty parameters to be tuned.The section develops a computational procedure for solving the MARCEDES optimization problem.Source: §4 · p. 5

4.1 Generalized Bayes framework

• A generalized Bayes prior models the p equation-specific sparsity parameters through a lower-dimensional hyperparameter θ.The row-wise ℓ1 term encourages sparsity, while h(B) softly penalizes deviations from acyclicity.Source: §4.1 · p. 5

4.2 Gradual enforcement of the DAGness penalty

• The DAG penalty is increased over a grid, with each stage initialized at the previous stage's solution.The gradual increase is intended to encourage broader exploration, including cyclic graphs and nearby DAGs.Source: §4.2 · pp. 5–6

4.3 Core optimization algorithm

• The core algorithm alternates over auxiliary variables, sparsity parameters, and B, using a variational ℓ1 formulation and smooth-gradient steps followed by row-wise shrinkage.For B, it takes an Adam step on the smooth component and applies a row-wise ISTA shrinkage step; diagonal entries are set to zero.Source: §4.3 · pp. 6–7

4.4 Cross-validation over the hyper-parameter

• K-fold cross-validation selects θ from a finite candidate set by minimizing the cross-validated risk, then computes the final estimator.The selected estimator uses the optimization scheme described in Sections 4.2 and 4.3.Source: §4.4 · p. 7

4.5 Empirical Bayes thresholding

• Empirical Bayes thresholding truncates small noisy off-diagonal coefficients to zero, followed by removal of the smallest-weight edges if cycles remain.Figure 1 presents a performance comparison for overall structure recovery.Source: §4.5 · p. 7, Figure 1 · p. 8

5 Numerical Experiments

• The simulations use ER-2 DAGs with Uniform, Laplace, and Student’s t errors, varying n across {50, 200} and p across {10, 20, 50, 100}; structure metrics include TPR, FDR, normalized SHD, and MCC.Parameter estimation accuracy is also assessed.Source: §5 · p. 8

5.1 Structure recovery

• MARCEDES improves structure recovery over competing methods; Gaussian-oriented GOLEM and PC have lower MCC and higher FDR than non-Gaussian methods.Among the non-Gaussian methods, MARCEDES demonstrates superior structure recovery.Source: §5.1 · pp. 8–9

5.2 Edge recovery

• Across small- and large-sample settings, MARCEDES consistently outperforms benchmarks on edge recovery, maintaining high TPR and favorable FDR.Figure 2 presents the edge-recovery comparison.Source: Figure 2 · p. 9, §5.2 · p. 9

5.1 Structure recovery

• With n = 50, MARCEDES improves MCC and normalized SHD; with n = 200, performance improves further, with very small SHD at large sample size and moderate dimension.The paper describes the latter as near-exact recovery of the underlying graph.Source: §5.1 · p. 9
• When p > n, some LiNGAM-based methods fail to produce estimates, while MARCEDES remains applicable and produces estimates.Source: §5.1 · p. 9

• MAR_{CEDES} combines mean absolute residual risk, row-specific sparsity penalties, and a soft DAG constraint in an unconstrained optimization framework.Source: §6 · p. 9
• Simulations show MARCEDES can substantially improve structure learning over existing methods.Source: §6 · p. 9
• Future work includes finding broader risk functions and establishing DAG-selection and parameter-estimation consistency guarantees.Source: §6 · p. 10
• Finite-sample estimates may contain spurious or missing edges and may be affected by latent confounding, measurement error, selection bias, or assumption violations; consequential uses warrant expert validation and additional evidence.Source: §6 · p. 10

• The likelihood is formed under independent Laplace(1) errors, with the joint distribution expressed using det(I − B).Source: §A.1 · p. 13
• Maximizing the Laplace log-likelihood over the scale parameter yields a criterion equal to a constant minus n times L_n(B).Source: §A.1 · p. 13

• The residual vector e is expressed as a transformation of ε: e = (I − B)X = (I − B)(I − B)−1ε = Wε, where W = (I − B)(I − B)−1.Source: §A.2 · p. 14
• The determinant expansion and independence are used to establish an expectation identity for the random matrix Σ.Source: §A.2 · p. 14
• The lemma’s equality condition requires W = PΔ; under DAG acyclicity, the next lemma sharpens this to W = I, equivalently e_j = ε_j for every j.Source: §A.2 · pp. 14–15
• The risk difference converges almost surely to a nonnegative quantity that equals zero only when the candidate matrix equals the true matrix.If det(I − B) = 0, the risk is infinite and the result holds trivially.Source: §A.2 · p. 16

B.1 Simulation setup

• The experiments compare MARCEDES with benchmark methods on data generated from weighted DAGs and linear non-Gaussian SEMs.The error distributions are Laplace(0.8), Student’s t with 5 degrees of freedom, and Uniform[−3, 3], assigned in equal proportions across structural equations.Source: §B.1 · p. 16
• Each experimental setting is repeated over 30 independent replications.Source: §B.1 · p. 16
• The study covers dimensions p ∈ {10, 20, 50, 100} and sample sizes n ∈ {50, 200}.These settings evaluate performance in both low- and high-dimensional regimes with small and large samples.Source: §B.1 · p. 16

B.2 Optimization method and implementation details

• The DAG penalty is introduced gradually using ΛD = {0, 25, 50}, rather than applied at full strength from the start.The stated target DAG-penalty parameter is λD = 50.Source: §B.2 · p. 16
• A preliminary global-penalty cross-validation suggests sparsity levels in (0, 1), typically close to zero, motivating a flexible Beta prior for sparsity parameters.The prior is reparameterized by its mean m and effective sample size τ, which provide interpretable controls for sparsity regularization.Source: §B.2 · pp. 16–17
• Initialization uses PC when p ≥ n and DirectLiNGAM otherwise, while allowing the empty graph as an alternative.The B-subproblem uses one Adam–ISTA inner iteration per outer iteration, with a maximum of 10^5 outer iterations and relative-change tolerance 10^-3.Source: §B.2 · p. 17
• Empirical Bayes thresholding is applied to off-diagonal estimates using the mean of 100 perturbation-based thresholds.The perturbations are Gaussian with standard deviation 0.1; the procedure uses a Laplace prior and median thresholding rule, then enforces a zero diagonal.Source: §B.2 · p. 17
• If thresholding leaves cycles, the algorithm repeatedly deletes the smallest-magnitude edge involved until the estimate is a DAG.This is described as making the smallest possible changes in magnitude to the thresholded coefficient matrix.Source: §B.2 · p. 17