随机化Nyström方法在谱聚类（Spectral Clustering）近似中的应用
题目描述：
在谱聚类（Spectral Clustering）中，核心步骤是计算一个大型相似矩阵（或亲和矩阵）\(W \in \mathbb{R}^{n \times n}\) 的 \(k\) 个最小特征值对应的特征向量（或拉普拉斯矩阵 \(L\) 的前 \(k\) 个特征向量）。由于 \(n\) 很大（例如数十万），对 \(W\) 或 \(L\) 进行精确特征分解（例如通过QR算法）代价过高（\(O(n^3)\)）。随机化Nyström方法（Randomized Nyström Method）通过对矩阵进行列采样，构造一个低秩近似，从而高效地近似这些特征向量。本题将详细介绍随机化Nyström方法在谱聚类近似中的算法步骤、数学原理、误差来源，并给出一个具体的计算流程示例。
解题过程循序渐进讲解：
第一步：问题形式化与背景知识
- 谱聚类的核心计算问题：
给定 \(n\) 个数据点，定义相似矩阵 \(W \in \mathbb{R}^{n \times n}\)，其中 \(W_{ij}\) 表示点 \(i\) 与点 \(j\) 的相似度（例如高斯核函数）。归一化拉普拉斯矩阵定义为：
\[ L = I - D^{-1/2} W D^{-1/2} \]
其中 \(D\) 是对角矩阵，\(D_{ii} = \sum_{j=1}^n W_{ij}\)。谱聚类的关键步骤是计算 \(L\) 的最小的 \(k\) 个特征值对应的特征向量（假设我们需要将数据分成 \(k\) 个簇）。精确计算 \(L\) 的特征分解需要 \(O(n^3)\) 的计算量和 \(O(n^2)\) 的存储空间，当 \(n\) 很大时不可行。
- 随机化Nyström方法的核心思想：
Nyström方法最初用于积分方程的数值解。在矩阵近似中，它通过对矩阵的列进行采样，利用采样到的列和相应的行来重建整个矩阵的低秩近似。对于对称半正定矩阵（如 \(W\)），随机化Nyström方法可以高效地近似其特征值和特征向量，而无需处理整个矩阵。
第二步：随机化Nyström方法的基本步骤
- 列采样：
从 \(n\) 个列中随机选择 \(l\) 个列（\(l \ll n\)），记下标集合为 \(I \subset \{1,2,\dots,n\}\)，大小为 \(l\)。将矩阵 \(W\) 分块为：
\[ W = \begin{bmatrix} W_{II} & W_{I\bar{I}} \\ W_{\bar{I}I} & W_{\bar{I}\bar{I}} \end{bmatrix} \]
其中 \(W_{II} \in \mathbb{R}^{l \times l}\) 是采样列之间的子矩阵，\(W_{I\bar{I}} \in \mathbb{R}^{l \times (n-l)}\) 是采样列与未采样列之间的子矩阵，且由于对称性 \(W_{\bar{I}I} = W_{I\bar{I}}^T\)。我们只计算（或已知）\(W_{II}\) 和 \(W_{I\bar{I}}\)，而 \(W_{\bar{I}\bar{I}}\) 未知。
- 构造近似：
定义矩阵 \(C = W_{:,I} \in \mathbb{R}^{n \times l}\)，即 \(W\) 中所有行、但仅包含被采样列 \(I\) 的子矩阵。定义 \(A = W_{II} \in \mathbb{R}^{l \times l}\)。经典的Nyström近似为：
\[ W \approx \tilde{W} = C A^{\dagger} C^T \]
其中 \(A^{\dagger}\) 是 \(A\) 的Moore-Penrose伪逆。如果 \(A\) 是满秩的，则 \(A^{\dagger} = A^{-1}\)。
-
对 \(A\) 进行特征分解：
计算 \(A\) 的特征分解：\(A = U_A \Lambda_A U_A^T\)，其中 \(U_A \in \mathbb{R}^{l \times l}\) 是正交特征向量矩阵，\(\Lambda_A = \text{diag}(\lambda_1, \dots, \lambda_l)\) 是特征值矩阵（按从大到小排列）。通常 \(A\) 是半正定的，特征值非负。 -
构造近似特征向量：
定义 \(V = C U_A \Lambda_A^{-1/2} \in \mathbb{R}^{n \times l}\)，其中 \(\Lambda_A^{-1/2}\) 是对角矩阵，元素为 \(1/\sqrt{\lambda_i}\)（对于非常小的特征值可以截断）。可以验证：
\[ \tilde{W} = V \Lambda_A V^T \]
这表明 \(V\) 的列近似是 \(W\) 的特征向量（尽管不一定正交），而 \(\Lambda_A\) 近似是特征值。但更常用的是正交化 \(V\) 以得到正交近似特征向量。
-
正交化：
对 \(V\) 进行QR分解：\(V = Q R\)，其中 \(Q \in \mathbb{R}^{n \times l}\) 的列标准正交，\(R \in \mathbb{R}^{l \times l}\) 是上三角矩阵。然后计算 \(B = Q^T W Q \approx Q^T (C A^{\dagger} C^T) Q\)。由于 \(C = W_{:,I}\) 和 \(A = W_{II}\)，实际上我们可以通过采样部分计算 \(B\)，而无需构造整个 \(W\)。 -
对 \(B\) 进行特征分解：
计算 \(B\) 的特征分解：\(B = U_B \Lambda_B U_B^T\)，其中 \(U_B \in \mathbb{R}^{l \times l}\) 正交，\(\Lambda_B\) 是特征值矩阵。则 \(W\) 的特征向量近似为 \(\tilde{U} = Q U_B \in \mathbb{R}^{n \times l}\)，特征值近似为 \(\Lambda_B\)。 -
截取前 \(k\) 个特征对：
由于我们只需要前 \(k\) 个最小的特征值对应的特征向量（对于拉普拉斯矩阵 \(L\)），我们从 \(\Lambda_B\) 中选取最小的 \(k\) 个特征值（实际上对于 \(W\)，我们可能需要最大的 \(k\) 个，因为 \(W\) 与 \(L\) 的特征值顺序相反，但原理相同），并取出对应的 \(\tilde{U}\) 中的 \(k\) 列。
第三步：在谱聚类中的具体应用步骤
-
输入：
- 数据点集 \(X = \{x_1, \dots, x_n\}\)，相似度函数 \(s(x_i, x_j)\)（如高斯核）。
- 聚类数 \(k\)。
- 采样列数 \(l\)（通常 \(l \approx 2k\) 或稍大，例如 \(l = k + 10\)）。
- 目标：近似计算归一化拉普拉斯矩阵 \(L\) 的最小的 \(k\) 个特征值对应的特征向量。
-
算法步骤：
a. 构造相似矩阵 \(W\)：
计算 \(W_{ij} = s(x_i, x_j)\)，但注意：我们不需要完整计算整个 \(W\)，只需在采样步骤中需要时计算部分。b. 随机采样 \(l\) 个列：
从 \(\{1,\dots,n\}\) 中均匀随机（或有偏，如按行范数采样）选取索引集 \(I\)。c. 计算子矩阵 \(C\) 和 \(A\)：
计算 \(C = W_{:,I}\)，即 \(C\) 的第 \(j\) 列是 \(W\) 的第 \(I_j\) 列。计算 \(A = W_{II} = C[I, :] \in \mathbb{R}^{l \times l}\)。d. 计算度矩阵 \(D\) 的近似：
计算 \(d_i = \sum_{j=1}^n W_{ij}\) 的对角矩阵 \(D\)。由于 \(W\) 没有完全计算，我们可以用 \(C\) 来近似：\(\tilde{d}_i = \sum_{j=1}^l C_{ij} \times (n/l)\)（基于采样列的缩放估计），或者更精确地，我们可以用 Nyström 近似 \(\tilde{W}\) 的行和来估计 \(d_i\)。e. 构造归一化拉普拉斯矩阵 \(L\) 的近似：
归一化拉普拉斯矩阵为 \(L = I - D^{-1/2} W D^{-1/2}\)。我们近似为：
\[ \tilde{L} = I - \tilde{D}^{-1/2} \tilde{W} \tilde{D}^{-1/2} \]
其中 $ \tilde{W} = C A^{\dagger} C^T $，$ \tilde{D} $ 是 $ \tilde{W} $ 的行和构成的对角矩阵。
f. 用 Nyström 方法近似 \(\tilde{L}\) 的特征分解：
对 \(\tilde{L}\) 应用上述 Nyström 步骤。但由于 \(\tilde{L}\) 是归一化的，我们可以直接在 \(W\) 上做近似，然后归一化。一个更稳定的方法是：
- 计算 \(A\) 的特征分解 \(A = U_A \Lambda_A U_A^T\)。
- 令 \(V = C U_A \Lambda_A^{-1/2}\)。
- 对 \(V\) 的每一行进行归一化：设 \(v_i\) 是 \(V\) 的第 \(i\) 行，计算 \(\bar{v}_i = v_i / \|v_i\|_2\)，形成矩阵 \(\bar{V}\)。
- 对 \(\bar{V}\) 进行QR分解：\(\bar{V} = Q R\)。
- 计算 \(B = Q^T L Q\)。由于 \(L\) 未知，但我们可以用 \(\tilde{L}\) 的近似形式计算。实际上，可以证明 \(B\) 可以用 \(A, C, \tilde{D}\) 表示，无需显式构造 \(L\)。
- 对 \(B\) 进行特征分解，得到特征值矩阵 \(\Lambda_B\) 和特征向量矩阵 \(U_B\)。
- 近似特征向量为 \(\tilde{U} = Q U_B\)。
g. 选取最小的 \(k\) 个特征值对应的特征向量：
从 \(\Lambda_B\) 中找出最小的 \(k\) 个特征值，取出对应的 \(\tilde{U}\) 的列，得到 \(U_k \in \mathbb{R}^{n \times k}\)。
h. 对 \(U_k\) 的行进行k-means聚类：
将 \(U_k\) 的每一行视为一个 \(k\) 维点，对这些点运行k-means算法，得到聚类标签。
第四步：误差分析与关键点
-
近似误差来源：
- 采样误差：由于只采样了 \(l\) 列，信息不完整。误差取决于矩阵的谱衰减特性。如果矩阵是快速低秩衰减的，则近似效果好。
- 数值误差：在计算 \(A^{\dagger}\) 或特征分解时的数值不稳定。可以通过在 \(\Lambda_A^{-1/2}\) 中截断小特征值（设置阈值）来稳定。
- 归一化误差：用 \(\tilde{D}\) 近似 \(D\) 引入的误差。
-
计算复杂度：
- 采样和计算 \(C\) 和 \(A\)：需要计算 \(O(nl)\) 个相似度（而不是 \(O(n^2)\)）。
- 对 \(A\) 进行特征分解：\(O(l^3)\)。
- QR分解和计算 \(B\)：\(O(nl^2)\)。
- 对 \(B\) 特征分解：\(O(l^3)\)。
- 总体：\(O(nl^2 + l^3)\)，当 \(l \ll n\) 时，远低于 \(O(n^3)\)。
-
改进采样策略：
- 均匀随机采样可能不是最好的。可以使用基于列范数的非均匀采样（重要性采样），以更高概率选择“重要”的列（即与其他列相关性高的列）。
- 也可以使用k-means++类似的策略，或基于杠杆得分的采样。
第五步：示例
假设 \(n=10000\), \(k=5\), 我们设置 \(l=20\)。
- 随机选取20个列索引 \(I\)。
- 计算相似度矩阵的对应部分：计算 \(C \in \mathbb{R}^{10000 \times 20}\)（每一列是第 \(I_j\) 列与所有点的相似度），以及 \(A = C[I, :] \in \mathbb{R}^{20 \times 20}\)。
- 计算 \(A\) 的特征分解 \(A = U_A \Lambda_A U_A^T\)，截断 \(\Lambda_A\) 中小于 \(10^{-8} \cdot \max(\Lambda_A)\) 的特征值。
- 计算 \(V = C U_A \Lambda_A^{-1/2} \in \mathbb{R}^{10000 \times r}\)（其中 \(r \leq 20\) 是截断后的秩）。
- 对 \(V\) 的每一行归一化得到 \(\bar{V}\)，对 \(\bar{V}\) 进行QR分解：\(\bar{V} = Q R\)。
- 计算 \(B = Q^T L Q\) 的近似（通过 \(\tilde{L}\) 的表达式，用 \(C, A, \tilde{D}\) 表示）。
- 计算 \(B\) 的特征分解，得到 \(U_B\) 和 \(\Lambda_B\)。
- 取 \(\Lambda_B\) 中最小的5个特征值对应的 \(U_B\) 的列，乘以 \(Q\) 得到 \(U_5 \in \mathbb{R}^{10000 \times 5}\)。
- 对 \(U_5\) 的10000行运行k-means聚类，分成5类。
总结：
随机化Nyström方法通过列采样，将大规模矩阵的特征分解问题转化为小规模矩阵的特征分解，极大地降低了计算和存储成本。在谱聚类中，它使我们能够处理大规模数据集，同时得到与精确谱聚类相近的聚类效果。关键参数是采样列数 \(l\) 和采样策略，这些需要根据具体问题和矩阵特性来调整。