# How to use eigenvalue information to efficiently diagonalize matrices?

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)

• Just to be clear, you want to compute something like $U D U^T$ when $U$ is the orthogonal matrix of eigenvectors of your specific $K$ matrix, and $D$ is some arbitrary diagonal matrix? – k20 Jan 17 '14 at 21:52
• Yes, that's right. The diagonal matrix is related to the eigenvalues in the following way. The diagonal entries are $\sqrt{\lambda_i}$ and $\sqrt{\frac{1}{\lambda_i}}$ – Abhinav Jan 17 '14 at 22:30
• Is your question just how to compute $A$ quickly, or is your question about a downstream calculation (like computing some other thing using $A$)? – k20 Jan 17 '14 at 22:54
• By the way, I'm not an expert on signal processing but you might be able to use the SciPy discrete cosine transform functions in some way. – k20 Jan 17 '14 at 23:02
• You could probably use the latter with scipy.solveh_banded. But going through A will require like O(N^3) for the similarity transform on a diagonal matrix D, whereas you could possibly use signal processing tricks to do it more efficiently like O(N^2) or O(N^2 log(N)) or something. I suspect this, because your matrix K is so close to the discrete fft or dct matrix. www-math.mit.edu/~gs/PIX/cupcakematrixtxt.jpg Also if you are only looking for the similarity transformation of the sqrt of eigenvalues this is the matrix square root (sqrtm). – k20 Jan 18 '14 at 0:18

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