Convolution is a pretty misused term in recent times with the advent of CNN.
Short answer, Convolution is a linear operator (check here) but what you are defining in context of CNN is not convolution, it is cross-correlation which is also linear in case of images (dot product).
Convolution:
Computers work only on discrete values, so the discrete time convolution is defined similarly as:
$$
y[n] = \sum_{k=\infty}^{\infty} x[k] h[n-k]
$$
which has the nice property of
$$Y(e^{j\omega}) = X(e^{j \omega})H(e^{j \omega})$$
where each $A(e^{j\omega})$ is the Fourier Transform of $a(t)$
So this is the basic idea of discrete convolution. The point to illustrate here is that convolution operation has some nice properties and very helpful.
Now there are 2 main types of convolution (as far as I know). One is linear convolution another is circular convolution. In terms of $2D$ signal they are defined in the following ways:
$$y[m,n] = \sum_{j=-\infty}^{\infty}\sum_{k=-\infty}^{\infty} x[i,j]h[m-i,n-j]$$
Circular convolution is the same except that the input sequence is finite from $0$ to $N-1$ giving rise to periodicity in frequency domain. The convolution operation is a little bit different due to this periodicity.
So these are the actual definitions of convolution. It has linearity property (clearly obvious since it is used to calculate output of LTI system), and it also has the linearity in terms of matrix multiplications, since we want to do the computer to do these calculations for us. I don't know a lot but there are many clever manipulations e.g the FFT algorithm an indispensable tool in signal processing (used to convert signals to frequency domain for a certain sampling rate). Like this you can define convolutions in terms of Hermitian matrices and stuff (only if you want to process in the $n$ domain it is much easier to process in the frequency domain).
For example a $1D$ circular convolution between 2 signals in $n$ domain defined in matrix form, can be written as follows $y_n=h_n*x_n$:
where
$$y_t = \begin{bmatrix} y_0 \\
y_1 \\
y_2 \\
y_3
\end{bmatrix} = \begin{bmatrix}
h_0 & h_3 & h_2 & h_1\\
h_1 & h_0 & h_3 & h_2\\
h_2 & h_1 & h_0 & h_3 \\
h_3 & h_2 & h_1 & h_0
\end{bmatrix} \begin{bmatrix} x_0 \\
x_1 \\
x_2 \\
x_3
\end{bmatrix}$$
The same can be done quite easily by converting the functions into frequency domain, multiplying and converting back to time domain.
When we move to $2D$ domain similar equation is formed except (in $n$ domain)instead of $h_0$ we will have a Hermitian matrix $H_0$ and we will perform the Kronecker product (I don't know the justification or proofs of these equations, it probably satisfies the convolution theorem and is fast when ran on computers). Again much easier to do in frequency domain.
When we move to multiple dimensions, it is called Multidimensional discrete convolution as per this Wikipedia article. As the article suggests the property where everything is in frequency domain.
$$Y(k_1,k_2,...,k_M)=H(k_1,k_2,...,k_M)X(k_1,k_2,...,k_M)$$
still holds good. When we do convolutions in the $n$ domain things get tricky and we have to do clever matrix manipulations as shown in the example above. Whereas, as stated above things get much easier in the frequency domain.
Its counter-intuitive that a picture has a frequency domain, and in general its actually not frequency. BUT in DSP we use a lot of filters whose math properties are similar to those of filters in the traditional sense of frequency and hence has the same calculations as in a frequency setting as shown in the 1st example of $1D$ signal.
The point convolution by definition is a linear operation. Check here if the explanations are unconvincing. There are Linear Time Varying systems and its output maybe determined by convolution, but then it uses the equation at a certain point of time i.e:
$$Y(e^{j\omega}, t) = H(e^{j\omega}, t)X(e^{j\omega} )$$.
Now whether it can be represented by matrix products in $n$ domain or not I cannot say. But generalising, it should be possible, only it will involve increasingly complex matrix properties and operations.
Cross Correlation:
Cross correlation is nothing similar to convolution, it has a bunch of interesting properties on its own, hence I do not like the two being mixed together. In signal processing it is mainly related to finding energy (auto-correlation) and has other applications in communication channels.
While in statistics, it is related to finding correlation (misused term since it can take any value, but correlation coefficient takes value between -1 to 1) between 2 stochastic processes, at different time instants, where we have multiple samplings of the signal or ensembles of the signal based on which the expectation is taken (See here). In CNNs maybe a way to see this would be, the image dataset is a ensemble of stochastic process, while the filter is another process which is fixed. And moving the filter over the image is time delay.
In the context of digital image processing cross correlation is just a dot product. And hence it can be represented by matrices. Whether cross correlation is linear or not is difficult to tell (at least I don't know), but in the context of image processing it seems linear. Bottom-line is that it can be definitely implemented by matrices as it is a dot product, whether it is linear in true mathematical sense is doubtful.
The more intuitive way in context of CNNs would be to view filters as just a bunch of shared weights for better regularisation instead of cross correlation or convolution.
