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

Question

问题描述

对于方阵，可以得到SVD

X= USV'

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

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

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

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

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

Answer 1

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

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

1. 使用 DSYTRD() 将矩阵化简为三对角形式

2. 通过 DSTEDC() 使用分而治之算法计算三对角矩阵的特征向量

3. 使用 DORMTR() 应用 DSYTRD() 报告的 Householder 反射。

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

1. 使用 dgebrd() 将 A 双对角化
2. 使用 DBDSDC()
使用分而治之算法计算双对角矩阵的 SVD
3. 使用 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 https://stackoverflow.com/questions/41109122/fitting-a-curve-to-a-power-law-distribution-with-curve-fit-does-not-work
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。

也可以执行准确性检查。