diff --git a/src/openfermion/transforms/opconversions/remove_symmetry_qubits.py b/src/openfermion/transforms/opconversions/remove_symmetry_qubits.py
index 5356d8d54..9ff91361c 100644
--- a/src/openfermion/transforms/opconversions/remove_symmetry_qubits.py
+++ b/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):
diff --git a/src/openfermion/transforms/opconversions/remove_symmetry_qubits_test.py b/src/openfermion/transforms/opconversions/remove_symmetry_qubits_test.py
index 1a95e48f3..21257d561 100644
--- a/src/openfermion/transforms/opconversions/remove_symmetry_qubits_test.py
+++ b/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()
+            )
+        )