I apologize if this question in a more general form has been asked before. I have a tridiagonal Toeplitz matrix $K$, whose eigenvalues and eigenvectors are known analytically for any dimension $N$ [1]. Specifically, in the notation of Section 2.2 of this paper, the matrix I wish to diagonalise is such that $a_1 = a_2$ $= -2$, and $ b_1 = b_2=c_1=c_2=1$. Or $$K = \begin{pmatrix}-2 & 1 & 0 & 0&..&0 \\1 & -2 & 1 & 0&..&0 \\0 & 1 & -2 & 1&..&0 \\. \\.\end{pmatrix}$$Then theorem 2.1 gives its eigenvalues and eigenvectors.They are: $\lambda_i = -4 \cos^2\left({\frac{(i+1)\pi}{2(N+1)}}\right)$ and eigenvectors $v_i = (.,.,\left[U(\frac{j}{2}, \frac{(\lambda_i+2)^2-2}{2}) + U(\frac{j}{2}-1, \frac{(\lambda_i+2)^2-2}{2})\right],.,.)^T$, where $U(j,x)$ is the Chebyshev polynomial of the second kind of order $j$. What I further wish to do is to use the matrix of eigenvectors, say $A$, to implement similarity transforms. (Note that $A$ is not sparse at all in this example). How best can I do it computationally? In Python, for $ N = 10^4 $, the program takes a lot of time to even calculate $A$ and store the matrix, presumably because of the presence of the Chebyshev polynomial. On the other hand, using the sparse linear algebra module scipy.sparse.linalg.eigs finishes the task faster, even though I don't give it any information! So, my question is this- how do I use a combination of my knowledge of the true eigenvalues and eigenvectors and the efficient linear algebra routines in SciPy to implement the similarity transform on a diagonal matrix D?Even if I use just the information about eigenvalues, can I efficiently compute the eigenvectors? Are there routines that compute the eigenvectors after taking an input of eigenvalues if they are known?For instance, the following code was my trial in Python:import numpy as npfrom scipy.special import eval_chebyufrom scipy import sparse as spN = 10**3def Mat(x,j):if j%2==0:    if j==0:        return 1    else:        return eval_chebyu(j/2, ((x+2)**2-2)/2) + eval_chebyu(j/2-1, ((x+2)**2-2)/2)else:    return (x+2)*eval_chebyu((j-1)/2, ((x+2)**2-2)/2)def norm(A):A=A.transpose()/np.sqrt(np.einsum('...i,...i', A,A ))return Adef tridiag(a, b, c, N):return sp.diags(a*np.ones(N-1), -1) + sp.diags(b*np.ones(N), 0) +sp.diags(c*np.ones(N-1), 1)K = tridiag(1,-2,1,N)A = np.zeros([N,N])lam = np.zeros(N)for i in range(0,N):    lam[i] = -4*(np.cos((i+1)*np.pi/(2*(N+1))))**2    for j in range(0,N):        A[i,j]= Mat(lam[i],j)A = norm(A)Tq_new = sp.diags(1/(2*omega*np.sqrt(-lam)),0)Tq_old = A.transpose().dot(Tq_new.todense()).dot(A)
How to use eigenvalue information to efficiently diagonalize matrices?
linear algebra;matrices;eigenvalues