svd奇异值分解python代码

### 使用 NumPy 和 SciPy 实现奇异值分解 (SVD) 在 Python 中，`numpy` 和 `scipy` 是常用的科学计算库，它们均提供了用于执行奇异值分解（Singular Value Decomposition, SVD）的功能。 #### 利用 NumPy 的 SVD 函数 以下是通过 `numpy.linalg.svd` 方法实现的一个简单例子： ```python import numpy as np # 定义一个矩阵 A A = np.array([[0, 1], [1, 1], [1, 0]]) # 执行 SVD 分解 u, s, vt = np.linalg.svd(A, full_matrices=True) print("U 矩阵:\n", u, "\n形状:", u.shape) print("\nΣ 向量:\n", s, "\n形状:", s.shape) print("\nV 转置矩阵:\n", vt, "\n形状:", vt.shape) ``` 此代码片段展示了如何对给定的矩阵 \( A \) 进行 SVD 分解，并打印出三个组成部分：\( U \)，\( Σ \)，以及 \( V^T \)[^1]。 #### 利用 SciPy 的 SVD 函数 同样可以借助 `scipy.linalg.svd` 来完成相同的操作。以下是一个示例： ```python from scipy import linalg import numpy as np # 定义一个矩阵 B B = np.random.rand(3, 2) # 执行 SVD 分解 u, s, vh = linalg.svd(B, full_matrices=False) print("U 矩阵:\n", u, "\n形状:", u.shape) print("\nΣ 向量:\n", s, "\n形状:", s.shape) print("\nVH 矩阵:\n", vh, "\n形状:", vh.shape) ``` 这段代码说明了当设置参数 `full_matrices=False` 时，返回的是缩减形式的 \( U \) 和 \( VH \) 矩阵[^3]。 #### 性能比较 对于大规模数据集而言，选择合适的工具至关重要。例如，在处理大小为 \( 1000 \times 1000 \) 的随机矩阵时，`scipy.linalg.svd` 显著优于 `numpy.linalg.svd`，尤其是在指定 LAPACK 驱动程序的情况下[^2]。 ```python from timeit import timeit from scipy import linalg as sl import numpy as np a = np.random.rand(1000, 1000) t_numpy = timeit(lambda: np.linalg.svd(a), number=10) t_scipy_default = timeit(lambda: sl.svd(a), number=10) t_scipy_gesvd = timeit(lambda: sl.svd(a, lapack_driver='gesvd'), number=10) print(f"Numpy svd 时间: {t_numpy:.4f} 秒") print(f"Scipy 默认 svd 时间: {t_scipy_default:.4f} 秒") print(f"Scipy gesvd svd 时间: {t_scipy_gesvd:.4f} 秒") ``` 以上代码可用于评估不同方法之间的效率差异[^2]。 --- ### 图像压缩应用案例 作为实际应用场景之一，可以通过保留较少数量的主要特征来近似原始图像，从而达到降低存储需求的目的。如下所示： ```python from scipy import linalg import numpy as np import cv2 as cv import matplotlib.pyplot as plt # 加载并转换灰度图 img = cv.imread('example_image.jpg', cv.IMREAD_GRAYSCALE) # 计算 SVD u, s, vh = linalg.svd(img, full_matrices=False) def reconstruct(k): """重建图像""" return (u[:, :k] @ np.diag(s[:k]) @ vh[:k, :]).clip(0, 255).astype(np.uint8) plt.figure(figsize=(12, 6)) for i, k in enumerate([5, 50, 100]): recon_img = reconstruct(k) plt.subplot(1, 3, i+1) plt.title(f'K={k}') plt.imshow(recon_img, cmap='gray') plt.axis('off') plt.tight_layout() plt.show() ``` 该脚本演示了利用不同程度的截断重构后的视觉效果变化情况[^3]。 ---