Why is the approximation valid, in the formula provided by Brian Karis?
"the information that $f(l_k,v)$ carries is low-frequency enough": As IneQuation explains, low-frequency was used to refer to the detail of the brdf function. I did actually mean that $f_r$ was low frequency though (which is the case with diffuse lighting), not $L_i$.

"with respect to $L_i$": what this means it that, since there aren't any large peaks in $f_r$, no part of $L_i$ will be wheighted much more heavily than any other in the integral.

If each function was, say, reciprocal of the other (meaning in this case that the peaks in each one would cancel out when multiplied by the other) the correct integral would be very different to the separated integral where this cancelation does not occurr. but since the two functions are not correlated (or they shouldn't be) this difference is likely to be fairly small.

Also, if you use a really high frequency BRDF that takes into acount only incident light around the normal and you separate it, you will get something that's way off. in this example the estimate using separation is 4x larger that the correct result: https://www.desmos.com/calculator/c4zgypthrc

Experimental results can be generated easily. here is an example: https://www.desmos.com/calculator/poyswojjsa it shows a low frequency brdf $f_r$ and a high frequency incident light $L_i$ function being multiplied before and after integration. the ratio between the correct result and the one obtained using separation is of about 1.6

Addendum regarding your latest comment on my other answer:

$$\frac{1}{N}\sum^{N}_{k=1}f(L_k) = \sum^{N}_{k=1}(f(L_k)\cdot \frac{1}{N})$$

Here we have "moved" or distributed $\frac{1}{N}$ into the sum. This formulation is equivalent but looks a bit more similar to an integral. if we set $\Delta L = \frac{1}{N}$ then, we can get

$$\sum^{N}_{k=1}f(L_k)\Delta L\approx \int_{L}f(L)dL$$

keeping in mind that the separation that occurs is, conceptually, separation of integrals, and that every integral "needs" its own $dL$. it follows that every montecarlo estimator needs its own $1/N$

I hope this helps you understand where the extra$\frac{1}{N}$ comes from