Fix #880: return simplified QubitOperator from symmetry_conserving_bravyi_kitaev (#1440)

Closes #880.

## Problem

`get_sparse_operator` raises on the operator produced by
`symmetry_conserving_bravyi_kitaev`:

```python
from openfermion import get_sparse_operator, FermionOperator, symmetry_conserving_bravyi_kitaev

fermion_op = FermionOperator("0^ 1^")
qubit_op = symmetry_conserving_bravyi_kitaev(fermion_op, 4, 2)
get_sparse_operator(qubit_op)   # ValueError: axis 0 index 7 exceeds matrix dimension 4
```

## Root cause

`symmetry_conserving_bravyi_kitaev` finishes by calling
`remove_indices`, which shifts qubit indices. When two qubits are mapped
onto the **same** new index, a term ends up with multiple Paulis acting
on one qubit, e.g.:

```
((0, X), (1, Y), (1, X))
```

`remove_indices` writes these terms straight into the operator’s
`.terms` dict, so they bypass the simplification `QubitOperator`
normally performs on construction. `qubit_operator_sparse` assumes each
qubit appears at most once per term (it grows the tensor product one
factor per Pauli), so a repeated qubit makes the per-term matrix larger
than the `n_qubits` Hilbert space and the assembly raises.

## Fix

Rebuild the operator after `remove_indices` so every term is routed back
through `QubitOperator`’s simplification, restoring canonical
one-Pauli-per-qubit form (e.g. `((0, X), (1, Y), (1, X)) → ((0, X), (1,
Z))`). This keeps the operator mathematically identical — the existing
eigenspectrum-based tests still pass — while satisfying the invariant
that `QubitOperator`s are simplified, which `get_sparse_operator` and
other consumers rely on.

## Tests

Added `test_output_is_simplified_qubit_operator`, a regression test
built from the reporter's reproducer. It now asserts the transform
equals an independently worked-out, fully-simplified `QubitOperator` for
that input (rather than re-deriving the expectation from the output),
checks every term is canonical (one Pauli per qubit), and checks
`get_sparse_operator` no longer raises. The existing eigenspectrum-based
tests continue to pass, confirming the operator is unchanged
mathematically.

Verified with the repo's check scripts: `check/pytest -m "not slow"
src/openfermion/transforms/opconversions/remove_symmetry_qubits_test.py`
passes, and `check/format-incremental` and `check/pylint-changed-files`
are clean. `check/mypy` reports no issues in the changed module.

---------

Co-authored-by: stark256-spec <stark256-spec@users.noreply.github.com>

Author: stark256-spec

src/openfermion/transforms/opconversions/remove_symmetry_qubits.py

@@ -16,7 +16,7 @@
 
 import copy
 
-from openfermion.ops.operators import FermionOperator
+from openfermion.ops.operators import FermionOperator, QubitOperator
 from openfermion.transforms.opconversions.bravyi_kitaev_tree import bravyi_kitaev_tree
 from openfermion.transforms.opconversions.term_reordering import reorder
 from openfermion.utils.indexing import up_then_down
@@ -95,11 +95,21 @@ def symmetry_conserving_bravyi_kitaev(fermion_hamiltonian, active_orbitals, acti
         qubit_hamiltonian, active_orbitals, parity_final_orb
     )
     qubit_hamiltonian = edit_hamiltonian_for_spin(
-        qubit_hamiltonian, active_orbitals / 2, parity_middle_orb
+        qubit_hamiltonian, active_orbitals // 2, parity_middle_orb
     )
-    qubit_hamiltonian = remove_indices(qubit_hamiltonian, (active_orbitals / 2, active_orbitals))
+    qubit_hamiltonian = remove_indices(qubit_hamiltonian, (active_orbitals // 2, active_orbitals))
+
+    # remove_indices() shifts qubit indices and can map two qubits onto the
+    # same index, producing terms with multiple Paulis acting on one qubit
+    # (e.g. ((0, 'X'), (1, 'Y'), (1, 'X'))). Rebuilding the operator routes
+    # each term back through QubitOperator's simplification, restoring the
+    # canonical one-Pauli-per-qubit form that consumers such as
+    # get_sparse_operator require. See issue #880.
+    simplified_hamiltonian = QubitOperator()
+    for term, coefficient in qubit_hamiltonian.terms.items():
+        simplified_hamiltonian += QubitOperator(term, coefficient)
 
-    return qubit_hamiltonian
+    return simplified_hamiltonian
 
 
 def edit_hamiltonian_for_spin(qubit_hamiltonian, spin_orbital, orbital_parity):

src/openfermion/transforms/opconversions/remove_symmetry_qubits_test.py

@@ -16,6 +16,7 @@
 
 import unittest
 
+import numpy
 import pytest
 
 from openfermion.hamiltonians import fermi_hubbard
@@ -26,7 +27,7 @@
     jw_get_ground_state_at_particle_number,
 )
 from openfermion.linalg import eigenspectrum
-from openfermion.ops.operators import FermionOperator
+from openfermion.ops.operators import FermionOperator, QubitOperator
 
 from openfermion.transforms.opconversions.remove_symmetry_qubits import (
     symmetry_conserving_bravyi_kitaev,
@@ -152,3 +153,36 @@ def test_single_operator(self):
         e_trafo = eigenspectrum(trafo_op)
         # Check eigenvalues
         self.assertSequenceEqual(e_op.tolist(), e_trafo.tolist())
+
+    def test_output_is_simplified_qubit_operator(self):
+        # Regression test for issue #880: symmetry_conserving_bravyi_kitaev
+        # used to return QubitOperators with un-simplified terms (multiple
+        # Paulis acting on the same qubit, e.g. ((0, 'X'), (1, 'Y'), (1, 'X'))
+        # left behind when remove_indices maps two qubit indices onto one),
+        # which made get_sparse_operator raise a ValueError.
+        op = FermionOperator("0^ 1^")
+        trafo_op = symmetry_conserving_bravyi_kitaev(op, active_orbitals=4, active_fermions=2)
+
+        # The result must equal the known, fully-simplified operator for this
+        # input, worked out independently (rather than re-derived from
+        # trafo_op, which would be a no-op now that the output is simplified).
+        expected_op = (
+            QubitOperator(((0, "X"), (1, "Z")), -0.25)
+            + QubitOperator(((0, "X"),), -0.25)
+            + QubitOperator(((0, "Y"), (1, "Z")), 0.25j)
+            + QubitOperator(((0, "Y"),), 0.25j)
+        )
+        self.assertEqual(trafo_op, expected_op)
+
+        # Every term is canonical: at most one Pauli per qubit.
+        for term in trafo_op.terms:
+            qubits = [qubit for qubit, _ in term]
+            self.assertEqual(len(qubits), len(set(qubits)))
+
+        # get_sparse_operator must no longer raise on the result, and must
+        # produce the sparse matrix of the independently-specified operator.
+        self.assertTrue(
+            numpy.allclose(
+                get_sparse_operator(trafo_op).toarray(), get_sparse_operator(expected_op).toarray()
+            )
+        )