Paper for the approximation formula provided by Brian Karis In these slides, specifically page 11, the following formula is reported:

$$
\frac{1}{N} \sum_{k=1}^N \frac{L_i(l_k)f(l_k,v)\cos(\theta_{l_k})}{p(l_k,v)} \approx \left( \frac{1}{N} \sum_{k=1}^N L_i(l_k) \right) \left( \frac{1}{N} \sum_{k=1}^N  \frac{f(l_k,v)\cos(\theta_{l_k})}{p(l_k,v)} \right)
$$

The formula above should be an approximation of the rendering equations. I'm looking for a paper that explains how that formula is derived, I cannot manage to find neither the title of a paper (other than the name "Dimitar" mentioned in the same page I pointed out).

I've done some research about some possible method that would justify. This is the closest thing I managed to find, however there's some hypothesis that doesn't seem to me fits.
I do not have a derivation but i can explain the reasoning behind its usage.

On the left side is a (monte carlo) estimator for the rendering equation. It states the very same thing that the rendering equation states but on a finite number of discrete samples.

on the right side we have the $L_i(l_k)$ factor separated and "integrated" on its own.
This approximation holds true when $L_i(l_k)$ is a constant, since $\sum_{k}f(k)\times n = n \times\sum_{k}f(k)$ - and, likewise, $\int_{k}f(k)\times n\space dk = n\times\int_{k}f(k)dk$ - when $n$ is a constant.

$L_i(l_k)$ can be though of the colour of the incoming light and $\int_{k}f(l_k,v)dk$ can be thought of as the intensity of the reflected light. multiplying those together you get that now you can precompute $\int_{k}L_i(l_k)dk$, saving a lot of time in either texture lookups (in the case of sky / environment lighting / screen space ray marching) or raytracing calculations.

However, since $L_i(l_k)$ is usually not a constant, this integral is not separable. Yet, in some cases (like when dealing with diffuse lighting or, more specifically, diffuse ambient lighting), the information that $f(l_k,v)$ carries is low-frequency enough that it can reasonably be approximated as a constant with respect to the distribution of $L_i(l_k)$. Experimentally, it can be shown that the final result is not too different from what you would expect.