Friday, June 19, 2009

Circulant matrices

Circulant matrices are freakin' awesome.

Numerical analysts like me like to use matrix multiplications to take derivatives. When we solve differential equations (especially partial differential equations) we have to invert these matrices. The thing that makes my job suck, as a numerical analyst, is that I often have to invert huge matrices; in my code, if I were to do this directly, I'd have to invert a 32768x32768 matrix for the most computationally inexpensive case. Some models require inverting matrices of tens of millions of elements. This would take forever to do!

The matrix that represents differentiation for uniform periodic domains is a circulant matrix. (Circulant matrices are evidently a type of Toeplitz matrix.) Multiplying a circulant matrix by a vector is like taking the convolution of the first row of the matrix with that vector. It just so happens that F(u*v)=F(u)F(v), where F is the discrete Fourier transform, u and v are vectors of equal length, and * is the convolution operator. The Fourier representation of a single row of my difference matrix is easy to find analytically.

This gives me solutions to Navier-Stokes equation in O(n log n) time instead of O(N^3) time (Gaussian elimination) or O(N^2) time (oodles of basic iterative solvers).

No comments:

Post a Comment