  The large $N$ limit of SU($N$) gauge theories is well understood in
perturbation theory. Also non-perturbative lattice studies have yielded
important positive evidence that 't Hooft's predictions are valid. We go far
beyond the statistical and systematic precision of previous studies by making
use of the Yang-Mills gradient flow and detailed Monte Carlo simulations of
SU($N$) pure gauge theories in 4 dimensions. With results for $N=3,4,5,6,8$ we
study the limit and the approach to it. We pay particular attention to
observables which test the expected factorization in the large $N$ limit. The
investigations are carried out both in the continuum limit and at finite
lattice spacing. Large $N$ scaling is verified non-perturbatively and with high
precision; in particular, factorization is confirmed. For quantities which only
probe distances below the typical confinement length scale, the coefficients of
the $1/N$ expansion are of $\mathrm{O}(1)$, but we found that large (smoothed)
Wilson loops have rather large $\mathrm{O}(1/N^2)$ corrections. The exact size
of such corrections does, of course, also depend on what is kept fixed when the
limit is taken.
