Элементы матрицы обратной для плохо обусловленной матрицы

Question

Я пытаюсь найти элементы матрицы, обратной для плохо обусловленной матрицы

Рассмотрим комплексную неэрмитову матрицу M, я знаю, что эта матрица имеет одно нулевое собственное значение и поэтому является единственной,Однако мне нужно найти сумму матричных элементов: v@f(M)@u, где u и v - оба вектора, а f (x) = 1 / x (фактически обратная матрица). Я знаю, что нулевое собственное значение не влияет на эту сумму, поэтому с сингулярностью нет явной проблемы. Тем не менее, мой код очень численно нестабилен, и я полагаю, что это является следствием ошибки при поиске собственных значений системы.

Начнем с построения предварительных матриц:

import numpy as np
import scipy as sc

g0 = np.array([0,0,1])
g1 = np.array([0,1,0])
e0 = np.array([1,0,0])

sm = np.outer(g0, e0)
sp = np.outer(e0, g0)



def spre(op):
    return np.kron(np.eye(op.shape[0]),op)

def spost(op):
    return np.kron(op.T,np.eye(op.shape[0]))

def sprepost(op1,op2):
    return np.kron(op1.T,op2)

sm_reg = spre(sm)
sp_reg = spre(sp)

spsm_reg=spre(sp@sm)

hil_dim = int(g0.shape[0])
cav_proj= np.eye(hil_dim).reshape(hil_dim**2,)

rho0 =(np.outer(e0,e0)).reshape(hil_dim**2,)




def ham(g):
    return g * (np.outer(g1,e0) + np.outer(e0, g1))


def lind_op(A):
    L = 2 * sprepost(A,A.conj().T) - spre(A.conj().T @ A)
    L += - spost(A.conj().T @ A)
    return L



def JC_lio(g, kappa, gamma):
    unit = -1j * (spre(ham(g)) - spost(ham(g)))
    lind = gamma * lind_op(np.outer(g0 , e0)) + kappa * lind_op(np.outer(g0 , g1))
    return unit + lind

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

def power_int(g, kappa, gamma):

    # Construct the non-Hermitian matrix of interest
    lio = JC_lio(g,kappa,gamma)

    #Find its left and right eigenvectors:
    ev, left, right = scipy.linalg.eig(lio, left=True,right=True)

    # Find the appropriate normalisation factors
    norm = np.array([(left.conj().T[ii]).dot(right.conj().T[ii]) for ii in range(len(ev))])

    #Find the similarity transformation for the problem
    P = right
    Pinv = (left/norm).conj().T


    #find the projectors for the Eigenbasis
    Proj = [np.outer(P.conj().T[ii],Pinv[ii]) for ii in range(len(ev))]


    #Find the relevant matrix elements between the Eigenbasis and the projectors --- this is where the zero eigenvector gets removed
    PowList = [(spsm_reg@ Proj[ii] @ rho0).dot(cav_proj) for ii in range(len(ev))]

    #apply the function
    Pow = 0

    for ii in range(len(ev)):
        if PowList[ii]==0:
            Pow = Pow
        else:
            Pow += PowList[ii]/ev[ii]


    return -np.pi * np.real(Pow)


#example run:

grange = np.linspace(0.001,10,40)

dat = np.array([power_int(g, 1, 1) for g in grange])

Запуск этого кода приводит к чрезвычайно осциллирующим результатам, где я ожидаю плавную кривую. Я подозреваю, что эта ошибка связана с плохой точностью определения собственных векторов, но я не могу найти какую-либо документацию по этому вопросу. Любые идеи будут приветствоваться.