fft matrix-vector multiplication

Question

I have to solve in MATLAB a linear system of equations A*x=B where A is symmetric and its elements depend on the difference of the indices: Aij=f(i-j).

I use iterative solvers because the size of A is say 40000x40000. The iterative solvers require to determine the product A*x where x is the test solution. The evaluation of this product turns out to be a convolution and therefore can be done dy means of fast fourier transforms (cputime ~ Nlog(N) instead of N^2). I have the following questions to this problem:

1. is this convolution circular? Because if it is circular I think that I have to use a specific indexing for the new matrices to take the fft. Is that right?

2. I find difficult to program the routine for the fft because I cannot understand the indexing I should use. Is there any ready routine which I can use to evaluate by fft directly the product A*x and not the convolution? Actually, the matrix A is constructed of 3x3 blocks and is symmetric. A ready routine for the product A*x would be the best solution for me.

3. In case that there is no ready routine, could you give me an idea by example how I could construct this routine to evaluate a matrix-vector product by fft?

Thank you in advance,

Panos

Answer 1

Very good and interesting question! :) For certain special matrix structures, the Ax = b problem can be solved very quickly.

Circulant matrices.

Matrices corresponding to cyclic convolution Ax = h*x (* - is convolution symbol) are diagonalized in the Fourier domain, and can be solved by:

x = ifft(fft(b)./fft(h));

Triangular and banded.

Triangular matrices and diagonally-dominant banded matrices are solved efficiently by sparse LU factorization:

[L,U] = lu(sparse(A)); x = U\(L\b);

Poisson problem.

If A is a finite difference approximation of the Laplacian, the problem is efficiently solved by multigrid methods (e.g., web search for "matlab multigrid").

Answer 2

Interesting question!

The convolution is not circular in your case, unless you impose additional conditions. For example, A(1,3) should equal A(2,1), etc.

You could do it with conv (retaining only the non-zero-padded part with option valid), which probably is also N*log(N). For example, let

A = [a b c d
     e a b c
     f e a b
     g f e a];

Then A*x is the same as

conv(fliplr([g f e a b c d]),x,'valid').'

Or more generally, A*x is the same as

conv(fliplr([A(end,1:end-1) A(1,:)]),x,'valid').'

Answer 3

I'd like to add some comments on Pio_Koon's answer.

First of all, I wouldn't advise to follow the suggestion for triangular and banded matrices. The time taken by a call to Matlab's lu() procedure on a large sparse matrix massively overshadows any benefits gained by solving the linear system as x=U\(L\b).

Second, in the Poisson problem you end up with a circulant matrix, therefore you can solve it using the FFT as described. In this specific case, your convolution mask h is a Laplacian, i.e., h=[0 -0.25 0; -0.25 1 -0.25; 0 -0.25 0].