Why can we add counterterms? I'm having a hard time understanding why renormalized perturbation theory works. Why is it permissible to add counterterms to the Lagrangian? Terms which are often divergent themselves and carry unphysical parameters. How is the Lagrangian still well-defined after renormalization?
 A: We can never "add counterterms" to a Lagrangian. The Lagrangian is just the Lagrangian; if you change it by adding terms, then you've got a different theory. There's a similar issue when textbooks talk about "adding" a term for "gauge fixing", which is equally nonsensical. Such language only exists because of historical confusion.
It all much sense in the framework of Wilsonian renormalization, in which every theory has a cutoff scale. Suppose a theory is described by a Lagrangian $\mathcal{L}_{\text{fund}}$ at a high cutoff scale $\Lambda$. Then it is equivalent to some theory with an effective Lagrangian $\mathcal{L}_{\text{eff}}$ at low cutoff scale $\mu$, and we define
$$\mathcal{L}_{\text{fund}} = \mathcal{L}_{\text{eff}} + \Delta \mathcal{L}.$$
Now, since the cutoff scale $\mu$ is low, quantum corrections in the effective theory are not as important. In particular, the full results in the effective theory can be roughly extracted by working at tree level.
This is the modern approach to the subject, but historically the logic went almost completely backwards. First, people computed only at tree level, and found they could fit physical observations using $\mathcal{L}_{\text{eff}}$. Then they tried to compute loop corrections with an infinite cutoff. If one does this with $\mathcal{L}_{\text{eff}}$, then one gets completely incorrect results, so textbooks speak of "adding the counterterm" $\Delta \mathcal{L}$. This is conceptually wrong because, if we wanted to work with a high cutoff, then we should have included $\Delta \mathcal{L}$ in the Lagrangian from the start. The way most textbooks phrase it, it looks like an ad hoc correction pasted on halfway through a calculation, but it really needs to be there the whole time. (For more detail about how the old textbook approach fits onto the modern approach, see here.)
Incidentally, the infinities are another red herring. In the modern view, we know that all cutoffs should be finite, in which case all predictions and parameters are finite. The trouble only occurs when one tries to take $\Lambda \to \infty$, which we have no reason to do. For example, textbooks will often tell you the QED coupling diverges at high energies, but it remains perturbative for any finite cutoff below the Planck scale.
A: The counterterms are a prickly subject, I recently got used to it, this will be a good opportunity to explain it if I've understood it well. 
Actually the counterterms are not added to the Lagrangian, they arise by splitting up the Lagrangian in its renormalized part and the rest , the counterterms. 
Upon the computation of loop diagrams in perturbation theory many -- not all -- turn out to be infinite. The physicists understood in the meantime (from the first confrontation to this problem until today) that the important parameters (typically mass and coupling of the particle to its interactions) of a quantum field theory (QFT) which appear in the Lagrangian are not the same than the parameters they actually measure. This is actually not so difficult to understand, you probably heard that the electron when exposed to electromagnetic (EM) interaction changes its charge because it is surrounded by a cloud of virtual electron-positron pairs. So if measured from a large distance the charge of the electron is shielded by this cloud and is not the real electron charge which is usually called the bare charge of the electron or particle. The same happens with the mass of the electron: if the self-interaction of the electron's field is taken into account, its mass changes. 
Everthing would be rather easy to understand if these parameter shifts due to interaction were finite, but unfortunately they are not. This observation is certainly related with the fact that we consider electrons as point-like particles. Doing the contrary usually also leads to contradictions, so in usual QFT electrons and other particles are considered as point-like. 
As computations with infinities are senseless, the first step is to carry out a regularisation to make the corresponding expressions finite which allows us to keep on with the calculations. The final goal is of course to get rid of the infinite or potentially infinite == regularized expressions. 
So in renormalized perturbation theory the following strategy is pursued: As we physicists know that infinite expressions appear in the loop computations, their appearence is already anticipated in the establishment of the Lagrangian. The bare parameters of the corresponding Lagrangian are assumed to be infinite and then split up in 2 parts: its renormalized part and the singular part, the latter corresponds to the counterterms. 
Then the not yet fixed counterterms are considered as part of the interactive Lagrangian and included into the computation of the loop diagrams and finally adapted (i.e. fixed) in order to compensate the infinite expressions in the loop diagram evaluation. This is done meticulously order by order in perturbation theoy. In this way the resulting expression become finite and independent of the chosen regularisation. 
I think, it is actually more intuitive to use bare perturbation theory where the anticipation of the splitting of the Lagrangian is not done and only at the very end of the computation after the regularisation the renormalized parameters are set in the expressions (with the final consequence that the dependence of the renormalized parameters and the bare parameters shows that by doing so the bare parameter will be infinite once the finite cutoff used for the regularisation is taken to infinite.).
But this latter procedure is technically more complicated in particular in perturbation theory of higher order and above all has the disadvantage that the perturbation series is developed in the bare coupling parameter which is supposed to be small but is actually not small at all. But renormalized perturbation theory where the development of the series is done in the renormalized coupling parameter has not this problem. 
This post only sketches the strategy applied in regularisation and renormalisation in QFT, many details should be still said, I limit myself to the following:
One has to be used to the fact that the bare Lagrangian contains formally many "infinite" (after the regularisation cutoff is taken to infinite) expressions, however all the symmetries a Lagrangian has to comply are usually preserved, so the Lagrangian has still all the properties which are required from a Lagrangian. And one has to bear in mind that a Lagrangian and "its contents" cannot be measured. This is actually what is now called "effective field theory", it is not considered to be a fundamental theory, a fundamental theory should be completely finite, so not exposed to infinities. This fundamental theory is still to be found, and many physicists are hard working on it. That also explains the interest of the physics community in string theory which is a possible answer to this problem .. but there are still many open questions. 
A: This point used to confuse me a lot. I think the following toy example may clear up much confusion. 
Imagine a theory with only one coupling constant which I will call $g$. We think Wilsonian, where our theory has a finite cutoff at a scale $\Lambda$ way beyond whatever scale of interest we are discussing. We know that if we integrate out modes progressively, one $d\Lambda$ slice at a time, we generate a smooth flow in the space of couplings. Suppose the one-loop beta function is $$\Lambda \frac{dg}{d\Lambda} = - \frac{g^2}{2\pi}$$ Again, in principle this can be derived totally from the Wilsonian picture where we integrate out thin shells in momentum space and work iteratively: no infinities in, no infinities out. This describes an asymptotically free theory for which perturbation theory becomes uniformly valid in the high energy regime. We can integrate this equation to find $$g(\Lambda) = \frac{g(\mu)}{1 + \frac{g(\mu)}{2\pi} \ln(\Lambda/\mu)}$$ As we take $\Lambda \to \infty$, we have $g \to 0$, a totally well defined limit. Here, $g(\Lambda)$ is the bare coupling and $g(\mu)$ is the effective (renormalized) coupling at some scale $\mu$. This describes a totally smooth (and logarithmically slow) flow down to whatever the scale $\mu$ is. The numerical difference between $g(\Lambda)$ and $g(\mu)$ will not be too serious (that is, it will not be infinite) provided we stay within the regime of perturbation theory, that is, don't go to too low energy. 
However, if we insist on doing everything in a power series in the renormalized coupling $g(\mu)$ (and this is the situation we consider in counterterms, where we expand the bare coupling as renormalized coupling plus perturbative corrections), we find $$g(\Lambda) = g(\mu) - \frac{g^2(\mu)}{2\pi} \ln(\Lambda/\mu) + \frac{g^3(\mu)}{4\pi^2} (\ln(\Lambda/\mu))^2 + ...$$
Even though the coupling has a totally regular $\Lambda \to \infty$ limit, there are "ultraviolet divergences" at each order of the expansion. The $g^2$ term neutralizes the one loop singularity which generated the flow in the first place, and the higher order terms are the leading logs from higher orders. The "divergences" in perturbation theory just come from a long RG flow, but disappear once we use the effective couplings, which is the general statement of universality (renormalizability). No couplings are truly singular (unless something drastic like a Landau pole happens), and renormalization is in this sense a small correction: after the bookkeeping, it just gives small logarithmic corrections to scaling. What makes it look like wizardry is using the perturbative technique without explaining the Wilsonian origin, but usually the perturbative technique is easier to calculate with. 
Hope this helped.
