How Does a Convolution Can Be Expressed as a Matrix Multiplication (Matrix Form)? I know this question may not be very relevant to programming, but if I don't understand the theory behind image processing I'll never be able to implement something in practice.

If I got it right Gaussian filters are convolved with an image for noise reduction since they compute a weighed average of a pixel's neighborhood and they are very useful in edge-detection, since you can apply a blur and derive the image at the same time by simply convolving with the derivative of a Gaussian function. 

But can anyone explain me, or give me some references on how are they computed?

E.g. Canny's edge detector talks about a 5x5 Gaussian filter, but how did they get those particular numbers? And how did they go from a continuous convolution to a Matrix multiplication?
For this operation to work, you need to imagine that your image is reshaped as a vector.
Then, this vector is multiplied on its left by the convolution matrix in order to obtain the blurred image.
Note that the result is also a vector the same size as the input, i.e., an image of the same size.
Each row of the convolution matrix corresponds to one pixel in the input image.
It contains the weight of the contributions of all the other pixels in the image to the blurred counterpart of the considered pixel.
Let's take an example: box blur of size $3 \times 3$ pixels on an image of size $6 \times 6$ pixels.
The reshaped image is a column of 36 elects, while the blur matrix has size $36 \times 36$.

Let's initialize this matrix to 0 everywhere.
Now, consider the pixel of coordinates $(i,j)$ in the input image (not on its border for simplicity). Its blurred counterpart is obtained by applying a weight of $1/9$ to itself and each of its neighbours at the positions $(i-1,j-1); (i-1,j), (i-1,j+1),\ldots,(i+1,j+1)$.
In the column vector, the pixel $(i,j)$ has the position $6*i + j$ (assuming lexicographic ordering). we report the weight $1/9$ in the $(6i+j)$-th line of the blur matrix.
Do the same with all other pixels.

A visual illustration of a closely related process (convolution + subtraction) can be found on my personal blog Computers Don't See (Yet) - Linearized model of LBPs.