怎么会呢纽比.利纳格.埃与纳比·利纳格·svd？

Question

问题描述

对于一个方阵，可以得到奇异值分解

X= USV'

分解，只需使用numpy.linalg.svd公司

^{pr2}$

常规或纽比.利纳格.埃，计算Hermitian矩阵X'X和XX'

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

在速度方面有什么不同吗？稳定性呢？在

Answer 1

实际上，numpy.linalg.svd和{}并不调用相同的Lapack例程。一方面，^{}表示LAPACK的dsyevd()，而{}使用LAPACK的{}。在

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

1. 使用DSYTRD（）将矩阵缩减为三对角形式

2. 通过DSTEDC（）

3. 使用DORMTR（）应用DSYTRD（）报告的Householder反射。

相反，为了计算SVD，^{}执行以下步骤，在job==A（需要U和VT）的情况下：

1. 使用dgebrd()对A进行双对角化
2. 用DBDSDC()分治算法计算双对角矩阵的奇异值分解
3. 使用dgebrd()两次应用dormbr()返回的矩阵P和Q来恢复双对角化，一次用于U，一次用于V

虽然LAPACK执行的实际操作非常不同，但策略在全局上是相似的。这可能源于这样一个事实：计算一般矩阵a的奇异值分解类似于对对称矩阵a^T.a进行特征分解

关于lapack分治SVD的精度和性能，请参见This survey of SVD methods：

• 他们经常获得基于QR的奇异值分解（SVD）的精度，尽管它没有被证明。在
• 如果没有发生通货紧缩，最坏的情况是O（n^3），但事实往往比这更好。在
• 内存需求是矩阵大小的8倍，这可能会让人望而却步。在

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

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

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，但实际复杂度相似，略低于n3，类似于n2.5。

也可以检查精确度。在

Answer 2

不，他们不使用相同的算法，因为他们做不同的事情。他们有点相关，但也有很大的不同。让我们从这样一个事实开始，您可以对m x n矩阵进行SVD，其中m和{}不需要相同。在

取决于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头文件中指定得非常好。在很多情况下，找出你面前的矩阵是什么样的是有意义的，这样才能做出正确的选择。矩阵是对称的还是厄米特的？是上对角线形式吗？它是半正定的吗。。。在

在运行时有雌雄同体的差异。但根据经验，EIG比SVD便宜。但这也取决于收敛速度，而收敛速度又很大程度上取决于矩阵的条件数，换句话说，矩阵有多不适定。。。在

奇异值分解通常是非常稳健和缓慢的算法，通常用于反演、通过截断进行速度优化、主成分分析，你所处理的矩阵只是一堆垃圾行；）