Ошибочные результаты для собственных разложений numpy/scipy

Question

Я считаю, что scipy.linalg.eig иногда дает противоречивые результаты. Но не каждый раз.

>>> import numpy as np
>>> import scipy.linalg as lin
>>> modmat=np.random.random((150,150))
>>> modmat=modmat+modmat.T  # the data i am interested in is described by real symmetric matrices
>>> d,v=lin.eig(modmat)
>>> dx=d.copy()
>>> vx=v.copy()
>>> d,v=lin.eig(modmat)
>>> np.all(d==dx)
False
>>> np.all(v==vx)
False
>>> e,w=lin.eigh(modmat)
>>> ex=e.copy()
>>> wx=w.copy()
>>> e,w=lin.eigh(modmat)
>>> np.all(e==ex)
True
>>> e,w=lin.eigh(modmat)
>>> np.all(e==ex)
False

Хотя я не лучший мастер по линейной алгебре, я понимаю, что собственное разложение по своей природе подвержено странным ошибкам округления, но я не понимаю, почему повторение вычислений приведет к другому значению. Но мои результаты и воспроизводимость меняются.

Какова именно природа проблемы - ну, иногда результаты приемлемо отличаются, а иногда нет. Вот некоторые примеры:

>>> d[1]
(9.8986888573772465+0j)
>>> dx[1]
(9.8986888573772092+0j)

Разница выше ~3e-13 не кажется невероятно большой. Вместо этого, реальная проблема (по крайней мере, для моего нынешнего проекта) состоит в том, что некоторые из собственных значений не могут согласиться с правильным знаком.

>>> np.all(np.sign(d)==np.sign(dx))
False
>>> np.nonzero(np.sign(d)!=np.sign(dx))
(array([ 38,  39,  40,  41,  42,  45,  46,  47,  79,  80,  81,  82,  83,
    84, 109, 112]),)
>>> d[38]
(-6.4011617320002525+0j)
>>> dx[38]
(6.1888785138080209+0j)

Подобный код в MATLAB, похоже, не имеет этой проблемы.

Answer 1

Разложения по собственным значениям удовлетворяют A V = V лямбда, что является всем, что гарантировано - например, порядок собственных значений - нет.

Ответьте на вторую часть вашего вопроса:

Современные компиляторы / библиотеки линейной алгебры производят / содержат код, который выполняет разные функции в зависимости от того, выровнены ли данные в памяти по (например) 16-байтовым границам. Это влияет на ошибку округления в вычислениях, поскольку операции с плавающей запятой выполняются в другом порядке. Небольшие изменения в ошибке округления могут затем повлиять на такие вещи, как упорядочение собственных значений, если алгоритм (здесь LAPACK/xGEEV) не гарантирует числовую стабильность в этом отношении.

(Если ваш код чувствителен к таким вещам, он некорректен! Запуск, например, на другой платформе или в другой версии библиотеки, может привести к аналогичной проблеме.)

Результаты обычно являются квази-детерминированными - например, вы получаете один из 2 возможных результатов, в зависимости от того, выровнен ли массив в памяти или нет. Если вас интересует выравнивание, проверьте A.__array_interface__['data'][0] % 16,

См. http://www.nccs.nasa.gov/images/FloatingPoint_consistency.pdf для получения дополнительной информации.

Answer 2

Я думаю, что ваша проблема в том, что вы ожидаете, что собственные значения будут возвращены в определенном порядке, и они не всегда будут одинаковыми. Сортируйте их, и вы будете в пути. Если я запускаю ваш код для генерации d а также dx с eig Я получаю следующее:

>>> np.max(d - dx)
(19.275224236664116+0j)

Но...

>>> d_i = np.argsort(d)
>>> dx_i = np.argsort(dx)
>>> np.max(d[d_i] - dx[dx_i])
(1.1368683772161603e-13+0j)