numpy.linalg.eigh vs numpy.linalg.svd如何？

Question

问题描述

对于方阵，可以得到SVD

X= USV'

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

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

例程或numpy.linalg.eigh，以计算埃尔米特矩阵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反射。

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

1. 使用Bidiagonalize A dgebrd()
2. 使用分治法计算对角矩阵的SVD。 DBDSDC()
3. 使用通过dgebrd()应用dormbr()两次而返回的矩阵P和Q来还原对角线化，一次用于U，一次用于V。

尽管LAPACK执行的实际操作有很大不同，但是这些策略在全局上是相似的。这可能源于以下事实：计算通用矩阵A的SVD与执行对称矩阵A ^ TA的特征分解相似

关于lapack分治法SVD的准确性和性能，请参阅此SVD方法调查：

• 他们经常实现基于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 /sf/ask/2877638571/
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。。

也可以执行准确性检查。