From graphs to DAGs: a low-complexity model and a scalable algorithm

The learning of directed acyclic graphs (DAGs) is an important problem for probabilistic and causal inference [Pea09, PJS17] with important applications in social sciences [MW15], genome research [SB09] and machine learning itself [PBM16, ABGLP19, SG21]. Through the development of probabilistic graphical models [Pea09, BPE14], DAGs are a most natural mathematical object to describe the causal relations among a number of variables. In today’s many application domains, the estimation of DAGs faces intractability issues as an ever growing number $d$ of variables is considered, due to the fact that estimating DAGs is NP-hard [Chi96]. The difficulty lies in how to enforce the acyclicity of graphs. Shimizu et al. [SHHK06] combined independent component analysis with the combinatorial linear assignment method to optimize a linear causal model (LiNGAM) and later proposed a direct and sequential algorithm [SIS⁺11] guaranteeing global optimum of the LiNGAM, for $O(d^{4})$ complexities.

Recently, Zheng et al. [ZARX18] proposed an optimization approach to learning DAGs. The breakthrough in this work, called NoTears, comes with the characterization of the DAG matrices by the zero set of a real-valued differentiable function on $\mathbb{R}^{d\times d}$, which shows that an $d\times d$ matrix $A$ is the adjacency matrix of a DAG if and only if the exponential trace satisfies

and thus the learning of DAG matrices can be cast as a continuous optimization problem subject to the constraint $h(A)=d$. The NoTears approach broadens the way of learning complex causal relations and provides promising perspectives to tackling large-scale inference problems [KGG⁺18, YCGY19, ZDA⁺20, NGZ20]. However, NoTears is still not suitable for large-scale applications as the complexity of computing the exponential trace and its gradient is $O(d^{3})$. More recently, Fang et al. [FZZ⁺20] proposed to represent DAGs by low-rank matrices with both theoretical and empirical validation of the low-rank assumption for a range of graph models. However, the adaptation of the NoTears framework [ZARX18] to low-rank model still yields a complexity of $O(d^{3})$ due to the DAG characteristic function in (1).

The contribution of the paper is to propose a new computational framework to tackle the scalability issues faced by the low-rank modeling of DAGs. We notice that the Hadamard product $\odot$ in characteristic functions as in (1) poses real obstacles to scaling up the optimization of NoTears [ZARX18] and NoTears-low-rank [FZZ⁺20]. To address these difficulties, we present a low-complexity model, named Low-Rank Additive Model (LoRAM), which is a composition of low-rank matrix factorization with sparsification, and then propose a novel approximation method compatible with LoRAM to compute the gradients of the exponential trace in (1). Formally, the gradient approximation—consisting of matrix computation of the form $(A,C,B)\to(\exp(A)\odot C)B$, where $A,C\in\mathbb{R}^{d\times d}$ and $B$ is a thin low-rank matrix—is inspired from the numerical analysis of [AMH11] for the matrix action of $\exp(A)$. We apply the new method to the computation of projections from graphs to DAGs through optimization with the differentiable DAG constraint.

Empirical evidence is presented to identify the cost and the benefits of the approximation method combined with Nesterov’s accelerated gradient descent [Nes83], depending on the considered range of problem parameters (number of nodes, rank approximation, sparsity of the target graph).

The main contributions are summarized as follows:

• •

  The LoRAM model, combining a low-rank structure with a flexible sparsification mechanism, is proposed to represent DAG matrices, together with a DAG characteristic function generalizing the exponential trace of NoTears [ZARX18].

• •

  An efficient gradient approximation method, exploiting the low-rank and sparse nature of the LoRAM model, is proposed. Under the low-rank assumption ($r\leq C\ll d$), the complexity of the proposed method is quadratic ($O(d^{2})$) instead of $O(d^{3})$ as shown in Table 1. Large efficiency gains, with insignificant loss of accuracy in some cases, are demonstrated experimentally in the considered range of application.

  	Search space	Memory req.	Cost for $\nabla h$
  NoTears [ZARX18]	$\mathbb{R}^{d\times d}$	$O(d^{2})$	$O(d^{3})$
  NoTears-low-rank [FZZ+20]	$\mathbb{R}^{d\times r}\times\mathbb{R}^{d\times r}$	$O(dr)$	$O(d^{3})$
  LoRAM (ours)	$\mathbb{R}^{d\times r}\times\mathbb{R}^{d\times r}$	$O(dr)$	$O(d^{2}r)$

A graph on $d$ nodes is defined and denoted as a pair $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $|\mathcal{V}|=d$ and $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$. By default, a directed graph is simply referred to as a graph. The adjacency matrix of a graph $\mathcal{G}$, denoted as $\mathbb{A}(\mathcal{G})$, is defined as the matrix such that $[\mathbb{A}(\mathcal{G})]_{ij}=1$ if $(i,j)\in\mathcal{E}$ and $0$ otherwise. Let $A\in\mathbb{R}^{d\times d}$ be any weighted adjacency matrix of a graph $\mathcal{G}$ on $d$ nodes, then by definition, the adjacency matrix $\mathbb{A}(\mathcal{G})$ indicates the nonzeros of $A$; the adjacency matrix $\mathbb{A}(\mathcal{G})$ is also the support of $A$, denoted as $\operatorname*{supp}(A)$, i.e., $[\operatorname*{supp}(A)]_{ij}=1$ if $A_{ij}\neq 0$ and $0$ otherwise. The number of nonzeros of $A$ is denoted as $\|A\|_{0}$. The matrix $A$ is called a DAG matrix if $\operatorname*{supp}(A)$ is the adjacency matrix of a directed acyclic graph (DAG). We define, by convention, the set of DAG matrices as follows: $\mathcal{D}_{d\times d}=\{A\in\mathbb{R}^{d\times d}:\operatorname*{supp}(A)\text{~{}defines a DAG}\}$.

We recall the following theorem that characterizes acyclic graphs using the matrix exponential trace—$\operatorname*{tr}(\exp(\cdot))$—where $\exp(\cdot)$ denotes the matrix exponential function. The matrix exponential will be denoted as $e^{\cdot}$ and $\exp(\cdot)$ indifferently. The operator $\odot$ denotes the matrix Hadamard product that acts on two matrices of the same size by elementwise multiplications.

The following corollary is a straightforward extension of the theorem above:

In view of the property above, we will refer to the composition of $\operatorname*{tr}(\exp(\cdot))$ and the operator $\sigma$ as a DAG characteristic function. Next, we show some more properties (proof in Appendix A) of the exponential trace.

In this section, we describe a low-complexity matrix representation of the adjacency matrices of directed graphs, and then a generalized DAG characteristic function for the new matrix model.

In the spirit of searching for best low-rank singular value decompositions and taking inspiration from [FZZ⁺20], the search of a full $d\times d$ (DAG) matrix $A$ is replaced by the search of a pair of thin factor matrices $(X,Y)$, in $\mathbb{R}^{d\times r}\times\mathbb{R}^{d\times r}$ for $1\leq r<d$, and the $d\times d$ candidate graph matrix is represented by the product $X{Y}^{\mathrm{T}}$. This matrix product has a rank bounded by $r$, with number of parameters $2dr\leq d^{2}$. However, the low-rank representation $(X,Y)\to X{Y}^{\mathrm{T}}$ generally gives a dense $d\times d$ matrix. Since in many scenarios the sought graph (or Bayesian network) is usually sparse, we apply a sparsification operator on $X{Y}^{\mathrm{T}}$ in order to trim abundant entries in $X{Y}^{\mathrm{T}}$. Accordingly, we combine the two operations and introduce the following model.

The candidate set $\Omega$ is to be fixed according to the specific problem. In the case of projection from a given graph to the set of DAGs, $\Omega$ can be fixed as the index set of the given graph’s edges.

The DAG characteristic function on the LoRAM search space is defined as follows:

Note that operators $\sigma_{2}$ and $\sigma_{\mathrm{abs}}$ (4) are two natural choices that meet the conditions (i)–(ii) of Corollary 2, since they both produce a nonnegative surrogate matrix of the $d\times d$ matrix $A_{\Omega}$ while preserving the support of $A_{\Omega}$.

Given a (generally non-acyclic) graph matrix $Z_{0}\in\mathbb{R}^{d\times d}$, let us consider the projection of $Z_{0}$ onto the feasible set (5):

where $A_{\Omega}(X,Y)$ is the LoRAM matrix (2), $h$ is given by Definition 5, and $\delta_{\epsilon}>0$ is a tolerance parameter. Based on Proposition 7 and given the objective function of (6), the solution to (6) provides a quasi DAG matrix closest to $Z_{0}$ and thus enables finding a projection of $Z_{0}$ onto the DAG matrix space $\mathcal{D}_{d\times d}$. More precisely, we tackle problem (6) using the penalty method and focus on the primal problem, for a given penalty parameter $\lambda>0$, followed by elementwise hard thresholding:

where $\mathbb{T}_{\epsilon^{\star}_{r,\Omega}}(z)=z\delta_{|z|\geq\epsilon^{\star}_{r,\Omega}}$ is the elementwise hard thresholding operator. The choice of $\lambda$ and $\epsilon^{\star}_{r,\Omega}$ is discussed in Section C.1 in the context of problem (6).

Problems (6) and (7) are nonconvex due to: (i) the matrix exponential trace $A\to\operatorname*{tr}(\exp(A))$ in $h$ (3) is nonconvex (as in [ZARX18]), and (ii) the matrix product $(X,Y)\to X{Y}^{\mathrm{T}}$ in $A_{\Omega}(X,Y)$ (2) is nonconvex. In view of the DAG characteristic constraint of (6), the augmented Lagrangian algorithms of NoTears [ZARX18] and NoTears-low-rank [FZZ⁺20] can be applied for the same objective as problem (6); it suffices for NoTears and NoTears-low-rank to replace the LoRAM matrix $A_{\Omega}(X,Y)$ in (6) by the $d\times d$ matrix variable and the dense matrix product $X{Y}^{\mathrm{T}}$ respectively.

However, the NoTears-based methods of [ZARX18, FZZ⁺20] have an $O(d^{3})$ complexity due to the composition of the elementwise operations (the Hadamard product $\odot$) with the matrix exponential in the $h$ function (1) and (3). We elaborate this argument in the next subsection and then propose a new computational method for computations involving the gradient of the DAG characteristic function $h$.

This section investigates the performance (computational gains and accuracy loss) of the proposed gradient approximation method (Algorithm 1) and thereafter reports on the performance of the LoRAM projection (6), compared to NoTears [ZARX18]. Sensitivity to the rank parameter $r$ of the proposed method is also investigated.

The implementation is available at https://github.com/shuyu-d/loram-exp.

This paper tackles the projection of matrices on DAG matrices, motivated by the identification of linear causal graphs. The line of research built upon the LiNGAM algorithms [SHHK06, SIS⁺11] has recently been revisited through the formal characterization of DAGness in terms of a continuously differentiable constraint by [ZARX18]. The NoTears approach of [ZARX18] however suffers from an $O(d^{3})$ complexity in the number $d$ of variables, precluding its usage for large-scale problems.

Unfortunately, this difficulty is not related to the complexity of the model (number of parameters of the model): the low-rank approach investigated by NoTears-low-rank [FZZ⁺20] also suffers from the same $O(d^{3})$ complexity, incurred in the gradient-based optimization phase.

The present paper addresses this difficulty by combining a sparsification mechanism with the low-rank model and using a new approximated gradient computation. This approximated gradient takes inspiration from the approach of [AMH11] for computing the action of exponential matrices based on truncated Taylor expansion. This approximation eventually yields a complexity of $O(d^{2}r)$, where the rank parameter is small ($r\leq C\ll d$). The experimental validation of the approach shows that the approximated gradient entails no significant error with respect to the exact gradient, for LoRAM matrices with a bounded norm, in the considered range of graph sizes ($d$) and sparsity levels. The proposed algorithm combining the approximated gradient with a Nesterov’s acceleration method [Nes83, Nes04] yields gains of orders of magnitude in computation time compared to NoTears on the same artificial benchmark problems. The approximation performance indicators reveal almost no performance loss for the projection problem in the setting of case (b) (where the matrix to be projected is perturbed with anti-causal links), while it incurs minor losses in terms of false discovery rate (FDR) in the setting of case (a) (with random additive spurious links).

Further developments aim to extend the approach and address the identification of causal DAG matrices from observational data. A longer term perspective is to extend LoRAM to the non-linear case, building upon the introduction of latent causal variables and taking inspiration from the non-linear independent component analysis and generalized contrastive losses [HST19]. Another perspective relies on the use of auto-encoders to yield a compressed representation of high-dimensional data, while constraining the structure of the encoder and decoder modules to enforce the acyclic property.

The authors warmly thank Fujitsu Laboratories LTD who funded the first author, and in particular Hiroyuki Higuchi and Koji Maruhashi for many discussions.