редкие eigs Сципи дают неправильные собственные векторы

Question

Я немного сбит с толку в отношении следующей проблемы:

Я вычисляю фиксированные точки квантового канала, что означает, что я хочу вычислить ведущий собственный вектор конкретной матрицы.Матрица такова, что ее размерность равна n ^ 2 xn ^ 2 и определена таким образом, что ведущее собственное значение, преобразованное в матрицу с формой nxn, является положительной матрицей (самосопряженной с положительными собственными значениями).

ЕслиЯ делаю это с scipy.sparse.linalg.eigs, но я получаю неправильные результаты.Однако точное вычисление (с использованием scipy.linalg.eig) работает нормально.Я пытался поиграть с аргументами для k и ncv для решателя, но я не смог правильно работать, если не установил k=n**2, в этом случае eigs просто ссылается на eig.Это, однако, не сработает в том случае, если я на самом деле имею в виду, где канал (super_op в приведенном ниже скрипте) фактически кодируется как LinearOperator.Поэтому я полагаюсь на использование eigs: /

Кто-нибудь знает, как правильно запустить этот прогон?

Спасибо всем заранее!

import numpy as np
from numpy.random import rand
from numpy import tensordot as td
from scipy.sparse.linalg import eigs
from scipy.linalg import eig


n = 16
d = 3

kraus_op = .5 - rand(n, d, n) + 1j * (.5 - rand(n, d, n))
super_op = td(kraus_op, kraus_op.conj(), [[1], [1]]).transpose(0, 2, 1, 3)

########
# Sparse
########

vals, vecs = eigs(super_op.reshape(n**2, n**2), k=n*(n-1), which='LM')

rho = vecs[:,0].reshape(n, n)

print('is self adjoint: ', np.allclose(rho, rho.conj().T))

super_op_times_rho = td(super_op, rho, [[2, 3], [0, 1]])

print('super_op(rho) == lambda * rho :', np.allclose(rho, super_op_times_rho/vals[0]))

########
# Exact
########

vals, vecs = eig(super_op.reshape(n**2, n**2))

rho = vecs[:,0].reshape(n, n)

print('is self adjoint: ', np.allclose(rho, rho.conj().T))

super_op_times_rho = td(super_op, rho, [[2, 3], [0, 1]])

print('super_op(rho) == lambda * rho :', np.allclose(rho, super_op_times_rho/vals[0]))

результат:

is self adjoint:  False
super_op(rho) == lambda * rho : True
is self adjoint:  True
super_op(rho) == lambda * rho : True

Для полноты:

Python 3.5.2 numpy 1.16.1
scipy 1.2.1

Answer 1

В конце концов, я нашел решение с помощью моих коллег:

В то время как eig дает собственные векторы (1), отсортированные (по величине) и (2), так, что первая запись действительна, eigs имеет другую сортировку, которая не должна быть такой же, как в eig, а также не регуляризирует комплексную фазу собственного вектора.Исправление фазы легко сделать, разделив тензор на фазу первого элемента (что соответствует тому, что первый диагональный элемент является действительным, чтобы избавиться от свободы выбора сложной фазы собственного вектора и сделать возможной эрмитность ...),

Таким образом, исправленный фрагмент кода для разреженного случая:

vals, vecs = eigs(super_op.reshape(chi**2, chi**2), k=chi*(chi-1), which='LM')
# find the index corresponding to the largest eigenvalue
arg = np.argmax(np.abs(vals))
rho = vecs[:,arg].reshape(chi, chi)
# regularize the output array 
rho *= np.abs(rho[0, 0])/rho[0, 0]