利用奇异值分解（SVD）求解线性方程组的稳健方法

import numpy as np
from scipy import linalg

np.random.seed(123)
v = np.random.rand(4)
A = v[:,None] * v[None,:] # A is a rank-1 matrix, leading to small singular values
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, b):
  U, S, Vt = linalg.svd(A, full_matrices=False)
  # 原始问题代码，直接计算伪逆
  # x_hat = Vt.T @ linalg.inv(np.diag(S)) @ U.T @ b # 错误写法，应为 S 的倒数
  # 更准确的伪逆计算应为 (U.T @ b) / S
  x_hat = Vt.T @ ((U.T @ b) / S) # 即使这样，仍可能因S中极小值导致不稳定
  return x_hat

x_svd_problematic = direct_ls_svd_problematic(A, b)
l2_svd_problematic = linalg.norm(A.dot(x_svd_problematic) - b)
print("svd (problematic implementation): ", l2_svd_problematic)

# 4. 使用 scipy.linalg.solve (针对方阵的精确解，此处用于正规方程)
x_solve = linalg.solve(A.T@A, A.T@b)
l2_solve = linalg.norm(A.dot(x_solve) - b)
print("scipy.linalg.solve (normal equations): ", l2_solve)

print("\nComparison of L2 norms:")
print(f"Manual (normal equations): {l2_manual}")
print(f"scipy.linalg.lstsq: {l2_lstsq}")
print(f"SVD (problematic): {l2_svd_problematic}")
print(f"scipy.linalg.solve (normal equations): {l2_solve}")

# 示例输出可能如下：
# manually (normal equations):  2.9751344995811313
# scipy.linalg.lstsq:  2.9286130558050654
# svd (problematic implementation):  6.830550019041984
# scipy.linalg.solve (normal equations):  2.928613055805065

def direct_ls_svd(A, b, rcond=1e-7):
  """
  使用SVD稳健地求解线性最小二乘问题 Ax = b。

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

  返回:
    np.ndarray: 最小二乘解 x_hat。
  """
  U, S, Vt = linalg.svd(A, full_matrices=False)

  # 创建一个掩码，过滤掉接近零的奇异值
  # 任何小于 rcond * max(abs(S)) 的奇异值都将被视为零
  m = (abs(S) / np.max(abs(S))) > rcond

  # 根据掩码裁剪 U, S, Vt
  # 仅保留对解有贡献的奇异值和对应的向量
  U_filtered, S_filtered, Vt_filtered = U[:, m], S[m], Vt[m, :]

  # 验证过滤后的矩阵是否能近似重构原始矩阵 A
  # assert np.allclose(U_filtered @ np.diag(S_filtered) @ Vt_filtered, A, atol=1e-5) # 注意：对于病态矩阵，可能无法精确重构

  # 求解 x_hat = V_filtered.T @ diag(1/S_filtered) @ U_filtered.T @ b
  # 等价于 Vt_filtered.T @ ((U_filtered.T @ b) / S_filtered)
  x_hat = Vt_filtered.T @ ((U_filtered.T @ b) / S_filtered)

  return x_hat

# 使用修正后的 SVD 求解器
x_svd_corrected = direct_ls_svd(A, b)
l2_svd_corrected = linalg.norm(A.dot(x_svd_corrected) - b)
print("svd (corrected implementation): ", l2_svd_corrected)

print("\nComparison after correction:")
print(f"scipy.linalg.lstsq: {l2_lstsq}")
print(f"SVD (corrected): {l2_svd_corrected}")
print(f"Are corrected SVD and scipy.linalg.lstsq L2 norms close? {np.allclose(l2_lstsq, l2_svd_corrected, rtol=1e-6)}")

# 示例输出可能如下：
# svd (corrected implementation):  2.928613055805065
#
# Comparison after correction:
# scipy.linalg.lstsq: 2.9286130558050654
# SVD (corrected): 2.928613055805065
# Are corrected SVD and scipy.linalg.lstsq L2 norms close? True