#  dotears: Scalable, consistent DAG estimation using observational and
interventional data

Albert Xue1 Jingyou Rao1 Sriram Sankararaman1∗ Harold Pimentel1∗
1University of California, Los Angeles
{asxue, roserao<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
(February 15, 2024)

###### Abstract

New biological assays like Perturb-seq link highly parallel CRISPR
interventions to a high-dimensional transcriptomic readout, providing insight
into gene regulatory networks. Causal gene regulatory networks can be
represented by directed acyclic graph (DAGs), but learning DAGs from
observational data is complicated by lack of identifiability and a
combinatorial solution space. Score-based structure learning improves
practical scalability of inferring DAGs. Previous score-based methods are
sensitive to error variance structure; on the other hand, estimation of error
variance is difficult without prior knowledge of structure. Accordingly, we
present dotears [doo-tairs], a continuous optimization framework which
leverages observational and interventional data to infer a single causal
structure, assuming a linear Structural Equation Model (SEM). dotears exploits
structural consequences of hard interventions to give a marginal estimate of
exogenous error structure, bypassing the circular estimation problem. We show
that dotears is a provably consistent estimator of the true DAG under mild
assumptions. dotears outperforms other methods in varied simulations, and in
real data infers edges that validate with higher precision and recall than
state-of-the-art methods through differential expression tests and high-
confidence protein-protein interactions.

**footnotetext: These authors contributed equally to this work.

##  1 Introduction

Understanding gene regulatory networks can identify mechanisms and pathways
linking GWAS significant variants to phenotype. Recent efforts to map
regulatory networks through trans-eQTLs are partly limited by power; for
example, the GTEx project finds only 143 trans-eQTLs in 838 individuals [1].
On the other hand, Võsa et al. detect almost 60,000 trans-eQTLs in
$\sim$31,000 individuals, across more than a third of trait-associated
variants [2]. These results imply that trans-regulatory relationships are
pervasive, but our ability to detect small eQTL effects is often limited by
small sample size regimes, observational data, and a high multiple testing
burden. Importantly, since rare tissues are unlikely to be sampled at
sufficient sample sizes, capturing gene regulation events across a wide array
of cell types and tissues requires another experimental method.

High-throughput genomic technologies like Perturb-seq provide a natural
alternative for learning gene regulatory networks. Perturb-seq links high-
dimensional transcriptomic readouts to known, highly parallel CRISPR
interventions, allows direct interrogation of causal regulatory relationships,
and has scaled genome-wide [3, 4, 5]. In particular, the effects of CRISPR
gene interventions are large in comparison to QTL effects, faciliating
inference of downstream regulatory relationships. Notably, analogous
experiments have already mapped gene-gene networks in yeast [6, 7, 8].

The inference of gene regulatory networks can be treated as a causal structure
learning problem, which considers learning relationships between variables in
the form of a Directed Acyclic Graph (DAG). Identifiability and scalability
are the primary difficulties in learning DAGs from data. For identifiability,
distinct DAGs may contain the same conditional independence relationships in
observational data, and DAGs are only identifiable up to Markov equivalence
[9, 10, 11]. For scalability, DAGs with NO TEARS introduces a continuous,
differentiable acyclicity constraint, which avoids combinatorial
characterizations of DAGs and allows for continuous optimization over
structure [12]. However, NO TEARS and related methods infer DAGs whose
topological order follows increasing marginal variance, and re-scaling data
can change or reverse their inferences [13].

Fundamentally, this is still an issue of identifiability. Because NO TEARS
uses observational data, it must choose a single member of a class of Markov
equivalent DAGs. The “tiebreaker” is then a function of the variance. However,
we show that interventional data can correct for variance sensitivity in NO
TEARS and related methods, and further that this correction is sufficient for
consistent estimation of structure.

Explicitly, exogenous error in the linear SEM drives variance sensitivity. Let
$X$ be a $p$-dimensional random vector (e.g. the distribution of gene
expression across $p$ genes), $W_{0}\in\mathbb{R}^{p\times p}$ the weighted
adjacency matrix of the true DAG, and $\epsilon$ a $p$-dimensional random
vector specifying the exogenous error. The linear SEM gives the autoregressive
formulation

$X=XW_{0}+\epsilon,$

where
$\Omega_{0}\coloneqq\text{Cov}\left(\epsilon\right)=\text{diag}\left(\sigma_{1}^{2},\dots,\sigma_{p}^{2}\right)$.
For any node $j$ we then have

$\text{Var}\left(X_{j}\right)=\sum_{i=1}^{p}w_{ij}^{2}\text{Var}\left(X_{i}\right)+\sigma_{j}^{2}.$
(1)

$\text{Var}\left(X_{i}\right)$ is itself a linear function of
$\sigma_{1}^{2},\dots,\sigma_{p}^{2}$. Critically, parental error variances
thus propagate to downstream nodes to provide a signal of structure. However,
in observational data the exact relationship between $W_{0}$, $\Omega_{0}$,
and $\text{Var}(X)$ is unidentifiable.

Consequently, given $\Omega_{0}$ then $W_{0}$ is recoverable even in
observational data [14]. However, estimation of $\Omega_{0}$ is difficult
without a priori knowledge of $W_{0}$. Previous methods either ignore
exogenous variance structure or use the conditional estimate
$\hat{\Omega}_{0}\mid W$ [12, 15]. We show that both procedures are sensitive
to exogenous variance structure even in the simplest two node DAG.

Hard interventions [16, 17] remove upstream variance in the linear SEM to
allow marginal estimation of $\Omega_{0}$ (Figure 1). We show that naive
incorporation of interventional data into the NO TEARS framework, without
estimation of $\Omega_{0}$, is insufficient for structural recovery. Finally,
we show correction by this marginal estimate is sufficient for structural
recovery.

Figure 1: Hard interventions allow marginal estimation of the error variances
$\Omega_{0}\coloneqq\text{diag}\left(\sigma_{1}^{2},\sigma_{2}^{2}\right)$.
(a) The true observational DAG $X_{1}\overset{w}{\to}X_{2}$, with
corresponding distributions. (b) Hard interventions remove incoming edges to
the target. The marginal variance shrinks due to removal of upstream variance
and effects from the intervention.

Accordingly, we present dotears, a novel optimization framework for structure
learning. dotears uses 1.) a novel marginal estimation procedure for
$\Omega_{0}$ using the structural consequences of interventions, and 2.) joint
estimation of the causal DAG from observational and interventional data, given
the estimated $\hat{\Omega}_{0}$. dotears provides a simple model that we
show, by extending results from [14], is a provably consistent estimator of
the true DAG under mild assumptions. In simulations dotears corrects for
exogenous variance structure and is robust to reasonable violations of its
modeling assumptions. We also apply dotears to the Perturb-seq experiment in
Replogle et al. [5]. dotears infers a sparse set of edges that validate with
high precision in differential expression tests and in an orthogonal set of
protein-protein interactions. In both simulations and real data, dotears
outperforms all other tested methods in all used metrics.

##  2 Model

###  2.1 Formulation

Let $G=\left([p],\,\mathcal{E}\right)$ a DAG on $p$ nodes with node set
$[p]\coloneqq\\{1\dots p\\}$ and edge set $\mathcal{E}$. We represent
$\mathcal{E}$ with the weighted adjacency matrix
$W\in\mathcal{D}\subset\mathbb{R}^{p\times p}$, where
$\mathcal{D}\subset\mathbb{R}^{p\times p}$ is the set of weighted adjacency
matrices on $p$ nodes whose support is a DAG. We denote the parent set of node
$i$ in the observational setting as $\text{Pa}(i)$. For $w_{ij}$ the $i,j$
entry in $W$, $|w_{ij}|>0$ indicates an edge $i\rightarrow j$ with weight
$w_{ij}$, equivalently denoted $i\overset{w_{ij}}{\rightarrow}j$.
$k=0,1,\dots,p$ indexes the intervention, where $k=0$ is reserved for the
observational system, and $k\neq 0$ denotes intervention on node $k$.
Similarly, we denote $\left(\cdot\right)^{(0)}$ for observational quantities,
and $\left(\cdot\right)^{(k)}$ for quantities under intervention on node $k$.
For brevity, a variable without a superscript is assumed to be observational;
for example, $X\equiv X^{(0)}$. $\boldsymbol{X}$ (bolded) denotes $n_{0}$
samples drawn from the $p$-dimensional random vector $X$ (unbolded), and
$\boldsymbol{\epsilon}$ (bolded) denotes $n_{0}$ samples drawn from the
$p$-dimensional random vector $\epsilon$ (unbolded). Similarly, if $X^{(k)}$
is a $p$-dimensional random vector, then
$\boldsymbol{X}^{(k)}\in\mathbb{R}^{n_{k}\times p}$ represents $n_{k}$
observations of $X$. We denote the total sample size
$n\coloneqq\sum_{k=0}^{p}n_{k}$, and the true weighted adjacency matrix
$W_{0}$.

The linear SEM is an autoregressive representation of $X^{(k)}$ and weighted
adjacency matrix $W_{0}^{(k)}$,

$X^{(k)}=X^{(k)}W_{0}^{(k)}+\epsilon^{(k)},\quad k=0,\dots,p.$ (2)

Here, $W_{0}^{(k)}$ is permutation-similar to a strictly upper triangular
matrix, representing the constraint $W_{0}^{(k)}\in\mathcal{D}$. For each $k$,
$\epsilon^{(k)}$ is a $p$-dimensional random vector such that
$\mathbb{E}\epsilon^{(k)}=0_{p}$, and
$\Omega_{0}^{(k)}\coloneqq\text{Cov}(\epsilon^{(k)})$. Denote $\epsilon_{i}$
as the $i$th element of $\epsilon$. Then $\epsilon_{i}^{(k)}$ is the exogenous
error on node $i$, such that $\mathbb{E}\epsilon_{i}^{(k)}=0$ and
$\epsilon_{i}^{(k)}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\epsilon_{j}^{(k)}$ for $i\neq j$. We further
define

$\displaystyle\Omega_{0}$ $\displaystyle\coloneqq\text{Cov}(\epsilon^{(0)})$
$\displaystyle=\text{diag}\left(\sigma_{1}^{2},\sigma_{2}^{2},\cdots,\sigma_{p}^{2}\right).$

###  2.2 Interventions

Motivated by recent work on a genome-wide screen that performs known single
interventions on all protein-coding genes [5], we consider the linear SEM with
known single interventions on all $p$ nodes. Accordingly, we obtain a system
of $p+1$ structural equations

$\displaystyle X^{(0)}$ $\displaystyle=X^{(0)}W^{(0)}+\epsilon^{(0)}$
$\displaystyle\vdots$ $\displaystyle X^{(p)}$
$\displaystyle=X^{(p)}W^{(p)}+\epsilon^{(p)}.$

In this setting we have complicated our problem. Before, with the single data
matrix $\boldsymbol{X}$, we inferred a single $W$; now, with the $p+1$ data
matrices
$\boldsymbol{X}^{(0)},\boldsymbol{X}^{(1)},\dots,\boldsymbol{X}^{(p)}$, it
seems that we must infer $p+1$ adjacency matrices
$W^{(0)},W^{(1)},\dots,W^{(p)}$.

Our model assumes hard interventions, i.e. that an intervention on node $k$
removes causal influences from observational parents of $k$ [16, 17]. Hard
interventions relax the do operation $\texttt{do}\left(X^{(k)}_{k}=0\right)$
to allow for residual noise, modeling the limit of interventional efficacy
combined with a noisy readout.

######  Assumption 1.

For an intervention $k\neq 0$, if in the observational setting

$X_{k}^{(0)}=\sum_{i\in\text{Pa}(k)}W_{ik}X_{i}^{(0)}+\epsilon_{k}^{(0)},$

then upon intervention on $k$

$X_{k}^{(k)}=\epsilon_{k}^{(k)}.$

Under Assumption 1, we can relate $W^{(0)}$ to $W^{(k)}$ by setting the $k$th
column of $W^{(0)}$ to $\vec{0}_{p}$, giving

$W^{(k)}_{ij}\coloneqq\begin{cases}W_{ij}&j\neq k\\\ 0&j=k\end{cases}.$ (3)

This gives $p+1$ data matrices $\boldsymbol{X}^{(k)}$ to jointly infer a
single weighted adjacency matrix $W$.

We now characterize the $p+1$ exogenous variance structures
$\Omega_{0}^{(k)}$. We assume that the exogenous variance of non-targeted
nodes $j\neq k$ is invariant.

######  Assumption 2.

For an intervention $k\neq 0$,
$\text{Var}\left(\epsilon_{j}^{(k)}\right)=\text{Var}\left(\epsilon_{j}^{(0)}\right)=\sigma_{j}^{2}$
for $j\neq k$.

We allow interventions to have effects on the error variance of the target,
but require the effect to be uniform across targets.

######  Assumption 3.

Let unknown $\alpha\in\mathbb{R}$. Then $\forall\,\,k$, if
$\text{Var}\left(\epsilon_{k}^{(0)}\right)=\sigma_{k}^{2}$ then
$\text{Var}\left(\epsilon_{k}^{(k)}\right)=\frac{\sigma_{k}^{2}}{\alpha^{2}}$.

Here, $\alpha$ is shared across interventions. As $\alpha\to\infty$, this
interventional model is equivalent to $\texttt{do}(X_{k}=0)$ [17]. Under these
assumptions, the variance of the target is

$\text{Var}\left(X_{k}^{(k)}\right)=\frac{\sigma_{k}^{2}}{\alpha^{2}}.$

Let $\widehat{\text{Var}}$ denote the unbiased sample variance. We then obtain
the estimator

$\hat{\Omega}_{0}=\text{diag}\left(\widehat{\text{Var}}\left(X_{1}^{(1)}\right),\dots,\widehat{\text{Var}}\left(X_{p}^{(p)}\right)\right).$
(4)

##  3 dotears

DAGs with NO TEARS [12] transforms the combinatorial constraint
$W\in\mathcal{D}$ into the continuous constraint $h(W)=0$, where $\circ$
denotes the Hadamard product and

$h(W)=\textrm{tr}\left[\exp\left({W\circ W}\right)\right]-p.$

Define $\|\cdot\|_{1}$ as the vector $\ell_{1}$ norm on
$\text{vec}\left(W\right)$, i.e. $\|\text{vec}\left(W\right)\|_{1}$, and
$\left|\left|{\cdot}\right|\right|_{F}$ the Frobenius norm. For some loss
function $\mathcal{F}$, the differentiability of $h$ allows for the
optimization framework

$\displaystyle\min_{W}$
$\displaystyle\,\,\mathcal{F}\left(W,\boldsymbol{X}\right)+\lambda\|W\|_{1}$
(5) s.t. $\displaystyle\,\,h(W)=0.$

We present dotears, a consistent, intervention-aware joint estimation
procedure for structure learning. Loh and Buhlmann (2014) showed that the
Mahalanobis norm is a consistent estimator of $W_{0}$, and is uniquely
minimized in expectation by $W_{0}$ given $\Omega_{0}$, but give no estimation
procedure for $\Omega_{0}$ [14, 18]. Note the Mahalanobis norm’s
characterization as inverse-variance-weighted by $\Omega_{0}$,

$\left|\left|\left(\boldsymbol{X}-\boldsymbol{X}W\right)\Omega_{0}^{-\frac{1}{2}}\right|\right|_{F}^{2}=\sum_{i=1}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}-\boldsymbol{X}W\right)_{i}\right|\right|_{F}^{2}.$
(6)

dotears solves the following optimization problem:

$\displaystyle\min_{W}$
$\displaystyle\,\,\frac{1}{p}\sum_{k=0}^{p}\frac{1}{2n_{k}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)\hat{\Omega}_{0}^{-\frac{1}{2}}\right|\right|_{F}^{2}+\lambda\|W\|_{1}$
(7) s.t. $\displaystyle h(W)=0,$

where

$W_{ij}^{(k)}=\begin{cases}W_{ij}&j\neq k\\\ 0&j=k\end{cases}.$

dotears retains the continuous DAG constraint and $\ell_{1}$ regularization of
$W$ from NO TEARS [12], but incorporates exogenous variance structure through
$\hat{\Omega}_{0}$ as well as interventional data $(k=1,\dots,p)$.

###  3.1 dotears successfully corrects for exogenous variance structure

We show that dotears is robust to exogenous variance structure, and motivate
the necessity of marginal estimation of $\hat{\Omega}_{0}$ using the simplest
non-trivial DAG $X_{1}\overset{w}{\rightarrow}X_{2}$. Assume
$X_{1}\overset{w}{\to}X_{2}$ has true weighted adjacency matrix
$W_{0}\coloneqq\begin{pmatrix}0&w\\\ 0&0\end{pmatrix}$ and SEM

$\displaystyle X_{1}$ $\displaystyle=\epsilon_{1},$
$\displaystyle\text{Var}\left(\epsilon_{1}\right)=\sigma_{1}^{2},$ (8)
$\displaystyle X_{2}$ $\displaystyle=wX_{1}+\epsilon_{2},$
$\displaystyle\text{Var}\left(\epsilon_{2}\right)=\sigma_{2}^{2}.$

Let $\gamma\coloneqq\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}$, such that
$\Omega_{0}=\begin{pmatrix}\gamma&0\\\ 0&1\end{pmatrix}$ and
$W_{0}=\begin{pmatrix}0&w\\\ 0&0\end{pmatrix}$. The least squares loss used by
NO TEARS is minimized in expectation if and only if

$|w|\geq\sqrt{1-\frac{1}{\gamma}},$ (9)

which is true if and only if the topological ordering of the DAG follows
increasing marginal variance, or equivalently a varsortable DAG. For the full
proof, see Supplementary Material 7.1, or [13, 18].

In Figure 2, we examine the performance of four different strategies of
correcting for $\Omega_{0}$ in simulations. NO TEARS (black) uses the least
squares loss, which ignores $\Omega_{0}$, while GOLEM-NV (orange) uses a
likelihood loss that estimates $\hat{\Omega}_{0}\mid W$ [12, 15]. We also
include the scenario when $\hat{\Omega}_{0}$ is set to $I_{p}$ in Eq. 7. Call
this NO TEARS interventional (green), the simplest extension of NO TEARS to
interventional data. NO TEARS interventional is aware of hard interventions
through the structure $W^{(k)}$, but ignores $\Omega_{0}$. NO TEARS
interventional is thus an ablation study on $\hat{\Omega}_{0}$.

For each set $(w,\gamma)\in\\{0.1,\dots,1.5\\}\times\\{1,2,4,10,100\\}$, we
draw 25 simulations of observational and interventional data, with sample size
$n=(p+1)*1000=3000$. For observational data, this is $3000$ observations from
the observational system; for interventional data, this is $1000$ observations
from each system $k=0,1,2$. For observational data, we draw Gaussian data
under the SEM in Eq. 8. For interventional methods, we draw Gaussian data
under the system of SEMs in Supplementary Material 7.4. To isolate the
behavior of the loss, we remove $\ell_{1}$ regularization. Full simulation
details and results are given in Supplementary Material 7.4.

NO TEARS does not correct for $\Omega_{0}$ and uses only observational data.
As a result, in Figure 2 it estimates correctly only on varsortable pairs of
$w,\gamma$. Note that the grey dashed line represents the theoretical
varsortability cutoff $|w|\leq\sqrt{1-\frac{1}{\gamma}}$ given in Eq. 9.

GOLEM-NV uses the maximum likelihood estimate $\hat{\Omega}\mid W$ under a
Gaussian model. Subsequently, joint estimation is performed over
$W,\hat{\Omega}\mid W$ using the Gaussian negative log likelihood and profile
likelihood for $\Omega_{0}$, which simplifies to

$\frac{1}{2}\sum_{i=1}^{p}\log\left(\left|\left|(X-XW)_{i}\right|\right|_{F}^{2}\right).$
(10)

This profile likelihood is insufficient to correct for exogenous variance
structure, and only infers varsortable structures in Figure 2. This evidence
qualitatively holds in simulations on three node topologies (Supplementary
Material 7.5), where the behavior of GOLEM-NV remains deterministic in
$w,\gamma$. Joint estimation of $W,\hat{\Omega}_{0}\mid W$ is thus still
sensitive to exogenous variance structure.

Through NO TEARS interventional, we also see that interventional data alone,
without correction for $\Omega_{0}$, is insufficient to infer the two node
structure. The NO TEARS interventional estimate is also deterministic in
$w,\gamma$ in Figure 2, and behaves almost identically to NO TEARS and GOLEM.
Note that since NO TEARS interventional is given interventional data for all
nodes, it operates in a fully identifiable setting (where NO TEARS and GOLEM-
NV do not).

Thus, neither interventional data nor correction by $\hat{\Omega}_{0}\mid W$
are alone sufficient to infer structure. dotears combines the two to give the
marginal estimate of $\Omega_{0}$ in Eq. 4 and a robust estimate of $W$.

We do not imply that other methods using interventional data cannot infer the
two node case under interventional data; in fact, many do successfully (see
Supplementary Material 7.4). Rather, we use 1.) an observational procedure
that ignores $\Omega_{0}$, 2.) an observational procedure that corrects for an
estimated $\hat{\Omega}_{0}$, 3.) an interventional procedure that ignores
$\Omega_{0}$, and 4.) dotears, an interventional procedure that corrects for
an estimated $\hat{\Omega}_{0}$, to motivate dotears as most parsimonious
model robust to exogenous variance structure under this line of thought.

(a) Under the equal variance assumption $\gamma=1$, correction for
$\Omega_{0}$ is unneeded. All methods are sufficient for structure recovery.

(b) At $\gamma=2$, NO TEARS (black), GOLEM-NV (orange), and NO TEARS
interventional (green) infer correctly on varsortable $w$.

(c) As $\gamma$ grows large, the varsortability bound approaches $|w|\geq 1$.
Once again, only dotears estimates correctly for all $w$.

Figure 2: Comparison of $\ell_{1}$ distance (lower is better) between true
structure and estimates from NO TEARS (black), GOLEM-NV (orange), NO TEARS
interventional (green), and dotears (blue). Each method corrects differently
for $\Omega_{0}$. For each $w=0.1,0.2,\dots,1.0$ and $\gamma=1,2,100$, we
generate Gaussian data from the structure $X_{1}\overset{w}{\to}X_{2}$ such
that $\sigma_{1}^{2}=\gamma\sigma_{2}^{2}$. For each pair $w,\gamma$, we draw
25 simulations at a sample size of $n=3000$. The dashed grey line represents
the varsortability bound $|w|\geq\sqrt{1-\frac{1}{\gamma}}$. Bars represent
standard errors; some standard errors are too small to see.

###  3.2 Optimization and consistency

We now wish to show that dotears is a consistent estimator of the true DAG.
However, two natural problems arise from our usage of $\hat{\Omega}_{0}$.
First, we have provided no estimation procedure for $\alpha$, but
$\mathbb{E}\hat{\Omega}_{0}\neq\Omega_{0}$ for $\alpha\neq 1$. However,
$\mathbb{E}\hat{\Omega}_{0}\propto\Omega_{0}$ for all $\alpha$, and constant
scalings of $\Omega$ are rescalings of the loss [14]. $\hat{\Omega}_{0}$ is
therefore well-specified for inference on observational data $k=0$.

However, if $\alpha\neq 1$ then $\hat{\Omega}_{0}$ is still misspecified for
interventional data $k\neq 0$. Under Assumption 3

$\text{Cov}\left(\epsilon^{(k)}\right)=\text{diag}\left(\sigma_{1}^{2},\sigma_{2}^{2},\dots,\frac{\sigma_{k}^{2}}{\alpha^{2}},\dots,\sigma_{p}^{2}\right),$

and thus
$\mathbb{E}\hat{\Omega}_{0}\not\propto\text{Cov}\left(\epsilon^{(k)}\right)$.
A naive approach might estimate $\hat{\Omega}_{0}$ from interventional data
only and then estimate $\hat{W}$ from observational data only. However, this
approach ignores a majority of our data and performs substantially worse in
simulations (Supplementary Material 7.6.6). We show that $\hat{\Omega}_{0}$ is
well-specified even for $k\neq 0$, and estimation of $\alpha$ is unnecessary
(see Corollary 2).

Under a sub-Gaussian assumption, Loh and Buhlmann show consistency of the
Mahalanobis norm on observational data given $\Omega_{0}$ [14]. We extend
these results, using $\hat{\Omega}_{0}\overset{p}{\to}\Omega_{0}$, to show

$\operatorname*{arg\,min}_{W}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)\hat{\Omega}_{0}^{-\frac{1}{2}}\right|\right|_{F}^{2}$

is a consistent estimator of $W_{0}^{(k)}$ for each $k$, where
$\overset{p}{\to}$ denotes convergence in probability. A full proof is given
in Supplementary Material 7.2.

##  4 Simulations

We evaluate structure learning methods across a range of DAG topologies,
effects distributions, and generative models. We benchmark methods that
leverage interventional data (dotears, GIES, IGSP, UT-IGSP, DCDI) and methods
using only observational data (NO TEARS, sortnregress, GOLEM-EV, GOLEM-NV,
DirectLiNGAM) [10, 12, 13, 15, 19, 20, 21, 22, 23]. dotears outperforms all
tested methods in DAG estimation and is robust to reasonable violations of the
model.

Figure 3: Method performance on large random graphs ($p$ = 40) using
Structural Hamming Distance (lower is better). Rows index Erdős-Rényi or Scale
Free topologies. Columns index parameterizations of edge density and weight,
ordered in increasing difficulty. For details see Supplementary Material
7.6.1. 10 simulations were drawn for each parameterization with sample size
$(p+1)*100=4100$. * indicates cross-validated methods. Methods are sorted by
average performance.

Some methods come with important caveats for evaluation. sortnregress is not
intended as a “true” structure learning method, but benchmarks the data’s
varsortability [13]. UT-IGSP can infer structure with unknown interventional
targets, but we constrain to known targets for fairness [21]. Most simulations
use Gaussian data, but dotears, NO TEARS, sortnregress, GIES, IGSP, and UT-
IGSP do not assume Gaussianity. The non-Gaussianity assumption is violated for
Direct-LiNGAM, and the equal variance assumption for GOLEM-EV [15, 23].

###  4.1 Large random graph simulations

We simulate synthetic data from large Erdős-Rényi and Scale-Free DAGs [24, 25]
($p=40$), with 10 replicates each. We simulate under four parameterizations:
$\\{\text{Low Density, High Density}\\}\times\\{\text{Weak Effects, Strong
Effects}\\}$. Observational and interventional data have matched sample size
$n=(p+1)\times 100=4100$. Methods using 5-fold cross-validation for
hyperparameter tuning are denoted by *. For non-binary methods, edge weights
are thresholded. Methods are robust to threshold choice (Supplementary
Material 7.6.3) [13]. For simulation, cross-validation, and thresholding
details, see Supplementary Material 7.6.1; for runtime and memory usage, see
Supplementary Material 7.6.2. We note that the memory usage of DCDI-G was
extreme even with no hidden layers.

On average, dotears outperforms all tested methods in structural recovery
(Structural Hamming Distance (SHD), Figure 3) and edge weight recovery
($\ell_{1}$ distance, Supplementary Figure 13). Furthermore, dotears
outperforms all other methods in most parameterizations. GIES matches dotears
in “Low Density” simulations, but performs substantially worse in “High
Density” simulations. We hypothesize that the greedy nature of GIES hinders
performance in more complex DAGs.

The primary modeling assumptions of dotears are 1.) hard interventions
(Assumption 1), 2.) shared $\alpha$ across interventions (Assumption 3), and
3.) linearity of the SEM. In Supplementary Material 7.7.1 and 7.7.2 we assess
the sensitivity of dotears to each assumption. In addition, in Supplementary
Material 7.7.1 we also consider 4.) simulations under different interventional
models.

We find that dotears performance can change under violations of the hard
intervention assumption, but is robust to violations of its interventional
model and linearity. Surprisingly, dotears is the second best performing
method under a mean-shift intervention model. Moreover, dotears outperforms
the neural network method DCDI-DSF in nonlinear simulations, even with
“imperfect” interventions. Under “imperfect” interventions the primary
determining factor is not nonlinearity, but the denseness of the DAG.

##  5 Results

Figure 4: a) dotears-inferred network. Edges with magnitude less than 0.2, and
genes without inferred edges, were removed. b) Precision-recall curves across
differential expression calls made by DESeq2. Dashed red lines indicate recall
of dotears at thresholds of $|w|<0.2,0.1,$ and $0.05$ respectively. c)
Precision-recall curves across high confidence protein-protein interactions
nominated by STRING. Dashed red lines indicate recall of dotears at thresholds
of $|w|<0.2,0.1,$ and $0.05$ respectively. d) dotears infers HSP90AB1
$\rightarrow$ HSP90AA1. HSP90AB1 knockdown increases expression of HSP90AA1,
but HSP90AA1 knockdown does not change HSP90AB1 expression. d) dotears
inferred edges show correlated gene expression in hold-out observational data.

We apply all benchmarked methods in Section 4 to a genome-wide Perturb-seq
experiment from Replogle et al. [5]. We validate inferred edges through 1.)
differential expression tests in the training data using DESeq2 [26, 27] and
2.) an orthogonal set of high-confidence protein-protein interactions from the
STRING database [28]. We also examine gene-gene correlations in held-out
observational expression data. High-confidence edges inferred by dotears show
differential expression and/or protein-protein interactions more frequently
than those found by other methods, and dotears outperforms all other methods
in precision and recall under reasonable thresholding.

Replogle et al. provide both normalized and raw data. We select the top 100
most variable genes in the raw observational data. We then benchmark all
methods on the normalized, feature-selected data. Cross-validation is not
performed due to low sample sizes in some knockdowns; instead, the $\ell_{1}$
penalty is arbitrarily set to 0.1 for all methods where appropriate. Figure 4a
shows the network inferred by dotears, thresholded at $|w|<0.2$. For full
details, see Supplementary Material 7.10.

We use DESeq2 differential expression calls and high-confidence protein-
protein interactions from the STRING database to validate inferred edges [26,
27, 28]. For differential expression, we call an edge $i\rightarrow j$ a true
positive if either gene shows differential expression under knockdown of the
other. This is because all methods struggle equally with predicting
directionality (see Supplementary Table 7). For protein-protein interactions,
we take high-confidence physical interactions as true positives, where here
“high-confidence” is defined by STRING as having a confidence level of over
70% [28]. Figures 4b and 4c show precision and recall across thresholds for
differential expression calls and the protein-protein interactions,
respectively. Vertical lines indicate different thresholding regimes.

dotears shows much higher precision than all other methods at equivalent
recall. Over 65% of inferred edges validated by either differential expression
or high-confidence protein-protein interactions. GIES, IGSP, UT-IGSP, and DCDI
are excluded because they infer binary edges. These methods inferred 3038,
3064, 3075, and 2039 out of a possible 4950 edges, respectively; no other
method predicted more than 700 even at a weight threshold of .05. Accordingly,
they had almost random precision - see Supplementary Material 7.10.3 for
detailed results for all methods at multiple thresholds. These results
reinforce concerns about the scalability of GIES to more complex scenarios.
For DCDI, we report intermediate results, since convergence was not obtained
under 24 hours on GPU training. CPU training attempts ran out of memory even
after allocation of 180GB; GPU training attempts also repeatedly ran out of
memory.

We re-run dotears on the same data excluding the observational data, and use
the held-out observational data to validate inferred edges in the knockdown
data. For each pair of genes $i,j$, Figure 4e shows the inferred edge weight
against the $r^{2}$ of $i,j$ in the observational data. Note that we take the
largest magnitude weight between $w_{ij}$ and $w_{ji}$. Figure 4e shows a
clear relationship between inferred edge weight in knockdown data and the
$r^{2}$ in observational expression data.

##  6 Discussion

We present dotears, a structure learning framework that uses interventional
data to estimate exogenous variance structure and subsequently leverages
observational and interventional data to learn the causal graph. We showed
that dotears is appropriate for Peturb-seq data analysis and can recover high-
confidence gene regulation events. Additionally, we verify that dotears is
robust to exogenous variance structure, prove that its loss function provides
consistent DAG estimation, and show that dotears outperforms all tested state-
of-the-art methods in simulations.

In simulations simple methods generally outperform complex methods in
structure recovery, even under complex generating mechanisms. In particular,
dotears and sortnregress regularly outperform more complex methods, including
neural network methods, even under modeling assumption violations and
nonlinear data. Their strong performance also points to the efficacy of using
patterns in variance to infer structure.

In real data, dotears infers DAGs with edges supported by knockdown expression
in training data, orthogonal high-confidence protein-protein interactions, and
gene expression correlations in held out observational data. dotears-inferred
edges validate with higher precision than any other method without sacrificing
power. In particular, edges with high inferred edge weights - in other words,
the edges dotears infers with highest “confidence” - are validated more
frequently than other methods’ inferences, regardless of validation scheme. We
value performance in these “high-confidence” edges most highly, since these
edges are the top candidates for follow-up experimental validation.

## References

  * [1] GTEx Consortium.  The gtex consortium atlas of genetic regulatory effects across human tissues.  Science, 369(6509):1318–1330, 2020.
  * [2] Urmo Võsa, Annique Claringbould, Harm-Jan Westra, Marc Jan Bonder, Patrick Deelen, Biao Zeng, Holger Kirsten, Ashis Saha, Roman Kreuzhuber, Seyhan Yazar, et al.  Large-scale cis-and trans-eqtl analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression.  Nature genetics, 53(9):1300–1310, 2021.
  * [3] Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles P Fulco, Livnat Jerby-Arnon, Nemanja D Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, et al.  Perturb-seq: dissecting molecular circuits with scalable single-cell rna profiling of pooled genetic screens.  cell, 167(7):1853–1866, 2016.
  * [4] Thomas M Norman, Max A Horlbeck, Joseph M Replogle, Alex Y Ge, Albert Xu, Marco Jost, Luke A Gilbert, and Jonathan S Weissman.  Exploring genetic interaction manifolds constructed from rich single-cell phenotypes.  Science, 365(6455):786–793, 2019.
  * [5] Joseph M Replogle, Reuben A Saunders, Angela N Pogson, Jeffrey A Hussmann, Alexander Lenail, Alina Guna, Lauren Mascibroda, Eric J Wagner, Karen Adelman, Gila Lithwick-Yanai, et al.  Mapping information-rich genotype-phenotype landscapes with genome-scale perturb-seq.  Cell, 185(14):2559–2575, 2022.
  * [6] Charles Boone, Howard Bussey, and Brenda J Andrews.  Exploring genetic interactions and networks with yeast.  Nature Reviews Genetics, 8(6):437–449, 2007.
  * [7] Michael Costanzo, Anastasia Baryshnikova, Jeremy Bellay, Yungil Kim, Eric D Spear, Carolyn S Sevier, Huiming Ding, Judice LY Koh, Kiana Toufighi, Sara Mostafavi, et al.  The genetic landscape of a cell.  science, 327(5964):425–431, 2010.
  * [8] Michael Costanzo, Benjamin VanderSluis, Elizabeth N Koch, Anastasia Baryshnikova, Carles Pons, Guihong Tan, Wen Wang, Matej Usaj, Julia Hanchard, Susan D Lee, et al.  A global genetic interaction network maps a wiring diagram of cellular function.  Science, 353(6306):aaf1420, 2016.
  * [9] Thomas Verma and Judea Pearl.  On equivalence of causal models.  In Proceedings of the Sixth Conference Annual Conference on Uncertainty in Artificial Intelligence (UAI-90), pages 220–227, 1990.
  * [10] Alain Hauser and Peter Bühlmann.  Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs.  The Journal of Machine Learning Research, 13(1):2409–2464, 2012.
  * [11] Chandler Squires and Caroline Uhler.  Causal structure learning: a combinatorial perspective.  Foundations of Computational Mathematics, pages 1–35, 2022.
  * [12] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  Advances in neural information processing systems, 31, 2018.
  * [13] Alexander Reisach, Christof Seiler, and Sebastian Weichwald.  Beware of the simulated dag! causal discovery benchmarks may be easy to game.  Advances in Neural Information Processing Systems, 34:27772–27784, 2021.
  * [14] Po-Ling Loh and Peter Bühlmann.  High-dimensional learning of linear causal networks via inverse covariance estimation.  The Journal of Machine Learning Research, 15(1):3065–3105, 2014.
  * [15] Ignavier Ng, AmirEmad Ghassami, and Kun Zhang.  On the role of sparsity and dag constraints for learning linear dags.  Advances in Neural Information Processing Systems, 33:17943–17954, 2020.
  * [16] Frederick Eberhardt and Richard Scheines.  Interventions and causal inference.  Philosophy of Science, 74(5):981–995, 2007.
  * [17] Judea Pearl.  Causality.  Cambridge university press, 2009.
  * [18] Marcus Kaiser and Maksim Sipos.  Unsuitability of notears for causal graph discovery when dealing with dimensional quantities.  Neural Processing Letters, 54(3):1587–1595, 2022.
  * [19] Karren Yang, Abigail Katcoff, and Caroline Uhler.  Characterizing and learning equivalence classes of causal dags under interventions.  In International Conference on Machine Learning, pages 5541–5550. PMLR, 2018.
  * [20] Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler.  Permutation-based causal inference algorithms with interventions.  Advances in Neural Information Processing Systems, 30, 2017.
  * [21] Chandler Squires, Yuhao Wang, and Caroline Uhler.  Permutation-based causal structure learning with unknown intervention targets.  In Conference on Uncertainty in Artificial Intelligence, pages 1039–1048. PMLR, 2020.
  * [22] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin.  Differentiable causal discovery from interventional data.  Advances in Neural Information Processing Systems, 33:21865–21877, 2020.
  * [23] Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvarinen, Yoshinobu Kawahara, Takashi Washio, Patrik O Hoyer, Kenneth Bollen, and Patrik Hoyer.  Directlingam: A direct method for learning a linear non-gaussian structural equation model.  Journal of Machine Learning Research-JMLR, 12(Apr):1225–1248, 2011.
  * [24] Paul Erdős, Alfréd Rényi, et al.  On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci, 5(1):17–60, 1960.
  * [25] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  science, 286(5439):509–512, 1999.
  * [26] Michael I Love, Wolfgang Huber, and Simon Anders.  Moderated estimation of fold change and dispersion for rna-seq data with deseq2.  Genome biology, 15(12):1–21, 2014.
  * [27] Constantin Ahlmann-Eltze and Wolfgang Huber.  glmgampoi: fitting gamma-poisson generalized linear models on single cell count data.  Bioinformatics, 36(24):5701–5702, 2020.
  * [28] Damian Szklarczyk, Rebecca Kirsch, Mikaela Koutrouli, Katerina Nastou, Farrokh Mehryary, Radja Hachilif, Annika L Gable, Tao Fang, Nadezhda T Doncheva, Sampo Pyysalo, et al.  The string database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest.  Nucleic acids research, 51(D1):D638–D646, 2023.
  * [29] Yuhao Wang, Chandler Squires, Anastasiya Belyaeva, and Caroline Uhler.  Direct estimation of differences in causal graphs.  Advances in neural information processing systems, 31, 2018.
  * [30] Johannes Köster and Sven Rahmann.  Snakemake—a scalable bioinformatics workflow engine.  Bioinformatics, 28(19):2520–2522, 2012.
  * [31] Peter Bühlmann, Jonas Peters, and Jan Ernest.  Cam: Causal additive models, high-dimensional order search and penalized regression.  2014\.
  * [32] David Maxwell Chickering.  Optimal structure identification with greedy search.  Journal of machine learning research, 3(Nov):507–554, 2002.

##  7 Supplementary Material

###  7.1 Least squares infers varsortable structures in expectation

We first derive an explicit cutoff for varsortability, matching the cutoff
given in Eq. 9.

######  Lemma 1.

Let $\gamma\coloneqq\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}$. The system
$X_{1}\overset{w}{\rightarrow}X_{2}$ is varsortable if and only if
$|w|\geq\sqrt{1-\frac{1}{\gamma}}$.

For the proof, see Supplementary Material 7.1.1. We verify that the least
squares loss infers varsortable structures by comparing the SEM in Eq. 8
against the false but Markov equivalent model
$X_{1}\overset{\delta}{\leftarrow}X_{2}$, with weighted adjacency matrix
$W_{\delta}\coloneqq\begin{pmatrix}0&0\\\ \delta&0\end{pmatrix}$ and SEM

$\displaystyle X_{1}$ $\displaystyle=\delta X_{2}+\epsilon_{1},$ (11)
$\displaystyle X_{2}$ $\displaystyle=\epsilon_{2}.$

Let $\left(\cdot\right)_{i}$ represent the $i$th column vector, and denote the
least squares loss as

$\ell\left(W,\boldsymbol{X}\right)\coloneqq\frac{1}{2n_{0}}\sum_{i=1}^{p}\left|\left|\left(\boldsymbol{X}-\boldsymbol{X}W\right)_{i}\right|\right|_{F}^{2}.$
(12)

######  Theorem 1.

Let $\gamma\coloneqq\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}$.
$\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)\leq\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)$
for all $\delta$ if and only if $|w|\geq\sqrt{1-\frac{1}{\gamma}}$.

######  Corollary 1.

In the system $X_{1}\overset{w}{\rightarrow}X_{2}$, the least squares loss is
uniquely minimized in expectation by the true DAG if and only if the system is
varsortable.

For the proof of Theorem 1, see Supplementary Material 7.1.2.

####  7.1.1 Proof of Lemma 1

###### Proof.

Under the generative SEM in Eq. 8, we have

$\begin{aligned} \text{Var}\left(X_{1}\right)&=\sigma_{1}^{2}\\\
\text{Var}\left(X_{2}\right)&=w^{2}\sigma_{1}^{2}+\sigma_{2}^{2}\end{aligned}.$
(13)

The system is therefore varsortable if and only if
$w^{2}\sigma_{1}^{2}+\sigma_{2}^{2}\geq\sigma_{1}^{2}$. Using the substitution
$\sigma_{1}^{2}=\gamma\sigma_{2}^{2}$, we obtain

$\displaystyle w^{2}\sigma_{1}^{2}+\sigma_{2}^{2}$
$\displaystyle\geq\sigma_{1}^{2}$ (14) $\displaystyle
w^{2}\gamma\sigma_{2}^{2}+\sigma_{2}^{2}$
$\displaystyle\geq\gamma\sigma_{2}^{2}$ $\displaystyle w^{2}\gamma+1$
$\displaystyle\geq\gamma$ $\displaystyle|w|$
$\displaystyle\geq\sqrt{1-\frac{1}{\gamma}}$

∎

####  7.1.2 Proof of Theorem 1

###### Proof.

From Eq. 12, we can calculate both
$\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)$ and
$\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)$ component-wise. Under
$W_{0}$ we obtain component-wise

$\displaystyle\left({X}-{X}W_{0}\right)_{1}$ $\displaystyle={X_{1}}$ (15)
$\displaystyle={\epsilon_{1}}$ by Eq. 8
$\displaystyle\left({X}-{X}W_{0}\right)_{2}$ $\displaystyle={X_{2}}-w{X_{1}}$
$\displaystyle={\epsilon_{2}}$ by Eq. 8

As a result,

$\displaystyle\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)$
$\displaystyle=\frac{1}{2n_{0}}\sum_{i=1}^{p}\mathbb{E}\left|\left|\left({\boldsymbol{X}}-{\boldsymbol{X}}W_{0}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle=\frac{1}{2n_{0}}\left(\mathbb{E}\left(\left|\left|{\boldsymbol{\epsilon}_{1}}\right|\right|_{F}^{2}\right)+\mathbb{E}\left(\left|\left|\boldsymbol{\epsilon}_{2}\right|\right|_{F}^{2}\right)\right)$
$\displaystyle=\frac{1}{2}\left(\sigma_{1}^{2}+\sigma_{2}^{2}\right)$

Here we note that the expected loss under the generative DAG is
$\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)=\frac{1}{2}\sum_{i=1}^{p}\text{Var}\left({\epsilon_{i}}\right)$.

Similarly, under $W_{\delta}$ we have

$\displaystyle\left(X-XW_{\delta}\right)_{1}$ $\displaystyle=X_{1}-\delta
X_{2}$ (16)
$\displaystyle=\epsilon_{1}-\delta\left(w\epsilon_{1}+\epsilon_{2}\right)$ by
Eq. 8 $\displaystyle=(1-\delta w)\epsilon_{1}-\delta\epsilon_{2}$
$\displaystyle\left(X-XW_{\delta}\right)_{2}$ $\displaystyle=X_{2}$
$\displaystyle=w\epsilon_{1}+\epsilon_{2}$ by Eq. 8

and therefore $\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)$ is

$\frac{1}{2}\left(\sigma_{1}^{2}\left(w^{2}+\left(1-\delta
w\right)^{2}\right)+\sigma_{2}^{2}\left(1+\delta^{2}\right)\right).$

We can now ask when
$\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)\leq\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)$,
or equivalently,

$\displaystyle\sigma_{1}^{2}\left(w^{2}+\left(1-\delta
w\right)^{2}\right)+\sigma_{2}^{2}\left(1+\delta^{2}\right)$
$\displaystyle\leq\sigma_{1}^{2}+\sigma_{2}^{2}$ (17)
$\displaystyle\gamma\sigma_{2}^{2}\left(w^{2}+\left(1-\delta
w\right)^{2}\right)+\sigma_{2}^{2}\left(1+\delta^{2}\right)$
$\displaystyle\leq\gamma\sigma_{2}^{2}+\sigma_{2}^{2}$
$\displaystyle\gamma\left(w^{2}+1-2\delta
w+\delta^{2}w^{2}\right)+1+\delta^{2}$ $\displaystyle\leq\gamma+1$
$\displaystyle\gamma\left(w^{2}-2\delta w+\delta^{2}w^{2}\right)+\delta^{2}$
$\displaystyle\leq 0$ $\displaystyle\delta^{2}\left(1+\gamma
w^{2}\right)-\delta\left(2\gamma w\right)+\gamma w^{2}$ $\displaystyle\leq 0$

The inequality becomes a quadratic in $\delta$, which when solved gives us
roots at

$\displaystyle\delta$ $\displaystyle=\frac{2\gamma w\pm\sqrt{(-2\gamma
w)^{2}-4\left(1+\gamma w^{2}\right)\left(\gamma w^{2}\right)}}{2\left(1+\gamma
w^{2}\right)}$ (18) $\displaystyle=\frac{2\gamma w\pm 2\gamma
w\sqrt{1-\frac{1}{\gamma}\left(1+\gamma w^{2}\right)}}{2\left(1+\gamma
w^{2}\right)}$ $\displaystyle=\frac{\gamma
w\left(1\pm\sqrt{1-w^{2}-\frac{1}{\gamma}}\right)}{1+\gamma w^{2}}.$

Note that for the quadratic in Eq. 17, $\gamma w^{2}\geq 0$. Then
$\exists\delta$ such that
$\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)\leq\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)$
when the term in the root exists in Eq. 18. Equivalently, we can say that
$\mathbb{E}\ell\left(W_{0},\boldsymbol{X}\right)\leq\mathbb{E}\ell\left(W_{\delta},\boldsymbol{X}\right)$
is guaranteed for all $\delta$ if and only if

$|w|\geq\sqrt{1-\frac{1}{\gamma}}.$ (19)

∎

###  7.2 Proof of Consistency

For a diagonal matrix $\Omega$, define

$\displaystyle\mathcal{L}_{\Omega}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)$
$\displaystyle\coloneqq\frac{1}{n_{k}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)\Omega^{-\frac{1}{2}}\right|\right|_{F}^{2}$
(20)
$\displaystyle=\frac{1}{n_{k}}\sum_{j=1}^{p}\frac{1}{\Omega_{jj}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{j}\right|\right|_{F}^{2},$

where $\Omega_{jj}$ is the $j,j$th entry of $\Omega$. Further, recall the
definition

$\Omega_{0}^{(k)}\coloneqq\text{Cov}\left(\epsilon^{(k)}\right).$

Denote convergence in probability by $\overset{p}{\rightarrow}.$ We show that
under a sub-Gaussian assumption, the dotears loss is a consistent estimator of
the true DAG.

######  Assumption 4.

For any $k$, we assume $X^{(k)}$ is a sub-Gaussian random vector with
parameter $\sigma^{2}$.

######  Lemma 2.

For any intervention $k$,

$\mathcal{L}_{{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=\mathcal{L}_{\Omega_{0}^{(k)}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)+\mathcal{O}_{W}(1)$

For the proof, see Supplementary Material 7.2.1.

######  Corollary 2.

$\arg\min_{W}\mathcal{L}_{\Omega_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=\arg\min_{W}\mathcal{L}_{\Omega_{0}^{(k)}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)$

For any intervention $k$ and any $\alpha$, Lemma 2 and Corollary 2 show that
true exogenous variance structure
$\Omega_{0}=\alpha^{2}\mathbb{E}\hat{\Omega}_{0}$ is sufficient for structure
recovery. Results from Loh and Buhlmann (2014) are then sufficient to
establish consistency of the estimator
$\operatorname*{arg\,min}_{W}\mathcal{L}_{\Omega_{0}}\left(W^{(k)},X^{(k)}\right)$.
We wish to show consistency of
$\operatorname*{arg\,min}_{W}\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)$,
where we have an estimated $\hat{\Omega}_{0}$ rather than $\Omega_{0}$.

For simplicity, in the following we assume $\alpha=1$ without loss of
generality, which is justified by Corollary 2 and Lemma 2. Note then that
$\Omega_{0}=\Omega_{0}^{(k)}$, which allows us to drop the (k) notation.

######  Remark 1.

For consistency under an estimated $\hat{\Omega}_{0}$, we have two cases: the
observational case $k=0$ and the interventional case $k\neq 0$. The
observational case is simplest - in particular, when $k=0$, $\hat{\Omega}_{0}$
is independent of $\boldsymbol{X}^{(0)}$, and we may therefore freely estimate
$W_{0}^{(0)}$ from
$\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(0)},\boldsymbol{X}^{(0)}\right)$.
When $k\neq 0$, we can assume independence of $\hat{\Omega}_{0}$ and
$\boldsymbol{X}^{(k)}$ under a data splitting framework, where some fraction
of our $n_{k}$ samples are reserved for estimation of $\hat{\sigma}_{k}^{2}$,
and the remaining samples are given in $\boldsymbol{X}^{(k)}$.

######  Assumption 5.

Let

$\hat{\sigma}_{k}^{2}\coloneqq\left(\hat{\Omega}_{0}\right)_{kk}.$

Then

$\hat{\sigma}_{k}^{2}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{X}^{(k)}_{k}.$

However, in the presented applications we use the full $n_{k}$ samples of
$\boldsymbol{X}_{k}^{(k)}$ for both estimates.

Note that we abuse notation to write $\boldsymbol{X}^{(k)}$ as the set of
samples used to estimate $\hat{W}_{0}^{(k)}$ given $\hat{\Omega}_{0}$.
Further, for convenience we abuse notation to write $n_{k}$ as the sample size
of $\boldsymbol{X}^{(k)}$ AND the sample size of $\hat{\sigma}_{k}^{2}$, since
the respective sample sizes are equal to $n_{k}$ up to a multiplicative
constant which can be ignored as $n_{k}\to\infty$.

We further define

$\mathbb{E}\mathcal{L}_{\Omega}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)\\\
\coloneqq\frac{1}{n_{k}}\sum_{j=1}^{p}\frac{1}{\Omega_{jj}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{j}\right|\right|_{F}^{2}.$

Here, we abuse notation to remove $\Omega_{jj}$ from the expectation.
Moreover, for the random vector $X$ we write

$\mathbb{E}\mathcal{L}_{\Omega}\left(W^{(k)},X^{(k)}\right)\\\
\coloneqq\sum_{j=1}^{p}\frac{1}{\Omega_{jj}}\mathbb{E}\left[\left(X^{(k)}-X^{(k)}W^{(k)}\right)_{j}^{2}\right],$

and note that the two definitions are equivalent under expectation. Note
further that

$\text{score}_{\Omega}(W)\equiv\mathbb{E}\mathcal{L}_{\Omega}(W,X),$

where score is defined in [14].

######  Remark 2.

Loh and Buhlmann differentiate between the weighted adjacency matrix $W$ and
the binary DAG $G$. In particular, they define

$\text{score}_{\Omega}\left(G,X\right)\coloneqq\min_{W\in\mathcal{U}_{G}}\left\\{\text{score}_{\Omega}\left(W,X\right)\right\\},$

where $\mathcal{U}_{G}$ is the set of weighted adjacency matrices with support
in $G$. However, for the true binary DAG $G_{0}$,

$\text{score}_{\Omega}\left(G,X\right)=\text{score}_{\Omega}\left(G_{0},X\right)\quad\forall
G\supseteq G_{0},$

which leads to identifiability issues (see Lemma 6 and discussion of Lemma 19
[14]).

We define everything directly on the weighted adjacency matrix $W$, and ignore
the binary DAG $G$. This discrepancy is essentially due to methodological
differences. Loh and Buhlmann first find a superset of the binary structure
$G_{0}$, and subsequently rely on sparse regression to recover edge weights
[14]. However, dotears searches directly in the space of weighted adjacency
matrices to find $W_{0}$. While the end result is similar, we use this remark
to explain discrepancies in notation and proof structure for readers following
with [14].

We start by restating results from [14]. Let $\Omega_{1}\in\mathbb{R}^{p\times
p}$ an arbitrary diagonal weight matrix. Define

$a_{\text{max}}\coloneqq\lambda_{\text{max}}\left(\Omega_{0}\Omega_{1}^{-1}\right),\,\,a_{\text{min}}\coloneqq\lambda_{\text{min}}\left(\Omega_{0}\Omega_{1}^{-1}\right)$

to be the maximum and minimum ratios between the corresponding diagonal
entries of $\Omega_{0}$ and $\Omega_{1}$, and further define

$\displaystyle\xi^{(k)}_{\Omega}\coloneqq\min_{\begin{subarray}{c}W\in\mathcal{D}\\\
W\neq
W_{0}\end{subarray}}\mathbb{E}\mathcal{L}_{\Omega}\left(W^{(k)},X^{(k)}\right)-\mathbb{E}\mathcal{L}_{\Omega}\left(W_{0}^{(k)},X^{(k)}\right)$
(21)

to be the additive gap between the expected loss of the true weighted
adjacency matrix $W_{0}$ and the expected loss of the next-best weighted
adjacency matrix under an arbitrary diagonal matrix $\Omega$. Note that
$\xi_{\Omega_{0}}>0$ (see Theorem 7 in [14]). Then we have the following
theorem, whose proof is given in [14]:

######  Theorem 2.

Suppose

$\frac{a_{\text{max}}}{a_{\text{min}}}\leq
1+\frac{\xi_{\Omega_{0}}^{(k)}}{p}.$ (22)

Then
$W_{0}^{(k)}\in\operatorname*{arg\,min}_{W\in\mathcal{D}}\left\\{\mathbb{E}\mathcal{L}_{\Omega_{1}}\left(W,X^{(k)}\right)\right\\}$.
If Inequality 22 is strict, then $W_{0}^{(k)}$ is the unique minimizer of
$\mathbb{E}\mathcal{L}_{\Omega_{1}}\left(W,X^{(k)}\right)$.

######  Theorem 3.

For any intervention $k$, as $n_{j}\to\infty$ for all $j=1,\dots,p$, with high
probability $W_{0}^{(k)}$ is the unique minimizer of
$\mathbb{E}\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},X^{(k)}\right)$, and
thus $\xi_{\hat{\Omega}_{0}}^{(k)}>0$.

For the proof, see Supplementary Material 7.2.2.

######  Lemma 3.

For any intervention $k$, suppose $X^{(k)}$ is a sub-Gaussian random vector
with parameter $\sigma^{2}$. Then $X^{(k)}_{k}$ is a sub-Gaussian random
variable with parameter $\sigma^{2}$, and for all $t\geq 0$, we have

$\left|\frac{1}{\hat{\sigma}_{k}^{2}}-\frac{1}{\sigma_{k}^{2}}\right|\leq\frac{2\sigma^{2}}{\sigma_{k}^{2}}\cdot\max\\{\delta,\delta^{2}\\}$

with probability at least $1-2\exp(-cn_{k}t^{2})$, where
$\delta=c^{\prime}\sqrt{\frac{1}{n_{k}}}+c^{\prime\prime}t$ and

$\frac{\sigma^{2}}{\sigma_{k}^{2}}\cdot\max\\{\delta,\delta^{2}\\}\leq\frac{1}{2}.$
(23)

###### Proof.

Lemma 23 is a re-statement of Lemma 29 in [14] (Appendix E) for the degenerate
case $p=1$. ∎

Following [14], we introduce the new notation $f^{(k)}_{\sigma_{j}}$. For a
node $j$ and a set $S\subseteq\\{1\dots p\\}\backslash\\{j\\}$, define

$f^{(k)}_{\sigma_{j}}(S)\coloneqq\frac{1}{\sigma_{j}^{2}}\mathbb{E}\left[\left(X^{(k)}_{j}-X^{(k)}_{S}w_{j}\right)^{2}\right],$

where $X^{(k)}_{j}$ is the $j$th column vector of $X^{(k)}$, and
$X^{(k)}_{S}w_{j}$ is the best linear predictor for $X^{(k)}_{j}$ regressed
upon $X_{S}^{(k)}$. Similarly,

$\hat{f}_{\sigma_{j}}^{(k)}(S)\coloneqq\frac{1}{\sigma_{j}^{2}}\frac{1}{n_{k}}\left|\left|\boldsymbol{X}^{(k)}_{j}-\boldsymbol{X}^{(k)}_{S}\hat{w}_{j}\right|\right|_{F}^{2},$

where $\hat{w}_{j}$ is the ordinary least squares solution for linear
regression of $\boldsymbol{X}^{(k)}_{j}$ upon $\boldsymbol{X}^{(k)}_{S}$, i.e.

$\displaystyle\hat{w}_{j}$
$\displaystyle\coloneqq\left(\left(\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\boldsymbol{X}^{(k)}_{S}\right)^{-1}\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\boldsymbol{X}^{(k)}_{j}\right)$
(24)
$\displaystyle=\left(\left(\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\boldsymbol{X}^{(k)}_{S}\right)^{-1}\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\left(\boldsymbol{X}^{(k)}_{S}w_{j}+\boldsymbol{e}^{(k)}_{j}\right)\right)$
$\displaystyle=w_{j}+\left(\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\boldsymbol{X}^{(k)}_{S}\right)^{-1}\left(\boldsymbol{X}^{(k)}_{S}\right)^{\textrm{T}}\boldsymbol{e}^{(k)}_{j},$

where note that

$\boldsymbol{X}_{j}^{(k)}=\boldsymbol{X}_{S}^{(k)}w_{j}+\boldsymbol{e}^{(k)}_{j}.$

Lastly, denote the vector $\ell_{2}$ norm as $\|\|\cdot\|\|_{2}$.

Loh and Buhlmann achieve a high-dimensional consistency result by first
conditioning on the support of a sparse precision matrix $\Theta^{(k)}$, since
the support of $\Theta^{(k)}\coloneqq\text{Cov}\left(X^{(k)}\right)^{-1}$
defines the moralized graph, an undirected graph obtained from a DAG by adding
edges between all parents with a shared child node [14, 29]. Since the edge
set of the moralized graph is a superset of the edge set of the true DAG, for
a given node $j$ one may condition on the maximum size of the putative
neighbor set, $N_{\Theta^{(k)}}(j)$. Indeed, in [14] Loh and Buhlmann restrict
$|N_{\Theta^{(k)}}(j)|\leq d$ for all $j$, with the only restriction on $d$
being that $d\leq n$.

dotears does not condition on $\Theta^{(k)}$, and therefore does not condition
on the moralized graph or $N_{\Theta}(j)$. To maintain the validity of our
consistency proof, we let $d=p-1\leq n$. Under this restriction, we no longer
have a high-dimensional result, but maintain consistency of the loss function
of dotears for low dimensionality of $p$.

######  Assumption 6.

For all $k=0,\dots,p$, let $d=p-1\leq n_{k}$.

######  Lemma 4.

For any intervention $k$, suppose $X^{(k)}$ is sub-Gaussian with parameter
$\sigma^{2}$. Then $\forall j\text{ and }S\subseteq\\{1\dots
p\\}\backslash\\{j\\}$,

$\sigma_{j}^{2}\left|\hat{f}^{(k)}_{\hat{\sigma}_{j}}(S)-\hat{f}^{(k)}_{\sigma_{j}}(S)\right|\leq\sigma^{2}\cdot\max\\{\delta,\delta^{2}\\}\cdot\frac{C}{n_{k}}\left|\left|\boldsymbol{e}^{(k)}_{j}\right|\right|_{2}^{2}$
(25)

with probability $1-c_{1}\exp(-c_{2}\log p)$, where
$\delta=c^{\prime}\sqrt{\frac{1}{n_{j}}}+c^{\prime\prime}\sqrt{\frac{\log
p}{n_{j}}}$, if

$\frac{\sigma^{2}}{\sigma_{j}^{2}}\cdot\delta\leq\frac{1}{2}.$ (26)

For the proof, see Supplementary Material 7.2.3.

######  Lemma 5.

For any intervention $k$, suppose $X^{(k)}$ is sub-Gaussian with parameter
$\sigma^{2}$. Then as $n_{i}\to\infty$ for all $i=1\dots p$, $\forall j\text{
and }S\subseteq\\{1\dots p\\}\backslash\\{j\\}$ it is true that

$\left|\hat{f}^{(k)}_{\hat{\sigma}_{j}}(S)-f^{(k)}_{\sigma_{j}}(S)\right|\leq
c_{0}\sigma^{4}\sqrt{\frac{\log p}{n_{j}}}+c_{1}\sigma^{2}\sqrt{\frac{\log
p}{n_{k}}}+c_{2}\frac{p}{n_{k}}$ (27)

with probability at least $1-c_{1}\exp(-c_{2}\log p)$.

For the proof, see Supplementary Material 7.2.4.

######  Theorem 4.

For any intervention $k$, suppose Inequality 27 holds, and suppose

$\left(c_{0}\sigma^{4}\sqrt{\frac{\log
p}{n_{j}}}+c_{1}\sigma^{2}\sqrt{\frac{\log
p}{n_{k}}}+c_{2}\frac{p}{n_{k}}\right)\cdot\sum_{j=1}^{p}\frac{1}{\sigma_{j}^{2}}<\frac{\xi_{{\Omega}_{0}}^{(k)}}{2}.$
(28)

Then

$\mathcal{L}_{\hat{\Omega}_{0}}\left(W_{0}^{(k)},\boldsymbol{X}^{(k)}\right)<\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)$

$\forall W\in\mathcal{D},$ and the estimator

$\hat{W}^{(k)}\coloneqq\operatorname*{arg\,min}_{W}\mathcal{L}_{\hat{\Omega}_{0}}\left(W,\boldsymbol{X}^{(k)}\right)$
(29)

is therefore consistent as $n_{j}\to\infty$ for all $j=1\dots p$.

For the proof, see Supplementary Material 7.2.5.

####  7.2.1 Proof of Lemma 2

###### Proof.

Note that

$\displaystyle\mathcal{L}_{\Omega_{0}^{(k)}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)\coloneqq$
$\displaystyle\frac{1}{n_{k}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)\left(\Omega_{0}^{(k)}\right)^{-\frac{1}{2}}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\sum_{i=1}^{p}\frac{1}{\text{Var}\left(\epsilon_{i}^{(k)}\right)}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\frac{\alpha^{2}}{\sigma_{k}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{k}\right|\right|_{F}^{2}+\frac{1}{n_{k}}\sum_{\begin{subarray}{c}i=1\\\
i\neq
k\end{subarray}}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle\mathcal{L}_{\Omega_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)\coloneqq$
$\displaystyle\frac{1}{n_{k}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)\Omega_{0}^{-\frac{1}{2}}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\sum_{i=1}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}.$

Let $I_{0}^{k}=\text{diag}(1,1,\dots,0,\dots,1)$ a modifier on the $p\times p$
identity matrix, where the diagonal entries are all $1$ except for the $k,k$
entry, which is $0$.

Then in the interventional system,

$W^{(k)}=WI_{0}^{k}$

and thus

$\mathcal{L}_{\Omega_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=\frac{1}{n_{k}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}WI_{0}^{k}\right)\Omega_{0}^{-\frac{1}{2}}\right|\right|_{F}^{2}$

We expand $\mathcal{L}_{\Omega_{0}}$ to obtain

$\displaystyle\mathcal{L}_{\Omega_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=$
$\displaystyle\frac{1}{n_{k}}\sum_{i=1}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\frac{1}{\sigma_{k}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{k}\right|\right|_{F}^{2}+\frac{1}{n_{k}}\sum_{\begin{subarray}{c}i=1\\\
i\neq
k\end{subarray}}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$

and note that

$\displaystyle\left(\boldsymbol{X}^{(k)}W^{(k)}\right)_{k}$
$\displaystyle=\boldsymbol{X}^{(k)}\left(WI_{0}^{k}\right)_{k}$
$\displaystyle=\boldsymbol{X}^{(k)}\vec{0}_{p}$
$\displaystyle=\vec{0}_{n_{k}}.$

As a result,

$\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{k}\right|\right|_{F}^{2}=\left|\left|\boldsymbol{X}_{k}^{(k)}\right|\right|_{F}^{2}=\mathcal{O}_{W}(1)$

is constant in $W$, and

$\displaystyle\mathcal{L}_{\Omega_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=$
$\displaystyle\frac{1}{n_{k}}\frac{1}{\sigma_{k}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{k}\right|\right|_{F}^{2}+\frac{1}{n_{k}}\sum_{\begin{subarray}{c}i=1\\\
i\neq
k\end{subarray}}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\frac{1}{\sigma_{k}^{2}}\left|\left|\boldsymbol{X}_{k}^{(k)}\right|\right|_{F}^{2}+\frac{1}{n_{k}}\sum_{\begin{subarray}{c}i=1\\\
i\neq
k\end{subarray}}^{p}\frac{1}{\sigma_{i}^{2}}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\mathcal{L}_{\Omega_{0}^{(k)}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)+\mathcal{O}_{W}(1)$

∎

####  7.2.2 Proof of Theorem 3

###### Proof.

Let
$a_{\text{max}}\coloneqq\lambda_{\text{max}}\left(\Omega_{0}\hat{\Omega}_{0}^{-1}\right)$,
and
$a_{\text{min}}\coloneqq\lambda_{\text{min}}\left(\Omega_{0}\hat{\Omega}_{0}^{-1}\right)$
similarly. We first prove that
$a_{\text{max}}\overset{p}{\rightarrow}\alpha^{2}$ without loss of generality.
For all $k$, note first that

$\frac{\sigma_{k}^{2}}{\hat{\sigma}_{k}^{2}}\overset{p}{\rightarrow}\alpha^{2}.$

Then by Continuous Mapping Theorem we have

$a_{\text{max}}\overset{p}{\rightarrow}{\alpha^{2}}.$

Similarly, $a_{\text{min}}\overset{p}{\rightarrow}{\alpha^{2}}$, and therefore

$\frac{a_{\text{max}}}{a_{\text{min}}}\overset{p}{\rightarrow}1.$

We note that $\frac{a_{\text{max}}}{a_{\text{min}}}=1$ is almost impossible
empirically, but can be controlled with high probability to $1$ with arbitrary
precision. Then by Theorem 2, with high probability $W_{0}^{(k)}$ is the
unique minimizer of
$\mathbb{E}\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},X^{(k)}\right)$ as
$n_{j}\to\infty$ for all $j=1\dots p$. Further,
$\xi_{\hat{\Omega}_{0}}^{(k)}>0$. ∎

####  7.2.3 Proof of Lemma 26

###### Proof.

Let the projection matrix of $\boldsymbol{X}_{S}^{(k)}$ be defined as
$P_{\boldsymbol{X}_{S}^{(k)}}\coloneqq\boldsymbol{X}_{S}^{(k)}\left(\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}\boldsymbol{X}_{S}^{(k)}\right)^{-1}\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}$.

Given the least squares estimate in Eq. 24, we may write

$\displaystyle\sigma_{j}^{2}\cdot\hat{f}^{(k)}_{\sigma_{j}}(S)$
$\displaystyle=\frac{1}{n_{k}}\left|\left|\boldsymbol{X}^{(k)}_{j}-\boldsymbol{X}_{S}^{(k)}\hat{w}_{j}\right|\right|_{F}^{2}$
(30)
$\displaystyle=\frac{1}{n_{k}}\left|\left|\boldsymbol{X}_{j}^{(k)}\left(w_{j}-\hat{w}_{j}\right)+\boldsymbol{e}^{(k)}_{j}\right|\right|_{F}^{2}$
$\displaystyle=\frac{1}{n_{k}}\left|\left|\left(I-P_{\boldsymbol{X}_{S}^{(k)}}\right)\boldsymbol{e}_{j}^{(k)}\right|\right|_{F}^{2}.$

Note by Triangle Inequality, we have for the $\ell_{2}$ norm
$\left|\left|\cdot\right|\right|_{2}$ and the spectral norm
$\left|\left|\left|\cdot\right|\right|\right|_{2}$

$\displaystyle\left|\left|\left|\left(I-P_{\boldsymbol{X}_{S}^{(k)}}\right)\boldsymbol{e}_{j}^{(k)}\right|\right|_{F}-\left|\left|\boldsymbol{e}_{j}^{(k)}\right|\right|_{2}\right|\leq$
$\displaystyle\left|\left|P_{\boldsymbol{X}_{S}^{(k)}}\boldsymbol{e}_{j}^{(k)}\right|\right|_{F}$
$\displaystyle\leq$
$\displaystyle\left|\left|\left|P_{\boldsymbol{X}_{S}^{(k)}}\right|\right|\right|_{2}\left|\left|\boldsymbol{e}_{j}^{(k)}\right|\right|_{2}$
$\displaystyle\leq$
$\displaystyle\left|\left|\boldsymbol{e}_{j}^{(k)}\right|\right|_{2},$

where the spectral norm of the projection matrix
$\left|\left|\left|P_{\boldsymbol{X}_{S}^{(k)}}\right|\right|\right|_{2}$ is
1. Then

$\sigma_{j}^{2}\cdot\hat{f}^{(k)}_{\sigma_{j}}(S)=\frac{1}{n_{k}}\left|\left|\boldsymbol{X}^{(k)}_{j}-\boldsymbol{X}_{S}^{(k)}\hat{w}_{j}\right|\right|_{F}^{2}\leq\frac{2}{n_{k}}\left|\left|\boldsymbol{e}_{j}^{(k)}\right|\right|_{2}^{2}.$

Expanding the left hand side of Eq. 25, we obtain

$\displaystyle\sigma_{j}^{2}\left|\hat{f}^{(k)}_{\hat{\sigma}_{j}}(S)-\hat{f}^{(k)}_{\sigma_{j}}(S)\right|=$
$\displaystyle\sigma_{j}^{2}\left|\frac{1}{\hat{\sigma}_{j}^{2}}-\frac{1}{\sigma_{j}^{2}}\right|\frac{1}{n_{k}}\left|\left|\boldsymbol{X}_{j}^{(k)}-\boldsymbol{X}_{S}^{(k)}\hat{w}_{j}\right|\right|_{F}^{2}$
$\displaystyle\leq$ $\displaystyle\sigma_{j}^{2}\cdot
2\frac{\sigma^{2}}{\sigma_{j}^{2}}\cdot\max\\{\delta,\delta^{2}\\}\cdot\frac{1}{n_{k}}\left|\left|\boldsymbol{X}_{j}^{(k)}-\boldsymbol{X}_{S}^{(k)}\hat{w}_{j}\right|\right|_{F}^{2}$
$\displaystyle\leq$ $\displaystyle
2\sigma^{2}\cdot\max\\{\delta,\delta^{2}\\}\cdot\frac{2}{n_{k}}\left|\left|\boldsymbol{e}^{(k)}_{j}\right|\right|_{2}^{2}$

with probability $\min\\{1-c_{1}\exp(-c_{2}\log p),1-2\exp(-cn_{j}t^{2})\\}$,
where $\delta=c^{\prime}\sqrt{\frac{1}{n_{j}}}+c^{\prime\prime}t$, if
$\frac{\sigma^{2}}{\sigma_{j}^{2}}\cdot\max\\{\delta,\delta^{2}\\}\leq\frac{1}{2}$.
The first inequality is by Lemma 23. Set $t=\sqrt{\frac{\log p}{n_{j}}}$. Then

$\displaystyle\delta$
$\displaystyle=c^{\prime}\sqrt{\frac{1}{n_{j}}}+c^{\prime\prime}\sqrt{\frac{\log
p}{n_{j}}},$ $\displaystyle\delta^{2}$
$\displaystyle=(c^{\prime})^{2}\frac{1}{n_{j}}+(c^{\prime\prime})^{2}\frac{\log
p}{n_{j}}+2c^{\prime}c^{\prime\prime}\frac{\sqrt{\log p}}{n_{j}}.$

By Assumption 6, $\log p<p-1<n_{j}$. Then
$\max\\{\delta,\delta^{2}\\}=\delta$, and
$1-2\exp(-cn_{j}t^{2})=1-2\exp(-c\log p)$. ∎

####  7.2.4 Proof of Lemma 5

###### Proof.

Give the singular value decomposition on the $n_{k}$ by $p$ matrix
$\boldsymbol{X}^{(k)}_{S}$ as

$\boldsymbol{X}^{(k)}_{S}=U\Sigma V^{\textrm{T}},$

where $U$ is an $n_{k}\times p$ matrix, $\Sigma$ is a diagonal $p\times p$
matrix, $V$ is a $p\times p$ matrix, and $U^{\textrm{T}}U=V^{\textrm{T}}V=I$.
Then

$\displaystyle\boldsymbol{X}_{S}^{(k)}\left(\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}\boldsymbol{X}_{S}^{(k)}\right)^{-1}\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}=$
$\displaystyle U\Sigma
V^{\textrm{T}}\left(V\Sigma^{2}V^{\textrm{T}}\right)^{-1}V\Sigma
U^{\textrm{T}}$ $\displaystyle=$ $\displaystyle UU^{\textrm{T}}.$

We substitute into the expansion in Equation 30 to obtain

$\displaystyle\sigma_{j}^{2}\cdot\hat{f}^{(k)}_{\sigma_{j}}(S)=$
$\displaystyle\frac{1}{n_{k}}\left|\left|\left(I-\boldsymbol{X}_{S}^{(k)}\left(\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}\boldsymbol{X}_{S}^{(k)}\right)^{-1}\left(\boldsymbol{X}_{S}^{(k)}\right)^{\textrm{T}}\right)\boldsymbol{e}_{j}^{(k)}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\left|\left|\left(I-UU^{\textrm{T}}\right)\boldsymbol{e}_{j}^{(k)}\right|\right|_{F}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\textrm{tr}\left[\left(\boldsymbol{\epsilon}_{j}^{(k)}\right)^{\textrm{T}}\left(I-UU^{\textrm{T}}\right)\left(I-UU^{\textrm{T}}\right)\left(\boldsymbol{e}_{j}^{(k)}\right)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\textrm{tr}\left[\left(\boldsymbol{e}_{j}^{(k)}\right)^{\textrm{T}}\left(I-UU^{\textrm{T}}\right)\left(\boldsymbol{e}_{j}^{(k)}\right)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n_{k}}\textrm{tr}\left[\left(\boldsymbol{e}_{j}^{(k)}\right)^{\textrm{T}}\left(\boldsymbol{e}_{j}^{(k)}\right)-\left(U^{\textrm{T}}\left(\boldsymbol{e}_{j}^{(k)}\right)\right)^{\textrm{T}}U^{\textrm{T}}\left(\boldsymbol{e}_{j}^{(k)}\right)\right].$

Let

$Y\coloneqq U^{\textrm{T}}\boldsymbol{e}_{j}^{(k)}\in\mathbb{R}^{p}.$

Then

$\displaystyle\sigma_{j}^{2}\cdot\hat{f}^{(k)}_{\sigma_{j}}(S)$
$\displaystyle=\frac{1}{n_{k}}\textrm{tr}\left[\left(\boldsymbol{e}_{j}^{(k)}\right)^{\textrm{T}}\left(\boldsymbol{e}_{j}^{(k)}\right)-Y^{\textrm{T}}Y\right]$
$\displaystyle=\frac{1}{n_{k}}\left[\sum_{i=1}^{n_{k}}\left(\boldsymbol{e}_{i,j}^{(k)}\right)^{2}-\sum_{l=1}^{p}Y_{l}^{2}\right].$

Here, $Y$ is a random vector with expectation 0 and covariance

$\displaystyle\text{Cov}(Y)$ $\displaystyle=\mathbb{E}(YY^{\textrm{T}})$
$\displaystyle=\mathbb{E}\left[U^{\textrm{T}}\left(\boldsymbol{e}_{j}^{(k)}\right)\left(\boldsymbol{e}_{j}^{(k)}\right)^{\textrm{T}}U\right]$
$\displaystyle=U^{\textrm{T}}\mathbb{E}\left[\left(\boldsymbol{e}_{j}^{(k)}\right)\left(\boldsymbol{e}_{j}^{(k)}\right)^{\textrm{T}}\right]U$
$\displaystyle=U^{\textrm{T}}\left(\sigma_{j}^{2}I_{n_{k}}\right)U$
$\displaystyle=\sigma_{j}^{2}U^{\textrm{T}}U$
$\displaystyle=\sigma_{j}^{2}I_{p}.$

As a result, for fixed $p$, $\sum_{l=1}^{p}Y_{l}^{2}=\mathcal{O}_{p}(1)$, and
as $n_{k}\to\infty$ we have
$\frac{1}{n_{k}}\sum_{l=1}^{p}Y_{l}^{2}=\mathcal{O}_{p}(1)\overset{p}{\to}0$.
Further, since $e_{i,j}^{(k)}$ are i.i.d sub-Gaussian with parameter at most
$c\sigma^{2}$, for $t\geq 0$ we can apply the sub-Gaussian tail bound

$\mathbb{P}\left(\frac{1}{n_{k}}\left|\sum_{i=1}^{n_{k}}\left(\boldsymbol{e}_{i,j}^{(k)}\right)^{2}-\mathbb{E}\left[\left(e_{i,j}^{(k)}\right)^{2}\right]\right|\geq
c\sigma^{2}t\right)\leq\,c_{1}\exp\left(-c_{2}n_{k}t^{2}\right).$

Note that

$\frac{1}{n_{k}}\mathbb{E}\left[\left(e_{i,j}^{(k)}\right)^{2}\right]=\sigma_{j}^{2}\cdot
f^{(k)}_{\sigma_{j}}(S)$

By setting $t=\sqrt{\frac{\log p}{n_{k}}}$, we therefore obtain the bound

$\displaystyle\sigma_{j}^{2}\left|\hat{f}^{(k)}_{\sigma_{j}}(S)-f^{(k)}_{\sigma_{j}}(S)\right|$
$\displaystyle\leq c\sigma^{2}\sqrt{\frac{\log
p}{n_{k}}}+\left|\frac{1}{n_{k}}\sum_{j=1}^{p}Y_{j}^{2}\right|$

with probability at least $c_{1}\exp\left(-c_{2}\log p\right)$.

We now use Lemma 26. We set $t^{\prime}=\sqrt{\frac{\log p}{n_{j}}}$, which
for $\delta=c^{\prime}\sqrt{\frac{1}{n_{j}}}+c^{\prime\prime}t^{\prime}$ gives

$\displaystyle\sigma_{j}^{2}\left|\hat{f}^{(k)}_{\hat{\sigma}_{j}}(S)-f^{(k)}_{\sigma_{j}}(S)\right|$
$\displaystyle\leq\sigma_{j}^{2}\left(\left|\hat{f}^{(k)}_{\hat{\sigma}_{j}}(S)-\hat{f}^{(k)}_{\sigma_{j}}(S)\right|+\left|\hat{f}^{(k)}_{\sigma_{j}}(S)-f^{(k)}_{\sigma_{j}}(S)\right|\right)$
$\displaystyle\leq\sigma^{2}\cdot\max\\{\delta,\delta^{2}\\}\cdot\frac{C}{n_{k}}\mathbb{E}\left[\left|\left|\boldsymbol{e}^{(k)}_{j}\right|\right|_{2}^{2}\right]$
$\displaystyle+c_{0}\sigma^{2}\sqrt{\frac{\log
p}{n_{k}}}+\left|\frac{1}{n_{k}}\sum_{j=1}^{p}Y_{j}^{2}\right|$
$\displaystyle\leq c_{0}^{\prime}\sigma^{4}\sqrt{\frac{\log
p}{n_{j}}}+c_{1}^{\prime}\sigma^{2}\sqrt{\frac{\log
p}{n_{k}}}+c_{2}^{\prime}\frac{p}{n_{k}}$

if

$\frac{\sigma^{2}}{\sigma_{j}^{2}}\cdot\delta\leq\frac{1}{2}.$

with probability at least $1-c_{1}\exp(-c_{2}\log p)$. ∎

####  7.2.5 Proof of Theorem 4

###### Proof.

Suppose the gap $\xi_{\hat{\Omega}_{0}}^{(k)}$ is nonzero, which we guarantee
with high probability by Theorem 3. Then the following inequality is valid, by
Eq. 27 and Eq. 28:

$\displaystyle\left|\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)-\mathbb{E}\mathcal{L}_{\Omega_{0}}\left(W^{(k)},X^{(k)}\right)\right|$
$\displaystyle\leq\sum_{j=1}^{p}\left|\hat{f}_{\hat{\sigma}_{j}}^{(k)}\left(\\{1\dots
p\\}\backslash\\{j\\}\right)-f_{\sigma_{j}}\left(\\{1\dots
p\\}\backslash\\{j\\}\right)\right|$ (31)
$\displaystyle<\frac{\xi^{(k)}_{\hat{\Omega}_{0}}}{2}$

for all $W\in\mathcal{D}$. Then for all $W\in\mathcal{D}$, $W\neq W_{0}$,

$\displaystyle\mathcal{L}_{\hat{\Omega}_{0}}\left(W_{0}^{(k)},\boldsymbol{X}^{(k)}\right)$
$\displaystyle<\mathbb{E}\mathcal{L}_{\Omega_{0}}\left(W_{0}^{(k)},X^{(k)}\right)+\frac{\xi_{\hat{\Omega}_{0}}^{(k)}}{2}$
$\displaystyle\leq\left(\mathbb{E}\mathcal{L}_{\Omega_{0}}\left(W^{(k)},X^{(k)}\right)-\xi_{\Omega_{0}}^{(k)}\right)+\frac{\xi_{\hat{\Omega}_{0}}^{(k)}}{2}$
$\displaystyle<\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right),$

where the first and third inequalities come from Eq. 31 and the second
inequality comes from the definition of the gap $\xi^{(k)}_{\Omega_{0}}$. ∎

###  7.3 Two node system - dotears

We return to the two-node system described by Eq. 8, and re-examine the system
under the dotears loss $\mathcal{L}_{\Omega_{0}}$, defined in Eq. 20. As
before, we define the ground-truth weighted adjacency matrix
$W_{0}\coloneqq\begin{pmatrix}0&w\\\ 0&0\end{pmatrix}$ and false weighted
adjacency matrix $W_{\delta}\coloneqq\begin{pmatrix}0&0\\\
\delta&0\end{pmatrix}$. For simplicity, we assume we are given
$\Omega_{\alpha}=\mathbb{E}\hat{\Omega}_{0}={\alpha^{2}}\Omega_{0}$, the
expected value of the estimator $\hat{\Omega}_{0}$. We show that
$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(k)},\boldsymbol{X}^{(k)}\right)<\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(k)},\boldsymbol{X}^{(k)}\right)$
for all $k=0,1,2$.

####  7.3.1 Observational system

In the observational system, we retain the generative SEM in Eq. 8:

$\displaystyle X_{1}^{(0)}$ $\displaystyle=\epsilon_{1}^{(0)}$ (32)
$\displaystyle X_{2}^{(0)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{2}^{(0)}$
$\displaystyle=w\epsilon_{1}^{(0)}+\epsilon_{2}^{(0)}.$

To calculate
$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(0)},\boldsymbol{X}^{(0)}\right)$,
we decompose the loss $\mathcal{L}$ component-wise by noting that

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)=\frac{1}{n_{k}}\sum_{i=0}^{p}\frac{\alpha^{2}}{\sigma_{i}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(k)}-\boldsymbol{X}^{(k)}W^{(k)}\right)_{i}\right|\right|_{F}^{2},$

to obtain

$\displaystyle\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{0}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(0)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(0)}$
$\displaystyle\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{0}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(0)}-w\boldsymbol{X}_{1}^{(0)}$
$\displaystyle=w\boldsymbol{\epsilon}_{1}^{(0)}+\boldsymbol{\epsilon}_{2}^{(0)}-w\boldsymbol{\epsilon}_{1}^{(0)}$
$\displaystyle=\boldsymbol{\epsilon}_{2}^{(0)}$

and therefore

$\displaystyle\frac{1}{n_{0}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{0}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\left(\frac{1}{n_{0}}\mathbb{E}\left|\left|\boldsymbol{\epsilon}_{1}^{(0)}\right|\right|_{F}^{2}\right)$
$\displaystyle=\left(\frac{\alpha^{2}}{\sigma_{1}^{2}}\right)\sigma_{1}^{2}$
$\displaystyle=\alpha^{2}$
$\displaystyle\frac{1}{n_{0}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{0}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\left(\frac{1}{n_{0}}\mathbb{E}\left|\left|\boldsymbol{\epsilon}_{2}^{(0)}\right|\right|_{F}^{2}\right)$
$\displaystyle=\left(\frac{\alpha^{2}}{\sigma_{2}^{2}}\right)\sigma_{2}^{2}$
$\displaystyle=\alpha^{2}$

As a result,

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(0)},\boldsymbol{X}^{(0)}\right)=2a^{2}.$

Similarly, we calculate
$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(0)},\boldsymbol{X}^{(0)}\right)$
component-wise:

$\displaystyle\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{\delta}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(0)}-\delta\boldsymbol{X}_{2}^{(0)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(0)}-\delta\left(w\boldsymbol{\epsilon}_{1}^{(0)}+\boldsymbol{\epsilon}_{2}^{(0)}\right)$
$\displaystyle=\left(1-\delta
w\right)\boldsymbol{\epsilon}_{1}^{(0)}-\delta\boldsymbol{\epsilon}_{2}^{(0)}$
$\displaystyle\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{\delta}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(0)}$
$\displaystyle=w\boldsymbol{\epsilon}_{1}^{(0)}+\boldsymbol{\epsilon}_{2}^{(0)}$

Let $\gamma\in\mathbb{R}^{+}$, such that
$\sigma_{1}^{2}=\gamma\sigma_{2}^{2}$. Then we calculate expected loss
component-wise as

$\displaystyle\frac{1}{n_{0}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{\delta}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\left((1-\delta
w)^{2}\sigma_{1}^{2}+\delta^{2}\sigma_{2}^{2}\right)$
$\displaystyle=\alpha^{2}\left((1-\delta
w)^{2}+\frac{\delta^{2}}{\gamma}\right)$
$\displaystyle\frac{1}{n_{0}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(0)}-\boldsymbol{X}^{(0)}W_{\delta}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\left(w^{2}\sigma_{1}^{2}+\sigma_{2}^{2}\right)$
$\displaystyle=\alpha^{2}\left(1+w^{2}\gamma\right)$

As a result,

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(0)},\boldsymbol{X}^{(0)}\right)=\alpha^{2}\left((1-\delta
w)^{2}+\frac{\delta^{2}}{\gamma}+1+w^{2}\gamma\right).$

We can now ask whether
$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(0)},\boldsymbol{X}^{(0)}\right)\geq\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(0)},\boldsymbol{X}^{(0)}\right)$
for all $w,\delta,\gamma$.

$\displaystyle\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(0)},\boldsymbol{X}^{(0)}\right)$
$\displaystyle\overset{?}{\geq}\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(0)},\boldsymbol{X}^{(0)}\right)$
$\displaystyle\alpha^{2}\left((1-\delta
w)^{2}+\frac{\delta^{2}}{\gamma}+1+w^{2}\gamma\right)$
$\displaystyle\overset{?}{\geq}2\alpha^{2}$ $\displaystyle(1-\delta
w)^{2}+\frac{\delta^{2}}{\gamma}+1+w^{2}\gamma$
$\displaystyle\overset{?}{\geq}2$ $\displaystyle 1-2\delta
w+\delta^{2}w^{2}+\frac{\delta^{2}}{\gamma}+w^{2}$
$\displaystyle\overset{?}{\geq}1$

In the end, we obtain the inequality

$w^{2}\left(\delta^{2}+\gamma\right)-w(2\delta)+\frac{\delta^{2}}{\gamma}\overset{?}{\geq}0$
(33)

which is also a quadratic in $w$. We can solve for the roots of $w$:

$\displaystyle w$
$\displaystyle=\frac{2\delta\pm\sqrt{4\delta^{2}-4\left(\frac{\delta^{2}}{\gamma}\right)\left(\delta^{2}+\gamma\right)}}{2\left(\delta^{2}+\gamma\right)}$
$\displaystyle=\frac{2\delta\pm
2\sqrt{\delta^{2}-\frac{\delta^{4}}{\gamma}-\delta^{2}}}{2\left(\delta^{2}+\gamma\right)}$
$\displaystyle=\frac{2\delta\pm
2\sqrt{-\frac{\delta^{4}}{\gamma}}}{2\left(\delta^{2}+\gamma\right)}$

$\gamma\in\mathbb{R}^{+}$ shows that $w$ has no solution in $\mathbb{R}$.
Moreover, the intercept term $\frac{\delta^{2}}{\gamma}>0$ in Eq. 33, proving
the result for all $w,\delta,\gamma$.

####  7.3.2 Intervention on node 1

We proceed similarly for the interventional cases. Upon intervention on node
1, the generative structure does not change, i.e.
$W_{0}^{(1)}=\begin{pmatrix}0&w\\\ 0&0\end{pmatrix}$, with SEM

$\begin{aligned} X_{1}^{(1)}&=\epsilon_{1}^{(1)}\\\
X_{2}^{(1)}&=wX_{1}^{(1)}+\epsilon_{2}^{(1)}\\\
&=w\epsilon_{1}^{(1)}+\epsilon_{2}^{(1)}\end{aligned}.$

Note that
$\text{Var}\left(\epsilon_{i}^{(1)}\right)=\frac{\sigma_{1}^{2}}{\alpha^{2}}=\frac{1}{\alpha^{2}}\text{Var}\left(\epsilon_{i}^{(0)}\right)$,
in accordance with Assumption 3. Component-wise under $W_{0}^{(1)}$ we have

$\displaystyle\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{0}^{(1)}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(1)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(1)}$
$\displaystyle\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{0}^{(1)}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(1)}-w\boldsymbol{X}_{1}^{(1)}$
$\displaystyle=\boldsymbol{X}_{2}^{(1)}-w\boldsymbol{X}_{1}^{(1)}$
$\displaystyle=w\boldsymbol{\epsilon}_{1}^{(1)}+\boldsymbol{\epsilon}_{2}^{(1)}-w\boldsymbol{\epsilon}_{1}^{(1)}$
$\displaystyle=\boldsymbol{\epsilon}_{2}^{(1)}$

which gives us the expected component-wise losses

$\displaystyle\frac{1}{n_{1}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{0}^{(1)}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\frac{\sigma_{1}^{2}}{\alpha^{2}}$
$\displaystyle=1$
$\displaystyle\frac{1}{n_{1}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{0}^{(1)}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\sigma_{2}^{2}$
$\displaystyle=\alpha^{2}$

and thus

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(1)},\boldsymbol{X}^{(1)}\right)=1+\alpha^{2}.$

Note that under intervention on $X_{1}$,
$W_{\delta}^{(1)}=\begin{pmatrix}0&0\\\ 0&0\end{pmatrix}$. Component-wise, we
then have

$\displaystyle\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{\delta}^{(1)}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(1)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(1)}$
$\displaystyle\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{\delta}^{(1)}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(1)}$
$\displaystyle=w\boldsymbol{\epsilon}_{1}^{(1)}+\boldsymbol{\epsilon}_{2}^{(1)}$

and the expected component-wise loss

$\displaystyle\frac{1}{n_{1}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{\delta}^{(1)}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\frac{\sigma_{1}^{2}}{\alpha^{2}}$
$\displaystyle=1$
$\displaystyle\frac{1}{n_{1}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(1)}-\boldsymbol{X}^{(1)}W_{\delta}^{(1)}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\left(w^{2}\frac{\sigma_{1}^{2}}{\alpha^{2}}+\sigma_{2}^{2}\right)$
$\displaystyle=\alpha^{2}\left(w^{2}\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}+1\right).$

Thus,

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(1)},\boldsymbol{X}^{(1)}\right)=1+\alpha^{2}+\alpha^{2}w^{2}\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}>\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(1)},\boldsymbol{X}^{(1)}\right),$

which proves the result.

####  7.3.3 Intervention on node 2

Under intervention on $X_{2}$, $W_{0}^{(2)}=\begin{pmatrix}0&0\\\
0&0\end{pmatrix}$, with corresponding SEM

$\begin{aligned} X_{1}^{(2)}&=\epsilon_{1}^{(2)}\\\
X_{2}^{(2)}&=\epsilon_{2}^{(2)}\\\ \end{aligned}.$

Component-wise, we obtain the terms

$\displaystyle\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{0}^{(2)}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(2)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(2)}$
$\displaystyle\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{0}^{(2)}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(2)}$
$\displaystyle=\boldsymbol{\epsilon}_{2}^{(2)}$

Note that
$\text{Var}\left(\epsilon_{2}^{(2)}\right)=\frac{\sigma_{2}}{\alpha^{2}}$ in
accordance with Assumption 3. Then the expected loss component-wise is

$\displaystyle\frac{1}{n_{2}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{0}^{(2)}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\sigma_{1}^{2}$
$\displaystyle=\alpha^{2}$
$\displaystyle\frac{1}{n_{2}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{0}^{(2)}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\frac{\sigma_{2}^{2}}{\alpha^{2}}$
$\displaystyle=1,$

giving

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{0}^{(2)},\boldsymbol{X}^{(2)}\right)=1+\alpha^{2}.$

Under $W_{\delta}^{(2)}=\begin{pmatrix}0&0\\\ \delta&0\end{pmatrix}$, we
obtain the terms component-wise

$\displaystyle\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{\delta}^{(2)}\right)_{1}$
$\displaystyle=\boldsymbol{X}_{1}^{(2)}-\delta\boldsymbol{X}_{2}^{(2)}$
$\displaystyle=\boldsymbol{\epsilon}_{1}^{(2)}-\delta\boldsymbol{\epsilon}_{2}^{(2)}$
$\displaystyle\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{\delta}^{(2)}\right)_{2}$
$\displaystyle=\boldsymbol{X}_{2}^{(2)}$
$\displaystyle=\boldsymbol{\epsilon}_{2}^{(2)}$

which in expectation gives the component-wise losses

$\displaystyle\frac{1}{n_{2}}\frac{\alpha^{2}}{\sigma_{1}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{\delta}^{(2)}\right)_{1}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{1}^{2}}\left(\sigma_{1}^{2}+\delta^{2}\frac{\sigma_{2}^{2}}{\alpha^{2}}\right)$
$\displaystyle=\alpha^{2}+\delta^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}}$
$\displaystyle\frac{1}{n_{2}}\frac{\alpha^{2}}{\sigma_{2}^{2}}\mathbb{E}\left|\left|\left(\boldsymbol{X}^{(2)}-\boldsymbol{X}^{(2)}W_{\delta}^{(2)}\right)_{2}\right|\right|_{F}^{2}$
$\displaystyle=\frac{\alpha^{2}}{\sigma_{2}^{2}}\frac{\sigma_{2}^{2}}{\alpha^{2}},$

which gives

$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(2)},\boldsymbol{X}^{(2)}\right)=1+\alpha^{2}+\delta^{2}\frac{\sigma_{2}^{2}}{\sigma_{1}^{2}}$

which is trivially greater than
$\mathbb{E}\mathcal{L}_{\Omega_{\alpha}}\left(W_{\delta}^{(2)},\boldsymbol{X}^{(2)}\right)$
for all $\alpha,\delta,\sigma_{1}^{2},\sigma_{2}^{2}$.

###  7.4 Full two node simulation details

In Section 3.1, we simulate from the system of SEMs

$\displaystyle X^{(0)}$ $\displaystyle=X^{(0)}W+\epsilon^{(0)}$ $\displaystyle
X^{(1)}$ $\displaystyle=X^{(1)}W^{(1)}+\epsilon^{(1)}$ $\displaystyle X^{(2)}$
$\displaystyle=X^{(2)}W^{(2)}+\epsilon^{(2)},$

such that in the observational system

$\displaystyle X_{1}^{(0)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(0)}\sim\mathcal{N}\left(0,\gamma\right)$ (34)
$\displaystyle X_{2}^{(0)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{2}^{(0)}$
$\displaystyle\epsilon_{2}^{(0)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 1 we have

$\displaystyle X_{1}^{(1)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(1)}\sim\mathcal{N}\left(0,\frac{\gamma}{\alpha^{2}}\right)$
(35) $\displaystyle X_{2}^{(1)}$
$\displaystyle=wX_{1}^{(1)}+\epsilon_{2}^{(1)}$
$\displaystyle\epsilon_{2}^{(1)}\sim\mathcal{N}\left(0,1\right),$

and under intervention on node 2

$\displaystyle X_{1}^{(2)}$ $\displaystyle=\epsilon_{1}^{(2)}$
$\displaystyle\epsilon_{1}^{(2)}\sim\mathcal{N}\left(0,\gamma\right)$ (36)
$\displaystyle X_{2}^{(2)}$ $\displaystyle=\epsilon_{2}^{(2)}$
$\displaystyle\epsilon_{2}^{(2)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right).$

For observational data, we simulate only from the SEM given in Eq. 34, with a
sample size $n=3000$. For interventional data, we draw a sample size of
$n_{k}=1000$ from each intervention $k$. This gives a total sample size of
$n=(p+1)\times 3000$ for interventional data, which matches the sample size of
observational data. We set $\alpha=4$. In Section 3.1, we presented empirical
results on two node simulations for
$(w,\gamma)\in\\{0.3,0.5,0.7,0.8,1.0\\}\times\\{1,2,100\\}$. Figure 5 shows
the $\ell_{1}$ distance between the ground truth DAG $W_{0}$ and the inferred
DAG for each method, for all
$(w,\gamma)\in\\{0.1,0.2,\dots,1.5\\}\times\\{1,2,4,10,100\\}$. For each
parameter combination of $(w,\gamma)$, we draw 25 instances of simulated data
for both the interventional and strictly observational case, at the described
sample size.

For dotears, sortnregress, and NO TEARS, we set the regularization parameter
$\lambda$ to 0, to isolate the performance of the loss function. In GOLEM-EV
and GOLEM-NV we set the $\ell_{1}$ regularization parameter $\lambda_{1}$ to
0, but let $\lambda_{2}=5$ to enforce DAG-ness, as recommended by the authors.

See Figure 5 for the full set of simulated $w,\gamma$.

Figure 5: Full simulations for the two-node DAG for all
$(w,\gamma)\in\\{0.1,0.2,\dots,1.5\\}\times\\{1,2,4,10,100\\}$, where
$X_{1}\overset{w}{\rightarrow}X_{2}$ and
$\sigma_{1}^{2}=\gamma\sigma_{2}^{2}$. Methods are compared using SHD (lower
is better).

###  7.5 Full three node simulations

Simulations in the two-node DAG showed that NO TEARS and GOLEM-NV are
deterministic functions of the true structure $W$ and the true exogenous
variance structure $\Omega_{0}$ (Sections 3.1, 7.4). Simulations in more
complex three-node topologies verify determinism of NO TEARS and GOLEM-NV in
$W$ and $\Omega_{0}$, but show they are distinct deterministic functions.

We simulate observational and interventional data under the three node chain
$X_{1}\overset{w}{\rightarrow}X_{2}\overset{w}{\rightarrow}X_{3}$, the three
node collider
$X_{1}\overset{w}{\rightarrow}X_{3}\overset{w}{\leftarrow}X_{2}$, and the
three node fork
$X_{2}\overset{w}{\leftarrow}X_{1}\overset{w}{\rightarrow}X_{3}$. In each
topology, we let $\Omega_{0}\coloneqq\begin{pmatrix}\gamma&0&0\\\ 0&1&0\\\
0&0&1\end{pmatrix}$, so that $\gamma$ is the exogenous variance of a source
node. For simplicity, we constrain the edge weights to be equal, and simulate
data under Gaussian exogenous variance for
$(w,\gamma)\in\\{0.1,\dots,1.5\\}\times\\{1,2,4,10,100\\}$. For interventional
data with interventions $i\in\\{0,1,2,3\\}$, we draw $n_{i}=1000$ observations
in all simulations, giving us a total sample size of $n=4000$ that is matched
in observational data. As in section 3.1, we set $\alpha=4$, and remove
$\ell_{1}$ regularization where appropriate. We benchmark dotears, NO TEARS,
sortnregress, GOLEM-EV, GOLEM-NV, DirectLingam, GIES, IGSP, UT-IGSP, and
DCDI-G [10, 12, 13, 15, 19, 20, 21, 22, 23].

We evaluate each method by the SHD between the ground truth DAG $W$ and the
method inferred DAG. In the two node case, the least squares loss performed
identically to the likelihood loss of GOLEM-NV; here, their performances
diverge, but remain essentially deterministic functions in $w$ and $\gamma$.

####  7.5.1 Chain

In the chain, we have the true structure $W_{0}=\begin{pmatrix}0&w&0\\\
0&0&w\\\ 0&0&0\end{pmatrix}$. We simulate under the system of SEMs

$\displaystyle X^{(0)}$ $\displaystyle=X^{(0)}W_{0}+\epsilon^{(0)}$
$\displaystyle X^{(1)}$ $\displaystyle=X^{(1)}W_{0}^{(1)}+\epsilon^{(1)}$
$\displaystyle X^{(2)}$ $\displaystyle=X^{(2)}W_{0}^{(2)}+\epsilon^{(2)},$
$\displaystyle X^{(3)}$ $\displaystyle=X^{(3)}W_{0}^{(3)}+\epsilon^{(3)},$

such that in the observational system

$\displaystyle X_{1}^{(0)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(0)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(0)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{2}^{(0)}$
$\displaystyle\epsilon_{2}^{(0)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(0)}$ $\displaystyle=wX_{2}^{(0)}+\epsilon_{3}^{(0)}$
$\displaystyle\epsilon_{3}^{(0)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 1 we have

$\displaystyle X_{1}^{(1)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(1)}\sim\mathcal{N}\left(0,\frac{\gamma}{\alpha^{2}}\right)$
$\displaystyle X_{2}^{(1)}$ $\displaystyle=wX_{1}^{(1)}+\epsilon_{2}^{(1)}$
$\displaystyle\epsilon_{2}^{(1)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$ $\displaystyle=wX_{2}^{(1)}+\epsilon_{3}^{(1)}$
$\displaystyle\epsilon_{3}^{(1)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 2

$\displaystyle X_{1}^{(2)}$ $\displaystyle=\epsilon_{1}^{(2)}$
$\displaystyle\epsilon_{1}^{(2)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(2)}$ $\displaystyle=\epsilon_{2}^{(2)}$
$\displaystyle\epsilon_{2}^{(2)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right),$
$\displaystyle X_{3}^{(2)}$ $\displaystyle=wX_{2}^{(2)}+\epsilon_{3}^{(2)}$
$\displaystyle\epsilon_{3}^{(2)}\sim\mathcal{N}\left(0,1\right),$

and under intervention on node 3

$\displaystyle X_{1}^{(3)}$ $\displaystyle=\epsilon_{1}^{(3)}$
$\displaystyle\epsilon_{1}^{(3)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(3)}$ $\displaystyle=wX_{1}^{(3)}+\epsilon_{2}^{(3)}$
$\displaystyle\epsilon_{2}^{(3)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$ $\displaystyle=\epsilon_{3}^{(3)}$
$\displaystyle\epsilon_{3}^{(3)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right).$

Figure 6 shows the full results.

(a) Results for $w$ in $0.1,\dots,0.8$, inclusive.

(b) Results for $w$ in $0.8,\dots,1.5$, inclusive.

Figure 6: Full simulations for the chain topology
$X_{1}\overset{w}{\rightarrow}X_{2}\overset{w}{\rightarrow}X_{3}$, where
$\Omega_{0}=\begin{pmatrix}\gamma&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{pmatrix}$. We
iterate across all
$(w,\gamma)\in\\{0.1,0.2,\dots,1.5\\}\times\\{1,2,4,10,100\\}$. Methods are
compared using SHD (lower is better). (a) shows results for $w$ in
$0.1,\dots,0.8$, inclusive, while (b) shows results for $w$ in
$0.8,\dots,1.5$, inclusive. $w=0.8$ is repeated across subpanels.

####  7.5.2 Collider

In the collider, we have the true structure $W_{0}=\begin{pmatrix}0&0&w\\\
0&0&w\\\ 0&0&0\end{pmatrix}$. We simulate under the system of SEMs

$\displaystyle X^{(0)}$ $\displaystyle=X^{(0)}W_{0}+\epsilon^{(0)}$
$\displaystyle X^{(1)}$ $\displaystyle=X^{(1)}W_{0}^{(1)}+\epsilon^{(1)}$
$\displaystyle X^{(2)}$ $\displaystyle=X^{(2)}W_{0}^{(2)}+\epsilon^{(2)},$
$\displaystyle X^{(3)}$ $\displaystyle=X^{(3)}W_{0}^{(3)}+\epsilon^{(3)},$

such that in the observational system

$\displaystyle X_{1}^{(0)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(0)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(0)}$ $\displaystyle=\epsilon_{2}^{(0)}$
$\displaystyle\epsilon_{2}^{(0)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(0)}$
$\displaystyle=wX_{1}^{(0)}+wX_{2}^{(0)}+\epsilon_{3}^{(0)}$
$\displaystyle\epsilon_{3}^{(0)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 1 we have

$\displaystyle X_{1}^{(1)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(1)}\sim\mathcal{N}\left(0,\frac{\gamma}{\alpha^{2}}\right)$
$\displaystyle X_{2}^{(1)}$ $\displaystyle=\epsilon_{2}^{(1)}$
$\displaystyle\epsilon_{2}^{(1)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$
$\displaystyle=wX_{1}^{(1)}+wX_{2}^{(1)}+\epsilon_{3}^{(1)}$
$\displaystyle\epsilon_{3}^{(1)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 2

$\displaystyle X_{1}^{(2)}$ $\displaystyle=\epsilon_{1}^{(2)}$
$\displaystyle\epsilon_{1}^{(2)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(2)}$ $\displaystyle=\epsilon_{2}^{(2)}$
$\displaystyle\epsilon_{2}^{(2)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right),$
$\displaystyle X_{3}^{(2)}$
$\displaystyle=wX_{1}^{(2)}+wX_{2}^{(2)}+\epsilon_{3}^{(2)}$
$\displaystyle\epsilon_{3}^{(2)}\sim\mathcal{N}\left(0,1\right),$

and under intervention on node 3

$\displaystyle X_{1}^{(3)}$ $\displaystyle=\epsilon_{1}^{(3)}$
$\displaystyle\epsilon_{1}^{(3)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(3)}$ $\displaystyle=\epsilon_{2}^{(3)}$
$\displaystyle\epsilon_{2}^{(3)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$ $\displaystyle=\epsilon_{3}^{(3)}$
$\displaystyle\epsilon_{3}^{(3)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right).$

(a) Results for $w$ in $0.1,\dots,0.8$, inclusive.

(b) Results for $w$ in $0.8,\dots,1.5$, inclusive.

Figure 7: Full simulations for the collider topology
$X_{1}\overset{w}{\rightarrow}X_{3}\overset{w}{\leftarrow}X_{2}$, where
$\Omega_{0}=\begin{pmatrix}\gamma&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{pmatrix}$. We
iterate across all
$(w,\gamma)\in\\{0.1,0.2,\dots,1.5\\}\times\\{1,2,4,10,100\\}$. Methods are
compared using SHD (lower is better). (a) shows results for $w$ in
$0.1,\dots,0.8$, inclusive, while (b) shows results for $w$ in
$0.8,\dots,1.5$, inclusive. $w=0.8$ is repeated across subpanels.

####  7.5.3 Fork

In the collider, we have the true structure $W_{0}=\begin{pmatrix}0&w&w\\\
0&0&0\\\ 0&0&0\end{pmatrix}$. We simulate under the system of SEMs

$\displaystyle X^{(0)}$ $\displaystyle=X^{(0)}W_{0}+\epsilon^{(0)}$
$\displaystyle X^{(1)}$ $\displaystyle=X^{(1)}W_{0}^{(1)}+\epsilon^{(1)}$
$\displaystyle X^{(2)}$ $\displaystyle=X^{(2)}W_{0}^{(2)}+\epsilon^{(2)},$
$\displaystyle X^{(3)}$ $\displaystyle=X^{(3)}W_{0}^{(3)}+\epsilon^{(3)},$

such that in the observational system

$\displaystyle X_{1}^{(0)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(0)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(0)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{2}^{(0)}$
$\displaystyle\epsilon_{2}^{(0)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(0)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{3}^{(0)}$
$\displaystyle\epsilon_{3}^{(0)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 1 we have

$\displaystyle X_{1}^{(1)}$ $\displaystyle=\epsilon_{1}^{(1)}$
$\displaystyle\epsilon_{1}^{(1)}\sim\mathcal{N}\left(0,\frac{\gamma}{\alpha^{2}}\right)$
$\displaystyle X_{2}^{(1)}$ $\displaystyle=wX_{1}^{(0)}+\epsilon_{2}^{(1)}$
$\displaystyle\epsilon_{2}^{(1)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$ $\displaystyle=wX_{1}^{(1)}+\epsilon_{3}^{(1)}$
$\displaystyle\epsilon_{3}^{(1)}\sim\mathcal{N}\left(0,1\right),$

under intervention on node 2

$\displaystyle X_{1}^{(2)}$ $\displaystyle=\epsilon_{1}^{(2)}$
$\displaystyle\epsilon_{1}^{(2)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(2)}$ $\displaystyle=\epsilon_{2}^{(2)}$
$\displaystyle\epsilon_{2}^{(2)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right),$
$\displaystyle X_{3}^{(2)}$ $\displaystyle=wX_{1}^{(2)}+\epsilon_{3}^{(2)}$
$\displaystyle\epsilon_{3}^{(2)}\sim\mathcal{N}\left(0,1\right),$

and under intervention on node 3

$\displaystyle X_{1}^{(3)}$ $\displaystyle=\epsilon_{1}^{(3)}$
$\displaystyle\epsilon_{1}^{(3)}\sim\mathcal{N}\left(0,\gamma\right)$
$\displaystyle X_{2}^{(3)}$ $\displaystyle=wX_{1}^{(3)}+\epsilon_{2}^{(3)}$
$\displaystyle\epsilon_{2}^{(3)}\sim\mathcal{N}\left(0,1\right),$
$\displaystyle X_{3}^{(1)}$ $\displaystyle=\epsilon_{3}^{(3)}$
$\displaystyle\epsilon_{3}^{(3)}\sim\mathcal{N}\left(0,\frac{1}{\alpha^{2}}\right).$

(a) Results for $w$ in $0.1,\dots,0.8$, inclusive.

(b) Results for $w$ in $0.8,\dots,1.5$, inclusive.

Figure 8: Full simulations for the fork topology
$X_{2}\overset{w}{\leftarrow}X_{1}\overset{w}{\rightarrow}X_{3}$, where
$\Omega_{0}=\begin{pmatrix}\gamma&0&0\\\ 0&1&0\\\ 0&0&1\\\ \end{pmatrix}$. We
iterate across all
$(w,\gamma)\in\\{0.1,0.2,\dots,1.5\\}\times\\{1,2,4,10,100\\}$. Methods are
compared using SHD (lower is better). (a) shows results for $w$ in
$0.1,\dots,0.8$, inclusive, while (b) shows results for $w$ in
$0.8,\dots,1.5$, inclusive. $w=0.8$ is repeated across subpanels.

###  7.6 Large random graph simulations

####  7.6.1 Data generation

For evaluation, data is generated for DAGs with $p=40$ nodes. Structures are
drawn from both Erdős-Rényi (ER) and Scale-Free (SF) DAGs [24, 25]. We
parameterize ER graphs as ER-$r$ and SF graphs as SF-$z$, where $r\in[0,1]$
represents the probability of assignment of each individual edge, whereas
$z\in\mathbb{Z}$ is the integer number of edges assigned per node.

We simulate under two parameterizations $\\{\text{Low Density, High
Density}\\}$ of the edge densities $(r,z)$. In the “Low Density”
parameterization, we give $(r,z)=(0.1,2)$. To evaluate performance on higher
density topologies, we also give “High Density” parameterizations, where
$(r,z)=(0.2,4)$.

Given an edge density scenario, we simulate under two parameterizations of the
edge weights. In the “Strong Effects” parameterization,
$w\sim\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$. We also give the
“Weak Effects” parameterization, which edge weights are drawn from
$w\sim\text{Unif}\left([-1.0,0.3]\cup[0.3,1.0]\right)$ such that $|w|\leq 1$
is guaranteed. Here, there is no guarantee any two nodes will be varsortable
(although they may be in practice, as a function of $\Omega_{0}$) [13, 18].
Table 1 summarizes all four possible simulation parameterizations.

For each node $i$, we draw $\sigma_{i}\sim\text{Unif}\left([0.5,2.0]\right)$,
and draw $n_{0}$ observations from
$\boldsymbol{\epsilon}_{i}^{(0)}\sim\mathcal{N}\left(0,\sigma_{i}^{2}\right)$.
For each DAG we generate an instance of observational data, where
$n_{0}=(p+1)*n_{k}=4100$, and an instance of interventional data, where
$n_{k}=100$ for all $k=0\dots p$ to match sample size. We set the distribution
of $\epsilon_{i}^{(k)}$ according to Assumptions 3 and 2 for $\alpha=4$.

Simulation Type | Topology | $w$
---|---|---
Low Density | , Strong Effects | ER-0.1 | OR SF-2 | $\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$
Low Density | , Weak Effects | ER-0.1 | OR SF-2 | $\text{Unif}\left([-1.0,-0.3]\cup[0.3,1.0]\right)$
High Density | , Strong Effects | ER-0.2 | OR SF-4 | $\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$
High Density | , Weak Effects | ER-0.2 | OR SF-4 | $\text{Unif}\left([-1.0,-0.3]\cup[0.3,1.0]\right)$
Table 1: For each simulation type for $p=40$ DAG simulations, we provide the
parameterizations of the Erdős-Rényi and Scale-Free topologies and the
distribution on the drawn edge weights.

For dotears, NO TEARS, sortnregress, GOLEM-EV, GOLEM-NV, and DCDI-G 5-fold
cross-validation was performed to select the regularization parameters. For
each drawn DAG, a separate data instance of interventional and observational
data was re-drawn from the same distribution specifically for cross-
validation. After choosing a $\lambda$ (or for GOLEM-EV and GOLEM-NV, the set
$(\lambda_{1},\lambda_{2})$) from the data for cross-validation, the methods
were evaluated on the original simulated data. For dotears, NO TEARS,
sortnregress, and DCDI-G, 5-fold cross-validation was performed across the
grid $\lambda\in\\{.001,.01,.1,1,10,100\\}$. For GOLEM-EV and GOLEM-NV, 5-fold
cross-validation was performed across the grid
$\lambda_{1}\times\lambda_{2}\in\\{.001,.01,.1,1,10,100\\}\times\\{.05,.5,5,50\\}$.
We threshold at 0.2 for “Weak Effects” simulations, and at 0.3 for “Strong
Effects” simulations. Methods are generally robust to thresholding choice.
Results for precision and recall on thresholded edges are shown in Figures 11
and 12, respectively.

####  7.6.2 Benchmarking

In Figure 9, we benchmark wallclock time and memory usage for all methods in
$p=40$ simulations [30] on the UCLA hoffman2 cluster. All continuous
optimization methods (dotears, NO TEARS, GOLEM-EV, and GOLEM-NV) have
significantly higher average runtimes than other methods, which is partially
explained by cross-validation procedures. dotears has relatively light memory
usage, outperformed only by NO TEARS and sortnregress.

DCDI-G has enormous memory requirements, especially relative to other methods.
We note that for DCDI-G, the reported benchmarks do not include memory usage
or runtime from cross-validation folds. Instead, we report only the “main” run
of DCDI-G, and thus DCDI-G is not denoted with the * indicating cross-
validation.

(a) Wallclock time across all $p=40$ simulations in seconds. Y-axis is in log-
scale.

(b) Memory usage across all $p=40$ simulations, defined as the maximum unique
set size used by the process, or max_uss (see psutil documentation).

Figure 9: Comparison of mean runtime (a) and mean memory usage (b) between
methods across all $p=40$ simulations. Error bars are standard errors.

####  7.6.3 Thresholding on large simulation results

Edge thresholding for weighted adjacency matrices is necessary for accurate
evaluation using SHD, but the choice of threshold can feel arbitrary. We find
that methods are generally robust to thresholding choice, following similar
results from [13].

Figure 10 examines the effect on thresholding of small weights in $W$ in large
random DAG simulations, for methods that infer weighted adjacency matrices
(dotears, NO TEARS, sortnregress, DirectLiNGAM, GOLEM-NV, and GOLEM-EV) [12,
13, 15, 23]. For simplicity, we summarize simulations in terms of the
generative edge distribution. Simulations with “Weak Effects”,
$w\sim\text{Unif}\left([-1.0,-0.3]\cup[0.3,1.0]\right)$, are shown in Figure
10(a); simulations with “Strong Effects”,
$w\sim\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$ are shown in Figure
10(b). We compare the SHD between the ground truth adjacency matrix and the
inferred adjacency matrix for each method at all thresholds between 0 and the
absolute lower bound of the true edge weight distribution (0.3 for “Weak
Effects”, 0.5 for “Strong Effects”). Any edge whose magnitude is below the
chosen threshold is set to 0 for SHD evaluation.

(a) Structural Hamming Distance (lower is better) as a function of $\hat{w}$
threshold for “Weak Effects” scenarios, where
$w\sim\text{Unif}\left([-1.0,-0.3]\cup[0.3,1.0]\right)$.

(b) Structural Hamming Distance (lower is better) as a function of $\hat{w}$
threshold for “Strong Effects” scenarios, where
$w\sim\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$.

Figure 10: Structural Hamming Distance (lower is better) as a function of
thresholding choice for $p=40$ simulations detailed in 7.6.1.
Parameterizations with “Weak Effects”,
$w\sim\text{Unif}\left([-1.0,-0.3]\cup[0.3,1.0]\right)$, are consolidated in
(a), while parameterizations with “Strong Effects”,
$w\sim\text{Unif}\left([-2.0,-0.5]\cup[0.5,2.0]\right)$, are consolidated in
(b).

Figure 10 shows that thresholding is necessary for evaluation of SHD between
weighted adjacency matrices, but also that methods are generally robust to
thresholding choice. Without thresholding (equivalently, a threshold of 0),
SHD results are inflated, but recover even at low thresholds.

Figure 11: Methodological comparison of structure learning methods on large
random graphs ($p=40$) using precision (higher is better). Each row represents
simulations drawn with a topology in {Erdős-Rényi, Scale Free}, respectively,
and each column represents a different parameterization of edge density and
edge weights ordered in terms of increasing difficulty. For details, see Table
1. 10 simulations were drawn for each set of parameters with sample size
$(p+1)*100=4100$, and * indicates methods performed with 5-fold cross-
validation on a separate draw from the same distribution. Y-axis scales are
not shared between figures. Methods are sorted, left to right, by strongest
average performance across all simulations. GIES, IGSP, and UT-IGSP infer
binary edges and are therefore omitted for fairness. Figure 12: Methodological
comparison of structure learning methods on large random graphs ($p=40$) using
recall (higher is better). Each row represents simulations drawn with a
topology in {Erdős-Rényi, Scale Free}, respectively, and each column
represents a different parameterization of edge density and edge weights
ordered in terms of increasing difficulty. For details, see Table 1. 10
simulations were drawn for each set of parameters with sample size
$(p+1)*100=4100$, and * indicates methods performed with 5-fold cross-
validation on a separate draw from the same distribution. Y-axis scales are
not shared between figures. Methods are sorted, left to right, by strongest
average performance across all simulations. GIES, IGSP, and UT-IGSP infer
binary edges and are therefore omitted for fairness.

####  7.6.4 Edge weight estimation for large random simulations

Accurate estimation of edge weights is important for understanding structure.
Figure 3 gives results on structural recovery through SHD, but does not inform
edge weight recovery. To measure edge weight recovery we use $\ell_{1}$
distance, defined as the vector $\ell_{1}$ norm between the flattened true
weighted adjacency matrix and the flattened inferred weighted adjacency
matrix. $\ell_{1}$ distance gives information on both structure recovery and
edge weight estimation simultaneously. In Figure 13, we benchmark the recovery
of edge weights for methods that return a weighted adjacency matrix (dotears,
sortnregress, NO TEARS, GOLEM-EV, GOLEM-NV, DirectLiNGAM [12, 13, 15, 23]) for
simulations given in 7.6.1. For fairness, we exclude methods that only return
a binary adjacency matrix. Methods are thresholded in the same manner as in
Figure 3.

Figure 13: Methodological comparison of structure learning methods on large
random graphs ($p=40$) using $\ell_{1}$ distance (lower is better). Each row
represents simulations drawn with a topology in {Erdős-Rényi, Scale Free},
respectively, and each column represents a different parameterization of edge
density and edge weights ordered in terms of increasing difficulty. For
details, see Table 1. 10 simulations were drawn for each set of parameters
with sample size $(p+1)*100=4100$, and * indicates methods performed with
5-fold cross-validation on a separate draw from the same distribution. Y-axis
scales are not shared between figures. Methods are sorted, left to right, by
strongest average performance across all simulations. GIES, IGSP, and UT-IGSP
infer binary edges and are therefore omitted for fairness.

dotears outperforms all other methods in terms of effect size recovery. In
addition, the relative ordering of the methods stays consistent with Figure 3.

####  7.6.5 Statistical significance of SHD distribution difference

For each method, and across all simulations in Section 7.6.1, we compare the
SHD distributions, marginalized across all simulation parameterizations,
against that of dotears. For each method, we perform a one-sided Mann-Whitney
U test to test difference between the method’s SHD distribution and the SHD
distribution of dotears. The resulting p-values are reported in Table 2, and
are significant for all methods.

Method | p-value
---|---
GIES | 2.08e-7
sortnregress | 1.1e-21
NO TEARS | 9.27e-23
DCDI-G | 1.04e-24
GOLEM-NV | 8.45e-28
GOLEM-EV | 4.47e-28
DirectLiNGAM | 4.47e-28
UT-IGSP | 4.467e-28
IGSP | 4.47e-28
Table 2: p-value from asymptotic one-sided Mann-Whitney U test comparing
distribution of indicated method SHD against distribution of dotears SHD.

####  7.6.6 Incorporation of interventional loss

We note that inferring $W$ from strictly observational data, given
$\hat{\Omega}_{0}$ from interventional data, has convenient theoretical
properties but ignores the majority of our data. Figure 14 compares the
performance, in Structural Hamming Distance, of dotears using both
observational and interventional data
$\left(\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)\text{
for }k=0\dots p\right)$ with dotears using only observational data
$\left(\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)\text{
for }k=0\right)$. We use the simulations drawn in Section 7.6.1, aggregated
across all scenarios in Table 1. Note that $\hat{\Omega}_{0}$ is estimated
identically for both.

Figure 14: Comparison of dotears performance when using only
$\mathcal{L}_{\hat{\Omega}_{0}}\left(W^{(k)},\boldsymbol{X}^{(k)}\right)$ for
observational data ($k=0$), right, and when including interventional data
($k=0\dots p$), left. Methods are run on simulation data generated in Section
7.6.1, averaged across parameterizations, and compared using Structural
Hamming Distance.

As expected, restricting dotears to inference only using the observational
data (Figure 14, right) drastically decreases performance when compared to
inference using both observational and interventional data (Figure 14, left).
This motivates including $k=1\dots p$ into the loss and proving consistency
for all $k$.

###  7.7 Sensitivity Analysis

####  7.7.1 Sensitivity under linear SEM

Figure 15: Method performance under nonlinear simulations in $p=10$ random
graphs, in Structural Hamming Distance (lower is better). Data are generated
by a neural network with Additive Noise Model (ANM) or a neural network with
non-additive noise (NN) model. See [22] and Supplementary Material 7.9 for
details. Scenarios are sorted by increasing difficulty. Methods are sorted by
average performance across scenarios.

We performed sensitivity analyses under a linear SEM on Erdős-Rényi graphs
($p=20$) with two edge density parameterizations crossed with two edge weight
parameterizations (we note that the specific density and edge weight
parameterizations differ from Section 4.1). We report results averaged over
$10$ replicates with a sample size of $n_{k}=100$. We report SHD in the main
text. For precision/recall for edge recovery, parameterization-specific
results, details on data generation, and modelling specifics, see
Supplementary Material 7.8.

We assessed the sensitivity of our results to violation of Assumption 1 by
allowing for $\pi$ percent of the observational parental variance into the
marginal variance of the target . Here, $\pi=0$ is a hard intervention. The
accuracy of dotears decreases as $\pi$ increases, or equivalently as the
interventional signal decreases (Figure 16(a)). This effect is dependent on
the magnitude of the edge weights and the density of the true DAG: in “High
Density” and “Strong Effects” scenarios, the decrease in accuracy is
noticeable due to higher incoming parental variance (Supplementary Figure 17).
Despite this decrease in accuracy, dotears remains the second best performing
method across parameterizations, even when excluding the $\pi=0$ scenario.

To assess sensitivity to violation of Assumption 3, we randomly perturb each
$\alpha$ per node to allow non-uniform interventional effects on error
variance across targets. dotears is robust to violations of Assumption 3 in
these simulations (Figure 16(b)).

Finally, our model assumes that the distribution of the target $X_{k}^{(k)}$
depends on that of $\epsilon_{k}^{(0)}$. In Figure 16(c), we explore
interventions with a fixed distribution,
$X_{k}^{(k)}\overset{iid}{\sim}\mathcal{N}(2,1)$ following the model used in
[22] and similar to that considered by [10]. In our model, this violates
Assumption 3.

In Figure 16(c), dotears is outperformed on average only by GIES in SHD, and
outperforms DCDI under their simulation model. Supplementary Figure 20 shows
that simulations with weak effects, where $0.3\leq|w|\leq 0.5$ drive this
performance dip; in simulations with strong effects, where $0.8\leq|w|\leq
1.0$, dotears again vastly outperforms other methods, including GIES.
Sufficient signal-to-noise ratio can thus overcome severe misspecification of
$\Omega_{0}$ under hard interventions.

####  7.7.2 Sensitivity under nonlinear SEM

To assess sensitivity to a nonlinear SEM, we generate nonlinear data under
models that were used to evaluate DCDI [22]. We simulated $10$ Erdős-Rényi
DAGs with $p=10$ and $n=10000$ total samples. For the non-linear models, we
considered, in turn, a neural network additive noise model (ANM) [31] and a
neural network with non-additive noise model (NN) [32]. We test both hard
interventions and “imperfect” interventions. Imperfect interventions add a
random vector drawn from $\mathcal{N}(0,1)$ to the last layer of the neural
network. See [22] and Supplementary Material 7.9 for details. Figure 15 shows
results in SHD.

With hard interventions, dotears maintains performance in DAG estimation
across increasing levels of non-linearity, and on average outperforms all
methods. dotears performs worse in the “imperfect” intervention setting. In
particular, dotears has average performance in “Low Density” scenarios, but
performs much worse in “High Density” simulations. This suggests that the
imperfect intervention model, rather than non-linearity, drives the relatively
worse performance of dotears. As in Figure 16(a), dense DAGs exacerbate
difficulties with imperfect interventions. However, on average dotears
outperforms the neural network DCDI even under “imperfect” interventions.

(a) Performance under soft interventions, where $\pi$ is the proportion of
allowed parental variance to the target.

(b) Performance under random per-node effects by interventions on target error
variance.

(c) Performance under interventions which fix the distribution of the target
to $\mathcal{N}(2,1)$.

Figure 16: Method performance on sensitivity analyses aggregated across
parameterizations of $p=20$ random graphs, in Structural Hamming Distance
(lower is better). See Supplementary Material 7.8 for details.

###  7.8 Sensitivity under linear SEM

For evaluation of dotears sensitivity to violations of the modeling
assumptions, data is generated for DAGs with $p=20$ nodes. Structures are
drawn from both Erdős-Rényi (ER) and Scale-Free (SF) DAGs [24, 25]. We
simulate under two parameterizations $\\{\text{Low Density, High Density}\\}$
of the edge densities $(r,z)$. In the “Low Density” parameterization, give
$(r,z)=(0.2,2)$. To evaluate performance on higher density topologies, we also
give “High Density” parameterizations, where $(r,z)=(0.5,6)$.

Given an edge density scenario, we simulate under two parameterizations of the
edge weights. In the “Strong Effects” parameterization,
$w\sim\text{Unif}\left([-1.0,-0.8]\cup[0.8,1.0]\right)$. We also give the
“Weak Effects” parameterization, which edge weights are drawn from
$w\sim\text{Unif}\left([-0.3,-0.5]\cup[0.3,0.5]\right)$. Table 1 summarizes
all four possible simulation parameterizations.

For each node $i$, we draw $\sigma_{i}\sim\text{Unif}\left([0.5,2.0]\right)$,
and draw $n_{0}$ observations from
$\boldsymbol{\epsilon}_{i}^{(0)}\sim\mathcal{N}\left(0,\sigma_{i}^{2}\right)$.
For each DAG we generate an instance of observational data, where
$n_{0}=(p+1)*n_{k}=2100$, and an instance of interventional data, where
$n_{k}=100$ for all $k=0\dots p$ to match sample size. We set the distribution
of $\epsilon_{i}^{(k)}$ according to Assumptions 3 and 2 for $\alpha=4$.
Cross-validation was performed as in Supplementary Material 7.6.1.

Simulation Type | Topology | $w$
---|---|---
Low Density | , Strong Effects | ER-0.2 | OR SF-2 | $\text{Unif}\left([-1.0,-0.8]\cup[0.8,1.0]\right)$
Low Density | , Weak Effects | ER-0.2 | OR SF-2 | $\text{Unif}\left([-0.5,-0.3]\cup[0.3,0.5]\right)$
High Density | , Strong Effects | ER-0.5 | OR SF-6 | $\text{Unif}\left([-1.0,-0.8]\cup[0.8,1.0]\right)$
High Density | , Weak Effects | ER-0.5 | OR SF-6 | $\text{Unif}\left([-0.5,-0.3]\cup[0.3,0.5]\right)$
Table 3: For each simulation type for $p=20$ sensitivity analyses, we provide
the parameterizations of the Erdős-Rényi and Scale-Free topologies and the
distribution on the drawn edge weights.

####  7.8.1 Hard Interventions

Figure 17: Method performance under soft or imperfect interventions, where
$\pi$ is the proportion of allow parental variance to the target. Each row is
a different parameterization of DAG density and edge weights; each column is a
different value of $\pi$. Results are shown in Structural Hamming Distance
(SHD), and methods are sorted by average performance in SHD. Figure 18: Method
performance under soft or imperfect interventions, where $\pi$ is the
proportion of allow parental variance to the target. Each row is a different
parameterization of DAG density and edge weights; each column is a different
value of $\pi$. Results are shown in precision of binary edge recovery, and
methods are sorted by average performance in SHD. Figure 19: Method
performance under soft or imperfect interventions, where $\pi$ is the
proportion of allow parental variance to the target. Each row is a different
parameterization of DAG density and edge weights; each column is a different
value of $\pi$. Results are shown in recall of binary edge recovery, and
methods are sorted by average performance in SHD.

Assumption 1 assumes that interventions remove all edges incoming to the
target, and therefore that the marginal variance of the target has no incoming
parental variance. We examine the model where

$X_{k}^{(k)}=\sqrt{\pi}\sum_{i=0}^{p}w_{ij}X_{i}^{(k)}+\epsilon_{k}^{(k)},$

where interventions downweight incoming parental edges by $\sqrt{\pi}$.
Accordingly, this means that

$\text{Var}\left(X_{k}^{(k)}\right)=\pi\sum_{i=0}^{p}w_{ij}^{2}\text{Var}\left(X_{i}^{(k)}\right)+\frac{\sigma_{k}^{2}}{\alpha^{2}},$

where for $\pi=0$ we recover hard interventions. We call these interventions
either soft or imperfect interventions, and examine method performance when
$\pi=\\{0,0.1,0.25,0.5\\}$. Note that other than the imperfect intervention,
the model and generating processes are the same as in Supplementary Material
7.6.1.

For each $\pi\in\\{0,0.1,0.25,0.5\\}$ and each parameterization in Table 3,
ten Erdős-Rényi DAGs and ten Scale-Free DAGs are drawn, with sample size
matched as described above. Results in Structural Hamming Distance are shown
in Supplementary Figure 17 per parameterization, and are shown averaged across
parameterizations in Figure 16(a). Results in binary precision and recall are
given in Supplementary Figure 18 and Supplementary Figure 19, respectively.

####  7.8.2 $\alpha$ perturbation

Assumption 3 assumes that for each node under intervention, the error variance
is reduced by a shared factor $\alpha$ across interventional targets, i.e.

$\text{Var}\left(\epsilon_{k}^{(k)}\right)=\frac{\sigma_{k}^{2}}{\alpha^{2}}$

for all $k$.

To test the sensitivity of dotears to this assumption, we draw data for which

$\text{Var}\left(\epsilon_{k}^{(k)}\right)=\frac{\sigma_{k}^{2}}{c_{k}^{2}\alpha^{2}},\quad
c_{k}\sim\text{Unif}([0.8,1.2]).$

$c_{k}$ is target-specific, and thus under this model $\Omega_{0}$ will be
misspecified.

Ten Erdős-Rényi DAGs and ten Scale-Free DAGs are drawn, with sample size
matched as described above. Results in Structural Hamming Distance, averaged
across parameterizations, are shown in Figure 16(b). Results for each
parameterization in Structural Hamming Distance are shown in Supplementary
Figure 20 (top row). Results in binary edge precision are shown in
Supplementary Figure 21. Results in binary edge recall are shown in
Supplementary Figure 22.

Figure 20: Method performance in linear sensitivity analyses. Each row is a
different sensitivity analysis; each column is a parameterization of DAG
density and edge weights. Results are shown in Structural Hamming Distance
(SHD), and methods are sorted by average performance in SHD. Figure 21: Method
performance under linear sensitivity analyses. Each row is a different
sensitivity analysis; each column is a parameterization of DAG density and
edge weights. Results are shown in precision of binary edge recovery, and
methods are sorted by average performance in SHD. Figure 22: Method
performance under linear sensitivity analyses. Each row is a different
sensitivity analysis; each column is a parameterization of DAG density and
edge weights. Results are shown in recall of binary edge recovery, and methods
are sorted by average performance in SHD.

####  7.8.3 Fixed Interventions

Our model assumes that the distribution of node $k$ under intervention on $k$
is dependent on the distribution of $\epsilon_{k}^{(0)}$. We also examine an
interventional model in which intervention on a node $k$ sets the distribution
with a fixed mean shift and error variance, i.e.

$X_{k}^{(k)}\sim\mathcal{N}(2,1),$

or equivalently $\epsilon_{k}^{(k)}$. This model is used in linear simulations
in [22] and is similar to that of [10].

###  7.9 Sensitivity to nonlinear SEM

To evaluate the performance of dotears on nonlinear data, we generate data
from the generative models in [22]. Specifically, structures are Erdős-Rényi
random DAGs with $p=10$ with single-node interventions on every node. The
total sample size is $n=10000$ for all draws and scenarios, evenly split
across the interventional scenarios (10 interventions and 1 observational) for
$n_{k}\approx 909$.

Data is generated either from the additive noise model (ANM) or the non-linear
with non-additive noise model (NN). For a node $j$, let
$N_{j}\sim\mathcal{N}(0,\sigma_{j}^{2})$, where
$\sigma_{j}^{2}\sim\text{Unif}[1,2]$. In the additive noise model, for a node
$j$ we have

$X_{j}\coloneqq f_{j}\left(\text{Pa}(j)\right)+0.4\cdot N_{j},$ (37)

where “$f_{j}$ are fully connected neural networks with one hidden layer of 10
units and leaky ReLU with a negative slope of 0.25 as nonlinearities” [22].
Imperfect interventions add a random vector of $\mathcal{N}(0,1)$ to the last
layer; for details, see [22].

In the nonlinear with non-additive noise model, for a node $j$ we have

$X_{j}\coloneqq f_{j}\left(\text{Pa}(X_{j}),N_{j}\right),$ (38)

where $f_{j}$ are “fully connected neural networks with one hidden layer of 20
units and tanh as nonlinearities” [22]. Imperfect interventions are obtained
in the same fashion as the additive noise model; for details, see [22].

We give two parameterizations for each setting: “Low Density”, in which the
expected number of edges is 10, and “High Density”, in which the expected
number of edges is 40. Results in SHD are given in Figure 15. Results in
precision and recall are given in Figure 23 and Figure 24, respectively.

Figure 23: Method performance in nonlinear sensitivity analyses. Each row is a
different sensitivity analysis; each column is a parameterization of DAG
density and edge weights. Results are shown in precision of binary edge
recovery, and methods are sorted by average performance in SHD across
parameterizations. Figure 24: Method performance in nonlinear sensitivity
analyses. Each row is a different sensitivity analysis; each column is a
parameterization of DAG density and edge weights. Results are shown in recall
of binary edge recovery, and methods are sorted by average performance in SHD
across parameterizations.

###  7.10 Genome-Wide Perturb-Seq

We benchmark all methods on a single-cell Perturb-seq experiment from Replogle
et al. [5]. Replogle et al. provide both raw count data as well as normalized
data, where normalized here means a z-scoring relative to the mean and
standard deviation of the control dataset, accounting for batch. Here, the
control dataset indicates a set of pre-selected control cells, and not simply
cells with non-targeting guides. For details, see [5].

We perform feature selection in the raw data. We select the top 100 most
variable genes in the raw observational data, excluding 1.) genes coding for
ribosomal proteins and 2.) genes without knockdown data. Here, the
observational data indicates cells incorporating non-targeting control guides.
We then take the normalized expression data for these selected 100 genes in
1.) the observational data, and 2.) each of the 100 knockdowns. This forms our
training set in Figures 4a-d.

In Figure 4e, we perform the same procedure as above, but exclude 1.) the
observational data. In other words, our training set is formed exclusively
from expression data in the 100 knockdowns.

Cross-validation is not performed due to low sample sizes in some knockdowns;
instead, the L1 penalty is arbitrarily set to 0.1 for all methods where
appropriate. Otherwise, method settings are the same as simulations.

For DCDI, we report intermediate results, since convergence was not obtained
under 24 hours even on GPU training. CPU training was attempted with 30 cores
and 180 GB memory, but was killed by the operating system. GPU training was
attempted, but was also repeatedly killed by the operating system for memory
constraints. DCDI is run as DCDI-G, since DCDI-DSF would not fit in memory,
and is also run with imperfect interventions and known interventions. DCDI
reported 2982 total edges, but many of these are both causal and anti-causal;
in total, 2039 “interactions” were reported.

####  7.10.1 Differential Expression Testing

We use DESeq2 to test for differential expression [26, 27]. We calculate size
factors across the raw feature-selected data, and use these in all downstream
tests. Next, for each knockdown $\text{ko}(i)$, we test for differential
expression for all genes $j$ compared to the observational data. For single-