Критерии сходимости для scipy.eigvalsh

Question

Я использую python (scipy) для вычисления собственных значений симметричной вещественной матрицы. Я в настоящее время с помощью функции

scipy.linalg.eigvalsh 

для вычисления собственных значений (http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigvalsh.html#scipy.linalg.eigvalsh). Глядя на исходный код для eigvalsh, кажется, что python делает вызов пакета fortran. Он также упоминает в документации, что ошибка будет выбрана при вычислении, не сходится.

Мой вопрос: каковы критерии конвергенции? и могу ли я изменить его (относительно легко)?

В моем конкретном приложении я вычисляю собственные значения последовательности матриц и замечаю сильную корреляцию между несколькими собственными значениями. Я хочу знать, не коррелирует ли корреляция исключительно по численным причинам. Если я смогу укрепить критерии конвергенции, тогда я смогу увидеть, увеличивается ли зависимость.

Answer 1

Если я прочитал исходный код справа, используется функция LAPACK dsyevr(). Если я это правильно понимаю, возиться со своими параметрами не обязательно приведет к более высокой точности. Если вам нужна высокая точность, вы можете попробовать mpmath:

import numpy as np 
from mpmath import mp 

print("*** Scipy calculations: ***") 
# Generate matrix: 
n = 25 
AA = np.random.randn(n, n) 
HH = np.dot(AA, AA.T) 

# Calculate eigenvalues and -vectors: 
w, VV = eigh(HH) # eigvalsh() calls also eigh() 

# Check Result: 
HH2 = np.dot(VV, np.diag(w).dot(VV.T)) 
dHH = HH - HH2 
elem_diff_max = np.abs(HH-HH2).max() 
print("Elements differ by maximally {}".format(np.abs(dHH).max())) 
print("Froebenius norm: {}".format(np.linalg.norm(HH-HH2,'fro'))) 

print("") 
print("*** Mpmath calculations (very slow): *** ") 
mp.dps = 40 # number of precision digits for mpmath 
mHH = mp.matrix(HH) # take previous atrix 
mw, mVV = mp.eigh(mHH) # and do eigem decomposition 

# Check rsults: 
mHH2 = mVV*mp.diag(mw)*mVV.T 
mdHH = mHH-mHH2 
#Curiously I could not figure out how to determine abs(mdHH).max(), 
hmax = mp.mpf(0) 
for r in mdHH.tolist(): 
    for c in r: 
     mc = c if c >= 0 else -c 
     hmax = mc if mc > hmax else hmax 

print("Elements differ by maximally {}".format(hmax)) 
print("Froebenius norm: {}".format(mp.norm(mdHH))) 
# Sample output (differs because of randn()): 
# 
# *** Scipy calculations: *** 
# Elements differ by maximally 6.48370246381e-14 
# Froebenius norm: 4.90996840307e-13 

# *** Mpmath calculations (very slow): *** 
# Elements differ by maximally 5.510129769479472693603452518229276614775e-39 
# Froebenius norm: 3.772588954060141733111171961647528674136e-38