Solving linear system $Ax=b$ via FFT

Question

On Wikipedia, I read that it is possible to solve $Ax=b$ when the matrix $A$ is circulant via the Fast Fourier Transform (FFT). For example, I have

$$\begin{bmatrix} 1 & 0 & 0 & -1 \\ -1 & 1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} | \\ x \\ | \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ -1 \end{bmatrix} $$

where $A$ is a circulant matrix. How do I solve this system using the FFT? I had a hard time understanding the concept. Any help is greatly appreciated.

Answer 1

We can use something called a similarity transform $A = F^{-1}DF$

For circulant matrices we can be sure that we can choose $F$ to be a DFT matrix and $D$ matrix will be a diagonal matrix with the Fourier coefficients of one row (or column).

This $F$ matrix is dense so so far we have not saved much by doing this transformation.

The FFT in matrix language is nothing else but a factorization of this $F$ matrix.

We can write $F = F_1F_2\cdots F_N$ . A product of matrices.

Now the neat part is that we can get these $F_k$ matrices to be sparse. Only two non-zero values every row no matter the size of the matrix.

We can show that $N$ here will be $\log_2(n)$ where $n$ is the side of $A$. Quite few matrices compared to the size. So we will get away with all in all something like $n+4n\log(n)$ multiplications instead of $n^2$ multiplications if we are to do a matrix-vector multiplication.

Answer 2

If $x$ has coordinates $x_1,x_2,x_3,x_4$, your matrix-vector equation amounts to:

$$[1,-1,0,0]\star[x_1,x_2,x_3,x_4]=[1,0,0,-1]\tag{1}$$

This is a deconvolution issue.

Here is the (classical) way to treat it using FFT (one should use the term DFT, Discrete Fourier Transform, but usage is to blend the two terms).

Let us take the FFT of both sides of (1):

$$FFT([1,-1,0,0]) \cdot FFT([x_1,x_2,x_3,x_4])=FFT([1,0,0,-1])$$

which is equivalent to:

$$FFT([x_1,x_2,x_3,x_4])=\dfrac{FFT([1,0,0,-1])}{FFT([1,-1,0,0])}\tag{2}$$

(division in the RHS being understood componentwise)

It remains now to apply the IFFT (Inverse FFT) to (2) to get $[x_1,x_2,x_3,x_4].$

Remark: I met deconvolution in image processing, in connection with "filters" giving rise to Wiener deconvolution.

Answer 3

As Jean-Marie explains, such matrix-vector product implements a circular convolution.

The vector $\mathbf{x}$ contains the filter coefficients (so I will wrote it $\mathbf{h}$ in the following).

We may interpret $\mathbf{b}$ as the output of a finite-impulse response (FIR) filter applied to the 4-periodic signal $\mathbf{a}=[1;-1;0;0]$. When filtering, we must consider $\mathbf{a}_{per}=[-1;0;0;\ 1;-1;0;0]$ where the 3 samples to the left have been added because of the periodicity..

Filtering the vector $\mathbf{a}_{per}$ with the (unknown) filter $\mathbf{h}$ yields the output $\mathbf{b}$ $$ \left\lbrace \begin{array}{ccc} b_1 = (-1)*h_4+0*h_3+0*h_2+1*h_1 = h_1 - h_4 \\ b_2 = 0*h_4+0*h_3+1*h_2+(-1)*h_1 = h_2 - h_1 \\ b_3 = 0*h_4+1*h_3+(-1)*h_2+0*h_1 = h_3 - h_2 \\ b_4 = 1*h_4+(-1)*h_3+0*h_2+0*h_1 = h_4 - h_3 \end{array} \right. $$

2)Solving for $\mathbf{h}$ is easy if you recognize that $\mathbf{A}$ is a circulant matrix. Indeed the DFT matrix $\mathbf{W}$, wiki will diagonalize $\mathbf{A}$. $$ \mathbf{A} = \mathbf{W}^{-1}\mathrm{diag}(\mathbf{f}_a)\mathbf{W} $$ with $\mathbf{f}_a$ the FFT of the signal $\mathbf{a}$

The matrix-vector product $\mathbf{W}\mathbf{z}$ is implemented using fft(z) while the matrix-vector product $\mathbf{W}^{-1} \mathbf{z}$ is implemented using ifft(z)

The solution writes $$ \mathbf{x} = \mathbf{W}^{-1} \mathrm{diag}^{-1}(\mathbf{f}_a) \mathbf{W} \mathbf{b} $$ In Matlab, the solution is found using z=fft(b)./fft(a); xsol = ifft(z)

In your specific case, there is a problem since we are dealing with zero-mean vectors (you have a divide by 0 operation). So before IFFT we force the DC to be null...

z=fft(b)./fft(a);z(1)=0;xsol = ifft(z)

The exact solution is [1;1;1;-3]/4