«上一篇 下一篇»
  计算机工程  2022, Vol. 48 Issue (4): 197-205  DOI: 10.19678/j.issn.1000-3428.0060601
0

引用本文  

陈璐瑶, 刘奇龙, 许云霞, 等. 基于超图正则化非负Tucker分解的图像聚类算法[J]. 计算机工程, 2022, 48(4), 197-205. DOI: 10.19678/j.issn.1000-3428.0060601.
CHEN Luyao, LIU Qilong, XU Yunxia, et al. Image Clustering Algorithm Based on Hypergraph Regularized Nonnegative Tucker Decomposition[J]. Computer Engineering, 2022, 48(4), 197-205. DOI: 10.19678/j.issn.1000-3428.0060601.

基金项目

国家自然科学基金(12061025);贵州省科学技术基金项目(黔科合基础[2020]1Z002);贵州省教育厅自然科学研究项目(黔教合KY字[2020]298)

通信作者

刘奇龙(通信作者),副教授、博士

作者简介

陈璐瑶(1995—),女,硕士研究生,主研方向为张量分解;
许云霞,讲师、博士研究生;
陈震,教授、博士生导师

文章历史

收稿日期:2021-01-15
修回日期:2021-03-24
基于超图正则化非负Tucker分解的图像聚类算法
陈璐瑶 , 刘奇龙 , 许云霞 , 陈震     
贵州师范大学 数学科学学院, 贵阳 550025
摘要:针对非负张量分解应用于图像聚类时忽略了高维数据内部几何结构的问题,在经典的张量非负Tucker分解的基础上,添加超图正则项以尽可能多地保留原始数据的内在几何结构信息,提出一种基于超图正则化非负Tucker分解模型HGNTD。通过构造超图刻画数据内部样本间的高阶关系,提高几何结构描述的准确性,针对超图正则化非负张量分解模型,基于交替非负最小二乘法,设计快速有效的超图正则化非负Tucker分解算法求解所给模型,证明算法在非负的条件下是收敛的,最终将算法应用于图像聚类。在Yale和COIL两个常用公开数据集上的实验结果表明,相对于k-means、非负矩阵分解、图正则化非负矩阵分解、非负Tucker分解和图正则化非负Tucker分解等算法,超图正则化非负Tucker分解算法聚类准确度提升了8.6%~11.4%,归一化互信息提升了2.0%~7.5%,具有更好的聚类效果。
关键词非负张量分解    Tucker分解    超图学习    交替非负最小二乘法    聚类分析    
Image Clustering Algorithm Based on Hypergraph Regularized Nonnegative Tucker Decomposition
CHEN Luyao , LIU Qilong , XU Yunxia , CHEN Zhen     
School of Mathematical Sciences, Guizhou Normal University, Guiyang 550025, China
Abstract: The internal geometry structure of high-dimensional data is ignored when nonnegative tensor decomposition is applied to image clustering.To solve this problem, we propose a Hypergraph regularized Nonnegative Tucker Decomposition(HGNTD) model by adding a hypergraph regularization term to preserve the intrinsic geometric structure information of original data as much as possible.This method is based on the classical nonnegative Tucker decomposition of tensors.Specifically, relying on the hypergraph regularization non-negative tensor decomposition model and based on the Alternating Nonnegative Least Squares(ANLS), a fast and effective hypergraph regularized nonnegative tucker decomposition algorithm is designed to solve the given model.The new algorithm is proven to converge under nonnegative conditions, thereby, it is applied to image clustering.The experimental results on Yale and COIL-100 data sets show that, compared with k-means, Nonnegative Matrix Factorization(NMF), Graph regularized Nonnegative Matrix Factorization(GNMF), Nonnegative Tucker Decomposition(NTD), and Graph regularized Nonnegative Tucker Decomposition(GNTD) algorithms, the clustering accuracy of the new algorithm is improved by approximately 8.6%~11.4%, and the Normalized Mutual Information(NMI) is improved by approximately 2.0%~7.5%.Therefore, it can be said to have a better clustering effect.
Key words: nonnegative tensor decomposition    Tucker decomposition    hypergraph learning    alternating nonnegative least squares    clustering analysis    

开放科学(资源服务)标志码(OSID):

0 概述

基于矩阵分解的数据处理方法近年来在计算机视觉、模式识别、信号处理、生物医学工程和图像工程等[1-3]研究领域中得到了成功的应用,如奇异值分解、主成分分析、独立成分分析、线性判别分析等方法。这些方法都是将原始矩阵近似分解为多个低秩矩阵的乘积[4]。分解的因子矩阵的元素可能是正数也可能是负数,但在实际应用问题中负值元素往往没有意义或无法解释,如图像数据的像素值不可能是负值,文档的数目及词汇的个数也不可能是负值。鉴于此,文献[5]提出一种新的矩阵分解方法——非负矩阵分解(Nonnegative Matrix Factorization,NMF),与传统的矩阵分解方法相比,非负矩阵分解所得结果均为非负值,这使得非负矩阵分解在实际问题中得到了广泛应用。文献[6]提出已有的非负矩阵分解方法通常是将二维的图像变为一维向量处理,这种做法破坏了图像像素点之间的几何结构关系,在数据降维和特征提取过程中容易导致信息丢失、维数灾难和小样例问题,文献[7]提出张量表示能在一定程度上避免上述问题。张量是矩阵和向量的高阶推广,其中0阶张量是标量、1阶张量是向量、2阶张量是矩阵。张量分解的两种经典模型为PARAFAC分解和Tucker分解[8-10]。PARAFAC分解是将一个N阶张量分解为若干个秩一张量和的形式;Tucker分解是将原张量分解成核张量和因子矩阵乘积的形式。文献[11]提出在非负性约束下,Tucker分解可以像矩阵分解一样提取基于局部特征的表示。因此,非负Tucker分解(Nonnegative Tucker Decomposition,NTD)是处理非负张量数据的一种常用方法。

文献[12-14]研究数据之间的几何结构,流形学习能够发现嵌入在高维数据空间中观测数据的低维本质结构,以取得更好的低维表示[15]。目前已有一些研究将流形学习理论与非负张量分解理论结合。例如,文献[16]考虑到图像空间的流形结构,将局部线性嵌入正则化项添加到非负平行因子分析(NPARAFAC)模型,研究了邻域保持的非负张量分解方法。另外,文献[17]提出拉普拉斯正则化非负张量分解方法(Laplace Regularized Nonnegative Tensor Decomposition,LRNTD)并用于图像聚类。在每种模式下,NPARAFAC通常比NTD需要更多的分量[18],导致NPARAFAC的数据表示效果不如NTD。传统的Tucker分解只涉及属性信息,而没有考虑图像之间的成对相似信息。文献[19]考虑属性信息与图像间局部表示的相似程度,提出一种图拉普拉斯Tucker分解(Graph Laplace Tucker Decomposition,GLTD)方法。文献[20]将流形正则化项引入非负Tucker分解中,提出一种流形正则化非负Tucker分解(Manifold Regularized Nonnegative Tucker Decomposition,MRNTD)方法。文献[21]在低维表示的前提下,为了同时保持内部的多线性结构和几何信息,将图正则化引入到非负Tucker分解中,提出一种图正则化的非负Tucker分解(Graph regularized Nonnegative Tucker Decomposition,GNTD)方法。

然而上述研究只考虑到了样本间的成对关系,忽略了样本间的高阶关系[22-24],文献[25-27]则利用超图可以解决这一不足。本文结合超图学习[22]和NTD算法,提出一种基于超图正则化非负Tucker分解算法(HyperGraph Regularized Nonnegative Tucker Decomposition,HGNTD)。利用K-最邻近(K-Nearest Neighbor,KNN)方法建立超图,将超图正则项添加到非负Tucker分解中,并基于交替非负最小二乘法(Alternating Nonnegative Least Squares,ANLS),设计超图正则化非负Tucker分解算法(HGNTD)求解模型。

1 相关工作 1.1 非负Tucker分解

张量分解有两种经典模型:PARAFAC分解[8]和Tucker分解[9]。这两种模型都是矩阵奇异值分解的高阶推广,也可看作是主成分分析(PCA)的高阶形式,且前者是后者的一种特例。

非负Tucker分解(NTD)[18]是把一个$ N $阶非负张量y$ \in {\mathbb{R}}_{+}^{{I}_{1}\times {I}_{2}\times \cdots \times {I}_{N}} $分解为核张量与因子矩阵的乘积,即:

$ \boldsymbol{y}\approx \boldsymbol{S}{\times }_{n\in \mathcal{N}}{\boldsymbol{A}}^{\left(n\right)} $

其中:$ \boldsymbol{S}\in {R}_{+}^{{J}_{1}\times {J}_{2}\times \cdots \times {J}_{N}} $为核张量;$ {\boldsymbol{A}}^{\left(n\right)}\in {\mathbb{R}}_{+}^{{I}_{n}\times {J}_{n}} $为因子矩阵$ ({J}_{n}\le {I}_{n}, \mathcal{N}=\{\mathrm{1, 2}, \cdots , N\left\}\right) $

非负Tucker分解可转化为如下最优化问题:

$ \begin{array}{l}\underset{\boldsymbol{s}, {\boldsymbol{A}}^{\left(n\right)}, n\in \mathcal{N}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\boldsymbol{y}-\boldsymbol{s}{\times }_{n\in \mathcal{N}}{\boldsymbol{A}}^{\left(n\right)}‖}_{F}^{2}\\ \begin{array}{cc}\mathrm{s}.\mathrm{t}.\boldsymbol{S}\ge 0, {\boldsymbol{A}}^{\left(n\right)}\ge 0, n\in \mathcal{N}& \end{array}\end{array} $

按张量模-$ N $展开,非负Tucker分解可表示为矩阵形式:

$ \begin{array}{l}\underset{{\boldsymbol{S}}_{\left(N\right)}, {\boldsymbol{A}}^{\left(N\right)}, {\boldsymbol{H}}_{N}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\boldsymbol{Y}}_{\left(N\right)}-{\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{S}}_{\left(N\right)}{\boldsymbol{H}}_{N}^{\mathrm{T}}‖}_{F}^{2}\\ \mathrm{s}.\mathrm{t}.{\boldsymbol{S}}_{\left(N\right)}\ge 0, {\boldsymbol{A}}^{\left(N\right)}\ge 0, {\boldsymbol{H}}_{N}\ge 0\end{array} $ (1)

其中:$ {\boldsymbol{Y}}_{\left(N\right)}\in {\mathbb{R}}_{+}^{{I}_{N}\times {I}_{1}{I}_{2}\cdots {I}_{N-1}} $$ {\boldsymbol{S}}_{\left(N\right)}\in {\mathbb{R}}_{+}^{{J}_{N}\times {J}_{1}{J}_{2}\cdots {J}_{N-1}} $分别为张量YS的模-$ N $展开矩阵[10]$ {\boldsymbol{H}}_{N}={\boldsymbol{A}}^{(N-1)} \otimes {\boldsymbol{A}}^{(N-2)}\cdots \otimes {\boldsymbol{A}}^{\left(1\right)} $$ \otimes $表示矩阵Kronecker积。将模-$ N $展开,$ {\boldsymbol{A}}^{\left(N\right)} $是与基矩阵$ {\boldsymbol{S}}_{\left(N\right)}{\boldsymbol{H}}_{N}^{\mathrm{T}} $相对应的系数矩阵。

NTD可以看作是NMF的多线性推广,NTD的矩阵形式(式(1))和NMF之间的一个最大的区别是前者的基矩阵是核张量与因子矩阵从模-1到模-$ \left(N-1\right) $的乘积,这有利于改善基矩阵的稀疏性。

1.2 超图学习

超图是普通图的推广,在普通图中,一个数据对象表示一个顶点,两个顶点之间的边用来描绘顶点之间的关系。超图中的超边连接3个或者更多的顶点可以更好地表示数据的高维信息,超图已经广泛应用于分类、聚类中。下文简单地介绍超图的一些基本理论[22, 26]

假设$ G=\left(V, E, \boldsymbol{W}\right) $,其中$ V=\left\{{v}_{1}, {v}_{2}, \cdots , {v}_{N}\right\} $是顶点集,集合$ V $中的每一个元素称为顶点;$ E=\left\{{e}_{1}, {e}_{2}, \cdots , {e}_{M}\right\} $是超边集,集合$ E $中的每一个元素称为超边,每条超边是顶点集$ V $的子集;$ \boldsymbol{W} $为超图的权重矩阵,是由超边的权重构成的对角矩阵,其中每条超边$ e $的权重是事先赋的一个正值$ w\left(e\right) $。如果以下两个条件成立:$ {e}_{j}\notin \mathrm{\varnothing }, \forall j\in \mathrm{1, 2}, \cdots , M $$ {e}_{1}\bigcup {e}_{2}\bigcup \cdots \bigcup {e}_{M}=V $,则$ G $是定义在集合$ V $上的超图。

超图可使用一个$ \left|V\right|\times \left|E\right| $的关联矩阵$ \boldsymbol{H} $来表示:

$ \boldsymbol{H}\left(v, e\right)=\left\{\begin{array}{c}1, \begin{array}{cc}v\in e& \end{array}\\ 0, \begin{array}{cc}v\notin e& \end{array}\end{array}\right. $ (2)

下文用几何图形解释超图,如图 1所示,其对应的顶点与超边的关联矩阵$ \boldsymbol{H} $表 1所示。

Download:
图 1 超图G的几何示意图 Fig. 1 Geometric schematic diagram of hypergraph G
下载CSV 表 1 图 1对应的关联矩阵H Table 1 The correlation matrix H corresponding to Fig. 1

图 1给出一个超图$ G=\left(V, E\right) $,实心节点表示为顶点,由椭圆虚线标记的集合表示超边。表 1表示为其对应的关联矩阵$ \boldsymbol{H} $,其中,$ V=\left\{{v}_{1}, {v}_{2}, \cdots , {v}_{7}\right\} $$ E=\left\{{e}_{1}, {e}_{2}, {e}_{3}\right\} $$ {e}_{1}=\left\{{v}_{1}, {v}_{2}, {v}_{3}\right\} $$ {e}_{2}=\left\{{v}_{3}, {v}_{4}, {v}_{5}\right\} $$ {e}_{3}=\left\{{v}_{5}, {v}_{6}, {v}_{7}\right\} $

对于顶点$ v\in V $$ v $所属的超边的权重之和称为$ v $的阶(或度),记为$ d\left(v\right)=\sum \limits_{\left\{\boldsymbol{e}\in E\left|v\in \boldsymbol{e}\right.\right\}}w\left(e\right) $;对于超边$ e $$ \in E $$ e $所包含的顶点数称为$ e $的阶(或度),记为$ \delta \left(e\right)=\left|e\right| $。因此:

$ d\left(v\right)=\sum \limits_{j=1}^{M}w({e}_{j})\boldsymbol{H}\left(v, {e}_{j}\right) \text{,} \delta \left(e\right)=\sum \limits_{j=1}^{N}\boldsymbol{H}\left({v}_{j}, e\right) $

分别用$ {\boldsymbol{D}}_{v}\mathrm{、}{\boldsymbol{D}}_{e}\mathrm{、}\boldsymbol{W} $表示顶点的度、超边的度和超边的权重所形成的对角矩阵,可分别表示为:

$ {\boldsymbol{D}}_{v}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[d\left({v}_{1}\right), d\left({v}_{2}\right), \cdots , d\left({v}_{N}\right)\right] $
$ {\boldsymbol{D}}_{e}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\delta \left({e}_{1}\right), \delta \left({e}_{2}\right), \cdots , \delta \left({e}_{M}\right)\right] $
$ \boldsymbol{W}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[w\left({e}_{1}\right), w\left({e}_{2}\right), \cdots , w\left({e}_{M}\right)\right] $

根据文献[22]定义,非标准化超图拉普拉斯矩阵为:

$ {\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}={\boldsymbol{D}}_{v}-\boldsymbol{B} $ (3)

其中:$ \boldsymbol{B}=\boldsymbol{H}\boldsymbol{W}{\boldsymbol{D}}_{e}^{-1}{\boldsymbol{H}}^{\mathrm{T}} $

根据超边的定义,允许包含任意数量的顶点。为了方便起见,本文采用在每条超边中包含相同顶点个数的一致超图。

2 超图正则非负Tucker分解

本文在NTD和超图学习基础上,提出张量数据超图的构造方法,并将超图的正则化项引入NTD目标函数中,最后给出求解HGNTD模型的有效算法。

2.1 超图构造

为了构造超图,可以把数据集中的每个数据点看作超图的顶点,将当前顶点和其最邻近的顶点看作一条超边,通过改变最邻近半径大小定义一组超边。根据超图的不同构造方法可以生成大量超边,使超图更好地描述数据结构。此外,超边权重的选择也很重要,目前最常用的加权方案有3种:0-1加权方案,热核(Heat Kernel)加权方案和点乘加权方案。不同的加权方案适用于不同的应用问题,如点乘权重通常适用于文本匹配问题,而热核加权方案常用于处理图像数据。

本文中采用的是热核加权方案,基于给定的图像空间邻域构造一个超图$ G=\left(V, E, \boldsymbol{W}\right) $。首先,将$ Y=\left\{{\boldsymbol{y}}_{1}, {\boldsymbol{y}}_{2}, \cdots , {\boldsymbol{y}}_{N}\right\} $中的每张图像$ {\boldsymbol{y}}_{i} $看作超图的一个顶点$ {v}_{i} $。然后,以每个顶点$ {v}_{i} $作为空间区域的“质心”节点,并使用K-最邻近(KNN)方法构造相应的超边$ {e}_{i}\in E $,因此,超图$ G $由构造在不同空间区域上的$ N $个超边组成。基于这些超边,通过式(2)获得关联矩阵$ \boldsymbol{H}\in {\mathbb{R}}^{\left|V\right|\times \left|E\right|} $。最后,使用热核加权方案,赋予每条超边一个正权重[28],即:

$ w\left({e}_{i}\right)=\sum \limits_{{\boldsymbol{y}}_{j}\in {e}_{i}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{‖{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}‖}^{2}}{{\sigma }^{2}}\right) $

其中:$ \sigma =\frac{1}{kN}\sum \limits_{i=1}^{N}\sum \limits_{\left\{{\boldsymbol{y}}_{j}\in {\boldsymbol{e}}_{i}\right\}}‖{\boldsymbol{y}}_{j}-{\boldsymbol{y}}_{i}‖ $表示所有顶点之间的平均距离。

根据上述超图的构造方法,每个顶点与至少一条超边连接,并且每条超边与权重相关联。对角权重矩阵$ \boldsymbol{W}\in {\mathbb{R}}^{N\times N} $由下式给出:

$ {\left(\boldsymbol{W}\right)}_{ij}=\left\{\begin{array}{l}w\left({e}_{i}\right), \begin{array}{cc}i=j& \end{array}\\ 0, \begin{array}{ccc}i\ne j& & \end{array}\end{array}\right. $

权重矩阵$ \boldsymbol{W} $的权值与两点间的欧氏距离成反比。两个点之间的欧氏距离越小,它们的相似性越高。因此,使用下式来刻画低维数据的光滑性:

$ C=\sum \limits_{e\in E}\sum \limits_{\left(i, j\right)\in e}\frac{w\left(e\right)}{\delta \left(e\right)}{‖{\boldsymbol{A}}_{i:}^{\left(N\right)}-{\boldsymbol{A}}_{j:}^{\left(N\right)}‖}^{2} $

在上式中,$ {\boldsymbol{A}}_{i:}^{\left(N\right)} $$ {\boldsymbol{A}}^{(N)} $的第$ i $行,可以重新表达如下:

$ C=\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right) $ (4)

其中:$ \mathrm{T}\mathrm{r}\left(\cdot \right) $表示矩阵的迹;$ {\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}} $是式(3)所定义的非标准化超图$ G $的超图拉普拉斯矩阵。

2.2 目标函数

将超图正则化项式(4)加入到原始NTD的目标函数中,得到HGNTD目标函数如下:

$ \begin{array}{l}\underset{\boldsymbol{s}, {\boldsymbol{A}}^{\left(n\right)}, n\in \mathcal{N}}{\mathrm{m}\mathrm{i}\mathrm{n}}{O}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}=\frac{1}{2}{‖\boldsymbol{y}-\boldsymbol{s}{\times }_{n\in \mathcal{N}}{\boldsymbol{A}}^{\left(n\right)}‖}_{F}^{2}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\lambda }{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right)\end{array} $
$ \mathrm{s}.\mathrm{t}.\ge 0, {\boldsymbol{A}}^{\left(n\right)}\ge 0, n\in N $ (5)

其中:$ \lambda $是一个非负参数,用于平衡超图正则化项和重构误差项;$ {\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}} $是式(3)所定义的非标准化超图$ G $的超图拉普拉斯矩阵。

采用拉格朗日乘子法,并考虑模-$ n $展开形式,然后确定基于式(5)的拉格朗日函数:

$ {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}=\frac{1}{2}{‖{\boldsymbol{Y}}_{\left(n\right)}-{\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}{)}^{\mathrm{T}}‖}_{F}^{2}+\\ \frac{\lambda }{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right)+ \mathrm{T}\mathrm{r}\left({\boldsymbol{\varPhi }}_{n}{\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)+\sum \limits_{l=1}^{N}Tr\left({\boldsymbol{\varPsi }}_{l}{{\boldsymbol{A}}^{\left(l\right)}}^{{}^{T}}\right) $

其中:$ \underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}={\boldsymbol{A}}^{\left(N\right)}\otimes \cdots \otimes {\boldsymbol{A}}^{(n+1)}\otimes {\boldsymbol{A}}^{(n-1)}\otimes \cdots \otimes {\boldsymbol{A}}^{\left(1\right)} $$ {\boldsymbol{\varPhi }}_{n} $表示$ {\boldsymbol{S}}_{\left(n\right)} $的拉格朗日乘子矩阵;$ {\boldsymbol{\varPsi }}_{l} $表示$ {\boldsymbol{A}}^{\left(l\right)}\left(l\in \mathcal{N}\right) $的拉格朗日乘子矩阵。

可将目标函数重新写为:

$ \begin{array}{l}{\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}=\frac{1}{2}\mathrm{T}\mathrm{r}\left(\left({\boldsymbol{Y}}_{\left(n\right)}-{\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}{\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right)}^{\mathrm{T}}\right){\left({\boldsymbol{Y}}_{\left(n\right)}-{\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}{\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right)}^{\mathrm{T}}\right)}^{\mathrm{T}}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\lambda }{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right)+\mathrm{T}\mathrm{r}\left({\boldsymbol{\varPhi }}_{n}{\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)+\sum \limits_{l=1}^{N}\mathrm{T}\mathrm{r}\left({\boldsymbol{\varPsi }}_{l}{{\boldsymbol{A}}^{\left(l\right)}}^{{}^{\mathrm{T}}}\right)=\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{Y}}_{\left(n\right)}{{\boldsymbol{Y}}_{\left(n\right)}}^{\mathrm{T}}\right)-\mathrm{T}\mathrm{r}\left({\boldsymbol{Y}}_{\left(n\right)}{\left({\boldsymbol{S}}_{\left(n\right)}{\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right)}^{\mathrm{T}}\right)}^{\mathrm{T}}{{\boldsymbol{A}}^{\left(n\right)}}^{{}^{\mathrm{T}}}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(n\right)}\left({\boldsymbol{S}}_{\left(n\right)}{\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right)}^{\mathrm{T}}\right){\left({\boldsymbol{S}}_{\left(n\right)}{\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right)}^{\mathrm{T}}\right)}^{\mathrm{T}}{{\boldsymbol{A}}^{\left(n\right)}}^{{}^{\mathrm{T}}}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{\lambda }{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right)+\mathrm{T}\mathrm{r}\left({\boldsymbol{\varPhi }}_{n}{\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)+\sum \limits_{l=1}^{N}\mathrm{T}\mathrm{r}\left({\boldsymbol{\varPsi }}_{l}{{\boldsymbol{A}}^{\left(l\right)}}^{{}^{\mathrm{T}}}\right)\end{array} $ (6)
2.3 迭代规则更新

为了求解式(6),采用交替非负最小二乘迭代方法,即每次只更新核心张量或一个因子矩阵,同时固定其他因子矩阵。

2.3.1 因子矩阵$ {\boldsymbol{A}}^{\left(n\right)}(n\ne N) $更新

$ {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}} $对于$ {\boldsymbol{A}}^{\left(n\right)}, n\ne N $的偏导数为:

$ \frac{\partial {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}}{\partial {\boldsymbol{A}}^{\left(n\right)}}={\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}- \\ \;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{Y}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}+{\boldsymbol{\varPsi }}_{n}$

由KKT条件,令$ \frac{\partial {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}}{\partial {\boldsymbol{A}}^{\left(n\right)}}=0 $$ {\left({\boldsymbol{\varPsi }}_{n}\right)}_{ij}{\boldsymbol{A}}_{ij}^{\left(n\right)}=0 $,其中$ {\left({\boldsymbol{\varPsi }}_{n}\right)}_{ij} $$ {\boldsymbol{A}}_{ij}^{\left(n\right)} $表示$ {\boldsymbol{\varPsi }}_{n} $$ {\boldsymbol{A}}^{\left(n\right)} $中的$ \left(i, j\right) $元素,从而得到等式:

$ {\left({\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}-{\boldsymbol{Y}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}\cdot {\boldsymbol{A}}_{ij}^{\left(n\right)}=0 $

因此,得到以下更新规则:

$ {\boldsymbol{A}}_{ij}^{\left(n\right)}\leftarrow {\boldsymbol{A}}_{ij}^{\left(n\right)}\frac{{\boldsymbol{P}}_{+}\left({\left({\boldsymbol{Y}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}\right)}{{\left({\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}} $ (7)

其中:$ {\boldsymbol{P}}_{+}\left(\xi \right)=\mathrm{m}\mathrm{a}\mathrm{x}(0, \xi ) $

2.3.2 因子矩阵$ {\boldsymbol{A}}^{\left(N\right)} $更新

$ {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}} $关于$ {\boldsymbol{A}}^{\left(N\right)} $的偏导数为:

$ \frac{\partial {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}}{\partial {\boldsymbol{A}}^{\left(N\right)}}={\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{S}}_{\left(N\right)}\left(\underset{p\ne N}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(N\right)}^{\mathrm{T}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{Y}}_{\left(N\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(N\right)}^{\mathrm{T}}+\lambda {\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}+{\boldsymbol{\varPsi }}_{N}$

同理,由KKT条件可以得到以下更新规则:

$ {\boldsymbol{A}}_{ij}^{{}^{\left(N\right)}}\leftarrow {\boldsymbol{A}}_{ij}^{{}^{\left(N\right)}}\frac{{\boldsymbol{P}}_{+}\left({\left({\boldsymbol{Y}}_{\left(N\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(N\right)}^{\mathrm{T}}+\lambda {\boldsymbol{A}}^{\left(N\right)}\bf{H}\bf{W}{\boldsymbol{D}}_{e}^{-1}{\boldsymbol{H}}^{\mathrm{T}}\right)}_{ij}\right)}{{\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{S}}_{\left(N\right)}\left(\underset{p\ne N}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(N\right)}^{\mathrm{T}}+\lambda {\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{D}}_{v}\right)}_{ij}} $ (8)
2.3.3 核张量S更新

基于式(5)采用拉格朗日乘子法,并考虑其向量化形式:

$ {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}=\frac{1}{2}{‖\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{Y}\right)-\boldsymbol{Q}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)‖}_{F}^{2}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; \frac{\lambda }{2}\mathrm{T}\mathrm{r}\left({\boldsymbol{A}}^{\left(N\right)}{\boldsymbol{L}}^{\mathrm{H}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}{\boldsymbol{A}}^{{\left(N\right)}^{\mathrm{T}}}\right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathrm{v}\mathrm{e}\mathrm{c}{\left(\boldsymbol{S}\right)}^{\mathrm{T}}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{\varPhi }\right)+\sum \limits_{l=1}^{N}{\rm{Tr}}\left({\boldsymbol{\varPsi }}_{l}{{\boldsymbol{A}}^{\left(l\right)}}^{{}^{T}}\right) $

其中:$ \boldsymbol{Q}={\boldsymbol{A}}^{\left(N\right)}\otimes \cdots \otimes {\boldsymbol{A}}^{\left(2\right)}\otimes {\boldsymbol{A}}^{\left(1\right)}\in {\mathbb{R}}_{+}^{{I}_{1}{I}_{2}\cdots {I}_{N}\times {J}_{1}{J}_{2}\cdots {J}_{N}} $$ \mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{\varPhi }\right) $$ \mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right) $的拉格朗日乘子。

$ {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}} $对于$ \mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right) $的偏导数为:

$ \frac{\partial {\mathcal{L}}_{\mathrm{H}\mathrm{G}\mathrm{N}\mathrm{T}\mathrm{D}}}{\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)}={\boldsymbol{Q}}^{\mathrm{T}}\boldsymbol{Q}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)-{\boldsymbol{Q}}^{\mathrm{T}}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{Y}\right)+\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{\varPhi }\right) $

同上,由KKT条件从而得到以下更新规则:

$ {\left(\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)\right)}_{i}\leftarrow {\left(\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)\right)}_{i}\frac{{P}_{+}\left({\left({\boldsymbol{Q}}^{\mathrm{T}}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{Y}\right)\right)}_{i}\right)}{{\left({\boldsymbol{Q}}^{\mathrm{T}}\boldsymbol{Q}\mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)\right)}_{i}} $ (9)

根据以上分析,给出本文算法。

算法1   HGNTD算法

输入  $ \boldsymbol{y}\in {\mathbb{R}}_{+}^{{I}_{1}\times {I}_{2}\times \cdots \times {I}_{N}} $,超图正则化参数$ \lambda > 0 $,邻域大小$ k > 0 $

输出  因子矩阵$ {\boldsymbol{A}}^{\left(n\right)}\left(n\in N\right) $与核张量

1.随机初始化因子矩阵$ {\mathrm{A}}^{\left(\mathrm{n}\right)}\left(\mathrm{n}\in \mathrm{N}\right) $

2.随机初始化核张量$ \mathrm{S} $

3.根据式(7)更新因子矩阵$ {\mathrm{A}}^{\left(\mathrm{n}\right)}\left(\mathrm{n}\in \mathrm{N}, \mathrm{n}\ne {\mathrm{N}}_{\mathrm{n}}\right) $

4.根据式(8)更新因子矩阵$ {\mathrm{A}}^{\left(\mathrm{N}\right)} $

5.根据式(9)核张量$ \mathrm{S} $

6.直到满足收敛条件时停止。

2.4 收敛性分析

定理1  若$ {\boldsymbol{A}}^{\left(n\right)}\ge 0\left(n\in N\right), \mathrm{v}\mathrm{e}\mathrm{c}\left(\boldsymbol{S}\right)\ge 0 $,则目标函数式(5)在更新迭代规则式(7)、式(8)和式(9)下是非增的。

为证明该定理,首先引入辅助函数的定义:

定义1  若$ G\left(x, {x}^{\text{'}}\right) $满足$ G\left(x, x\right)=F\left(x\right) $$ G\left(x, {x}^{\text{'}}\right)\ge $ $ F\left(x\right) $,则$ G\left(x, {x}^{\text{'}}\right) $$ F\left(x\right) $的辅助函数[29]

引理1  若$ G\left(x, {x}^{\text{'}}\right) $$ F\left(x\right) $的辅助函数,则$ F\left(x\right) $在更新迭代规则[29]

$ {x}^{t+1}=\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\;G\left(x, {x}^{t}\right) $ (10)

下是非增的。

证明  根据定义1和更新迭代规则式(10),有:

$ F\left({x}^{t+1}\right)\le G\left({x}^{t+1}, {x}^{t}\right)\le G\left({x}^{t}, {x}^{t}\right)=F\left({x}^{t}\right) $

为证明目标函数式(5)在更新规则式(7)下是非增的,先构造关于矩阵$ {\boldsymbol{A}}^{\left(n\right)} $$ (i, j) $元素$ {\boldsymbol{A}}_{ij}^{\left(n\right)} $的辅助函数。对固定的行列指标$ (i, j) $,令$ {\boldsymbol{A}}^{\left(n\right)}\left(x\right) $是通过将矩阵$ {\boldsymbol{A}}^{\left(n\right)} $$ (i, j) $元素替换为变量$ x $,其余元素不变得到的矩阵函数。

引理2  函数:

$ G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)={F}_{ij}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)+{F}_{ij}^{\mathrm{\text{'}}}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\left(x-{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)+\\ \frac{1}{2}\cdot \frac{{\left(\left({\boldsymbol{A}}^{\left(n\right)}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\right){\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}}{{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}}{\left(x-{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)}^{2} $ (11)

$ {F}_{ij}\left(x\right) $的辅助函数,其中:

$ {F}_{ij}\left(x\right)=\frac{1}{2}{‖{\left({\boldsymbol{Y}}_{\left(n\right)}\right)}_{i:}-{\left(\left({\boldsymbol{A}}^{\left(n\right)}\left(x\right)\right){\boldsymbol{S}}_{\left(n\right)}(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}{)}^{\mathrm{T}}\right)}_{i:}‖}_{F}^{2} $ (12)

证明  显然,$ G\left(x, x\right)={F}_{ij}\left(x\right) $。要证明$ G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right) $$ {F}_{ij}\left(x\right) $的辅助函数,只需证明$ G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\ge {F}_{ij}\left(x\right) $即可。由式(12)得:

$ {F}_{ij}^{\mathrm{\text{'}}}\left(x\right)={\left(\left({\boldsymbol{A}}^{\left(n\right)}\left(x\right)\right){\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\boldsymbol{Y}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)_{ij}} $
$ {F}_{ij}^{\mathrm{\text{'}}\mathrm{\text{'}}}\left(x\right)={\left({\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{jj} $

将函数$ {F}_{ij}\left(x\right) $$ {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}} $处泰勒展开得:

$ {F}_{ij}\left(x\right)={F}_{ij}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)+{F}_{ij}^{\mathrm{\text{'}}}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\left(x-{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)+ \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{F}_{ij}^{\mathrm{\text{'}}\mathrm{\text{'}}}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right){\left(x-{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)}^{2} $ (13)

因为:

$ {\left(\left({\boldsymbol{A}}^{\left(n\right)}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\right){\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}=\\ \sum \limits_{l=1}^{{J}_{n}}{\boldsymbol{A}}_{il}^{\left(n\right)}{\left({\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{lj}\ge \\ {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}{\left({\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{jj} $

所以:

$ \frac{{\left(\left({\boldsymbol{A}}^{\left(n\right)}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right)\right){\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{{}^{\mathrm{T}}}\right)}_{ij}}{{\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}}\ge \\{\left({\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{{}^{\mathrm{T}}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{{}^{\mathrm{T}}}\right)}_{jj}={F}_{ij}^{\mathrm{\text{'}}\mathrm{\text{'}}}\left({\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right) $

比较式(11)、式(13),得$ G\left(x, {A}_{ij}^{{\left(n\right)}^{t}}\right)\ge {F}_{ij}\left(x\right) $,因此,$ G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right) $$ {F}_{ij}\left(x\right) $的辅助函数。

下面证明定理1。

证明  用式(11)中的$ G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right) $代替式(10)中的$ G\left(x, {x}^{t}\right) $,即取$ {\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t+1}}}=\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\;G\left(x, {\boldsymbol{A}}_{ij}^{{\left(n\right)}^{t}}\right) $

因此:

$ {\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t+1}}}={\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t}}}-{\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t}}}\frac{{F}^{\text{'}}\left({\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t}}}\right)}{{\left({\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}}= \\ \;\;\;\;\;\;\;\;\;{\boldsymbol{A}}_{ij}^{{}^{{\left(n\right)}^{t}}}\frac{{\left({\boldsymbol{Y}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }{\boldsymbol{A}}^{\left(p\right)}\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}}{{\left({\boldsymbol{A}}^{\left(n\right)}{\boldsymbol{S}}_{\left(n\right)}\left(\underset{p\ne n}{\otimes }\left({{\boldsymbol{A}}^{\left(p\right)}}^{{}^{\mathrm{T}}}{\boldsymbol{A}}^{\left(p\right)}\right)\right){\boldsymbol{S}}_{\left(n\right)}^{\mathrm{T}}\right)}_{ij}} $ (14)

由于式(11)是$ F\left({\boldsymbol{A}}_{ij}^{\left(n\right)}\right) $的辅助函数,因此在更新迭代规则式(14)(即式(7))下目标函数式(5)是非增的。同理可证,在更新迭代规则式(8)与式(9)下目标函数式(5)也是非增的。故而定理1成立。

由定理1知算法1是局部收敛的。

3 数值实验

现实中的高维图像数据依据人们能否获得先验信息可以分成两类:一类是可以通过人工标记或可以直接获得先验标签信息的数据;另一类是没有任何类别标签信息的数据。

图像数据分析与处理的方法分为有标签数据的图像分类和无标签数据的图像聚类。它们的区别是:分类是事先定义好类别,类别数不变。分类器需要由人工标注的分类训练语料训练得到,属于有监督学习。而聚类是没有事先规定类别,类别数不确定。聚类不需要人工标注和预先训练分类器,类别在聚类过程中自动生成,属于无监督学习。

本文进行的是无监督的聚类实验,因此无需训练集和测试集,数据集的标签只用于和聚类实验结果进行对比。

为验证本文算法的有效性,选择如下代表性算法进行比较:NMF[5],GNMF[14],NTD[18],GNTD[21]和K-means[30]。本文分别在两个常用公开的数据集Yale和COIL-100上进行聚类实验,实验环境为:Matlab 2019a,Intel® Core(TM) i5-7500 CPU@ 3.40 GHz,3.41 GHz处理器,8.00 GB内存,64位的Win10操作系统。

3.1 数据描述

实验中所使用数据集如下:

1)Yale。该数据集是由耶鲁大学计算视觉与控制中心创建,包含15位志愿者,每个采集志愿者包含光照、表情和姿态的变化以及遮挡脸部(不戴眼镜、戴普通眼睛和戴墨镜)状态下的11张样本,总共有165幅分辨率为100×100像素的人脸图像。在实验中,Yale数据集的所有图像都被调整为32×32像素的大小,Yale数据集的部分图片如图 2所示。

Download:
图 2 Yale数据集的部分图片 Fig. 2 Some images of Yale dataset

2)COIL-100。该数据集是哥伦比亚大学采集的(COIL-100)图像数据集。COIL-100由7 200幅尺寸为128×128像素的彩色图像组成,包含100个目标,每个目标有72幅不同角度的图像。在实验中,每个彩色图像被调整为大小32×32像素的灰度图像。COIL-100数据集的部分图片如图 3所示。

Download:
图 3 COIL-100数据集的部分图片 Fig. 3 Some images of COIL-100 dataset

传统的非负矩阵分解方法将每幅图像向量化后构成相应的1 024×165和1 024×7 200的矩阵再做分解,破坏了图像像素几何结构关系。非负张量分解方法直接对32×32×11K和32×32×72K(其中K是簇数)两个张量进行分解,不会破坏图像原有的结构。

3.2 评估指标

为评估聚类效果,使用两个流行的评估指标来评估聚类效果:一种度量标准是聚类准确度(Accuracy);另一种度量标准是归一化互信息(Normalized Mutual Information,NMI)。

聚类准确度(Accuracy,ACC)定义如下[31]

$ {A}_{\mathrm{A}\mathrm{C}\mathrm{C}}=\frac{\sum \limits_{i=1}^{n}\delta \left({C}_{i}, \mathrm{m}\mathrm{a}\mathrm{p}\left({C}_{i}^{\text{'}}\right)\right)}{n} $

其中:$ {C}_{i} $是真实类标签集合;$ {C}_{i}^{\text{'}} $是算法得到的聚类标签集合;$ n $是数据样本总数;$ \delta \left(x, y\right) $是一个kronecker函数,如果$ x=y $,则$ \delta \left(x, y\right)=1 $,否则$ \delta \left(x, y\right)=0 $,并且$ \mathrm{m}\mathrm{a}\mathrm{p}\left({C}_{i}^{\text{'}}\right) $是将每个聚类标签$ {C}_{i}^{\text{'}} $映射到真实类标签集$ {C}_{i} $中等效标签的映射函数。为得到最佳映射,通常使用Kuhn-Munkres算法。

给定两个图片聚类集$ {C}_{i} $$ {C}_{i}^{\text{'}} $,它们之间的互信息度量M$ {}_{\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right) $如下:

$ {M}_{\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right)=\sum \limits_{{c}_{i}\in C, {c}_{j}^{\text{'}}\in {C}^{\text{'}}}p\left({c}_{i}, {c}_{j}^{\text{'}}\right)\cdot \mathrm{l}\mathrm{b}\frac{p\left({c}_{i}, {c}_{j}^{\text{'}}\right)}{p\left({c}_{i}\right)p\left({c}_{j}^{\text{'}}\right)} $

其中:$ p\left({c}_{i}\right) $$ p\left({c}_{j}^{\text{'}}\right) $分别表示随机选取一个样本分别属于$ {c}_{i} $$ {c}_{j}^{\text{'}} $的概率;$ p\left({c}_{i}, {c}_{j}^{\text{'}}\right) $是随机选取的样本同时属于$ {c}_{i} $$ {c}_{j}^{\text{'}} $的联合概率;$ {M}_{\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right) $的值在0和$ \mathrm{m}\mathrm{a}\mathrm{x}\left(H\left(C\right), H\left({C}^{\text{'}}\right)\right) $之间,$ H\left(C\right) $$ H\left({C}^{\text{'}}\right) $分别是$ C $$ {C}^{\text{'}} $的熵。为了简化对不同集合之间的比较,使用以下归一化度量标准来代替$ {M}_{\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right) $。归一化互信息$ {N}_{\mathrm{N}\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right) $取0到1之间的值,当两个聚类集合完全独立时是0,当两个聚类集合完美匹配时是1,即两个样本点相似度越高,NMI的值越接近于1。

$ {N}_{\mathrm{N}\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right)=\frac{{M}_{\mathrm{M}\mathrm{I}}\left({C}_{i}, {C}_{i}^{\text{'}}\right)}{\mathrm{m}\mathrm{a}\mathrm{x}\left(H\left(C\right), H\left({C}^{\text{'}}\right)\right)} $
3.3 对比算法与参数设置

为验证本文所提算法的有效性,将HGNTD分别与NMF[5]、GNMF[14]、NTD[18]、GNTD[21]和k-means[30]算法进行比较,具体参数如下:

1)基于NMF算法的聚类。目标函数中距离使用Frobenius范数度量,基矩阵的行数$ r=10 $,即低维子空间的维数为$ r $

2)基于GNMF算法的聚类。目标函数中距离使用Frobenius范数度量,并采用热核(Heat Kernel)加权方案,使用k-最邻近方法构造权重矩阵$ \boldsymbol{W} $,其中$ k=3 $,图的正则化参数$ \lambda $取0.5。

3)基于NTD算法的聚类。目标函数中距离使用Frobenius范数度量,核张量的规模为$ \frac{1}{2}{I}_{1}\times \frac{1}{2}{I}_{2}\times \cdots \times r $

4)基于GNTD算法的聚类。目标函数中距离使用Frobenius范数度量,核张量的规模同NTD算法。采用热核加权方案,GNTD的其他参数设置与GNMF的设置相同。

5)基于HGNTD算法的聚类。目标函数中使用Frobenius范数度量距离。核张量的规模同NTD算法。采用热核加权方案,使用$ k $-最邻近方法构造权重矩阵W,其中$ k=3 $,超图正则化参数$ \lambda =0.5 $

3.4 实验与结果分析

具体实验步骤如下:

1)为使实验更具有说服力,对聚类实验配置进行随机化处理,即在不同的类别个数进行聚类实验。每次随机选择数据库中K个类别(COIL-100:K=2,4,…,20;Yale:K=3,5,…,15)。

2)对于k-means算法,作为一种基线(Baseline)算法,只在原始数据空间中进行聚类。对于NMF、GNMF、NTD、GNTD和HGNTD算法,在进行分解后,可以得到各个对比算法对应的低维子空间表示,那么聚类则在这些低维子空间中进行。然后,再将k-means应用于这些重构表征上进行聚类。

3)计算两个聚类指标精度(ACC)和归一化互信息(NMI)。

4)重复上述步骤20次,并记录各个聚类算法的平均性能以及上下偏差情况。

Yale数据集的聚类准确度和归一化互信息的结果如表 2表 3所示,COIL-100数据集的聚类准确度和归一化互信息的结果如表 4表 5所示,其中粗体为最优值。

下载CSV 表 2 Yale数据集聚类准确度 Table 2 Clustering accuracy of Yale dataset  
下载CSV 表 3 Yale数据集聚类归一化互信息 Table 3 Clustering normalized mutual information of Yale dataset  
下载CSV 表 4 COIL-100数据集聚类准确度 Table 4 Clustering accuracy of COIL-100 dataset  
下载CSV 表 5 COIL-100数据集聚类归一化互信息 Table 5 Clustering normalized mutual information of COIL-100 dataset  

表 2~表 5可以看出:在保留几何信息和内部结构的情况下,HGNTD与NMF、GNMF、NTD、GNTD和k-means算法相比,聚类性能更好。本文方法比其他方法的聚类准确度提升了8.6~11.4%,归一化互信息提升了2.0~7.5%,这也表明本文提出算法对图像聚类是有效的。

表 6表 7给出Yale数据集上NMF、GNMF、NTD、GNTD和HGNTD算法分别在维数$ r=4 $$ r=5 $下的对比结果,其中粗体为最优值。

下载CSV 表 6 Yale数据集上不同算法的聚类准确度 Table 6 Clustering accuracy of different algorithms on Yale dataset  
下载CSV 表 7 Yale数据集上不同算法聚类归一化互信息 Table 7 Clustering normalized mutual information of different algorithms on Yale dataset  

表 6表 7可以看出:在不同维数下,本文算法的聚类结果要比其他算法的聚类结果更好。

图 4所示为NMF、GNMF、NTD、GNTD和HGNTD算法在Yale数据集上的人脸重构效果,分别给出在Yale数据集上的部分图像,其中,第1行为原始图像,第2行~第6行分别为HGNTD算法、GNTD算法、GNTD算法、GNMF算法和NMF算法的样本重构图像。其中,HGNTD算法的重构图像比其他算法下的样本重构图像更加清晰,复原效果更好。

Download:
图 4 Yale数据集部分样本重构图像 Fig. 4 Reconstructed images from partial samples of Yale dataset
4 结束语

本文提出一种基于超图正则化非负Tucker分解算法(HGNTD)。利用$ k $-最邻近方法构造超图,将超图正则项添加到非负Tucker分解中,建立基于超图正则化非负Tucker分解模型,基于交替非负最小二乘法算法,给出超图正则化非负Tucker分解算法求解此模型,并证明算法是收敛的。在公开数据集上的实验结果表明,与k-means、NMF、GNMF、NTD、GNTD等算法相比,该算法具有更好的聚类效果。由于在本文实验中发现图像只需少数基底表示即可,因此下一步将继续对基于稀疏约束的超图正则化非负张量分解算法进行研究,以提高模型的适用性。

参考文献
[1]
VASWANI N, BOUWMANS T, JAVED S, et al. Robust subspace learning: robust PCA, robust subspace tracking, and robust subspace recovery[EB/OL]. [2020-12-10]. https://arxiv.org/pdf/1711.09492.pdf.
[2]
LIU J, OSHER S. Block matching local SVD operator based sparsity and TV regularization for image denoising[J]. Journal of Scientific Computing, 2019, 78(1): 607-624. DOI:10.1007/s10915-018-0785-8
[3]
SAFONT G, SALAZAR A, VERGARA L, et al. Probabilistic distance for mixtures of independent component analyzers[J]. IEEE Transactions on Neural Networks and Learning Systems, 2018, 29(4): 1161-1173. DOI:10.1109/TNNLS.2017.2663843
[4]
FU X, HUANG K, SIDIROPOULOS N D, et al. Nonnegative matrix factorization for signal and data analytics: identifiability, algorithms, and applications[EB/OL]. [2020-12-10]. https://arxiv.org/pdf/1803.01257.pdf.
[5]
LEE D D, SEUNG H S. Learning the parts of objects by non-negative matrix factorization[J]. Nature, 1999, 401(6755): 788-791. DOI:10.1038/44565
[6]
YAN S, XU D, YANG Q, et al. Multilinear discriminant analysis for face recognition[J]. IEEE Transactions on Image Processing, 2007, 16(1): 212-220. DOI:10.1109/TIP.2006.884929
[7]
CICHOCKI A, MANDIC D, DE LATHAUWER L, et al. Tensor decompositions for signal processing applications from two-way to multiway component analysis[J]. IEEE Signal Processing Magazine, 2015, 32(2): 145-163. DOI:10.1109/MSP.2013.2297439
[8]
HARSHMAN R A. Foundations of the parafac procedure: models and conditions for an "explanatory" multimodal factor analysis[EB/OL]. [2020-12-10]. https://psychology.uwo.ca/faculty/harshman/wpppfac0.pdf.
[9]
TUCKER L R. Some mathematical notes on three-mode factor analysis[J]. Psychometrika, 1966, 31(3): 279-311. DOI:10.1007/BF02289464
[10]
KOLDA T G, TAMARA G, BADER, et al. Tensor decompositions and applications[J]. SIAM Review, 2009, 51(3): 455-500. DOI:10.1137/07070111X
[11]
ZHOU G, CICHOCKI A, ZHAO Q, et al. Efficient nonnegative tucker decompositions: algorithms and uniqueness[J]. IEEE Transactions on Image Processing, 2015, 24(12): 4990-5003. DOI:10.1109/TIP.2015.2478396
[12]
LI S Z, HOU X W, ZHANG H J, et al. Learning spatially localized, parts-based representation[C]//Proceedings of 2001 IEEE Computer Society Conference on Computer Vision & Pattern Recognition. Kauai, USA: IEEE Press, 2001: 1582-1598.
[13]
CAI D, HE X, WU X, et al. Non-negative matrix factorization on manifold[C]//Proceedings of the 8th IEEE International Conference on Data Mining. Washington D. C., USA: IEEE Press, 2008: 63-72.
[14]
CAI D, HE X, HAN J, et al. Graph regularized nonnegative matrix factorization for data representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2011, 33(8): 1548-1560. DOI:10.1109/TPAMI.2010.231
[15]
LENG C, ZHANG H, CAI G, et al. Graph regularized LP smooth non-negative matrix factorization for data representation[J]. IEEE/CAA Journal of Automatica Sinica, 2019, 6(2): 584-595. DOI:10.1109/JAS.2019.1911417
[16]
WANG Y, GUI L, ZHANG Y. Neighborhood preserving non-negative tensor factorization for image representation[C]//Proceedings of IEEE International Conference on Acoustics. Washington D. C., USA: IEEE Press, 2012: 3889-3392.
[17]
WANG C, HE X, BU J, et al. Image representation using Laplacian regularized nonnegative tensor factorization[J]. Pattern Recognition, 2011, 44(10/11): 2516-2526.
[18]
KIM Y D, CHOI S. Nonnegative tucker decomposition[C]// Proceedings of IEEE Conference on Computer Vision & Pattern Recognition. Washington D. C., USA: IEEE Press, 2008: 1-8.
[19]
JIANG B, DING C, TANG J, et al. Image representation and learning with graph-Laplacian tucker tensor decomposition[J]. IEEE Transactions on Cybernetics, 2019, 49(4): 1417-1426. DOI:10.1109/TCYB.2018.2802934
[20]
LI X, NG M K, CONG G, et al. MR-NTD: manifold regularization nonnegative tucker decomposition for tensor data dimension reduction and representation[J]. IEEE Transactions on Neural Networks and Learning Systems, 2017, 28(8): 1787-1800. DOI:10.1109/TNNLS.2016.2545400
[21]
QIU Y, ZHOU G, ZHANG Y, et al. Graph regularized nonnegative tucker decomposition for tensor data representation[C]//Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing. Washington D. C., USA: IEEE Press, 2019: 8613-8617.
[22]
ZHOU D, HUANG J, SCHLKOPF B. Learning with hypergraphs: clustering, classification, and embedding[C]//Proceedings of the 20th Annual Conference on Neural Information Processing Systems. Vancouver, Canada: MIT Press, 2006: 134-146.
[23]
AGARWAL S, BRANSON K, BELONGIE S. Higher order learning with graphs[C]//Proceedings of the 23rd International Conference on Machine Learning. Pittsburgh, USA: [s. n.], 2006: 121-128.
[24]
LIU M, ZHANG J, GUO X, et al. Hypergraph regularized sparse feature learning[J]. Neurocomputing, 2017, 237(10): 185-192.
[25]
ZENG K, YU J, LI C, et al. Image clustering by hyper-graph regularized non-negative matrix factorization[J]. Neurocomputing, 2014, 138(22): 209-217.
[26]
WANG W, QIAN Y, TANG Y Y. Hypergraph-regularized sparse NMF for hyperspectral unmixing[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(2): 681-694. DOI:10.1109/JSTARS.2015.2508448
[27]
罗永恩, 胡继承, 徐茜. 基于超图的多模态关联特征处理方法[J]. 计算机工程, 2017, 43(1): 226-230.
LUO Y E, HU J C, XU Q. Multimodal correlation feature processing method based on hypergraph[J]. Computer Engineering, 2017, 43(1): 226-230. (in Chinese)
[28]
BELKIN M, NIYOGI P. Laplacian eigenmaps and spectral techniques for embedding and clustering[J]. Advances in Neural Information Processing Systems, 2001, 14(6): 585-591.
[29]
LEE D D, SEUNG H S. Algorithms for non-negative matrix factorization[C]//Proceedings of International Conference on Neural Information Processing Systems. [S. 1.]: MIT Press, 2000: 135-153.
[30]
MACQUEEN J B. Some methods for classification and analysis of multivariate observations[C]//Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability. Los Angeles, USA: University of California Press, 1967: 238-247.
[31]
XU W, LIU X, GONG Y. Document clustering based on non-negative matrix factorization[C]//Proceedings of the 26th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval. New York, USA: ACM Press, 2003: 267-273.