Solving Ax=b by FFT I read in Wiki that it is possible to solve Ax=b via Fast Fourier Transform given that A is a circulant matrix. For example, I have  $\begin{bmatrix}
1 & 0 & 0 & -1 \\
-1 & 1 & 0 & 0 \\
0 & -1 & 1 & 0 \\
0 & 0 & -1 & 1 
\end{bmatrix}x=$$ \begin{bmatrix}
1 \\
0 \\
0 \\
-1 
\end{bmatrix}  $
where A is a circulant matrix.
How do I solve this system using FFT? I had a hard time understanding the concept. Any help is greatly appreciated. Thank you for your understanding.
 A: 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
A: 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.
A: 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.
