Are there mathematical reasons for convolution in neural networks beyond expediency? In convolutional neural networks (CNN) the matrix of weights at each step gets its rows and columns flipped to obtain the kernel matrix, before proceeding with the convolution. This is explained on a series of videos by Hugo Larochelle here:

Computing the hidden maps would correspond to doing a discrete
  convolution with a channel from the previous layer, using a kernel
  matrix [...], and that kernel is computed from the hidden weights
  matrix $W_{ij}$, where we flip the rows and the columns.


If we were to compare the reduced steps of a convolution to regular matrix multiplication as in other types of NN, expediency would be a clear explanation. However, this might not be the most pertinent comparison...
In digital imaging processing the application of convolution of a filter to an image (this is a great youtube video for a practical intuition) seems related to: 


*

*The fact that convolution is associative while (cross-)correlation is not. 

*The possibility to apply filters in the frequency domain of the image as multiplications, since convolution in the time domain is equivalent to multiplication in the frequency domain (convolution theorem).



In this particular technical environment of DSP correlation is defined as:
$$F\circ I(x,y)=\sum_{j=-N}^{N}\sum_{i=-N}^N\, F(i,j)\,I(x+i, y+j)$$
which is essentially the sum of all the cells in a Hadamard product:
$$\small F\circ I(x,y)=\Tiny\begin{bmatrix}F[-N,-N]\,I[x-N,y-N]&\cdots&F[-N,0]\,I[x-N,y-N]&\cdots& F[-N,N]\,I[x-N,y+N]\\
\vdots&\ddots&\vdots&\ddots&\vdots\\
F[0,-N]\,I[x,y-N]&\cdots&F[0,0]\,I[x,y]&\cdots& F[0,N]\,I[x,y+N]\\
\vdots&\ddots&\vdots&\ddots&\vdots\\
F[N,-N]\,I[x+N,y-N]&\cdots&F[N,0]\,I[x+N,y]&\cdots& F[N,N]\,I[x+N,y+N]\\
\end{bmatrix}$$
where $F(i,j)$ is a filter function (expressed as a matrix), and $I(x,y)$ is the pixel value of an image at location $(x,y)$:

The objective of cross-correlation is to assess how similar is a probe image to a test image. The calculation of a cross-correlation map relies on the convolution theorem.

On the other hand, convolution is defined as:
$$F* I(x,y)=\sum_{j=-N}^{N}\sum_{i=-N}^N\, F(i,j)\,I(x-i, y-j)$$
which as long as the filter is symmetric, it is the same as a correlation operation with the rows and columns of the filter flipped:
$$\small F* I(x,y)=\Tiny\begin{bmatrix}F[N,N]\,I[x-N,y-N]&\cdots&F[N,0]\,I[x-N,y-N]&\cdots& F[N,-N]\,I[x-N,y+N]\\
\vdots&\ddots&\vdots&\ddots&\vdots\\
F[0,N]\,I[x,y-N]&\cdots&F[0,0]\,I[x,y]&\cdots& F[0,-N]\,I[x,y+N]\\
\vdots&\ddots&\vdots&\ddots&\vdots\\
F[-N,-N]\,I[x+N,y-N]&\cdots&F[-N,0]\,I[x+N,y]&\cdots& F[-N,-N]\,I[x+N,y+N]\\
\end{bmatrix}$$


Convolution in DSP is meant to apply filters to the image (e.g. smoothing, sharpening). As an example, after convolving Joseph Fourier's face with a Gaussian convolution filter: $\small\begin{bmatrix}    1&4&7&4&1\\
    4&16&26&16&4\\
    7&26&41&26&7\\
    4&16&26&16&4\\
    1&4&7&4&1\end{bmatrix}$ the edges on his face are fuzzier:


Computationally, both operations are a Frobenius inner product, amounting to calculating the trace of a matrix multiplication. 

Questions (reformulating after comments and first answer):


*

*Is the use of convolutions in CNN linked to FFT?

From what I gather so far the answer is no. FFTs have been used to
  speed up GPU implementations of
  convolutions. However, FFT are not usually part of the structure or activation functions in CNNs, despite the use of convolutions in the pre-activation steps. 



*Is convolution and cross-correlation in CNN equivalent?

Yes, they are equivalent.



*If it is a simple as "there is no difference", what is the point of flipping the weights into the kernel matrix?

Neither the associativity of convolution (useful in math proofs), nor any considerations regarding FTs and the convolution theorem are applicable. In fact, it seems as though the flipping doesn't even take place (cross-correlation being simply mislabeled as convolution) (?).  

 A: There are no differences in what neural networks can do when they use convolution or correlation. This is because the filters are learned and if a CNN can learn to do a particular task using convolution operation, it can also learn to do the same task using correlation operation (It would learn the rotated version of each filter).
To find more details about the reasons that people sometimes find it more intuitive to think about convolution than correlation, this post may be useful.
There remains this question that if there is no difference between convolution and cross-correlation, what is the point of flipping the weights into the kernel matrix? I would like to include some sentences from the Deep learning book by Ian Goodfellow et al. to answer this question:

The only reason to flip the kernel is to obtain the commutative property. While the commutative property is useful for writing proofs, it is not usually an important property of a neural network implementation... Many machine learning libraries implement cross-correlation but call it convolution.

The takeaway is that although convolution is a favorite operation in classic machine vision applications, it is replaced by correlation in many of the implementations of the convolutional neural networks.
A: There's a practical reason for the link between FFTs and convolution.
Convolution is slow in the time/image domain. Applying an $n \times n$ filter to one pixel requires $O(n^2)$ multiplications and additions. Applying it to every pixel in an $N \times N$ image thus requires $n^2N^2$ operations. This grows quickly, and the large number of operations not only requires extra time, but also introduces more numerical error. 
The Convolution Theorem says that convolution in the time domain is equivalent to pointwise multiplication in the frequency domain. FFTs are fast: they have good asymptotic performance $O(N^2 \log N^2)$ and the actual implementations are often highly optimized. Switching to the Fourier domain thus lets you perform a convolution in $O(N^2)$ time (which is dominated by the pointwise multiplication), instead of $O(n^2N^2)$. This can provide a considerable speedup, even though it seems much more complicated to go down the FFT -> multiplication -> inverse FFT route. More here
