Sampling the light vs. sampling the BRDF: difference in intensity
Shoot me... After the fixes I posted above, I started debugging Mitsuba to see whether it does anything differently on a large scale. Of course it didn't; my code is completely correct, the dimness is due to this macro def:

#define EPSILON   1e-5


This is used in divisions to avoid dividing with zero. Unfortunately, it accumulates to a large value during the calculations. Changing it to 1e-9 fixes everything:



For future generations I leave the source code up (just remember to change this value).

UPDATE:

I found another mistake in my calculations, but it isn't apparent in this scene. Microsurface-based BRDFs are usually sampled according to the microfacet normal distribution function (D(h, a), in this case GGX).

Don't forget however, that you also have a shadowing-masking function (G(l, v, a)) which has to be evaluated according to the chosen microfacet distribution. This means that I can't just mix-and-match f.e. GGX sampling and the Smith-Schlick evaluation; I have to use the Smith-GGX solution:

float G_Smith_GGX(float ndotl, float ndotv, float roughness)
{
    float a = roughness * roughness;
    float a2 = a * a;

    float lambda_l = sqrt(a2 + (1.0 - a2) * ndotl * ndotl);
    float lambda_v = sqrt(a2 + (1.0 - a2) * ndotv * ndotv);

    return (2.0 * ndotl * ndotv) / (ndotv * lambda_l + ndotl * lambda_v);
}


Which explains the UE4 approach, as:

vec3 F = F_Schlick(F0, ldoth);
float G = G_Smith_GGX(ndotl, ndotv, roughness);

// PDF = (D(h) * dot(n, h)) / (4 * dot(v, h))
return (F * G * vdoth) / (ndotv * ndoth + EPSILON);


See this article for the derivation.