I read here that "a randomly initialized filter acts very much like an edge detector!". I want to know if there are any papers describing and explaining this phenomenon.
Intuition, for a small case
Why? Let's look at the simplest possible case, where the kernel is 1x2 (i.e., two pixels wide and one pixel high).
Here's a kernel matrix for an edge detector that detects vertical edges:
$$E_1 = [-1 \; +1]$$
Here's another matrix for an edge detector that also detects vertical edges:
$$E_2 = [+1 \; -1]$$
In particular, convolution with $E_1$ responds strongly to edges that are vertical and where the right side of the edge has higher intensity and the left side has lower intensity (whiter on the right, blacker on the left). $E_2$ responds strongly to vertical edges with the opposite intensity variation (whiter on the left, blacker on the right).
We can see that the same is true for other examples. For instance, consider
$$E_3 = [+1.7 \; -1.7]$$
This behaves almost identically to $E_2$, except for rescaling at its output.
For simplicity, let's focus on convolution kernels that have been normalized so their entries sum to zero. (Many convolution kernels we use in practice have this form, because they have nice properties.)
What about a random 1x2 matrix? Well, if it has been normalized, then it will necessarily have the form
$$M = [+\alpha \; -\alpha]$$
for some constant $\alpha$ (where $\alpha$ is random). Regardless of what $\alpha$ is, we see that this behaves like either $E_1$ or $E_2$.
Conclusion: a random, normalized 1x2 convolution filter behaves like a vertical edge detector, with high probability.
This can be generalized to other sizes of kernels. For a random, normalized 2x1 kernel, we'll basically get a horizontal edge detector. A random, normalized 2x2 kernel will likely respond to edges in some direction (maybe at some angle; the specific direction will depend on the entries of the matrix).
As the matrix gets larger and larger, this effect diminishes (I think). However, typically in practice we use a relatively small convolution kernel (a fairly small receptive field), so we really care mostly about the behavior with small convolution kernels.
What about normalization? So far I've been talking about matrices that were generated randomly and then normalized to sum to zero. What if we skip the normalization step?
Well, this doesn't change very much. Without normalization, a $m \times n$ random matrix $M$ can be decomposed into $M=M' + c \cdot A$ where $M'$ is normalized so its entries sum to zero, $c$ is a constant, $A$ is the matrix containing $1/(mn)$ in every entry.
Convolving the input image with $M'$ has a good chance of being basically applying an edge detector filter to the input image (as discussed above). Convolving the input image with $c \cdot A$ will basically blur the image and then multiply by $c$. Thus, the output image is the (weighted) sum of these two: the (weighted) sum of an edge detector filter plus a blurred version of the image.
For a random matrix, there's a decent chance that $c$ will be relatively small (since the mean of the sum of a bunch of random variables is often fairly small), so multiplying by $c$ makes the output image be determined more by the edge detector than by the blurring. If $c$ is small, we can basically ignore the second term of the sum and approximate the output image as the result of an edge detector filter.
That's why convolution with a random matrix has a decent chance of looking like applying some kind of edge detector to the input image.
Looking at the example images on those pages, to my eyes the first one looks more or less like an edge detector, and the second one looks more or less like a blur operator. This is consistent with the analysis above.