利用SVD解决线性方程组：数值稳定性的关键优化

import numpy as np
from scipy import linalg

np.random.seed(123)
v = np.random.rand(4)
A = v[:,None] * v[None,:] # 生成一个秩为1的矩阵，因此会有多个接近0的奇异值
b = np.random.randn(4)

# 方法1: 使用正规方程组（通常不推荐，数值不稳定）
x_manual = linalg.inv(A.T.dot(A)).dot(A.T).dot(b)
l2_manual = linalg.norm(A.dot(x_manual) - b)
print("manually (Normal Equations): ", l2_manual)

# 方法2: 使用scipy.linalg.lstsq (推荐)
x_lstsq = linalg.lstsq(A, b)[0]
l2_lstsq = linalg.norm(A.dot(x_lstsq) - b)
print("scipy.linalg.lstsq: ", l2_lstsq)

# 方法3: 初始自定义SVD实现 (存在问题)
def direct_ls_svd_problematic(A_matrix, b_vector):
  # 注意：此函数在原始问题中期望x是输入，y是输出，但这里我们将其调整为A, b
  # calculate the economy SVD for the data matrix A_matrix
  U,S,Vt = linalg.svd(A_matrix, full_matrices=False)

  # 尝试直接计算伪逆，但未处理接近零的奇异值
  # x_hat = Vt.T @ linalg.inv(np.diag(S)) @ U.T @ b_vector # 这种方式对S=0的值会报错
  # 更常见的SVD解法形式
  S_inv_diag = np.diag(1/S) # 如果S中有0或接近0的值，这里会出问题
  x_hat = Vt.T @ S_inv_diag @ U.T @ b_vector
  return x_hat

# 运行问题代码
# x_svd_problematic = direct_ls_svd_problematic(A, b) # 可能会因除以零而失败
# 为了演示问题，我们直接使用原始问题中的SVD代码，它没有直接计算伪逆，但仍会受到小奇异值影响
# 原始问题中的 direct_ls_svd 函数返回的是残差，这里需要修改以返回x_hat
def direct_ls_svd_original(A_matrix, b_vector):
    U, S, Vt = linalg.svd(A_matrix, full_matrices=False)
    # 原始代码中直接使用 S 参与计算，但未过滤
    # x_hat = Vt.T @ linalg.inv(np.diag(S)) @ U.T @ b_vector # 原始问题中的实现
    # 调整为更常见的SVD最小二乘解形式
    S_inv = np.diag(1.0 / S) # 这里是潜在的数值问题来源
    x_hat = Vt.T @ S_inv @ U.T @ b_vector
    return x_hat

try:
    x_svd_original = direct_ls_svd_original(A, b)
    l2_svd_original = linalg.norm(A.dot(x_svd_original) - b)
    print("svd (original problematic): ", l2_svd_original)
except np.linalg.LinAlgError as e:
    print(f"svd (original problematic) failed: {e}")
except RuntimeWarning as e:
    print(f"svd (original problematic) warning: {e}")


# 方法4: 使用scipy.linalg.solve (针对A.T@A可逆的情况)
x_solve = linalg.solve(A.T@A, A.T@b)
l2_solve = linalg.norm(A.dot(x_solve) - b)
print("scipy.linalg.solve: ", l2_solve)

print("\n--- 原始代码运行结果 ---")
print("manually (Normal Equations): ", l2_manual)
print("scipy.linalg.lstsq: ", l2_lstsq)
# 假设 direct_ls_svd_original 运行成功，这里打印其结果
# print("svd (original problematic): ", l2_svd_original) # 如果运行失败则不打印
print("scipy.linalg.solve: ", l2_solve)

# 比较l2_manual和l2_lstsq
print("np.allclose(l2_manual, l2_lstsq, rtol=1.3e-1):", np.allclose(l2_manual, l2_lstsq, rtol=1.3e-1))

def direct_ls_svd_corrected(A_matrix, b_vector, rcond=1e-7):
  """
  使用SVD计算线性最小二乘解，并进行奇异值过滤以提高数值稳定性。

  参数:
    A_matrix (np.ndarray): 系数矩阵 A。
    b_vector (np.ndarray): 右侧向量 b。
    rcond (float): 相对条件数。小于 rcond * max(S) 的奇异值将被视为零。

  返回:
    np.ndarray: 最小二乘解 x_hat。
  """
  # 计算经济型SVD
  U, S, Vt = linalg.svd(A_matrix, full_matrices=False)

  # 构建一个掩码，过滤掉小于 rcond * max(S) 的奇异值
  # 这些奇异值被认为是数值上的零，它们的倒数会导致不稳定
  mask = (abs(S) / np.max(abs(S))) > rcond

  # 仅保留有效的奇异值及其对应的U和Vt部分
  U_filtered, S_filtered, Vt_filtered = U[:, mask], S[mask], Vt[mask, :]

  # 验证重构的A是否接近原始A（可选，用于调试）
  # assert np.allclose(U_filtered @ np.diag(S_filtered) @ Vt_filtered, A_matrix)

  # 计算最小二乘解 x_hat = V_filtered.T @ Sigma_filtered_inv @ U_filtered.T @ b_vector
  # 这里使用更数值稳定的形式：x_hat = V_filtered.T @ ((U_filtered.T @ b_vector) / S_filtered)
  x_hat = Vt_filtered.T @ ((U_filtered.T @ b_vector) / S_filtered)

  return x_hat

# 使用修正后的SVD函数
x_svd_corrected = direct_ls_svd_corrected(A, b)
l2_svd_corrected = linalg.norm(A.dot(x_svd_corrected) - b)
print("\nsvd (corrected with filtering): ", l2_svd_corrected)

# 再次比较
print("np.allclose(l2_lstsq, l2_svd_corrected, rtol=1.3e-7):", np.allclose(l2_lstsq, l2_svd_corrected, rtol=1.3e-7))