Perturbative proof of unitarity of $S$-matrix in QED In any standard textbook on QFT I know it is claimed that the $S$-matrix in QED is a unitary operator. I have never seen any proof of it. This should be compared with the analogous property of $S$-matrix in the non-relativistic scattering theory when the unitarity of the $S$-matrix is a difficult theorem which is proved only under certain assumptions on the form of the interaction potential (e.g. the Coulomb potential is excluded from the theory). 

I would expect that there should be a way to check the unitarity of the $S$-matrix in QED perturbatively, in each order of the perturbation theory. Was it done somewhere? A reference would be very helpful.

Remark. In the known to me literature, the elements of $S$-matrix are always computed between two (initial and final) collections of free particles states. The unitarity of $S$-matrix in this form implies that a bunch of free particles cannot form after an interaction a collection of bound states or a single bound state (since the probability of this process is zero). Thus the unitarity of $S$-matrix has very physical consequences. 
 A: This is not a full answer but hopefully it is helpful to you.
The S-matrix is defined as a change of basis between the 'in' and 'out' states, both span the same Hilbert space. So, by definition, if an S-matrix exists, it must be unitary. It remains to show that the 'in' and 'out' states can be constructed for a given theory. Certain restrictions on the interaction are required for that construction to be possible. Formal construction of the 'in' and 'out' states are part of the LSZ formalism so you can read more on that. For QED, the problem is a bit more complicated because of the long-range interaction due to the massless photon.
$H_0$, $H$ and bound states
On your remarks about the spectra of $H_0$, $H$ and bound states, to quote Weinberg's words:

Note that $H_0$ is assumed here to have the same spectrum as the full
  Hamiltonian $H$. This requires that the masses appearing in $H_0$ be
  the physical masses that are actually measured, which are not
  necessarily the same as the 'bare' mass terms appearing in $H$; the
  difference if there is any must be included in the interaction $V$,
  not $H_0$. Also, any relevant bound states in the spectrum of $H$
  should be introduced into $H_0$ as if they were elementary particles.

My interpretation of this in the context of QED is that you need to put in the physical mass, renormalized fields to pick out the free Hamiltonian $H_0$. The rest (including counterterms) will be considered parts of the interaction. Now, what you include in $H_0$ depend on the type of scattering you wish to consider. 
If you want scattering between photons and electrons, then $H_0$ is just the textbook free QED Hamiltonian. You want protons scattering as well? No problem, just put in another Dirac field with the correct mass for the protons.
But if you need to consider scattering with a Hydrogen atom bound state, the above free Hamiltonian $H_0$ won't do. You will need to include new terms to allow for the free Hydrogen atom. If these is no interaction, there is certainly no bound state. The easiest way to do this, as Weinberg pointed out, is to treat the bound state as an elementary particle and include a new field and appropriate interaction for it. Protons are also composite particles but did you think twice about including a free Dirac field for it?
Sidney Coleman provided the following (non-relativistic) example to shred light on this topic. Consider three particles (same mass for simplicity) with central potentials between them:
$$
H=\sum_{i=1}^3\frac{p_i^2}{2m}+V_{12}(|\mathbf{x}_1-\mathbf{x}_2|)+V_{23}(|\mathbf{x}_2-\mathbf{x}_3|)+V_{13}(|\mathbf{x}_1-\mathbf{x}_3|)
$$
Here $V_{12}$ is assumed to be strong enough to make a bound state by itself.
Now, if you want to consider scattering between the three types of free particles, the 'in'/'out' states would be $|\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3\rangle^{\text{in,out}}$, and the free Hamiltonian is simply:
$$
H_0=\sum_{i=1}^3\frac{p_i^2}{2m}
$$
But if you want to consider $|\mathbf{p},\mathbf{p}_3\rangle^{\text{in,out}}$, the scattering between a 1-2 pair bound state with combined momentum $\mathbf{p}$ and a type-3 particle, the free Hamiltonian is different:
$$
H_0=\frac{p^2}{2\mu}+\frac{p^2_{cm}}{4m}+V_{12}(r)+\frac{p_3^2}{2m}
$$
Here, we need to include the $V_{12}(r)$ interaction to allow for the free 1-2 bound state. I can imagine it is possible to apply a similar approach to a field theory with a bound state but introducing a new field to treat each bound state as a new elementary particle is certainly much simpler.
