Actually, the given pipeline was used in the old days of Graph Neural Networks.
Canonical paper on the subject is Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering.
You start from arbitrary graph with adjacency matrix $A_{ij}$ (let us assume that graph is undirected), such that:
$$
A_{ij} =
\begin{cases}
1 & \text{if there and edge between vertices $i$ and $j$} \\
0 & \text{otherwise}
\end{cases}
$$
Then one constructs graph Laplacian:
$$
L = D-A
$$
$D$ is the degree matrix (number of edges entering the given vertex). There are also different normalizations in the literature, like $I - D^{-1/2} A D^{-1/2}$ is called normalized Laplacian.
This matrix has several properties:
- It is symmetric
- Non-negative definite
From the first statement it follows, that the matrix can be diagonalized due to the Spectral theorem,
Therefore, it makes sense to perform the eigendecomposition of this operator. And the eigenvectors form the graph Fourier basis.
Note, that in the special case, when the graph is a regular square grid, graph Laplacian just the discrete Laplace operator:
$$
\begin{pmatrix}
0 & -1 & 0 \\
-1 & 4 & -1 \\
0 & -1 & 0 \\
\end{pmatrix}
\qquad
(\text{for 2d case})
$$
And the Fourier basis consists of plane waves:
$$
\sim e^{i (k_i i + k_j j)}
$$
The next important fact is the Convolution theorem that states that convolution of two signals can be done as inverse Fourier transform of the dot product of their Fourier transforms:
$$
f * g = \mathcal{F}^{-1} [\mathcal{F}[f] \cdot \mathcal{F}[g]]
$$
These operations correspond to steps 3 and 4 in the pipeline in the OP.
I am not aware of the simple geometrical intuition in this setting. But in the following research, full eigendecomposition was truncated to Chebyshev polynomials, and then up to the first term in the decomposition, which gave rise to Graph Convolutional Networks by T.Kipf.
They allow for more intuitive and visual interpretation. Given adjacency matrix $A$, input feature map on the graph $H^{(l)}$ one defines the output of the layer to be:
$$
f(H^{(l)}, A) = \sigma(\tilde{D}^{-1/2} A \tilde{D}^{-1/2} H^{(l)} W^{(l)})
$$
$\tilde{D}$ is the graph Laplacian for adjacency matrix with self-loops (arrows $i \rightarrow i$). $W^{(l)}$ is the matrix of learnable parameters.
In essence, one transforms the feature vectors $ H^{(l)}$ by a pointwise linear transformation, and aggregates the information from the neighborhood with some coefficients. Equivalently the function above can be rewritten as;
$$H_{i}^{l+1}=\eta\left(\frac{1}{\hat{d}_{i}} \sum_{j \in N_{i}} \hat{\boldsymbol{A}}_{i j}{W}^{l} H_{j}^{l}\right)$$
From the formula above one can see, that graph convolution is the summation of the features in the neighborhood vertices which are transformed by the weight matrix ${W}^{l}$ and division by the normalizing constant.
