Watch out: Your problem arises because this second-quantized Hamiltonian contains terms that do not conserve the particle number (e.g. $d b$ in the "G" term). As a result, when "diagonalizing" such a Hamiltonian, you actually need a Bogoliubov transformation (instead of a simple unitary transformation) to obtain new normal modes. The requirement for this transformation, which mixes creation and annihilation operators, is that the new operators $d'$ and $b'$ fulfill the same bosonic commutation relations as the old ones (e.g. $[d',d'^{\dagger}]=1$).

While it is formally possible to write everything in terms of a matrix $M$ (as you have done), finding the Bogoliubov transformation is not equivalent to diagonalizing this matrix in the usual way. You need to fulfill the condition $[X_j,X_k^{\dagger}]=\delta_{jk} \sigma_k$, where $\sigma_k=+1$ for the $k=1,2$ and $\sigma_k=-1$ for $k=3,4$ in your example (the minus sign comes about because $X_3^{\dagger}=X_1$ and so on). Setting your Bogoliubov transformation to be $X_j=S_{jn} Y_n$ (summation over repeated indices), this implies $S \Sigma S^{\dagger} = \Sigma$, where $\Sigma$ would be the diagonal matrix with entries given by $\sigma_k$ ($1,1,-1,-1$). This is different from the usual unitary transformation, and indeed is known as a symplectic transformation.