【NeurIPS】神经成像时空回归模型的有效分层贝叶斯推理 Efficient Hierarchical Bayesian Inference

#NeurIPS #

今天分享的是NeurIPS 2021的一篇论文《神经成像时空回归模型的有效分层贝叶斯推理》

原文链接：https://arxiv.org/abs/2111.01692

摘要

神经影像学及其他领域的一些问题需要对多任务稀疏层次回归模型的参数进行推断，例如M/EEG逆问题、基于任务的功能磁共振成像分析的神经编码模型、气候科学。在这些领域中，需要推断的模型参数和测量噪声都可能表现出复杂的时空结构。现有的工作要么忽略了时间结构，要么推理方案需要的计算量很大。为了克服这些限制，作者设计了一个新的、灵活的层次贝叶斯框架，在其中模型参数、噪声的时空动力学被建模为具有克罗内克乘积的协方差结构。作者框架中的推理基于优化-最小化优化算法，并保证了收敛性。作者提出的高效算法针对时间自协方差矩阵利用了黎曼几何原理。对于由托普利茨矩阵描述的平稳动力学，采用了循环嵌入理论。作者证明了凸边界的性质，并推导了所得算法的更新规则。通过对来自M/EEG的合成和真实的神经数据上进行实验，作者证明了他们的方法可以提高性能。

1.介绍

在神经科学、神经成像等领域，经常对回归问题使用概率图形模型，其中模型参数和测量噪声都可能具有复杂的时空相关结构。本文的重点是脑源成像(BSI)，即利用非侵入性磁电或脑电图(M/EEG)传感器重建脑源活动。虽然之前在BSI回归问题中已经提出了许多包含时空相关性的生成模型，但推理解决方案要么实施了特定的简化和模型约束，要么总体上忽略了时间相关性，因此没有完全解决其固有的时空问题结构。

2.贡献

• 作者提出了一种新的和有效的算法框架，使用Type-II贝叶斯回归，明确地考虑了模型系数和测量噪声的时空协方差结构。
• 作者重点研究了脑源重建中的M/EEG逆问题，采用了一种多任务回归的方法，并将脑源重建问题表述为一个概率生成模型。
• 作者对具有克罗内克结构的协方差矩阵，利用黎曼流形上测地线凸性的概念，推导出了健壮、快速和高效的优化-最小化(MM)优化算法，用于具有可证明收敛保证的模型推理。
• 作者除了推导出了针对正定时间协方差矩阵的更新规则外，还假设时间协方差矩阵具有Toeplitz结构，这样可以在所提出的MM框架内产生计算效率较高的更新规则。

3.时空生成模型

作者重点研究脑源成像（BSI）问题，采用了多任务回归的方法，一个多任务线性回归问题的数学表述为$Y_g=L{X}_g+E_{g}forg=1,...,G(1)$

其中，

• $L\in R^{M\times N}$表示前向矩阵

• $X_g\in R^{N\times T}$表示一组脑源活动

• $Y_g\in R^{M\times T}$表示M/EEG测量结果

• $E_g\in R^{M\times T}$表示噪声

所以本文中BSI的目标是从给定L的M/EEG测量结果中推断出潜在的大脑活动。

作者考虑将多任务线性回归问题表示为概率图形模型，给出了两个假设：

• 假设1：脑源的时空结构可以用高斯分布来建模，则$p(X_g|\Gamma ,B)\sim N(0,{\Sigma }_0)$,其中${\Sigma }_0=\Gamma \otimes B$。$\Gamma$和B分别表示空间和时间协方差矩阵，并设置$\Gamma$是对角矩阵，B为大小T×T的正定矩阵。

• 假设2：噪声的时空结构可以用高斯分布来建模，则协方差矩阵${\Sigma }_e=\Lambda \otimes {\rm Y}$。$\Lambda$和${\rm Y}$分别表示空间噪声和时间噪声协方差矩阵，设置$\Lambda$是对角矩阵，并在这里假设噪声和脑源具有相同的时间结构，所以${\rm Y}=B$。

其中，⊗表示克罗内克乘积，指两个任意大小的矩阵间的运算。

则本文的概率图形模型如下：

• $\Lambda$是空间噪声的协方差矩阵，它和时间协方差B的克罗内克积形成了${\Sigma }_e$；$\Gamma$是空间协方差，和时间协方差B的克罗内克积形成了${\Sigma }_0$。
• 而${\Sigma }_e$是独立的高斯噪声$E_g$的方差，${\Sigma }_0$是脑源活动$X_g$的方差。
• 最后加上矩阵L，形成多任务回归问题的$Y_g$。
• 通过这样的方式把多任务回归问题表示成概率模型，同时把$X_g$和$E_g$划分成了时间与空间的结构。

4.算法实现

4.1 重要公式

时空生成模型中包含三个未知参数$\Gamma ,\Lambda ,B$，本文提出的算法可以在交替迭代过程中对其进行优化。在此给出两个重要公式。

脑源的后验分布$x_g$是高斯分布，即$p(x_g|y_g,\Gamma ,\Lambda ,B)\sim N({\overline{x}}_g,{\Sigma }_x)$

${\overline{x}}_g$是它的均值，${\Sigma }_x$是它的方差。公式如下：

$${\overline{x}}_g={\Sigma }_0{D}^T(\Lambda \otimes B+D{\Sigma }_0{D}^T{)}^{-1}{y}_g={\Sigma }_0{D}^T{\tilde{\Sigma }}_y^{-1}{y}_g(2)$$

$${\Sigma }_x={\Sigma }_0-{\Sigma }_0{D}^T{\tilde{\Sigma }}_y^{-1}D{\Sigma }_0(3)$$

4.2 流形学习

针对参数B的更新，文章给出了直观的几何表示：

$B\in S_{++}$,这是一个正定矩阵集，对B的更新采用了流形学习的思想，图中展示了{$B^k,M_{time}^k$}这个矩阵对之间的测地线路径，通过计算这个矩阵对之间的几何平均值来更新$B^{k+1}$。

其中流形学习指的是，假设数据在高维空间的分布位于某一更低维的流形上，基于这个假设来进行数据的分析。

4.3 full Dugh

文章给出了算法full Dugh的具体实现

该算法中涉及公式6，公式9，公式7，公式10、11、12，如下：

$$M_{time}^k:=\frac{1}{MG}\sum _{g=1}^G{Y}_g^T({\Sigma }_y^k{)}^{-1}{Y}_g(6)$$

$$M_{SN}^k:=H^k{\Phi }^T({\Sigma }_y^k{)}^{-1}{M}_{space}^k({\Sigma }_y^k{)}^{-1}\Phi H^k,with(9)$$

$$M_{space}^k:=\frac{1}{TG}\sum _{g=1}^G{Y}_g(B^k{)}^{-1}{Y}_g^T(9)$$

$$B^{k+1}\leftarrow (B^k{)}^{\frac{1}{2}}((B^k{)}^{\frac{-1}{2}}{M}_{time}^k(B^k{)}^{\frac{-1}{2}}{)}^{\frac{1}{2}}(B^k{)}^{\frac{1}{2}}(7)$$

$$H^{k+1}=diag(h^{k+1}),{h_i}^{k+1}\leftarrow \sqrt{\frac{{g_i}^k}{{z_i}^k}}fori=1,...,N+M,where(10)$$

$$g:=diag(M_{SN}^k)(11)$$

$$z:=diag({\Phi }^T({\Sigma }_y^k{)}^{-1}\Phi )(12)$$

• 公式6的作用是对$M_{time}^k$进行定义
• 公式9的作用是对$M_{SN}^k$进行定义
• 公式7的作用是对B进行更新
• 公式10、11、12的作用是对H进行更新

所以，对后验分布的最终估计${\overline{x}}_g$可以从模型参数H、B的任何初始猜测开始得到，在更新公式(7)(10、11、12)和(2)之间进行迭代。

4.4 thin Dugh

文章给出了算法thin Dugh的具体实现

该算法中涉及公式16、17、18，公式19，如下：

$${p_l}^{k+1}\leftarrow \sqrt{\frac{{{\hat{g}}_l}^k}{{{\hat{z}}_l}^k}}forl=1,...,L,where(16)$$

$$\hat{g}:=diag(P^k{Q}^H(B^k{)}^{-1}{M}_{time}^k(B^k{)}^{-1}Q{P}^k)(17)$$

$$\hat{z}:=diag(Q^H(B^k{)}^{-1}Q)(18)$$

$${\overline{x}}_g=(\Gamma \otimes B)D^T{\tilde{\Sigma }}_y^{-1}{y}_g=tr(QP(\Pi \odot Q^H{Y_g}^T{U}_x)({U_x}^{T}L{\Gamma }^T))(19)$$

• 该算法中的B是正定托普利茨矩阵
• 公式19的作用是给出对${\overline{x}}_g$的计算
• 公式10、11、12的作用是对H进行更新
• 公式16、17、18的作用是对p进行更新

所以，对后验的计算通过在更新公式(16、17、18)(10、11、12)和(19)之间迭代得到。

5.实验对比

四种不同脑源成像方案的比较：

在宽范围的信噪比、不同的时间样本数量、不同的AR模型顺序三种情况下对四种脑源成像的算法进行了性能评估，EMD（陆地移动距离）和TCE（时程相关误差）是进行性能评价的两个度量标准。

可以看到，Dugh在EMD度量方面取得了优越的性能，而在TCE方面则具有竞争力。

6.总结

• 针对脑源成像问题，采用多任务线性回归，提出了脑源和噪声的时空生成模型
• 提出了Full Dugh和Thin Dugh两个算法
• 在真实脑磁图和脑电图数据中进行了实验，表现出稳健和具有竞争力的性能
• 可以用于其他领域，如高斯过程回归中的核宽度、字典学习问题中的字典元素需要从数据中推断出来