# Dirac: a command-line $\gamma$-matrix calculator

Sergii Kutnii

###### Abstract

A software for simplification of Dirac matrix polynomials that arise in
particle physics problems is implemented.

††journal: Computer Physics Communications

m m m⟨#1|#2|#3⟩ m m⟨#1|#2⟩ mD#1 > , m m ∫#1exp{i∫d^4x[#2]} > , m ∫#1

##  1 Introduction

Many problems in high-energy physics require simplification of polynomials of
Dirac matrices. As an example, consider the problem of classifying possible
self-interactions of relativistic fermions. Interaction lagrangian has to be a
Lorentz scalar. Let (Itzykson and Zuber, 2012, 48–54)

$\begin{split}&\Gamma_{1}=1\\\ &\Gamma^{\mu}_{2}=\gamma^{\mu}\\\
&\Gamma^{\mu_{1}\mu_{2}}_{3}=\sigma^{\mu_{1}\mu_{2}}\\\
&\Gamma^{\mu}_{4}=\gamma^{5}\gamma^{\mu}\\\ &\Gamma_{5}=\gamma^{5}\end{split}$
(1)

Then a generic quartic self-interaction term can be written as

$\mathcal{L}_{4Aij}=t_{A\bar{\mu}\bar{\nu}}\bar{\psi}\Gamma_{i}^{\bar{\mu}}\psi\bar{\psi}\Gamma_{j}^{\bar{\nu}}\psi,$
(2)

where $\bar{\mu},\bar{\nu}$ are Lorentz multi-indices, i.e. combinations of
zero to two indices, and $t_{Aij\bar{\mu}\bar{\nu}}$ is some Lorentz-invariant
tensor with up to four indices. By virtue of Weyl’s theorem on invariants of
orthogonal groups (Weyl, 1946, 53), such tensors can only be built from the
metric and Levi-Civita symbol. If parity symmetry is required, then only
metric is allowed. Not all terms of the form (2) are independent though. The
matrices (1) form a complete basis in the algebra of $4\times 4$ complex
matrices $gl(4,\mathbb{C})$. It is easy to show that for the Grassmann field
$\psi$

$\psi\bar{\psi}=-\frac{1}{4}\bar{\psi}\psi-\frac{1}{4}\gamma_{\mu}\bar{\psi}\gamma^{\mu}\psi-\frac{1}{8}\sigma_{\mu\nu}\bar{\psi}\sigma^{\mu\nu}\psi+\frac{1}{4}\gamma^{5}\gamma_{\mu}\bar{\psi}\gamma^{5}\gamma^{\mu}\psi-\frac{1}{4}\gamma^{5}\bar{\psi}\gamma^{5}\psi.$
(3)

Therefore,

$\bar{\psi}\Gamma_{i}^{\bar{\lambda}}\psi\bar{\psi}\Gamma_{k}^{\bar{\nu}}\psi=\Lambda_{j}\bar{\psi}\Gamma_{i}^{\bar{\lambda}}\Gamma_{j\bar{\mu}}\Gamma_{k}^{\bar{\nu}}\psi\bar{\psi}\Gamma_{j}^{\bar{\mu}}\psi,$
(4)

where summation over repeated indices is implied, and $\Lambda_{j}$ are
coefficients in the expansion (3). But completeness also implies the
decomposition

$\Gamma_{i}^{\bar{\lambda}}\Gamma_{j\bar{\mu}}\Gamma_{k}^{\bar{\nu}}=K_{ijkl\bar{\mu}\bar{\kappa}}^{\bar{\lambda}\bar{\nu}}\Gamma_{l}^{\bar{\kappa}},$
(5)

for certain numeric coefficients
$K_{ijkl\bar{\mu}\bar{\kappa}}^{\bar{\lambda}\bar{\nu}}$. Therefore, quartic
scalar combinations of fermions are subject to identities

$\bar{\psi}\Gamma_{i}^{\bar{\kappa}}\psi\bar{\psi}\Gamma_{j}^{\bar{\lambda}}\psi={f_{ijkl}}^{\bar{\kappa}\bar{\lambda}}_{\bar{\mu}\bar{\nu}}\bar{\psi}\Gamma_{k}^{\bar{\mu}}\psi\bar{\psi}\Gamma_{l}^{\bar{\nu}}\psi.$
(6)

These are the well-known Fierz identities (Fierz, 1937), (Itzykson and Zuber,
2012, 161–162). They can be derived manually with relative ease.

But completeness also implies an analogous set of identities for sixth-order
terms

$I_{ijk}^{\bar{\lambda}\bar{\mu}\bar{\nu}}=\bar{\psi}\Gamma_{i}^{\bar{\lambda}}\psi\bar{\psi}\Gamma_{j}^{\bar{\mu}}\psi\bar{\psi}\Gamma_{k}^{\bar{\nu}}\psi.$
(7)

The need to consider such terms arose in a study of bosonization in a Nambu-
Jona-Lasinio-like model with additional color symmetry (Kutnii, 2023b). The
idea of bosonization, proposed by Kikkawa (Kikkawa, 1976), Ebert, and
Reinhardt (Ebert and Reinhardt, 1978), (Reinhardt, 1978), is to utilize the
relation

$\exp\left[\lambda{}\int\,A^{2}(x)d^{d}x\right]\sim\funint{B}\exp\left\\{\int\left[\frac{1}{4\lambda}B(x)^{2}+A(x)B(x)\right]d^{d}x\right\\},$
(8)

a functional analog of Fourier or Laplace transform, to represent a generic
four-fermion interaction term as

$\begin{split}&\exp\left[-i\int{}G_{\bar{i}}\bar{\psi}\tau_{\bar{i}\bar{\alpha}}\psi\bar{\psi}\tau_{\bar{i}}^{\bar{\alpha}}\psi{}d^{d}x\right]\sim\\\
&\sim\funint{\Sigma}\exp\left\\{-i\int\left[-\frac{1}{4G_{\bar{i}}}\Sigma_{\bar{i}\bar{\alpha}}\Sigma_{\bar{i}}^{\bar{\alpha}}+\Sigma_{\bar{i}\bar{\alpha}}\bar{\psi}\tau_{\bar{i}}^{\bar{\alpha}}\psi\right]d^{d}x\right\\}.\end{split}$
(9)

Here

$\begin{split}&\tau_{i0}^{\bar{\alpha}}=\Gamma_{i}^{\bar{\alpha}}\\\
&\tau_{i1}^{\bar{\alpha}}=T^{a}\Gamma_{i}^{\bar{\beta}},\end{split}$ (10)

where $T^{a}$ are the generators of $su(N)$ algebra, and shorthand notation
$\bar{i}$ is used for index pairs $i0,i1$. The generating functional of the
model can thus be expressed as

$\begin{split}&Z\left[\bar{\eta}\eta\right]\equiv\funint{\bar{\psi},\psi}\exp\left\\{-i\int\left[\bar{\psi}\widehat{\partial}\psi+G_{\bar{i}}\bar{\psi}\tau_{\bar{i}\bar{\alpha}}\psi\bar{\psi}\tau_{\bar{i}}^{\bar{\alpha}}\psi+\bar{\eta}\psi+\bar{\psi}\eta\right]d^{4}x\right\\}\sim\\\
&\sim\funint{\bar{\psi},\psi,\Sigma}\exp\left\\{-i\int\left[\bar{\psi}\left(\widehat{\partial}+\widehat{\Sigma}\right)\psi-\frac{1}{4G_{\bar{i}}}\Sigma_{\bar{i}\bar{\alpha}}\Sigma_{\bar{i}}^{\bar{\alpha}}+\bar{\eta}\psi+\bar{\psi}\eta\right]d^{4}x\right\\}\end{split}$
(11)

where
$\widehat{\Sigma}=\tau_{\bar{i}\bar{\alpha}}\Sigma_{\bar{i}}^{\bar{\alpha}}$.
The quartic Fierz identities (6) together with first fundamental theorem of
invariant theory for $su(N)$ matrices (Kraft and Procesi, 1996, 21–22) have
two major implications: that the most general massless parity-preserving NJL-
like lagrangian with global $U(N)$ symmetry is

$\begin{split}\mathcal{L}(\bar{\psi},\psi)&=\bar{\psi}i\widehat{\partial}\psi+G_{1}\bar{\psi}\psi\bar{\psi}\psi+G_{2}\bar{\psi}\gamma_{\mu}\psi\bar{\psi}\gamma^{\mu}\psi+\\\
&+G_{3}\bar{\psi}\sigma_{\mu\nu}\psi\bar{\psi}\sigma^{\mu\nu}\psi+G_{4}\bar{\psi}\gamma^{5}\gamma_{\mu}\psi\bar{\psi}\gamma^{5}\gamma^{\mu}\psi+G_{5}\bar{\psi}\gamma^{5}\psi\bar{\psi}\gamma^{5}\psi,\end{split}$
(12)

and that bosonization is not uniquely determined by interaction. The
lagrangian

$\begin{split}&\mathcal{L}(\bar{\psi},\psi,\widehat{\Sigma})=\bar{\psi}\left(i\widehat{\partial}+\widehat{\Sigma}\right)\psi-\frac{1}{4\lambda_{10}}\Sigma_{10}^{2}-\frac{1}{4\lambda_{11}}\Sigma_{11}^{a}\Sigma_{11}^{a}-\\\
&-\frac{1}{4\lambda_{20}}\Sigma_{20\,\mu}\Sigma^{\mu}_{20}-\frac{1}{4\lambda_{21}}\Sigma^{a}_{21\,\mu}\Sigma^{a\mu}_{21}-\frac{1}{4\lambda_{30}}\Sigma_{30\,\mu\nu}\Sigma^{\mu\nu}_{30}-\frac{1}{4\lambda_{31}}\Sigma^{a}_{31\,\mu\nu}\Sigma^{a\mu\nu}_{31}-\\\
&-\frac{1}{4\lambda_{40}}\Sigma_{40\,\mu}\Sigma^{\mu}_{40}-\frac{1}{4\lambda_{41}}\Sigma^{a}_{41\,\mu}\Sigma^{a\mu}_{41}-\frac{1}{4\lambda_{50}}\Sigma_{50}^{2}-\frac{1}{4\lambda_{51}}\Sigma^{a}_{51}\Sigma^{a}_{51}\\\
\widehat{\Sigma}&=\Sigma_{10}+T^{a}\Sigma^{a}_{11}+\gamma_{\mu}\Sigma^{\mu}_{20}+T^{a}\gamma_{\mu}\Sigma^{a\mu}_{21}+\\\
&+\sigma_{\mu\nu}\Sigma^{\mu\nu}_{30}+T^{a}\sigma_{\mu\nu}\Sigma^{a\mu\nu}_{31}+\gamma^{5}\gamma_{\mu}\Sigma^{\mu}_{40}+T^{a}\gamma^{5}\gamma_{\mu}\Sigma^{a\mu}_{41}+\gamma^{5}\Sigma_{50}+T^{a}\gamma^{5}\Sigma^{a}_{51}\end{split}$
(13)

defines bosonization of (12), provided that

$\begin{bmatrix}G_{1}\\\ G_{2}\\\ G_{3}\\\ G_{4}\\\
G_{5}\end{bmatrix}=\begin{bmatrix}\lambda_{10}\\\ \lambda_{20}\\\
\lambda_{30}\\\ \lambda_{40}\\\
\lambda_{50}\end{bmatrix}+\begin{bmatrix}-\frac{N+4}{8N}&-\frac{1}{2}&-\frac{3}{2}&\frac{1}{2}&-\frac{1}{8}\\\
-\frac{1}{2}&\frac{N-2}{4N}&0&\frac{1}{4}&\frac{1}{8}\\\
-\frac{1}{16}&0&\frac{N-2}{4N}&0&-\frac{1}{16}\\\
\frac{1}{8}&\frac{1}{4}&0&\frac{N-2}{4N}&-\frac{1}{8}\\\
-\frac{1}{8}&\frac{1}{2}&-\frac{3}{2}&-\frac{1}{2}&-\frac{N+4}{8N}\end{bmatrix}\begin{bmatrix}\lambda_{11}\\\
\lambda_{21}\\\ \lambda_{31}\\\ \lambda_{41}\\\ \lambda_{51}\end{bmatrix}.$
(14)

It is clear that half of $\lambda_{\bar{i}}$ are free parameters and the
system of equations (14) has nontrivial solutions even if all $G_{i}=0$. This
suggests a possibility to consider a generalized bosonization: to insert

$\begin{split}\text{const}&=\funint{\Sigma}\exp\left\\{-i\int\left[\frac{\lambda_{11\bar{i}}^{-\frac{1}{2}}}{2}\Sigma_{11\bar{i}\bar{\alpha}}-\lambda_{11\bar{i}}^{\frac{1}{2}}\bar{\psi}\tau_{\bar{i}\bar{\alpha}}\psi\right]\times\right.\\\
&\left.\times\left[\frac{\lambda_{11\bar{i}}^{-\frac{1}{2}}}{2}\Sigma_{11\bar{i}}^{\bar{\alpha}}-\lambda_{11\bar{i}}^{\frac{1}{2}}\bar{\psi}\tau_{\bar{i}}^{\bar{\alpha}}\psi\right]d^{4}x-\right.\\\
&\left.-i\int\left[\varkappa_{21\bar{i}}^{-\frac{1}{2}}\bar{\Sigma}_{21\bar{i}\bar{\alpha}}-\left(\varkappa_{21\bar{i}}^{\frac{1}{2}}\bar{\psi}\tau_{\bar{i}\bar{\alpha}}+\lambda_{21\bar{i}}^{\frac{1}{2}}\bar{\psi}\tau_{\bar{i}\bar{\alpha}}\psi\bar{\psi}\right)\right]\times\right.\\\
&\left.\times\left[\varkappa_{21\bar{i}}^{-\frac{1}{2}}\Sigma_{21\bar{i}}^{\bar{\alpha}}-\left(\varkappa_{21\bar{i}}^{\frac{1}{2}}\tau_{\bar{i}}^{\bar{\alpha}}\psi+\lambda_{21\bar{i}}^{\frac{1}{2}}\psi\bar{\psi}\tau_{\bar{i}}^{\bar{\alpha}}\psi\right)\right]d^{4}x\right\\}\end{split}$
(15)

into path integral expression for the generating functional. This would
naturally produce terms of type (7) (properties of $su(N)$ generators imply
that all terms can be reduced to that form). Fierz identities for such terms
will lead to cancellation conditions analogous to (14). The main obstacle to
deriving them is combinatorial explosion of complexity. There are, up to
permutations, 35 terms of the form (7) (basis matrices can occur more than
once within a single term), each of which has two $\psi\bar{\psi}$ factors
that can be expanded along the lines of (3). Thus there are 875 products of
five basis matrices, each producing up to several dozen terms.

Only Lorentz-invariant sixth-order combinations of fermions can be present in
the interaction lagrangian. The even-parity ones can be built as contractions
of $\psi\bar{\psi}$ with eleven Lorentz-invariant tensor products

$\begin{split}&\mathbf{\theta}_{111}=\mathbf{1}\otimes\mathbf{1}\otimes\mathbf{1}\\\
&\mathbf{\theta}_{122}=\mathbf{1}\otimes\gamma_{\mu}\otimes\gamma^{\mu}\\\
&\mathbf{\theta}_{133}=\mathbf{1}\otimes\sigma_{\mu\nu}\otimes\sigma^{\mu\nu}\\\
&\mathbf{\theta}_{144}=\mathbf{1}\otimes\gamma^{5}\gamma_{\mu}\otimes\gamma^{5}\gamma^{\mu}\\\
&\mathbf{\theta}_{155}=\mathbf{1}\otimes\gamma^{5}\otimes\gamma^{5}\\\
{}&\mathbf{\theta}_{523}=\gamma^{5}\otimes\gamma_{\mu}\otimes\gamma^{5}\gamma^{\mu}\\\
&\mathbf{\theta}_{223}=\gamma_{\mu}\otimes\gamma_{\nu}\otimes\sigma^{\mu\nu}\\\
&\mathbf{\theta}_{443}=\gamma^{5}\gamma_{\mu}\otimes\gamma^{5}\gamma_{\nu}\otimes\sigma^{\mu\nu}\\\
&\mathbf{\theta}_{243}=\epsilon_{\kappa\lambda\mu\nu}\gamma^{\kappa}\otimes\gamma^{5}\gamma^{\lambda}\otimes\sigma^{\mu\nu}\\\
&\mathbf{\theta}_{533}=\epsilon_{\kappa\lambda\mu\nu}\gamma^{5}\otimes\sigma^{\kappa\lambda}\otimes\sigma^{\mu\nu}\\\
&\mathbf{\theta}_{333}=\sigma_{\kappa}^{\lambda}\otimes\sigma_{\lambda}^{\mu}\otimes\sigma_{\mu}^{\kappa}\end{split}$
(16)

These lead to 275 matrix products to simplify, some of them being quite
complex, such as the product of five $\sigma$-matrices. The need for an
automated computation tool becomes clear. However, all such tools are either
proprietary with prohibitive cost or their support of tensor and Dirac matrix
algebra is tailored for different kinds of problems. Cadabra (Peeters, 2007)
could be the best existing open-source tool for the task, but it uses
different matrix basis and is designed for computations in arbitrary number of
spatio-temporal dimensions. The need to support arbitrary dimension numbers
and arbitrary metrics makes it hard to introduce $\gamma^{5}$ since the latter
is an essentially four-dimensional object not easily generalizable to
arbitrary dimension count. But simplification of Dirac matrix products
required to produce Fierz identities has to take into account identities
involving $\gamma^{5}$ too. e.g.

$\gamma_{\lambda}\gamma_{\mu}\gamma_{\nu}=\eta_{\lambda\mu}\gamma_{\nu}-\eta_{\lambda\nu}\gamma_{\mu}+\eta_{\mu\nu}\gamma_{\lambda}+i\epsilon_{\lambda\mu\nu\rho}\gamma^{5}\gamma^{\rho}$
(17)

(Pal, 2007) which is valid for four-dimensional Lorentz metric only. Cadabra,
for example, failed to produce any analog of (17). GiNaC and sympy lack full
support of the basis (1) having no predefined construct for $\sigma_{\mu\nu}$
so even if they can simplify Dirac matrix polynomials, they would make it
necessary to write some basis mapping code to present the results in the
required form. The philosophy behind dirac is different from that of the tools
mentioned above: instead of trying to solve a general problem with a general-
purpose computer algebra system, dirac simplifies $\gamma$-matrix polynomials
in four-dimensional Minkowski spacetime and produces output suitable for copy-
pasting into LaTeXdocuments. The goal was to produce the most efficient
solution for one highly specific use case, rather than competing with much
older and more feature-rich computer algebra systems.

##  2 Representation and algorithms

Completeness of the basis (1) implies that

$\Gamma_{i}^{\bar{\alpha}}\Gamma_{j}^{\bar{\beta}}={C_{ij}^{k}}^{\bar{\alpha}\bar{\beta}}_{\bar{\gamma}}\Gamma^{\bar{\gamma}}_{k}.$
(18)

The multiplication structure constants
${C_{ij}^{k}}^{\bar{\alpha}\bar{\beta}}_{\bar{\gamma}}$ can be grouped into
pseudo-matrices. Let

$\mathbf{\Gamma}^{\bar{\mu}}=\begin{bmatrix}1&\gamma^{\mu}&\sigma^{\mu_{1}\mu_{2}}&\gamma^{5}\gamma^{\mu}&\gamma^{5}\end{bmatrix}$
(19)

Then (18) is equivalent to

$\Gamma_{i}^{\bar{\mu}}\mathbf{\Gamma}^{\bar{\nu}}=\mathbf{\Gamma}^{\bar{\lambda}}\left(\mathbf{C}_{i}^{\bar{\mu}}\right)_{\bar{\lambda}}^{\bar{\nu}}.$
(20)

The nice thing about pseudo-matrices $\mathbf{C_{i}^{\bar{\mu}}}$ is that they
form a representation of $\Gamma_{i}^{\bar{\mu}}$: it is easy to verify that
multiplication of $\Gamma$ corresponds to multiplication of the respective
pseudo-matrices $\mathbf{C}$ in the same order plus contraction of matching
tensor indices. Thus, all pseudo-matrices can be constructed recursively from
basic multiplicative identities, a good reference for which can be found in
(Pal, 2007). Let the metric be denoted with $\eta$, and imaginary unit with
$I$. The pseudo-matrix counterparts to $\gamma^{\mu}$ are

$\begin{split}&\mathbf{C}_{2}^{\mu}=\begin{bmatrix}0&\eta^{\mu\nu}&0&0&0\\\
\delta^{\mu}_{\lambda}&0&I\left(\eta^{\mu\nu_{1}}\delta_{\lambda}^{\nu_{2}}-\eta^{\mu\nu_{2}}\delta_{\lambda}^{\nu_{1}}\right)&0&0\\\
0&-\frac{I}{2}\left(\delta_{\lambda_{1}}^{\mu}\delta_{\lambda_{2}}^{\nu}-\delta_{\lambda_{2}}^{\mu}\delta_{\lambda_{1}}^{\nu}\right)&0&\frac{1}{2}{\epsilon^{\mu\nu}}_{\lambda_{1}\lambda_{2}}&0\\\
0&0&-{\epsilon^{\mu\nu_{1}\nu_{2}}}_{\lambda}&0&-\delta^{\mu}_{\lambda}\\\
0&0&0&-\eta^{\mu\nu}&0\end{bmatrix}\end{split}$ (21)

and $\gamma^{5}$ is represented with

$\begin{split}&\mathbf{C}_{5}=\begin{bmatrix}0&0&0&0&1\\\
0&0&0&\delta_{\lambda}^{\nu}&0\\\
0&0&-\frac{I}{2}{\epsilon^{\nu_{1}\nu_{2}}}_{\lambda_{1}\lambda_{2}}&0&0\\\
0&\delta_{\lambda}^{\nu}&0&0&0\\\ 1&0&0&0&0\end{bmatrix}\end{split}$ (22)

Given that, one can compute

$\left(\mathbf{C}_{4}^{\mu}\right)^{\bar{\nu}}_{\bar{\lambda}}=\left(\mathbf{C}_{5}\right)^{\bar{\kappa}}_{\bar{\lambda}}\left(\mathbf{C}_{2}^{\mu}\right)_{\bar{\kappa}}^{\bar{\nu}}=\leavevmode\resizebox{216.81pt}{}{$\begin{bmatrix}0&0&0&-\eta^{\mu\nu}&0\\\
0&0&-{\epsilon^{\mu\nu_{1}\nu_{2}}}_{\lambda}&0&-\delta^{\mu}_{\lambda}\\\
0&-\frac{1}{2}{\epsilon^{\nu_{1}\nu_{2}}}_{\lambda_{1}\lambda_{2}}&0&-\frac{I}{2}\left(\delta_{\lambda_{1}}^{\mu}\delta_{\lambda_{2}}^{\nu}-\delta_{\lambda_{2}}^{\mu}\delta_{\lambda_{1}}^{\nu}\right)&0\\\
\delta^{\mu}_{\lambda}&0&I\left(\eta^{\mu\nu_{1}}\delta_{\lambda}^{\nu_{2}}-\eta^{\mu\nu_{2}}\delta_{\lambda}^{\nu_{1}}\right)&0&0\\\
0&\eta^{\mu\nu}&0&0&0\end{bmatrix}$}$ (23)

and the most complicated of all

$\begin{split}&\left(\mathbf{C}_{3}^{\mu_{1}\mu_{2}}\right)^{\bar{\nu}}_{\bar{\lambda}}=\frac{I}{2}\left[\left(\mathbf{C}_{2}^{\mu_{1}}\right)_{\bar{\lambda}}^{\bar{\kappa}}\left(\mathbf{C}_{2}^{\mu_{2}}\right)_{\bar{\kappa}}^{\bar{\nu}}-\left(\mathbf{C}_{2}^{\mu_{2}}\right)_{\bar{\lambda}}^{\bar{\kappa}}\left(\mathbf{C}_{2}^{\mu_{1}}\right)_{\bar{\kappa}}^{\bar{\nu}}\right]=\\\
&=\leavevmode\resizebox{346.89731pt}{}{$\begin{bmatrix}0&0&\eta^{\mu_{1}\nu_{1}}\eta^{\mu_{2}\nu_{2}}-\eta^{\mu_{2}\nu_{1}}\eta^{\mu_{1}\nu_{2}}&0&0\\\
0&I\left(\delta^{\mu_{1}}_{\lambda}\eta^{\mu_{2}\nu}-\delta^{\mu_{2}}_{\lambda}\eta^{\mu_{1}\nu}\right)&0&-{\epsilon^{\mu_{1}\mu_{2}\nu}}_{\lambda}&0\\\
\frac{1}{2}\left(\delta^{\mu_{1}}_{\lambda_{1}}\delta^{\mu_{2}}_{\lambda_{2}}-\delta^{\mu_{2}}_{\lambda_{1}}\delta^{\mu_{1}}_{\lambda_{2}}\right)&0&\begin{array}[]{l}\frac{I}{2}\left[\delta^{\mu_{1}}_{\lambda_{1}}\delta^{\nu_{2}}_{\lambda_{2}}\eta^{\mu_{2}\nu_{1}}-\delta^{\mu_{2}}_{\lambda_{1}}\delta^{\nu_{2}}_{\lambda_{2}}\eta^{\mu_{1}\nu_{1}}-\right.\\\
-\delta^{\mu_{1}}_{\lambda_{1}}\delta^{\nu_{1}}_{\lambda_{2}}\eta^{\mu_{2}\nu_{2}}+\delta^{\mu_{2}}_{\lambda_{1}}\delta^{\nu_{1}}_{\lambda_{2}}\eta^{\mu_{1}\nu_{2}}-\\\
-\delta^{\mu_{1}}_{\lambda_{2}}\delta^{\nu_{2}}_{\lambda_{1}}\eta^{\mu_{2}\nu_{1}}+\delta^{\mu_{2}}_{\lambda_{2}}\delta^{\nu_{2}}_{\lambda_{1}}\eta^{\mu_{1}\nu_{1}}+\\\
\left.+\delta^{\mu_{1}}_{\lambda_{2}}\delta^{\nu_{1}}_{\lambda_{1}}\eta^{\mu_{2}\nu_{2}}-\delta^{\mu_{2}}_{\lambda_{2}}\delta^{\nu_{1}}_{\lambda_{1}}\eta^{\mu_{1}\nu_{2}}\right]\end{array}&0&-\frac{I}{2}{\epsilon^{\mu_{1}\mu_{2}}}_{\lambda_{1}\lambda_{2}}\\\
0&-{\epsilon^{\mu_{1}\mu_{2}\nu}}_{\lambda}&0&I\left(\eta^{\mu_{2}\nu}\delta^{\mu_{1}}_{\lambda}-\eta^{\mu_{1}\nu}\delta^{\mu_{2}}_{\lambda}\right)&0\\\
0&0&-I\epsilon^{\mu_{1}\mu_{2}\nu_{1}\nu_{2}}&0&0\end{bmatrix}$}\end{split}$
(24)

Then any product of $\Gamma$-matrices can be represented as follows:

$\Gamma_{i_{1}}^{\bar{\mu}_{1}}\ldots\Gamma_{i_{k}}^{\bar{\mu}_{k}}\ldots\Gamma_{i_{n}}^{\bar{\mu}_{n}}\rightarrow{\left(\mathbf{C}_{i_{1}}^{\mu_{1}}\right)}^{\bar{\lambda}_{2}}_{\bar{\lambda}_{1}}\ldots{\left(\mathbf{C}_{i_{k}}^{\mu_{k}}\right)}_{\bar{\lambda}_{k}}^{\bar{\lambda}_{k+1}}\ldots{\left(\mathbf{C}_{i_{n}}^{\mu_{n}}\right)}_{\bar{\lambda}_{n}}^{\bar{\lambda}_{n+1}}.$
(25)

On the other hand, any matrix in $gl(4,\mathbb{C})$ can be written as product
of the basis row (19) and a coefficient column vector. In particular,

$1\rightarrow\begin{bmatrix}1\\\ 0\\\ 0\\\ 0\\\ 0\end{bmatrix}.$ (26)

Any matrix can be multiplied by $1$. This implies that column vector
representation of the left hand side in (25) can be obtained by multiplying
its pseudo-matrix counterpart with the right hand side of (26) on the right,
which is equivalent to taking the first column of the rightmost pseudo-matrix
and multiplying it by the remaining pseudo-matrices on the left in the right
to left order.

All that remains is to simplify the components of the resulting column vector
which are polynomials in the metric, Kronecker, and Levi-Civita symbols with
complex coefficients. This can be done in three stages. First, Levi-Civita
powers can be expanded using

$\epsilon_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}=-\sum_{P}\text{sgn}(P)\eta_{\mu_{1}P\nu_{1}}\eta_{\mu_{2}P\nu_{2}}\eta_{\mu_{3}P\nu_{3}}\eta_{\mu_{4}P\nu_{4}},$
(27)

where $P$ labels all permutations of $\nu_{1},\ldots,\nu_{4}$:

$P\left\\{\mu_{1},\mu_{2},\mu_{3},\mu_{4}\right\\}=\left\\{P\mu_{1},P\mu_{2},P\mu_{3},P\mu_{4}\right\\},$
(28)

so if
$P\left\\{\mu_{1},\mu_{2},\mu_{3},\mu_{4}\right\\}=\left\\{\mu_{2},\mu_{3},\mu_{4},\mu_{1}\right\\}$
then $P\mu_{1}=\mu_{2},P\mu_{2}=\mu_{3}$ and so on. Then all contractible
indices can be contracted. using the fact that contraction of any tensor with
the metric or Kronecker amounts to simple index replacement, and Levi-Civita
symbol with any index pair contracted is zero. Contraction guarantees that all
remaining indices are free. Then all terms having same tensor composition up
to permutations of indices in constituent tensors can be collected, taking
into account Levi-Civita index permutation signs.

The only remaining major part of a $\Gamma$-matrix calculator is expression
parser. A version of shunting yard algorithm was implemented. Input syntax is
described in the next section.

##  3 The software

The dirac software is available on github (Kutnii, 2023a). It is implemented
in C++20 and uses CMake (min version 3) as its build system. Ruby 3 language
interpreter, ‘‘open3’’ Ruby gem, and LaTeX installation with `pdflatex` are
required to run example scripts. The standard CMake build procedure (from the
project root)

⬇

mkdir build

cd build

cmake ..

make

produces ‘‘dirac’’ executable in the ‘‘build’’ folder.

There are two main ways to invoke the executable. The first is to provide the
expression to simplify with `-e` key. In this case the executable simplifies
the expression, prints the result to the console and exits. Example:

⬇

./dirac -e "\gamma_\mu\gamma_\nu"

\eta_{\nu\mu} -
I\eta_{\mu\omega_{1}}\eta_{\nu\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

(the quotes around the expression prevent the terminal from eating the
backslashes). This mode is mainly useful for scripting. The ‘‘examples’’
folder contains Ruby scripts that demonstrate how scripting can be used to
batch computations and generate LaTeX output from the results.

There are other command line keys as well. They affect input processing and
output formatting.

If no expression is provided via command line, the application starts an
interactive shell. The shell accepts three primary types of input.

Quit-expression is the single word `quit` which exits the shell.

Set-expression:

⬇

dirac:> #set <var-name> <var-value>

sets the variable name identified by `var-name` to `var-value`. The variables
are documented below.

mode: arithmetic mode. Possible values: `float` and `rational`. Default is
`rational`. Equivalent command line option `-m`. Example:

⬇

dirac:> #set mode float

dirac:> \gamma5\sigma^{\mu\nu}

-0.500000I{\epsilon^{\mu\nu}}_{\omega_{1}\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

dirac:> #set mode rational

dirac:> \gamma5\sigma^{\mu\nu}

- \frac{I}{2}{\epsilon^{\mu\nu}}_{\omega_{1}\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

or

⬇

./dirac -m float

The `mode` variable also affects input parsing. Floating point numeric values
are acceptable in float mode and are considered errors when the mode is
rational.

The rational arithmetic implementation is the simplest possible. A rational
number is implemented as a pair of two largest system integer types. Since

$\frac{a}{b}+\frac{c}{d}=\frac{ad+bc}{bd},$ (29)

longer rational expressions will tend to result in larger numerators and
denominators. Thus, the implementation is prone to integer overflows in very
large $\gamma$-matrix expressions, when the result of an operation on integers
is larger than the largest integer representable by the machine. It is not
clear, however, whether such overflows will happen in any practical example.
Therefore a conscious decision was made to retain the naive rational
implementation to avoid either adding library dependencies or implementing
arbitrary precision arithmetic from scratch. A more robust implementation may
be added in the future if there is user demand for it.

line_terms: number of terms per output line. Possible values: integers or
`inf` (meaning ‘‘infinity’’). Default: `inf`. When this variable is set to a
nonzero integer constant, LaTeX line breaks (as if inside ‘split’ environment)
are inserted after each `line_terms` terms. Example:

⬇

dirac:> #set line_terms 2

dirac:> \gamma_\kappa\gamma_\lambda\gamma_\mu\gamma_\nu

&\eta_{\lambda\kappa}\eta_{\nu\mu} + \eta_{\mu\lambda}\eta_{\nu\kappa} -\\\

&-\eta_{\mu\kappa}\eta_{\nu\lambda} + \left[ -
I\eta_{\kappa\omega_{1}}\eta_{\lambda\omega_{2}}\eta_{\nu\mu} -\right.\\\

&\left. - I\eta_{\kappa\omega_{1}}\eta_{\mu\lambda}\eta_{\nu\omega_{2}} +
I\eta_{\kappa\omega_{1}}\eta_{\mu\omega_{2}}\eta_{\nu\lambda} + \right.\\\

&\left.+I\eta_{\kappa\lambda}\eta_{\omega_{1}\nu}\eta_{\omega_{2}\mu} -
I\eta_{\kappa\mu}\eta_{\omega_{1}\nu}\eta_{\omega_{2}\lambda} + \right.\\\

&\left.+I\eta_{\kappa\nu}\eta_{\omega_{1}\mu}\eta_{\omega_{2}\lambda}\right]\sigma^{\omega_{1}\omega_{2}}
- I\epsilon_{\lambda\mu\nu\kappa}\gamma^5

dirac:> #set line_terms inf

dirac:> \gamma_\kappa\gamma_\lambda\gamma_\mu\gamma_\nu

\eta_{\lambda\kappa}\eta_{\nu\mu} + \eta_{\mu\lambda}\eta_{\nu\kappa}
-\eta_{\mu\kappa}\eta_{\nu\lambda} + \left[ -
I\eta_{\kappa\omega_{1}}\eta_{\lambda\omega_{2}}\eta_{\nu\mu} -
I\eta_{\kappa\omega_{1}}\eta_{\mu\lambda}\eta_{\nu\omega_{2}} +
I\eta_{\kappa\omega_{1}}\eta_{\mu\omega_{2}}\eta_{\nu\lambda} +
I\eta_{\kappa\lambda}\eta_{\omega_{1}\nu}\eta_{\omega_{2}\mu} -
I\eta_{\kappa\mu}\eta_{\omega_{1}\nu}\eta_{\omega_{2}\lambda} +
I\eta_{\kappa\nu}\eta_{\omega_{1}\mu}\eta_{\omega_{2}\lambda}\right]\sigma^{\omega_{1}\omega_{2}}
- I\epsilon_{\lambda\mu\nu\kappa}\gamma^5

Equivalent command line option: `-l`:

⬇

./dirac -l 2

dummy: dummy index template. Possible values: string literals. Greek letters
in LaTeX notation are good choice. Default is `\omega`. Equivalent command
line option: `-d`. Example:

⬇

dirac:> \gamma_\mu\gamma_\nu

\eta_{\nu\mu} -
I\eta_{\mu\omega_{1}}\eta_{\nu\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

dirac:> #set dummy \sigma

dirac:> \gamma_\mu\gamma_\nu

\eta_{\nu\mu} -
I\eta_{\mu\sigma_{1}}\eta_{\nu\sigma_{2}}\sigma^{\sigma_{1}\sigma_{2}}

or

⬇

./dirac -d "\sigma"

apply_symmetry: controls whether the coefficient terms at $\sigma^{\mu\nu}$
are merged in the output by taking into account $\sigma$’s antisymmetry.
Possible values: `true` or `false`. Default is `true`. Command line
equivalent: `-s`. Example:

⬇

dirac:> #set apply_symmetry false

dirac:> \gamma_\mu\gamma_\nu

\eta_{\nu\mu} + \left[ - \frac{I}{2}\eta_{\mu\omega_{1}}\eta_{\nu\omega_{2}} +
\frac{I}{2}\eta_{\mu\omega_{2}}\eta_{\nu\omega_{1}}\right]\sigma^{\omega_{1}\omega_{2}}

dirac:> #set apply_symmetry true

dirac:> \gamma_\mu\gamma_\nu

\eta_{\nu\mu} -
I\eta_{\mu\omega_{1}}\eta_{\nu\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

or

⬇

./dirac -s false

###  3.1 Math expressions

All input lines that are neither quit-expressions nor set-expressions are
considered computable math. The dirac application tries to parse and compute
them. Math syntax is LaTeX-like with some differences. A valid expression
consists of

  * 1.

literals: alphanumeric sequences preceded by `\`, e.g. `\gamma`;

  * 2.

integer or floating-point numbers;

  * 3.

arithmetic operators `+,-,*,/`;

  * 4.

subscript `_`;

  * 5.

superscript `^`;

  * 6.

brackets `{...}` (round or square brackets are not recognized for the sake of
implementation simplicity).

Subscript `[head]_[tail]` and superscript `[head]^[tail]` are tensorial
expressions. Unlike LaTeX, multiple non-bracketed subscripts and superscripts
to a single head are possible, but multi-level are not: `\eta_\mu_\nu` is
valid input while `\gamma^{\eta_{\mu\nu}}` is not.

Literals’ interpretation depends on their position in the input. Literals
inside the tail of a tensorial expression are interpreted verbatim as tensor
index labels. Otherwise, only a limited number of special literals is
recognized:

  * 1.

`\I` \- imaginary unit;

  * 2.

`\gamma` \- Dirac gamma-matrix;

  * 3.

`\sigma` \- Dirac sigma-matrix;

  * 4.

`\gamma5` \- $\gamma^{5}$ matrix;

  * 5.

`\eta` \- Minkowski metric;

  * 6.

`\delta` \- Kronecker delta;

  * 7.

`\epsilon` \- Levi-Civita symbol.

The app does not perform any validation of tensorial expression consistency
save for checking that all basic tensors have correct index counts at
computation stage. Nevertheless, mathematically valid input results in
mathematically valid output by design of the pseudo-matrix representation.

Some examples:

⬇

dirac:> \gamma_\lambda\gamma_\mu\gamma_\nu\gamma^\lambda

4\eta_{\mu\nu}

dirac:> \gamma_\kappa\gamma_\lambda\gamma_\mu\gamma_\nu\gamma^\kappa

\left[-2\eta_{\mu\lambda}\eta_{\omega_{1}\nu}
-2\eta_{\nu\mu}\eta_{\omega_{1}\lambda} +
2\eta_{\nu\lambda}\eta_{\omega_{1}\mu}\right]\gamma^{\omega_{1}} +
2I\epsilon_{\nu\lambda\mu\omega_{1}}\gamma^5\gamma^{\omega_{1}}

dirac:> \gamma_\mu\gamma_\nu - \gamma_\nu\gamma_\mu

- 2I\eta_{\mu\omega_{1}}\eta_{\nu\omega_{2}}\sigma^{\omega_{1}\omega_{2}}

###  3.2 Dirac as a library

It is also possible to use dirac as a library. Reusable C++ classes and
routines are defined in the `src/algebra` directory. `CMakeLists.txt` defines
a library target dirac_common that includes those components. Top-level
constructs are defined in `src/algebra/Gamma.hpp`. These are

  * 1.

`GammaPolynomial` \- gamma-ring element i.e. a polynomial of Dirac matrices
and Lorentz-invariant symbols with complex coefficients;

  * 2.

`CanonicalExpr` \- canonical gamma-expression i.e. a linear combination of
$1,\gamma^{\mu},\sigma^{\mu\nu},\gamma^{5}\gamma^{\mu},\gamma^{5}$ with
Lorentz-invariant coefficients;

  * 3.

`reduceGamma` \- main simplification routine that converts a `GammaPolynomial`
to a `CanonicalExpr`.

In the `examples/fierz_gen` directory dirac_common library is applied to the
problem of computing sixth-order Fierz identities stated in the introduction.

##  4 Conclusions

The dirac software is a fully functional command line calculator for
$\gamma$-matrix polynomials. Scripting can be used to process multiple
expressions in a batch ang generate LaTeX documents with the output.

## References

  * Ebert and Reinhardt (1978) Ebert, D., Reinhardt, H., 1978\.  Functional approach to nuclear field theory: A schematic model with pairing and particle-hole forces.  Nuclear Physics A 298, 60–76.  URL: https://www.sciencedirect.com/science/article/pii/0375947478900076, doi:https://doi.org/10.1016/0375-9474(78)90007-6.
  * Fierz (1937) Fierz, M., 1937.  Zur fermischen theorie des $\beta$-zerfalls.  Zeitschrift für Physik 104, 553–565.
  * Itzykson and Zuber (2012) Itzykson, C., Zuber, J., 2012\.  Quantum Field Theory.  Dover Books on Physics, Dover Publications.
  * Kikkawa (1976) Kikkawa, K., 1976.  Quantum corrections in superconductor models.  Progress of Theoretical Physics 56, 947–955.  URL: https://doi.org/10.1143/PTP.56.947, doi:10.1143/PTP.56.947, arXiv:https://academic.oup.com/ptp/article-pdf/56/3/947/5455249/56-3-947.pdf.
  * Kraft and Procesi (1996) Kraft, H., Procesi, C., 1996\.  A Primer of Classical Invariant Theory: Preliminary Version.  Preprint.
  * Kutnii (2023a) Kutnii, S., 2023a.  URL: https://github.com/skutnii/dirac.
  * Kutnii (2023b) Kutnii, S., 2023b.  Bosonization, effective action, and R-operation in a generalized Nambu-Jona-Lasinio model.  arXiv:2304.07118.
  * Pal (2007) Pal, P.B., 2007.  Representation-independent manipulations with dirac matrices and spinors.  URL: https://arxiv.org/abs/physics/0703214, doi:10.48550/ARXIV.PHYSICS/0703214.
  * Peeters (2007) Peeters, K., 2007.  Cadabra: a field-theory motivated symbolic computer algebra system.  Computer Physics Communications 176, 550–558.  URL: https://www.sciencedirect.com/science/article/pii/S0010465507000318, doi:https://doi.org/10.1016/j.cpc.2007.01.003.
  * Reinhardt (1978) Reinhardt, H., 1978.  Nuclear field theory.  Nuclear Physics A 298, 77–92.  URL: https://www.sciencedirect.com/science/article/pii/0375947478900088, doi:https://doi.org/10.1016/0375-9474(78)90008-8.
  * Weyl (1946) Weyl, H., 1946.  The classical groups: their invariants and representations.  1, Princeton university press.


6 September, 2023 23 December, 2023 26 December, 2023 2 January, 2024

Corresponding author: Petar M. Djurić<EMAIL_ADDRESS>This work
was supported by the National Science Foundation under Awards 2212506 and
2021002.

#  Dagma-DCE: Interpretable, Non-Parametric Differentiable Causal Discovery

Daniel Waxman ${}^{\textbf{1}}$ 0009-0004-0168-5547 Graduate Student Member
IEEE  Kurt Butler ${}^{\textbf{1}}$ 0000-0002-1520-4909 Graduate Student
Member  IEEE  and Petar M. Djurić ${}^{\textbf{1}}$ 0000-0001-7791-3199 Life
Fellow  IEEE  Department of Electrical and Computer Engineering, Stony Brook
University, Stony Brook, NY 11794 USA

###### Abstract

We introduce Dagma-DCE, an interpretable and model-agnostic scheme for
differentiable causal discovery. Current non- or over-parametric methods in
differentiable causal discovery use opaque proxies of “independence” to
justify the inclusion or exclusion of a causal relationship. We show
theoretically and empirically that these proxies may be arbitrarily different
than the actual causal strength. Juxtaposed to existing differentiable causal
discovery algorithms, Dagma-DCE uses an interpretable measure of causal
strength to define weighted adjacency matrices. In a number of simulated
datasets, we show our method achieves state-of-the-art level performance. We
additionally show that Dagma-DCE allows for principled thresholding and
sparsity penalties by domain-experts. The code for our method is available
open-source at https://github.com/DanWaxman/DAGMA-DCE, and can easily be
adapted to arbitrary differentiable models.

Causal discovery, interpretable machine learning, differential causal effect

##  1 INTRODUCTION

The discovery of causal relationships is the subject of many scientific
questions, in medicine, economics, the social sciences, and more, and is
therefore a question of utmost importance in science and engineering. Given an
observational signal, one may therefore be interested in learning the causal
structure underlying the data. Under Pearl’s formalism of causality [1], the
causal structure is captured by a _causal graph_ , which is a directed acyclic
graph (DAG) with edges denoting direct causation of one variable to another.

A core principle of causal inference is the notion of _intervention_ , where
one can interact with the system that produces signals. This is in contrast to
causal inference from _observational_ data, where we can only observe the
outputs of a system, as is typically the case in signal processing or machine
learning. Unfortunately, causal discovery using observational data is
generally ill-posed and requires strong assumptions about the system we
observe [2, pp. 135]. Different methods may use different assumptions, where
the validity of the assumptions may be extremely domain- or problem-specific.
Even when assumptions are made to ensure (at least partial) identifiability of
the causal graph, causal discovery can be extremely expensive, in terms of
both the computation and the amount of required data.

One major class of methods, referred to as _constraint-based_ methods, is
based on the observation that a causal graph encodes a set of conditional
independencies (CIs). By assuming that all CIs are encoded in the graph, an
assumption known as _causal faithfulness_ , constraint-based methods utilize
statistical CI tests to determine which causal graphs are consistent with the
data.

While some of the earliest successful causal discovery algorithms are
constrained-based, such as the PC Algorithm and Fast Causal Inference [3],
they can be very costly. Namely, they are both expensive computationally,
scaling exponentially with the number of variables in the worst case, and they
are expensive in terms of data. Additionally, CI tests tend not to be very
powerful, and constrained-based causal discovery methods require conditional
data for many different values and combinations of random variables. Moreover,
causal faithfulness is not a testable assumption in general, though some
weaker versions are [4].

An alternative approach to causal discovery using observational data is to use
_score-based_ methods. These methods attempt to find a causal model
$\mathcal{M}$ and its associated causal graph $\mathcal{G}$, which maximizes a
score function $\mathcal{S}(\mathcal{M},\mathcal{D})$ for observational data
$\mathcal{D}=\\{{\boldsymbol{\mathbf{x}}}_{1},{\boldsymbol{\mathbf{x}}}_{2},\dots,{\boldsymbol{\mathbf{x}}}_{N}\\}$
[2, pp. 148]. For example, one can search for the causal model which maximizes
the likelihood or minimizes the mean square error (MSE). In practice, one
typically optimizes a version of these objectives which is penalized for model
complexity, such as the Bayesian information criterion. For certain classes of
causal models, this approach is simple and effective — indeed, early attempts
at score-based causal discovery were made contemporaneously to early
constraint-based approaches [5]. Despite their simplicity, score-based methods
suffer from the super-exponential search space of DAGs with $d$ variables, as
well as identifiability issues [2, pp. 149].

Fundamentally troubling to score-based methods is the difficulty of
optimization in discrete spaces, which tend to be combinatorially more
difficult than continuous objectives. To address this, Notears (short for _N
on-combinatorial Optimization via Trace Exponential and Augmented lagRangian
for Structure learning_) [6] proposes a novel algebraic constraint on
adjacency matrices, such that the constraint is satisfied if and only if the
adjacency matrix $\mathbf{A}$ represents a DAG. Notears then performs
constrained continuous optimization on a penalized MSE objective to learn a
causal graph without explicitly exploring the discrete space of all DAGs.

Notears was originally defined for linear models, with the adjacency matrix
entries being coefficients of the corresponding linear functions. Notears+ [7]
extends Notears to nonlinear functions and achieves state-of-the-art
performance on simulated data. In lieu of linear coefficients, Notears+ uses
opaque proxies to define $\mathbf{A}$, providing advantages in optimization
but at the cost of interpretability and generality of the results.

Several methods have provided alternative acyclicity constraints which are
faster to compute or are otherwise preferable in optimization. One striking
example is Dagma (short for _DAG s via $\mathbf{M}$-matrices with Acyclicity_)
[8], which in practice can optimize up to an order of magnitude faster than
Notears.

Our contribution is to provide a novel learning objective that results in
obtaining interpretable weighted causal graphs. Our method has several
important strengths. It is model-agnostic in nature, with the ability to use
any twice-differentiable model trainable with gradient descent. Furthermore,
the sparsity penalty used in the scoring function $\mathcal{S}(\cdot,\cdot)$
is also model-agnostic. Lastly, we can interpret the causal graphs with the
differentiable causal effect (DCE) metric [9]. We base our optimization on the
method of Dagma, leading us to name our method Dagma-DCE.

The rest of this paper is structured as follows: in Section II, we review
score-based causal discovery and the current state-of-the-art in
differentiable causal discovery. In Section III, we show that the adjacency
matrix returned by Notears and Dagma is not interpretable, first by a
theoretical argument and then by a toy example. We propose a new graph
definition and optimization problem in Section IV, forming Dagma-DCE, and
discuss its interpretability in Section V. Finally, we run experiments on
synthetic data in Section VI showing that Dagma-DCE performs competitively
with other methods and returns legitimately different adjacency matrices than
Dagma. Concluding remarks are provided in Section VII.

##  2 BACKGROUND

We begin by reviewing some concepts of causal inference, namely structural
causal models (SCMs), which formalize causality. We then introduce
differentiable causal discovery methods, which utilize continuous optimization
for score-based causal discovery.

###  2.1 Structural Causal Models and Causal Discovery

SCMs provide a formalism for causality, where causal relationships are
described by a DAG and accompanying functional assignments. To make this more
clear, let ${\boldsymbol{\mathbf{x}}}=[x_{1},\dots,x_{d}]$ be $d$-dimensional
random vector, and $\mathcal{G}$ be a DAG whose vertices are the components of
${\boldsymbol{\mathbf{x}}}$. Then, according to $\mathcal{G}$, an SCM
describes each component $x_{j}$ as a function of its parents
$\operatorname{Pa}(x_{j})$ and independent noise $z_{j}$,
$f_{j}(\operatorname{Pa}(x_{j}),z_{j})$, i.e.,

$x_{j}\coloneqq f_{j}\left(\operatorname{Pa}(x_{j}),z_{j}\right),\quad
j=1,\dots,d.$ (1)

We use the notation $\coloneqq$ to reinforce the somewhat philosophical point
that these are _assignments_ , and not merely equivalence [2]. We will
frequently refer to the SCM of a system as simply the “causal model”
$\mathcal{M}$.

The goal of structure learning is then to identify the SCM of a system, given
data $\mathcal{D}$. If one cannot interact with the system, then we must make
assumptions on the form of each $f_{j}$ and the distribution of the noise
$z_{j}$ to uniquely determine the SCM [2, pp. 139]. These assumptions can be
parametric, for example, assuming each $f_{j}$ is linear and each $z_{j}$ is
non-Gaussian [10], or they can be non-parametric, for example, assuming the
$f_{j}$ are strictly nonlinear with the $z_{j}$ being additive and Gaussian
[11, Corollary 31]. Once we make such a restriction to our model, we call it
_identifiable_.

After establishing identifiability, there are several families of techniques
for learning the SCM from observational data. Here, we focus on score-based
methods, though the interested reader is pointed to [12] for a recent,
comprehensive review. In score-based causal discovery, we find a causal model
$\mathcal{M}$ which maximizes a score function
$\mathcal{S}(\mathcal{M},\mathcal{D})$. In practice, we parameterize models by
some vector ${{\boldsymbol{\mathbf{\theta}}}}$ and then find the model
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$ which optimizes
$\mathcal{S}(\cdot,\cdot)$:

$\hat{{{\boldsymbol{\mathbf{\theta}}}}}=\arg\max_{{\boldsymbol{\mathbf{\theta}}}}\mathcal{S}(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}},\mathcal{D})\textrm{
s.t. }\mathcal{G}_{{\boldsymbol{\mathbf{\theta}}}}\text{ is a DAG.}$ (2)

The main drawback of naïve score-based methods is that they have to search
over the space of DAGs with $d$ nodes, which grows extremely quickly. For
$d=5$, this number is already $29281$; for $d=14$, it is already greater than
$10^{36}$ [13, A003024]. Therefore, for even small graphs, one must resort to
greedy or approximate techniques [2, pp. 150].

###  2.2 Differentiable Causal Discovery

Differentiable causal discovery methods are score-based methods that recast
Eq. 2 as a constrained continuous optimization problem; as before, let
$\mathcal{S}(\cdot,\cdot)$ be a score function that reflects how well a model
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$ explains the data
$\mathcal{D}$, and
$\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\in\mathbb{R}_{\geq 0}^{d\times
d}$ be a weighted adjacency matrix derived from
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$. Then, using a function
$h(\cdot)$, which is $0$ if and only if
$\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}$ represents a DAG, the
optimization problem Eq. 2 becomes

$\displaystyle\max_{{{\boldsymbol{\mathbf{\theta}}}}\in{{\boldsymbol{\mathbf{\theta}}}}}$
$\displaystyle\mathcal{S}(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}},\mathcal{D}),$
(3) s.t. $\displaystyle h(\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}})=0.$

The earliest differentiable causal discovery algorithm with a DAG constraint
is Notears [6], which uses a trace of a nonlinear function on the adjacency
matrix $\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}$, i.e.,

$h_{\textsc{Notears}}(\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}})=\operatorname{tr}\left(\exp(\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\odot\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}})\right)-d,$

where $\odot$ is the Hadamard (elementwise) product. The motivation for this
constraint is that powers of the adjacency matrix
$\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}^{k}$ reflect how many closed
walks of length $k$ are present in the corresponding graph
$\mathcal{G}_{{\boldsymbol{\mathbf{\theta}}}}$. Since DAGs have no closed
walks of any length, $h(\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}})$ is $0$
precisely when $\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}$ is a DAG.

Quickly, several simplifications on the same constraint were derived. For
example, one need not calculate the entire matrix exponential, with a
polynomial sufficing instead [14].

The alternative constraint that we use is Dagma [8]; with $s>0$, the Dagma
constraint function is

$h_{\textsc{Dagma}}(\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}})=-\log\det\left(s\mathbb{1}_{d}-\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\odot\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\right)+d\log
s,$ (4)

where $\mathbb{1}_{d}$ is the $d\times d$ identity matrix. This constraint is
motivated by the fact that $\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}$
describes a DAG if and only if it is nilpotent. While the constraint given by
Eq. 4 is, in general, a relaxation of nilpotency, the authors of Dagma
characterize the set of matrices for which zeroness of Eq. 4 is equivalent to
acyclicity, called $\mathbf{M}$-matrices. They then show that the set of
$\mathbf{M}$-matrices includes all DAGs, and that the gradient of
$h_{\textsc{Dagma}}(\cdot)$ is sufficiently well-behaved, to enable
constrained optimization within $\mathbf{M}$-matrices.

It remains to specify the model
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$, the score function
$\mathcal{S}$, and how to derive
$\mathcal{G}_{{\boldsymbol{\mathbf{\theta}}}}$ from
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$. The approach introduced in
Notears+ and used in Dagma111In general, Dagma refers to both linear and
nonlinear models using Eq. 4. To avoid confusion, Dagma exclusively refers to
the nonlinear case here. is to define

${{\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}}_{ij}}\triangleq\lVert\partial_{i}f_{j}\rVert_{L^{2}}=\sqrt{\int_{\mathbb{R}^{d}}\left(\partial_{i}f_{j}({\boldsymbol{\mathbf{x}}})\right)^{2}\,d{\boldsymbol{\mathbf{x}}}},$
(5)

where $f_{j}$ is the function ${\boldsymbol{\mathbf{x}}}\mapsto x_{j}$
according to $\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$. To guarantee Eq.
5 is well-defined, it is required that $\lVert\partial_{i}f_{j}\rVert_{L^{2}}$
is finite. Accordingly, the Notears+ authors restrict
$\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$ to a space of functions called
the Sobolev space $H^{1}(\mathbb{R}^{d})$.

However, the Notears+ authors argue that for causal discovery, the zeroness of
${{\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}}_{ij}}$ is what matters, and so
a condition equivalent to
${{\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}}_{ij}}=0$ is derived for
several classes of models. For example, in a multilayer perception (MLP)
network, the $L^{2}$ norm of the first layer in the neural network is chosen.
The score function is defined to be a sparsity-penalized MSE222Dagma uses the
log-likelihood instead for benchmarks, but we use the MSE to simplify
comparisons, and avoid assumptions on the data-generating process.,

$\mathcal{S}(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}},\mathcal{D})=-(\operatorname{MSE}(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}},\mathcal{D})+\lambda_{1}g(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}})),$
(6)

where $\lambda_{1}$ is a hyperparameter and
$g(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}})$ defines a sparsity penalty.
For MLPs, Notears+ and Dagma use the $\ell_{1}$ norm of the first layer MLP
weights.

In practice, it is very difficult to exactly follow the constraint in Eq. 3.
For this reason, as well as to avoid false positives, an additional
thresholding step is performed, where $\mathbf{A}_{ij}$ is set to $0$ if it is
less than some predetermined hyperparameter $C$.

##  3 NON-INTERPRETABILITY OF Dagma

A key motivation for Dagma-DCE is that existing methods utilize a non-
interpretable weighted adjacency matrix
$\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}$. In this section, we argue
theoretically that the adjacency matrix of Dagma is not interpretable and
subsequently show empirically that Dagma fails to return intuitively-
interpretable results on linear SCMs.

For the purposes of this section, we restrict ourselves to models where
$f_{j}$ is modeled with an MLP network, though our arguments can be adjusted
for other cases. The theoretical section can safely be skipped, but shows
explicitly that adjacency matrices derived from the first layer of a neural
network may be quite different than those defined by derivatives.

###  3.1 Theoretical Argument

For the purposes of this section, let $f_{j}\colon\mathbb{R}^{d}\to\mathbb{R}$
be an MLP with weight matrices $\mathbf{A}^{(1)},\dots,\mathbf{A}^{(M)}$ and
activation function $\sigma(\cdot)$. Then the $(M-1)$-hidden layer MLP $f_{j}$
can be expressed as

$f_{j}({\boldsymbol{\mathbf{x}}})=\mathbf{A}^{(M)}\sigma\left(\mathbf{A}^{(M-1)}\sigma\left(\cdots\sigma\left(\mathbf{A}^{(1)}{\boldsymbol{\mathbf{x}}}\right)\right)\right).$
(7)

Despite the exact analytic result that

$\lVert\partial_{i}f_{j}\rVert_{L^{2}}=0\iff\lVert\mathbf{A}^{(1)}_{i\cdot}\rVert_{L^{2}}=0,$
(8)

differentiable causal discovery methods rely on thresholding small, nonzero
values. The following lemma shows that when using sigmoid activation
functions, as is done in Notears+ and Dagma, the two norms in Eq. 8 may be
arbitrarily different.

######  Lemma 1.

Let $\sigma(\cdot)$ denote the sigmoid activation function. Then for any
$\delta,\epsilon>0$, there exists an MLP $f_{j}$ with weight matrices
$\mathbf{A}^{(1)},\dots,\mathbf{A}^{(M)}$ such that
$\lVert\mathbf{A}^{(1)}_{i\cdot}\rVert_{L^{2}}<\epsilon$ but
$\lVert\partial_{i}f_{j}\rVert_{L^{2}}>\delta$.

As (7) is a repeated composition of smooth functions, it suffices to show the
result for the $1$-hidden layer MLP

$g({\boldsymbol{\mathbf{x}}})=\mathbf{A}^{(2)}\sigma\left(\mathbf{A}^{(1)}{\boldsymbol{\mathbf{x}}}\right).$

If $\mathbf{A}^{(2)}$ is a $K\times H$ matrix, then

$[g(x)]_{k}=\sum_{h}^{H}\mathbf{A}^{(2)}_{kh}\sigma\left(\sum_{j=1}^{d}\mathbf{A}^{(1)}_{hj}x_{j}\right).$

For notational clarity, let

$z_{h}=\sigma\left(\sum_{j=1}^{d}\mathbf{A}^{(1)}_{hj}x_{j}\right).$

Additionally, note that $\sigma^{\prime}(\cdot)>0$. Then so long as
$\mathbf{A}^{(1)}_{hi}$ is nonzero for some $h$, we have by chain rule that

$\left\lvert\frac{\partial z_{h}}{\partial
x_{i}}\right\rvert>\epsilon^{\prime}$

for some $\epsilon^{\prime}>0$. Then by directly calculating the norm,

$\displaystyle\lVert\partial_{i}f_{j}\rVert_{L^{2}}^{2}$
$\displaystyle=\int\sum_{k=1}^{K}\left(\frac{\partial g_{k}}{\partial
x_{i}}\right)^{2}\,dx_{i}$
$\displaystyle=\int\sum_{k=1}^{K}\left(\sum_{h=1}^{H}\mathbf{A}^{(2)}_{kh}\frac{\partial
z_{h}}{\partial x_{i}}\right)^{2}\,dx_{i}$
$\displaystyle\geq\int\sum_{k=1}^{K}\left(\sum_{h=1}^{H}\mathbf{A}^{(2)}_{kh}\epsilon^{\prime}\right)^{2}\,dx_{i}.$

we conclude that it is unbounded. Therefore, we can choose $\mathbf{A}^{(2)}$
such that the function norm is arbitrarily large, and in particular, larger
than $\delta$.

While this result may not be particularly “generic,” in the sense that the
resulting function $f_{j}$ may be unlikely to be learned in practice, it does
show that $\partial_{i}f_{j}$ can be arbitrarily different than any nonzero
$\mathbf{A}^{(1)}$. This situation only gets worse for the more popular
rectified linear unit (ReLU) activation function, where
$\lVert\partial_{i}f_{j}\rVert_{L^{2}}$ is arbitrary for a _fixed_ $f_{j}$.
This is formalized in the following lemma:

######  Lemma 2.

Let $\sigma(\cdot)$ denote the ReLU activation function. Then for any $s>0$
and $\mathbf{A}^{(1)},\dots,\mathbf{A}^{(M)}$, there exist matrices
$\mathbf{B}^{(1)},\dots,\mathbf{B}^{(M)}$ such that

$\displaystyle f_{j}({\boldsymbol{\mathbf{x}}})$
$\displaystyle=\mathbf{A}^{(M)}\sigma\left(\mathbf{A}^{(M-1)}\sigma\left(\cdots\sigma\left(\mathbf{A}^{(1)}{\boldsymbol{\mathbf{x}}}\right)\right)\right)$
$\displaystyle=\mathbf{B}^{(M)}\sigma\left(\mathbf{B}^{(M-1)}\sigma\left(\cdots\sigma\left(\mathbf{B}^{(1)}{\boldsymbol{\mathbf{x}}}\right)\right)\right),$

and $\lVert\mathbf{B}^{(1)}_{i\cdot}\rVert_{L^{2}}=s$.

Note that for any real $z$, the ReLU function is invariant to positive
rescaling: $\sigma(z)=\sigma(az)/a$ for any $a>0$. This implies, more
generally, that
$\mathbf{S}^{-1}\mathbf{A}^{(2)}\sigma(\mathbf{S}\mathbf{A}^{(1)}{\boldsymbol{\mathbf{x}}})=\mathbf{A}^{(2)}\sigma(\mathbf{A}^{(1)}{\boldsymbol{\mathbf{x}}})$,
where $\mathbf{S}$ is a diagonal matrix [15]. The lemma then immediately
follows by setting $\mathbf{B}^{(1)}=\mathbf{S}\mathbf{A}^{(1)}$ and
$\mathbf{B}^{(2)}=\mathbf{S}^{-1}\mathbf{A}^{(2)}$, with the appropriate
choice of $\mathbf{S}$.

###  3.2 Empirical Result with a Linear SCM

(a) Dagma

(b) Dagma-DCE

Figure 1: The difference between the magnitude of the true derivatives in a
linear causal model to the magnitude of the weighted graph in Dagma for a
random $10\times 10$ Erdös-Rényi directed graph with $20$ expected edges. Gray
boxes surrounding each cell denote the magnitude of the ground-truth linear
coefficient.

In the lemmas above, we have shown that the adjacency matrix of Notears+ and
Dagma may be arbitrarily different than the $L^{2}$ function norm. Here, we
illustrate this fact in practice by applying Dagma to a linear structural
equation model (SEM), and showing that the returned adjacency matrix does not
correspond to the intuitive notion of causal strength.

For this purpose, we simulated an Erdös-Rényi (ER) graph with $10$ random
variables and an expectation of $20$ edges, creating a DAG as in [6]. Given a
DAG, a SEM assigns functional relationships between a variable and its
parents. In this case, a linear SEM is assigned by randomly sampling linear
coefficients from $\mathcal{U}\left((-2.0,-0.5)\cup(0.5,2.0)\right)$.

The use of a linear SEM allows for an obvious ground-truth weighted graph to
compare to, namely the linear coefficients. We then ran Dagma with each
$f_{j}$ being modeled as an MLP with one hidden layer of $10$ units and
compared the resulting weighted adjacency matrix to the ground truth in Fig.
1(a). For the purposes of this demonstration, we tested with
$\lambda_{1}\in\\{0,10^{-4},10^{-3},10^{-2},10^{-1}\\}$ and chose the value
which minimized the Frobenius norm
$\lVert\mathbf{A}_{\text{est}}-\mathbf{A}_{\text{true}}\rVert_{F}$ to give a
favorable comparison. This resulted in using $\lambda_{1}=10^{-2}$, with all
other hyperparameters set to their default values.

Two things are visually apparent: first, Dagma learns the graph topology quite
well, including one erroneous edge in an otherwise perfect prediction. Second,
however, is that the weights recovered by Dagma do not obviously correspond to
the interpretable, linear coefficients, i.e., it cannot easily be interpreted
as a “strength of causation.” Even worse, it is not easy to predict when Dagma
is over- or under-estimating the strength of an edge. While we defer a
definition of Dagma-DCE and discussion of Fig. 1(b) to the subsequent
sections, a view of Fig. 1(b) shows a much closer resemblance to linear
coefficients. It is important to note that this does not indicate “errors” in
Dagma, but instead that the magnitude of the weighted graph _is not defined_
to be the true “causal strength” — their very definition has lost meaning
compared to the linear case.

##  4 PROPOSED SOLUTION

In this section, we provide a definition of the learning objective for Dagma-
DCE, and subsequently discuss practical details of implementing such a scheme.
In particular, we discuss how Dagma-DCE is model-agnostic in both formulation
and implementation, and we also explore simple methods for accelerating and
pre-training models.

###  4.1 Definition of the Optimization Problem

An initial solution would be to estimate Eq. 5 directly; however, we argue
$H^{1}(\mathbb{R}^{d})$ is the wrong space to define models over. Not only are
many interesting functions outside of $H^{1}(\mathbb{R}^{d})$, but also simple
functions, for example linear functions, do not have a finite Sobolev norm
over all of $\mathbb{R}^{d}$. Additionally, from a data-centric perspective,
there is little reason to care about the behavior of the model where there is
a vanishingly small probability mass.

Instead, we choose to define the adjacency matrix as the $L^{2}$ norm of the
derivative $\partial_{i}f_{j}$ with respect to the probability distribution of
${\boldsymbol{\mathbf{X}}}$, which we denote
$\mathbb{P}^{\boldsymbol{\mathbf{X}}}$. While the theory of Sobolev spaces
plays no further role in our work, to adapt the language of [7], models are
defined over the _weighted Sobolev space
$H^{1}(\mathbb{R}^{d},\mathbb{P}^{\boldsymbol{\mathbf{X}}})$_ [16]. This is
simply to say that the models’ function and first (weak) derivatives are
square-integrable with respect to $\mathbb{P}^{\boldsymbol{\mathbf{X}}}$. We
then offer the alternative definition and Monte Carlo approximation

$\displaystyle\mathbf{A}_{ij}\triangleq\lVert\partial_{i}f_{j}\rVert_{L^{2}(\mathbb{P}^{\boldsymbol{\mathbf{X}}})}$
$\displaystyle=\sqrt{\int_{\mathbb{R}^{d}}\left(\partial_{i}f_{j}({\boldsymbol{\mathbf{x}}})\right)^{2}p_{{\boldsymbol{\mathbf{X}}}}({\boldsymbol{\mathbf{x}}})\,d{\boldsymbol{\mathbf{x}}}}$
$\displaystyle\approx\sqrt{\frac{1}{N}\sum_{n=1}^{N}\left(\partial_{i}f_{j}({\boldsymbol{\mathbf{x}}}_{n})\right)^{2}},$
(9)

where $p_{\boldsymbol{\mathbf{X}}}({\boldsymbol{\mathbf{x}}})$ is the density
of $\mathbb{P}^{\boldsymbol{\mathbf{X}}}$ with respect to the Lebesgue
measure, and ${\boldsymbol{\mathbf{x}}}_{n}\sim
p_{\boldsymbol{\mathbf{X}}}({\boldsymbol{\mathbf{x}}})$. By noting that the
entries of $\mathbf{A}$ are the root-mean-square entries of the Jacobian
matrix, it is then straightforward to calculate $\mathbf{A}_{ij}$ for any
differentiable model, with the above being a Monte Carlo approximation over
empirical Jacobian matrices.

This definition of $\mathbf{A}$ also allows for a principled sparsity penalty
— indeed, if $\partial_{i}f_{j}$ is square-integrable with respect to
$\mathbb{P}^{\boldsymbol{\mathbf{X}}}$, then Jensen’s inequality implies that
it is also integrable, which allows us to define and approximate

$g(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}})=\lVert\partial_{i}f_{j}\rVert_{L^{1}(\mathbb{P}^{\boldsymbol{\mathbf{X}}})}\approx\frac{1}{N}\sum_{n=1}^{N}\lvert\partial_{i}f_{j}({\boldsymbol{\mathbf{x}}}_{n})\rvert.$
(10)

Combining Eqs. 3, 4, 6, 4.1 and 10, we arrive at the optimization problem

$\displaystyle\min_{{{\boldsymbol{\mathbf{\theta}}}}\in{{\boldsymbol{\mathbf{\theta}}}}}$
$\displaystyle\operatorname{MSE}(\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}},\mathcal{D})+\lambda_{1}\lVert\partial_{i}f_{j}\rVert_{L^{1}(\mathbb{P}^{\boldsymbol{\mathbf{X}}})}$
(11) s.t.
$\displaystyle-\log\det\left(s\mathbb{1}_{d}-\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\odot\mathbf{A}_{{\boldsymbol{\mathbf{\theta}}}}\right)+d\log
s=0.$

Similar to Dagma, we use a central path optimization technique to solve the
constrained optimization problem.

Our approach has the same convergence properties as Notears and Dagma, using
similar assumptions. Namely, we assume causal sufficiency, i.e., a lack of
latent confounders, and identifiability, which ensures that the optima of the
score are the true causal graph. There are many options for ensuring
identifiability, for example, thrice-differentiable strictly-nonlinear
additive noise models suffice [11]. Additionally, while it is not necessary
for the optimization problem to be well-defined, the interpretability of
Dagma-DCE relies on once-differentiability of the underlying $f_{j}$. Namely,
this excludes discrete or mixed-type data, since how a model interpolates
functions between discrete labels is arbitrary.

###  4.2 Dagma-DCE is Model Agnostic

The derivatives in Section 4.1 may be available analytically for many models.
For the simplest case, take $\mathcal{M}_{{\boldsymbol{\mathbf{\theta}}}}$ to
be linear, i.e.,

$x_{j}\coloneqq\sum_{i=1}^{d}\beta_{ij}x_{i}.$

Then the linear coefficients $\beta_{ij}$ are the derivatives used in Section
4.1. Note that Eq. 11 then becomes the problem of linear Dagma, meaning Dagma-
DCE is a legitimate generalization.

Simple (semi-)analytic expressions are also available for additive models:

$x_{j}\coloneqq\sum_{k=1}^{K}\beta_{jk}g_{k}({\boldsymbol{\mathbf{x}}}),$ (12)

in which case

$\partial_{i}f_{j}=\sum_{k=1}^{K}\beta_{jk}\frac{\partial g_{k}}{\partial
x_{i}}.$

Analytic expressions of coefficients are even available for several classes of
non-parametric models. For example, by the representer theorem [17, pp. 132],
a Gaussian process (GP) posterior can be written in the form of Eq. 12, where
$g({\boldsymbol{\mathbf{x}}})$ depends on the kernel $\kappa(\cdot,\cdot)$ and
$\beta_{j1},\dots,\beta_{jn}$ are easily calculated:

$\begin{split}f_{j}&\sim\mathcal{GP}({\boldsymbol{\mathbf{m}}},\kappa);\\\
x_{j}&\coloneqq
f_{j}({\boldsymbol{\mathbf{x}}}),\end{split}\;\Longleftrightarrow\;\begin{split}x_{j}\coloneqq\sum_{n=1}^{N}\beta_{n}\kappa({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}_{n})\end{split}.$
(13)

Closed-form expressions for several common kernels, including the automatic
relevance detection squared exponential (ARD-SE) and Matérn-3/2 kernels, are
derived in [9].

Furthermore, even when analytic expressions become cumbersome to derive by
hand, the derivative can often be obtained via automatic differentiation (AD)
tools. Among many other differentiable models, a notable case is with MLPs. It
is also the case for extensions of Notears+-like objectives to time series,
for example the one-dimensional convolutional neural networks used in NTS-
Notears [18].

###  4.3 Practical Considerations

Several general comments are made in Section 4.1 of [8] which also hold for
Dagma-DCE, particularly concerning the optimization details. For example, the
hyperparameter $s$ can be changed freely, even between optimization
iterations. This can affect performance, as $s$ controls the volume of the
feasible region. The authors of [8] recommend starting with $s=1$ and slowly
reducing its value. Additionally, different optimizers can be used to solve
sub-problems along the central path; Dagma uses Adam [19], which we also use.
Finally, the use of central path methods also creates additional
hyperparameters, including central path coefficients and weight
initializations, which are discussed.

While backward passes over the Jacobian is indeed expensive and a drawback of
the Dagma-DCE, the Monte Carlo approximation can be effectively parallelized
with respect to $N$, particularly on GPUs. Therefore, while scaling with $d$
is inherently difficult, scaling with $N$ is quite inexpensive.

In practice, we additionally see benefits from pre-training for Dagma-DCE,
namely by training a single iteration of Dagma. Intuitively, while Dagma and
Dagma-DCE are fundamentally different optimization problems, their underlying
task is the same. Therefore, a single iteration of Dagma serves as a good
starting point, and saves costly backward passes over the Jacobian.

Finally, Dagma-DCE has several hyperparameters, including the strength of
$L_{1}$ regularization $\lambda_{1}$ and thresholds. In the next section, we
discuss how the formulation of Dagma-DCE makes the threshold hyperparameter
interpretable. However, the issue of setting $\lambda_{1}$ remains. One
option, explored in [7] and [18], is to optimize $\lambda_{1}$ with respect to
some metric using cross-validation if labeled data/causal graph pairs exist.
How to choose $\lambda_{1}$ when these labeled pairs are not available is an
interesting task for future research.

##  5 INTERPRETABILITY OF Dagma-DCE

The interpretability of the Dagma-DCE approach stems from the observation that
the partial derivatives, $\partial_{i}f_{j}$, can be interpreted through the
lens of causal strength. We begin this section with general notes on causal
strengths, followed by the differential causal effect, which forms the basis
of interpretability for Dagma-DCE. We revisit the empirical example of Section
III-3.2 and show that Dagma-DCE recovers the intuitive result. Finally, we
discuss how thresholding in Dagma-DCE is a fundamentally interpretable
operation, allowing the incorporation of prior knowledge from experts and
decision-makers.

###  5.1 Measuring Causal Strength

The quantification of causal strength is a highly nontrivial task, with
several competing definitions [20]. Despite this, we find it helpful to ground
our reasoning in the linear case, where an unambiguous notion of causal
strength is available. In particular, the linear model coefficients quantify
exactly how a change in one variable will induce changes in its children.
Consider a linear SCM with coefficients given by a matrix
$[\mathbf{B}]_{ij}=\beta_{ij}$ and additive noise $z_{j}$:

$x_{j}:=\sum_{i=1}^{d}\beta_{ij}x_{i}+z_{j},\qquad j=1,...,d.$

Then the entry $\beta_{ij}$ represents the direct influence of $x_{i}$ to
$x_{j}$. Note that ${\mathbf{B}}$ can be viewed as a causal graph, provided
that it is acyclic, or equivalently, nilpotent.

Occasionally, these direct coefficients summarize all causal information that
$x_{i}$ holds about $x_{j}$. Indeed, if $\mathbf{B}$ is viewed as a graph,
this is the case if the only path $x_{i}$ and $x_{j}$ is the direct one, in
which case $\mathbf{B}_{ij}$ is the total causal strength. Otherwise, in
general, we must sum across all paths in the graph, and multiply the
coefficients along each path. This can be expressed compactly by a geometric
sum,

$\mathbf{T}=\sum_{\ell=1}^{d}\mathbf{B}^{\ell}.$

Causal sufficiency of the system then tells us that the entries
$\mathbf{T}_{ij}$ of the matrix $\mathbf{T}$ are sufficient to describe how
any change to $x_{i}$ will affect $x_{j}$, and thus, $\mathbf{T}$ can be
useful to ask counterfactual questions or to understand the impact of various
causes to a single effect.

###  5.2 Differential Causal Effect and Dagma-DCE

The DCE provides a nonlinear version of the linear model causal strength by
considering the local behavior of the model [9]. The Jacobian matrix,
$\mathcal{J}_{f}$, replaces the role of ${\mathbf{B}}$, describing how changes
to a parent affect the child node when there are no other paths in the system.
Similarly, the geometric sum of the Jacobian also defines the total causal
effect from source to effect, for small changes.

Relating to Dagma-DCE, the adjacency matrix defined in (4.1) can be
interpreted as the root-mean-squared DCE — as a result, $\mathbf{A}_{ij}\neq
0$ if and only if the DCE, $\partial_{i}f_{j}$, is not identically zero. The
magnitude of $\mathbf{A}_{ij}$ quantifies how strong the interaction is, in
terms of the energy of the DCE functions.

###  5.3 Empirical Results with a Linear SCM

We quickly revisit the experiment of Section III-3.2, where Dagma was fit to a
linear SCM, and its weighted adjacency matrix compared to the ground truth.
The setup is nearly identical, with the same parameterization of each function
$f_{j}$. The main difference is that we do not tune the sparsity penalty
hyperparameter, instead we simply set $\lambda_{1}=0$. Like Dagma, the Dagma-
DCE method obtains an accurate estimate of the graph topology, recovering the
underlying graph perfectly. Unlike Dagma, we additionally recover the
“intuitive” notion of causal strength.

###  5.4 Thresholding in Dagma-DCE

An ad-hoc component of Notears+-like methodologies is the thresholding step.
Indeed, across different datasets, the thresholding hyperparameter can have
large consequences to accuracy, as shown in [7]. However, many conceivable
causal inference tasks will not have any validation set to calibrate the
thresholding parameter with, meaning the user will have to make a heuristic
estimate.

On the contrary, Dagma-DCE’s weighted adjacency matrix entries are
interpretable, so the choice of a threshold carries direct meaning to the
modeler or expert. This decision can even be made on an edge-by-edge basis, or
across (a possibly rescaled version of) $\mathbf{A}$.

For example, one may decide to normalize the entry $A_{ij}$ of by the marginal
variance of $x_{j}$ for the purposes of thresholding, reminiscent of ANOVA
[21]. An alternative scheme could be normalizing by the sum of the $j$th
column; while subadditivity of norms means this does not actually scale by the
total derivative $\partial f_{j}$, it can be viewed as scaling by the total
causal effect. Finally, since the entries of $\mathbf{A}$ are now
unconstrained, one could attempt thresholding methods developed for Granger
causality, for example via clustering [22] or denoising [23]. How this
thresholding step should be interpreted and executed forms an interesting
question for future work.

(a)

(b)

Figure 2: Resulting SID (top left), SHD (top right), $F_{1}$ Score (bottom
left), and time elapsed (bottom right) for random data generated from (2(a))
the ER-4 GP-additive model and (2(b)) the ER-4 MLP model, as detailed in
Section 6.2. Boxes show the median and quartiles across $T=10$ trials for
Dagma and Dagma-DCE, and $T=5$ trials for Notears, with whiskers showing the
minimum and maximum values.

##  6 EXPERIMENTS

We perform several experiments similar to those conducted by the authors of
Dagma [8], providing comparisons to Dagma and Notears+. We begin with a brief
discussion of metrics used to evaluate results, followed by experimental
details, results, and discussion.

###  6.1 Metrics

Evaluating results in causal discovery is complex in the sense that several
different metrics exist with different goals and calculations. The two most
basic metrics we report are the structural Hamming distance (SHD), which
calculates the number of incorrect edges a graph contains, and the $F_{1}$
score, which is the geometric mean of precision and recall.

Additionally, we include the structural intervention distance (SID) [24],
which is arguably more well-suited to causal inference. The SID is motivated
by the idea that, when calculating causal effects via “adjustment sets,” some
errors are more severe than others. For example, including spurious relations
does not affect causal effects, while reversing the order of an edge does. The
SID calculates the number of edges $i\to j$ such that the adjustment set
corresponding to the calculation of the causal effect of $x_{i}$ on $x_{j}$ is
not valid.

Finally, we report the time elapsed in minutes in each of our synthetic
benchmarks.

###  6.2 Experimental Details

We replicate a common synthetic dataset in the causal discovery literature,
generating ER-$4$ DAGs according to the procedure of [7], with functional
relationships assigned in one of two ways. The first way is using additive
GPs:

$f_{j}({\boldsymbol{\mathbf{x}}})=\sum_{i\in\operatorname{Pa}(x_{j})}g_{ij}(x_{i}),$

where each $g_{ij}$ is sampled from an RBF GP with lengthscale $1$. The other
is with a random MLP network with hidden size $100$ and a sigmoid activation.
Crucially, all weights are samples from
$\mathcal{U}\left((-2.0,-0.5)\cup(0.5,2.0)\right)$.

As noted in Section IV-4.3, the hyperparameters of each method may generally
be tuned with cross-validation. Nonetheless, to facilitate a fair comparison
between methods, we use a set of “defaults” as in [7, 8, 18]. For Dagma and
Notears+, we use the hyperparameters noted in Appendix C.2.2 of [8]. However,
since the sparsity penalty of Dagma-DCE is fundamentally different than that
of other methods, it is difficult to find a “default” setting for a fair
comparison. We chose $\lambda_{1}=3.5\times 10^{-2}$, which provided a
sparsity penalty of similar magnitude in the MLP experiments. This is
imperfect, but allows for the initial score functions of Dagma and Dagma-DCE
to be similar without unfairly optimizing $\lambda_{1}$. Similarly, we used
$0.25$ as thresholds for comparison, which provides a similar true positive
rate in the MLP experiments. Additionally, due to the added cost of Dagma-DCE,
we set the maximum number of iterations to $10\%$ of the corresponding values
in Dagma.

All experiments were run on a NVIDIA Titan RTX GPU. The authors’
implementations of Notears+333https://github.com/xunzheng/notears and
Dagma444https://github.com/kevinsbello/dagma were used, with slight
modifications to allow GPU acceleration. All methods (including Dagma-DCE)
used the same MLP architecture used in the experiments of [7, 8].

###  6.3 Results & Discussion

Our results show that Dagma-DCE achieves the state-of-the-art on common
synthetic datasets: in additive GPs (Fig. 2(a)), the median SID is universally
lower than Dagma or Notears+. Similarly, the $F_{1}$ score is higher for all
but $D=50$ nodes, and the SHD is significantly lower for $d\leq 30$ nodes.
While Dagma-DCE is indeed more computationally expensive than Dagma, it is
still substantially faster than Notears across all graph sizes.

In the MLP experiments (Fig. 2(b)), our results are somewhat less impressive,
but still competitive. This is despite what we argue is a generating process
that gives a distinct advantage to Dagma and Notears+: by recalling that these
methods perform thresholding on the norm of layer weights, sampling weights
from $\mathcal{U}\left((-2.0,-0.5)\cup(0.5,2.0)\right)$ builds-in
identifiability in a perhaps unrealistic manner.

Finally, without asserting what the “correct” magnitude of weighted adjacency
matrices is, it is interesting to note that Dagma and Dagma-DCE return
adjacency matrices in which the relative magnitudes differ considerably. We
measure this with Kendall’s $\tau$-b and the Spearman correlation coefficient,
which measure rank correlation. These both take on values in $[-1,1]$, where a
value of $\pm 1$ indicates complete agreement or disagreement, and a value of
$0$ indicates that orderings are unrelated. The resulting values for our
synthetic experiments can be found in Table 1.

Table 1: Rank correlations between adjacency matrix entries of Dagma and Dagma-DCE across all trials for synthetic datasets, given as the mean plus-minus one standard deviation. Dataset | Kendall’s $\tau\text{-b}$ | Spearman’s $\rho$
---|---|---
Additive GPs | $0.40\pm 0.09$ | $0.53\pm 0.11$
MLP | $0.55\pm 0.06$ | $0.74\pm 0.07$

##  7 CONCLUSION

We reformulated the objectives of Notears+ and Dagma in order to extract a
principled, interpretable weighted adjacency matrix. With this comes
principled, interpretable notions of an $\ell_{1}$ sparsity penalty, and
thresholding the resulting graph.

A drawback of this method is its computational complexity, requiring a
backward pass over the Jacobian; however, empirical results indicate it is
still much faster than some other competing methods that use an augmented
Lagrangian scheme for optimization. Furthermore, we were able to achieve
state-of-the-art performance on synthetic datasets with small-to-medium-sized
graphs. Interesting directions for future work include applications to real
datasets and a formal workflow for determining sparsity and thresholding
hyperparameters.

## References

  * [1] Judea Pearl,  Causality,  Cambridge University Press, 2009.
  * [2] J. Peters, D. Janzing, and B. Scholkopf,  Elements of Causal Inference: Foundations and Learning Algorithms,  Adaptive Computation and Machine Learning Series. MIT Press, 2017.
  * [3] Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman,  Causation, Prediction, and Search,  MIT Press, 2000.
  * [4] Jiji Zhang and Peter Spirtes,  “Detection of unfaithfulness and robust causal inference,”  Minds and Machines, vol. 18, pp. 239–271, 2008.
  * [5] David Maxwell Chickering,  “Optimal structure identification with greedy search,”  Journal of Machine Learning Research, vol. 3, pp. 507–554, 2002.
  * [6] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing,  “DAGs with NO TEARS: Continuous optimization for structure learning,”  Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [7] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing,  “Learning sparse nonparametric DAGs,”  in International Conference on Artificial Intelligence and Statistics. PMLR, 2020, pp. 3414–3425.
  * [8] Kevin Bello, Bryon Aragam, and Pradeep Ravikumar,  “DAGMA: Learning DAGs via M-matrices and a log-determinant acyclicity characterization,”  Advances in Neural Information Processing Systems, vol. 35, pp. 8226–8239, 2022.
  * [9] Kurt Butler, Guanchao Feng, and Petar M Djurić,  “A differential measure of the strength of causation,”  IEEE Signal Processing Letters, vol. 29, pp. 2208–2212, 2022.
  * [10] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan,  “A linear non-Gaussian acyclic model for causal discovery.,”  Journal of Machine Learning Research, vol. 7, no. 10, 2006.
  * [11] Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf,  “Causal discovery with continuous additive noise models,”  Journal of Machine Learning Research, vol. 15, pp. 2009–2053, 2014.
  * [12] Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden,  “D’ya like DAGS? A survey on structure learning and causal discovery,”  ACM Computing Surveys, vol. 55, no. 4, pp. 1–36, 2022.
  * [13] OEIS Foundation Inc.,  “The On-Line Encyclopedia of Integer Sequences,” 2023,  Published electronically at http://oeis.org.
  * [14] Dennis Wei, Tian Gao, and Yue Yu,  “DAGs with No Fears: a closer look at continuous optimization for learning Bayesian networks,”  Advances in Neural Information Processing Systems, vol. 33, pp. 3895–3906, 2020.
  * [15] Joachim Bona-Pellissier, François Bachoc, and François Malgouyres,  “Parameter identifiability of a deep feedforward ReLU neural network,”  arXiv preprint arXiv:2112.12982, 2021.
  * [16] Alois Kufner and Bohumír Opic,  “How to define reasonably weighted Sobolev spaces,”  Commentationes Mathematicae Universitatis Carolinae, vol. 025, no. 3, pp. 537–554, 1984.
  * [17] C.E. Rasmussen and C.K.I. Williams,  Gaussian Processes for Machine Learning,  Adaptive Computation and Machine Learning series. MIT Press, 2005.
  * [18] Xiangyu Sun, Oliver Schulte, Guiliang Liu, and Pascal Poupart,  “NTS-NOTEARS: Learning nonparametric DBNs with prior knowledge,”  in International Conference on Artificial Intelligence and Statistics. PMLR, 2023, pp. 1942–1964.
  * [19] Diederik Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  in International Conference on Learning Representations (ICLR), San Diega, CA, USA, 2015.
  * [20] Dominik Janzing, David Balduzzi, Moritz Grosse-Wentrup, and Bernhard Schölkopf,  “Quantifying causal influences,”  The Annals of Statistics, vol. 41, no. 5, pp. 2324–2358, 2013.
  * [21] Elliott Sober,  “Apportioning causal responsibility,”  The Journal of Philosophy, vol. 85, no. 6, pp. 303–318, 1988.
  * [22] Nattaporn Plub-in and Jitkomut Songsiri,  “Estimation of Granger causality of state-space models using a clustering with Gaussian mixture model,”  in 2019 IEEE International Conference on Systems, Man and Cybernetics (SMC). IEEE, 2019, pp. 3853–3858.
  * [23] Khashaye Bayati, Karthikeyan Umapathy, and Soosan Beheshti,  “Granger causality pattern learning equipped with noise invalidation soft thresholding,”  in 2023 31st European Signal Processing Conference (EUSIPCO). IEEE, 2023, pp. 2028–2032.
  * [24] Jonas Peters and Peter Bühlmann,  “Structural intervention distance for evaluating causal graphs,”  Neural Computation, vol. 27, no. 3, pp. 771–799, 2015.