The way I typically look at this is by relating the maximal difference to the magnitude of the results (0.000168/139) in order to avoid the problem that relative difference has with cancellation (or you could look at only positive inputs). This gives about 1.2e-6, which would seem standard fare for float32 computation.
A thing to “know” – well, let’s say anecdotally be aware of – is that efficient methods to compute convolutions - prominently the Winograd method - may be more prone to see numerical precision issues than the naïve counterparts (and im2col is just rearranging the inputs to the defining formula). Intuitively, this isn’t surprising because it cleverly exploits algebraic relations between partial results which only works approximately in floating point arithmetic.
Part of the price you pay for im2col in is the materialization of very large matrices, which is costly in terms of memory but also time (because of limited memory bandwidth).
If you’re very determined, you could try to get a good kernel for direct convolution from TVM or use the algorithm selection of the MIOpen / CuDNN libraries to your advantage.
One way I always use to verify that everything is all right in general is to use float64 as a reference. (Provided both use the same code path.)
Best regards
Thomas