Bogoliubov transformation for bosons (matrix calculation) I'd like to know if there is a general numerical method of diagonalizing the bosonic quadratic Hamiltonian below
$$H=\sum_{i,j=1}^NT_{ij}b_i^\dagger b_j+\frac{1}{2}\sum_{i,j=1}^N\left(U_{ij}b_i^\dagger b_j^\dagger+U_{ij}^*b_ib_j\right)\!,$$
using the Bogoliubov transformation, where $T$ is an $N\times N$ Hermitian matrix and $U$ is $N\times N$ symmetric but can in general be complex. One may rewrite $H$ in $2N\times 2N$ matrix form
$$H=\frac{1}{2}\begin{pmatrix}
b^\dagger & b
\end{pmatrix}\begin{bmatrix}
T & U\\
U^* & T^*
\end{bmatrix}\begin{pmatrix}
b\\ b^\dagger
\end{pmatrix}+\mathrm{const},$$
where $b=(b_1,b_2,\ldots,b_N)$ is the vector of bosonic annihilation operators. The main challenge is that the Bogoliubov transformation for bosons
$$\begin{pmatrix}
b\\b^\dagger
\end{pmatrix}=\begin{bmatrix}
M & N\\
N^* & M^*
\end{bmatrix}\begin{pmatrix}
\xi \\ \xi^\dagger
\end{pmatrix}$$
is not unitary, but simplectic, as has been discussed in many similar questions. To preserve the bosonic commutation relations of $[b,b^\dagger]$ and $[\xi,\xi^\dagger]$, the transformation matrix satisfies
$$\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}\begin{pmatrix}
M^\dagger & N^T\\
N^\dagger & M^T
\end{pmatrix}=\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}\!.$$
For fermions, $-I\,$ becomes $I$ and the transformation matrix is unitary. So is there a general numerical procedure to find the bosonic Bogoliubov transformation that diagonalizes $H$ into
$$H=\frac{1}{2}\begin{pmatrix}
\xi^\dagger & \xi
\end{pmatrix}\begin{bmatrix}
\Gamma & 0\\
0 & \Gamma
\end{bmatrix}\begin{pmatrix}
\xi \\ \xi^\dagger
\end{pmatrix}+\mathrm{const},$$
where $\Gamma$ is real-diagonal?
 A: Self answer:
We start with a simple identity
$$\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}\begin{pmatrix}
M^\dagger & N^T\\
N^\dagger & M^T
\end{pmatrix}\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}=\begin{pmatrix}
M^\dagger & -N^T\\
-N^\dagger & M^T
\end{pmatrix}\!.\quad(*)$$
Then the symplectic condition that the canonical transformation matrix satisfies can be rewritten as
$$\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}\begin{pmatrix}
M^\dagger & -N^T\\
-N^\dagger & M^T
\end{pmatrix}=\begin{bmatrix}
I & 0\\
0 & I
\end{bmatrix}\!.$$
Therefore, we see the transformation matrix and its inverse multiplied together to get the identity matrix. Plugging the transformation into the Hamiltonian, we obtain
$$\begin{pmatrix}
M^\dagger & N^T\\
N^\dagger & M^T
\end{pmatrix}\begin{bmatrix}
T & U\\
U^* & T^*
\end{bmatrix}\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}=\begin{bmatrix}
\Gamma & 0\\
0 & \Gamma
\end{bmatrix}\!.$$
We now use again the equation $(*)$ to obtain
$$\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}\begin{pmatrix}
M^\dagger & -N^T\\
-N^\dagger & M^T
\end{pmatrix}\begin{bmatrix}
I & 0\\
0 & -I
\end{bmatrix}\begin{bmatrix}
T & U\\
U^* & T^*
\end{bmatrix}\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}=\begin{bmatrix}
\Gamma & 0\\
0 & \Gamma
\end{bmatrix}\!,$$
which then simplifies to
$$\begin{pmatrix}
M^\dagger & -N^T\\
-N^\dagger & M^T
\end{pmatrix}\begin{bmatrix}
T & U\\
-U^* & -T^*
\end{bmatrix}\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}=\begin{bmatrix}
\Gamma & 0\\
0 & -\Gamma
\end{bmatrix}\!,$$
and hence we have
$$\begin{bmatrix}
T & U\\
-U^* & -T^*
\end{bmatrix}=\begin{pmatrix}
M & N\\
N^* & M^*
\end{pmatrix}\begin{bmatrix}
\Gamma & 0\\
0 & -\Gamma
\end{bmatrix}\begin{pmatrix}
M^\dagger & -N^T\\
-N^\dagger & M^T
\end{pmatrix}\!.$$
By diagonalizing the non-Hermitian matrix on the left-hand side (which equals $\mathrm{diag}(I,-I)$ times the Hamiltonian), we obtain the transformation matrix whose columns are the eigenvectors of the non-Hermitian matrix.
A: I would suggest the following reference: A. G. D. Maestro   and M. J. Gingras, J. Phys. Condens. Matter 16, 3339 (2004). It has an appendix that fully describes how to do those steps (in k-space). 
Please note that it is not enough to diagonalize the matrix [T, U;-U*,-T*], the eigenvectors of this matrix are not necessarily eigenvectors of [T, U;U*,T*] as explained in the reference as you the Bogoliubov matrix must satisfy three conditions at the same time. Letting H be the square matrix that you want to diagonalize, Z is the Bogoliubov transformation, and G=diag[I,-I], the conditions that Z needs to satisfy are
(1) $Z^{\dagger}HZ=\rm diag[e1,..,en,e1,..en]$ 
where ek is the energy dispersion depending on the system in hands. For example, in rare earth pyrochlores you have 4 bands and thus n=4. 
(2) $ZGZ^{\dagger}=G$; this is basically the commutator relations in compact form [i.e. you want to impose that your new operators are bosonic as well]. 
(3) $Z^{-1} GM Z=-\rm diag[e1,...,en,-e1,..,-en]$ 
The requirement in (3) is actually comes from combining (1,2). Surprisingly, solving (3) is not always enough! Considering the energy dispersions for bosons you can get them from (3). For the Bogoliubov matrix you need first to construct a matrix $\tilde{Z}$ from the eigenvectors of GM, then impose that 
\begin{equation}
Z=\tilde{Z}P
\end{equation}
where P is block diagonal matrix. Using this expression in (2), we find that $PGP^{\dagger}=W$, where $W=(\tilde{Z}^{\dagger}G\tilde{Z})^{-1}$ is block diagonal. The procedure is to evaluate W using the result obtained earlier for $\tilde{Z}$, then using the formula $PGP^{\dagger}=W$, we conclude that the ith block in P is related to the ith block in W as $\pm P_{i}P_{i}^{\dagger}=W_{i}$. Now, since W is Hermitian by construction (check the definition of W above), then we can use linear algebra to say that the block $W_{i}$ is Hermitian and diagonalizable, i.e. there exist an invertible matrix $X_{i}$ such that $W_{i}=X_{i}D_{i}X_{i}^{-1}$, where $D_{i}$ is a diagonal matrix with eigenvalues of $W_{i}$ on the diagonal. Using this expression of $W_{i}$ together with $\pm P_{i}P_{i}^{\dagger}=W_{i}$, one can easily verify that 
\begin{equation}
P_{i}=X_{i}\sqrt{\pm D_{i}}X_{i}^{-1}
\end{equation} 
This way you construct the matrix P and thus the Bogoliubov matrix becomes 
\begin{equation}
Z=\tilde{Z}P
\end{equation}
This is basically the summary of the numerical (could be sometimes exact depending on the problem in hand) method for diagonalizing quadratic bosonic Hamiltonian in $\vec{k}$-space. 
