Intuitive Idiot-Proof Validation of Importance Sampling for Monte Carlo Computations (video)

Do not close the tab if you're not into Raytracing, the method I present is applicable to every algorithm using Importance Sampling even beyond Monte Carlo methods, such as Gaussian Quadratures with IS.

So you have an Importance Sampling routine (i.e. a sample transform from the uniform distribution) and you wish to validate that its actually producing samples according to the Probability Density Function that you (or whoever implemented the routine) claim it matches.

How can you make sure that this is actually the case, beyond just looking at the results of a simulation with your new routine and seeing if they seem "off"?

Well there exist t-tests, Chi Squared test, or simply histogramming the resulting samples into discrete bins and hoping that with wide enough bins and enough samples the average occupancy will resemble the PDF.

However the t-test and Chi Squared will only give you a single metric, with confidence or perceived "error" being highly dependant on numerical methods and precision. While binning your samples will take a really long time to converge to something resembling your PDF, especially if it's multidimensional!

There exists another way, ridiculously simple, easy and intuitive to grasp, stemming from a little application of elementary Multivariable Calculus. Yet for some reason academic literature sidesteps the issue of testing Importance Sampling implementations (sure they give proofs that you can work out by hand, but there's not a lot of interest in automated testing, which is strange because you can bet on a first implementation containing a mistake).

See, while in the process of producing samples (be they random or LDS for QMC) we deal with discrete values and points, the actual functions transforming samples from the (usually) uniform distribution are (usually) continuous.

Long story short, the absolute determinant of the Jacobian of your Importance Sampling's sample transform function multiplied by the pdf evaluated at the transformed sample, should be equal to 1.0

The most epic thing is that it works with all possible mappings (like VNDF sampling in Rendering), not only ones that have been developed via the Inversion Method of conditional probabilities.

Now you can evaluate that Jacobian:

1. Numerically (i.e. Finite Differences)
2. Programmatically and Automatically with a framework like Enoki
3. Using a symbolic manipulation package such as Maple, Mathematica or Sage Math
4. By hand (tedious and error prone)

We chose the first method because it allows for the implementation to be tested in-situ.

This is a very Idiot-Proof method (myself being one), and it is very uniform and robust, as long as you're not completely changing domains the code stays the same for testing importance sampling of many unrelated functions. In our case of Monte Carlo Rendering we can use the same code to validate GGX, Beckmann, and other BRDF's VNDF sampling routines. This is naturally a 2D problem, but you can apply the same logic to an arbitrarily dimensional Importance Sampling (however I recommend you run your test at locations given by a Low Discrepancy Sequence and not a uniform grid).

In other words absolute value of the Jacobian's determinant is inversely proportional to the PDF, however for IEEE754 robustness reasons I'd rather have you remember the first formulation. These reasons are outlined in the video.

That sort of insight is incredibly useful, especially if you're building new Importance Sampling routines that have discontinuities, because the mapping could fail only in one specific region. With a Chi Squared or other test that gives a single number unified metric of PDF and Importance Sampling agreement, it's impossible to see where your mapping works and where it doesn't.

I've initially wanted to write this down in a form of a format article/paper with some LaTeX, but because I needed to conduct a training session for my Mid/Senior Developer on the subject anyway, I've decided to record it instead of doubling my work.

P.S. I've made one claim in the video, that the Sample Transform needs to be a bijection, it does not. There are certain, but weaker constraints, mainly that the image of any element under the Importance Sampling function is continuous, and the function itself has a finite number of discontinuities (weaker condition that "monotonically increasing" or "monotonically decreasing). So this will continue to work with PDFs such as delta distributions. It will also work with PDFs that have zeros because your Importance Sampling should never generate a sample that has a PDF of 0, if it does then you have an error in the implementation and so you'll be comparing with a FLT_MAX and seeing a natural failure.