The assumption underlying such model is the same as lots of other models for skin rendering; the subsurface scattering can be approximated as a diffusion phenomenon. This is good because in highly scattering media, the distribution of light loses dependency from the angle and tends to isotropy.
The dipole approximation is a formulation for the resolution of such diffusion problem in an analytical fashion.
Basically they start by approximating the BSSRDF as a multiple scattering and single scattering component. The multiple scattering is then defined as:

Where $F_t$ are Fresnel terms and $R$ is the diffusion profile expressed as function of the distance between the entry and exit point.
This $R$ is referred to as diffusion profile and they formulate this profile via a dipole approximation. The contribution of the incoming light ray is considered to be the one of two virtual sources: one negative beneath the surface and one positive above it (that's why dipole)

Here in the picture r is the $\|x_i - x_o\|$ above. The contribution of those light sources is dependent on various factors such as the distance of the light from the surface, scattering coefficient etc. (See below for a more detailed description of the formula itself).
This model only account for multiple scattering events, but that's fine enough for skin. It must be noticed though that for some translucent materials (e.g. smoke and marble) the single scattering is fundamental. That paper propose a single scattering formulation, but is expensive.
The diffusion profile is usually approximated for real-time application as a series of gaussian blurs (like in the seminal works of D'Eon et al. in GPU Gems 3 then used for the Jimenez's SSSSS) so to make it practical for real time scenarios. In this wonderful paper there are details on such approximation.
A picture from that paper show actually how good is this formulation:

As side note the dipole approximation assumes that the material is semi-infinite, however this assumption doesn’t hold with thin slabs and multi-layered material such as the skin. Building on the dipole work, Donner and Jensen [2005] proposed the multi-pole approximation that accounts for the dipole problems.
With this model instead of a single dipole, the authors use a set of them to describe the scattering phenomenon. In such formulation the reflectance and transmittance profiles can be obtained by summing up the contribution of the different dipoles involved
EDIT:
I am putting here the answers to a couple of @NathanReed 's questions in the comment section:
Even with the diffusion profile approximation, the BSSRDF model still requires integrating over a radius of nearby points on the surface to gather incoming light, correct? How is that accomplished in, say, a path tracer? Do you have to build some data structure so you can sample points on the surface nearby a given point?
The BSSRDF approximation still need to be integrated over a certain area, yes.
In the paper linked they used a Montecarlo ray-tracer randomly sampling around a point with a density defined as:
$\sigma_{tr}e^{-\sigma_{tr}d}$
Where that sigma value is the effective extinction coefficient defined below ( it is dependent on scattering and absorbing coefficient, which are properties of the material) and d is the distance to the point we are evaluating. This density is this way defined because the diffusion term has exponential fall-off.
In [Jensen and Buhler 2002] they proposed an acceleration technique. One of the main concepts was to decouple the sampling from the evaluation of the diffusion term. This way they perform a hierarchical evaluation of the information computed during the sampling phase to cluster together distant samples when it comes to evaluating the diffusion. The implementation described in the paper uses an octree as structure.
This technique, according to the paper, is order of magnitude faster than the full Monte Carlo integration.
Unfortunately I never got myself into an off-line implementation, so I can't help more than this.
In the real-time sum-of-Gaussians approximations the correct radius is implicitly set when defining the variance of the Gaussian blurs that need to be applied.
Why one positive and one negative light? Is the goal for them to cancel each other in some way?
Yes, the dipole source method (which dates way before Jensen's paper) is such defined to satisfy some boundary condition. Specifically the fluence must be zero at a certain extrapolated boundary that has a distance from the surface of $2AD$ where

Being $F_{dr}$ the fresnel reflectivity of the slab considered and that sigma value is the reduced extinction coefficient described below.
EDIT2: I have expanded (a tiny bit) some of the concepts in this answer in a blog post: http://bit.ly/1Q82rqT
For those who are not scared by lots of greek letters in a formula, here's an extract from my thesis where the reflectance profile is briefly described in each term:
