了解 Numpy 如何进行 SVD

计算科学 线性代数 Python 拉帕克
2021-12-22 15:44:01

我一直在使用不同的方法来计算矩阵的秩和矩阵方程组的解。我遇到了函数 linalg.svd。将此与我自己用高斯消元法求解系统的努力进行比较,它似乎更快更精确。我试图了解这是怎么可能的。

据我所知, linalg.svd 函数使用 QR 算法来计算矩阵的特征值。我知道这在数学上是如何工作的,但我不知道 Numpy 是如何做到如此快速且不损失太多精度的。

所以我的问题是:numpy.svd 函数是如何工作的,更具体地说,它如何快速准确地完成它(与高斯消除相比)?

3个回答

你的问题有很多问题。

不要使用高斯消元法(LU 分解)来计算矩阵的数值秩。 LU 因式分解在浮点运算中不可靠。相反,使用揭示秩的 QR 分解(例如xGEQPXxGEPQY在 LAPACK 中,其中 x 是 C、D、S 或 Z,尽管这些例程很难追踪;请参阅JedBrown 对相关问题的回答),或使用 SVD (奇异值分解,例如xGESDDor xGESVD,其中 x 又是 C、D、S 或 Z)。SVD 是一种更准确、更可靠的数值秩确定算法,但它需要更多的浮点运算。

但是,对于求解线性系统,LU 分解(带有部分旋转,这是 LAPACK 中的标准实现)在实践中非常可靠。在某些病理情况下,带有部分枢轴的 LU 分解是不稳定的(参见数值线性代数中的第 22 讲由 Trefethen 和 Bau 了解详细信息)。QR 分解是一种用于求解线性系统的更稳定的数值算法,这可能就是它为您提供如此精确结果的原因。但是,对于方阵,它比 LU 因式分解需要更多的浮点运算(我相信;JackPoulson 可能会纠正我)。对于矩形系统,QR 分解是更好的选择,因为它会为超定线性系统产生最小二乘解。SVD 也可用于求解线性系统,但它会比 QR 分解更昂贵。

janneb 是正确的,numpy.linalg.svd 是xGESDDLAPACK 中的包装器。奇异值分解分两个阶段进行。首先,将要分解的矩阵简化为双对角形式。LAPACK 中用于简化为双对角形式的算法可能是 Lawson-Hanson-Chan 算法,它确实在某一点使用了 QR 分解。Trefethen 和 Bau 在数值线性代数中的第 31 讲概述了这一过程。然后,xGESDD使用分治算法从双对角矩阵中计算奇异值和左右奇异向量。要了解此步骤的背景知识,您需要参考Golub 和 Van Loan 的Matrix Computations或 Jim Demmel 的Applied Numerical Linear Algebra

最后,您不应将奇异值与特征值混淆这两组数量并不相同。SVD 计算矩阵的奇异值。Cleve Moler 的Numerical Computing with MATLAB很好地概述了奇异值和特征值之间的差异通常,给定矩阵的奇异值与其特征值之间没有明显的关系,除非在正规矩阵的情况下,其中奇异值是特征值的绝对值。

由于您问题的措辞,我假设您的矩阵是正方形的。LAPACK 的 SVD 例程,例如zgesvd,对于方阵基本上分三个阶段进行:

  1. 隐式计算酉矩阵UAVA,作为 Householder 变换的乘积,使得一般矩阵A简化为实数上双对角矩阵B:=UAHAVA. 在该子程序退出时,Householder 向量(标准化,以便它们的第一个条目隐含地为一个)UAVA分别存储在的部分B在主对角线和超对角线的下方和右侧。这一步需要O(n3)工作。
  2. 用于计算实对称三对角矩阵的特征值分解的算法的变体用于计算双对角 SVD。尽管据我所知,最知名的实对称三对角 EVP (MRRR) 算法对于双对角 SVD 来说还不是稳定的,但这里有一个有趣的讨论。LAPACK 目前对双对角 SVD使用分而治之的方法。双对角 SVD 产量{UB,VB,Σ}这样B=UBΣVBH. 这一步需要O(n2)当 SVD 的 MRRR 稳定时工作,但目前需要尽可能多的O(n3)工作。
  3. 这一步很简单。自从UABVAH=A,A=(UAUB)Σ(VAVB)H,因此可以在应用隐式定义的 Householder 变换后显式形成 SVDUAVAUBVB, 分别。这一步需要O(n3)工作。

numpy.linalg.svd 是来自 LAPACK 的 {Z,D}GESDD 的包装器。反过来,LAPACK 是由一些世界上最重要的数值线性代数专家精心编写的。事实上,如果一个对该领域不熟悉的人能够成功击败 LAPACK(无论是速度还是准确性),这将是非常令人惊讶的。

至于为什么 QR 比高斯消除更好,那可能更适合https://scicomp.stackexchange.com/