Неявные ошибки алгоритма Фрэнсиса

Question

Я получаю очень противоречивые результаты от реализации имплицитного алгоритма Фрэнсиса для собственных пар для матрицы nxn произвольного размера. Преобладающая проблема заключается в том, что, казалось бы, наугад, мой алгоритм возвращает все «нан» числа для собственных векторов. Другой заключается в том, что решения не соответствуют классическим критериям «Ax = lambda * x».

import numpy as np
import scipy.linalg as sla
import matplotlib.pyplot as plt
import numpy.linalg as la

def QR(A):
    m,n = A.shape
    Qmaster = np.identity(m)
    for j in range(n-1):
        for i in range(j+1,m):
            Q = np.identity(m)
            x_1 = A[j,j] #c
            x_2 = A[i,j] #s
            beta = np.maximum(x_1,x_2) #normalize refelctors by the 
     #maximum of the two to avoid over/underflow\
            x_1 /= beta
            x_2 /= beta
            delta = np.sqrt(x_1**2 + x_2**2)
            x_1 /= delta
            x_2 /= delta
            Q[i,i] = x_1
            Q[j,j] = x_1
            Q[j,i] = (-1)*x_2
            Q[i,j] = x_2
            Qmaster = Qmaster@Q #continuously updates the Q for 
                                 #return
            A = Q.T@A
    return Qmaster, A

def reflector(x,matdim=None):
    # for a given vector x, return the matrix that returns a vector 
    #[tau, 0,0,...]
# When using to reduce a matrix to hessenberg, a tuple argument can 
#be given to pad the matrix with zeros
    x = x.astype(np.float64)
    tau = np.linalg.norm(x)*np.sign(x[0])
    x[0] += tau 
    gamma = x[0]/tau
    for i in range(1,len(x)):
        x[i] /= x[0]
    x[0] = 1
    x = x[np.newaxis,:].T
    Q = np.eye(len(x)) - (gamma)*(x@x.T)
    if matdim != None:
        offset = np.abs(Q.shape[1] - matdim)
        I = np.eye(matdim)
        I[offset:,offset:] = Q
        Q = I

    return Q

def hessenberg(A):
    # converts a matrix into hessenberg form
    #m,n = A.shape
    for i in range(n-2):
        vec = A[i+1:,i]
        Q = reflector(vec,matdim=n)
        A = Q@A
        A = A@Q.T
    return A

def rayleighiteration(A,rho):
    m,n = A.shape
    I = np.eye(n)
    residual = []
    x = np.ones(n)
    iters = 500
    for _ in range(iters):
        x = la.inv(A - rho*I)@x
        x += 1e-13
        x /=la.norm(x)
        rho = x.T@A@x
     return x

def franciseig(A,iters=1000):
    B = hessenberg(A)
    _,n = B.shape
    eigs = np.zeros(n)
    eigvecs = np.zeros((n,n))
    for i in range(iters):
        Q,R = QR(B)
        B = R@Q
    for j in range(n):
        eigs[j] = B[j,j]
        x = rayleighiteration(A,eigs[j])
        eigvecs[:,j] = np.squeeze(x)
    return eigs,eigvecs

np.random.seed(13)
n = 7
X = np.random.rand(n,n)
w,v = sla.eig(X)
eig,vec = franciseig(X)

for i in range(n):
    print(la.norm(X@vec[:,i] - eig[i]*vec[:,i]))

Я ожидаю, что в последней строке кода все напечатанные значения должны быть численно 0 (около 1e-15).