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.