Python中scipy.linalg.eig计算随机矩阵特征值不稳定原因

假设存在数据矩阵X，样本量num_samples = 1600，数据维度dim_data = 2，可通过RBF核构建规模为1600*1600的相似度矩阵S。对矩阵逐行做归一化处理：将行内所有元素乘以1 / 该行元素总和，即可得到方阵形式的右随机矩阵，理论上该矩阵存在值为1的特征值，对应元素全为1的常数特征向量。

将全1向量与该矩阵做乘积即可验证其确实为特征向量，但使用scipy.linalg.eig计算得到的、对应特征值1的特征向量仅为分段常数。
对随机生成的同规模数据转换得到的随机矩阵使用scipy.linalg.eig计算，均能稳定得到对应特征值1的常数特征向量。

核心问题：使用scipy.linalg.eig计算随机矩阵特征值时，哪些因素会引发数值不稳定问题？

RBF核计算函数 #

def kernel(sigma,X):
  """
  param sigma: variance
  param X: (num_samples,data_dim)
  """
  squared_norm = np.expand_dims(np.sum(X**2,axis=1),axis=1) + np.expand_dims(np.sum(X**2,axis=1),axis=0)-2*np.einsum('ni,mi->nm',X,X)
  return np.exp(-0.5*squared_norm/sigma**2)

行归一化生成随机矩阵函数 #

def normalize(array):
  degrees = []
  M = array.shape[0]
  for i in range(M):
    norm = sum(array[i,:])
    degrees.append(norm)
  degrees_matrix = np.diag(np.array(degrees))
  P = np.matmul(np.linalg.inv(degrees_matrix),array)
  return P

测试流程代码 #

#generate the data
points = np.linspace(0,4*np.pi,1600)
Z = np.zeros((1600,2))
Z[0:800,:] = np.array([2.2*np.cos(points[0:800]),2.2*np.sin(points[0:800])]).T
Z[800:,:] = np.array([4*np.cos(points[0:800]),4*np.sin(points[0:800])]).T

X = np.zeros((1600,2))
X[:,0] = np.where(Z[:,1] >= 0, Z[:,0] + .8 + params[1], Z[:,0] - .8 + params[2])
X[:,1] = Z[:,1] + params[0]

#create the stochastic matrix P
P = normalize(kernel(.05,X))

#inspect the eigenvectors
e,v = scipy.linalg.eig(P)
p = np.flip(np.argsort(e))
e = e[p]
v = v[:,p]

plot_array(v[:,0])

#check on synthetic data:
Y = np.random.normal(size=(1600,2))
P = normalize(kernel(Y))

#inspect the eigenvectors
e,v = scipy.linalg.eig(P)
p = np.flip(np.argsort(e))
e = e[p]
v = v[:,p]

plot_array(v[:,0])

对计算得到的特征向量和理论常数特征向量做偏差统计，结果如下：

[-1.36116641e-05 -1.36116641e-05 -1.36116641e-05 ...  5.44472888e-06
  5.44472888e-06  5.44472888e-06]
norm =  0.9999999999999999
max difference =  0.04986484253966891
max difference / element value -3663.3906291852545

构建核矩阵时sigma取值越小，排序后的特征值衰减越平缓：

• 当sigma=0.05时，scipy.linalg.eig输出的前4个特征值四舍五入后均为1，会出现特征向量计算精度不足的问题，得到的前5个特征向量为分段常数形态：
  
• 当sigma提升至0.5时，即可得到正常的常数特征向量，对应的前5个特征向量形态如下：
  

这个数值不稳定问题本质是近重复特征值引发的特征向量空间混淆：

• 当sigma取值极小（比如0.05）时，RBF核的衰减速度极快，只有距离非常近的点之间才有非零相似度。构造的两个同心圆被上下偏移后，四个分段区域的点跨区域距离远大于sigma，跨区域相似度几乎为0，整个相似度矩阵近似为分块对角矩阵，每个子块都是独立的右随机矩阵，各自对应一个值为1的特征值。浮点数精度下这几个特征值的差值小于计算误差阈值，scipy.linalg.eig无法区分这几个几乎相等的特征值，返回的是这几个特征向量张成空间中的任意一组正交基，也就是观测到的分段常数向量。
• 随机生成的高斯数据点在空间中均匀分布，不存在完全割裂的分块结构，最大特征值1和其余特征值之间有明显的谱间隔，因此特征值分解结果稳定，能正确返回常数特征向量。
• 当sigma提升到0.5时，核函数的作用范围足够覆盖跨区域的点距离，整个图结构完全连通，仅存在唯一的最大特征值1，因此可以正常计算得到常数特征向量。

小sigma场景下的解决方法：

1. 右随机矩阵对应特征值1的特征向量恒为全1向量，这是定义直接给出的性质，不需要通过通用特征值分解计算，直接使用即可。
2. 如果需要计算其余特征向量，可以换用针对大规模稀疏矩阵设计的特征值求解器（比如scipy.sparse.linalg.eigsh），指定只求解最大的k个特征对，避免多个近邻特征值引发的基混淆问题。

内容的提问来源于stack exchange，提问作者kiyopi