Bug in QuadraticHamiltonian diagonalization

[kevinsung]

The following code raises an error:

import numpy as np
import openfermion as of

hermitian_part = np.array(
    [[0.0, 1.0, 0.0],
     [1.0, 0.0, 1.0],
     [0.0, 1.0, 0.0]],
    dtype=complex
)
antisymmetric_part = np.array(
    [[0.0, 1.0j, 0.0],
     [-1.0j, 0.0, 1.0j],
     [0.0, -1.0j, 0.0]],
    dtype=complex
)

hamiltonian_quad = of.QuadraticHamiltonian(hermitian_part, antisymmetric_part)
hamiltonian = of.get_fermion_operator(hamiltonian_quad)
hamiltonian_jw = of.get_sparse_operator(of.jordan_wigner(hamiltonian))

orbital_energies, transformation_matrix, constant = hamiltonian_quad.diagonalizing_bogoliubov_transform()
occupied_orbitals = (1,)
eig = np.sum(orbital_energies[list(occupied_orbitals)]) + constant

circuit = of.prepare_gaussian_state(
    cirq.LineQubit.range(3),
    hamiltonian_quad,
    occupied_orbitals=occupied_orbitals)
state = cirq.final_state_vector(circuit, dtype=np.complex128)

np.testing.assert_allclose(hamiltonian_jw @ state, eig * state, atol=1e-8)

This is due to an incorrect assumption in the QuadraticHamiltonian code regarding the format of the Schur decomposition.