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.


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,


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


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.

  • 2
    $\begingroup$ "... which makes me expect that the solutions to this problem won't be substantially faster than the solution of the problem for general matrices" You are mixing theoretical and practical questions here! In theory, it's easy to get an $O(n^2\log n)$ algorithm using Lancoz tridiagonalization based on FFT and a QR like diagonalization of the tridiagonal matrix, but see scicomp.stackexchange.com/questions/2975/… and scicomp.stackexchange.com/questions/10842/… $\endgroup$ – Thomas Klimpel Jan 14 '15 at 9:33
  • $\begingroup$ One reason why I am mixing them up is that I am actually inclined to think that the OP primarily cares about the practical, not theoretical, questions. $\endgroup$ – Luboš Motl Jan 14 '15 at 15:21
  • $\begingroup$ I don't think that solving the system with Levinson's algorithm could work, because the algorithm is based on extending the solution for the smaller problem with one less row and column, but in that smaller problem that eigen values will probably differ, and so the system will have only the trivial solution of a zero-vector. $\endgroup$ – yairchu Feb 6 '17 at 20:04

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.


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.