Python奇异值分解

原创已于 2023-07-29 13:15:39 修改 · 1.4k 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#python #numpy #scipy #svd #奇异值分解

于 2023-02-17 08:00:00 首次发布

scipy 专栏收录该内容

71 篇文章

订阅专栏

文章介绍了奇异值分解的基本原理，包括特征值和奇异值的关系，以及非方阵的SVD表示。在numpy和scipy库中，都提供了SVD的实现。通过对比测试，展示了两者在计算效率上的差异，numpy.linalg.svd比scipy.linalg.svd慢，而后者可以通过设置lapack_driver参数优化性能。

文章目录

基本原理

当 $A$ 是方阵时，可以很容易地进行特征分解： $A=W\Sigma W^{-1}$ ，其中 $\Sigma$ 是 $A$ 的特征值组成的对角矩阵。如果 $W$ 由标准正交基组成，则 $W^{-1}=W^T$ ，特征分解可进一步写成 $W^T\Sigma W$ 。

然而，当 $A$ 不是方阵时，情况大不一样了，但仍然可以将 $A$ 表示成 $A=U\Sigma V^T$ 的形式，其中 $\Sigma$ 也是对角矩阵，对角线上的每个元素被称作奇异值。

奇异值的求解过程和特征值息息相关，因为把 $A$ 变成方阵很简单，只要乘以转置就行。故令 $L=AA^T$ ， $R=A^TA$ ，则 $L, R$ 都可以求特征值 $\lambda_i$ 和特征向量，其中 $L$ 的特征向量为 $A$ 的左奇异向量， $R$ 的特征向量为右奇异向量。对应的奇异值 $\sigma_i=\sqrt{\lambda_i}$ 。

numpy的实现

numpy.linalg中提供了奇异值分解函数svd，参数为

svd(a, full_matrices=True, compute_uv=True, hermitian=False)

其中

a 待分解矩阵，维度为**(M, N)**
full_matrices 若为True，则U, Vh分别为**(M, M)和(N, N)；否则分别为(M, K), (K, N)**，K为M, N中较小的那个
compute_uv 如果为False则不计算U, Vh
hermitian 为True时，表示处理的是实对称矩阵

scipy实现

scipy.linalg中也提供了奇异值分解函数svd，其参数为

svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd')

其中与numpy.linalg相同的参数，其意义也相同，不相同的部分，各参数含义如下

overwrite_a 如果为True，则直接对a进行修改
check_finite 如果为True，则进行有限性检查
lapack_driver SVD分解的方法，有两个选择
- ‘gesdd’ 效率更高
- ‘gesvd’ 此为Matlab和R中使用的方法

其返回值即 $\Sigma, V^T$ 。

scipy.linalg还提供了两个和SVD相关的函数，svdvals(a)用于求a的奇异值；diagsvd(s, M, N)通过s, M, N，创建一个 $\Sigma$ 矩阵。

对比测试

下面测试一下svd

import numpy as np
import scipy.linalg as sl
a = np.random.rand(5,5)
u1, s1, vh1 = sl.svd(a)
u2, s2, vh2 = np.linalg.svd(a)
print(s1)
# [2.63698545 0.94063722 0.36159198 0.21052102 0.19014115]
print(s1-s2)
# [ 0.0  0.0  1.11022302e-16 -2.77555756e-17 0.0]

numpy和scipy的结果是几乎相同的，下面测试一下不同方法进行奇异值分解的时间

from timeit import timeit
a = np.random.rand(1000,1000)
timeit(lambda:sl.svd(a), number=10)
# 1.870287900000001
timeit(lambda:np.linalg.svd(a), number=10)
# 13.355788999999998
timeit(lambda:sl.svd(a, lapack_driver='gesvd'), number=10)
# 3.873418600000001