Renormalization Group approaches to many-body problems: An algebraic perspective

This is related to a question I asked a few months ago here. I haven’t looked at it again until recently because I was busy with something very different, but now I am coming back at it and wanted to ask something related with greater detail.

The method

Let me look at a general (Numerical) Renormalization Group sort of approach in a purely algebraic form: Let us say we have two vector spaces, $E$ and $F$, of dimensions $N$ and $M$ respectively, and let us denote by ${|e_i rangle }_{i=1}^{N}$ and ${ | f_j rangle }_{i=1}^{M}$ basis vectors for $E$ and $F$ respectively. Each space an be thought as the Hilbert space of the possible states of a system. For example, if $E$ is the space of states of an electron in a single site, $E = { |{ uparrow} rangle, |{ downarrow} rangle } $. Say for $E$ and $F$ we have linear operators
$$A : E rightarrow E $$
$$B : F rightarrow F $$
that we now how to diagonalize. These operators are of course the hamiltonians of the "blocks" defining vector spaces $E$ and $F$.

When we connect the two spaces $E$ and $F$, the linear operator $D : E otimes F rightarrow E otimes F $ describing the hamiltonian will look like
$$D = A otimes 1_{F}+ 1_{E}otimes B + underbrace{sum_{i i^prime j j^prime} C_{i i^prime j j^prime} |e_i otimes f_j rangle langle e_{i^prime}otimes f_{j^prime} |}_{C},$$
which is an $N times M$ matrix (I am using $|e_i otimes f_j rangle$ as an abbreviation for $|e_i rangle otimes | f_j rangle$). $E otimes F$ is an abbreviation for the vector space $ { sum_{i,j} a_{i j} | e_i otimes f_j rangle }.$
The way I see it, when implementing numerically a RG strategy, what one does abstractly boils down to this: One takes an orthonormal set of vectors ${| v_alpha rangle }_{alpha=1}^n$ in $E$ and another orthogonal set of vectors ${| u_beta rangle }_{beta=1}^m$ in $F$, and call $tilde{E} = text{Span}left( { v_alpha }_{alpha=1}^n right)$ and $tilde{F} = text{Span}left( { u_beta }_{beta=1}^m right)$ (In Wilson’s attack of the Kondo problem, I think, these would be some low-lying eigenstates of the hamiltonian; in White’s DMRG solution of the spin chains, these are eigenvectors of a density matrix extracted from the ground state of a larger system). Now, in terms of these sets we can construct approximations for $A$ and $B$ like
$$tilde{A} := text{Pr}_{tilde{E}} := sum_{alpha alpha^prime} langle v_alpha | A | v_{alpha^prime} rangle |v_alpha rangle langle v_{alpha^prime} |$$ and
$$tilde{B} := text{Pr}_{tilde{F}} := sum_{beta beta^prime} langle u_beta | B | u_{beta^prime} rangle |u_beta rangle langle u_{beta^prime} |,$$
($text{Pr}_X$ means projection over set $X$) and by expanding
$$| v_alpha rangle = sum_i v_{alpha i} | e_{i} rangle,$$
$$| u_beta rangle = sum_j u_{beta j} | f_{j} rangle$$
and calling
$$tilde{C}_{alpha alpha^prime beta beta^prime} = sum_{i i^prime j j^prime} v^ast_{alpha i}u^ast_{beta j} C_{i i^prime j j^prime} v_{alpha^prime i^prime}u_{beta^prime j^prime} $$
it is clear that we can write an approximation for $D$ which looks like
$$tilde{D} = tilde{A}otimes 1_{tilde{F}} + 1_{tilde{E}}otimes tilde{B} + underbrace{sum_{alpha alpha^prime beta beta^prime} tilde{C}_{alpha alpha^prime beta beta^prime} | v_alpha otimes u_beta rangle langle v_{alpha^prime} otimes u_{beta^prime}|}_{tilde{C}},$$
which is now only an $n times m$ matrix.
The process is to be iterated in some way, of course, to increase the size. For example, $tilde{E} otimes tilde{F}$ will be renamed $E$ for the next step. Maybe (as in the standard DMRG), $F$ will be always the state space of a single site. In this case the inputs for each new step are prepared as:
$$E rightarrow tilde{E} otimes tilde{F} $$
$$F rightarrow F$$
$$A rightarrow tilde{D} $$
$$B rightarrow B.$$

The Question(s)

Now well, the DMRG solution to spin chains is relatively easy to understand from this point of view, but I am blocked when I try to use this to diagonalize, for instance, very large Hubbard chains. The problem gets even harder if there is hopping to second, third or higher-order neighbors. The problem is how to keep track of how the connecting terms between the block (contained originally in $C$) change under repeated applications of the projections to the reduced sets. The first iteration is authomatic, but then, starting with the second one, it gets extraordinarily messy to rewrite the connections in the next reduced basis set (in the $n-$th step it will involve projecting $n$ times the hopping terms that one has defined in the full many-body basis). The second question is how to deal with the fluctuations in the number of electrons between the blocks. The way I have written it, it is clear that in principle each matrix is a block matrix with the structure of a direct sum corresponding to different number of electrons. It must be like this so that the matrices $D$ and its approximation $tilde{D}$ accounts for all possible fluctuations of charge between the blocks being put together, but after a few iterations this makes the matrices become too large (when probably one is only interested in building approximate matrices for the subspace in which the number of electrons in the "final" big system of real interest is fixed, e.g. N electrons in a lattice of N sites, etc)

How are the technicalities dealt with in real calculations? I have extensively looked for detailed accounts of the algorithms employed, but I have not encounter any full explanation except for spin Heisenberg-like systems (where DMRG was originally applied), where the connectivity terms are particularly simple and there is no fluctuation of particles between the blocks, so the scheme above works just fine.