Firstly, I highly suggest reading Eric Heitz's paper "Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs", which covers the full derivation of microfacet-based BRDFs.
The $\frac{1}{4(N \cdot V)(N \cdot L)}$ term is a side effect of the derivation of the BRDF for specular microfacets. Specifically, it comes from the Jacobian of the reflection transformation. See the paper and/or Walter et. al's 2007 paper for more details.
As for your concern for a divide by zero, the definition of the rendering equation prevents it. Let me explain:
For $(N \cdot V)$ to equal zero, they must be orthogonal. In this case, the visibility term of the rendering equation will cull the cases at / below the horizon.