I finally found some more time to work on my symbolic regression library (tentatively named CuSR) by simply procrastinating on my other ongoing projects. My previous post on this can be found here. Since adding features is much more rewarding when applying them to real problems, I decided to search for alternatives to the Trowbridge-Reitz (GGX) distribution function, originally introduced in the paper Average irregularity representation of a rough surface for ray reflection and later popularized as the GGX distribution in the paper Microfacet Models for Refraction through Rough Surfaces. Now it’s being used throughout computer graphics.
A commonly used form for the Trowbridge-Reitz distribution for $\theta \leq \pi/2$ is:
$$ D(\theta) = \frac{\alpha_g^2}{\pi (1+\cos^2\theta(\alpha_g^2-1))^2} \tag{1}$$
$\theta$ is the angle between the normal vector $\vec{n}$ and the halfway vector between the view and light vector $\vec{h}$.
$\alpha_g$ controls the shape of the distribution and can be considered a roughness parameter when the distribution is used in BRDFs.
Alternatives to this function should follow some rules to be viable (for any $\alpha_g \in [0,1]$):
$$D(\theta) \geq 0 \tag{2}$$
$$\text{Monotonically non-increasing} \tag{3}$$
$$\int_\Omega cos\theta_hD(\theta_h) d\omega_h = 1 \tag{4}$$
$$D(\theta) = C\text{ for }\alpha_g = 1 \tag{5}$$
$$\alpha_1 < \alpha_2 \Rightarrow \int_0^{\pi/2} \theta_hD_{\alpha_1}(\theta_h) d\theta_h < \int_0^{\pi/2} \theta_hD_{\alpha_2}(\theta_h) d\theta_h\ \tag{6}$$
Equation (4) ensures that the distribution integrates to 1 over the hemisphere $\Omega$ when weighted by $cos\theta_h$. This is equal to Equation (4) in the GGX paper and represents only a special case, the general case was manually verified. Equation (5) is not strictly necessary, but it ensures that for fully rough surfaces ($\alpha_g = 1$), the distribution becomes a constant, similar to a Lambertian distribution. This property aligns with the behavior of the Trowbridge-Reitz distribution. Equation (6) ensures that an increasing $\alpha_g$ results in a broader lobe.
I implemented numerical checks for these constraints and evaluated approximately 420 million candidate functions using the operators $+,-,*,/,^2,^3$. Some care was taken to not generate redundant expressions. On my RTX 4070 Ti, I can evaluate approximately 20 million candidates per second on 512 data points (128 different $\theta$ and 4 different $\alpha_g$). This results in a runtime of just 21 seconds and I would have liked to let it run longer, but for now the amount of expressions per run is still limited by the main memory.
To remove duplicates, I first used a SymPy-based script, then processed the remaining expressions with ChatGPT o3-mini, followed by a manual review. Ultimately, only two strikingly similar functions remained: the already-known Trowbridge-Reitz distribution and a single novel alternative that also satisfies all necessary criteria:
$$ D(\theta) = \frac{\alpha_g^4}{\pi (1+\cos\theta(\alpha_g^2-1))^3} \tag{7}$$
Although the overall yield of this exercise was lower than expected, this function appears to be a viable alternative to existing distributions. I have not found it documented in the literature – if you are aware of any references, please let me know. Its shape lies between the Trowbridge–Reitz and Beckmann distributions. A qualitative comparison is beyond the scope of this post and will depend on the specific dataset being approximated.

Its CDF and inverse CDF are (where $\xi$ is a uniform random number in $[0,1]$):
$$CDF(\theta) =\frac{\frac{\alpha_g^4}{(1+(\alpha_g^2-1)\cos\theta)^2}-1}{\alpha_g^4-1}$$
$$\theta =arccos\left(\frac{1}{\alpha_g^2-1}\left(\frac{1}{\sqrt{\xi\left(1-\frac{1}{\alpha_g^4}\right)+\frac{1}{\alpha_g^4}}}-1\right)\right)$$
Similarly to the “Generalized-Trowbridge-Reitz (GTR)” introduced in Physically Based Shading at Disney, the above distribution function can also be used in a generalized form:
$$ D(\theta) = \frac{C}{(1+\cos\theta(\alpha_g^2-1))^\gamma} \tag{8}$$
$$ 1/C = \frac{2\pi}{(\alpha_g^2-1)^2}\left(\frac{\alpha_g^{4-2\gamma}-1}{2-\gamma}-\frac{\alpha_g^{2-2\gamma}-1}{1-\gamma}\right) \tag{9}$$
Note that, similarly to GTR, a higher $\gamma$ will lead to a smaller lobe for the same roughness values.