Polynomial filter on symmetric tridiagonal matrices

During my research, I encountered a problem involving symmetric tridiagonal matrices and their eigenvectors, which I have been attempting to solve since then. So far, I have only managed to obtain a solution in a special case. The problem concerns two real-valued symmetric tridiagonal matrices $A$ and $M$ of size $(l+1) \times (l+1)$.

Definitions:

Matrix $A$ is hollow (zeros on its main diagonal) and is defined as follows:

\begin{equation} A := A^{m, l} := \begin{bmatrix} 0 & a_1 & 0 & \ldots & 0 & 0 \\ a_1 & 0 & a_2 & \ldots & 0 & 0 \\ 0 & a_2 & 0 & \ldots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & 0 & a_l \\ 0 & 0 & 0 & \ldots & a_l & 0 \end{bmatrix} \label{eq:A_matrix} \end{equation}

with $[a_1, \dots, a_l]$ in the off-diagonal given as $a_k=\sqrt{k(m-k+1)}$. The parameter $m$ is an integer and has to be larger or equal $l$, i.e., $l\leq m$.

Matrix $M$ shares the same off-diagonal elements as $A$ but includes additional nonzero diagonal components, \begin{equation} M := A + D, \end{equation} where $D$ is a diagonal matrix. The entries of $D$ are parameterized by a variable $0 < s \leq \tfrac{1}{2}$ and are given by \begin{equation} D = D(s) := \mathrm{diag}(0, d, 2d, \dots, l d), \label{eq:D_matrix} \end{equation} where the shift $d := d(s)$ is defined in terms of $s$ as \begin{equation} d(s) := -\frac{2b(s)}{a(s)}, \qquad a(s) := 2\sqrt{s(1-s)}, \quad b(s) := 2s - 1. \end{equation}

I am interested in a polynomial $P(A)$ defined over matrix $A$. The polynomial is constructed as follows: Let $\lambda_0< \dots < \lambda_l$ denote the eigenvalues of the matrix $M$, and $v_0,\dots ,v_l$ the corresponding eigenvectors. Since matrix $M$ is a real symmetric tridiagonal matrix (a Jacobi matrix), its eigenvalues are guaranteed to be real and nondegenerate. The polynomial $P$ is defined to have roots at the values \begin{equation} \tilde\lambda_i = \lambda_i a+m b \end{equation} for all $i$ except for the largest eigenvalue $\lambda_l$, i.e, $i\neq l$.
The polynomial $P$ thus has $l$ roots and therefore is of degree $l$. This defines $P$ uniquely up to an overall multiplicative constant, which is irrelevant for our purposes.

Furthermore, we define a upper triangular matrix $K \in \mathbb{R}^{(l+1)\times (l+1)}$ via \begin{equation} K_{i,j} := \sqrt{\frac{\binom{m}{i}}{\binom{m}{j}}}\binom{m-i}{j-i}a^ib^{j-i},\quad \text{for } i\leq j \end{equation} and zero elsewhere.

Claim to be proven:

Let $e_1$ denote the first basis state of $\mathbb{R}^{l+1}$, i.e., $e_1 = [1,0,\dots,0]^t$. I want to prove: \begin{equation} \boxed{ K\cdot P(A) e_1 \propto v_{l}. } \end{equation} In words, applying the polynomial $P(A)$ to the initial vector $e_1$ and subsequently acting with a correction matrix $K$ yields the principal eigenvector of $M=A+D$ (up to normalization). This relation holds for all $l,m\in \mathbb{N}$, with $l \leq m$ and $0<s\leq 1/2$.

Comments:

(C1): What I managed to prove so far is the special case of $s=1/2$. There \begin{equation} a=1, b = 0, d=0, K = I \end{equation} with $I$ denoting the identity matrix. This also results into $M = A$ and that $v_l$ is the principal vector of $A$ itself. Furthermore we have $\tilde\lambda_i = \lambda_i$. Hence the polynom $P$ has roots on all eigenvalues of matrix $A$ except $\lambda_l$ Thus, up to a global factor

\begin{equation} P(x) = \prod_{i=0}^{l-1}(x-\lambda_i). \end{equation}

To see that $P(A)$ acts as a filter for the principal eigenvector one decomposes $e_1$ into the eigenstates of matrix $A$ via $e_1 = \sum_{i=0}^lc_1 v_i$. One gets:

\begin{equation} P(A)e_1 = \sum_{i=0}^lc_i P(A)v_i = \sum_{i=0}^lc_i P(\lambda_i)v_i = \sum_{i=0}^lc_i \prod_{j=0}^{l-1}(\lambda_i-\lambda_j)v_i. \end{equation}

It is easy to see that

\begin{equation} \prod_{j=0}^{l-1}(\lambda_i-\lambda_j) \propto \delta_{il}. \end{equation}

Hence, only one component related to index $l$ survives resulting into

$P(A)e_1\propto v_l$.

(C2): If $l=m$ matrix $A$ resembles the angular momentum operators $J_x$ known in quantum mechanics and $D$ can be related to $J_z$, s.t. \begin{equation} M = \frac{2}{a}(aJ_x+bJ_z)+d\frac{m}{2} I \end{equation} is related to a tilted field Hamiltonian. If $l<m$ it is a truncated version of it.

(C3): Matrix $A$ and $D$ do not commute if $d\neq 0$ i.e., $s\neq1/2$.

(C4): Powers of $A$ applied onto the initial vector $e_1$ form a Krylov subspace $\{A^le_1, \dots Ae_1, e_1\}$

(C4): It can be easily checked that $a^2+b^2 = 1$. Thus the parameters $a$ and $b$ appearing in $K$ and $M$ can be related to a rotation angle $\theta$ as $a = \cos \theta, b = \sin \theta$.

(C5): The inverse of $K$ is structurally similar to $K$. It is again a upper triangular matrix with entries \begin{equation} K_{i,j} := \sqrt{\frac{\binom{m}{i}}{\binom{m}{j}}}\binom{m-i}{j-i}a^{-i}(-b)^{j-i},\quad \text{for } i\leq j \end{equation} and zero elsewhere.

(C6): Eigenstates of symmetric tridiagonal matrices are related to orthogonal polynomials.

(C7): I append a short Python script which numerically supports the claim.

import numpy as np
from math import sqrt
from numpy.linalg import matrix_power
from scipy.special import binom
from numpy.polynomial import Polynomial


if __name__ == '__main__':
    # Parameters for matrix (l <= m) m,l integers
    m = 10
    l = 6  # l+1 = matrix size

    # Parameter for shifts (0 < s <= 0.5)
    s = 1/2

    # meta parameters
    a = 2 * sqrt(s*(1-s))
    b = 2 * s - 1

    # shift defining matrix M
    d = -2 * b/a


    # Build tri-diagonal matrices
    off_diag = [np.sqrt(i * (m - i + 1)) for i in range(1, l+1)]
    diag = np.arange(l + 1) * d

    A = np.diag(off_diag, 1) +  np.diag(off_diag, -1)
    M = A + np.diag(diag)

    print('Matrix M =\n', M, '\n')

    # Calculate spectra and principal eigenvector
    evals, evecs = np.linalg.eig(M)
    principal_vector = evecs[:,np.argmax(evals)]

    def shift(x):
        return x * a + m * b

    # Eigenvalues shifted without the largest
    evals_shifted = sorted([shift(e) for e in evals])[:l]


    ''' 
        Polynomial P(A)
    '''

    # Convert from roots to coefficients of polynomial
    P = Polynomial.fromroots(evals_shifted)
    coefs = P.coef

    # Define polynomial
    PolyA = sum(
        [
            coefs[i] * (matrix_power(A,i))
            for i in range(l+1)
        ]
    )

    # Correction matrix
    K = np.zeros((l + 1, l + 1))
    for row in range(l + 1):
        for col in range(row, l + 1):
            K[row, col] = (
                   sqrt(binom(m, row) / binom(m, col)) *
                   binom(m - row, col - row) *
                   (a ** row) *
                   (b) ** (col - row)
            )


    # Initial vector
    e1 = np.array([1]+[0]*l)

    PAe1 = (K  @ PolyA)  @ e1
    PAe1 /= np.linalg.norm(PAe1)

    print(
        'Pricipal eigenvector of M =\n\t',
        principal_vector
    )
    print(
        'P(A)e1 = \n\t',
        PAe1
    )
    
    # note that the normalized vector still can have a global minus sign in front
    error  = min(
        np.linalg.norm(principal_vector - PAe1),
        np.linalg.norm(-principal_vector - PAe1)
    )
    print('error = ', error)

    ''' 
    # Inverse of K
    Kinv = np.zeros((l+1,l+1))
    for row in range(l+1):
        for col in range(row, l+1):
            Kinv[row, col] = (
                    sqrt(binom(m, row) / binom(m, col)) *
                    binom(m-row,col-row) *
                    (a ** (-col)) *
                    (-b) ** (col-row)
            )
    
    '''

```