scipy eigh дает отрицательные собственные значения для положительно полуопределенной матрицы

У меня возникают проблемы с функцией eigh scipy, возвращающей отрицательные собственные значения для положительных полуопределенных матриц. Ниже представлен MWE.

Функция hess_R возвращает положительную полуопределенную матрицу (это сумма матрицы первого ранга и диагональной матрицы, обе с неотрицательными элементами).

import numpy as np
from scipy import linalg as LA

def hess_R(x):
    d = len(x)
    H = np.ones(d*d).reshape(d,d) / (1 - np.sum(x))**2
    H = H + np.diag(1 / (x**2))
    return H.astype(np.float64)

x = np.array([  9.98510710e-02 ,  9.00148922e-01 ,  4.41547488e-10])
H = hess_R(x)
w,v = LA.eigh(H)
print w

Напечатаны собственные значения

[ -6.74055241e-271   4.62855397e+016   5.15260753e+018]

Если я заменю np.float64 на np.float32 в операторе возврата hess_R, я получу

[ -5.42905303e+10   4.62854925e+16   5.15260506e+18]

вместо этого, поэтому я предполагаю, что это какая-то проблема с точностью.

Есть ли способ исправить это? Технически мне не нужно использовать восемь, но я думаю, что это основная проблема с другими моими ошибками (извлечение квадратных корней из этих матриц, получение NaN и т. д.).


python numpy scipy precision eigenvalue
person angryavian    schedule 24.04.2016    source источник
comment
Если я использую LA.eig вместо LA.eigh, я получаю другие собственные значения: [ 5.15260753e+18+0.j 3.22785571e+01+0.j 4.62855397e+16+0.j]   -  person Peaceful    schedule 24.04.2016
comment
ИМХО, ваша функция Hess_R не возвращает фактическую матрицу Гессе. поэтому eigh возвращает ложный результат в вашем случае.   -  person B. M.    schedule 24.04.2016
comment
@Б.М. Не могли бы вы еще пояснить, что вы имеете в виду? Что вместо этого возвращает функция?   -  person angryavian    schedule 24.04.2016
comment
Вы можете увидеть разные результаты в зависимости от того, с какой реализацией LAPACK связана ваша версия numpy. Используя OpenBLAS v2.18, я получаю собственные значения 3.14000000e+02, 4.62855397e+16, 5.15260753e+18 с scipy.linalg.eigh и 1.06862038e+02, 4.62855397e+16, 5.15260753e+18 с numpy.linalg.eigh для вашего примера H (разница также может быть частично связана с потерей точности при выводе значений с плавающей запятой в x в виде строк).   -  person ali_m    schedule 24.04.2016
comment
К вашему сведению: вы можете заменить np.ones(d*d).reshape(d,d) на np.ones((d, d)).   -  person Warren Weckesser    schedule 24.04.2016

Ответы (1)


arrow_upward
4
arrow_downward

Я думаю, проблема в том, что вы достигли предела точности с плавающей запятой. Хорошее эмпирическое правило для результатов линейной алгебры состоит в том, что они хороши примерно до одной части 10 ^ 8 для float32 и примерно до одной части 10 ^ 16 для float 64. Похоже, что отношение вашего наименьшего к наибольшему собственное значение здесь меньше 10^-16. Из-за этого возвращаемому значению нельзя доверять, и оно будет зависеть от деталей используемой вами реализации собственного значения.

Например, вот четыре разных решателя, которые должны быть у вас в наличии; посмотрите на их результаты:

# using the 64-bit version
for impl in [np.linalg.eig, np.linalg.eigh, LA.eig, LA.eigh]:
    w, v = impl(H)
    print(np.sort(w))
    reconstructed = np.dot(v * w, v.conj().T)
    print("Allclose:", np.allclose(reconstructed, H), '\n')

Выход:

[ -3.01441754e+02   4.62855397e+16   5.15260753e+18]
Allclose: True 

[  3.66099625e+02   4.62855397e+16   5.15260753e+18]
Allclose: True 

[ -3.01441754e+02+0.j   4.62855397e+16+0.j   5.15260753e+18+0.j]
Allclose: True 

[  3.83999999e+02   4.62855397e+16   5.15260753e+18]
Allclose: True

Обратите внимание, что все они согласны с двумя большими собственными значениями, но значение наименьшего собственного значения меняется от реализации к реализации. Тем не менее, во всех четырех случаях входная матрица может быть восстановлена ​​с точностью до 64 бит: это означает, что алгоритмы работают так, как ожидалось, с доступной им точностью.

person jakevdp    schedule 24.04.2016