摘要部分
传统的基于图的聚类算法在大规模高光谱图像聚类问题中效果不好,主要是因为较高的计算复杂度。
本文中提出新的方法,称为基于锚图的快速谱聚类。fast spectral clustering with anchor graph (FSCAG)
具体来说,在构造锚图中考虑了高光谱的频谱和空间特性。
首先构造锚图,然后对锚图进行谱聚类。可以将复杂度降低到 O ( n d m ) O(ndm) O(ndm),传统的基于图的聚类方法至少需要 O ( n 2 d ) O(n^2d) O(n2d)
n , d , m n,d,m n,d,m分别是样本数,特征数,锚点数
介绍
高光谱图像HSI:
数据可看作三维立方体,有复杂空间结构,其中每一个像素点抽取出来都是一条连续的光谱曲线。
X,Y轴表示空间信息,Z轴表示光谱信息。
HSI聚类的目的:
将给定的图像分成组,使得同一组中的像素尽可能彼此相似,而分配给不同组的像素不相似。
HSI聚类的复杂性:
- 相似矩阵的构造, O ( n 2 d ) O(n^2d) O(n2d)时间复杂度
- 拉普拉斯矩阵的特征值分解
算法实现部分
首先结合HSI的空间信息构建锚图。然后,设计了一种对该图进行谱分析的有效方法。
锚图构造
设:
X
=
[
x
1
,
…
,
x
n
]
T
∈
R
n
×
d
X=[x_1,\ldots,x_n]^T\in\mathbb{R}^{n\times d}
X=[x1,…,xn]T∈Rn×d
表示数据矩阵
其中,n是数据点的数量,d是特征维度,每个数据点 x i ∈ R d x_i\in\mathbb{R}^d xi∈Rd 属于c个类之一。
给定数据矩阵 X X X,每个数据点 x i x_i xi表示为afinity graph(节点是数据样本,边代表数据之间的相似度)上的一个顶点,每个边表示一对顶点的相似关系。
x i x_i xi和 x j x_j xj之间边的权重定义为 a i j a_{ij} aij, A = { a i j } ∈ R n × n , ∀ i , j ∈ 1 , … n A=\{a_{ij}\}\in\mathbb{R}^{n\times n},\forall i,j\in1,\ldots n A={aij}∈Rn×n,∀i,j∈1,…n 表示afinity graph的相似矩阵。
最近的研究采用基于锚的策略来构建相似矩阵A。基于锚的战略首先从数据点生成m个锚,其中 m ≪ n m\ll n m≪n,然后设计矩阵 Z ∈ R n × m Z\in\mathbb{R}^{n\times m} Z∈Rn×m,测量数据点和锚之间的相似性。
基于锚的策略的第一步是锚的生成,可以通过随机选择或使用k-means方法来实现。
-
随机选择通过从计算复杂度 O(1) 的数据点中随机抽样来选择 m 个锚点
虽然随机选择不能保证选择的 m 个 anchor 总是好的,但是对于大规模 HSI 聚类来说是极快的。
-
k-means 方法使用 m 个聚类中心作为锚点。 聚类中心是比较有代表性的anchors,但是k-means的计算复杂度是 O ( n d m t ) O(ndmt) O(ndmt),其中t是迭代次数,大规模HSI聚类不可能。
设:
U
=
[
u
1
,
…
,
u
m
]
T
∈
R
m
×
d
U=[u_1,\ldots,u_m]^T \in \mathbb{R}^{m\times d}
U=[u1,…,um]T∈Rm×d
表示生成的锚。
Z中的第(i,j)个元素定义为:
z
i
j
=
K
(
x
i
,
u
j
)
∑
s
∈
Φ
i
K
(
x
i
,
u
s
)
,
∀
j
∈
Φ
i
(1)
z_{ij} = \frac{K(x_i,u_j)}{\sum_{s\in\Phi_i}^{} {K(x_i,u_s)}},\forall j\in \Phi_i \tag {1}
zij=∑s∈ΦiK(xi,us)K(xi,uj),∀j∈Φi(1)
其中 Φ i ⊂ 1 , … , m \Phi_i\subset{1,\ldots,m} Φi⊂1,…,m表示 U U U中 x i x_i xi的k近邻的索引, K ( ) K() K()是给定的核函数[14],
高斯核函数:
K
(
x
i
,
u
j
)
=
e
x
p
(
−
∣
∣
x
i
−
u
j
∣
∣
2
2
2
σ
2
)
K(x_i,u_j) = exp(\frac{-||x_i-u_j||_2^2}{2\sigma ^2})
K(xi,uj)=exp(2σ2−∣∣xi−uj∣∣22)
基于内核的方法总是会带来额外的参数,例如
σ
\sigma
σ,通过求解以下问题获得
Z
Z
Z的第i行[18]:
min
z
i
T
1
=
1
,
z
i
j
≥
0
∑
j
=
1
m
∣
∣
x
i
−
u
j
∣
∣
2
2
z
i
j
+
γ
z
i
j
2
(2)
\mathop{\min}_{z_i^T\textbf{1}=1,z_{ij}\geq0} \sum_{j=1}^{m} {||x_i-u_j||_2^2z_{ij}}+\gamma z_{ij}^2 \tag{2}
minziT1=1,zij≥0j=1∑m∣∣xi−uj∣∣22zij+γzij2(2)
z
i
T
表
示
的
是
Z
的
第
i
行
,
γ
是
正
则
化
参
数
z_i^T表示的是Z的第i行,\gamma是正则化参数
ziT表示的是Z的第i行,γ是正则化参数,等式(2)是无参数的
但是它没有考虑HSI的空间信息,这可能导致由于噪声、离群值或混合像素的存在,聚类HSI图像中出现一些孤立像素。
受[6]和[20]的启发,建议通过直接修改(2)中的目标函数来合并空间信息:
min
z
i
T
1
=
1
,
z
i
j
≥
0
∑
j
=
1
m
∣
∣
x
i
−
u
j
∣
∣
2
2
z
i
j
+
α
∣
∣
x
i
‾
−
u
j
∣
∣
2
2
z
i
j
+
γ
z
i
j
2
(3)
\mathop{\min}_{z_i^T\textbf{1}=1,z_{ij}\geq0} \sum_{j=1}^{m} {||x_i-u_j||_2^2z_{ij}}+\alpha ||\overline{x_i}-u_j||_2^2z_{ij} +\gamma z_{ij}^2 \tag{3}
minziT1=1,zij≥0j=1∑m∣∣xi−uj∣∣22zij+α∣∣xi−uj∣∣22zij+γzij2(3)
x
i
‾
\overline{x_i}
xi是位于
x
i
x_i
xi周围窗口内相邻像素的平均值,可以提前计算。参数
α
\alpha
α控制原始HSI和相应平均滤波HSI之间的平衡(折衷)。
d
i
j
x
=
∣
∣
x
i
−
u
j
∣
∣
2
2
,
d
i
j
x
‾
=
∣
∣
x
i
‾
−
u
j
∣
∣
2
2
,
d_{ij}^x=||x_i-u_j||_2^2,\quad d_{ij}^{\overline{x}}=||\overline{x_i} -u_j||_2^2,\quad
dijx=∣∣xi−uj∣∣22,dijx=∣∣xi−uj∣∣22,
(3)写成向量的形式:
min
z
i
∣
∣
z
i
+
1
2
γ
d
i
∣
∣
2
2
s.t.
z
i
T
1
=
1
,
z
i
j
≥
0
(4)
\mathop{\min}_{z_i} {||z_i+\frac{1}{2\gamma }d_i||_2^2} \quad \textbf{s.t.} \quad z_i^T\textbf{1}=1 ,z_{ij}\geq0\tag{4}
minzi∣∣zi+2γ1di∣∣22s.t.ziT1=1,zij≥0(4)
(4)的推导:
min
z
i
T
1
=
1
,
z
i
j
≥
0
∑
j
=
1
m
∣
∣
x
i
−
u
j
∣
∣
2
2
z
i
j
+
α
∣
∣
x
i
‾
−
u
j
∣
∣
2
2
z
i
j
+
γ
z
i
j
2
=
min
∑
i
=
1
m
∑
j
=
1
m
(
d
i
j
x
z
i
j
+
α
d
i
j
x
‾
z
i
j
)
+
γ
∣
∣
z
∣
∣
F
2
=
min
∑
i
=
1
m
∑
j
=
1
m
z
i
j
(
d
i
j
x
+
α
d
i
j
x
‾
)
+
γ
∑
i
=
1
m
z
i
T
z
i
=
min
∑
i
=
1
m
∑
j
=
1
m
z
i
j
d
i
j
+
γ
∑
i
=
1
m
z
i
T
z
i
=
min
∑
i
=
1
m
z
i
T
d
i
+
γ
∑
i
=
1
m
z
i
T
z
i
=
γ
min
∑
i
=
1
m
(
1
γ
z
i
T
d
i
+
z
i
T
z
i
)
=
γ
min
(
z
i
T
z
i
+
2
1
2
γ
z
i
T
d
i
+
d
i
T
d
i
4
γ
2
−
d
i
T
d
i
4
γ
2
)
最
后
一
项
常
数
舍
去
=
γ
min
∣
∣
z
i
+
1
2
γ
d
i
∣
∣
2
2
\begin{aligned} & \mathop{\min}_{z_i^T\textbf{1}=1,z_{ij}\geq0} \sum_{j=1}^{m} {||x_i-u_j||_2^2z_{ij}}+\alpha || \overline{x_i}-u_j||_2^2z_{ij} +\gamma z_{ij}^2 \\ & = \mathop{\min} \sum_{i=1}^m \sum_{j=1}^m (d_{ij}^x z_{ij} + \alpha d_{ij}^{\overline{x}} z_{ij} )+\gamma||z||_F^2 \\ & =\mathop{\min} \sum_{i=1}^m \sum_{j=1}^m z_{ij}(d_{ij}^x + \alpha d_{ij}^{\overline{x}}) + \gamma \sum_{i=1}^m z_i^Tz_i \\ & =\mathop{\min} \sum_{i=1}^m \sum_{j=1}^m z_{ij} d_{ij}+\gamma \sum_{i=1}^m z_i^Tz_i \\ & = \mathop{\min} \sum_{i=1}^m z_i^Td_i + \gamma \sum_{i=1}^m z_i^Tz_i \\ & = \gamma \mathop{\min} \sum_{i=1}^m (\frac{1}{\gamma}z_i^T d_i+ z_i^Tz_i) \\ & = \gamma \mathop{\min}( z_i^Tz_i + 2 \frac{1}{2\gamma} z_i^Td_i + \frac{d_i^Td_i}{4\gamma^2}-\frac{d_i^Td_i}{4\gamma^2}) \qquad最后一项常数舍去 \\ & = \gamma \mathop{\min} ||z_i+ \frac{1}{2\gamma} d_i||_2^2 \end{aligned}
minziT1=1,zij≥0j=1∑m∣∣xi−uj∣∣22zij+α∣∣xi−uj∣∣22zij+γzij2=mini=1∑mj=1∑m(dijxzij+αdijxzij)+γ∣∣z∣∣F2=mini=1∑mj=1∑mzij(dijx+αdijx)+γi=1∑mziTzi=mini=1∑mj=1∑mzijdij+γi=1∑mziTzi=mini=1∑mziTdi+γi=1∑mziTzi=γmini=1∑m(γ1ziTdi+ziTzi)=γmin(ziTzi+22γ1ziTdi+4γ2diTdi−4γ2diTdi)最后一项常数舍去=γmin∣∣zi+2γ1di∣∣22
依据[19]中的工作,最好学习具有精确k个非零值的稀疏
z
i
z_i
zi。因此,
Z
Z
Z是稀疏的,并且可以大大减轻后续频谱分析的计算负担,参数
γ
\gamma
γ可以设置为
γ
=
(
k
/
2
)
d
i
,
k
+
1
−
(
1
/
2
)
∑
j
=
1
k
d
i
,
j
\gamma=(k/2)d_{i,k+1}-(1/2) \sum_{j=1}^{k}d_{i,j}
γ=(k/2)di,k+1−(1/2)∑j=1kdi,j,(4)的解是:
z
i
j
=
d
i
,
k
+
1
−
d
i
,
j
k
d
i
,
k
+
1
−
∑
j
=
1
k
d
i
j
(5)
z_{ij}=\frac{d_{i,k+1}-d_{i,j}}{kd_{i,k+1}-\sum_{j=1}^kd_{ij}} \tag{5}
zij=kdi,k+1−∑j=1kdijdi,k+1−di,j(5)
详细计算过程参照[19],因此只需要
O
(
n
d
m
)
O(ndm)
O(ndm)来计算矩阵
Z
Z
Z,相似矩阵
A
A
A可通过[14]得到:
A
=
Z
Λ
−
1
Z
T
(6)
A=Z \Lambda ^{-1} Z^T \tag{6}
A=ZΛ−1ZT(6)
对角矩阵
Λ
∈
R
m
×
m
\Lambda\in \mathbb{R}^{m\times m}
Λ∈Rm×m 定义为
Λ
j
j
=
∑
i
=
1
n
z
j
j
\Lambda_{jj} =\sum_{i=1}^n z_{jj}
Λjj=∑i=1nzjj
锚图的谱分析
SC的目标函数:
min
F
T
F
=
I
t
r
(
F
T
L
F
)
(7)
\mathop{\min}_{F^TF=I}tr(F^TLF) \tag{7}
minFTF=Itr(FTLF)(7)
F
∈
B
n
×
c
F\in\mathbb{B}^{n\times c}
F∈Bn×c 是聚类指标矩阵,c是类数
L = D − A L=D-A L=D−A在图论中是拉普拉斯矩阵,度矩阵 D ∈ R n × n D\in \mathbb{R}^{n\times n} D∈Rn×n 定义为对角矩阵,第i个对角元素为 ∑ j = 1 n a i j \sum_{j=1}^naij ∑j=1naij
注意,F的元素被约束为离散值,这使得(7)很难求解。这个问题众所周知的解决方案是将矩阵F从离散值松弛为连续值。通过对L进行特征值分解,我们得到了由对应于最小c特征值的特征向量组成的松弛连续解。然后,可以使用k均值聚类来计算离散解。
拉普拉斯矩阵的特征值分解的时间复杂度是
O
(
n
2
c
)
O(n^2c)
O(n2c),并不适合大型HSI聚类。从(6)中得知
A
A
A的第
(
i
,
j
)
(i,j)
(i,j)个元素是
a
i
j
=
z
i
T
Λ
−
1
z
j
a_{ij}=z_i^T \Lambda^{-1} z_j
aij=ziTΛ−1zj,满足
a
i
j
T
=
a
j
i
a_{ij}^T=a_{ji}
aijT=aji,
A
A
A是对称阵。注意:
z
i
j
≥
0
z_{ij} \geq0
zij≥0并且
z
i
T
1
=
1
z_i^T \textbf{1}=1
ziT1=1,可以验证
A
A
A是满足
a
i
j
a_{ij}
aij的双随机矩阵。并且:
∑
j
=
1
n
a
i
j
=
∑
j
=
1
n
z
i
T
Λ
−
1
z
j
=
z
i
T
∑
j
=
1
n
Λ
−
1
z
j
=
z
i
T
1
=
1
(8)
\sum_{j=1}^{n}a_{ij}=\sum_{j=1}^{n}z_i^T \Lambda^{-1} z_j=z_i^T\sum_{j=1}^{n}\Lambda^{-1} z_j=z_i^T \textbf{1}=1 \tag{8}
j=1∑naij=j=1∑nziTΛ−1zj=ziTj=1∑nΛ−1zj=ziT1=1(8)
因此,相似矩阵
A
A
A被自动归一化,度矩阵
D
=
I
D=I
D=I,
I
I
I表示单位矩阵。
L
=
I
−
A
L=I-A
L=I−A,(7)等价于以下问题:
max
F
T
F
=
I
t
r
(
F
T
A
F
)
(9)
\mathop{\max}_{F^TF=I}tr(F^TAF) \tag{9}
maxFTF=Itr(FTAF)(9)
根据(6)
A
A
A可以写为
A
=
B
B
T
A=BB^T
A=BBT,其中
B
=
Z
Λ
−
(
1
/
2
)
B=Z\Lambda^{-(1/2)}
B=ZΛ−(1/2),
B
B
B的奇异值分解可以写为:
B
=
U
∑
V
T
B=U\sum V^T
B=U∑VT
其中 $U\in \mathbb{R}^{n\times n},\sum \in \mathbb{R}^{n \times m},V \in \mathbb{R}^{m\times m} $ 分别是左奇异向量矩阵、奇异值矩阵,右奇异向量矩阵。
很容易验证 U U U的列向量是相似矩阵 A = B B T A=BB^T A=BBT的特征向量。我们不直接对A进行特征值分解,而是对 B B B进行奇异值分解,以获得 F F F的松弛连续解。计算复杂度可以减少到 O ( n m c + m 2 c ) O(nmc+m^2c) O(nmc+m2c)
然后,使用k均值聚类来计算离散解。
计算复杂性分析
假设我们有一个数据矩阵 X ∈ R n × d X∈ \mathbb{R}^{n×d} X∈Rn×d,并从 X X X生成 m m m个锚。FSCAG的步骤和相应的计算复杂度总结如下:
- 需要 O ( 1 ) O(1) O(1)通过随机选择获得 m m m个锚。
- 需要 O ( n d m ) O(ndm) O(ndm)来获得矩阵 Z Z Z来构造锚图。
- 通过对矩阵B进行奇异值分解,需要 O ( n m c + m 2 c ) O(nmc+m^2c) O(nmc+m2c)来获得 F F F的松弛连续解。
- 需要 O ( n m c t ) O(nmct) O(nmct)对最终聚类结果的松弛连续解执行 k m e a n s kmeans kmeans,其中 t t t是迭代数。
考虑到 m ≪ n , c ≪ d m\ll n,c \ll d m≪n,c≪d,通常 t t t很小,FSCAG的总体计算复杂度为 O ( n d m ) O(ndm) O(ndm)。
在算法1中总结了我们算法的细节.
算法1 FSCAG算法
输入:HSI数据矩阵 X ∈ R n × d X∈ R^{n×d} X∈Rn×d,类的数量c,锚的数量m,折衷参数α,邻居数量k
- 通过随机选择生成m个锚。
- 通过(5)获得矩阵 Z Z Z。
- 通过对矩阵 B B B进行奇异值分解,得到 F F F的松弛连续解,其中 B = Z Λ − 1 2 B=Z\Lambda^{-\frac{1}{2}} B=ZΛ−21。
- 对F的松弛连续解执行 k − m e a n s k-means k−means,以获得最终聚类结果。
**输出:**c类
实验结果
对三个广泛使用的HSI数据集执行FSCAG。选择k均值[4]、FCM[5]、FCM_ S1[6]和SC[11]作为基准。k-均值、FCM和FCM_ S1的计算复杂度为 O ( n d c t ) O(ndct) O(ndct),而SC的计算复杂程度为 O ( n 2 d ) O(n^2d) O(n2d)
三个HSI数据集的定量评估 (OM ERROR):
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GMeNRATO-1665046344250)(三个HSI数据集的定量评估.png)]
三个HSI数据集上的运行时间:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-yX7v0hd1-1665046344252)(三个HSI数据集上的运行时间.png)]
聚类结果及分析
- 在印度松数据集上进行聚类实验。FSCAG中使用的参数分别设置为 m = 500 、 k = 5 和 α = 0.6 m=500、k=5和\alpha=0.6 m=500、k=5和α=0.6。k-means、FCM和SC获得较差的聚类性能,但相比之下,FCM_S1和FSCAG提高了聚类性能,以获得更平滑的聚类图,这清楚地反映了合并空间信息的重要性。
- 在Salinas数据集上进行聚类实验。参数分别设置为 m = 1000 、 k = 5 和 α = 0.8 m=1000、k=5和α=0.8 m=1000、k=5和α=0.8。FSCAG获得了最高的精度,最佳OA为74.63%,kappa系数为0.7173。与其他算法相比,FSCAG产生了更多的同质区域和更好的聚类图。k均值、FCM、FCM_S1和FSCAG的运行时间具有相同的数量级。Salinas数据集属于大规模HSI数据集,样本总数为111 104。专注于两种基于图的聚类方法,FSCAG只需要11.9秒,但由于“内存不足(OM)错误”,SC无法处理Salinas的数据集。
- 在Pavia中心数据集上进行聚类实验。参数分别设置为m=1000、k=5和α=0.3。FSCAG获得最高精度,最佳OA为81.55%,kappa系数为0.7441。我们看到FSCAG比其他算法产生更多的同质区域和更好的聚类图。此外,FCM_ S1和FSCAG通过考虑空间信息获得更好的聚类结果和更平滑的聚类图。
参数敏感度
有三个参数需要调整: m , k , α m,k,α m,k,α。
从第II-C节中,可以看到计算复杂度主要由参数m决定。参数k与矩阵Z的稀疏性相关,对计算复杂度几乎没有影响。此外,参数α与计算复杂度无关。
在Salinas数据集上进行参数敏感性实验,Salinas数据集选择较小的m=1000和k=5是合理的。
……
结论
摘要部分已给出……
_ S1和FSCAG通过考虑空间信息获得更好的聚类结果和更平滑的聚类图。
参数敏感度
有三个参数需要调整: m , k , α m,k,α m,k,α。
从第II-C节中,可以看到计算复杂度主要由参数m决定。参数k与矩阵Z的稀疏性相关,对计算复杂度几乎没有影响。此外,参数α与计算复杂度无关。
在Salinas数据集上进行参数敏感性实验,Salinas数据集选择较小的m=1000和k=5是合理的。
……
结论
摘要部分已给出……