scipy.linalg.sparse.eigsh возвращает отрицательные и несогласованные собственные значения для положительной полуопределенной матрицы

Question

Я пытаюсь использовать scipy.linalg.sparse.eigsh (назовем это методом 1: M1) для вычисления наименьших собственных значений матрицы Лапласа вещественной симметричной полуопределенной матрицы W.

В качестве теста я выполнил вычисление для scipy.linalg.eigh (метод 2: M2). Кажется, M1 возвращает разные собственные значения от M2, и, кроме того, эти собственные значения кажутся неправильными.

Вот код:

# Initialization
N=10
np.random.seed(0) # for reproducibility 
Nvp=4

# Create a real symmetric semi-definite matrix W, precision float64
X=np.random.random((N,N))
W=np.dot(X, X.T) 
W=np.array(W, dtype=np.float64) 

# Compute its Laplacian matrix A, precision float64
d=np.array([sum(W[:][j]) for j in range(N)], dtype=np.float64)
D=np.diag(d)
A=D-W
for i in range(N):
    for j in range(N):
        A[i][j]/=np.sqrt(d[i], dtype=np.float64)*np.sqrt(d[j], dtype=np.float64)
A=np.array(A, dtype=np.float64)

Давайте проверим, правильно ли отформатирован A:

>>> A.dtype
dtype('float64')
>>> np.allclose(A, A.T)
True

Теперь давайте запустим несколько тестов:

## 1) Compute A's smallest eigenvalues by 2 different means
wA2, vA2 = la.eigh(A)
wA1, vA1 = sparse.eigsh(A, k=Nvp, sigma=0, which='LM')

## 2) Compute W's smallest eigenvalues by 2 different means
wW2, vW2 = la.eigh(W)
wW1, vW1 = sparse.eigsh(W, k=Nvp, sigma=0, which='LM')

# Output computed eigenvalues
print(wA2[:Nvp])
print(wA1[:Nvp])
print(wW2[:Nvp])
print(wW1[:Nvp])

Вот вывод:

>>>[-1.88737914e-15  9.03999970e-01  9.23513143e-01  9.52678423e-01]

[-4.93242313e-01 -8.14381749e-17  9.22235466e-01  9.44848237e-01]

[0.00575077 0.04770667 0.08565863 0.16319798]

[0.00575077 0.04770667 0.08565863 0.16319798]

Этот первый вывод показывает, что M1 вычисляет отрицательное собственное значение для A, что математически невозможно. Более того, если один проверяет вычисления других, скажем, третий:

>>> np.dot(A, vA2[:,2])-wA2[2]*vA2[:,2]
array([-0.01183104, -0.25059123,  0.47372558, -0.31358618, -0.2968036 ,
   -0.04832199,  0.40973325, -0.01369472,  0.33267859, -0.27122678])

Это даже не близко к нулю. Я должен добавить, что вычисленные собственные значения отличаются каждый раз. Для меня это связано с вектором инициализации. Я бы сказал, что scipy.linalg.sparse.eigsh недостаточно повторяется, чтобы приблизиться к реальному результату, но установка maxiter=1000000 не влияет на странные результаты. Что касается отрицательного собственного значения, я, к сожалению, не в курсе.

Я бегу:

Python 3.7.0 (по умолчанию, 28 июня 2018 г., 13:15:42)

[GCC 7.2.0] :: Anaconda, Inc. на linux

NumPy и SciPy созданы для Intel MKL.

Может ли кто-нибудь просветить меня? Заранее спасибо за ваше время.

Answer 1

Ваша матрица A является единственной, то есть 0 является собственным значением A. Значение sigma в методе инверсии сдвига само по себе не должно быть собственным значением. Взгляните на раздел «Режим сдвига и инверсии спектрального преобразования» документации ARPACK.

Те же формулы приведены в учебнике SciPy . Формула для метода сдвига-инверсии:

inv(A - sigma*M) @ M @ x = nu*x

( Argh , мне бы очень хотелось, чтобы в stackoverflow была разметка LaTeX!) Когда sigma равен 0, A - sigma*M это просто A, а inv(A) не существует.