Diagonalizing a linearized optomechanical Hamiltonian in the paper Nonlinear Interaction Effects in a Strongly Driven Optomechanical Cavity, the authors diagonalize the Hamiltonian (equation (2) in the paper)
$$H_1=−Δd^†d+ω_Mb^†b+G(d+d^†)(b+b^†)$$
in the appendix (equations S1-S3). First they define the vector $\vec X =[b\,d\,b^†\,d^†]^T$, and then say it can be done "by standard means". I assumed they would mean writing the Hamiltonian as
$$H=\vec X^TM\vec X,\quad M=\left(\begin{matrix}\omega_M&0&0&0\\G&-\Delta&G&0\\0&0&0&0\\G&0&G&0\end{matrix}\right)$$
but M is not diagonalisable according to WolframAlpha. 
What am I missing here?
 A: Actually the form you are looking for is
$$
H = X^\dagger M X = \left[\begin{array}{cccc}b^\dagger & d^\dagger & b & d\end{array}\right]\left(\begin{array}{cccc}\alpha_1 & \beta_1 & \gamma_1 & \delta_1 \\ \alpha_2 & \beta_2 & \gamma_2 & \delta_2 \\ \alpha_3 & \beta_3 & \gamma_3 & \delta_3 \\ \alpha_4 & \beta_4 & \gamma_4 & \delta_4 \end{array}\right)\left[\begin{array}{c}b\\ d\\ b^\dagger\\d^\dagger\end{array}\right] 
$$
For the matrix $M$ try 
$$
M = \frac{1}{2}\left(\begin{array}{cccc}\omega_M & G & 0 & G\\ G & -\Delta & G & 0 \\ 0  & G & \omega_M &G \\ G & 0 & G & -\Delta \end{array}\right)
$$
Indeed, 
$$
H = \frac{1}{2}\left[\begin{array}{cccc}b^\dagger & d^\dagger & b & d\end{array}\right] \left(\begin{array}{cccc}\omega_M & G & 0 & G\\ G & -\Delta & G & 0 \\ 0  & G & \omega_M &G \\ G & 0 & G & -\Delta \end{array}\right)\left[\begin{array}{c}b\\ d\\ b^\dagger\\d^\dagger\end{array}\right] = \\\; \\
= \frac{1}{2}\left[\begin{array}{cccc}b^\dagger & d^\dagger & b & d\end{array}\right]\left[\begin{array}{c}\omega_M b + Gd + Gd^\dagger\\ Gb  - \Delta d + Gb^\dagger\\ Gd + \omega_Mb^\dagger +  Gd^\dagger\\Gb + Gb^\dagger - \Delta d^\dagger \end{array}\right] = \\\; \\
= \frac{1}{2}\left[ \omega_Mb^\dagger b + Gb^\dagger d + Gb^\dagger d^\dagger + Gd^\dagger b  - \Delta d^\dagger d + G d^\dagger b^\dagger +\\
Gbd + \omega_M bb^\dagger + Gbd^\dagger + Gdb + Gdb^\dagger - \Delta dd^\dagger + Gb^\dagger d \right] 
$$
Summing up and rearranging,
$$
H = \frac{\omega_M}{2}\left( b^\dagger b + bb^\dagger\right) - \frac{\Delta}{2}\left( d^\dagger d + dd^\dagger\right) + G\left( b^\dagger d + b^\dagger d^\dagger + bd + bd^\dagger \right)
$$
and finally,
$$
H = \frac{\omega_M - \Delta}{2} - \Delta d^\dagger d + \omega_M b^\dagger b + G\left(b + b^\dagger \right)\left(d + d^\dagger \right) 
$$
where I assumed the canonical commutation relations (bosons).
All looks good and matrix $M$ is symmetric and diagonalizable.
