一种基于可调系数M矩阵逆的因果发现方法及系统

技术领域

本发明涉及机器学习及数据挖掘技术领域，特别地涉及一种基于可调系数M矩阵逆的因果发现方法及系统。

背景技术

随着数据获取手段的多样化和数据的日益丰富，社会学、生物医学、经济学等领域积累了海量的观测数据，这些数据是对事实或现象直接观测的结果，其中蕴含了数据背后真实的分布规律。近年来，如何直接从一组变量的观测数据中发现变量间的因果关系已成为统计学和人工智能等领域的研究热点之一。有向无环图(DAG)被用于描述多个变量之间的因果关系，其中节点表示随机变量，有向边表示变量间的因果关系。因果结构发现的目标即为从观测的样本数据中学习一个能准确刻画变量间因果关系的DAG，通常被称为DAG结构学习或贝叶斯网络结构学习问题。

基于评分函数从观测数据中学习变量间因果结构的算法大致可归为以下两类：基于组合优化的评分函数算法和基于连续优化的评分函数算法，前者是一个NP-难问题，只适用于节点数较小的图；后者提出一个刻画DAG结构特征的连续函数作为无环性约束，进而利用连续优化求解真实DAG结构。基于连续优化的评分函数算法虽然具有优美的数学形式和性质，但关于无环性约束中高阶项的系数分配出现了两种不同的观点：一方认为高阶项需要“小”系数以避免数值爆炸，另一方认为高阶项需要“大”系数以避免梯度消失。双方观点皆有依据，因此，给高阶项分配固定的系数(“小”系数或“大”系数)不可能同时解决数值爆炸和梯度消失这两个问题，进而导致对应的因果发现算法的准确率和效率较低。

发明内容

为了解决现有技术的不足，本发明通过引入M矩阵推导出一个新的无环性约束连续函数，该函数能自适应调整高阶项的系数，并在此无环性约束函数基础上，基于增广拉格朗日框架提出了一种新的基于连续优化的因果发现算法。具体来说，在优化的初始阶段，候选图与目标DAG之间差异较大(即候选图的谱半径较大而目标DAG的谱半径为零)，此时给高阶项分配较“小”的系数，以避免数值爆炸问题。随着优化的进行，候选图逐渐趋近于目标图，即谱半径逐渐趋近于零，此时给高阶项分配较“大”的系数，以避免高阶项的环路信息和梯度信息消失问题。基于上述分析，本发明同时避免了优化过程中可能存在的数值爆炸和梯度消失问题，使得算法能更好地学习真实DAG结构，进而提升因果发现的准确率和效率。

为了实现上述目的，本发明采用如下技术方案：

一种基于可调系数M矩阵逆的因果发现方法，所述方法至少包括：

S1.获取观测数据样本并进行数据预处理，构建数据矩阵；

S2.初始化噪声变量矩阵，基于加权邻接矩阵并采用线性结构方程表示变量间的因果关系，确定优化的评分函数，并结合M矩阵的特性提出一个可调系数的无环性约束函数；

S3.利用增广拉格朗日框架将原带约束的优化问题转化为无约束优化问题；

S4.初始化加权邻接矩阵、拉格朗日乘子、二次惩罚项系数、无环性约束的参数，设定L1正则项系数、算法精度、权重阈值；

S5.计算加权邻接矩阵与其自身哈达玛积的谱半径及无环性约束的系数；

S6.将二次惩罚项系数增大十倍，并把加权邻接矩阵代入无约束优化问题，解出对应的拉格朗日乘子；

S7.将S6步骤解出的拉格朗日乘子代入无约束优化问题，解出对应的加权邻接矩阵；

S8.判断加权邻接矩阵是否符合算法精度要求，若符合则继续向下执行；否则重新执行步骤S5；

S9.对加权邻接矩阵进行阈值处理；

S10.输出最优的加权邻接矩阵对应的因果图。

本发明还提供一种基于可调系数M矩阵逆的因果发现系统，所述系统用于所述的基于可调系数M矩阵逆的因果发现方法，所属系统包括：

数据集预处理模块，用于获取观测数据样本并进行数据预处理，构建数据矩阵；

参数设定模块，用于设定噪声变量矩阵、加权邻接矩阵、评分函数、无环性约束函数、拉格朗日乘子、二次惩罚项系数、无环性约束的系数、L1正则项系数、算法精度、权重阈值；

迭代优化求解模块，用于求解将原带约束的优化问题利用增广拉格朗日框架转化后的无约束优化问题；

阈值处理模块，用于对加权邻接矩阵进行阈值处理，即将加权邻接矩阵中所有绝对值小于权重阈值的元素设为零；

可视化输出模块，用于输出最优的加权邻接矩阵对应的因果图。

与现有技术相比，本发明技术方案的有益效果是：

本发明在步骤S2中通过引入M矩阵并利用其特性提出了一个新的可调系数的无环性约束函数，将新无环性约束应用于因果发现算法可有效解决以往工作中存在的数值爆炸和梯度消失等问题，使算法可以更快地收敛且趋近真实DAG结构，进而提升算法的准确率和效率。

附图说明

图1表示本发明实施例中提出的基于可调系数M矩阵逆的因果发现方法的流程图；

图2表示本发明实施例中提出的基于可调系数M矩阵逆的因果发现系统的结构图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

对于本领域技术人员来说，附图中某些公知内容说明可能省略是可以理解的；

附图中描述的位置关系仅用于示例性说明，不能理解为对本专利的限制；

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例：

注意：本实施例中用小写斜体表示标量，小写加粗斜体表示向量，大写加粗斜体表示矩阵。

如图1所示的基于可调系数M矩阵逆的因果发现方法的流程图，参见图1，所述方法包括：

S1.获取观测数据样本并进行数据预处理，构建数据矩阵；

在本实施例中，设初始数据集为d维随机变量

S2.初始化噪声变量矩阵，基于加权邻接矩阵并采用线性结构方程表示变量间的因果关系，确定优化的评分函数，并结合M矩阵的特性提出一个可调系数的无环性约束函数；

具体过程为：

(1)噪声变量矩阵表示为

(2)加权邻接矩阵为

(3)表示变量间的因果关系的线性结构方程为X＝XW+N；

(4)所述评分函数Q(W,X)由均方误差损失函数与L1正则项构成，表达式为

(5)无环性约束函数是结合M矩阵逆和矩阵迹函数得到的函数h(W)，其表达式为

其中上述(5)细分为以下过程以说明本发明提出的无环性约束函数的优势：

(5a)由矩阵的幂级数可得，

(5b)h(W)幂级数的系数1/s

(5c)具体来说，在优化初始阶段，

S3.利用增广拉格朗日框架将原带约束的优化问题转化为无约束优化问题进行求解；

在本实施例中，原带约束的优化问题表述为：

其中

其中σ＞0为二次惩罚项的系数，

S4.初始化加权邻接矩阵、拉格朗日乘子、二次惩罚项系数、无环性约束的参数，设定L1正则项系数、算法精度、权重阈值；

在本实施例中，初始化加权邻接矩阵W的每一个元素服从均匀分布并控制在区间[-2.0,-0.5]∪[0.5,2.0]内，初始化拉格朗日乘子

S5.计算加权邻接矩阵与其自身哈达玛积的谱半径及无环性约束的参数；

在本实施例中，先计算出加权邻接矩阵与其自身哈达玛积的谱半径

S6.将二次惩罚项系数增大十倍，并把加权邻接矩阵代入无约束优化问题，解出对应的拉格朗日乘子；

在本实施例中，由于仅当h(W)＝0时公式(2)才与公式(1)等价，因此在优化过程中采用逐步增大二次惩罚项系数σ的策略，使其作为惩罚因子以确保h(W)＝0，即在每一次迭代优化过程中将二次惩罚项系数σ增大十倍。将加权邻接矩阵代入公式(2)后，此时公式(2)中只有拉格朗日乘子

S7.将S6步骤解出的拉格朗日乘子代入无约束优化问题，解出对应的加权邻接矩阵；

在本实施例中，将S6中求出的拉格朗日乘子

S8.判断加权邻接矩阵是否符合算法精度要求，若符合则继续向下执行；否则重新执行步骤S5；

在本实施例中，判断此时加权邻接矩阵W对应的无环性约束函数h(W)是否小于S4中设定的算法精度τ＝10

S9.对加权邻接矩阵进行阈值处理；

在本实施例中，步骤S9为将步骤S8中求出的加权邻接矩阵W进行阈值处理，即将W中所有绝对值小于权重阈值θ的元素设为零，以排除可能诱导环的虚假关联，经阈值处理后可以减少图中假因果边的出现并保持稀疏性，从而提高算法的准确性。

S10.输出最优的加权邻接矩阵对应的因果图。

在本实施例中，经过步骤S9得到加权邻接矩阵的最优解W

如图2所示，本发明还提供一种基于可调系数M矩阵逆的因果发现系统，所述系统用于所述的基于可调系数M矩阵逆的因果发现方法，所属系统包括：

数据集预处理模块，用于获取观测数据样本并进行数据预处理，构建数据矩阵；

参数设定模块，用于设定噪声变量矩阵、加权邻接矩阵、评分函数、无环性约束函数、拉格朗日乘子、二次惩罚项系数、无环性约束的系数、L1正则项系数、算法精度、权重阈值；

迭代优化求解模块，用于求解将原带约束的优化问题利用增广拉格朗日框架转化后的无约束优化问题；

阈值处理模块，用于对加权邻接矩阵进行阈值处理，即将加权邻接矩阵中所有绝对值小于权重阈值的元素设为零；

可视化输出模块，用于输出最优的加权邻接矩阵对应的因果图。

显然，本发明的上述实施例仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定了。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其他不同形式的变化或变动。这里无需也无法对所有实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。