经典多维尺度分析（CMDS）是一种降维技术，旨在将高维数据投影到低维空间，同时尽可能保留数据点之间的距离关系。然而，在实际应用中，当输入数据为距离矩阵时，可能会遇到矩阵中包含无穷值（inf）的情况。这种情况通常发生在图中存在不连通的点时，这些点之间的距离被设置为无穷大。原始的CMDS算法在处理包含无穷值的距离矩阵时会报错，因为后续的计算（例如，计算中心矩阵和特征值分解）无法处理无穷值。
为了解决这个问题，我们需要在CMDS算法中添加一个预处理步骤，将距离矩阵中的无穷值替换为一个足够大的有限值。这样做可以保证后续的计算能够正常进行，并且不会对最终的降维结果产生显著的影响。
下面是修改后的CMDS算法的Python代码：
import numpy as np
from sklearn.metrics 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':
Xt = X.T
D = euclidean_distances(Xt,Xt)
# Check for inf values in the distance matrix
if np.any(np.isinf(D)):
# Replace inf values with a large but finite value
D[np.isinf(D)] = np.finfo(D.dtype).max
# Centering matrix
H = np.eye(D.shape[0]) - np.ones(D.shape) / D.shape[0]
# Double-center the distance matrix
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 = evecs[:, :n_dim]
# Projecting data to the new space
Y = np.sqrt(np.diag(evals[:n_dim])) @ evecs.T
return Y, evals, evecs代码解释：
- 导入必要的库: numpy 用于数值计算，sklearn.metrics.euclidean_distances 用于计算欧氏距离（如果输入类型为原始数据）。
-
cmds 函数: 该函数实现了CMDS算法。
-
输入参数:
- X: 输入数据，可以是原始数据矩阵 (d, n) 或距离矩阵 (n, n)。
- n_dim: 目标空间的维度。
- input_type: 指示输入数据类型，可以是 'raw'（原始数据）或 'distance'（距离矩阵）。
-
数据预处理:
- 如果 input_type 为 'raw'，则计算原始数据的欧氏距离矩阵。
- 关键修改： 检查距离矩阵 D 中是否存在无穷值 (np.isinf(D))。如果存在，则将这些无穷值替换为当前数据类型能表示的最大有限值 (np.finfo(D.dtype).max)。 np.finfo(D.dtype).max 会根据 D 的数据类型（例如 float64）返回该类型可以表示的最大值。
- 计算中心矩阵: H = np.eye(D.shape[0]) - np.ones(D.shape) / D.shape[0]。
- 双重中心化距离矩阵: B = -0.5 * H @ D**2 @ H。
- 特征值分解: evals, evecs = np.linalg.eigh(B。
- 排序特征值和特征向量: 按照特征值降序排列。
- 选择前 n_dim 个特征向量: 用于投影数据。
- 投影数据到新空间: Y = np.sqrt(np.diag(evals[:n_dim])) @ evecs.T。
- 返回结果: 投影后的嵌入 Y，特征值 evals，和特征向量 evecs。
-
输入参数:
使用示例：
import numpy as np
# 创建一个包含无穷值的距离矩阵
D = np.array([[0, 1, np.inf],
[1, 0, 2],
[np.inf, 2, 0]])
# 设置目标维度
n_dim = 2
# 使用修改后的CMDS算法进行降维
Y, evals, evecs = cmds(D, n_dim, input_type='distance')
# 打印结果
print("Projected embeddings:\n", Y)
print("Eigenvalues:\n", evals)
print("Eigenvectors:\n", evecs)注意事项：
- 将无穷值替换为一个足够大的有限值可能会对降维结果产生一定的影响，尤其是在图中存在大量不连通的点时。因此，在实际应用中，需要根据具体情况选择合适的替换值。一种更高级的策略是使用图论方法来处理不连通的图，例如，将不连通的图分解成多个连通的子图，然后分别对每个子图进行降维。
- 该代码使用了numpy和scikit-learn库。请确保已经安装了这些库。
总结：
通过在CMDS算法中添加一个预处理步骤，将距离矩阵中的无穷值替换为一个足够大的有限值，可以有效地解决CMDS算法在处理包含无穷值的距离矩阵时遇到的问题。这种方法简单易行，并且可以保证算法的正常运行。在实际应用中，需要根据具体情况选择合适的替换值，以获得最佳的降维效果。
