Eigenvectors of a symmetric positive definite Toeplitz matrix

Question

I wish to efficiently compute the eigenvectors of an n x n symmetric positive definite Toeplitz matrix K. A full eigendecomposition would be even better.

Although I assumed this would be a well addressed problem in the numerical linear algebra literature, I have found surprisingly little on this topic, despite extensive searching.

It is possible to solve linear systems involving Toeplitz efficiently using a variety of techniques. For example, one can embed a Toeplitz matrix into a circulant matrix to efficiently perform Toeplitz matrix vector products using fast Fourier transforms. Then one can solve linear systems using linear conjugate gradients, which only involves matrix vector products. This procedure would require O(n log n) time and O(n) space, for convergence to within machine precision.

There are also algorithms which will compute the exact inverse and determinant of K in O(n^2) time and O(n^2) space [O(n) space if only the determinant is required]. Moreover, slide 108 of http://www.math.cinvestav.mx/~grudsky/Talks/Talk_11.pdf has some nice pointers about how to approximate the eigenvalues of a Hermitian Toeplitz matrix at the cost of a single FFT.

However, I am struggling to find a good approach for a full eigendecomposition of K.

Any pointers, or thoughts greatly appreciated! It would be embarrassing if I had missed something obvious in the literature, which is possible. At the same time, I would still be happy to know there is a good solution to this problem.

As a minor caveat, I am not troubled by (small) approximations, e.g., algorithms that use fast Fourier transforms (FFTs). Ideally I would also like an algorithm that does not just have good asymptotic complexity, but has for instance a break-even runtime with standard O(n^3) eigendecompositions for small n (e.g., n <= 500).

Thank you.

PS. Possible approaches:

• A first approach could be to approximate the eigenvalues as suggested in the Grudsky slides, and then use the inverse iteration method for the eigenvectors, with circulant embeddings for fast Toeplitz matrix vector products. But this would have a complexity of O(n^2 log n) and is approximate and would be quite slow for all but very large values of n. I feel the eigenvectors can be determined accurately and efficiently in O(n^2).

• Another more promising approach could be to treat a symmetric Toeplitz matrix as a low rank correction to a circulant matrix (for which eigenvectors are readily available). But I am not sure about the details.

Answer 1

First, a meta-comment. Numerical linear algebra is de facto a part of engineering or, more precisely, computer science (practical sciences studying how we do something and dealing with our practical relationships to mathematical objects and not the primary objects themselves) so servers about computer science similar to this one may be more helpful in solving such questions.

Now, my guess. The problem to find eigenstates of a Toeplitz matrix is not exactly solvable in the same sense as in the special case of the circulant which makes me expect that the solutions to this problem won't be substantially faster than the solution of the problem for general matrices.

But if you can find the eigenvalues by the method you mentioned, I believe that the problem to find eigenstates is a special case of a Toeplitz system,

https://en.wikipedia.org/wiki/Toeplitz_matrix#Solving_a_Toeplitz_system

although a singular one, and there exists an $O(n^2)$ algorithm for that, the Levinson algorithm

https://en.wikipedia.org/wiki/Levinson_recursion

The Toeplitz system is the problem to find $x$ obeying $Ax=b$ for a given Toeplitz matrix $A$ and a column $b$. If you know an eigenvalue $\lambda$ and you want to find the eigenvector $x$, you are solving $(A-\lambda\cdot {\bf 1})x=0$ which seems like a Toeplitz system, one with a shifted (still Toeplitz) new matrix $A-\lambda\cdot {\bf 1}$, albeit a system with $b=0$. Of course, the solution is not unique – the scaling is arbitrary – and my experience is non-existent and preventing me from saying whether the method breaks down for $b=0$ or $b\to 0$, which you may try if $b=0$ doesn't work, but it seems as an obvious attempt constructed from the pieces that are out there.

Answer 2

You might try searching for algorithms that find eigenvalues and eigenvectors of bisymmetric matrices.

That is, a symmetric Toeplitz matrix is bisymmetric: it is symmetric about both its diagonals. If $A$ is symmetric, and $J$ is the exchange matrix (reverse operator), then $A=A^T$ and $JA=AJ$. (You might look for symmetric centrosymmetric or or symmetric persymmetric algorithms.)

But where the inverse of a bisymmetric matrix is bisymmetric, the inverse of a Toeplitz matrix is not Toeplitz (in general). I mention this because the eigenvectors of the inverse of a non-singular matrix are the same as for the matrix itself (with eigenvalues the reciprocals of those for the matrix itself).

YMMV. I honestly don't know if there are specialized algorithms for the bisymmetric eigenvalue problem.

Answer 3

You say that you can compute the exact inverse of a Toeplitz matrix in O(n^2) time and space but you may know already that you only need O(n) space providing you can use each row of the inverse in the order it is generated and then overwrite it.