diff --git a/src/openfermion/transforms/opconversions/remove_symmetry_qubits.py b/src/openfermion/transforms/opconversions/remove_symmetry_qubits.py
index 5356d8d54..41ec02120 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
@@ -99,7 +99,17 @@ def symmetry_conserving_bravyi_kitaev(fermion_hamiltonian, active_orbitals, acti
     )
     qubit_hamiltonian = remove_indices(qubit_hamiltonian, (active_orbitals / 2, active_orbitals))
 
-    return qubit_hamiltonian
+    # 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 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..0f59d497c 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,27 @@ 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), 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)
+
+        # Every term must be in canonical form: 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 ...
+        sparse_op = get_sparse_operator(trafo_op)
+
+        # ... and must match the explicitly-simplified operator.
+        simplified = QubitOperator()
+        for term, coefficient in trafo_op.terms.items():
+            simplified += QubitOperator(term, coefficient)
+        expected = get_sparse_operator(simplified)
+
+        self.assertTrue(numpy.allclose(sparse_op.toarray(), expected.toarray()))