Действительно, numpy.linalg.svd и numpy.linalg.eigh не вызывают одну и ту же процедуру Лапака.С одной стороны, numpy.linalg.eigh относится к LAPACK dsyevd(), в то время как numpy.linalg.svd использует LAPACK dgesdd().
Общей точкой между этими подпрограммами является использование алгоритма «разделяй и властвуй» Куппена, впервые разработанного для решения трехдиагональных задач на собственные значения.Например, dsyevd() обрабатывает только эрмитову матрицу и выполняет следующие шаги, только если требуются собственные векторы:
Преобразование матрицы в трехдиагональную форму с помощью DSYTRD ()
Вычислить собственные векторы трехдиагональной матрицы с использованием алгоритма «разделяй и властвуй» с помощью DSTEDC ()
Применить отражение Домохозяина, о котором сообщает DSYTRD ()используя DORMTR ().
Напротив, для вычисления SVD dgesdd() выполняет следующие шаги, в случае задания == A (требуются U и VT):
- Бидиагонализация A с использованием
dgebrd() - Вычисление SVD бидиагональной матрицы с использованием алгоритма «разделяй и властвуй» с использованием
DBDSDC() - Отмена бидиагонализации с использованием матриц Pи Q возвращается
dgebrd(), применяя dormbr() дважды, один раз для U и один раз для V.
В то время как фактические операции, предусмотренные LAPACK, очень разные, стратегии во всем мире схожи.Это может быть обусловлено тем фактом, что вычисление SVD общей матрицы A аналогично выполнению собственного разложения симметричной матрицы A ^ TA
Относительно точности и характеристик разрыва и захвата SVD Лапака см. Thisобзор методов SVD :
- Они часто достигают точности SVD на основе QR, хотя это не доказано.
- Худший случай - 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()
Для таких одномерных диффузионных матриц восемь опережают SVD, но фактическая сложность аналогична, чуть ниже, чем n ^ 3, что-то вроде n ^ 2.5.,
Также можно выполнить проверку точности.