How to use eigenvalue information to efficiently diagonalize matrices?

Question

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 np
from scipy.special import eval_chebyu
from scipy import sparse as sp

N = 10**3

def 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 A

def 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)

Answer 1

Gilbert Strang explains the matrix square root of the second difference matrix here.

In particular,

import numpy
import scipy.linalg

def f(N, p):
    a = 2 * (-1)**(p-1)
    b = 1/numpy.tan((numpy.pi * (2*p - 1)) / (4 + 4*N))
    c = 1/numpy.tan((numpy.pi * (2*p + 1)) / (4 + 4*N))
    return a - b + c

N = 10000
K_tri = numpy.eye(N) - numpy.eye(N, k=1)
K = K_tri + K_tri.T

s = numpy.arange(N)
T = scipy.linalg.toeplitz(f(N, s))
H = scipy.linalg.hankel(f(N, s+2), f(N, N+1+s))
K_sqrt = (0.5 / (N+1)) * (T - H)

print(numpy.max(numpy.abs(numpy.dot(K_sqrt, K_sqrt) - K)))

2.14717132963e-13