Fix partial_trace when not removing subsystems (#7613)

* Fix `partial_trace` when not removing subsystems

There were places, such as in `quantum_info.partial_trace` where an
`OpShape` output of `OpShape.remove` was mutated.  In cases where
nothing was removed from the `OpShape`, it returned itself, and
consequently caused lingering issues with the object that owned the
shape.  This ensures a copy is always produced, and adds in a special
case to handle `DensityMatrix` with nothing traced out.

* Check qargs not shape equality

Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

qiskit/quantum_info/operators/op_shape.py

@@ -332,7 +332,7 @@ def subset(self, qargs=None, qargs_l=None, qargs_r=None):
         )
 
     def remove(self, qargs=None, qargs_l=None, qargs_r=None):
-        """Return the reduced OpShape with specified qargs removed"""
+        """Return a new :class:`OpShape` with the specified qargs removed"""
         if qargs:
             # Convert qargs to left and right qargs
             if qargs_l or qargs_r:
@@ -342,7 +342,7 @@ def remove(self, qargs=None, qargs_l=None, qargs_r=None):
             if self._num_qargs_r:
                 qargs_r = qargs
         if qargs_l is None and qargs_r is None:
-            return self
+            return self.copy()
 
         # Format integer qargs
         if isinstance(qargs_l, Integral):

qiskit/quantum_info/states/utils.py

@@ -60,6 +60,9 @@ def partial_trace(state, qargs):
         return DensityMatrix(rho, dims=traced_shape._dims_l)
 
     # Density matrix case
+    # Empty partial trace case.
+    if not qargs:
+        return state.copy()
     # Trace first subsystem to avoid coping whole density matrix
     dims = state.dims(qargs)
     tr_op = SuperOp(np.eye(dims[0]).reshape(1, dims[0] ** 2), input_dims=[dims[0]], output_dims=[1])

releasenotes/notes/fix-partial_trace-no-systems-0dc2df3007942eb6.yaml

@@ -0,0 +1,7 @@
+---
+fixes:
+  - |
+    If :func:`.quantum_info.partial_trace` is asked to trace out *no* subsystems,
+    it will now correctly return the
+    :class:`.DensityMatrix` of the input state with all dimensions remaining
+    rather than throwing an error.

test/python/quantum_info/states/test_utils.py

@@ -32,6 +32,7 @@ def test_statevector_partial_trace(self):
         self.assertEqual(partial_trace(psi, [0]), DensityMatrix.from_label("10"))
         self.assertEqual(partial_trace(psi, [1]), DensityMatrix.from_label("1+"))
         self.assertEqual(partial_trace(psi, [2]), DensityMatrix.from_label("0+"))
+        self.assertEqual(partial_trace(psi, []), DensityMatrix(psi))
 
     def test_density_matrix_partial_trace(self):
         """Test partial_trace function on density matrices"""
@@ -42,6 +43,7 @@ def test_density_matrix_partial_trace(self):
         self.assertEqual(partial_trace(rho, [0]), DensityMatrix.from_label("10"))
         self.assertEqual(partial_trace(rho, [1]), DensityMatrix.from_label("1+"))
         self.assertEqual(partial_trace(rho, [2]), DensityMatrix.from_label("0+"))
+        self.assertEqual(partial_trace(rho, []), rho)
 
     def test_shannon_entropy(self):
         """Test shannon_entropy function"""