Estimation Rates for Sparse Linear Cyclic Causal Models

TL;DR

This paper analyzes the statistical estimation rates for sparse linear cyclic causal models, introducing a new estimator and comparing it to existing methods, with theoretical and empirical insights.

Contribution

It presents a novel two-step penalized maximum likelihood estimator and establishes its near minimax optimality for sparse cyclic causal graphs.

Findings

The new estimator outperforms LLC in synthetic experiments.

Asymptotic near minimax optimality is proven for the maximum likelihood estimator.

Practical advantages are demonstrated through numerical experiments.

Abstract

Causal models are important tools to understand complex phenomena and predict the outcome of controlled experiments, also known as interventions. In this work, we present statistical rates of estimation for linear cyclic causal models under the assumption of homoscedastic Gaussian noise by analyzing both the LLC estimator introduced by Hyttinen, Eberhardt and Hoyer and a novel two-step penalized maximum likelihood estimator. We establish asymptotic near minimax optimality for the maximum likelihood estimator over a class of sparse causal graphs in the case of near-optimally chosen interventions. Moreover, we find evidence for practical advantages of this estimator compared to LLC in synthetic numerical experiments.

Equations

$$X=B^{\ast }X+Z,Z\sim \mathcal{N}(0,I),$$

$$d_H(x,y)=\{i\in [p]:x\ne y\}$$

$$\Vert B{\Vert }_{\infty }=\underset{i,j}{\max }|B_{i,j}|,\Vert B{\Vert }_1=\sum _{i,j}|B_{i,j}|,$$

$$\Vert B{\Vert }_{\infty ,\infty }=\underset{i\in [p_1]}{\max }\Vert B_{i,:}{\Vert }_1.$$

$$\mathcal{B}:=\{B\in \mathrm{I}{\mathrm{R}}^{p\times p}:B_{i,i}=0,\text{ for all }i=1,\dots ,p\}.$$

$$X=\sum _{k=0}^{\infty }(B^{\ast }{)}^{k}Z,$$

$$x_{t+1}=B^{\ast }{x}_t,x_0=Z.$$

$${\Sigma }^{\ast }=(I-B^{\ast }{)}^{-1}(I-B^{\ast }{)}^{-\top }.$$

$$d(B)=\underset{i\in [p]}{\max }|\{j:B_{i,j}\ne 0\}|,$$

$$[p]={\mathcal{U}}_e\dot{\cup }{\mathcal{J}}_e,$$

$$(U_e{)}_{i,j}=\{\begin{array}{rl}1,& i=j\text{ and }i\in {\mathcal{U}}_e\\ 0,& \text{otherwise},\end{array}(J_e{)}_{i,j}=\{\begin{array}{rl}1,& i=j\text{ and }i\in {\mathcal{J}}_e\\ 0,& \text{otherwise}.\end{array}$$

$${\Sigma }^{\ast ,e}=(I-U_e{B}^{\ast }{)}^{-1}(I-U_e{B}^{\ast }{)}^{-\top },$$

$${\Theta }^{\ast ,e}=({\Sigma }^{\ast ,e}{)}^{-1}=(I-U_e{B}^{\ast }{)}^{\top }(I-U_e{B}^{\ast }).$$

$$\mathcal{B}(p,d,\eta ):=\{B\in \mathrm{I}{\mathrm{R}}^{p\times p}:B_{i,i}=0\text{ for }i\in [p],\Vert B{\Vert }_{\mathrm{op}}\le 1-\eta ,d(B)\le d\},$$

$$\kappa =\kappa (\mathcal{E})=\underset{i\ne j\in [p]}{\max }|\{e\in \mathcal{E}:i\in {\mathcal{U}}_e,j\in {\mathcal{J}}_e\}|.$$

$$n\ge pd{E}^2\log \left(1+\frac{p}{4d}\right),$$

$$\inf_{\hat{B}}\sup_{B^{\ast}\in\mathcal{B}(p,d,\eta)}{\rm I}\kern-1.79993pt{\rm P}\big{(}\|\hat{B}-B^{\ast}\|_{F}^{2}\gtrsim\frac{pdE}{\kappa n}\log\big{(}1+\frac{p}{4d}\big{)}\big{)}\geq c,$$

$$(T_i^{\ast }{)}_{\ell ,:}={\mathfrak{e}}_j^{\top }{\Sigma }^{\ast ,e}{P}_i^{\top },(t_i^{\ast }{)}_{\ell }={\Sigma }_{j,i}^{\ast ,e}$$

$$(T_i^{\ast }{)}_{\ell ,:}=\left[\begin{array}{cccccc}0& \dots & 1& \dots & 0& {\Sigma }_{j,{\mathcal{U}}_e\setminus \{i\}}^{\ast ,e}\end{array}\right],$$

$${\hat{\Sigma }}^e=\frac{1}{n_e}\sum _{k=1}^{n_e}{X}_k^e(X_k^e{)}^{\top },$$

$${\hat{b}}_i={\mathrm{argmin}}_{b\in \mathrm{I}{\mathrm{R}}^{p-1}}\Vert {\hat{T}}_{i}b-{\hat{t}}_i{\Vert }_2^2+\lambda \Vert b{\Vert }_1,i=1,\dots ,p,$$

$$({\hat{B}}_{\mathrm{llc}}{)}_{i,-i}={\hat{b}}_i^{\top },({\hat{B}}_{\mathrm{llc}}{)}_{i,i}=0,i\in [p].$$

$$\mathcal{C}(d):=\{v\in \mathrm{I}{\mathrm{R}}^p:\text{for all }S\subseteq [p]\text{ with }|S|\le d,\Vert v_{S^c}{\Vert }_1\le 3\Vert v_S{\Vert }_1\}.$$

$$\rho (d)=$$

$$n\gtrsim$$

$$\lambda \asymp$$

$$\Vert {\hat{B}}_{\mathrm{llc}}-B^{\ast }{\Vert }_F^2\lesssim$$

$$\tilde{\rho }(d):=\underset{i\in [p]}{\min }\underset{\begin{array}{c}v\in \mathrm{I}{\mathrm{R}}^p,v\ne 0,\\ |\mathrm{supp}v|\le d\end{array}}{\inf }\frac{\Vert T_i^{\ast }v{\Vert }_2}{\Vert v{\Vert }_2}.$$

$$\ell (\Theta ,\hat{\Sigma })=\mathrm{Tr}(\hat{\Sigma }\Theta )-\log \det (\Theta )$$

$$\mathcal{L}(B)=\mathcal{L}(B,{\hat{\Sigma }}^1,\dots ,{\hat{\Sigma }}^E)=\sum _{e\in \mathcal{E}}\ell ({\Theta }^e(B),{\hat{\Sigma }}^e),$$

Full text

t2Supported by NSF awards IIS-1838071, DMS-1712596 and DMS-TRIPODS-1740751; ONR grant N00014-17- 1-2147 and grant 2018-182642 from the Chan Zuckerberg Initiative DAF.

1 Introduction

Directed graphical models [Pea09, SGS00] provide a useful framework for interpretation, inference and decision making in many areas of science such as biology, sociology, and environmental sciences [FLNP00, Dun66, KH88]. Unlike their undirected counterparts that merely encode the structure of probabilistic dependence between random variables directed graphical models reveal causal effects that are the basis of scientific discovery  [Pea09].

Most frequently, the model is assumed to be governed by a directed acyclic graph (DAG) $G=(V,E)$, where $V=\{X_{1},\dots,X_{p}\}$ are the variables of an observed system and $E$ is a set of edges such that there is no directed cycle in $G$. In such models, known as Bayes networks [Pea09], the variables follow a joint distribution that factorizes according to the graph $G$ in the sense that node $i$ is independent of other nodes conditionally on its parents. The absence of cycles allows for a direct interpretation of the causal structure between the variables $X_{1},\ldots,X_{p}$ whereby a directed edge corresponds to a causal effect. At the same time, most complex systems showcase feedback loops that can be both positive and negative and the need to extend Bayes networks to allow for cycles was recognized long ago.

A large body of work focuses on learning Bayes from observational data, that is, data drawn independently from the joint distribution of $(X_{1},\dots,X_{p})$. Observational data is rather abundant but even in the cyclic cases, it is known to lead to a severe lack of identifiability: Such data, even in infinite abundance, can only yield an equivalence class—the Markov equivalence class—of DAGs that are all compatible with the conditional independence relation in the given data. While a DAG in the Markov equivalence class can already yield decisive scientific insight [MKB09], searching over the space of DAGS is often computationally hard. Many algorithms have been proposed over the years such as the PC algorithm [SGS00] and Greedy Equivalence search [Chi02] and max-min hill climbing [TBA06], but all of them rely on the notion of faithfulness of the distribution, i.e., the assumption that all conditional dependence relations that could be compatible with the DAG $G$ are actually fulfilled by the distribution of $X$. In fact, for consistency of these algorithms, one needs to assume that these dependencies observe a signal-to-noise ratio that allows to detect them with high probability [KB07, LB14, vdGB13]. Extensions that allow certain kinds of cycles, [Ric96, RS96, SM09, IOS*+*10, LSRH12] have been proposed but at the expense of having an increased number of graphs in each equivalence class.

Recent breakneck advances in data collection processes such as the spread of A/B testing for online marketing or targeted gene editing with CRISPR-Cas9 are contributing to the proliferation of interventional data, the gold standard for causal inference. With unlimited interventions on any combination of nodes, learning a directed graphical model becomes a trivial task. However, exhaustively performing all interventions is a daunting and costly task and recent work has focused on finding a small number of interventions for several classes of DAGs [SKDV15, KDV17]. For graphs with cycles, [HEH12] have characterized the system of interventions necessary to learn a parametric linear structural equation model (SEM) [BH77, Bol83], in which all variables are real valued and the causal relationships given by the edges $E$ are linear. Formally, this model postulates that the following equation holds (in distribution) for observational samples from $X$:

[TABLE]

where we exclude explicit self-loops by assuming that the diagonal of $B^{\ast}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ is zero. By writing $X=(I-B^{\ast})^{-1}Z$ and assuming that the corresponding inverse matrix exists, this allows us to handle underlying graphs that are cyclic. This model has been extensively studied in [HEH12], where it is shown that if we have access to data from a sufficiently rich system of interventions, i.e., if enough variables are randomized and are thus made independent of the influence of their parents encoded in $B^{\ast}$, then on a population level, $B^{\ast}$ is identifiable by a method of moments type estimator that the authors call LLC (for linear, latent, causal).

In this paper, we present upper and lower bounds for the reconstruction of $B^{\ast}$ in Frobenius norm for classes of sparse $B^{\ast}$, corresponding to graphs with bounded in-degree, using multiple observations for each intervention setup. We also provide upper bounds for the original LLC estimator with $\ell_{1}$-penalization term as well as an $\ell_{1}$-penalized maximum likelihood estimator, all under the simplifying assumption that the noise or disturbance variables $Z$ are Gaussian, independent of each other, and have unit variance. Moreover, we provide numerical evidence that a non-convex ADMM type algorithm can be used to find a solution to this maximum likelihood problem, albeit without convergence guarantees.

1.1 Related work

It is known that several variants of the model (1.1) are identifiable from observational data, including non-linear SEMs [HJM*+*09] or non-Gaussian noise [SHHK06]. Linear SEMs with Gaussian noise can be identifiable under additional assumptions, for example when the components of the noise have equal variances and the underlying graph is a DAG [LB14, PB14] or when the underlying graph is random and sparse [AR18].

As for assumptions on the noise that guarantee identifiability from observational data, one example is the recovery of a linear structural equation model under non-Gaussian noise via Independent Component Analysis, [SHHK06], which under additional restrictive assumptions on the structure of the underlying graph can be extended to the Gaussian case, [AR18].

Moreover, many more approaches to dealing with cycles and/or interventions are known, such as convex regularizers in an exponential family model [SNM07, SM09], independence testing [IOS*+*10], Independent Component Analysis [LSRH12], and adapting Greedy Equivalence Search to handle interventional data [HB12, WSU18]. From the above, it seems that the linear Gaussian case is somewhat of a worst case example for identifiability of the ground truth matrix, especially when allowing cycles, and thus warrants the investigation of controlled interventions to eliminate ambiguity, which is the main contribution of [HEH12]. Similar models have been considered for applications, for example in computational biology, see [CBG13], where identifiability is not provided by controlled experiments on the variance, but rather by a mean shift of some variables.

Our work extends the results in [HEH12] by providing explicit upper bounds for their suggested method, as well as presenting an alternative estimator that leads to upper bounds independent of the conditioning of the experiments as explained in Section 3.3. In spirit, our results are similar to consistency guarantees obtained in [vdGB13] and [WSU18], but we focus on the case where enough interventions are performed to identify the ground truth structure matrix $B^{\ast}$, alleviating the need for additional assumptions on $B^{\ast}$.

1.2 Structure of the paper

The rest of the paper is structured as follows: In Section 2, we give an overview of the linear structural equation model we consider and the main assumptions we make. In Section 3, we present lower bounds, upper bounds for LLC, and upper bounds for a two-step maximum likelihood estimator. In Section 4, we derive a non-convex variant of ADMM to solve part of the numerical optimization problem for the penalized maximum likelihood estimator and explore its performance on synthetic and semi-synthetic data. The proofs of the main results are deferred to Sections C – E in the appendix, and we collect general lemmas used in all the proofs in Section F. Section A contains a short argument for why experimental data is necessary given our assumptions, and Section B provides a way of speeding up our numerical calculations.

1.3 Notation

We write $a\lesssim b$ for two quantities $a$ and $b$ if there exists an absolute constant $C>0$ such that $a\leq Cb$, and similarly for $a\gtrsim b$. For a natural number $p$, we denote by $[p]=\{1,\dots,p\}$. Given a set $S$, we write $|S|$ for its cardinality.

Let $x,y\in{\rm I}\kern-1.79993pt{\rm R}^{p}$. We write $\operatorname{supp}x$ for the indices of non-zero elements of $x$,

[TABLE]

for the Hamming distance between $x$ and $y$, and $\|x\|_{p}$ for the $\ell^{p}$ norm of $x$.

For two matrices $A,B\in{\rm I}\kern-1.79993pt{\rm R}^{p_{1}\times p_{2}}$, we abbreviate the $i$th row by $B_{i,:}$ and the $i$th column by $B_{:,i}$. Similarly, $B_{i,-j}$ denotes the $i$th row of $B$ where the $j$th element is omitted. Further, $\|B\|_{F}$ denotes the Frobenius norm, $\|B\|_{\mathrm{op}}$ the operator norm,

[TABLE]

and $\|B\|_{\infty,\infty}$ the operator norm of $B$ with respect to the $\ell^{\infty}$ norm, which is

[TABLE]

If $A$ is a square invertible matrix, we denote by $A^{-1}$ its inverse and by $A^{-\top}$ the transpose of $A^{-1}$. We denote the smallest and largest singular value of $A$ by $\sigma_{\min}(A)$ and $\sigma_{\max}(A)$, respectively. If $A$ and $B$ are symmetric, we write $A\prec B$ if $B-A$ is positive definite, and similarly for $A\succ B$. By $I\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$, we denote the identity matrix.

For a function $f:{\rm I}\kern-1.79993pt{\rm R}^{p_{1}}\to{\rm I}\kern-1.79993pt{\rm R}^{p_{2}}$, we denote its derivative at a point $x\in{\rm I}\kern-1.79993pt{\rm R}^{p_{1}}$ applied to a vector $h\in{\rm I}\kern-1.79993pt{\rm R}^{p_{1}}$ by $Df(x)[h]$. We write $\mathsf{subG}$ and $\mathsf{subE}$ to denote sub-Gaussian and sub-Exponential distributions as defined in Definition 24.

2 Model and assumptions

Before summarizing our explicit assumptions, we give a definition of observations under a linear cyclic structural equation model with and without interventions. We assume that a linear SEM on a random vector $X=(X_{1},\dots,X_{p})$ is given by a matrix $B^{\ast}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ without self-cycles, i.e., $B^{\ast}\in\mathcal{B}$ with

[TABLE]

Without any intervention, each observation is an independent copy of $X=(I-B^{\ast})^{-1}Z$, where $Z$ can in principle be any noise variable. Since non-Gaussian noise can lead to identifiability from observational data through exploiting this particular property [HJM*+*09, LSRH12], we focus on Gaussian noise, and make the simplifying assumption that $Z\sim\mathcal{N}(0,I)$. In order to guarantee that $(I-B^{\ast})^{-1}$ exists, we assume $\|B^{\ast}\|_{\mathrm{op}}<1$ which in particular allows us to write

[TABLE]

and $X$ can be interpreted as the steady state distribution of an auto-regressive process $\{x_{t}\}_{t\geq 0}$ governed by the dynamics

[TABLE]

Hence, $X$ is distributed according to $X\sim\mathcal{N}(0,\Sigma^{\ast})$ with

[TABLE]

In order to obtain results in the high-dimensional regime $p\asymp n$, we additionally assume that the in-degree of $B^{\ast}$ is bounded, resulting in a sparse matrix $B^{\ast}$. That is, if we denote the maximum in-degree of a matrix $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ by

[TABLE]

then we assume $d(B^{\ast})\ll p$.

Moreover, we assume that we have access to interventional, a.k.a. experimental, data, which is modeled as follows, keeping in line with the definition from [HEH12]. An experiment $e$ is given by a partition

[TABLE]

with associated projection matrices

[TABLE]

In effect, all nodes in $\mathcal{J}_{e}$ are intervened on, i.e., they are not influenced by their parents anymore. We assume that they follow a standard Gaussian distribution $\mathcal{N}(0,1)$, leading to a random variable $X^{e}\sim\mathcal{N}(0,\Sigma^{\ast,e})$ corresponding to experiment $e$ with covariance matrix

[TABLE]

and inverse covariance matrix (concentration matrix)

[TABLE]

[HEH12] provide the following criterion to identify $B^{\ast}$ from interventional data associated with $\mathcal{E}$.

Definition 1** (Completely separating system).**

The set of experiments $\mathcal{E}$ is a completey separating system if for every $i\neq j\in[p]$, there exists $e\in\mathcal{E}$ such that $i\in\mathcal{J}_{e}$ and $j\in\mathcal{U}_{e}$.

Note that [HEH12] call the separation condition for a pair $(i,j)\in[p]^{2}$ the pair condition. They show that Definition 1 guarantees identifiability of $B^{\ast}$ from observational data. Conversely, they show that if $\mathcal{E}$ is not separating, there exists a ground truth system that is not satisfied, albeit allowing a more general covariance structure on the error terms $Z^{e}_{k}$ for the latter construction than we do.

We are now in a position to state our assumptions.

A1****Structure matrix.

For any two positive integers $d\leq p$ and $\eta\in(0,1/2]$, let $\mathcal{B}(p,d,\eta)$ denote the set of sparse matrices defined by

[TABLE]

and assume $B^{\ast}\in\mathcal{B}(p,d,\eta)$.

A2****Interventions.

Let $\mathcal{E}$ be a set of experiments with associated partitions $\{(\mathcal{U}_{e},\mathcal{J}_{e})\}_{e\in\mathcal{E}}$ and projection matrices $\{(U_{e},J_{e})\}_{e\in\mathcal{E}}$ as in (2.4) and (2.5), respectively. Assume that $\mathcal{E}$ is separating in the sense of Definition 1.

A3****Noise.

*Assume $n\in\mathbb{N}$ is divisible by $E:=|\mathcal{E}|$, set $n_{e}=n/E$ for $e\in\mathcal{E}$, and for $k\in[n_{e}],e\in\mathcal{E}$, denote by $Z_{k}^{e}\sim\mathcal{N}(0,I)$ i.i.d. Gaussian random vectors. Then, we assume that we have access to observations of the form $X_{k}^{e}=(I-U_{e}B^{\ast})^{-1}Z_{k}^{e}$. *

A few remarks are in order.

A1. The bound $\|B^{\ast}\|_{\mathrm{op}}\leq 1-\eta$ guarantees invertibility of $I-UB^{\ast}$ for any projection matrix $U$ and stationarity of the process (2.2).

A2. As mentioned, this is the same assumption under which [HEH12] show identifiability of $B^{\ast}$ under more general assumptions than the ones presented here, in particular allowing more general noise variances and hidden variables. Note that their proof of necessity of this assumption does not exactly match our assumption because our noise variances are restricted, so in principle, identifiability from observational data could be possible under a weaker condition. However, we give evidence in Section A that at least observational data alone is not sufficient to recover a general $B^{\ast}$.

Intuitively, the fact that $\mathcal{E}$ is separating guarantees that $B^{\ast}$ can be recovered from submatrices of $\{\Sigma^{\ast,e}\}_{e\in\mathcal{E}}$ via solving a system of linear equations, a fact that is made more precise in Section 3.2. Since we are interested in recovering $B^{\ast}$ under otherwise minimal assumptions on $B^{\ast}$, this is the case we consider for the theoretical contributions of this paper. We do however investigate the behavior of the two estimators considered in Section 3 with respect to a violation of this assumption numerically in Section 4.

A3. The assumption of Gaussian noise is not critical for our analysis, and in fact all our proofs extend readily to sub-Gaussian noise. Similarly, the assumption $n_{e}=n/E$ can be replaced by $n_{e}\asymp n/E$, that is, the number of observations in all experiments is comparable. On the other hand, the assumption that ${\rm I}\kern-1.79993pt{\rm E}[Z_{k}^{e}]=0$, ${\rm I}\kern-1.79993pt{\rm E}[(Z_{k}^{e})^{2}]=1$ might be restrictive in practice. We conjecture that it might be relaxed while maintaining many of the guarantees we give in Section 3, but due to the notational burden associated with incorporating these additional factors into the estimation, we chose to leave this topic as the subject of future research. Note that while the assumption of equal variances implies identifiability from observational data in the case where $B^{\ast}$ is assumed to be acyclic [LB14, PB14], it does not in the cyclic case, see Section A. Hence, the assumptions as presented still lead to a class rich enough to require controlled experiments to estimate $B^{\ast}$.

Remark 2**.**

It was shown in [HEH13] that the minimum number of experiments necessary to obtain a completely separating system is of the order $\log(p)$, which can be seen by a simple binary coding argument. Hence, if we are able to pick the experiments in the most parsimonious way possible, $E=O(\log(p))$ only contributes a logarithmic factor to any of the rates presented in Section 3.

3 Main results

3.1 Lower bounds

First, we give lower bounds for the estimation of matrices $B^{\ast}\in\mathcal{B}(p,d,\eta)$. To that end, let $\kappa$ denote the redundancy of the experiments $\mathcal{E}$. It is defined as the maximum number of experiments that separate two variables,

[TABLE]

Theorem 3**.**

There exists a constant $c>0$ such that if $d\leq p/4$ and

[TABLE]

then,

[TABLE]

where the infimum is taken over all measurable functions of the data $\{X_{k}^{e}\}_{e\in\mathcal{E},\,k\in[n_{e}]}$.

The proof of Theorem 3 is deferred to Section C. We remark that there is a mismatch in the lower bound and the range of $n$ for which it is effective that is of order $E$. In the case of a minimal system of completely separating interventions, by Remark 2, this mismatch is of order $\log(p)$.

3.2 Upper bounds for the LLC estimator

Next, we give bounds on the performance of the LLC estimator introduced in [HEH12]. We briefly summarize the algorithm below, which can be seen as a moment estimator for $B^{\ast}$.

3.2.1 The LLC estimator

Denote by $b^{\ast}_{i}\in{\rm I}\kern-1.79993pt{\rm R}^{p-1}$ the $i$th row of $B^{\ast}$, where we omit the $i$th entry, which is assumed to be zero since $B^{\ast}\in\mathcal{B}$. Formally, $b^{\ast}_{i}=(P_{i}B_{i,:}^{\top})=B_{i,-i}^{\top}$, where $P_{i}\colon{\rm I}\kern-1.79993pt{\rm R}^{p}\to{\rm I}\kern-1.79993pt{\rm R}^{p-1}$ denotes the projection operator that omits the $i$th coordinate.

LLC is motivated by the observation that on the population level, each $b^{\ast}_{i}$ satisfies a linear system $T^{\ast}_{i}b^{\ast}_{i}=t^{\ast}_{i}$, where $T^{\ast}_{i}\in{\rm I}\kern-1.79993pt{\rm R}^{m_{i}\times(p-1)}$ and $t^{\ast}_{i}\in{\rm I}\kern-1.79993pt{\rm R}^{m_{i}}$ for some $m_{i}\geq 1$ are defined as follows. For $i=1,\dots,p$, define the matrix $T^{\ast}_{i}$ and the column vector $t^{\ast}_{i}$ row by row. For each experiment $e$ such that $i\in\mathcal{U}_{e}$ and each $j\in\mathcal{J}_{e}$, add a row to $T^{\ast}_{i}$ and to $t^{\ast}_{i}$, say with index $\ell=\ell(e,j)$, that is of the form

[TABLE]

where $\mathfrak{e}_{j}$ is the $j$th canonical vector of ${\rm I}\kern-1.79993pt{\rm R}^{p}$. To better visualize $(T^{\ast}_{i})_{\ell,:}$, one may rearrange the indices so that $\mathcal{J}_{e}=\{1,\dots,|\mathcal{J}_{e}|\}$, in which case we have

[TABLE]

where “$1$” appears in the $j$th coordinate. Let $m_{i}$ denote the total number of such rows obtained by scanning through all experiments $e$ such that $i\in\mathcal{U}_{e}$ and $j$ such that $j\in\mathcal{J}_{e}$.

When $\mathcal{E}$ is a completely separating system, $T^{\ast}_{i}b_{i}=t^{\ast}_{i}$ has the unique solution $b^{\ast}_{i}=(B^{\ast}_{i,-i})^{\top}$, [HEH12]. The LLC estimator is obtained by substituting $\Sigma^{\ast,e}$ in the above definitions with its empirical counterpart $\hat{\Sigma}^{e}$ defined by

[TABLE]

except for where the variances are known exactly due to the fact that an intervention is performed. This leads to a linear system of the form $\hat{T}_{i}b_{i}=\hat{t}_{i}$. Rather than solving the linear system exactly, the LLC estimator is obtained by minimizing a penalized least squares problem to promote sparsity in the resulting estimate:

[TABLE]

where $\lambda>0$ is a tuning parameter. The solutions to the above problems are assembled into the LLC estimator $\hat{B}_{\mathrm{llc}}$ by setting

[TABLE]

3.2.2 Statistical performance

The upper bounds we give for the performance of LLC depend on additional constants that are not directly controlled for an arbitrary $B^{\ast}\in\mathcal{B}(p,d,\eta)$. Loosely speaking, they pertain to the conditioning of the $\ell^{1}$-regularized least squares problems that are solved to obtain $\hat{B}_{\mathrm{llc}}$. These constants are defined as follows. Denote by

[TABLE]

Then, define

[TABLE]

We are now in a position to state the first rate of convergence for the LLC estimator.

Theorem 4** (Rates for LLC estimator).**

Let assumptions A1 – A3 hold and fix $\delta\in(0,1)$. Assume further that

[TABLE]

Then LLC estimator $\hat{B}_{\mathrm{llc}}$ defined in (3.7) with $\lambda$ chosen such that

[TABLE]

satisfies

[TABLE]

*with probability at least $1-\delta$. *

The proof is deferred to Section D. It uses standard arguments for the LASSO, together with perturbation results for regression with noisy design from [LW11] in Lemma 16 to handle the presence of noise in the matrices $\hat{T}_{i}$.

Remark 5**.**

Unfortunately, it is not clear whether the factors $\rho(d),R(d),\tilde{R}$ stay bounded with increasing $p,\,d,$ and $E$, uniformly over all possible ground truth matrices $B^{\ast}\in\mathcal{B}(p,d,\eta)$. Hence, even though the explicit dependence on $p,\,d,$ and $E$ in the upper bounds (3.11) matches the lower bounds (3.3), we can not claim this rate to be (near) minimax optimal.

Remark 6**.**

Comparing the definitions of $\rho(d)$ and $R(d)$, one might prefer an alternative definition of the former of the form

[TABLE]

In fact, these two quantities are related, albeit for different values of $d$, see [BLT18, Section 8]. We choose $\rho(d)$ instead of $\tilde{\rho}(d)$ for the sake of a simpler presentation.

In order to address the issues raised in the previous remark, we next give a penalized maximum likelihood estimator.

3.3 Upper bounds for two-step penalized likelihood

3.3.1 Two-step maximum likelihood estimator

One shortcoming in the rate for LLC for large $n$ in Theorem 4 are the constants $\rho(d)$ and $\tilde{R}$ which might actually grow with $p$, see Remark 5. Moreover, as a moment estimator, it does not naturally behave well with respect to model misspecification. This motivates a different estimator based on a penalized maximum likelihood approach.

Recall that the negative log-likelihood of a multivariate Gaussian with empirical covariance matrix $\hat{\Sigma}$ and precision matrix $\Theta$ is given by.

[TABLE]

Thus, the negative log-likelihood for the whole model is proportional to

[TABLE]

where $\Theta^{e}(B)=(I-U_{e}B)^{\top}(I-U_{e}B),$ and

[TABLE]

In order to exploit sparsity in the underlying matrix $B^{\ast}$, we need to penalize $\mathcal{L}(B)$ before minimizing it. However, due to the non-linear dependence of $\Sigma^{e}$ on $B$, a vanilla $\ell_{1}$-penalization term might be acting at the wrong scale globally. To overcome this limitation, we propose a two-step estimation procedure. First an initial guess $\hat{B}_{\mathrm{init}}$ is produced using a penalization acting on the scale of the concentration matrices. This initial guess is subsequently refined to $\hat{B}$ as the solution to the $\ell_{1}$-penalized log-likelihood restricted to a small ball around $\hat{B}_{\mathrm{init}}$.

In the first step, we employ penalization with a term resembling a graphical lasso penalty for each experiment,

[TABLE]

leading to the penalized log-likelihood

[TABLE]

The initialization estimator is then given by

[TABLE]

Note that this is not a convex optimization problem and it is hard to solve in general. However, we do give a local optimization algorithm in Section 4 that attempts to find a local minimum for (3.14).

In the second step, this estimator is refined by employing a different regularization term,

[TABLE]

[TABLE]

and the estimator is given by

[TABLE]

with a suitably chosen localization parameter $R_{\mathrm{loc}}>0$.

The loss function (3.16) is again non-convex and hence hard to optimize, but local optimization algorithms seem to produce good results, see Section 4.

3.3.2 Statistical performance

Assuming we have access to the global minima $\hat{B}_{\mathrm{init}}$ and $\hat{B}_{\mathrm{loc}}$, we show the following rates for $\hat{B}_{\mathrm{loc}}$:

Theorem 7**.**

Under assumptions A1 – A3, if

[TABLE]

and the parameters for the estimators $\hat{B}_{\mathrm{init}}$ and $\hat{B}_{\mathrm{loc}}$ are chosen such that

[TABLE]

then

[TABLE]

with probability at least $1-\delta$.

The proof is deferred to Section E. It is based on the one hand on convexity properties of the Gaussian log-likelihood function that were developed in the context of convex optimization problems for estimation of sparse concentration matrices in [RBLZ08] and [LW13], and on the other to new structural results on the difference $\Theta^{e}(B)-\Theta^{\ast,e}$ between concentration matrices expressed in terms of $B-B^{\ast}$; see Lemma 18.

Note that the upper bound (3.20) is worse by a factor of $E$ and a log factor than the lower bound (3.3) in Theorem 3. However, the completely separating system $\mathcal{E}$ can be chosen to be as small as $E\asymp\log(p)$, see [HEH13] and Remark 2, in which case this eventual rate is almost minimax optimal up to logarithmic terms.

4 Numerical experiments

Recall that $\ell(\Theta,\hat{\Sigma})=\mathop{\mathsf{Tr}}(\hat{\Sigma}\Theta)-\log\det(\Theta)$ and that we want to find solutions to the two regularized maximum likelihood problems,

[TABLE]

Both problems are non-convex and there is no obvious strategy for how to find global minima. However, since they are continuous, we can empirically study the performance of optimization algorithms designed for convex problems, hoping to obtain at least local minima. In the following, we describe how candidate solutions for both (4.1) and (4.2) can be found efficiently and demonstrate their performance based on experiments with synthetic data. Additionally, we give a low-rank update approach in Appendix B that can be used to speed up calculations when the number of experiments $E$ is large, but for each experiment, the number of controlled variables $|\mathcal{J}_{e}|$ is small.

4.1 Solving the initialization problem by non-convex ADMM

The difficulty in solving problem (4.1) is to handle the non-smooth penalty terms of non-linear transformations of $B$, $\|\Theta^{e}(B)\|_{1}$. We use a non-linear version of the Alternating Direction Method of Multipliers (ADMM) algorithm, which allows us to introduce additional variables $\Theta^{e}$, constrain them to fulfill $\Theta^{e}=\Theta^{e}(B)$, and keep the resulting dimensionality blowup manageable.

The ADMM algorithm [GM76, GM75, BPC*+*11] is a splitting algorithm intended to solve convex optimization problems of the form

[TABLE]

where $x\in{\rm I}\kern-1.79993pt{\rm R}^{m},y\in{\rm I}\kern-1.79993pt{\rm R}^{\ell}$, $f$ and $g$ are convex functions on ${\rm I}\kern-1.79993pt{\rm R}^{m}$ and ${\rm I}\kern-1.79993pt{\rm R}^{\ell}$, respectively, and $F\in{\rm I}\kern-1.79993pt{\rm R}^{m\times k}$, $G\in{\rm I}\kern-1.79993pt{\rm R}^{\ell\times k}$, $b\in{\rm I}\kern-1.79993pt{\rm R}^{k}$. Introducing the dual variable $u$, a step size $\rho>0$, and starting with an initialization $x^{0},y^{0},u^{0}$, its iterations are given by

[TABLE]

which is the so called scaled form of ADMM.

Note that while in the case of convex objective functions and linear constraints, there are well-established convergence results for ADMM, [Gab83, EB92], results about convergence to a stationary point for non-convex variants are scarce, requiring either linear constraints [WYZ15] or further modifications and additional assumptions [BKSV15].

In order to apply a non-convex ADMM variant, we rewrite problem (4.1) as

[TABLE]

Then, introducing dual variables $\Lambda^{e}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$, $e=1,\dots,E$, the outer iteration of our algorithm is given by

[TABLE]

Note that (4.6) is a convex problem, resembling the graphical LASSO [FHT08] or SPICE [RBLZ08] but with an additional quadratic penalty term. We can solve these subproblems with an extension of the QUIC algorithm [HDRS11] that employs coordinate descent to iteratively find Newton directions.

Problem (4.7) on the other hand is a non-convex problem, albeit without constraints. Hence, we can use any local optimization algorithm. For our experiments, we choose L-BFGS [LN89] to perform this approximate minimization, yielding a stationary point of the objective function.

In order to find a suitable step size parameter $\rho$, we allow varying $\rho^{k}$ and employ the dual-balancing strategy from [HYW00, WL01].

4.2 Solving local problem by Augmented Lagrangian Method

In order to find a local minimum of (4.2), we employ the Augmented Lagrangian Method [NW06] that transforms the inequality constraint into a box constraint and iteratively solves for the associated dual variable. It leads to the following iteration, where $u$ is a slack variable for the $\ell_{2}$ constraint and $\lambda$ is the associated dual variable.

[TABLE]

To solve (4.9), we use L-BFGS-B [ZBLN97], transforming the $\ell^{1}$-regularization into a linear term with additional non-negativity constraints,

[TABLE]

4.3 Experimental setup

We perform experiments with synthetic and semi-synthetic data to gauge the performance of the maximum likelihood procedure, comparing it to the LLC algorithm [HEH12].

For the synthetic benchmarks, we study two types of graph structures: (directed) random regular graphs and graphs composed of disconnected cliques. For the semi-synthetic benchmarks, we use a gene-regulatory network from [CBG13] consisting of 39 genes. Note that in [CBG13], the authors employ a model very similar to ours, but instead of allowing controlled experiments on certain nodes, they consider so-called expression quantitative trait loci (eQTL) as proxies for interventions, which changes their model compared to the one considered here. Nonetheless, part of the output of their estimator is a linear causal network, which is what we consider as ground truth to simulate data following the Gaussian model introduced in Section 2. Example ground truth matrices for the two random models and the semi-synthetic matrix from [CBG13] are given in Figure 1. There, we set $p=39$ and $d=3$ for the random models to coincide with $p$ and $d$ in the semi-synthetic case.

In the following two sections, we give more details about data generation and parameter tuning.

4.3.1 Models

Synthetic graphs:

The ground truth graphs are generated by first obtaining the (directed) adjacency matrix $B_{\mathrm{adj}}\in\{0,1\}^{p\times p}$, a matrix $B_{\mathrm{val}}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ containing edge values, and finally setting $B^{\ast}$ to be the Hadamard product of the two, normalized to have operator norm $1-\eta=0.5$,

[TABLE]

Here, $B_{\mathrm{val}}$ consists of independent standard Gaussian entries, and $B_{\mathrm{adj}}$ is the adjacency matrix of either a regular random graph or one composed of disconnected cliques.

Random regular graphs:

$\operatorname{supp}((B_{\mathrm{adj}})_{i,:})$ is constructed by sampling $d$ times uniformly at random without replacement from $\{1,\dots,p\}\setminus\{i\}$ and all elements in the support are assigned $1$.

Disconnected cliques:

$B_{\mathrm{adj}}$ is the adjacency matrix of a graph consisting of $\lfloor p/d\rfloor$ disconnected $d$-cliques and an additional disconnected $p-d\lfloor p/d\rfloor$ clique if $d$ does not divide $p$. This model is meant to illustrated clustered variables that operate in modules, akin to the ones arising in gene regulatory networks.

4.3.2 Tuning parameters

Choice of $\lambda$:

To keep the comparison simple, we use an oracle choice of $\lambda_{\mathrm{init}}$, $\lambda_{\mathrm{loc}}$ and $R_{\mathrm{loc}}$. For the first two, this means choosing them such that $\|\hat{B}_{\mathrm{init}}-B^{\ast}\|_{F}$ and $\|\hat{B}_{\mathrm{loc}}-B^{\ast}\|_{F}$ is minimal. For $R_{\mathrm{loc}}$, we choose $R_{\mathrm{loc}}=2\|\hat{B}_{\mathrm{init}}-B^{\ast}\|_{F}$. In practice, both parameters could be chosen by cross-validation.

Initialization of optimization algorithm:

We initialize the calculation of $\hat{B}_{\mathrm{init}}$ for the largest value of $\lambda_{\mathrm{init}}$ with the all zeros matrix and then warmstart the calculation with the output of the calculation for the next larger value of $\lambda_{\mathrm{init}}$. The calculation of $\hat{B}_{\mathrm{loc}}$ is initialized with the output of $\hat{B}_{\mathrm{init}}$. To investigate the dependence on the initialization, we also try initializing the calculation of $\hat{B}_{\mathrm{init}}$ with a strict triangular matrix whole upper elements consist of independent $\mathcal{N}(0,10)$ random variables, as well as running the likelihood optimization without constraints as described previously with the same initialization.

Systems of interventions:

We consider three choices for the experiments $\mathcal{E}$. The first one, which we call binary, consists of separating the nodes with a bisection approach similar to the construction given in [Dic69] that leads to $E=O(\log p)$. The second one, which we call bounded, is given by [Cai84] and produces experiments whose sizes $|\mathcal{J}_{e}|$ are bounded by $k$. In this case, $E=O(n/k)$. The third kind corresponds to $k=1$, taking $\mathcal{J}_{i}=\{i\}$ for $i\in[p]$, which we call single-node experiments.

Repetitions:

All errors are averaged over 32 random repetitions of sampling $B^{\ast}$ and the observations $X^{e}_{k}$.

4.4 Results

4.4.1 Performance

Random regular graphs with oracle choice:

In Figure 2, we collect comparisons for the estimation rates of $\hat{B}_{\mathrm{llc}}$, $\hat{B}_{\mathrm{init}}$, and $\hat{B}_{\mathrm{loc}}$, varying $n,p,d,$ and $E$, respectively, where the varying $E$ case is given by bounded experiments with a varying bound on the size $k$ of the experiments which, of course, governs the total number $E$ of experiments needed for separation. In all other cases, we consider binary experiments.

Figure 2(a) indicates that all three estimators exhibit a risk that scales as $1/n$ and displays a clear ordering in the performance of the three candidates where $\hat{B}_{\mathrm{llc}}$ performs worse than $\hat{B}_{\mathrm{init}}$, which in turn is worse than $\hat{B}_{\mathrm{loc}}$.

In Figure 2(b), we observe a scaling with respect to $p$ that is slightly worse than guaranteed by our theorems and could be due to the presence of log factors. In Figure 2(c), we in turn see that the scaling with respect to $d$ is slightly better than expected, hinting at good adaptation to the sparsity parameter $d$. Most interestingly, in Figure 2(d), we observe that the scaling with respect to $E$ when increasing the number of experiments appears to be better than predicted by our theory: about $E^{1/2}$ for $\hat{B}_{\mathrm{llc}}$ and $\hat{B}_{\mathrm{init}}$, about $E^{1/3}$ for $\hat{B}_{\mathrm{loc}}$. This different behavior is even more striking in Figure 3(b) where the performance of $\hat{B}_{\mathrm{loc}}$ appears to decay at most logarithmically in $E$.

Disconnected clique graphs:

In Figure 3(a), we plot the same experiment as in Figure 2(a), only this time with disconnected clusters instead of random regular graphs. We notice a similar behavior, with the key difference of the performance of $\hat{B}_{\mathrm{init}}$ surpassing that of $\hat{B}_{\mathrm{loc}}$. This could be explained by the fact that the penalization in the objective $\hat{B}_{\mathrm{init}}$ is particularly suited for the estimation of this kind of graphs since the sparsity of $(I-B^{\ast})^{\top}(I-B^{\ast})$ in this case almost coincides with the one of $B^{\ast}$, which can be seen from the argument that led to (E.22) in the proof of Theorem 7.

Semi-synthetic graph:

The performance of the three estimators on the semi-synthetic data built from the graph taken in [CBG13] appears in Figure 3(b). The LLC estimator $\hat{B}_{\mathrm{llc}}$ performs similarly to $\hat{B}_{\mathrm{init}}$ and both suffer in comparison to $\hat{B}_{\mathrm{loc}}$ either in terms of absolute performance and in terms of scaling with $E$.

4.4.2 Stability

Role of initialization:

In Figure 4(a), we show the same setup as in Figure 2(a), only this time, the calculation for $\hat{B}_{\mathrm{init}}$ is initialized with a random matrix as outlined in Section 4.3. Additionally, we plot the result of optimizing an unconstrained version of $\hat{B}_{\mathrm{loc}}$ with the same bad initialization, denoted by $\hat{B}_{\mathrm{unconstr}}$. We observe that the performance of the latter is very bad due to the non-convex nature of the objective together with the fact that a bad initialization point is chosen. However, even though $\hat{B}_{\mathrm{init}}$ is found through solving a non-convex objective as well, it seems to be robust enough to yield comparable performance and hence serve as a good initialization for calculating $\hat{B}_{\mathrm{loc}}$ even with a poor initial choice of $B$.

Missing experiments:

In Figure 4(b) we investigate the robustness to systems of interventions that do not fulfill the separability condition in Definition 1. For this, we consider single-node experiments and plot the number of experiments that are missing from a completely separating set of such experiments (in which case we would have $E=p$). The likelihood based approaches are much more robust in this case, and to a larger degree than degree of freedom calculations as in Appendix A would suggest.

Appendix A Non-identifiability in the cyclic case for equal variances

In this section, we give a brief argument to show that for generic matrices $B^{\ast}$, unlike the acyclic case considered in [LB14, PB14], having equal noise variance as required in Assumption A3 does not lead to identifiability from observational data.

The argument is based on counting dimensions of the null space of the non-linear maps

[TABLE]

One limitation of our argument is that it does not cover the potential identifiability of $B^{\ast}$ from observational data under the additional assumption of bounded in-degree $d(B^{\ast})\leq d$.

Proposition 8**.**

Define the integer

[TABLE]

*Then the matrix $B^{\ast}\in\mathcal{B}$ is not uniquely determined by $\Theta^{\ast,e}=\Theta^{e}(B^{\ast}),e\in\mathcal{E}$ whenever $m<p^{2}-p$. In particular, without interventions, this condition holds as soon as $p\geq 4$. *

Proof.

Consider the maps

[TABLE]

defined by stacking all $\Theta^{e}$ into one vector and accepting respectively matrices with zero diagonal and arbitrary diagonal. Similarly, denote by $\bar{\Theta}^{e}$ the map $\Theta^{e}$ when not restricted to matrices with zero-diagonal.

We show that the derivative of $\Theta$ has constant rank bounded above by $p^{2}-p-(m-p)$ at a point $B^{\ast}\in\mathcal{B}$. In turn, whenever $m<p$, this implies the existence of $\tilde{B}\neq B^{*}$ such that $\Theta(B^{\ast})=\Theta(\tilde{B})$ by the constant rank theorem

First, let $B\in\mathcal{B}$ be arbitrary and compute the derivative of $\bar{\Theta}$ at a point $B$ by computing the derivative $D\bar{\Theta}^{e}(B)$ of the individual maps $\bar{\Theta}^{e}:{\rm I}\kern-1.79993pt{\rm R}^{p\times p}\to{\rm I}\kern-1.79993pt{\rm R}$. For any $\bar{H}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$, it holds

[TABLE]

where we used the fact that $U_{e}^{2}=U_{e}$.

Next, we compute the dimension of the null space of $D\Theta(B)$. To that end, observe that for any $H\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ such that $D\Theta^{e}(B)[H]=0$, it holds $D\bar{\Theta}^{e}(B)[\bar{H}]=0$ and $(\bar{H})_{ii}=0$ for all $i\in[p]$. To characterize the dimensionality of the subspace of such matrices $H$, we first consider the null space of $D\bar{\Theta}^{e}(B)[\bar{H}]=0$ and then intersect it with the subspace given by $(\bar{H})_{ii}=0$.

Abbreviate $\bar{G}=\bar{H}A$. By (A.6), $D\bar{\Theta}^{e}(B)[\bar{H}]=0$ for all $B\in\mathcal{B}$ whenever $(U_{e}\bar{G})+(U_{e}\bar{G})^{\top}=0$ for all $e$. We permute the indices such that $\mathcal{J}_{e}=\{1,\dots,|\mathcal{J}_{e}|\}$ to write this equality in block form:

[TABLE]

For each $e\in\mathcal{E}$ the three nonzero blocks above translate into the following conditions:

[TABLE]

As a result $\bar{H}=\bar{G}(I-B)$ is the image of $(I-B)$ through the linear operator $\bar{G}$ that lives in the intersections of the orthogonal subspaces defined by the above constraints. Thus, each constraint contributes 1 to the codimension of the null space of $D\bar{\Theta}(B)$. Equivalently, each constraint contributes 1 to the rank $D\bar{\Theta}(B)$.

Next, we discuss how to deal with the fact that we need to compute $\mathrm{rank}(D\Theta(B))$ instead of $\mathrm{rank}(D\bar{\Theta}(B))$, where $B$ is restricted to lie in the subspace of matrices with zero-diagonal, thus we need to restrict $\bar{H}$ above accordingly. Intuitively, we want to say that the rank can increase by at most $p$, the number of additional linear constraints on the null space, but we need to further establish that there is a $B^{\ast}\in\mathcal{B}$ such that the rank of $D\Theta(B)$ is constant in a neighborhood of $B^{\ast}$. Adding to that the constraint that $\bar{H}$ has null diagonal, we get $\mathrm{rank}(D\Theta(B))\leq m$, where $m$ is defined in (A.2).

Next, we show that $\mathrm{rank}(D\Theta(B))$ is, in fact, constant and equal to some $r^{*}$ in a neighborhood of $B^{*}$ to apply the constant rank theorem.

To that end, let $B^{\ast}$ be such that $r^{\ast}:=\mathrm{rank}(D\Theta(B^{*}))\geq\mathrm{rank}(D\Theta(B))$ for all $B\in\mathcal{B}$. Considering $D\Theta(B^{\ast})$ a matrix let $S^{\ast}$ be a maximal principal minor and denote the restriction of $D\Theta(B^{\ast})$ to $S^{\ast}$ by $[D\Theta(B^{\ast})]_{S^{\ast}}$. By definition, we have

[TABLE]

Moreover, the map $B\mapsto f(B):=\det[D\Theta(B^{\ast})]_{S^{\ast}}$ is a polynomial in the elements of $B$ such that $f(B^{\ast})\neq 0$. By continuity, it also holds that $f(B)\neq 0$ in an open neighborhood of $B^{\ast}$ as well, and thus $\operatorname{rank}(D\Theta(B))\geq r^{\ast}$ in that neighborhood. But since $r^{\ast}$ is maximal, $\operatorname{rank}D\Theta(B)=r^{\ast}$ for $B$ in an open neighborhood of $B^{\ast}$.

The above means that we can apply the constant rank theorem [Boo86, Theorem II.7.1] to obtain diffeomorphisms $\varphi\colon{\rm I}\kern-1.79993pt{\rm R}^{p^{2}-p}\supset V_{1}\to U_{1}\subseteq{\rm I}\kern-1.79993pt{\rm R}^{p^{2}-p}$, $\psi\colon{\rm I}\kern-1.79993pt{\rm R}^{E\times p^{2}}\supset U_{2}\to V_{2}\subseteq{\rm I}\kern-1.79993pt{\rm R}^{E\times p^{2}}$, with $U_{j},V_{j}$ open sets for $j\in\{1,2\}$ such that

[TABLE]

If $r^{\ast}<p^{2}-p$, we obtain a continuum of pre-images of $\Theta(B^{\ast})$ as

[TABLE]

for all $(0,\dots,0,x_{p^{2}-p-r^{\ast}+1},\dots,x_{p^{2}-p})\in V_{1}$, which includes points other than [math] because $V_{1}$ is an open set.

To conclude, recall that $r^{*}\leq m$ so that $m<p^{2}-p$ is a sufficient condition for the failure of injectivity of $\Theta$. This completes the first part of the proof.

To obtain the conclusion without interventions, note that in this case

[TABLE]

so that $m<p^{2}-p$ whenever $p>3$. ∎

Appendix B Numerical speed-up

When many experiments are performed with a small number of nodes that are intervened on, say $|\mathcal{J}_{e}|\leq k$, calculating the log-likelihood term in the algorithms considered in Section 4 in a naive way takes $O(Ep^{3})$ operations: both calculating $\mathop{\mathsf{Tr}}(\Theta^{e}(B)\hat{\Sigma}^{e})$ and performing a Cholesky decomposition for each of the $E$ matrices $\Theta^{e}(B)=L^{e}(L^{e})^{\top}$ with $L^{e}$ lower triangular takes $O(p^{3})$ time. The Cholesky decomposition in turn is used to compute

[TABLE]

The computational complexity can be improved by using a low rank decomposition of $\Theta^{e}(B)$, both for computing the trace term and the Cholesky decomposition $(I-B)^{\top}(I-B)=LL^{\top}$. To see this, write $J_{e}=I-U_{e}$ and decompose

[TABLE]

Hence, the Cholesky decomposition $L^{e}$ can be computed by a rank $k$ update followed by a rank $k$ downdate of $L$, which takes $O(kp^{2})$ [See04]. Computation of the trace terms $\mathop{\mathsf{Tr}}(\hat{\Sigma}^{e}\Theta^{e}(B))$ can be sped up analogously, also taking $O(kp^{2})$ time.

Hence, the total time to compute the log-likelihood is $O(p^{3}+Ekp^{2})$. In a similar manner, computing the objective function for step (4.7) in the non-convex ADMM procedure can be done in $O(p^{3}+Ekp^{2})$ time, although one iteration takes $O(Ep^{3})$ time due to the complexity of performing step (4.6).

Appendix C Proof of lower bounds

C.1 Proof of Theorem 3

To begin, recall the definition of the redundancy factor

[TABLE]

The proof of Theorem 3 is based on standard techniques for minimax lower bounds [Tsy09].

Theorem 9** ([Tsy09, Theorem 2.5]).**

Denote by $\mathcal{G}\subseteq{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ a set of possible hypotheses with associated probablity measures $P_{B}$ for $B\in\mathcal{G}$.

Fix $M\geq 2,s>0,\alpha\in(0,1/8),$ and assume that there exists $B_{0},\dots,B_{M}\in\mathcal{G}$, such that

1. (i)

$\|B_{j}-B_{k}\|_{F}\geq 2s>0$* for all $0\leq j<k\leq M$;* 2. (ii)

$\mathsf{KL}(P_{j}|P_{0})\leq\alpha\log M$* for all $j=1,\dots,M$, where $P_{j}=P_{B_{j}}$.*

Then,

[TABLE]

where the infimum in (C.2) is taken over all measurable functions $\hat{B}$ on the observations.

Before proceeding with the proof, we first present two lemmas. Lemma 10 gives a way to upper bound the Kullback-Leibler divergence between two Gaussian distributions in terms of their concentration matrices, while Lemma 11 contains a version of the Varshamov-Gilbert lemma adopted to produce candidate matrices with the same row sparsity. Their proofs can be found in Section C.2 and Section C.3, respectively.

In the following, we denote by $d_{H}(A,B)$ the Hamming distance between two matrices $A,B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$. It is defined by $d_{H}(A,B)=|\{(i,j)\in[p]^{2}:A_{i,j}\neq B_{i,j}\}|$.

Lemma 10**.**

Let $\Theta_{1},\Theta_{2}\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ be two positive definite concentration matrices and $P_{1}=\mathcal{N}(0,\Theta_{1}^{-1})$ and $P_{2}=\mathcal{N}(0,\Theta_{2}^{-1})$ the associated Gaussian distributions. If

[TABLE]

then,

[TABLE]

Lemma 11**.**

Given $m\geq 1$ and $1\leq d\leq m/2$, there is a family $H_{1},\dots,H_{M}$ of matrices in $\{0,1\}^{m\times m}$ such that

1. (i)

Every row of $H_{i}$ is $d$-sparse for $i=1,\dots,M$; 2. (ii)

$\displaystyle d_{H}(H_{i},H_{j})\geq\frac{md}{2},\quad i\neq j$; 3. (iii)

$\displaystyle\log M\geq\frac{md}{16}\log\left(1+\frac{m}{2d}\right).$**

Taking the above lemmas as given, we proceed to prove Theorem 3.

We apply Theorem 9 by constructing an appropriate set of hypotheses $B_{0},\dots,B_{M}$. Without loss of generality, assume that $p$ is even. Set $B_{0}=0$ to be the all zeros matrix, and apply Lemma 11 with $m=p/2$ to obtain $M$ matrices $H_{1},\dots,H_{M}\in\{0,1\}^{p/2\times p/2}$ with pairwise Hamming distance at least $pd/4$ and with

[TABLE]

We define $B_{i}$, $i=1,\dots,M$ as block matrices by setting

[TABLE]

By construction, for every $i\in[M]$, every row of $B_{i}$ is $d$-sparse, and $B_{i}$ has zero-diagonal. Moreover, by $\kappa\geq 1$, $\alpha<1/8$, and by assumption

[TABLE]

so we get from the Gershgorin circle theorem that

[TABLE]

in light of $\eta\leq 1/2$. Hence, $B_{i}\in\mathcal{B}(p,d,\eta)$ for all $i\in[M]$.

Next, we can lower bound the pairwise distances by

[TABLE]

which yields the needed separation in Theorem 9(i).

We proceed to estimate the KL divergence between two distributions corresponding to matrices $B_{0}$ and any $B_{i}$, $i\geq 1$. Decompose the difference between the concentration matrices as

[TABLE]

Because $\|B_{i}\|_{\mathrm{op}}\leq 1/5$, (C.11) together with $\|U_{e}\|_{\mathrm{op}}\leq 1$ and the sub-multiplicativity of the operator norm implies

[TABLE]

so the hypothesis of Lemma 10 is satisfied for all pairs $(\Theta^{e}(B_{i}),\Theta^{e}(B_{0}))$. By Lemma 10, (C.12), and the tensorization property of the KL divergence, we obtain

[TABLE]

Since $B_{i}$ is defined to be anti-symmetric, we have

[TABLE]

Moreover,

[TABLE]

Combining (C.15), (C.16) and (C.17), we arrive at

[TABLE]

It remains to compute the Frobenius norm of each $B_{i}$,

[TABLE]

Hence, because

[TABLE]

we obtain

[TABLE]

Finally, we can pick $\alpha=\frac{1}{16}$ in Theorem 9 to conclude that

[TABLE]

for some constant $c>0$.

C.2 Proof of Lemma 10

The Kullback-Leibler divergence between two Gaussians $P_{1}=\mathcal{N}(0,\Theta_{1}^{-1})$ and $P_{2}=\mathcal{N}(0,\Theta_{2}^{-1})$ is given by

[TABLE]

Using the fact that the first derivative of $\Theta\mapsto-\log\det\Theta$ is $-\Theta^{-1}$ and the second derivative is $\Theta^{-1}\otimes\Theta^{-1}$, we employ a Taylor expansion about $\Theta_{1}$ (compare (E.31) in proof of Lemma 17) to obtain

[TABLE]

for some $\tilde{\Theta}=t\Theta_{1}+(1-t)\Theta_{2}$, $t\in[0,1]$. By considering the square root of $\tilde{\Theta}^{-1}$, this can be expressed in terms of the Frobenius norm of the difference $\Theta_{2}-\Theta_{1}$,

[TABLE]

By Weyl’s inequality, [Fra12, Section 6.7, Theorem 2],

[TABLE]

Hence, if $\|\Theta_{2}-\Theta_{1}\|_{\mathrm{op}}\leq\lambda_{\mathrm{min}}(\Theta_{2})/2$, then

[TABLE]

which together with (C.28) implies (C.4).

C.3 Proof of Lemma 11

We use the probabilistic method to show the existence of the family $H_{1},\dots,H_{M}$, modifying a standard argument that can be found in [Tsy09, Lemma 2.9].

Let $H_{1},\dots,H_{M}$ be $M$ independent random matrices $H_{k}$, where each row of $H_{k}$ is a zero-one-vector corresponding to a subset of $[p]$ with cardinality $d$ drawn uniformly at random. More precisely, for the $i$th row of the matrix $H_{k}$, draw $U_{1}^{i}$ uniformly from $\{1,\dots,m\}$, and $U_{j}^{i}$ conditioned on $U^{i}_{1},\dots,U^{i}_{j-1}$ uniformly from the set $\{1,\dots,m\}\setminus\{U^{i}_{1},\dots,U^{i}_{i-1}\}$, $j=2,\dots,d$. Then, set

[TABLE]

By a union bound, the probability that there exists a pair $k,\ell$ for which $d_{H}(H_{k},H_{\ell})\leq md/2$ can be bounded by the probability of this occurring for one draw of $H_{1}$, comparing to a fixed $H_{0}$ with d-sparse rows, say $(H_{0})_{kl}={\rm 1}\kern-2.40005pt{\rm I}_{\{k\leq d\}}$,

[TABLE]

because for each row, every $d$ sparse pattern is equally likely.

We can lower bound the Hamming distance by the number of elements in $\mathrm{supp}(H_{0})$ on which $H_{1}$ is one, i.e.,

[TABLE]

with $Z_{i,j}={\rm 1}\kern-2.40005pt{\rm I}(U^{i}_{j}\leq d)$. Then, $Z_{i,j}\sim\mathsf{Bern}(q_{i,j})$ with $q_{i,1}=\frac{d}{m}$ and, noting that $d\leq m/2$,

[TABLE]

From there, apply a Chernoff bound, that is, pick $\lambda>0$ and estimate

[TABLE]

For a Bernoulli distribution $Z\sim\mathsf{Bern}(q)$, the moment generating function is given by

[TABLE]

Together with the observation that $Z_{i,j}$ is stochastically dominated by a $\mathsf{Bern}(2d/m)$ distribution, we can estimate

[TABLE]

Setting $\lambda=\log\left(1+\frac{m}{2d}\right)$, we get that

[TABLE]

Combining (C.38), (C.35), (C.32), and the estimate ${M\choose 2}\leq M^{2}$, we obtain

[TABLE]

since $d\leq m/2$, provided we choose $M$ such that

[TABLE]

Setting,

[TABLE]

we have thus shown that there exists a family fulfilling the conditions of Lemma 11.

Appendix D Proof of LLC upper bounds

D.1 Notation and lemmas

We start by recalling the notation from Section 3.2. We denote the $i$th row of $B^{\ast}$ by $b^{\ast}_{i}\in{\rm I}\kern-1.79993pt{\rm R}^{p-1}$, omitting the diagonal element which is assumed to be zero. From empirical covariances of the performed experiments, for $i\in[p]$, we obtain estimators $\hat{T}_{i}$ for $T^{\ast}_{i}$ and $\hat{t}_{i}$ for $t^{\ast}_{i}$, where $T^{\ast}_{i}b^{\ast}_{i}=t^{\ast}_{i}$. Then, we solve the associated $\ell^{1}$-regularized least squares problem,

[TABLE]

and assemble its solutions into $\hat{B}_{\mathrm{llc}}$ as

[TABLE]

where $P_{i}\in{\rm I}\kern-1.79993pt{\rm R}^{(p-1)\times p}$ denotes the projection matrix that omits the $i$th coordinate.

In particular, the $T^{\ast}_{i}$ are defined row-wise, adding a row $\mathfrak{e}_{j}^{\top}\Sigma^{\ast,e}P_{i}^{\top}$. for each experiment $e$ such that $i\in\mathcal{U}_{e}$ and each entry $j\in\mathcal{J}_{e}$. Similarly, the vector $t^{\ast}_{i}$ is defined by appending the corresponding entries $\Sigma^{\ast,e}_{j,i}$. Estimators for $T^{a}st_{i}$ and $t^{\ast}_{i}$ in turn are given by row-wise assembling empirical counterparts of the above quantities, so that a generic $\ell$th row of $T^{\ast}_{i}$ and $\ell$th entry of $t^{\ast}_{i}$ are given by

[TABLE]

respectively.

Next, recall the following quantities that enter the rate:

[TABLE]

where

[TABLE]

We restate Theorem 4 for convenience.

Theorem 4.

Let assumptions A1 – A3 hold and fix $\delta\in(0,1)$. Assume further that

[TABLE]

Then, the LLC estimator $\hat{B}_{\mathrm{llc}}$ defined in (3.7) with $\lambda$ chosen such that

[TABLE]

satisfies

[TABLE]

with probability at least $1-\delta$.

The proof relies on the following key lemmas. Lemma 12 yields control on the stochastic error, while Lemma 13 ensures that the linear system we solve via $\ell^{1}$-regularization is well-conditioned for that purpose.

Lemma 12**.**

Under the assumptions of Theorem 4, writing

[TABLE]

for a fixed $C>0$, there is an event $\mathcal{A}$ such that ${\rm I}\kern-1.79993pt{\rm P}(\mathcal{A})\geq 1-\delta$, and on $\mathcal{A}$,

[TABLE]

Lemma 13**.**

Assume that the same hypotheses as in Theorem 4 hold. On the same event $\mathcal{A}$ as in Lemma 12, we have

[TABLE]

where $\mathcal{C}(d)$ is the set of vectors fulfilling the cone condition in (D.7).

D.2 Proof of Theorem 4

Since the experiments are completely separating, it follows from [HEH12] that

[TABLE]

Fix $i\in[p]$ and abbreviate $T^{\ast}=T^{\ast}_{i}$, $\hat{T}=\hat{T}_{i}$, $t^{\ast}=t^{\ast}_{i}$, $\hat{t}=\hat{t}_{i}$, $b^{\ast}=b^{\ast}_{i}$, and $\hat{b}=\hat{b}_{i}$.

On the event $\mathcal{A}$ from Lemma 12 the following holds. By definition of $\hat{b}$, we have

[TABLE]

Set $h=\hat{b}-b^{\ast}$ and rearrange to obtain

[TABLE]

By Hölder’s inequality

[TABLE]

By the assumptions on $n$ and Lemma 12, $\|\hat{T}^{\top}(\hat{T}b^{\ast}-\hat{t})\|_{\infty}\leq 4\tilde{R}\phi_{n}$. Denote by $S$ the support of $b^{\ast}$. By triangle inequality and splitting between $S$ and $S^{c}$, we can bound the regularization term by

[TABLE]

Add $\lambda\|h\|_{1}/2$ on both sides of (D.15) to obtain

[TABLE]

Now, assume $\lambda\geq 8\tilde{R}\phi_{n}$, which by Lemma 12 matches the assumed scaling of

[TABLE]

Together with (D.17), we get that $h$ fulfills the cone condition $\|h_{S^{c}}\|_{1}\leq 3\|h_{S}\|_{1}.$ In turn, by Lemma 13, taking into account that the assumptions on $n$ and $\phi_{n}$ are fulfilled by assumption, we obtain

[TABLE]

Moreover, by the Cauchy-Schwarz inequality $\|h_{S}\|_{1}\leq\sqrt{d}\|h\|_{2}$, so combined with (D.17) and (D.18), we have

[TABLE]

Re-introducing the index $i$ and summing the above over all $i\in[p]$, we get

[TABLE]

D.3 Proof of Lemma 12

The proof of Lemma 12 consists of two parts that correspond to Lemma 14 and Lemma 15 below.

Let $\phi_{n}>0$ and define the events $\mathcal{A}_{1},\mathcal{A}_{2},\mathcal{A}_{3}$ as follows:

[TABLE]

Lemma 14 gives an upper bound on $\|\hat{T}^{\top}(\hat{T}b^{\ast}-\hat{t})\|_{\infty}$ in terms of $\phi_{n}$, while Lemma 15 gives a high-probability bound on $\phi_{n}$. We give the proofs of both of these lemmas after finishing the proof of Lemma 12.

Lemma 14** (Trace term estimate).**

If

[TABLE]

then on the event $\mathcal{A}_{1}\cap\mathcal{A}_{2}\cap\mathcal{A}_{3}$,

[TABLE]

Lemma 15** (Control on stochastic error).**

Let $\delta\in(0,1)$. If $n\gtrsim E\log(e\kappa p/\delta)$,

[TABLE]

and we set

[TABLE]

for a fixed constant $C>0$, we have that

[TABLE]

Adjusting the constants in the requirement on $n$ (D.8) in Theorem 4, we can ensure the requirements (D.24) and (D.26) and thus Lemma 12 follows by setting $\mathcal{A}=\mathcal{A}_{1}\cap\mathcal{A}_{2}\cap\mathcal{A}_{3}$ and combining Lemma 14 and Lemma 15.

Proof of Lemma 14.

We fix $i\in[p]$ and, as before, omit it for notational convenience It holds

[TABLE]

where we used the fact that $T^{\ast}b^{\ast}=t^{\ast}$ and that for an arbitrary matrix $A\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ and vector $x\in{\rm I}\kern-1.79993pt{\rm R}^{p}$, $\|Ax\|_{\infty}\leq\|A\|_{\infty,\infty}\|x\|_{\infty}$. Since by the definition of $\mathcal{A}_{2}$,

[TABLE]

we have that combined with the definitons of $\mathcal{A}_{1}$, $\mathcal{A}_{3}$, and $\tilde{R}$,

[TABLE]

if $\phi_{n}\leq\tilde{R}/p$. ∎

Proof of Lemma 15.

For all three events, we write each element of the associated matrices or vectors as a sum over independent sub-exponential random variables and apply Bernstein’s inequality, Lemma 27.

We start by controlling $\max_{i}\|(\hat{T}_{i}-T^{\ast}_{i})b^{\ast}_{i}\|_{\infty}$. Let $i\in\{1,\dots,p\}$ and $\ell\in\{1,\dots,m_{i}\}$. The $\ell$th row of $\hat{T}_{i}-T^{\ast}_{i}$ corresponds to an experiment $e=e(i,\ell)$ such that $i\in\mathcal{U}_{e}$, and an index $j=j(\ell)\in\mathcal{J}_{e}$, which means we can write

[TABLE]

where we used that $B^{\ast}_{i,i}=0$, so that $P_{i}^{\top}b^{\ast}_{i}=(B^{\ast}_{i,:})^{\top}$. Moreover, with independent normal random vectors $Z^{e}_{k}$ for $e=1,\dots,E$ and $k=1,\dots,n/E$, $\hat{\Sigma}^{e}$ is of the form

[TABLE]

so that

[TABLE]

We proceed to control the $\ell^{2}$ norm of the vectors that are being multiplied with $Z_{k}^{e}$. Lemma 23 yields that

[TABLE]

and

[TABLE]

Hence, by Lemma 26,

[TABLE]

and by Lemma 25, $\mathfrak{e}_{\ell}^{\top}\hat{T}_{i}b^{\ast}_{i}\sim\mathsf{subE}(\eta^{-2})$. Now, Bernstein’s inequality in Lemma 27 and ${\rm I}\kern-1.79993pt{\rm E}[\hat{T}_{i}]=T^{\ast}_{i}$ allows us conclude that for $t_{1}>0$,

[TABLE]

A union bound over all indices $i\in[p]$ and all $\ell$, taking into account that there are at most $\kappa p$ rows in every $T^{\ast}_{i}$, then yields

[TABLE]

Similarly, for $j\in\{1,\dots,p-1\}$ and any row index $\ell$,

[TABLE]

and, as before,

[TABLE]

so that $\mathfrak{e}_{\ell}^{\top}\hat{T}_{i}\mathfrak{e}_{j}\sim\mathsf{subE}(\tilde{\eta}^{-2})$. Hence, by Bernstein’s inequality, for $t_{2}>0$,

[TABLE]

A union bound over all $i\in[p]$, $j\in[p-1]$, and row indices $\ell$ yields

[TABLE]

In particular, the union of the two events in (D.42) and (D.47) occurs with probability at most $\delta$ if

[TABLE]

Taking into account that all $\hat{T}_{i}$ and $\hat{t}_{i}$ are of the form we investigated in (D.44), we get the claim of the lemma if we choose

[TABLE]

for a suitable constant $C$ and assume $n\gtrsim E\log(e\kappa p/\delta)$. ∎

D.4 Proof of Lemma 13

To obtain the result, we employ the following lemma.

Lemma 16** ([LW11, Lemma 12]).**

If for a matrix $\Gamma\in{\rm I}\kern-1.79993pt{\rm R}^{k\times k}$, $k\in\mathbb{N}$ and an integer $s\geq 1$, it holds that

[TABLE]

then

[TABLE]

To this end, let $v\in{\rm I}\kern-1.79993pt{\rm R}^{p-1}$ be a $d$ sparse vector with $\|v\|_{2}=1$, as well as $i\in[p]$, and denote by $\hat{T}=\hat{T}_{i}$, $T^{\ast}=T^{\ast}_{i}$. Then,

[TABLE]

On the one hand, by the definition of $R(d)$, (D.5),

[TABLE]

On the other hand, if the event $\mathcal{A}_{2}$ occurred, then by definition

[TABLE]

Thus, denoting by $S$ the support of $v$, we can further estimate

[TABLE]

Combined, (D.56) and (D.58) yield

[TABLE]

Now, let $h\in{\rm I}\kern-1.79993pt{\rm R}^{p-1}$ be a vector that fulfills the cone condition of order $d$. That is, there is a set of indices $S\subseteq[p-1]$ with $|S|\leq d$ such that $\|h_{S^{c}}\|_{1}\leq 3\|h_{S}\|_{1}.$ This in turn implies that

[TABLE]

by the Cauchy-Schwarz inequality. By Lemma 16, (D.58), and the definition of $R(d)$ in (D.4), we have

[TABLE]

Combined, if

[TABLE]

which is guaranteed from the assumptions of Theorem 12, we get the claim,

[TABLE]

Appendix E Proof of upper bounds for penalized maximum likelihood estimator

E.1 Notation and lemmas

In the following section, we present the proof of Theorem 7, whereas the proofs of several key lemmas are deferred to later sections.

We begin by recalling the estimators and restating Theorem 7. The loss functions are given by

[TABLE]

where

[TABLE]

We consider the penalty terms

[TABLE]

leading to the objective functions

[TABLE]

and

[TABLE]

Finally, the estimators are defined as

[TABLE]

where $\lambda_{\mathrm{init}},\lambda_{\mathrm{loc}}$ and $R_{\mathrm{loc}}$ are tuning parameters that are to be determined.

Theorem 7.

Under assumptions A1 – A3, if

[TABLE]

and the parameters for the estimators $B_{\mathrm{init}}$ and $B_{\mathrm{loc}}$ are chosen such that

[TABLE]

then

[TABLE]

with probability at least $1-\delta$.

First, we present three key lemmas used in the proof of Theorem 7. Lemma 17 yields curvature estimates of the likelihood function in terms of the difference of the concentration matrices associated with a candidate matrix $B$ while Lemma 18 allows us to relate the difference of the concentration matrices to the difference in the underlying matrices, $B-B^{\ast}$. Finally, Lemma 19 gives bounds on a stochastic error term.

To facilitate the presentation, we present the lemmas with the following set of notations and assumptions. Let $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ be an arbitrary matrix and $\mathcal{E}$ a set of completely separating experiments as in assumption A2 with associated matrices $\{J_{e},U_{e}\}_{e\in\mathcal{E}}$. Moreover, assume that $B^{\ast}\in\mathcal{B}(p,d,\eta)$. Then, we denote by

[TABLE]

the concentration matrices associated with $B$ and $B^{\ast}$, respectively, as well as the associated differences between the structure matrices and the concentration matrices by

[TABLE]

respectively. We also abbreviate

[TABLE]

The first lemma follows from convexity arguments that also appear in [RBLZ08, LW13].

Lemma 17** (Lower bounds on Gaussian log-likelihood function, [RBLZ08, LW13]).**

With $\mathcal{L}$ defined as in (3.13), it holds for any $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ that

[TABLE]

*where $c_{1}=18^{-1}$. *

Lemma 18** (Upper and lower bounds on $\|\Delta\|_{F}$ in terms of $\|H\|_{F}$).**

If $B\in\mathcal{B}$, that is, $B$ has zero diagonal, we have

[TABLE]

Lemma 19** (Trace term estimates).**

Let $\delta\in(0,1)$. Denote by $\phi_{n}$ and $\psi_{n}$ the rates

[TABLE]

for an appropriately chosen constant $C>0$. If $n\gtrsim E\log(epE/\delta)$, then with probability at least $1-\delta$, it holds for any $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ that

[TABLE]

and

[TABLE]

For the proof of Theorem 7, we additionally introduce the following abbreviations. For $\diamond\in\{\mathrm{init},\mathrm{loc}\}$, let

[TABLE]

With this, we are ready to give the proof of 7.

E.2 Proof of Theorem 7

Proof sketch: The proof of Theorem 7 is split into two parts. First, we show that the initialization estimator $\hat{B}_{\mathrm{init}}$ performs well enough to allow us to choose $R_{\mathrm{loc}}$ sufficiently small, so that the log-likelihood in an $R_{\mathrm{loc}}$-neighborhood of $B_{\mathrm{init}}$ has large enough curvature. Second, we show that locally, $\hat{B}_{\mathrm{loc}}$ achieves the desired rate.

Both proofs are based on re-arranging the optimality condition for the penalized log-likelihood, bounding the occurring trace term with high-probability, and exploiting the curvature of the log-likelihood function.

Step 1, basic inequality: By definition of the estimator $\hat{B}_{\mathrm{init}}$,

[TABLE]

Comparing to the ground truth $B^{\ast}$ yields the basic inequality

[TABLE]

which implies

[TABLE]

Applying the lower bound on the negative log-likelihood (E.12) in Lemma 17 then yields

[TABLE]

Step 1, estimate error term: Next, we bound the trace term

[TABLE]

with high probability using Lemma 19. For the remainder of the proof, we place ourselves on the event of probability at least $1-\delta$ on which the statement of Lemma 19 holds. Thus, we can estimate the trace term in (E.20) by

[TABLE]

Denoting the support of $\Theta^{\ast,e}=(I-U_{e}B^{\ast})^{\top}(I-U_{e}B^{\ast})$ by $S^{e}_{\mathrm{init}}$, we have

[TABLE]

Moreover, by triangle inequality,

[TABLE]

Combined with the definition of the penalization term,

[TABLE]

Now, assume $\lambda_{\mathrm{init}}\geq\psi_{n}$, which matches the assumed scaling of $\lambda_{\mathrm{init}}$ to obtain

[TABLE]

Note that we can control the size of the support $|S^{e}_{\mathrm{init}}|$ by the in-degree of $B^{\ast}$. Namely, if we decompose

[TABLE]

which is a sum over the outer product of $d+1$ sparse vectors by the assumption that the in-degree of the underlying graph is bounded by $d$, and hence

[TABLE]

In turn, Hölder’s inequality yields

[TABLE]

Bounds on $\|\Delta_{\mathrm{init}}\|_{F}^{2}$: If $\|\Delta_{\mathrm{init}}\|_{F}\geq 1$, by (E.20) and (E.22), we have

[TABLE]

which yields a contradiction if

[TABLE]

By the assumption that $\lambda_{\mathrm{init}}\asymp\psi_{n}$ and the value of $\psi_{n}$ in 19, this holds if

[TABLE]

If $\|\Delta_{\mathrm{init}}\|_{F}\leq 1$, again by combining (E.20) and (E.22), we have

[TABLE]

Dividing by $\|\Delta_{\mathrm{init}}\|_{F}$ and squaring then implies

[TABLE]

By Lemma 19 and the choice of $\lambda_{\mathrm{init}}\asymp\psi_{n}$, this leads to

[TABLE]

Bounds on $\|H_{\mathrm{init}}\|_{F}$: In order to relate $\|\Delta_{\mathrm{init}}\|_{F}$ to $\|H_{\mathrm{init}}\|_{F}$ we appeal to Lemma 18. If $n$ is large enough for (E.23) to hold, then by the lower bound (E.13) in Lemma 18,

[TABLE]

and hence

[TABLE]

which concludes the analysis for the initialization estimator.

Step 2, basic inequality: We have

[TABLE]

Suppose $R_{\mathrm{loc}}\geq\|H_{\mathrm{init}}\|_{F}$, which we achieve by (E.25) and choosing $n$ large enough later, once $R_{\mathrm{loc}}$ has been chosen. Then, comparing to the ground truth $B^{\ast}$ yields the basic inequality

[TABLE]

which implies

[TABLE]

Applying the lower bound on the negative log-likelihood (E.12) in Lemma 17 yields

[TABLE]

Step 2, estimate error term: We resort to Lemma 19, this time in the form of (E.52), which yields

[TABLE]

First, we want to ensure $\|\Delta_{\mathrm{loc}}\|_{F}\leq 1$. The upper bound on $\Delta$ in Lemma 18, (E.15), achieves this if

[TABLE]

for a small enough constant $c_{2}\leq 1$. By the triangle inequality and (E.25), this is true if

[TABLE]

Second, since we want the bound (E.14) to be effective within the ball $\|B-\hat{B}_{\mathrm{init}}\|_{F}\leq R_{\mathrm{loc}}$ over which the optimization in step 2 is constrained, we choose $n$ large enough to guarantee

[TABLE]

This again follows from triangle inequality and (E.25) if

[TABLE]

Third, to control the $\|H_{\mathrm{loc}}\|_{\infty,\infty}$ term in (E.28), observe that

[TABLE]

by Hölder inequality. To guarantee $\|H_{\mathrm{loc}}\|_{\infty,\infty}\leq 2$, it is enough to ask for $\|H_{\mathrm{init}}\|_{F}^{4}\leq 1/p^{2}$ and $R_{\mathrm{loc}}\leq 1/\sqrt{p}$ by triangle inequality. By (E.25), the former is be satisfied if

[TABLE]

Combined, in addition to the assumptions made in step 1, if

[TABLE]

then

[TABLE]

In turn, from (E.14), we obtain

[TABLE]

Writing $S:=\operatorname{supp}(I-B^{\ast})$, we then see that

[TABLE]

and by triangle inequality,

[TABLE]

Together with (E.28) and observing that we can assume $\lambda_{\mathrm{loc}}\geq 4\phi_{n}$, it follows that

[TABLE]

Applying the Cauchy-Schwarz inequality gives

[TABLE]

Finally, we divide by $\|H_{\mathrm{loc}}\|_{F}$, take squares, observe that $|S|\leq p(d+1)$ use $\lambda_{\mathrm{loc}}\asymp\phi_{n}$, and plug in the value of $\phi_{n}$ in Lemma 19 to obtain

[TABLE]

which concludes the proof.

E.3 Proof of Lemma 17

Let $R_{1}>0$ and recall the notation

[TABLE]

for the negative log-likelihood of a centered multivariate Gaussian distribution. Let $\Theta^{\ast},\hat{\Sigma}$ be a positive definite matrix and a positive semi-definite matrix, respectively, and set $\Sigma^{\ast}=(\Theta^{\ast})^{-1}$. Noting that the first derivative of $\Theta\mapsto-\log\det\Theta$ is $-\Theta^{-1}$ and the second derivative is $\Theta^{-1}\otimes\Theta^{-1}$, by computing a Taylor expansion of $\ell$ with differential remainder term about $\Theta^{\ast}$, we have that

[TABLE]

for some $t\in[0,1]$ and $\tilde{\Theta}=\Theta^{\ast}+t(\Theta-\Theta^{\ast})$.

Denote the matrix square root of $\tilde{\Theta}^{-1}$ by $\tilde{\Theta}^{-1/2}$. Then, we can further lower bound the quadratic term by

[TABLE]

By the spectral theorem, we can express the smallest eigenvalue of $\tilde{\Theta}^{-1/2}$ in terms of the largest eigenvalue of $\tilde{\Theta}$,

[TABLE]

Now, recall

[TABLE]

where

[TABLE]

and introduce

[TABLE]

and denote by

[TABLE]

the Frobenius norm of the collection of $\tilde{\Delta}^{e}$ when viewed as a tensor. We now apply the expansion (E.31) and the estimate (E.32) to each of the summands, distinguishing two cases.

First, if $\|\Delta\|_{F}\leq R_{1}$, then also $\|\Theta^{e}-\Theta^{\ast,e}\|_{F}\leq R_{1}$ for all $e\in\mathcal{E}$ and we get

[TABLE]

Therefore, from (E.32) we get a lower bound of the form

[TABLE]

with $c_{1}=(\max_{e\in\mathcal{E}}\|\Theta^{\ast,e}\|_{\mathrm{op}}+R_{1})^{-2}/2$.

Second, if $\|\Delta\|_{F}>R_{1}$, we can leverage the convexity of $\Theta\mapsto-\log\det\Theta$ to again obtain lower bounds. Define $g(s)$ for $s\in[0,1]$ by

[TABLE]

Since $\ell$ is convex in $\Theta$, $g$ is convex in $s$, and we obtain

[TABLE]

Plugging in $t=R_{1}/\|\Delta\|_{F}$, we are in the first case that was discussed and can appeal to (E.33), which yields

[TABLE]

Combined, we get

[TABLE]

Finally, setting $R_{1}=1$ and observing that $\max_{e}\|\Theta^{\ast,e}\|_{\mathrm{op}}\leq 2$ by Lemma 23 yields the claim.

E.4 Proof of Lemma 18

In this section, we abbreviate

[TABLE]

We also need the following linear transformation of $H$, which we denote by $G$,

[TABLE]

First, we give a lemma that allows us to estimate the Frobenius norm of $G$ by its off-diagonal elements.

Lemma 20**.**

Let $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ and denote by $H$ and $G$ the matrices in (E.34) and (E.35), respectively. Moreover, write $G_{D}$ and $G_{D^{c}}$ for the restriction of the matrix $G$ to its diagonal indices and off-diagonal elements, respectively. If $\|B^{\ast}\|_{\mathrm{op}}<1$ and $H_{D}=0$, then

[TABLE]

Proof.

By the definition of $G$, we know that $H=GA^{-1}=G(I-B^{\ast})$. The restriction $H_{D}=0$ implies

[TABLE]

Since $B^{\ast}$ has zero diagonal, for each $i\in[p]$, we can solve for $G_{ii}$ and obtain

[TABLE]

By the Cauchy-Schwarz inequality,

[TABLE]

Finally, summing over all $i$ gives

[TABLE]

Since $\|B^{\ast}\|_{\mathrm{op}}<1$ and by Lemma 23,

[TABLE]

we have the claim,

[TABLE]

With this, we proceed to prove Lemma 18.

To start, let $e\in\mathcal{E}$. We have

[TABLE]

Since $U_{e}^{\top}U_{e}=U_{e}^{2}=U_{e}$, we can simplify the terms in the above expression as

[TABLE]

which leads to

[TABLE]

Hence, by Lemma 23,

[TABLE]

Write $G=HA$.

First, to further lower bound the above expression, consider the diagonal of the $\mathcal{J}_{e}\times\mathcal{J}_{e}$ block of the matrix. There, we have $(U_{e}G)_{\mathcal{J}_{e},\mathcal{J}_{e}}=0$ and $(U_{e}G)^{\top}_{\mathcal{J}_{e},\mathcal{J}_{e}}=0$, and thus

[TABLE]

where we used $\|h\|_{1}\leq\sqrt{p}\|h\|_{2}$ for a vector $h\in{\rm I}\kern-1.79993pt{\rm R}^{p}$, which follows from Hölder’s inquality. Summing over the experiments $\mathcal{E}$, together with the assumption of $\mathcal{E}$ being completely separating, Hölder’s inequality, Lemma 20, and Lemma 23, we get

[TABLE]

Second, focusing on the $\mathcal{U}_{e}\times\mathcal{J}_{e}$ block of the matrix

[TABLE]

we note that by the Cauchy-Schwarz inequality and the elementary inequality $(a+b)^{2}\geq\frac{1}{2}a^{2}-b^{2}$ for $a,b\in{\rm I}\kern-1.79993pt{\rm R}$,

[TABLE]

Summing over the experiments, taking into account that by symmetry the same estimate holds for the $\mathcal{J}_{e}\times\mathcal{U}_{e}$ block, and bounding maximum and minimum singular values by Lemma 23, we obtain a lower bound of

[TABLE]

Finally, we can upper bound $\|\Delta\|_{F}$ in terms of $\|H\|_{F}$, starting from (E.37), by

[TABLE]

E.5 Proof of Lemma 19

In this section, we abbreviate

[TABLE]

and introduce the events

[TABLE]

where the terms $T_{1},T_{2},T_{3}$ are upper bounded by rates $\phi_{n}$ and $\psi_{n}$ to be made precise in Lemma 22, while Lemma 21 shows how $\phi_{n}$ and $\psi_{n}$ can be used to estimate the trace term.

Lemma 21** (Trace term estimates).**

On the event $\mathcal{A}_{1}\cap\mathcal{A}_{2}\cap\mathcal{A}_{3}$, it holds for any $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ that

[TABLE]

and

[TABLE]

Lemma 22** (Control on stochastic error).**

Let $\delta\in(0,1)$. There exists an absolute constant $C$ such that if

[TABLE]

and

[TABLE]

then with probability at least $1-\delta$, it holds that

[TABLE]

where $\mathcal{A}_{1},\mathcal{A}_{2},\mathcal{A}_{3}$ are defined as in (E.50).

Combined, Lemmas 21 and 22 yield Lemma 19.

Proof of Lemma 21.

First, by Hölder’s inequality,

[TABLE]

Identifying the first term as $T_{1}$ and using the estimate $T_{1}\leq\psi_{n}$ yields (E.51).

Second, by the same calculation that led to (E.42), we decompose the trace term as

[TABLE]

The first term in (E.55) can be bounded by

[TABLE]

while the second term can be controlled by

[TABLE]

For each entry of the matrix on the right of (E.57), indexed by $i,j\in[p]$, we have

[TABLE]

so that

[TABLE]

Combined with the estimate (E.55), this yields the second claim, (E.52). ∎

Proof of Lemma 22.

To begin, recall the definition of $\hat{\Sigma}^{e}$, as a sum of $n/E$ *i.i.d. *samples, that is, for $i,j\in[p]$,

[TABLE]

By the definition of the sample distribution, we can write

[TABLE]

where the $Z^{e}_{\ell}$ follow a $\mathcal{N}(0,1)$ distribution and are i.i.d., and Lemma 26 ensures that both $Y_{1}$ and $Y_{2}$ are $\mathsf{subG}(\sigma_{\max}(I-U_{e}B^{\ast}))$ random variables. By Lemma 25, we obtain that $(X^{e}_{\ell})_{i}(X^{e}_{\ell})_{j}\sim\mathsf{subE}(\sigma_{\max}(I-U_{e}B^{\ast})^{2})$.

Similarly,

[TABLE]

Here, we have $\tilde{Y_{1}}\sim\mathsf{subG}(1)$ and $\tilde{Y_{2}}=\mathfrak{e}_{j}^{\top}(I-U_{e}B^{\ast})Z^{e}_{\ell}\sim\mathsf{subG}(\sigma_{\max}(I-U_{e}B^{\ast}))$. By again applying Lemma 25, this means that $(A_{e}^{-1}X^{e}_{\ell})_{i}X^{e}_{\ell,j}\sim\mathsf{subE}(\sigma_{\max}(I-U_{e}B^{\ast}))$.

Having established this, to obtain an estimate for $T_{1}$, we employ Bernstein’s inequality, Lemma 27 to the sum in (E.58) for each $e\in\mathcal{E},i,j\in[p]$ to see

[TABLE]

for $t_{1}>0$, with an absolute constant $c_{B}$ and $K_{1}=\max_{e}\sigma_{\max}(I-U_{e}B)^{2}$. Here, we made use of the fact that subtracting $\Sigma^{\ast,e}$ centers the variables in the sum and that there are $n/E$ independent summands in (E.42). By a union bound,

[TABLE]

To bound $T_{2}$, for $i,j\in[p]$, we write

[TABLE]

with

[TABLE]

By Bernstein’s inequality, Lemma 27, for $t_{2}>0$,

[TABLE]

where $c_{B}$ is an absolute constant and

[TABLE]

A union bound then yields

[TABLE]

To bound $T_{3}$, we proceed similarly. Using $(X^{e}_{\ell})_{i}(X^{e}_{\ell})_{j}\sim\mathsf{subE}(\sigma_{\max}(I-U_{e}B)^{2})$ instead of $(A_{e}^{-1}X^{e}_{\ell})_{i}X^{e}_{\ell,j}\sim\mathsf{subE}(\sigma_{\max}(I-U_{e}B))$, for $t_{3}>0$, we have

[TABLE]

where $K_{3}=\max_{e}\sigma_{\max}(I-U_{e}B^{\ast})^{2}$.

Combined, recalling that by Lemma 23, $\sigma_{\max}(I-U_{e}B^{\ast})\leq 2$ and applying a union bound, we see that the union of the events in (E.62), (E.64), and (E.66) occurs at most with probability $\delta$ if

[TABLE]

Restricting $n$ to be large enough so that the effective part of the bound is the square root term in both cases then yields the claim. ∎

Appendix F Technical lemmas

Lemma 23**.**

*If $B\in{\rm I}\kern-1.79993pt{\rm R}^{p\times p}$ is such that $\|B\|_{\mathrm{op}}\leq 1-\eta$ for some $\eta>0$, then we have *

[TABLE]

Moreover, for any diagonal matrix $U=\mathop{\mathsf{diag}}u$ with $u\in\{0,1\}^{p}$,

[TABLE]

Proof.

Let $B$ and $U$ as in the assumptions above. First, note that by its diagonal structure,

[TABLE]

Next, We can relate the maximum and minimum singular values to the operator norm and employ sub-additivity and sub-multiplicativity as follows:

[TABLE]

Moreover, for $i\in[p]$, denoting the standard unit vector with $1$ in the $i$th coordinate by $\mathfrak{e}_{i}$, we have

[TABLE]

and the same argument yields the bound for $\sum_{k=1}^{p}B_{i,k}^{2}$ by transposing the matrix and $\|B\|_{\mathrm{op}}=\|B^{\top}\|_{\mathrm{op}}$. ∎

Definition 24** (Sub-Gaussian and sub-Exponential random variables).**

We call a random variable $X$ sub-Gaussian with variance proxy $\sigma^{2}$, written $X\sim\mathsf{subG}(\sigma^{2})$, if

[TABLE]

We call a random variable sub-exponential with parameter $\lambda$, written $X\sim\mathsf{subE}(\lambda)$, if

[TABLE]

Lemma 25** (Product of $\mathsf{subG}$ random variables is $\mathsf{subE}$, [Ver18, Lemma 2.7.7]).**

If $X\sim\mathsf{subG}(\sigma_{X}^{2})$ and $Y\sim\mathsf{subG}(\sigma_{Y}^{2})$, then

[TABLE]

Lemma 26** (Sum of independent sub-Gaussian variables, [Ver18, Proposition 2.6.1]).**

If $X_{1},\dots,X_{n}$ are $n$ independent mean-zero random variables such that $X_{i}\sim\mathsf{subG}(\sigma_{i}^{2})$, then

[TABLE]

Lemma 27** (Bernstein’s inequality, [Ver18, Theorem 2.8.1]).**

Let $X_{1},\dots,X_{n}$ be $n$ independent mean-zero random variables such that $X_{i}\sim\mathsf{subE}(\lambda_{i})$. Then, there is an absolute constant $c_{B}$ such that for $t>0$,

[TABLE]