经典多维尺度变换（CMDS）概述
经典多维尺度变换（cmds），又称主坐标分析（principal coordinate analysis, pcoa），是一种常用的降维技术，旨在将高维数据点映射到低维空间，同时尽可能保留原始数据点之间的距离关系。其核心思想是通过对距离矩阵进行双重中心化，然后进行特征分解，从而找到数据在低维空间中的最优表示。cmds广泛应用于数据可视化、模式识别和机器学习等领域。
CMDS算法的基本步骤如下：
- 构建距离矩阵D：计算所有数据点对之间的距离。如果输入是原始数据，通常使用欧氏距离；如果输入已经是距离矩阵，则直接使用。
- 构建中心化矩阵H：$H = I - \frac{1}{n} \mathbf{1}\mathbf{1}^T$，其中$I$是单位矩阵，$n$是数据点数量，$\mathbf{1}$是全1向量。
- 双重中心化平方距离矩阵B：$B = -\frac{1}{2} H D^{(2)} H$，其中$D^{(2)}$表示距离矩阵$D$中每个元素平方后的矩阵。
- 特征分解：对矩阵$B$进行特征分解，得到特征值和特征向量。
- 降维投影：选择最大的$n_dim$个正特征值及其对应的特征向量，构建低维嵌入。
处理距离矩阵中的无穷大值
在某些应用场景中，特别是当数据点表示图中的节点，而距离表示它们之间的路径长度时，两个不连通的节点之间的距离通常会被标记为无穷大（inf）。这在图论中是一个完全有效的概念，但在CMDS的数值计算中，无穷大值会导致严重问题。
原始CMDS算法在计算双重中心化平方距离矩阵$B$时，涉及$D^2$的操作。如果$D$中包含inf，那么$D^2$中的对应元素将变为inf，进而导致矩阵乘法H @ D**2 @ H的结果中包含inf或NaN（Not a Number）。随后的特征分解np.linalg.eigh(B)将无法处理这些非有限值，从而引发运行时错误。
为了使CMDS算法能够鲁棒地处理包含不连通点（即距离为inf）的场景，我们需要在计算$B$之前对距离矩阵进行预处理。
解决方案：替换无穷大值
最直接且有效的方法是识别距离矩阵中的所有无穷大值，并将其替换为一个巨大的、但有限的数值。这个替换值应该足够大，以反映不连通点之间“无限远”的语义，但又不能真正是inf，以避免数值计算错误。Python的numpy库提供了np.finfo(D.dtype).max，它能返回给定数据类型所能表示的最大有限浮点数，这通常是一个理想的替换值。
实现步骤：
- 检查输入距离矩阵D中是否存在无穷大值。
- 如果存在，则将所有inf值替换为np.finfo(D.dtype).max。
通过这种方式，我们既保留了不连通点距离“非常大”的语义，又避免了数值计算的崩溃。
示例代码
以下是修改后的CMDS函数，它集成了处理无穷大距离值的功能：
import numpy as np
from sklearn.metrics.pairwise import euclidean_distances
def cmds(X, n_dim, input_type='raw'):
"""
Classical(linear) multidimensional scaling (MDS)
Parameters
----------
X: (d, n) array or (n,n) array
input data. The data are placed in column-major order.
That is, samples are placed in the matrix (X) as column vectors
d: dimension of points
n: number of points
n_dim: dimension of target space
input_type: it indicates whether data are raw or distance
- raw: raw data. (n,d) array.
- distance: precomputed distances between the data. (n,n) array.
Returns
-------
Y: (n_dim, n) array. projected embeddings.
evals: (n_dim) eigen values
evecs: corresponding eigen vectors in column vectors
"""
if input_type == 'distance':
D = X
elif input_type == 'raw':
# Transpose X to (n, d) for euclidean_distances
Xt = X.T
D = euclidean_distances(Xt, Xt)
else:
raise ValueError("input_type must be 'raw' or 'distance'")
# 检查距离矩阵中是否存在无穷大值，并进行替换
if np.any(np.isinf(D)):
# 将inf值替换为该数据类型所能表示的最大有限浮点数
# 这样做可以避免在后续计算中因inf值导致错误，同时保留其“非常远”的语义
D[np.isinf(D)] = np.finfo(D.dtype).max
# Centering matrix
n = D.shape[0]
H = np.eye(n) - np.ones((n, n)) / n
# Double-center the distance matrix
# 注意：这里D**2是元素级的平方操作
B = -0.5 * H @ (D**2) @ H
# Eigen decomposition
evals, evecs = np.linalg.eigh(B)
# Sorting eigenvalues and eigenvectors in decreasing order
sort_indices = np.argsort(evals)[::-1]
evals = evals[sort_indices]
evecs = evecs[:, sort_indices]
# Selecting top n_dim eigenvectors
evecs_selected = evecs[:, :n_dim]
# Projecting data to the new space
# 确保特征值非负，因为它们理论上应代表方差
# 实际应用中，由于数值精度或非欧氏距离，可能出现微小负特征值，
# 但对于CMDS，通常只考虑正特征值。这里假设前n_dim特征值是有效的。
valid_evals = np.maximum(0, evals[:n_dim]) # 避免负数开方
Y = np.sqrt(np.diag(valid_evals)) @ evecs_selected.T
return Y, evals, evecs
注意事项
- 替换值的选择： np.finfo(D.dtype).max是一个稳健的选择，因为它确保了替换值在当前浮点数类型下是最大的有限数。如果使用一个任意大的常数，需要确保它足够大以区分连通点，但又不能过大导致其他数值溢出（尽管这种情况在现代浮点数系统中不常见）。
- 语义影响： 将inf替换为有限大值会稍微改变不连通点之间的“真实”距离，但对于CMDS的目标（在低维空间中近似距离关系），这种近似是合理的。它允许算法继续运行并提供一个可解释的嵌入，即使这些点在原始空间中是完全不连通的。
- 适用场景： 这种处理方法特别适用于距离矩阵中inf表示“不可达”或“不连通”的情况，例如在图论或网络分析中。如果inf在你的应用中具有其他特定的、需要严格区分的语义，可能需要考虑更复杂的处理策略。
- 负特征值： 在CMDS中，理论上所有特征值都应非负。然而，由于数值精度问题或输入距离矩阵并非严格欧氏距离（例如，经过inf替换后），可能会出现微小的负特征值。在投影步骤np.sqrt(np.diag(evals[:n_dim]))中，如果evals包含负数，np.sqrt会产生复数。通常做法是取max(0, eval)来避免复数，如示例代码所示。
总结
通过在CMDS算法中引入对距离矩阵中无穷大值的检测和替换机制，我们显著提升了算法的鲁棒性。这种方法使得CMDS能够有效处理包含不连通关系的数据集，从而扩展了其在复杂网络和图结构数据分析中的应用范围。在实际应用中，理解并妥善处理数据中的特殊值（如inf或NaN）是构建稳定、可靠数据分析流程的关键一环。
