The answer to your question - Is convolution linear is that - 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]=∑k=∞∞x[k]h[n−k]
which has the nice property of
Y(ejω)=X(ejω)H(ejω)
where each A(ejω) 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 is 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]=∑j=−∞∞∑k=−∞∞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 the frequency domain. The convolution operation is a little bit different due to this periodicity.
So these are the actual definitions of convolution. It has a linearity property (clearly obvious since it is used to calculate output of an LTI system), and it also has linearity in terms of matrix multiplications, since we want 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 yn=hn∗xn: where
yt=⎡⎣⎢⎢⎢y0y1y2y3⎤⎦⎥⎥⎥=⎡⎣⎢⎢⎢h0h1h2h3h3h0h1h2h2h3h0h1h1h2h3h0⎤⎦⎥⎥⎥⎡⎣⎢⎢⎢x0x1x2x3⎤⎦⎥⎥⎥
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 h0 we will have a Hermitian matrix H0 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 the frequency domain.
