Fix bogoliubov_transform for spin-preserving Gaussian transformations (#1360)

Fixes #776.

`bogoliubov_transform` checks whether the transformation matrix is
spin-block-diagonal before it ever looks at the matrix shape, and that
check only makes sense for a square N x N matrix. For the N x 2N matrix
you get out of a general (non-particle-conserving) Bogoliubov
transformation it ends up inspecting the wrong blocks, decides that a
transform which doesn't actually mix spin is "block diagonal", and then
slices it into pieces of the wrong shape. The matrix from the issue dies
with:

```
ValueError: Bad shape for transformation_matrix. Expected (1, 1) or (1, 2) but got (1, 3).
```

### What I changed

I moved the block detection into a helper, `_spin_blocks`, that handles
both shapes. For the N x 2N case, write the matrix as `W = (W1 W2)` (it
acts on all the creation operators followed by all the annihilation
operators). The transform leaves the two spin sectors alone exactly when
both `W1` and `W2` are block diagonal, and in that case each sector is
its own (N/2) x N Gaussian transform built from the matching diagonal
blocks of `W1` and `W2`. The recursion that was already there handles
the rest. As a bonus, for spin-symmetric systems (BCS-style pairing
Hamiltonians, say) you now do two half-size Givens decompositions
instead of one full-size one.

### The part that took some digging

I didn't want to call this correct based on eigenstate energies alone,
so I checked the circuit at the operator level: does it actually give `U
a_p† U⁻¹ = b_p†` for every mode? That caught something the energy checks
can't.

Under Jordan-Wigner a spin-down ladder operator drags a Z string across
all of the spin-up qubits, so it anticommutes with the spin-up parity
operator. If the spin-up circuit flips spin-up parity (which happens
whenever its Gaussian decomposition uses an odd number of particle-hole
transforms), then splitting the sectors apart flips the sign of every
spin-down operator. What makes it sneaky is that it's invisible to
energy tests: starting from a fixed occupation state, that bad sign is
just a global phase, so the state still comes out as the correct
eigenstate with the correct energy even though the underlying unitary is
wrong.

So when the spin-up circuit flips parity I add a layer of Z gates on the
spin-down qubits to put the signs back (`_preserves_parity` decides
whether that's needed). With a fixed computational-basis initial state
it's only a global phase, so I skip it there.

I also ran into the general Gaussian decomposition not reliably
diagonalizing exactly-block-structured matrices (`U† H U` comes out
non-diagonal for some perfectly valid block inputs), so I went with
always splitting these instead of trying to detect when it's safe to
fall back. Splitting is the more correct path here, not just the faster
one.

### Tests

- `test_bogoliubov_transform_spin_block_gaussian_regression`: the actual
matrix from #776, normalized so it satisfies `W1 W1† + W2 W2† = I`. This
used to raise ValueError.
- `test_bogoliubov_transform_spin_block_operator_algebra`: checks `U
a_p† U⁻¹ = b_p†` on the full unitary, for both a parity-preserving and a
parity-flipping spin-up sector. This is the one that fails without the
sign fix.
- `test_bogoliubov_transform_spin_block_gaussian`: eigenstate energies
through the state-prep path.
- `test_bogoliubov_transform_spin_mixing_gaussian_not_split`: a
transform that genuinely mixes spin should still go through the general
path rather than getting split.

Outside the committed tests I ran a wider sweep (k = 2, 3, 4, real and
complex, a range of seeds) cross-checking three things that each catch a
sign error on their own: the operator algebra, whether `U† H U` is
diagonal, and the eigenstate energies. Everything passes, including 27
of the parity-flipping cases. The full `circuits/` suite passes and
black/pylint are clean.

Author: blazingphoenix7

src/openfermion/circuits/primitives/bogoliubov_transform.py

@@ -92,13 +92,14 @@ def bogoliubov_transform(
         initial_state = _occupied_orbitals(initial_state, n_qubits)
     initially_occupied_orbitals = cast(Optional[Sequence[int]], initial_state)
 
-    # If the transformation matrix is block diagonal with two blocks,
-    # do each block separately
-    if _is_spin_block_diagonal(transformation_matrix):
-        up_block = transformation_matrix[: n_qubits // 2, : n_qubits // 2]
+    # If the transformation matrix does not mix the two spin sectors, do each
+    # sector separately. This is both faster (two half-size decompositions) and
+    # more robust: the general fermionic Gaussian decomposition is numerically
+    # unreliable on exactly block-structured matrices. See issue #776.
+    spin_blocks = _spin_blocks(transformation_matrix)
+    if spin_blocks is not None:
+        up_block, down_block = spin_blocks
         up_qubits = qubits[: n_qubits // 2]
-
-        down_block = transformation_matrix[n_qubits // 2 :, n_qubits // 2 :]
         down_qubits = qubits[n_qubits // 2 :]
 
         if initially_occupied_orbitals is None:
@@ -110,6 +111,17 @@ def bogoliubov_transform(
                 i - n_qubits // 2 for i in initially_occupied_orbitals if i >= n_qubits // 2
             ]
 
+        # Spin-down ladder operators carry a Jordan-Wigner Z string over all
+        # spin-up qubits, so they anticommute with the spin-up parity operator.
+        # If the spin-up circuit flips spin-up parity (an odd number of
+        # particle-hole transformations) the factorized circuit would realize
+        # the spin-down operators with a flipped sign. A layer of Z gates on the
+        # spin-down qubits restores the documented action U a_p^dagger U^-1 =
+        # b_p^dagger. When an initial computational basis state is fixed this
+        # only contributes an irrelevant global phase, so it is skipped there.
+        if initially_occupied_orbitals is None and not _preserves_parity(up_block):
+            yield (cirq.Z(q) for q in down_qubits)
+
         yield bogoliubov_transform(up_qubits, up_block, initial_state=up_orbitals)
         yield bogoliubov_transform(down_qubits, down_block, initial_state=down_orbitals)
         return
@@ -122,6 +134,38 @@ def bogoliubov_transform(
         yield _gaussian_basis_change(qubits, transformation_matrix, initially_occupied_orbitals)
 
 
+def _spin_blocks(matrix: numpy.ndarray) -> Optional[Tuple[numpy.ndarray, numpy.ndarray]]:
+    r"""Split a transformation matrix into independent spin-sector transforms.
+
+    Returns a pair of transformation matrices (spin-up block, spin-down block)
+    if the transformation does not mix the two spin sectors, assuming modes are
+    ordered with all spin-up modes before all spin-down modes. Returns None if
+    the transformation mixes the sectors.
+
+    For an $N \times N$ matrix the blocks are the diagonal blocks. An
+    $N \times 2N$ matrix $W = (W_1 \; W_2)$ acts on the column vector of all
+    creation operators followed by all annihilation operators, so it preserves
+    the spin sectors if and only if both $W_1$ and $W_2$ are block diagonal;
+    the sector transforms are then $(N/2) \times N$ matrices formed from the
+    corresponding diagonal blocks of $W_1$ and $W_2$.
+    """
+    n, m = matrix.shape
+    if n % 2:
+        return None
+    k = n // 2
+    if m == n:
+        if _is_spin_block_diagonal(matrix):
+            return matrix[:k, :k], matrix[k:, k:]
+        return None
+    left_block = matrix[:, :n]
+    right_block = matrix[:, n:]
+    if _is_spin_block_diagonal(left_block) and _is_spin_block_diagonal(right_block):
+        up_block = numpy.concatenate([left_block[:k, :k], right_block[:k, :k]], axis=1)
+        down_block = numpy.concatenate([left_block[k:, k:], right_block[k:, k:]], axis=1)
+        return up_block, down_block
+    return None
+
+
 def _is_spin_block_diagonal(matrix) -> bool:
     n = matrix.shape[0]
     if n % 2:
@@ -131,6 +175,30 @@ def _is_spin_block_diagonal(matrix) -> bool:
     return numpy.isclose(max_upper_right, 0.0) and numpy.isclose(max_lower_left, 0.0)
 
 
+def _preserves_parity(transformation_matrix: numpy.ndarray) -> bool:
+    r"""Whether the Bogoliubov circuit for this block preserves fermion parity.
+
+    A particle-number-conserving (square) transformation uses only Givens
+    rotations and always preserves parity. A general $N \times 2N$ Gaussian
+    transformation preserves parity if and only if the orthogonal matrix that
+    implements it on the Majorana operators lies in the identity component of
+    O(2N), i.e. has determinant $+1$; the other component requires an odd
+    number of particle-hole transformations, each an X gate that flips the
+    parity. That determinant equals the determinant of the unitary
+    Bogoliubov-de Gennes extension of $W = (W_1 \; W_2)$, namely
+    $\begin{pmatrix} W_1 & W_2 \\ W_2^* & W_1^* \end{pmatrix}$, which is real and
+    equal to $\pm 1$, so the parity is read off directly without running the
+    full fermionic Gaussian decomposition.
+    """
+    n, m = transformation_matrix.shape
+    if m == n:
+        return True
+    w1 = transformation_matrix[:, :n]
+    w2 = transformation_matrix[:, n:]
+    bogoliubov_de_gennes = numpy.block([[w1, w2], [numpy.conjugate(w2), numpy.conjugate(w1)]])
+    return numpy.isclose(numpy.linalg.det(bogoliubov_de_gennes), 1.0)
+
+
 def _occupied_orbitals(computational_basis_state: int, n_qubits) -> List[int]:
     """Indices of ones in the binary expansion of an integer in big endian
     order. e.g. 010110 -> [1, 3, 4]"""

src/openfermion/circuits/primitives/bogoliubov_transform_test.py

@@ -116,6 +116,182 @@ def test_spin_symmetric_bogoliubov_transform(
     numpy.testing.assert_allclose(quad_ham_sparse.dot(state), energy * state, atol=atol)
 
 
+def _independent_spin_sector_transform(n_spatial_orbitals, seed_up, seed_down, real):
+    """Assemble an N x 2N Bogoliubov transformation that does not mix spin.
+
+    Builds two independent non-particle-conserving quadratic Hamiltonians, one
+    per spin sector, and places their diagonalizing transforms into a combined
+    matrix. The matrix acts on the column vector of all creation operators
+    followed by all annihilation operators, so each sector's creation
+    (annihilation) block lands in the left (right) half, offset by the sector's
+    position within that half.
+    """
+    quad_ham_up = random_quadratic_hamiltonian(
+        n_spatial_orbitals, conserves_particle_number=False, real=real, seed=seed_up
+    )
+    quad_ham_down = random_quadratic_hamiltonian(
+        n_spatial_orbitals, conserves_particle_number=False, real=real, seed=seed_down
+    )
+    up_energies, up_matrix, up_constant = quad_ham_up.diagonalizing_bogoliubov_transform()
+    down_energies, down_matrix, down_constant = quad_ham_down.diagonalizing_bogoliubov_transform()
+
+    n = n_spatial_orbitals
+    n_qubits = 2 * n
+    transformation_matrix = numpy.zeros((n_qubits, 2 * n_qubits), dtype=complex)
+    transformation_matrix[:n, :n] = up_matrix[:, :n]
+    transformation_matrix[:n, n_qubits : n_qubits + n] = up_matrix[:, n:]
+    transformation_matrix[n:, n:n_qubits] = down_matrix[:, :n]
+    transformation_matrix[n:, n_qubits + n :] = down_matrix[:, n:]
+    return (
+        transformation_matrix,
+        (quad_ham_up, up_energies, up_constant),
+        (quad_ham_down, down_energies, down_constant),
+    )
+
+
+def test_bogoliubov_transform_spin_block_gaussian_regression():
+    """Regression test for issue #776.
+
+    A non-square transformation matrix that does not mix spin sectors used to
+    be misidentified as a square spin-block-diagonal matrix, leading to bogus
+    slicing and a ValueError when constructing the circuit.
+    """
+    qubits = LineQubit.range(2)
+    phase = -9.57167901e-01 - 2.89533435e-01j
+    phase /= abs(phase)
+    transformation_matrix = numpy.array([[phase, 0, 0, 0], [0, 1, 0, 0]], dtype=complex)
+
+    circuit = cirq.Circuit(
+        bogoliubov_transform(qubits, transformation_matrix), cirq.I.on_each(*qubits)
+    )
+
+    # The transformation is a phase rotation of mode 0, so it should leave
+    # computational basis states invariant up to phase.
+    state = circuit.final_state_vector(initial_state=0, qubit_order=qubits, dtype=numpy.complex128)
+    expected_state = numpy.zeros(4, dtype=numpy.complex128)
+    expected_state[0] = 1
+    cirq.testing.assert_allclose_up_to_global_phase(state, expected_state, atol=1e-6)
+
+
+@pytest.mark.parametrize(
+    'n_spatial_orbitals, real, up_orbitals, down_orbitals',
+    [(2, True, [0], [0, 1]), (2, False, [1], []), (3, True, [0, 2], [1]), (3, False, [], [0, 2])],
+)
+def test_bogoliubov_transform_spin_block_gaussian(
+    n_spatial_orbitals, real, up_orbitals, down_orbitals, atol=5e-5
+):
+    """A non-particle-conserving transform that does not mix spin sectors must
+    prepare eigenstates of the combined Hamiltonian with the correct energy."""
+    n_qubits = 2 * n_spatial_orbitals
+    qubits = LineQubit.range(n_qubits)
+    sim = cirq.Simulator(dtype=numpy.complex128)
+
+    transformation_matrix, up_data, down_data = _independent_spin_sector_transform(
+        n_spatial_orbitals, 51624, 48201, real
+    )
+    quad_ham_up, up_energies, up_constant = up_data
+    quad_ham_down, down_energies, down_constant = down_data
+
+    # Combined Hamiltonian with spin-down modes shifted to come after spin-up
+    ferm_op = openfermion.get_fermion_operator(quad_ham_up)
+    for term, coefficient in openfermion.get_fermion_operator(quad_ham_down).terms.items():
+        shifted_term = tuple((index + n_spatial_orbitals, action) for index, action in term)
+        ferm_op += openfermion.FermionOperator(shifted_term, coefficient)
+    combined_ham_sparse = get_sparse_operator(ferm_op, n_qubits=n_qubits)
+
+    energy = (
+        sum(up_energies[i] for i in up_orbitals)
+        + sum(down_energies[i] for i in down_orbitals)
+        + up_constant
+        + down_constant
+    )
+    occupied_orbitals = list(up_orbitals) + [i + n_spatial_orbitals for i in down_orbitals]
+    initial_state = sum(2 ** (n_qubits - 1 - int(i)) for i in occupied_orbitals)
+
+    circuit = cirq.Circuit(
+        bogoliubov_transform(qubits, transformation_matrix, initial_state=initial_state)
+    )
+    state = sim.simulate(
+        circuit, initial_state=initial_state, qubit_order=qubits
+    ).final_state_vector
+
+    numpy.testing.assert_allclose(combined_ham_sparse.dot(state), energy * state, atol=atol)
+
+
+@pytest.mark.parametrize(
+    'n_spatial_orbitals, seed_up, seed_down, real',
+    [
+        # The (2, 1000, ...) case has a parity-preserving spin-up sector while
+        # the (2, 1001, ...) case has a parity-flipping one; the latter is the
+        # regression guard for the spin-down Jordan-Wigner sign correction.
+        (2, 1000, 2000, True),
+        (2, 1001, 2001, True),
+        (3, 1000, 2000, False),
+        (3, 1001, 2001, True),
+    ],
+)
+def test_bogoliubov_transform_spin_block_operator_algebra(
+    n_spatial_orbitals, seed_up, seed_down, real, atol=1e-6
+):
+    """The factorized circuit must implement the documented transformation
+    exactly: U a_p^dagger U^-1 = b_p^dagger for every mode p, including the
+    cross-sector Jordan-Wigner sign carried by spin-down operators. Eigenstate
+    energies cannot detect a wrong sign on b_p^dagger, so this checks the full
+    operator algebra of the circuit's unitary directly.
+    """
+    n_qubits = 2 * n_spatial_orbitals
+    qubits = LineQubit.range(n_qubits)
+    transformation_matrix, _, _ = _independent_spin_sector_transform(
+        n_spatial_orbitals, seed_up, seed_down, real
+    )
+
+    circuit = cirq.Circuit(
+        bogoliubov_transform(qubits, transformation_matrix), cirq.I.on_each(*qubits)
+    )
+    unitary = circuit.unitary(qubit_order=qubits)
+
+    def ladder(mode, action):
+        return get_sparse_operator(
+            openfermion.FermionOperator(((mode, action),)), n_qubits=n_qubits
+        ).toarray()
+
+    for p in range(n_qubits):
+        transformed = unitary @ ladder(p, 1) @ unitary.conj().T
+        expected = numpy.zeros_like(transformed)
+        for q in range(n_qubits):
+            expected += transformation_matrix[p, q] * ladder(q, 1)
+            expected += transformation_matrix[p, n_qubits + q] * ladder(q, 0)
+        numpy.testing.assert_allclose(transformed, expected, atol=atol)
+
+
+def test_bogoliubov_transform_spin_mixing_gaussian_not_split(atol=5e-5):
+    """A transform whose annihilation block mixes spin sectors must use the
+    general procedure even if its creation block is spin-block-diagonal."""
+    n_qubits = 4
+    qubits = LineQubit.range(n_qubits)
+    sim = cirq.Simulator(dtype=numpy.complex128)
+
+    # A pairing Hamiltonian whose antisymmetric part couples the two sectors
+    quad_ham = random_quadratic_hamiltonian(
+        n_qubits, conserves_particle_number=False, real=True, seed=63003
+    )
+    quad_ham_sparse = get_sparse_operator(quad_ham)
+    energies, transformation_matrix, constant = quad_ham.diagonalizing_bogoliubov_transform()
+
+    occupied_orbitals = [0, 2]
+    energy = sum(energies[i] for i in occupied_orbitals) + constant
+    initial_state = sum(2 ** (n_qubits - 1 - int(i)) for i in occupied_orbitals)
+
+    circuit = cirq.Circuit(
+        bogoliubov_transform(qubits, transformation_matrix, initial_state=initial_state)
+    )
+    state = sim.simulate(
+        circuit, initial_state=initial_state, qubit_order=qubits
+    ).final_state_vector
+
+    numpy.testing.assert_allclose(quad_ham_sparse.dot(state), energy * state, atol=atol)
+
+
 @pytest.mark.parametrize(
     'n_qubits, conserves_particle_number', [(4, True), (4, False), (5, True), (5, False)]
 )