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 indices: Aij=f(ij) .

I use iterative solvers because size A is 40000x40000. Iterative solvers need to determine the product A*x , where x is the test solution. Evaluation of this product is a convolution and therefore can be carried out using fast Fourier transforms (cputime ~ Nlog(N) instead of N^2 ). I have the following questions for this problem:

• Is this a circular bundle? Because, if it is circular, I believe that I need to use certain indexing for new matrices to take fft. It is right?

• I find it difficult to program the procedure for fft because I cannot figure out which indexing I should use. Is there any ready-made procedure that I can use to evaluate directly A*x product through fft, and not convolution? In fact, the matrix A is built of 3x3 blocks and is symmetrical. For me, the best solution would be a ready-made subprogramming for the product A*x .

• If there is no ready-made subroutine, could you give me an example on the example of how I could build this procedure for computing a matrix-vector product via fft?

Thanks in advance,

Panos

Answer 1

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

Circulating matrices.

The matrices corresponding to the cyclic convolution Ax = h*x (* is the convolution symbol) are diagonalized in the Fourier domain and can be solved using:

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

Triangular and striped.

triangular matrices and diagonally dominant striped matrices are solved efficiently using 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 effectively solved using multigrid methods (for example, web search for a "multilrid" matlab).

Answer 2

Interest Ask!

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 can do this with conv (keeping only the non-zero part with the valid option), which is probably also N * log (N). For example, let

 A = [abcd eabc feab gfea]; 

Then A*x coincides with

 conv(fliplr([gfeabcd]),x,'valid').' 

Or, as a rule, A*x matches

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

Answer 3

I would like to add some comments to Pio_Koon's answer.

First of all, I would not recommend following the proposal of triangular and ribbon matrices . The time taken to call the Matlab lu() procedure on a large sparse matrix massively outshines any of the benefits obtained by solving a linear system as x=U\(L\b) .

Secondly, in the Poisson problem you get a circulant matrix, so you can solve it using FFT, as described. In this particular case, your convolutional mask h is Laplacian, i.e. h=[0 -0.25 0; -0.25 1 -0.25; 0 -0.25 0] h=[0 -0.25 0; -0.25 1 -0.25; 0 -0.25 0] .