numpy.linalg.eigh 与 numpy.linalg.svd 相比如何？

Question

问题描述

对于方阵，可以得到SVD

X= USV'

分解，通过简单地使用numpy.linalg.svd

u,s,vh = numpy.linalg.svd(X)

例程或numpy.linalg.eigh，计算厄密矩阵X'X和上的eig分解XX'

他们使用相同的算法吗？调用相同的 Lapack 例程？

在速度方面有什么区别吗？和稳定性？

Answer 1

不，他们不使用相同的算法，因为他们做不同的事情。它们有些相关，但也非常不同。让我们从您可以在 m x n 矩阵上执行 SVD 这一事实开始，其中 m 和 n 不需要相同。

取决于numpy的版本，你在做什么。以下是 lapack 中用于双精度的特征值例程：
http://www.netlib.org/lapack/explore-html/d9/d8e/group__double_g_eeigen.html
以及相应的 SVD 例程：
http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing.html

套路有差异。差异很大。如果您关心细节，它们在 fortran headers 中指定得非常好。在许多情况下，了解您面前的矩阵类型对于做出良好的例程选择是有意义的。矩阵是symmetric/hermitian？它是上对角线形式吗？是半正定的吗？ ...

运行时存在巨大差异。但根据经验，EIG 比 SVD 便宜。但这也取决于收敛速度，而收敛速度又在很大程度上取决于矩阵的条件数，换句话说，矩阵的不适定性如何，...

SVD 通常是非常稳健和缓慢的算法，经常用于反演、通过截断进行速度优化、主成分分析，实际上你正在处理的矩阵只是一堆乱七八糟的行；）

Answer 2

的确，numpy.linalg.svd和numpy.linalg.eigh调用的不是同一个Lapack例程。一方面，numpy.linalg.eigh 引用了 LAPACK 的 dsyevd() 而 numpy.linalg.svd 使用了 LAPACK 的 dgesdd()。

这些例程之间的共同点是使用 Cuppen 的分而治之算法，该算法首先设计用于解决三对角特征值问题。例如，dsyevd() 仅处理 Hermitian 矩阵并执行以下步骤，仅当需要特征向量时：

使用 DSYTRD()
通过 DSTEDC()
应用 DSYTRD() 使用 DORMTR() 报告的 Householder 反射。

相反，要计算 SVD，dgesdd() 执行以下步骤，在 job==A 的情况下（需要 U 和 VT）：

使用 dgebrd()
使用 DBDSDC()
使用 dgebrd() 返回的矩阵 P 和 Q 应用 dormbr() 两次，一次用于 U，一次用于 V

虽然 LAPACK 执行的实际操作非常不同，但策略在全球范围内是相似的。这可能源于这样一个事实，即计算一般矩阵 A 的 SVD 类似于执行对称矩阵 A^T.A.

的特征分解。

关于 lapack 分治 SVD 的准确性和性能，请参阅 This survey of SVD methods:

他们经常达到基于 QR 的 SVD 的准确性，尽管尚未得到证实
如果没有发生通货紧缩，最坏的情况是 O(n^3)，但通常证明比那更好
内存要求是矩阵大小的 8 倍，这可能会让人望而却步

关于对称特征值问题，复杂度为 4/3n^3（但通常证明比这更好），内存占用约为 2n^2 加上矩阵的大小。因此，如果您的矩阵是对称的，最好的选择可能是 numpy.linalg.eigh。

可以使用以下代码计算特定矩阵的实际复杂度：

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

# see 
def func_powerlaw(x, m, c):
    return np.log(np.abs( x**m * c))

import time

start = time.time()
print("hello")
end = time.time()
print(end - start)

timeev=[]
timesvd=[]
size=[]

for n in range(10,600):
    print n
    size.append(n)
    A=np.zeros((n,n))
    #populate A, 1D diffusion. 
    for j in range(n):
        A[j,j]=2.
        if j>0:
            A[j-1,j]=-1.
        if j<n-1:
            A[j+1,j]=-1.

    #EIG
    Aev=A.copy()
    start = time.time()
    w,v=np.linalg.eigh(Aev,'L')
    end = time.time()
    timeev.append(end-start)
    Asvd=A.copy()
    start = time.time()
    u,s,vh=np.linalg.svd(Asvd)
    end = time.time()
    timesvd.append(end-start)


poptev, pcov = curve_fit(func_powerlaw, size[len(size)/2:], np.log(timeev[len(size)/2:]),p0=[2.1,1e-7],maxfev = 8000)
print poptev

poptsvd, pcov = curve_fit(func_powerlaw, size[len(size)/2:], np.log(timesvd[len(size)/2:]),p0=[2.1,1e-7],maxfev = 8000)
print poptsvd

plt.figure()

fig, ax = plt.subplots()

plt.plot(size,timeev,label="eigh")
plt.plot(size,[np.exp(func_powerlaw(x, poptev[0], poptev[1])) for x in size],label="eigh-adjusted complexity: "+str(poptev[0]))

plt.plot(size,timesvd,label="svd")
plt.plot(size,[np.exp(func_powerlaw(x, poptsvd[0], poptsvd[1])) for x in size],label="svd-adjusted complexity: "+str(poptsvd[0]))


ax.set_xlabel('n')
ax.set_ylabel('time, s')

#plt.legend(loc="upper left")

ax.legend(loc="lower right")
ax.set_yscale("log", nonposy='clip')

fig.tight_layout()



plt.savefig('eigh.jpg')
plt.show()

对于这样的一维扩散矩阵，eigh优于svd，但实际复杂度相似，略低于n^3，大约n^2.5。

也可以检查准确性。

numpy.linalg.eigh 与 numpy.linalg.svd 相比如何？

how does numpy.linalg.eigh vs numpy.linalg.svd?

python

numpy

lapack

numerical

问题描述