本章主要讲解一下在连续特征值的情况下如何进行因果图的挖掘。之前的介绍中节点的属性大部分以是否值为主。
之前讲的优化的框架是如下的公式表达。

$$max\ F(A) \\ s.t.G(A) \in DAGs \tag{1.1}$$

F(A)表示图的打分， A 为因果结构的邻接矩阵,G(A)为A 所对应的有向无环图，DAGS 是所有可能的DAG组成的离散空间.公式1.1能够在低维做一个优化，并且取得比较好的结果，但是本身公式1.1是一个非凸的优化问题，无法直接求解。为解决这个问题，学者们提出了一种连续优化方法口: 基于迹指数和增广拉格朗日的非组合优化结构学习( Non-combinatorial Optimization via TraceExponential and Augmented lagRangian for Structure learning，NOTEARS)。该方法采用数值解法，将 DAG 学习问题转化为约束的数值优化问题，下面我们将介绍该方法。

模型构造

假设数据 $D \in R^{n \times N}$包含 n 个独立同分布的观测样本，其中每个观测样本由N个变量组成，即 $V\{V_{1},..,V_{N}\}$,而数据对应的总体分布为 P(V)。结构学习的任务是从所有结构中找出服从数据分布 P(V)的结构。 连续优化方法引入结构等式模型 ( Structural Equation ModelSEM)的概念，以描述变量之间的关联及相互作用。形式上，SEM 由两个部分变量组成:

$$V_{j}=A_{j}^{T}V_{-j} + U_{j} \tag{2.1}$$

$V_{-j}$表示节点集合中除去$V_{j}$的所有节点，$U_{j}$表示偏置。至于这个表示方式一般使用权重邻接矩阵的方式，也就是$A_{ij}$表示邻接矩阵A中，i节点和j节点之间的权重。通过矩阵的表达， 简化了计算的方式。

$$V_{j}=A_{j}^{T}V + U_{j} \tag{2.2}$$

这里的V表示为矩阵的除了$V_{j}$以外的行向量，而$A_{j}^{T}$表示权重。很明显这里有一个可优化的方程。

$$l(A,D)=|\frac{1}{2n}(D-DA)|^{2} \tag{2.3}$$

也就是真实值与其他的特征与权重的乘积的差别，这就是一种损失函数。其实就是最小二乘法的损失函数，如果将这个最小二乘法损失降低到最小，那么我们就获得了真实的因果结构。此外对于一个回归性问题，G的分布应该是稀疏的，要加入一个稀疏正则项，从而得到最后的损失函数为

$$F(A)=l(A;D)+\lambda ||A||=\frac{1}{2n}|D-DA|^{2}+\lambda ||A|| \tag{2.4}$$

这个时候优化方程就变成了如下的样子

$$min\ F(A) \\ s.t. G(A) \in DAGs \tag{2.5}$$

很明显接下来就要解决的问题就是无环的约束问题啦， 怎么能通过量化的方式表达呢？

无环的约束方法

我们先来许愿， 希望这个无环的判断长成什么样子呢？ 我先来许愿

1. 有一个函数h(A)=0，就表示无环，也就是h(A)越大也能表示环越多。
2. 函数h的光滑的
3. 导数是便于计算的
   这里咱们直接给一个定理，如果$B= A \odot A$，A是权重矩阵，对应图的G位有向无环图，当且仅当$h(A)=tr(e^{B})-N=0$, 其中$\odot$表示哈达码积，$e^{B}$表示矩阵指数，tr表示矩阵的迹，并且有

$$tr(e^{B})=tr(I)+tr(B)+tr(\frac{B^{2}}{2!})+tr(\frac{B^{3}}{3!}) \tag{2.6}$$

可以进一步看上面这个公式2.6。tr(I)=N永远成立，也就是tr(B)=0才能保证上面的公式成立，这里有个小推论，矩阵中有长度为l的环必然有导致$B^{l}$的迹。
有了上面这个定理。公式2.5进一步可落地啦。

$$min\ F(A) \\ s.t. h(A)=0 \tag{2.7}$$

上面这个运筹方程的求解可以使用增广拉格朗日方法进行求解，使用对偶函数近似最优解。但是这一些都是基于线性的假设。
增广拉格朗日方法在原损失函数F(A)的基础上添加上h(A)的二次惩罚项，$\rho$表示惩罚系数，转化为如下的问题

$$Dual(\alpha)=min_{A \in R^{N \times N}}L^{\rho}(A,\alpha) \\ L^{\rho}(A,\alpha) = F(A)+\frac{\rho}{2} |h(A)|^{2}+\alpha h(A) \tag{2.8}$$

神经网络进行DAG学习

如果变量之间存在非线性关系，显然这是一个基础假设的崩塌，那么如何解决非线性的问题呢？
第一反应应该是神经网络类的方法，是的，神经网络在结构复杂到一定程度下能拟合任何的函数形式。那我们新的范式表达也十分简单。

$$V_{j}=MLP_{j}(V_{-j}, \theta) \tag{2.9}$$

这里需要注意一下，我们实际上是为每个$V_{j}$都有一个神经网络。也就是我们有多个MLP。整体结构如下图所示。

邻接矩阵表达

这里还要解决一个问题，有了神经网络网络的协助，最后图的结构还是需要一个邻接矩阵的表达，已经融入和神经网络中的参数，怎么表达邻接矩阵的关系呢？可以试想一些神经网络的传播过程，如果$MLP_{j}$表示$V_{j}$变量的生成网络， 如果参数的矩阵的第k列都是0，表达什么含义呢？说明$V_{k}$这个变量无法传递到下一层，是不是说明$V_{j}$和$V_{k}$是独立的，且没有边关系。这里公式2.8就可以通过神经网路的方式进行表示。

$$min_{\theta} (\frac{1}{n} \sum_{j}^{N}(l(V_{j},MLP_{j}(V_{-j},\theta_{j}))+\lambda|W_{j}|)+\frac{\rho}{2}|h(A(\theta))|^{2}+\alpha|h(A(\theta))|) \tag{2.10}$$

进一步给迭代的规则

$$\theta^{k}=argmin L_{\rho^{k}}(\theta,\alpha^{k}) \\ \alpha^{k+1}=\alpha^{k}+\rho^{k} h(A^{k}) \\ \rho^{k+1}=\eta \rho^{k}\ if|h(A^{k})|>\gamma |h(A^{k-1})| \\ \rho^{k+1}=\rho^{k}$$

可见F(A)可以通过神经网络拟合，通过神经网路的学习，从生成模型中提取出邻接矩阵，在进行阈值的筛选就可以得到对应的有向无环图。

DAG-GNN

之前的方式是通过MLP当时学习DAG，接下来的这一个章节学习的是通过GNN学习权重邻接矩阵，并将DAG的学习问题转化为无环约束的数字优化问题。

问题转化

首先还是使用结构等式框架（SEM）为基础，回顾线性表达,D为数据样本矩阵，U是噪声矩阵，A是邻接矩阵。

$$D=A^{T}D+U \tag{3.1}$$

当时数据对应的变量中所有的变量按照拓扑排序的方式，这权重矩阵A是一个严格上三角矩阵

$$D=(I-A^{T})^{-1} U \tag{3.2}$$

进一步可以表示为$D=f_{A}(U)$，以U为输入，返回D为输出的图神经网络的形式

$$D=f_{2}((I-A^{T})^{-1} f_{1}(U) ) \tag{3.3}$$

公式3.2中，$f_{1},f_{2}$表示参数函数，可以分别对D和U进行线性和非线性转换。公式3.3可以看做是一个解码器，学习DAG的任务变成了优化D的生成模型和学习权重邻接矩阵A的问题。
这个时候如果给定样本D和U的分布形式（高斯分布但是均值和方差未知），希望通过公式3.3得到D的生产模型，并通过优化生成模型得到邻接矩阵A。这个时候可以通过最大似然的方式求得D的分布形式。 也就是最大化log p(D),也就是最大化D的概率。因为还有U的参数，所以是最大化如下的公式表达。

$$\frac{1}{n} \sum_{k=1}^{n}log p(D^{k})=\frac{1}{n}\sum_{k=1}^{n} \int p(D^{k}|U)p(U) dU \tag{3.4}$$

这个时候发现上面的公式很难计算，这个可以采用变分编码器（VAE）的方式拟合生成机制，进而获得一个接近真实分布的变分后验。

如上图所示，变分解码器将输入的数据用解码器压缩到低维空间，其做法是将原始数据D经过编码器得到因变量U的分布所需的各项元素，得到因变量U，然后得到重构数据D‘。
在使用VAE的时候，给定一个样本$D^{k}$,编码器将其编码成$q(U|D^{k})$的隐变量U，然后经过解码器从U从重构出来$p(D^{k}|U)$的$D^{k}$，则编码器表示为如下形式(对应公式3.3)

$$U=f_{4}((I-A^{T})f_{3}(D)) \tag{3.5}$$

解码器表示为

$$D=f_{2}((I-A^{T})^{-1}f_{1}(D)) \tag{3.6}$$

令$f_{3}$为MLP，而$f_{4}$表示恒等映射，$f_{2}$为MLP，$f_{1}$表示恒等映射，如下图所示

上图的神经网络可以形式化如下的方式

$$[U均值|U标准差]=(I-A^{T}) MLP(D,W_{1},W_{2}) \\ MLP(D,W_{1},W_{2})=relu(DW_{1})W_{2} \\ [D均值|D标准差]=MLP((I-A^{T})^{-1}U, W_{3},W_{4})$$

这里还需进一步进入一个无环约束定理，
令$A \in R^{N \times N}$为一个有向图对应的邻接矩阵，对任意的$\alpha>0$，该有向无环图当且仅当如下的条件成立

$$tr[(I+\alpha A \circ A)^{N}]-N=0 \tag{3.7}$$

从而根据重构损失和无环约束就可以表示为如下的表达式

$$min_{A,\theta}\ f(A, \theta) \\ s.t. h(A)=tr[(I+\alpha A \circ A)^{N}]-N=0 \tag{3.8}$$

问题中未知的包括邻接矩阵A，$\theta=[W_{1},W_{2},W_{3},W_{4}]$,进一步使用增广矩阵的方式

$$L(A,\theta)=f(A,\theta)+\alpha h(A)+\frac{\rho}{2} |h(A)|^{2} \tag{3.9}$$

类似的更新公式

$$(A^{k}, \theta)=argmin L(A,\theta) \\ \alpha^{k+1}=\alpha^{k}+\rho^{k} h(A^{k}) \\ \rho^{k+1}=\eta \rho^{k}\ if|h(A^{k})|>\gamma|h(A^{k-1})| \\ \rho^{k+1}=\rho^{k} \tag{3.10}$$

优化完成后就可以得到权重矩阵A和相应的有向无环图。因为使用神经网络重构数据切显示的定义了邻接矩阵，所以此处就要更新矩阵权重，又要更新神经网络参数。

总而言之

上面介绍了一些新的方法进行DAG的学习，当然还有些方式也能起到类似的学习下效果，例如基于对抗网络的因果结构学习SAM，如果感兴趣可以查阅相关的资料。