Диагонализация плохо обусловленной матрицы и невозможность вычисления собственных векторов. Различные результаты с NumPy / Scipy

Question

Я имею дело с разреженной матрицей с очень маленькими элементами. Рассмотрим вектор:

vec=[-1.e-76 -1.e-72 -1.e-68 -1.e-64 -1.e-60 -1.e-56 -1.e-52 -1.e-48 -1.e-44
-1.e-40 -1.e-36 -1.e-32 -1.e-28 -1.e-24 -1.e-20 -1.e-16 -1.e-12 -1.e-08
-1.e-04 -1.e-02 -1.e-04 -1.e-08 -1.e-12 -1.e-16 -1.e-20 -1.e-24 -1.e-28
-1.e-32 -1.e-36 -1.e-40 -1.e-44 -1.e-48 -1.e-52 -1.e-56 -1.e-60 -1.e-64
-1.e-68 -1.e-72 -1.e-76]

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

H=sps.diags([vec,vec],[-1,1],dtype='f8')

Меня интересуют собственные значения, но еще больше - собственные векторы. , Насколько я знаю, есть два способа справиться с диагонализацией: scipy.linalg и numpy.linalg и первый лучше.

 denseHam=H.toarray()

Правильный спектр собственных значений задается всеми этими функциями:

import numpy as np
import scipy.linalg as la
s1= la.eigvalsh(denseHam)
s2= np.linalg.eigvalsh(denseHam)
s3= np.linalg.eigvals(denseHam) #I did not expect that!

Правильный спектр:

spectrum=[-3.16230928e-03 -3.16227766e-08 -3.16227766e-13 -3.16227766e-18
-3.16227766e-23 -3.16227766e-28 -3.16227766e-33 -3.16227766e-38
-3.16227766e-43 -3.16227766e-48 -3.16227766e-53 -3.16227766e-58
-3.16224604e-63  3.16224604e-63  3.16227766e-58  3.16227766e-53
 3.16227766e-48  3.16227766e-43  3.16227766e-38  3.16227766e-33
 3.16227766e-28  3.16227766e-23  3.16227766e-18  3.16227766e-13
 3.16227766e-08  3.16230928e-03]

Тем не менее, другие функции (которые также включают вычисление собственных векторов) не работают, и я не могу продолжать, потому что мне нужны собственные векторы.

Я должен сказать, что C ++ способен правильно вычислять и собственные векторы.

Итак, у меня два вопроса:

1. Почему функция np.linalg.eigh(denseHam) дает спектр, отличный от np.linalg.eigvalsh(denseHam)?
2. Есть ли способ правильно вычислить собственные векторы с помощью python?

Заранее большое спасибо!

--- ОБНОВЛЕНИЕ ------ Я привожу здесь минимальный полный пример. Обратите внимание на явное вырождение numpy.linalg.eigh:

import numpy as np
import scipy.sparse as sps

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

s1=np.linalg.eigvalsh(denseHam)
(s2,basis)=np.linalg.eigh(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem])

Answer 1

Краткий ответ: Попробуйте LAPACK's scipy.linalg.lapack.dsyev()!

# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

Функции np.linalg.eigvalsh() и np.linalg.eigh() оба вызывают LAPCK DSYEVD, как указано в их документации. Тем не менее, они предоставляют разные собственные значения: eigh() не удается. Причина, вероятно, выгравирована в исходном коде DSYEVD. Действительно, если собственные векторы не нужны, процедура масштабирует матрицу, сводит матрицу к тригональной форме (DSYTRD) и наконец звонит DSTERF. Эта последняя процедура вычисляет все собственные значения симметричной трехдиагональной матрицы, используя вариант Пал-Уокера-Кахана алгоритма QL или QR. Напротив, если необходимы собственные векторы, DSYEVD переключается на DSTEDC после масштабирования и приведения к треугольной форме. DSTEDC вычисляет все собственные значения и собственные векторы симметричной трехдиагональной матрицы с использованием метода «разделяй и властвуй».

1. Небольшая норма входной матрицы не имеет значения, поскольку в таком случае, вероятно, применяется масштабирование. Поскольку действительная симметричная матрица имеет собственные значения очень различной величины (от 3.16224604e-63 до 3.16230928e-03), плохо подготовлен. Точность большинства процедур линейной алгебры, включая вычисление собственных значений, существенно зависит от этого свойства. Предоставленная матрица уже имеет тригональную форму.

2. np.linalg.eigh() не удается. Вероятно, это связано с использованием стратегии «разделяй и властвуй».

3. np.linalg.eigvalsh(), кажется, работает. Поэтому похоже, что DSTERF сработало. Тем не менее, эта процедура предоставляет только собственные значения.

Поскольку LAPACK предоставляет различные алгоритмы для вычисления собственных значений и собственных векторов , ваш код C ++, вероятно, использует другой, такой как dsyev(). После масштабирования и сводя матрицу к тригональной форме, эта процедура сначала вызывает DORGTR для генерации ортогональной матрицы, а затем вызывает DSTEQR для получения собственных векторов. Надеемся, что эта функция может быть вызвана с помощью Низкоуровневых функций LAPACK от scipy (scipy.linalg.lapack)

Я добавил несколько строк в ваш код для вызова этой функции. scipy.linalg.lapack.dsyev() вычисляет собственные значения, подобные значениям np.linalg.eigvalsh() для этой плохо обусловленной матрицы.

import numpy as np
import scipy.sparse as sps
import scipy.linalg.lapack

vec=np.array([-1.e-76, -1.e-72, -1.e-68, -1.e-64, -1.e-60, -1.e-56, -1.e-52,
       -1.e-48, -1.e-44, -1.e-40, -1.e-36, -1.e-32, -1.e-28, -1.e-24,
       -1.e-20, -1.e-16, -1.e-12, -1.e-08, -1.e-04, -1.e-02, -1.e-04,
       -1.e-08, -1.e-12, -1.e-16, -1.e-20, -1.e-24, -1.e-28, -1.e-32,
       -1.e-36, -1.e-40, -1.e-44, -1.e-48, -1.e-52, -1.e-56, -1.e-60,
       -1.e-64, -1.e-68, -1.e-72, -1.e-76])
H=sps.diags([vec,vec],[-1,1],dtype='f8')
denseHam=H.toarray()

# eigenvalues only, boils down to LAPACK's DSTERF
s1=np.linalg.eigvalsh(denseHam)
# LAPACK's dsyevd, divide and conquer
(s2,basis)=np.linalg.eigh(denseHam)
# LAPACK's dsyev instead of LAPACK's dsyevd
(s3,basis3,error)=scipy.linalg.lapack.dsyev(denseHam)

print("Note the difference between the eigenvalues computed with eigvalsh (1stcolumn) and eigh (2nd column)")
for elem in range(len(s1)):
    print (s1[elem],"     ",s2[elem], s3[elem])

Вы также можете избежать сокращения до тригональной формы, поскольку ваша матрица уже является тригональной. Масштабирование, вероятно, требуется, чтобы избежать численных проблем. Тогда DORGTR и DSTEQR могут быть вызваны напрямую через функции LAPACK для Cython . Но это требует больше работы и этапа компиляции.

Answer 2

Несколько удивительно, что они дают разные результаты, так как я ожидаю, что и numpy, и scipy позвонят в LAPACK, но в целом это был и мой опыт.

Обратите внимание, что привязки scipy предлагают больше аргументов для игры; и numpy, вероятно, использует разные значения по умолчанию. Некоторые эксперименты кажутся необходимыми; Ваша проблема заключается не только в том, чтобы иметь очень маленькие элементы (которые могут быть решены простым масштабированием, если это приводит к недостаточному заполнению), но ваша проблема также очень «жесткая», с собственными значениями, охватывающими более 70 порядков. C ++ может дать вам собственные векторы, но я не удивлюсь, если они будут загрязнены числовым шумом до такой степени, что станут бесполезными.

Звучит как проблема, где было бы гораздо надежнее решить ее в каком-то преобразованном / предварительно подготовленном пространстве. Строка документации не говорит, могут ли функции LAPACK обрабатывать 128-битное число с плавающей запятой; в прошлый раз, когда я пытался, они этого не делали, но если они это сделают сейчас, обязательно используйте это как минимум вместо 64-битных.