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

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

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

作者重点研究脑源成像（BSI）问题采用多任务回归的方法，多任务线性回归的数学表达为 Y g = L X g E g ?? f o r ?? g = 1 , . . . , G ?? ( 1 ) Y_g=LX_g E_g\;for\;g=1,...,G\;(1) Yg=LXg Egforg=1,...,G(1)

其中，

• L ∈ R M × N L\in R^{M×N} L∈RM×N表示前向矩阵

• X g ∈ R N × T X_g \in R^{N×T} Xg​∈RN×T表示一组脑源活动

• Y g ∈ R M × T Y_g \in R^{M×T} Yg​∈RM×T表示M/EEG测量结果

• E g ∈ R M × T E_g \in R^{M×T} Eg​∈RM×T表示噪声

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

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

• 假设1：脑源的时空结构可以用高斯分布来建模，则 p ( X g ∣ Γ , B ) ∼ N ( 0 , Σ 0 ) p(X_g|\Gamma,B)\sim N(0,\Sigma_0) p(Xg​∣Γ,B)∼N(0,Σ0​),其中 Σ 0 = Γ ⊗ B \Sigma_0=\Gamma⊗B Σ0​=Γ⊗B。 Γ \Gamma Γ和B分别表示空间和时间协方差矩阵，并设置 Γ \Gamma Γ是对角矩阵，B为大小T×T的正定矩阵。

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

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

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

• Λ \Lambda Λ是空间噪声的协方差矩阵，它和时间协方差B的克罗内克积形成了 Σ e \Sigma_e Σe​； Γ \Gamma Γ是空间协方差，和时间协方差B的克罗内克积形成了 Σ 0 \Sigma_0 Σ0​。
• 而 Σ e \Sigma_e Σe​是独立的高斯噪声 E g E_g Eg​的方差， Σ 0 \Sigma_0 Σ0​是脑源活动 X g X_g Xg​的方差。
• 最后加上矩阵L，形成多任务回归问题的 Y g Y_g Yg​。
• 通过这样的方式把多任务回归问题表示成概率模型，同时把 X g X_g Xg​和 E g E_g Eg​划分成了时间与空间的结构。

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

脑源的后验分布 x g x_g xg​是高斯分布，即 p ( x g ∣ y g , Γ , Λ , B ) ∼ N ( x ‾ g , Σ x ) p(x_g|y_g,\Gamma,\Lambda,B)\sim N(\overline{x}_g,\Sigma_x) p(xg​∣yg​,Γ,Λ,B)∼N(xg​,Σx​)

x ‾ g \overline{x}_g xg​是它的均值， Σ x \Sigma_x Σx​是它的方差。公式如下：

x ‾ g = Σ 0 D T ( Λ ⊗ B + D Σ 0 D T ) − 1 y g = Σ 0 D T Σ ~ y − 1 y g    ( 2 ) \overline{x}_g=\Sigma_0D^T (\Lambda⊗B+D\Sigma_0 D^T)^{-1} y_g=\Sigma_0D^T \tilde{\Sigma}_y^{-1} y_g\;(2) xg​=Σ0​DT(Λ⊗B+DΣ0​DT)−1yg​=Σ0​DTΣ~y−1​yg​(2)

Σ x = Σ 0 − Σ 0 D T Σ ~ y − 1 D Σ 0    ( 3 ) \Sigma_x=\Sigma_0-\Sigma_0D^T\tilde{\Sigma}_y^{-1}D\Sigma_0\;(3) Σx​=Σ0​−Σ0​DTΣ~y−1​DΣ0​(3)

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

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

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

文章给出了算法full Dugh的具体实现 该算法中涉及公式6，公式9，公式7，公式10、11、12，如下：

M t i m e k : = 1 M G ∑ g = 1 G Y g T ( Σ y k ) − 1 Y g    ( 6 ) M^k_{time}:=\frac1{MG}\sum_{g=1}^GY_g^T(\Sigma_y^k)^{-1}Y_g\;(6) Mtimek​:=MG1​g=1∑G​YgT​(Σyk​)−1Yg​(6)