利用奇异值分解（SVD）求解线性最小二乘问题：数值稳定性与实践优化

import numpy as np
from scipy import linalg

# 1. 数据准备
np.random.seed(123)
# 创建一个接近秩亏损的矩阵A，以模拟小奇异值的情况
# 通过给一个秩为1的矩阵添加微小噪声，使其成为一个病态但满秩的矩阵
v = np.random.rand(4)
A = v[:,None] * v[None,:] + np.random.rand(4,4) * 1e-3 
b = np.random.randn(4)

print("--- 线性最小二乘问题求解对比 ---")

# 2. 对比方法一：通过正规方程组求解（可能存在数值不稳定）
try:
    x_normal_eq = linalg.inv(A.T @ A) @ A.T @ b
    l2_normal_eq = linalg.norm(A @ x_normal_eq - b)
    print(f"正规方程组 (手动实现) L2范数: {l2_normal_eq:.10f}")
except linalg.LinAlgError:
    print("正规方程组 (手动实现) 求解失败：矩阵奇异或接近奇异。")
    l2_normal_eq = np.nan

# 3. 对比方法二：使用scipy.linalg.lstsq（推荐的标准方法）
# 这是一个经过高度优化和数值稳定的实现，通常作为基准
x_lstsq = linalg.lstsq(A, b)[0]
l2_lstsq = linalg.norm(A @ x_lstsq - b)
print(f"scipy.linalg.lstsq L2范数: {l2_lstsq:.10f}")

# 4. 问题SVD实现：未处理小奇异值
# 这个函数直接对所有奇异值求逆，可能导致数值爆炸
def direct_ls_svd_problematic(A_matrix, b_vector):
    U, S, Vt = linalg.svd(A_matrix, full_matrices=False)
    # 直接对S中的每个奇异值取倒数，如果S中包含极小值，会产生巨大误差
    S_inv = np.diag(1/S)
    x_hat = Vt.T @ S_inv @ U.T @ b_vector
    return x_hat

x_svd_problematic = direct_ls_svd_problematic(A, b)
l2_svd_problematic = linalg.norm(A @ x_svd_problematic - b)
print(f"SVD (未处理小奇异值) L2范数: {l2_svd_problematic:.10f}")

# 原始SVD输出示例（可能因随机种子略有不同）：
# 正规方程组 (手动实现) L2范数: 2.9286130558
# scipy.linalg.lstsq L2范数: 2.9286130558
# SVD (未处理小奇异值) L2范数: 6.8305500190 (或更高)

def direct_ls_svd_optimized(A_matrix, b_vector, rcond=1e-7):
  """
  使用奇异值分解（SVD）求解线性最小二乘问题，并处理小奇异值以提高数值稳定性。

  参数:
  A_matrix (np.ndarray): 设计矩阵。
  b_vector (np.ndarray): 目标向量。
  rcond (float): 相对条件数阈值。小于 rcond * max(S) 的奇异值将被视为零。

  返回:
  np.ndarray: 最小二乘解向量 x_hat。
  """
  # 计算