Agreeing well with experiment is an understatement: the theoretical value matches the experimental value to 10 significant figures, making it the most accurate prediction in all of science.
The $g$-factor for elementary spin-$1/2$ particles is $g=2a+2$, where $a$ is (essentially defined as) the anomalous magnetic moment - the constant $2$ is the classical, tree-level result from the Dirac equation, which is modified by higher-order diagrams in the perturbative series.
By comparing to the non-relativistic approximation of the Dirac equation (i.e. the Hamiltonian for the Zeeman effect etc.) and finding the term that corresponds to the classical $\vec{B}\cdot\vec{S}$ interaction, you can show that contributions to the anomalous magnetic moment must come from terms in off-shell $\bar{u}_2\gamma^\mu u_1$ matrix elements that have the same structure as $\ \!\not D^2-D^2$, where $\ \!\not D$ and $D$ are the covariant derivatives from minimal coupling in spinor and scalar QED respectively. Hence the anomalous magnetic moment may be obtained perturbatively by considering diagrams of the form:

where the photon is strictly off-shell and the central blob indicates that you can throw in whatever loops you want (that are consistent with the theory, obviously). Eventually you'll find that the matrix element can be parameterised by two form factors:
$$
i\mathcal{M}^\mu=-ie\bar{u}_2\left(F_1\left(\frac{p^2}{m^2}\right)\gamma^\mu+\frac{i\sigma^{\mu\nu}}{2m}p_\nu F_2\left(\frac{p^2}{m^2}\right)\right)u_1
$$
[see e.g. Schwartz, Quantum Field Theory and the Standard Model, Eq. (17.10)] [1]
Finally, as only the $F_2$ term has the structure of $\ \!\not D^2-D^2$, we can identify $a$ with $F_2(0)$, where we have chosen an energy scale that matches experimental conditions. The $F_1$ term is purely electric, modifying the scale dependence of the electric charge, giving rise to effects like the Lamb shift (For a more explicit derivation of the form factors, see David Morgante's answer here).
From here, it's just a matter of combinatorially enumerating all the loop diagrams at successive orders and calculating them (more recent endeavours have used computer techniques to complete corrections at $\alpha^4$!). Here's an example of a 1-loop diagram that contributes to $a$ (this is actually the only contributing 1-loop diagram):

Dyson's argument tells us that the radius of convergence for QED perturbation series is zero, no matter how small the expansion parameter/coupling constant ($\alpha$, in this case), the series will never converge: the best we can do is an asymptotic but non-convergent series, typically over Feynman diagrams. So how can we expect to calculate anything perturbatively, even after renormalisation? How can we make any predictions using QED at all? Although asymptotic series are eventually infinite, you can truncate the series at a certain finite order and obtain a good approximation to the modeled function. The non-convergence is surprisingly a good thing: from modern beta-function arguments, the Landau pole arises as the QED coupling constant varies with energy, and the absolute convergence of the series would imply the instability of the vacuum state (see e.g. [3]) - so an asymptotic series is a very reasonable tradeoff. However, we will be unable to make predictions at finite energies by summing over all orders, since the series will diverge (and badly too, since eventually each successive term typically gets larger and larger).
Determining when to best truncate the series is, in general, a difficult problem known as optimal truncation - you might also try employing techniques like Borel resummation so as not to throw away the information encoded in the later terms (but this is thought to be not all that helpful for QED [3]). For QED and other physical field theories which utilise asymptotic series, however, reasonably justified arguments based on non-perturbative effects suggest terminating the series at $N = 1/g'$, where $g'$ is the coupling constant associated with the theory (at $p^2 = 0$, for our purposes here). As $\alpha(0)\approx\frac{1}{137}$, we can guess that $N=137$ would be a good place to stop - but given our current technological limitations, we might never hope to reach beyond $N=5$ or $6$. As the typical behaviour of asymptotic series is to "converge before they diverge", first- to fourth-order effects seem to contribute the lion's share of the true value of the magnetic moment, and hence the $g$-factor, which is why our predictions match so well with experiment.
Further Reading:
- [1]: Schwartz, Quantum Field Theory and the Standard Model
- [2]: arXiv:0806.2196
- [3]: G.Parisi, Phys. Lett. 76B (1978) 65