Regularized matrix data clustering and its application to image analysis

Xu Gao; Weining Shen; Liwen Zhang; Jianhua Hu; Norbert J. Fortin; Ron D. Frostig; Hernando Ombao

doi:10.1111/biom.13354

生物计量学。作者手稿；PMC 2022年9月1日提供。

以最终编辑形式发布为：

生物计量学。2021年9月；77(3): 890–902.

2020年8月24日在线发布。数字对象标识：10.1111/生物13354

预防性维修识别码：项目管理委员会7884484

尼姆斯：美国国家卫生研究院1635244

PMID：32799339

正则化矩阵数据聚类及其在图像分析中的应用

徐高,¹ 沈伟宁,¹ 张丽文,² 胡建华,^三诺伯特·J·福廷,⁴ 罗恩·弗罗斯蒂格,^4,⁵和埃尔南多·奥姆巴奥⁶

作者信息版权和许可信息 PMC免责声明

关联数据

补充资料: 附录。
NIHMS1635244-附录.pdf（540K）
GUID:57EC3349-3ACF-4CA7-B345-E2C6A744AE87

摘要

我们提出了一种新的正则化混合模型来聚类矩阵值数据。该方法假设每个簇具有可分离的协方差结构，并对每个簇的平均信号施加稀疏结构（例如，低秩、空间稀疏）。我们将该问题表示为带有正则化项的矩阵正态分布的有限混合模型，然后开发了一种用于高效计算的期望最大化算法。理论上，我们证明了所提出的估计对于各种惩罚函数的选择是强一致的。仿真和脑信号研究的两个应用证实了该方法的卓越性能，包括比竞争对手更好的预测精度和解决方案的科学解释性。

关键词：聚类、成像、矩阵正态分布、混合模型、正则化、时频分析

1 |. 简介

在过去十年中，生成具有复杂统计结构的大容量数据集的技术取得了巨大发展。其中，矩阵值数据通常出现在大脑图像和信号中，其中采样单元可以被视为二维阵列（即矩阵），例如脑电图（EEG）和局部场电位（LFP）。这些信号通常是高维的，具有复杂的结构，如空间/时间相关性、低秩和稀疏性(高等。，2019年;胡等。，2019年;王等。, 2019). 本文的主要目标是提供一种新的方法来聚类矩阵值数据，同时考虑到它们的复杂结构。

聚类是统计学和许多科学应用中的一个基本问题，包括研究大脑功能及其对刺激的反应的研究(金，2015). 本文的一个主要动机是由合著者Fortin的实验室进行的非空间工作记忆实验，以研究气味顺序的神经元学习过程(艾伦等。, 2016). 海马神经元时间编码的发现扩展了我们对情景记忆神经生物学的基本理解，从而为阐明记忆损伤的潜在神经机制提供了跨物种基础。在整个实验过程中，来自同一气味端口的五种气味（表示为ABCDE）被呈现给大鼠。每种气味呈现都是由一个鼻子戳开始的。对大鼠进行测试，以正确识别气味呈现的顺序是按顺序（ABCDE）还是按顺序（如AABDE、ABCDD），方法是将它们的鼻子放在端口中，直到信号被确认或在信号发出之前撤回。将12个微电极植入大鼠海马，记录LFP。这里感兴趣的主要科学问题是探索和理解LFP信号中的（潜在）特征，这些特征与发展序列气味记忆的神经机制有关。为了证明聚类分析的威力，我们进行了探索性分析。图1（此图在本文电子版中以彩色显示，任何提及的颜色均指该版本）显示了一只老鼠跨越12个微电极的平滑LFP，用于五种气味（ABCDE）及其聚集在气味上的相关平均信号。在每个图中，x个-轴表示按间隔[0，1]缩放的时间，而年-轴表示不同的微电极（通道）。很明显，不同序列的平均模式差异很大，当我们比较每种气味中不同河豚的信号时，有很强的空间依赖性。在所有六个热峰中，我们观察到，1到8和9到12的河豚在河豚身上分别形成了气味A、B和D的两个不同的范式。此外，这两个“范式”也随着时间的变化而演变，这表明，本研究中的聚类分析将有助于揭示LFP中的潜在模式/结构，从而对其与不同气味的关系提供更多的见解。具体而言，理想的聚类方法应能够通过以下方式发现LFP中的“潜在特征”：（a）应确定行向（四极内）和列向（四足间）相关性，从而提供关于空间和时间相关性的见解；（b）它应该能够通过引入正则化项来揭示真实的平均差，从而提高LFP信号的信噪比；（c）它应该能够揭示数据固有的稀疏性本质（例如，检测图像信号的边界）；（d）通过将所提出的方法应用于时域和频域，它应该能够识别出仅仅通过对这些信号的视觉检查不容易识别的不同“结构”。

保存图片、插图等的外部文件。对象名为nihms-1635244-f0002.jpg

在单独的窗口中打开

图1

前6项（内存编码实验）：气味的平滑LFP ABCDE及其平均值（聚合所有气味），其中x个-轴表示重新缩放的时间和年-轴代表不同的河豚（通道）；底部4（大鼠中风研究）：中风前后所有600个时期中微电极10和20的时频图。此图在本文的电子版中以彩色显示，任何提及颜色的地方都指的是该版本

在本文中，我们将重点放在使用有限混合模型进行聚类，因为统计推断可以以计算效率高、概念简单的方式进行，并且结果具有良好的概率解释。现有方法，如双聚类(芝加哥等。, 2017)，分层聚类(尤安等。，2018年,2019)，光谱聚类(Ng公司等。, 2002)，和（随机向量的）混合模型不直接适用于矩阵数据分析，因为输入协变量集被视为向量，其中没有考虑矩阵结构及其可解释性(莱斯和奥格登，2007年). 此外，通过对矩阵进行矢量化，得到的输入空间维数可以非常大，即第页×q个矩阵将转换为pq值-维向量，这给计算和理论带来了额外的挑战。

为了解决上述问题，我们提出了一种新的惩罚混合模型来聚类矩阵值数据。该框架受到了维罗利（2011）和高等。(2019)其中，每个混合成分由矩阵正态分布表示，其协方差矩阵可以分解为两个独立的列和行协方差矩阵的Knocker积(Dawid，1981年;Dutilleul，1999年). 可分离协方差结构既提供了计算便利性，因为它有效地减少了协方差参数的数量，又提供了有用的实际解释，因为它根据时间域和空间域分离变量。此外，我们考虑采用具有三种不同规范（即，ℓ₁,ℓ₂以及核范数），以使该方法在捕获平均结构中的不同特征时具有灵活性和鲁棒性，从而增强正确识别聚类的能力。例如，核范数的使用为真实图像提供了有用的低阶近似(周和李，2014)和的使用ℓ₁-范数有助于检测图像边界(王等。, 2017). 我们引入了一种新的期望最大化（EM）类型的算法，可以对所有三种惩罚准则进行有效计算。理论上，我们通过以下方法证明了所提出估计量的强一致性结果范和李（2001）我们必须对其进行修改以适应矩阵值数据。

注意，矩阵正态分布及其可分离协方差结构是我们提出的混合模型的构建块。虽然不在我们当前论文的范围内，但可以推广矩阵正态分布，并考虑更一般的协方差模型，如尖峰协方差模型(多诺霍等。, 2018). 我们期望我们提出的对矩阵平均信号应用不同正则化准则的想法在这些场景中仍然有用。

论文的其余部分组织如下。在第2节，我们简要回顾了矩阵正态分布，然后介绍了所提出的惩罚矩阵正态混合模型。我们提出了一种新的EM类型的计算算法，并讨论了一种相关的单步长（OSL）算法。在第3节，我们通过显示正则化估计的一致性结果来研究该方法的理论性质。我们通过中的模拟研究来评估有限样本的性能第4节最后，我们将所提出的方法应用于分析从气味记忆和中风实验的成像研究中获得的两个LFP数据集第5节和6支持信息中提供了其他数值结果、技术证明以及计算代码和数据集的描述。

2|。方法

我们首先回顾了矩阵正态分布及其有限混合第2.1节和2.2矩阵正态分布的混合模型及其在三向数据分类中的应用已经由维罗利（2011）。我们将讨论并强调我们在以下方面的主要贡献第2.3节.

2.1 |. 矩阵正态分布

我们简要回顾了矩阵正态分布。该分布将用于对矩阵值图像数据建模。我们定义了第页×第页随机矩阵Y具有平均值的矩阵正态分布M（M）和协方差矩阵U型和V（V），表示为明尼苏达州_{第页,第页}(M（M）,U型,V（V）)，如果其密度函数为

（f） (Y（Y） ∣ M（M）, U型, V（V）) = \frac{经验 (- \frac{1}{2} 信托收据 ({V（V）}^{- 1} {(Y（Y） - M（M）)}^{T型} {U型}^{- 1} (Y（Y） - M（M）)))}{{(2 π)}^{第页 第页 / 2} {| V（V） |}^{第页 / 2} {| U型 |}^{第页 / 2}},

(1)

哪里 $M（M） \in ℝ^{第页 \times 第页}$ , $U型 \in ℝ^{第页 \times 第页}$ , $V（V） \in ℝ^{第页 \times 第页}$ ，矩阵U型和V（V）被视为协方差矩阵中的介于和之间的矩阵，·和tr分别表示矩阵的行列式和迹。在Gupta和Nagar（1999年），的等效定义(1)由提供

v（v） e（电子） c (Y（Y）) ~ N个 (v（v） e（电子） c (M（M）), V（V） \otimes U型),

(2)

哪里血管内皮细胞是列矢量化操作，⊗是Kronecker乘积。矩阵正态分布为图像和时空数据建模提供了通常的多元正态分布的自然扩展。例如，在时空统计中，矩阵中的行和列方向分别对应于空间和时间属性。协方差结构的可分性在空间统计等应用中特别有用(哈斯1995;Cressie 2015年)和电生理数据(高等。, 2020). 据信，这种假设大大减少了参数的数量，从而为大规模时空数据提供了一种高效且稳健的计算过程(Genton，2007年).

矩阵正态分布的统计推断通常通过似然函数进行。给出i.i.d.观察结果Y（Y）₁,Y（Y）₂…，Y（Y）_n个~明尼苏达州_{第页,第页}(M（M）,U型,V（V）)，对数似然函数为

Ş (M（M）, U型, V（V）) = - \frac{n个 第页 第页}{2} 日志 2 π - \frac{n个 第页}{2} 日志 | V（V） | - \frac{n个 第页}{2} 日志 | U型 | - \frac{1}{2} \sum_{我 = 1}^{n个} 信托收据 {{V（V）}^{- 1} {({Y（Y）}_{我} - M（M）)}^{T型} {U型}^{- 1} ({Y（Y）}_{我} - M（M）)} .

(3)

得到的最大似然估计量（MLE）M（M）,U型、和V（V）具体如下：

\hat{M（M）} = \bar{Y（Y）}, \hat{U型} = \frac{1}{n个 第页} \sum_{我 = 1}^{n个} ({Y（Y）}_{我} - \bar{Y（Y）}) {\hat{V（V）}}^{- 1} {({Y（Y）}_{我} - \bar{Y（Y）})}^{T型}, \hat{V（V）} = \frac{1}{n个 第页} \sum_{我 = 1}^{n个} {({Y（Y）}_{我} - \bar{Y（Y）})}^{T型} {\hat{U型}}^{- 1} ({Y（Y）}_{我} - \bar{Y（Y）}) .

(4)

保存图片、插图等的外部文件。对象名为nihms-1635244-f0001.jpg

请注意，这两者U型和V（V）只能识别为常数倍数(Dutilleul，1999年)，但他们的Kronecker产品U型⊗V（V）保持不变。没有用于的封闭式解决方案 $\hat{U型}$ 和 $\hat{V（V）}$ 。可以使用迭代算法以数值方式获取其值，如算法1中所述。

2.2 |. 矩阵正态混合模型

为了进行概率聚类，我们考虑了一个矩阵正态混合模型及其使用EM算法的推理。提供身份证(第页×第页)-维度观测Y（Y）₁…，Y（Y）_n个来自以下混合物K（K）矩阵正态分布_j个= (M（M）_j个,U型_j个,V（V）_j个)和重量π₁…，π_K（K）属于K（K）-维单纯形，用Δ表示_K（K）.然后混合物密度可以写为

\sum_{k个 = 1}^{K（K）} π_{k个} （f） (Y（Y） ∣ Θ_{k个}) = \sum_{k个 = 1}^{K（K）} π_{k个} M（M） {N个}_{第页, 第页} ({M（M）}_{k个}, {U型}_{k个}, {V（V）}_{k个}),

(5)

哪里（f）是具有平均值的矩阵正态分布M（M）_k个和协方差U型_k个和V（V）_k个定义如下(1). 我们使用θ=（θ₁, …, Θ_K（K）;π₁…，π_K（K）)表示（5）中的参数集合。那么log-likelihood函数是

ℓ_{o（o） b条 秒} (Θ) = \sum_{我 = 1}^{n个} 日志 {\sum_{j个 = 1}^{K（K）} π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})} .

(6)

EM算法(登普斯特等。, 1977)已广泛应用于统计的所有领域，并已被证明适用于高斯混合模型。因此，我们还使用EM算法进行推理，这是一种由期望（E）和最大化（M）步骤组成的迭代方法。

在E阶中，观测的后验概率Y（Y）_我属于j个根据贝叶斯定理得出第个簇，如下所示：

α_{我 j个} = \frac{π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})}{\sum_{我 = 1}^{K（K）} π_{我} （f） ({Y（Y）}_{我} ∣ Θ_{我})} .

(7)

在M步中，通过求解无约束优化问题获得参数向量的估计

{\hat{Θ}}_{j个} = \underset{Θ_{j个}}{参数 最大值} \sum_{我 = 1}^{n个} \sum_{j个 = 1}^{K（K）} α_{我 j个} 日志 {π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})} .

通过一些代数，我们得到了以下关系式：

{\hat{π}}_{j个} = \frac{\sum_{我 = 1}^{n个} α_{我 j个}}{n个}, {\hat{M（M）}}_{j个} = \frac{\sum_{我 = 1}^{n个} α_{我 j个} {Y（Y）}_{我}}{\sum_{我 = 1}^{n个} α_{我 j个}}, {\hat{U型}}_{j个} = \frac{\sum_{我 = 1}^{n个} α_{我 j个} ({Y（Y）}_{我} - {\hat{M（M）}}_{j个}) {\hat{V（V）}}_{j个}^{- 1} ({Y（Y）}_{我} - {\hat{M（M）}}_{j个})'}{第页 \sum_{我 = 1}^{n个} α_{我 j个}}, {\hat{V（V）}}_{j个} = \frac{\sum_{我 = 1}^{n个} α_{我 j个} ({Y（Y）}_{我} - {\hat{M（M）}}_{j个})' {\hat{U型}}_{j个}^{- 1} ({Y（Y）}_{我} - {\hat{M（M）}}_{j个})}{第页 \sum_{我 = 1}^{n个} α_{我 j个}} .

(8)

请注意 ${\hat{U型}}_{j个}$ 和 ${\hat{V（V）}}_{j个}$ ,j个= 1, …,k个没有闭合形式的解，但可以按照算法1中的类似步骤从数值上获得。

2.3 |. 惩罚矩阵正态混合模型

矩阵正规模型的混合先前已由讨论维罗利（2011）用于对三向阵列数据进行分类。然而，对于许多成像研究，需要考虑潜在的空间（矩阵）结构（例如第1节). 通过对平均矩阵信号使用惩罚函数，例如低阶近似，可以有效地模拟这种结构(周和李，2014)或基于全变分形式的惩罚(王等。, 2017). 在本文中，我们提出了一种惩罚方法，即在矩阵正态混合模型中的每个混合成分的均值上包含一个惩罚项。惩罚函数的形式为ℓ₁,ℓ₂，或平均矩阵的核范数M（M）₁…，M（M）_k个.惩罚函数的选择取决于每个簇的平均结构的领域知识，如稀疏性、平滑性和低秩(绿色，1990年)，它给出了易于解释的结果，也减轻了计算负担。具体来说，我们考虑惩罚对数似然函数

问 (Θ; λ) = \sum_{我 = 1}^{n个} 日志 {\sum_{j个 = 1}^{K（K）} π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})} - λ \sum_{j个 = 1}^{K（K）} P（P） ({M（M）}_{j个}),

(9)

哪里P（P）（·）是一些惩罚函数，例如ℓ₁,ℓ₂和核规范，以及λ≥0是调谐参数。对于ℓ₁和ℓ₂-惩罚，范数定义在向量化矩阵平均值上M（M）_j个对于核规范惩罚，定义为以下奇异值之和M（M）_j个中提出的模型(9)是通用的，事实上，设置λ=0等价于在中引入的没有正则化的矩阵法线的混合模型维罗利（2011）.

类似第2.2节，我们提出了一种改进的EM算法来估计参数。E步骤的进行方式与公式（7）M步骤归结为解决优化问题，

\hat{Θ} = \underset{Θ}{参数 最大值} \sum_{我 = 1}^{n个} \sum_{j个 = 1}^{K（K）} α_{我 j个} 日志 {π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})} - λ \sum_{j个 = 1}^{K（K）} P（P） ({M（M）}_{j个}) .

(10)

注意，解决方案 $\hat{Θ}$ 可能没有显式形式。兰格（1995）提出了一种与EM算法相关的梯度方法。它通过对牛顿法进行一次迭代来代替M步。替代方法，例如代理函数(兰格等。, 2000)和超松弛EM算法(于，2012)文献中也介绍了。潘和沈（2007）介绍ℓ₁-对混合单变量正态模型的平均参数进行惩罚。他们使用次梯度方法获得了M步的显式解。绿色（1990）开发了可应用于更一般情况的“OSL”算法。受上述方法的启发，我们开发了一种次梯度方法，当ℓ₁-使用规范和OSL方法ℓ₂和核规范。

对于ℓ₁-规范处罚，遵循以下类似推导潘和沈（2007）,M（M）_j个可以由更新

\hat{{M（M）}_{j个}} = 签名 ({\tilde{M（M）}}_{j个}) {(| {\tilde{M（M）}}_{j个} | - \frac{λ}{{\sum^{​}}_{我 = 1}^{n个} α_{我, j个}} {U型}_{我} 1_{第页 \times 第页} {V（V）}_{我})}_{+}, j个 = 1, \dots, K（K）,

(11)

哪里 ${\tilde{M（M）}}_{j个} = \frac{\sum_{我 = 1}^{n个} α_{我, j个} {Y（Y）}_{我}}{\sum_{我 = 1}^{n个} α_{我, j个}}$ 是的更新M（M）_我没有惩罚，B类₊=最大值(B类，0），1_{第页×第页}是所有1、sign（）和（.）的矩阵₊都是组件式运算符。

对于ℓ₂-范数惩罚，目标函数的偏导数 $问_{ℓ_{2}} (Θ)$ 是

\frac{\partial 问_{ℓ_{2}} (Θ)}{\partial {M（M）}_{j个}} = {U型}_{j个}^{- 1} \sum_{我 = 1}^{n个} α_{我, j个} ({Y（Y）}_{我} - {M（M）}_{j个}) {V（V）}_{j个}^{- 1} - 2 λ {M（M）}_{j个}, j个 = 1, \dots, K（K） .

因此，M（M）_j个可以由更新

{\hat{M（M）}}_{j个} = {\tilde{M（M）}}_{j个} - \frac{2 λ}{\sum_{我 = 1}^{n个} α_{我 j个}} {U型}_{j个} {M（M）}_{j个} {V（V）}_{j个},

(12)

哪里U型_j个,M（M）_j个,V（V）_j个是上一步的更新。

核规范处罚‖·‖_*定义为矩阵奇异值之和，类似的推导得出

{\hat{M（M）}}_{j个} = {\tilde{M（M）}}_{j个} - \frac{λ}{\sum_{我 = 1}^{n个} α_{我 j个}} {U型}_{j个} Φ_{j个} Ω_{j个}^{T型} {V（V）}_{j个},

(13)

哪里M（M）_j个具有奇异值分解 ${M（M）}_{j个} = Φ_{j个} Λ_{j个} Ω_{j个}^{T型}$ .核规范正则化的使用本质上等同于L（左）₁-奇异值的正则化M（M）。有关派生的更多详细信息(11)至(13)支持信息第3.3节中给出。

总之，估计过程涉及初始化算法以及E步和M步之间的交替算法。在这里，我们提供更多细节。

I.（初始化）

选择一个固定的聚类数K（K）。我们首先对原始矩阵值观测值进行矢量化Y（Y）₁…，Y（Y）_n个并应用K（K）-表示实现初始集群成员身份值，写为S公司₁…，S公司_K（K），其中S公司_j个= {我|Y（Y）_我在里面j个第个簇}。请注意K（K）-在这一步中也可以使用平均值，例如，将观测值随机分配给不同的簇。然后，对于每个簇，θ的初始值_j个可以通过以下方法获得第2.1节和π_j个可以通过以下方式直接估算 ${\hat{π}}_{j个} = | {S公司}_{j个} | / n个$ ，其中|S公司_j个|是的基数S公司_j个.

二、。（E步骤）

我们通过更新后验隶属度

α_{我 j个} = \frac{π_{j个} （f） ({Y（Y）}_{我} ∣ Θ_{j个})}{\sum_{我 = 1}^{K（K）} π_{我} （f） ({Y（Y）}_{我} ∣ Θ_{我})} .

三、（M步）

更新平均参数M（M）_j个关于各种处罚方程式（11）至(13)分别是。其他参数π_j个,U型_j个,V（V）_j个可以更新以下内容方程式（8）和算法1。

四、（停止标准）

重复第二步和第三步，直到达到一定的迭代次数或平均参数估计值发生变化M（M）_j个（根据Frobenius范数）低于某些预先规定的截止值。

四、（选择簇数）

集群中的一个关键问题是确定集群的数量。灵感来自Smyth（2000），我们采用交叉验证惩罚似然（CVPL）作为关键度量来考虑预测标准。我们分割数据集Y（Y）= {Y（Y）₁…，Y（Y）_n个}分为训练和测试组，表示为Y（Y）_火车,Y（Y）_测试，然后安装k个-混合模型（用于k个=1，…）启用Y（Y）_火车然后使用估计参数获得惩罚log-likelihood函数Y（Y）_测试，表示为问(Y（Y）_测试|Y（Y）_火车,k个). CVPL的一个好特性是它的期望值是真惩罚似然函数和k个-混合惩罚似然加上一些常数。考虑到这个度量，我们可以通过首先划分来定义CVPLY（Y）= (Y（Y）¹…，Y（Y）^L（左）)相等地变成L（左）随机分配，然后考虑

问 (L（左）, k个) = {L（左）}^{- 1} \sum_{我 = 1}^{L（左）} 问 ({Y（Y）}^{我} ∣ {Y（Y）}^{- 我}, k个),

哪里Y（Y）^−我是数据Y（Y）不包括Y（Y）^我.然后我们选择 $\hat{k个} = {argmax（最大值）}_{k个} 问 (L（左）, k个)$ 对于预选的L（左）在数值分析中，我们选择L（左）=3，即三重交叉验证，以便于计算。为选择较大的值L（左）也有可能，未来将有兴趣研究CVPL作为模型选择标准在不同值下的稳定性L（左）。我们将在模拟研究中为计算提供更多细节。

3|。理论

在本节中，我们研究了所提出的矩阵正规混合模型的理论性质，假设簇的真实数量，表示为K（K），已知。对混合模型（特别是高斯混合模型）的理论理解受到了关注（特别是，见陈，2017). 在这里，我们首先给出了矩阵正态混合模型MLE的一致性结果。众所周知，如果参数空间没有约束，则一元正态分布混合的可能性是无界的。因此，我们考虑一个约束参数空间 $Ψ^{{d日}_{1}, {d日}_{2}}$ 作为

Ψ^{{d日}_{1}, {d日}_{2}} = {Θ : {(π_{1}, \dots, π_{k个})}^{T型} \in Δ_{K（K）}, {M（M）}_{1}, \dots, {M（M）}_{k个} \in ℝ^{第页 \times 第页} \times {V（V）}_{1}, \dots, {V（V）}_{K（K）} \in ℝ^{第页 \times 第页}, {U型}_{1}, \dots, {U型}_{k个} \in ℝ^{第页 \times 第页} : \times \underset{1 \leq 小时 \neq j个 \leq k个}{最小值} ρ ({U型}_{小时} {U型}_{j个}^{- 1}) \geq {d日}_{1}, \times \underset{1 \leq 小时' \neq j个' \leq k个}{最小值} ρ ({V（V）}_{小时'} {V（V）}_{j个'}^{- 1}) \geq {d日}_{2}, ρ ({U型}_{我}) > 0, ρ ({V（V）}_{我}) > 0 对于 我 = 1, \dots, K（K） .},

(14)

哪里d日₁,d日₂是在（0，1]，Δ上定义的两个固定常数_K（K）表示K（K）-一维单纯形，以及ρ（·）表示矩阵的最小特征值。我们首先给出了矩阵正态混合模型的一致性结果。证据见Web附录C.

定理1。

让Y₁, … ,Y（Y）_n个 是矩阵正态分布混合物的随机样本，定义见(5)用参数索引Θ属于这个空间的 $Ψ^{{d日}_{1}, {d日}_{2}}$ 由定义(14).假设真参数值 $Θ_{0} \in Ψ^{{d日}_{1}, {d日}_{2}}$ 并表示 ${\hat{Θ}}_{n个}$ 作为解决方案 $参数 {最大值}_{Ψ^{{d日}_{1}, {d日}_{2}}} ℓ_{o（o） b条秒} (Θ)$ .然后 ${\hat{Θ}}_{n个}$ 收敛到Θ₀ 几乎可以肯定.

这里，θ₀= (π₁₀…，π_K（K）0;M（M）₁₀…，M（M）_K（K）0;V（V）₁₀⊗U型₁₀…，V（V）_K（K）0⊗U型_K（K）0)是真实混合物模型参数的集合， ${\hat{Θ}}_{n个} = ({\hat{π}}_{1}, \dots, {\hat{π}}_{K（K）}; {\hat{M（M）}}_{1}, \dots, {\hat{M（M）}}_{K（K）}; {\hat{V（V）}}_{1} \otimes {\hat{U型}}_{1}, \dots, {\hat{V（V）}}_{K（K）} \otimes {\hat{U型}}_{K（K）})$ ，以及 ${\hat{Θ}}_{n个}$ 至θ₀表示θ中每个分量的收敛。参数所在的条件(14)在实践中不容易检查。一个充分条件是将所有特征值限定在给定的区间内(一,b条)用于数值稳定性(Ingrassia和Rocci，2007年).

接下来，我们证明了在温和的条件下(9)是一致的。它为在矩阵正规混合模型中使用正则化提供了理论依据。我们定义了一个新的参数空间 ${\bar{Ψ}}^{{d日}_{1}, {d日}_{2}}$ 作为

{\bar{Ψ}}^{{d日}_{1}, {d日}_{2}} = {(π_{1}, \dots, π_{K（K）}, {M（M）}_{1}, \dots, {M（M）}_{K（K）}, {V（V）}_{1} \otimes {U型}_{1}, \dots, {V（V）}_{K（K）} \otimes {U型}_{K（K）}) \in Ψ^{{d日}_{1}, {d日}_{2}} : \frac{σ_{我} ({U型}_{小时})}{σ_{我} ({V（V）}_{小时})} = c_{小时} 的 我 = 1, \dots, 最小值 {第页, 第页}, 和 小时 = 1, \dots, K（K）},

(15)

哪里σ_我(U型_小时)表示我矩阵的第th特征值U型_小时和c_小时是每一个的正常数小时= 1, …,K（K）.

定理2。

让Y₁, … ,Y（Y）_n个 是混合矩阵正态分布的随机样本(5).让 ${\hat{Θ}}_{λ}$ 成为受惩罚可能性的最大化者(9)其中惩罚函数的形式为ℓ₁, ℓ₂或核规范。假设真参数值 $Θ_{0} \in {\bar{Ψ}}^{{d日}_{1}, {d日}_{2}}$ ,和λ→ 0作为n→ ∞,然后 $‖ {\hat{Θ}}_{λ} - Θ_{0} ‖ = {o（o）}_{第页} (1)$ .

在这里， $‖ {\hat{Θ}}_{λ} - Θ_{0} ‖ = {\sum^{}}_{我 = 1}^{K（K）} {| {\hat{π}}_{我} - π_{我 0} | + {‖ {\hat{M（M）}}_{我} - {M（M）}_{我 0} ‖}_{F类} + {‖ {\hat{V（V）}}_{我} \otimes {\hat{U型}}_{我} - {V（V）}_{我 0} \otimes {U型}_{我 0} ‖}_{F类}}$ .可以使用范和李（2001）在更强的条件下。然而，这里的重点是介绍新方法及其在图像聚类中的应用；因此，我们选择在本文中不追求这一方向。

4 |. 模拟

我们首先评估提出的CVPL标准是否能够在不同场景下识别正确的簇数。所有模拟结果均基于200次蒙特卡洛复制；报告的CVPL和ARI值是这些复制的平均值。在方案一和方案二中，我们从两个比例相等的矩阵正态分布以及十字和矩形的平均结构的混合物中生成数据，如所示Web图1在这两种情况下，行和列协方差矩阵都遵循自回归设置，其中 $覆盖（cov） {{Y（Y）}_{{k个}_{1}, 我_{1}}, {Y（Y）}_{{k个}_{2}, 我_{2}}} = {0.9}^{| {k个}_{1} - {k个}_{2} | + | 我_{1} - 我_{2} |}$ , 1 ≤k个_我≤第页, 1 ≤我_我≤第页。在场景I中，我们设置样本大小n个=500和图像大小第页=第页= 60. 在场景II中，我们让n个＝100，第页=第页=80。

我们将建议的方法应用于ℓ₁,Ş₂和核规范处罚，并总结结果Web表1如前所述，当λ=0，我们的模型与下面提出的无正则化的混合模型等价维罗利（2011）可以看出，所提出的方法能够选择真实的簇数(K（K）=2）在大多数情况下，基于CVPL的最大值。在三种正则化方法中，通过比较CVPL值，核范数优于ℓ₁和ℓ₂规范和ℓ₂normal的表现略好于ℓ₁规范。这是意料之中的，因为真正的平均值结构具有较低的秩，但不是入门级稀疏结构。

为了获得真实的信噪比（SNR），在场景III中，我们使用了来自实际数据示例的估计参数第5.2节（两个集群下的Theta带时频分析，用于4至8 Hz之间的Theta活动振荡）。我们设定了样本量n个=500，并采用建议的方法。Web表2显示CVPL值有利于两个簇，这与我们的设置一致。在场景IV中，我们通过向场景II中的原始图像添加白噪声来生成合成信号，从而使协方差不再可分离。如所示Web表2，所提出的方法可以在大多数情况下识别正确数量的集群，并且ℓ₂正规化似乎是一个恰当的选择。与模拟场景I和II相比，鉴于噪声水平增加，场景III和IV中的聚类更具挑战性。因此，图像中的潜在稀疏结构变得更加难以检测，这也是为什么我们没有观察到CVPL值在λ=0（无正规化）和λ≠ 0. 对于场景I到场景IV，我们还基于最高CVPL值，给出了正确检测集群数量（超过200个Monte-Carlo复制）的百分比Web表8。我们发现，在大多数情况下K（K）可以以超过70%的概率进行识别。对于选择概率低于70%的情况，一种可能的解释是，存在另一种“竞争性”K（K）（通常靠近真实K（K）)其CVPL值与真实值的CVPL值非常接近K（K）例如，场景一没有惩罚。

我们还进行了深入的研究，以了解所提出的三种正则化方法在不同稀疏度水平上的表现。在场景V中，我们生成值为0和1的三角形和正方形的平均结构。在场景六中，我们遵循与V相同的设置，只是通过添加高斯噪声将三角形的所有零值更改为接近0的小数值。在场景VII中，我们考虑一个新的分离正方形的平均结构（参见Web图6). 然后，我们应用所提出的方法并计算平均参数估计的均方误差（MSE）。发件人Web表7可以看出，使用适当的正则化惩罚将导致更准确的平均参数估计。在场景V中，值为0和1，ℓ₁具有λ=0.5优于其他惩罚和设置。比较CVPL值后，结果也一致。此外，ℓ₂具有λ=0.2和核规范λ=1分别对场景VI和场景VII最有效。结果证实了直觉ℓ₁规范支持零，ℓ₂范数缩小了较小的值，而核范数推动了低阶结构（例如，场景七）。

然后，我们将拟议方法的性能与其他竞争方法进行比较：K（K）-意味着，双聚类(特纳等。, 2005)、和光谱聚类(Ng公司等。, 2002). 我们使用与场景I至场景VII相同的设置生成信号。为了评估性能，我们计算了调整后的随机指数（ARI）(米利根和库珀，1986年)将聚类结果与潜在事实进行比较。ARI是一个最大值为1的数字；它测量两个集群解决方案之间的一致性（即使两个集群方案具有不同数量的集群）。此外，我们考虑了所有方法的预测精度。Web表3到6总结了场景I至IV的结果。所提出的方法在所有场景下都优于所有三种竞争方法，这可能是因为我们的方法考虑了信号中潜在的稀疏结构。例如，在Web表4，提出的方法优于所有三种竞争方法，而双聚类的性能略优于K（K）-在场景IV中，即使在数据生成过程中违反了可分离协方差假设，我们的方法仍然比竞争方法具有更好的预测性能。

值得一提的是，我们在模拟中的主要目标是证明CVPL在以有效的方式指导调谐参数的选择和对图像聚类的正则化效应方面是有用的。如果有人对基于最高CVPL值选择“最佳”模型更感兴趣。理想情况下，可以通过对三个调谐参数（正则化范数、，λ和簇数）。在我们的模拟研究中，我们首先搜索两个值的广泛范围λ（例如{0，.2，.5，.8,1,1.5,2,5,10}）和K（K）（例如，1到6），然后缩小到一个较小的范围，以便于计算。

5 |. 大鼠神经生物学实验中气味记忆数据的分析

在本节中，我们使用提出的方法分析从非空间事件的内存编码实验中获得的LFP数据集(艾伦等。, 2016;胡等。, 2020). 在那次实验中，老鼠被训练在实验期间识别一系列五种气味。在大多数情况下，这五种气味的顺序相同（“按顺序“气味”），但有一些违规行为（“不合顺序的“气味）。例如，气味序列ABCDE公司是“按顺序”的ABBDE公司这是一种“外排”的气味。老鼠被要求用鼻子戳住端口，以正确识别气味是“内排”还是“外排“。在整个实验过程中，基于12个显示任务关键型单细胞活性的微电极收集了峰电位和LFP数据。LFP数据集包含247个试验，采样率为1000 HzT型=2000个时间点。Web图3给出了12只河豚LFP信号的快照。聚类分析对本研究很有用，因为它可以揭示LFP信号中的潜在模式信息，并深入了解它们与不同气味和输入/输出序列的关系。

5.1 |. 图像聚类的时域分析

作为第一步，我们将重点放在时域上，以研究原始多微电极信号与“序列内”或“序列外”模式之间的关联。我们在所有247个试验中对原始LFP信号实施了所提出的方法。表1总结了不同簇数和惩罚之间的CVPL值。基于最高CVPL值，我们的方法为所有惩罚准则选择了两个聚类，这两个聚类对应于输入和输出序列。此外，我们的方法与K（K）-意味着，这表明所提出的方法在检测与“输入”或“输出”序列相关的潜在结构方面具有期望的性能。

表1

气味记忆研究（时域）：针对不同数量的聚类和惩罚，交叉验证惩罚可能性（CVPL）和调整后的随机指数（ARI）值

		CVPL公司				急性呼吸道感染
处罚	λ	K（K）=1	K（K）= 2	K（K）=3	K（K）= 4	我们的方法	K（K）-手段	双星簇	光谱
	0	1.243	1.290^一	1.285	1.281	0.768
ℓ ₁	0.5	1.222	1.253^一	1.253^一	1.246	0.786	0.499	0.583	0.512
	1	1.217	1.243^一	1.206	1.204	0.768
	1.5	1.208	1.249^一	1.234	1.218	0.780
ℓ ₂	0.5	1.264	1.302^一	1.107	1.240	0.768
	1	1.267	1.301^一	1.027	1.202	0.774	0.499	0.583	0.512
	1.5	1.263	1.298^一	1.189	1.235	0.756
核能	0.5	1.253	1.309^一	1.299	1.274	0.756
	1	1.252	1.299^一	1.287	1.277	0.733	0.499	0.583	0.512
	1.5	1.249	1.290^一	1.286	1.214	0.711

在单独的窗口中打开

^美国每行的最高CVPL值（×10⁵).

作为进一步的研究，研究人员还希望了解LFP信号如何与实验中老鼠正确识别气味序列的能力相关。由于序列外试验的样本量较小，我们只关注序列内试验。换句话说，我们关注实验的“敏感性”（真阳性率）。Web表9总结了CVPL和ARI值。基于最高CVPL值，使用K（K）=2簇为首选。基于ARI值ℓ₁核规范规范化的表现优于ℓ₂范数，这可能表明信号中可能存在低秩和稀疏平均结构。对于簇的数量，K（K）=2是首选值，因为其CVPL值通常（10倍中的6倍）高于通过其他选择获得的值K（K）在使用不同的正则化规范下。与K（K）-手段和其他竞争方法。

5.2 |. 时频聚类分析

接下来，我们从时间-频率的角度研究潜在结构。艾伦等。(2016)表明两个特定的振荡带（θ：4到8 Hz和慢γ：20到40 Hz）产生强大的功率，并在检测输入/输出序列中发挥重要作用。Web图3表明时频图表明，低频θ带获得的功率远大于慢γ带。我们将该方法分别应用于θ和慢伽马谱带。表2给出了与θ带真实气味序列相比的聚类数和ARI结果。基于CVPL值，我们发现如果不使用惩罚，则首选三簇模型，并且在采用惩罚的模型中，具有核范数正则化的五簇模型是最佳的。然后我们深入研究了这两个模型，发现与气味A、B和D相关的时频信号在视觉上非常相似。这在某种程度上解释了两种方法之间的不同结果，并证明了我们提出的模型的明显优势。通过添加核范数惩罚，该方法能够在聚类过程中纳入更多信息，从而产生更多的聚类。ARI值也略有下降，有利于核规范规范化。然后我们研究了其他方法的结果，与之前的观察结果类似，气味C和E的结果最一致，而气味A、B和D导致不同的分组。我们的方法提供了一些证据，表明低频带（θ）和气味序列之间的联系。我们还将我们的方法应用于慢伽马波段，并获得了类似的结果（见支持信息，Web表10).

表2

气味记忆研究（时频域）：基于θ谱的“顺序”试验中获得的CVPL和ARI值

		CVPL公司				急性呼吸道感染
处罚	λ	K（K）=1	K（K）= 2	K（K）= 3	K（K）= 4	K（K）= 5	我们的	K（K）-手段	双星簇	光谱
	0	9.432	11.001	11.300^一	11.198	11.172	0.712
ℓ ₁	0.5	7.285	8.516	8.975	8.849	8.997^一	0.692	0.674	0.681	0.678
	1	7.345	8.650	8.632	8.725	8.745^一	0.703
	1.5	7.324	8.571	8.705^一	8.556	8.701	0.709
ℓ ₂	0.5	7.732	8.965^一	8.881	8.671	7.277	0.693
	1	7.634	8.719^一	8.388	8.544	7.616	0.686	0.674	0.681	0.678
	1.5	7.534	8.650	8.632	8.825^一	8.745	0.682
核能	0.5	7.430	9.034	9.196	9.183	9.259^一	0.707
	1	7.243	9.013	9.166	9.255	9.263^一	0.714	0.674	0.681	0.678
	1.5	7.143	8.571	9.040^一	8.995	8.969	0.712

在单独的窗口中打开

^美国每行的最大值（×10^三).

6 |. 大鼠中风数据分析

在本节中，我们将提出的方法应用于另一个LFP数据集，该数据集来自加州大学欧文分校弗罗斯蒂格实验室进行的大鼠中风实验，其中记录了中风前后的LFP。有32个微电极，每个大鼠植入四层。我们在这里考虑的数据是根据中风前后5分钟的信号收集的。采样率为1000 Hz，每个历元长度为1秒。这个实验中有趣的一个科学问题是识别中风前后变化的潜在模式。中的初步时频分析图1（底部4）显示了两个微电极的对数功率谱。这些结果是通过分别平均中风前后的所有时间得到的。大多数微电极的LFP在每个时代内都表现平稳，中风前后的LFP差异似乎很低。然而，对于某些微电极（尤其是10），存在一些不可忽视的动力学，并且在冲程前和冲程后存在明显差异。这些发现表明，在进行聚类分析以确定中风前后变化的潜在模式时，对所有通道进行平均或矢量化可能不是最佳方法。

我们还研究了中风前后所有32个微电极的动力学。图2显示了跨越微电极的β和慢γ频带的时频图。对数功率谱是通过对历元进行平均得到的。通过比较中风前后的曲线图，我们观察到两个频带的微电极之间存在强烈的相关性。这突出了在混合模型中引入正则化的重要性。此外，通过比较中风前后的曲线图，有明显的局部差异迹象。如果在进行聚类分析时不考虑矩阵结构而对原始信号进行天真的矢量化，那么这种差异很容易被忽略。

保存图片、插图等的外部文件。对象名为nihms-1635244-f0003.jpg

在单独的窗口中打开

图2

大鼠卒中研究：卒中前后所有微电极之间特定频带的时频图。这张图以彩色出现在本文的电子版中，任何提到彩色的都是指该版本

我们将所提出的方法应用于中风前后所有时期的时频图像。表3显示了不同数量的簇和惩罚标准的CVPL值。除了一个例外，所有场景都建议两个聚类，分别对应“正常”和“中风”两种状态。我们还将聚类结果与真值状态指数进行了比较，并总结了所提方法的ARI和K（K）-平均值表4。与K（K）-手段。通过正则化，我们的方法产生的ARI值增加了80%（相比之下λ=0，带λ≠ 0). 特别是，当核定额惩罚与λ= 2. 对于β波段也获得了类似的结果。这些发现与初步分析中的发现一致。

表3

大鼠中风研究：β和慢γ波段对数功率谱的CVPL值

		CVPL（慢伽马射线）					CVPL（β）
处罚	λ	K（K）=1	K（K）= 2	K（K）= 3	K（K）= 4	K（K）= 5	K（K）=1	K（K）= 2	K（K）= 3	K（K）= 4	K（K）= 5
	0	2.781	2.941	2.964^一	2.764	2.822	4.538	4.645^一	4.627	4.526	4.598
ℓ ₁	1	2.364	2.472^一	2.031	1.513	0.421	3.832	3.980^一	3.268	3.594	1.676
	2	2.043	2.106^一	1.370	1.288	0.621	4.032	4.167^一	3.373	3.227	3.277
ℓ ₂	0.5	2.564	2.688^一	2.474	2.306	2.184	4.129	4.245^一	4.179	4.036	3.424
	1	2.354	2.484^一	2.188	1.787	1.895	3.845	4.063^一	3.557	3.429	3.329
	2	2.245	2.338^一	2.163	1.539	1.733	3.743	4.024^一	3.699	2.972	3.206
核能	0.5	2.745	2.806^一	2.627	2.502	2.303	4.356	4.464^一	4.299	4.130	3.963
	1	2.434	2.556^一	2.362	1.946	1.720	4.022	4.191^一	3.977	3.618	3.371
	2	2.634	2.748^一	1.689	1.257	0.684	3.543	3.687^一	3.274	2.795	2.262

在单独的窗口中打开

^美国不同频带上每行的最高CVPL值（×10⁴).

表4

大鼠中风研究：慢伽马和β谱带的ARI与“中风”的关系

		ARI（慢伽马）				ARI（β）
处罚	λ	我们的方法	K（K）-手段	双星簇	光谱	我们的方法	K（K）-手段	双星簇	光谱
	0	0.507				0.887
ℓ ₁	0.5	0.981	0.751	0.783	0.764	0.942	0.715	0.741	0.724
	1	0.961				0.914
	2	0.951				0.861
ℓ ₂	0.5	0.951				0.941
	1	0.951	0.751	0.783	0.764	0.878	0.715	0.741	0.724
	2	0.961				0.787
核能	0.5	0.951				0.941
	1	0.960	0.751	0.783	0.764	0.942	0.715	0.741	0.724
	2	1				0.951

在单独的窗口中打开

7 |。总结性意见

本文提出了一种正则化概率聚类框架来分析矩阵数据。与现有的聚类方法相比，如K（K）-也就是说，优点如下：（a）通过直接处理矩阵数据，我们能够同时捕获行和列的相关性；（b）该框架能够揭示信号和图像固有的自然稀疏性；（c）通过使用CVPL，该方法能够找到适用于不同类型信号结构和稀疏程度的最优正则化方法和调谐参数；（d）提出的方法是基于理论基础的；提供直观的解释性；并且具有较低的计算成本（通过并行计算），因此适用于大数据集。

为了在图像中合并不同的稀疏结构，我们在此框架中提出了三种正则化方法，ℓ₂范数惩罚较大的矩阵元素值，ℓ₁范数以识别非零为目标，核范数鼓励低阶近似。在实践中，可以根据CVPL准则或关于图像结构的先验知识选择最佳正则化范数。我们的仿真结果证实了该方法在不同场景下使用适当的正则化方法在图像结构恢复和估计精度提高（如MSE）方面的优异性能。

尽管本文提供了一些有希望的结果，但在分析矩阵数据时仍然存在许多公开的问题。例如，在当前的工作中，选择簇的数量依赖于一些预先指定的度量（CVPL）。作为由引入的框架的扩展维罗利（2011）可以在聚类分析中引入贝叶斯框架，对聚类数进行贝叶斯推理。与在回归模型中使用弹性网类似，考虑两者的组合将是未来的兴趣所在ℓ₁和ℓ₂我们模型中的惩罚。在提出的混合模型中加入正则化协方差矩阵估计也将很有帮助，特别是对于分析大尺寸图像。另一个值得关注的未来工作方向是开发推理程序，用于量化聚类质心和图像信号之间的相似性（及其不确定性），以进行预测。

补充材料

附录

单击此处查看。^{（540K，pdf格式）}

确认

作者感谢主编、副主编和两位审稿人对本文早期版本的建设性和有益评论。沈的研究部分得到西蒙斯基金会奖512620和国家科学基金会拨款DMS 1509023的支持。胡锦涛的努力得到了国家卫生拨款研究所R01AI143886、R01CA219896和CCSG P30 CA013696的部分支持。福廷的研究得到了NIH拨款R01-MH115697、NSF拨款IOS-1150292和BCS-1439267以及白厅基金会拨款2010-05-84的支持。弗罗斯蒂格的研究部分得到了莱杜克基金会的支持。

资金筹措信息

西蒙斯基金会，赠款/奖项编号：512620；国家卫生研究所，拨款/奖励编号：R01AI143886、R01CA219896、CCSG P30 CA013696；NIH，授予/授予编号：R01-MH115697；NSF，授予/奖励编号：IOS-1150292，BCS-1439267；白厅基金会，拨款/奖项编号：2010-05-84

脚注

数据可用性声明

本文的支持信息中提供了支持本文研究结果的数据。

支持信息

第4-6节中引用的Web附录、表格和图表可在Wiley Online Library的生物统计学网站上获取。计算代码和数据集可从以下网站获得（有关详细说明，请参阅支持信息第4节）：https://github.com/XuGaoUCI/Regularized_Matrix_Clustering网站.

参考文献

Allen TA、Salz DM、McKenzie S和Fortin NJ（2016）ca1神经元的非空间序列编码.神经科学杂志,36, 1547–1563.[PMC免费文章][公共医学][谷歌学者]
陈杰（2017）混合模型下MLE的一致性.统计科学,32, 47–63.[谷歌学者]
Chi EC、Allen GI和Baraniuk RG（2017年）凸双簇.生物计量学,73, 10–19. [公共医学][谷歌学者]
Cressie N（2015）空间数据统计纽约：John Wiley&Sons。[谷歌学者]
Dawid AP（1981）一些矩阵变量分布理论：符号考虑和贝叶斯应用.生物特征,68, 265–274.[谷歌学者]
Dempster AP、Laird NM和Rubin DB（1977年）通过EM算法从不完整数据中获得最大似然.英国皇家统计学会杂志。B系列（方法学）,39, 1–38.[谷歌学者]
Donoho DL、Gavish M和Johnstone IM（2018年）尖峰协方差模型中特征值的最优收缩.统计年刊,46, 1742–1778.[PMC免费文章][公共医学][谷歌学者]
Dutilleul P（1999）矩阵正态分布的MLE算法.统计计算与模拟杂志,64, 105–123.[谷歌学者]
Euan C、Ombao H和Ortega J（2018）分层谱合并算法：一种新的时间序列聚类.分类杂志,35, 71–99.[谷歌学者]
Euan C、Sun Y和Ombao H（2019）基于相干性的脑连接可视化时间序列聚类.应用统计年鉴,13, 990–115.[谷歌学者]
范J、李锐（2001）基于非冲突惩罚似然的变量选择及其oracle性质.美国统计协会杂志,96, 1348–1360.[谷歌学者]
高X、沈伟、胡杰、福田N和Ombao H（2019）基于正则化矩阵数据聚类的局部场势建模2019年第9届国际IEEE/EMBS神经工程会议（NER），第597-602页。电气与电子工程师协会。[谷歌学者]
Gao X、Shen W、Shahbaba B、Fortin N和Ombao H（2020）演化状态空间模型及其在局部场电位时频分析中的应用.中国统计局,30, 1561–1582.[PMC免费文章][公共医学][谷歌学者]
Gao X、Shen W、Ting C-M、Cramer SC、Srinivasan R和Ombao H（2019）用copula高斯图形模型估计大脑连通性2019年IEEE第16届国际生物医学成像研讨会（ISBI 2019），第108–112页。电气与电子工程师协会。[谷歌学者]
Genton MG（2007）时空协方差矩阵的可分离逼近.环境计量学,18, 681–695.[谷歌学者]
Green PJ（1990）EM在惩罚似然估计中的应用.英国皇家统计学会杂志。B系列（方法学）,52, 443–452.[谷歌学者]
Gupta AK和Nagar DK（1999年）矩阵变量分布，音量104佛罗里达州博卡拉顿：CRC出版社。[谷歌学者]
哈斯TC（1995）应用于硫酸盐湿沉降的时空过程局部预测.美国统计协会杂志,90, 1189–1199.[谷歌学者]
Hu L、Fortin NJ和Ombao H（2019）高维多通道大脑信号建模.生物科学统计,11, 91–126.[谷歌学者]
Hu L、Guindani M、Fortin N和Ombao H（2020年）多试验脑信号差分连接的层次贝叶斯模型.计量经济与统计,15, 117–135.[PMC免费文章][公共医学][谷歌学者]
Ingrassia S和Rocci R（2007）多元高斯有限混合的约束单调EM算法。计算统计学& 数据分析,51, 5339–5351.[谷歌学者]
国王RS（2015）聚类分析与数据挖掘：简介弗吉尼亚：Stylus Publishing，LLC。[谷歌学者]
兰格·K（1995）与EM算法局部等价的梯度算法.英国皇家统计学会杂志。B系列（方法学）,57, 425–437.[谷歌学者]
Lange K、Hunter DR和Yang I（2000）使用替代目标函数的优化传递.计算与图形统计杂志,9, 1–20.[谷歌学者]
Milligan GW和Cooper MC（1986年）层次聚类分析外部标准的可比性研究.多元行为研究,21, 441–458. [公共医学][谷歌学者]
Ng AY、Jordan MI和Weiss Y（2002）关于谱聚类：分析和算法《神经信息处理系统进展会议》，第849-856页。[谷歌学者]
Pan W和Shen X（2007）基于惩罚模型的聚类及其在变量选择中的应用.机器学习研究杂志,8, 1145–1164.[谷歌学者]
Reiss PT和Ogden RT（2007）函数主成分回归与函数偏最小二乘.美国统计协会杂志,102, 984–996.[谷歌学者]
Smyth P（2000）基于交叉验证似然的概率聚类模型选择.统计与计算,10, 63–72.[谷歌学者]
Turner H、Bailey T和Krzanowski W（2005）通过系统性能测试证明了微阵列数据的改进双聚类.计算统计与数据分析,48，235–254。[谷歌学者]
维罗利C（2011）用于三向数据分类的矩阵正态分布的有限混合.统计与计算,21，511–522。[谷歌学者]
维罗利C（2011）基于模型的三向数据结构聚类.贝叶斯分析,6, 573–602.[谷歌学者]
王X、朱赫和ADNI（2017）基于全变分的广义标度图像回归模型.美国统计协会杂志,112, 1156–1168.[PMC免费文章][公共医学][谷歌学者]
王毅、丁C-M、高X和Ombao H（2019）低维嵌入脑信号的探索性分析2019年第9届国际IEEE/EMBS神经工程会议（NER），第997-1002页。电气与电子工程师协会。[谷歌学者]
于毅（2012）单调超松弛EM算法.计算与图形统计杂志,21, 518–537.[谷歌学者]
周浩、李磊（2014）正则化矩阵回归.英国皇家统计学会杂志。B系列（方法学）,76, 463–483.[PMC免费文章][公共医学][谷歌学者]

正则化矩阵数据聚类及其在图像分析中的应用

徐高

沈伟宁

张立文

胡建华

诺伯特·J·福廷

罗恩·弗罗斯蒂格

埃尔南多·奥姆巴奥

关联数据

摘要

1 |. 简介

2|。 方法

2.1 |. 矩阵正态分布

2.2 |. 矩阵正态混合模型

2.3 |. 惩罚矩阵正态混合模型

I.（初始化）

二、。（E步骤）

三、 （M步）

四、 （停止标准）

四、 （选择簇数）

3|。 理论

定理1。

定理2。

4 |. 模拟

5 |. 大鼠神经生物学实验中气味记忆数据的分析

5.1 |. 图像聚类的时域分析

表1

5.2 |. 时频聚类分析

表2

6 |. 大鼠中风数据分析

表3

表4

7 |。 总结性意见

补充材料

附录

确认

脚注

参考文献

2|。方法

三、（M步）

四、（停止标准）

四、（选择簇数）

3|。理论

7 |。总结性意见