Basis change matrices from spatial low rank are not unitary

[kevinsung]

Tested this with LiH. Are they supposed to be unitary?

import cirq
import openfermion
from openfermion import low_rank_spatial_two_body_decomposition, prepare_one_body_squared_evolution

def load_molecular_hamiltonian(
        geometry,
        basis,
        multiplicity,
        description,
        n_active_electrons,
        n_active_orbitals):

    molecule = openfermion.MolecularData(
            geometry, basis, multiplicity, description=description)
    molecule.load()

    n_core_orbitals = (molecule.n_electrons - n_active_electrons) // 2
    occupied_indices = list(range(n_core_orbitals))
    active_indices = list(range(n_core_orbitals,
                                n_core_orbitals + n_active_orbitals))

    return molecule.get_molecular_hamiltonian(
            occupied_indices=occupied_indices,
            active_indices=active_indices)

# 4-qubit LiH 2-2 with bond length 1.45
bond_length = 1.45
geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., bond_length))]
hamiltonian = load_molecular_hamiltonian(
        geometry, 'sto-3g', 1, format(bond_length), 2, 2)

# Perform the low rank decomposition of two-body operator.
eigenvalues, one_body_squares, _, one_body_correction = (
    low_rank_spatial_two_body_decomposition(
        hamiltonian.two_body_tensor,
        spin_basis=True))

# Get basis transformation matrices and print whether they are unitary
for j in range(len(eigenvalues)):
    _, basis_change_matrix = (
        prepare_one_body_squared_evolution(one_body_squares[j]))
    print(cirq.is_unitary(basis_change_matrix))

Output:

False
False
True
False

[babbush]

Yes, this is a problem. We'll look into it.