特征分解提取特征向量及后续SVD计算异常问题求助

问题背景 #

执行以下矩阵运算流程时，无法得到正确的正特征值对应U矩阵，且未提取到负特征向量：

1. 对m×n矩阵A做SVD分解，得到$A=USV^T$；
2. 构造对称矩阵$U'=\frac{1}{2}(U+U^T)$；
3. 对m×m的U'做特征分解；
4. 提取正特征值对应的特征向量构成矩阵X（m×k）；
5. 对$X^TX$做SVD分解，获取属于SO(n)的U矩阵；
6. 针对负特征值重复步骤4、5。

附实现代码：

import numpy as np

def svd(A):
    U, S, VT = np.linalg.svd(A, full_matrices=True)
    return U, S, VT

def symmetric(U):
    U_symmetric = 0.5 * (U + U.T)
    return U_symmetric

def eigenvalue_decomposition(U):
    eigenvalues, eigenvectors = np.linalg.eig(U)
    return eigenvalues, eigenvectors

def extract_positive_eigenvectors(eigenvalues, eigenvectors):
    positive_indices = np.where(eigenvalues > 0)
    if len(positive_indices[0]) == 0:
        return None
    X = eigenvectors[:, positive_indices]
    return X

def extract_negative_eigenvectors(eigenvalues, eigenvectors):
    negative_indices = np.where(eigenvalues < 0)
    if len(negative_indices[0]) == 0:
        return None
    X = eigenvectors[:, negative_indices]
    return X

A = np.array([[1, 2, 2], [0, 1, 2]])

U, _, _ = svd(A)
U_symmetric = symmetric(U)
eigenvalues, eigenvectors = eigenvalue_decomposition(U_symmetric)

# extract positive eigenvectors (if any)
X_positive = extract_positive_eigenvectors(eigenvalues, eigenvectors)

if X_positive is not None:
    X_positive_2d = X_positive.reshape(X_positive.shape[0], -1)
    U_positive, _, _ = svd(np.dot(X_positive_2d.T, X_positive_2d))
else:
    U_positive = None

print("Matrix U from SVD of XTX for positive eigenvectors:")
print(U_positive)

# extract negative eigenvectors (if any)
X_negative = extract_negative_eigenvectors(eigenvalues, eigenvectors)

if X_negative is not None:
    X_negative_2d = X_negative.reshape(X_negative.shape[0], -1)
    U_negative, _, _ = svd(np.dot(X_negative_2d.T, X_negative_2d))
else:
    U_negative = None

print("Matrix U from SVD of XTX for negative eigenvectors:")
print(U_negative)

排查方向与优化建议 #

一、数值精度与特征值判断问题 #

• 浮点数误差导致的误判：直接用eigenvalues > 0会忽略因数值计算误差产生的极小正/负特征值，甚至对称矩阵的特征值可能出现极小虚部。解决方式：
  1. 先对特征值取实部：eigenvalues = np.real(eigenvalues)，消除虚部干扰；
  2. 设置精度阈值，比如positive_indices = np.where(eigenvalues > 1e-10)，避免将接近0的极小值误判为0。
• 特征分解函数选择错误：np.linalg.eig是通用特征分解函数，对对称矩阵可能返回带虚部的结果。改用np.linalg.eigh——该函数专门针对对称/厄米矩阵优化，返回纯实特征值和特征向量，结果更稳定。

二、负特征向量未提取的原因 #

原矩阵U是正交矩阵，$U'=\frac{1}{2}(U+U^T)$的特征值范围在[-1,1]之间。若U非对称，U'的负特征值绝对值可能极小，被严格的<0判断过滤。加入精度阈值后（比如eigenvalues < -1e-10），可有效提取真实的负特征值对应向量。

三、$X^TX$的SVD结果不符合SO(n)要求 #

• SO(n)要求矩阵是行列式为1的正交矩阵，但np.linalg.svd返回的U矩阵行列式可能为-1。需额外判断并修正：

  def adjust_to_so(U):
      if np.linalg.det(U) < 0:
          U[:, -1] *= -1  # 翻转最后一列，将行列式转为1
      return U
  

• $X^TX$是对称半正定矩阵，其SVD等价于特征分解，直接用np.linalg.eigh获取正交特征向量矩阵，效率更高且结果更稳定，无需冗余的SVD计算。

四、代码优化版本 #

import numpy as np

def svd(A):
    U, S, VT = np.linalg.svd(A, full_matrices=True)
    return U, S, VT

def symmetric(U):
    return 0.5 * (U + U.T)

def eigenvalue_decomposition(U):
    # 对称矩阵专用特征分解，返回实值结果并按特征值绝对值降序排序
    eigenvalues, eigenvectors = np.linalg.eigh(U)
    idx = np.argsort(np.abs(eigenvalues))[::-1]
    return eigenvalues[idx], eigenvectors[:, idx]

def extract_positive_eigenvectors(eigenvalues, eigenvectors, tol=1e-10):
    eigenvalues = np.real(eigenvalues)
    positive_indices = np.where(eigenvalues > tol)
    if len(positive_indices[0]) == 0:
        return None
    return eigenvectors[:, positive_indices].reshape(eigenvectors.shape[0], -1)

def extract_negative_eigenvectors(eigenvalues, eigenvectors, tol=1e-10):
    eigenvalues = np.real(eigenvalues)
    negative_indices = np.where(eigenvalues < -tol)
    if len(negative_indices[0]) == 0:
        return None
    return eigenvectors[:, negative_indices].reshape(eigenvectors.shape[0], -1)

def get_so_matrix(XTX):
    _, U = np.linalg.eigh(XTX)
    if np.linalg.det(U) < 0:
        U[:, -1] *= -1
    return U

# 测试流程
A = np.array([[1, 2, 2], [0, 1, 2]])
U, _, _ = svd(A)
U_symmetric = symmetric(U)
eigenvalues, eigenvectors = eigenvalue_decomposition(U_symmetric)

# 正特征值处理
X_positive = extract_positive_eigenvectors(eigenvalues, eigenvectors)
U_positive = get_so_matrix(X_positive.T @ X_positive) if X_positive is not None else None
print("Matrix U from positive eigenvectors (SO(n)):")
print(U_positive)

# 负特征值处理
X_negative = extract_negative_eigenvectors(eigenvalues, eigenvectors)
U_negative = get_so_matrix(X_negative.T @ X_negative) if X_negative is not None else None
print("\nMatrix U from negative eigenvectors (SO(n)):")
print(U_negative)

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