Repository: qutip/qutip (python)

Issue:
Support measurments of observables with degenerate eigenvalues
Currently the `qutip.measurement.measure_observable` function always uses rank-1 projectors to determine the final state after the measurement. But for obervables with degenerate eigenvalues you might just want to project the state onto a higher-dimensional eigenspace. Eg. measuring the $Z \otimes Z$ operator on the Bell state $\ket{\psi} = 1/\sqrt{2} (\ket{00} + \ket{11})$ using qutip will return either $\ket{00}$ or $\ket{11}$ despite the state being an eigenstate of the observable.
```
import qutip as qt
import qutip.measurement as meas
import numpy as np

#Prepare Bell state
state = qt.tensor([qt.basis(2, 0)]*2) + qt.tensor([qt.basis(2, 1)]*2)
state = state.unit()

# ZZ operator
op = qt.tensor([qt.sigmaz()]*2)

result_eigenvalue, final_state = meas.measure(state, op)

#Check that the measurement didn't change the state (because it is an eigenstate of the measurement operator)
assert np.isclose(qt.fidelity(final_state, state), 1) #fails
```

I am aware that I can get my intended behaviour by explicity passing the list of projection operators to the `measure` function, but the following enhancement would be very useful:
An option `respect_degenerate_eigenvalues=False` in `qutip.measurement.measure_observable` that groups the eigenstate projectors based on eigenvalue degeneracy. In the above example the rank-2 projectors |00><00| + |11><11| and |01><01| + |10><10| would be used and the Bell state would be preseved during a ZZ measurement.

Pull request:
Makes measurement support degenerate states
**Description**
Measurement functions did not take into account degenerate states of the operator. (#1800, #2054)
Also the returned probabilities could be negative. (#2265)

I added a tolerance to these function for the minimal probability. Smaller / negative probabilities are rounded to `0`.
I also changed `measurement_statistics_observable` to support degenerate operators: It always return the projectors, with projectors and probabilities of degenerate sub-systems summed. It also only return the (eigenvalues, projector, probability) set for cases where the probability is greater than `tol`.

**Related issues or PRs**
fix #1800, #2054, #2029, #2265

Patch:
diff --git a/doc/changes/2331.misc b/doc/changes/2331.misc
deleted file mode 100644
index 4231d540f..000000000
--- a/doc/changes/2331.misc
+++ /dev/null
@@ -1,1 +0,0 @@
-Fix typos in `expect` documentation
diff --git a/doc/changes/2333.bugfix b/doc/changes/2333.bugfix
deleted file mode 100644
index c8a1bda03..000000000
--- a/doc/changes/2333.bugfix
+++ /dev/null
@@ -1,1 +0,0 @@
-Fixed two problems with the steady_state() solver in the HEOM method.
\ No newline at end of file
diff --git a/doc/changes/2342.misc b/doc/changes/2342.misc
new file mode 100644
index 000000000..1163fc45c
--- /dev/null
+++ b/doc/changes/2342.misc
@@ -0,0 +1,1 @@
+Allow measurement functions to support degenerate operators.
diff --git a/qutip/__init__.py b/qutip/__init__.py
index 19f373c8b..3fe38033a 100644
--- a/qutip/__init__.py
+++ b/qutip/__init__.py
@@ -50,7 +50,7 @@ from .entropy import *
 from .partial_transpose import *
 from .continuous_variables import *
 from .distributions import *
-
+from . import measurement
 
 # utilities
 from .utilities import *
diff --git a/qutip/core/expect.py b/qutip/core/expect.py
index bb72fccd9..9051c53f7 100644
--- a/qutip/core/expect.py
+++ b/qutip/core/expect.py
@@ -16,7 +16,7 @@ def expect(oper, state):
     Parameters
     ----------
     oper : qobj/array-like
-        A single or a `list` of operators for expectation value.
+        A single or a `list` or operators for expectation value.
 
     state : qobj/array-like
         A single or a `list` of quantum states or density matrices.
@@ -25,8 +25,8 @@ def expect(oper, state):
     -------
     expt : float/complex/array-like
         Expectation value.  ``real`` if ``oper`` is Hermitian, ``complex``
-        otherwise. A (nested) array of expectaction values if ``state`` or
-        ``oper`` are arrays.
+        otherwise. A (nested) array of expectaction values of state or operator
+        are arrays.
 
     Examples
     --------
diff --git a/qutip/core/qobj.py b/qutip/core/qobj.py
index c050143ec..2b4dcda24 100644
--- a/qutip/core/qobj.py
+++ b/qutip/core/qobj.py
@@ -147,10 +147,12 @@ class Qobj:
 
     Parameters
     ----------
-    inpt: array_like, data object or :obj:`.Qobj`
+    inpt: array_like
         Data for vector/matrix representation of the quantum object.
     dims: list
         Dimensions of object used for tensor products.
+    type: {'bra', 'ket', 'oper', 'operator-ket', 'operator-bra', 'super'}
+        The type of quantum object to be represented.
     shape: list
         Shape of underlying data structure (matrix shape).
     copy: bool
@@ -160,12 +162,8 @@ class Qobj:
 
     Attributes
     ----------
-    data : object
-        The data object storing the vector / matrix representation of the
-        `Qobj`.
-    dtype : type
-        The data-layer type used for storing the data. The possible types are
-        described in `Qobj.to <./classes.html#qutip.core.qobj.Qobj.to>`__.
+    data : array_like
+        Sparse matrix characterizing the quantum object.
     dims : list
         List of dimensions keeping track of the tensor structure.
     shape : list
@@ -175,8 +173,7 @@ class Qobj:
         'operator-bra', or 'super'.
     superrep : str
         Representation used if `type` is 'super'. One of 'super'
-        (Liouville form), 'choi' (Choi matrix with tr = dimension),
-        or 'chi' (chi-matrix representation).
+        (Liouville form) or 'choi' (Choi matrix with tr = dimension).
     isherm : bool
         Indicates if quantum object represents Hermitian operator.
     isunitary : bool
@@ -214,16 +211,10 @@ class Qobj:
         Create copy of Qobj
     conj()
         Conjugate of quantum object.
-    contract()
-        Contract subspaces of the tensor structure which are 1D.
     cosm()
         Cosine of quantum object.
     dag()
         Adjoint (dagger) of quantum object.
-    data_as(format, copy)
-        Vector / matrix representation of quantum object.
-    diag()
-        Diagonal elements of quantum object.
     dnorm()
         Diamond norm of quantum operator.
     dual_chan()
@@ -241,14 +232,10 @@ class Qobj:
         object.
     inv()
         Return a Qobj corresponding to the matrix inverse of the operator.
-    logm()
-        Matrix logarithm of quantum operator.
     matrix_element(bra, ket)
         Returns the matrix element of operator between `bra` and `ket` vectors.
     norm(norm='tr', sparse=False, tol=0, maxiter=100000)
         Returns norm of a ket or an operator.
-    overlap(other)
-        Overlap between two state vectors or two operators.
     permute(order)
         Returns composite qobj with indices reordered.
     proj()
@@ -256,8 +243,6 @@ class Qobj:
     ptrace(sel)
         Returns quantum object for selected dimensions after performing
         partial trace.
-    purity()
-        Calculates the purity of a quantum object.
     sinm()
         Sine of quantum object.
     sqrtm()
@@ -355,7 +340,7 @@ class Qobj:
             raise TypeError('Qobj data must be a data-layer format.')
         if self._dims.shape != data.shape:
             raise ValueError('Provided data do not match the dimensions: ' +
-                             f"{self._dims.shape} vs {data.shape}")
+                             f"{dims.shape} vs {self._data.shape}")
         self._data = data
 
     def to(self, data_type):
@@ -919,6 +904,7 @@ class Qobj:
         else:
             evals, evecs = _data.eigs(self.data, isherm=self._isherm)
 
+        numevals = len(evals)
         dV = _data.diag([np.sqrt(evals, dtype=complex)], 0)
         if self.isherm:
             spDv = _data.matmul(dV, evecs.conj().transpose())
diff --git a/qutip/measurement.py b/qutip/measurement.py
index da90387dd..5b29c4083 100644
--- a/qutip/measurement.py
+++ b/qutip/measurement.py
@@ -13,7 +13,7 @@ __all__ = [
 
 import numpy as np
 
-from . import Qobj, expect, identity, tensor
+from . import Qobj, expect, identity, tensor, settings
 
 
 def _verify_input(op, state):
@@ -36,7 +36,7 @@ def _verify_input(op, state):
         raise ValueError("state must be a ket or a density matrix")
 
 
-def _measurement_statistics_povm_ket(state, ops):
+def _measurement_statistics_povm_ket(state, ops, tol=None):
     r"""
     Returns measurement statistics (resultant states and probabilities)
     for a measurements specified by a set of positive operator valued
@@ -51,6 +51,9 @@ def _measurement_statistics_povm_ket(state, ops):
       List of measurement operators :math:`M_i` (specifying a POVM such that
       :math:`E_i = M_i^\dagger M_i`).
 
+    tol : float, optional
+        Smallest value for the probabilities. Smaller probabilities will be
+        rounded to ``0``.
 
     Returns
     -------
@@ -65,19 +68,23 @@ def _measurement_statistics_povm_ket(state, ops):
     """
     probabilities = []
     collapsed_states = []
+    if tol is None:
+        tol = settings.core["atol"]
 
     for i, op in enumerate(ops):
-        p = np.absolute((state.dag() * op.dag() * op * state))
-        probabilities.append(p)
-        if p != 0:
-            collapsed_states.append((op * state) / np.sqrt(p))
+        psi = op * state
+        p = np.absolute(psi.overlap(psi))
+        if p >= tol:
+            collapsed_states.append(psi / np.sqrt(p))
+            probabilities.append(p)
         else:
             collapsed_states.append(None)
+            probabilities.append(0.)
 
     return collapsed_states, probabilities
 
 
-def _measurement_statistics_povm_dm(density_mat, ops):
+def _measurement_statistics_povm_dm(density_mat, ops, tol=None):
     r"""
     Returns measurement statistics (resultant states and probabilities)
     for a measurements specified by a set of positive operator valued
@@ -92,6 +99,10 @@ def _measurement_statistics_povm_dm(density_mat, ops):
         List of measurement operators :math:`M_i` (specifying a POVM s.t.
         :mathm:`E_i = M_i^\dagger M_i`)
 
+    tol : float, optional
+        Smallest value for the probabilities. Smaller probabilities will be
+        rounded to ``0``.
+
     Returns
     -------
     collapsed_states : list of :class:`.Qobj`
@@ -105,20 +116,24 @@ def _measurement_statistics_povm_dm(density_mat, ops):
     """
     probabilities = []
     collapsed_states = []
+    if tol is None:
+        tol = settings.core["atol"]
 
     for i, op in enumerate(ops):
         st = op * density_mat * op.dag()
         p = st.tr()
-        probabilities.append(p)
-        if p != 0:
+        if p >= tol:
             collapsed_states.append(st/p)
+            probabilities.append(p)
         else:
             collapsed_states.append(None)
+            probabilities.append(0.)
+
 
     return collapsed_states, probabilities
 
 
-def measurement_statistics_povm(state, ops):
+def measurement_statistics_povm(state, ops, tol=None):
     r"""
     Returns measurement statistics (resultant states and probabilities) for a
     measurement specified by a set of positive operator valued measurements on
@@ -137,6 +152,10 @@ def measurement_statistics_povm(state, ops):
            projectors (s.t. :math:`E_i = M_i^\dagger = M_i`)
         3. kets (transformed to projectors)
 
+    tol : float, optional
+        Smallest value for the probabilities. Smaller probabilities will be
+        rounded to ``0``. Default is qutip's core settings' ``atol``.
+
     Returns
     -------
     collapsed_states : list of :class:`.Qobj`
@@ -160,12 +179,12 @@ def measurement_statistics_povm(state, ops):
         raise ValueError("measurement operators must sum to identity")
 
     if state.isket:
-        return _measurement_statistics_povm_ket(state, ops)
+        return _measurement_statistics_povm_ket(state, ops, tol)
     else:
-        return _measurement_statistics_povm_dm(state, ops)
+        return _measurement_statistics_povm_dm(state, ops, tol)
 
 
-def measurement_statistics_observable(state, op):
+def measurement_statistics_observable(state, op, tol=None):
     """
     Return the measurement eigenvalues, eigenstates (or projectors) and
     measurement probabilities for the given state and measurement operator.
@@ -178,33 +197,57 @@ def measurement_statistics_observable(state, op):
     op : :class:`.Qobj`
         The measurement operator.
 
+    tol : float, optional
+        Smallest value for the probabilities.
+        Default is qutip's core settings' ``atol``.
+
     Returns
     -------
     eigenvalues: list of float
         The list of eigenvalues of the measurement operator.
 
-    eigenstates_or_projectors: list of :class:`.Qobj`
-        If the state was a ket, return the eigenstates of the measurement
-        operator. Otherwise return the projectors onto the eigenstates.
+    projectors: list of :class:`.Qobj`
+        Return the projectors onto the eigenstates.
 
     probabilities: list of float
         The probability of measuring the state as being in the corresponding
         eigenstate (and the measurement result being the corresponding
         eigenvalue).
     """
+    probabilities = []
+    values = []
+    projectors = []
     _verify_input(op, state)
+    if tol is None:
+        tol = settings.core["atol"]
 
     eigenvalues, eigenstates = op.eigenstates()
-    if state.isket:
-        probabilities = [abs(e.overlap(state))**2 for e in eigenstates]
-        return eigenvalues, eigenstates, probabilities
-    else:
-        projectors = [e.proj() for e in eigenstates]
-        probabilities = [expect(v, state) for v in projectors]
-        return eigenvalues, projectors, probabilities
+    # Detect groups of eigenvalues within atol of each other.
+    # A group will be [False] * N + [True]
+    groups = np.append(np.diff(eigenvalues) >= tol, True)
+
+    present_group = []
+    for i in range(len(eigenvalues)):
+        present_group.append(i)
+        if not groups[i]:
+            continue
 
+        projector = 0
+        for j in present_group:
+            projector += eigenstates[j].proj()
+        probability = expect(projector, state)
 
-def measure_observable(state, op):
+        if probability >= tol:
+            probabilities.append(probability)
+            values.append(np.mean(eigenvalues[present_group]))
+            projectors.append(projector)
+
+        present_group = []
+
+    return values, projectors, probabilities
+
+
+def measure_observable(state, op, tol=None):
     """
     Perform a measurement specified by an operator on the given state.
 
@@ -221,6 +264,10 @@ def measure_observable(state, op):
     op : :class:`.Qobj`
         The measurement operator.
 
+    tol : float, optional
+        Smallest value for the probabilities.
+        Default is qutip's core settings' ``atol``.
+
     Returns
     -------
     measured_value : float
@@ -269,19 +316,17 @@ def measure_observable(state, op):
     The measurement result is the same, but the new state is returned as a
     density matrix.
     """
-    eigenvalues, eigenstates_or_projectors, probabilities = (
-        measurement_statistics_observable(state, op))
-    i = np.random.choice(range(len(eigenvalues)), p=probabilities)
+    eigenvalues, projectors, probabilities = (
+        measurement_statistics_observable(state, op, tol))
+    i = np.random.choice(len(eigenvalues), p=probabilities)
     if state.isket:
-        eigenstates = eigenstates_or_projectors
-        state = eigenstates[i]
+        state = (projectors[i] * state) / probabilities[i]**0.5
     else:
-        projectors = eigenstates_or_projectors
         state = (projectors[i] * state * projectors[i]) / probabilities[i]
     return eigenvalues[i], state
 
 
-def measure_povm(state, ops):
+def measure_povm(state, ops, tol=None):
     r"""
     Perform a measurement specified by list of POVMs.
 
@@ -303,6 +348,10 @@ def measure_povm(state, ops):
            :math:`E_i = M_i^\dagger = M_i`)
         3. kets (transformed to projectors)
 
+    tol : float, optional
+        Smallest value for the probabilities.
+        Default is qutip's core settings' ``atol``.
+
     Returns
     -------
     index : float
@@ -311,13 +360,14 @@ def measure_povm(state, ops):
     state : :class:`.Qobj`
         The new state (a ket if a ket was given, otherwise a density matrix).
     """
-    collapsed_states, probabilities = measurement_statistics_povm(state, ops)
-    index = np.random.choice(range(len(collapsed_states)), p=probabilities)
+    collapsed_states, probabilities = (
+        measurement_statistics_povm(state, ops, tol))
+    index = np.random.choice(len(collapsed_states), p=probabilities)
     state = collapsed_states[index]
     return index, state
 
 
-def measurement_statistics(state, ops):
+def measurement_statistics(state, ops, tol=None):
     r"""
     A dispatch method that provides measurement statistics handling both
     observable style measurements and projector style measurements(POVMs and
@@ -342,14 +392,18 @@ def measurement_statistics(state, ops):
           2. projection operators if ops correspond to projectors (s.t.
              :math:`E_i = M_i^\dagger = M_i`)
           3. kets (transformed to projectors)
+
+    tol : float, optional
+        Smallest value for the probabilities.
+        Default is qutip's core settings' ``atol``.
     """
     if isinstance(ops, list):
-        return measurement_statistics_povm(state, ops)
+        return measurement_statistics_povm(state, ops, tol)
     else:
-        return measurement_statistics_observable(state, ops)
+        return measurement_statistics_observable(state, ops, tol)
 
 
-def measure(state, ops):
+def measure(state, ops, tol=None):
     r"""
     A dispatch method that provides measurement results handling both
     observable style measurements and projector style measurements (POVMs and
@@ -374,8 +428,12 @@ def measure(state, ops):
           2. projection operators if ops correspond to projectors (s.t.
              :math:`E_i = M_i^\dagger = M_i`)
           3. kets (transformed to projectors)
+
+    tol : float, optional
+        Smallest value for the probabilities.
+        Default is qutip's core settings' ``atol``.
     """
     if isinstance(ops, list):
-        return measure_povm(state, ops)
+        return measure_povm(state, ops, tol)
     else:
-        return measure_observable(state, ops)
+        return measure_observable(state, ops, tol)
diff --git a/qutip/solver/heom/bofin_solvers.py b/qutip/solver/heom/bofin_solvers.py
index 1061c1c49..23b82f629 100644
--- a/qutip/solver/heom/bofin_solvers.py
+++ b/qutip/solver/heom/bofin_solvers.py
@@ -966,8 +966,8 @@ class HEOMSolver(Solver):
             L = L.tocsc()
             solution = spsolve(L, b_mat)
 
-        data = _data.Dense(solution[:n ** 2].reshape((n, n), order='F'))
-        data = _data.mul(_data.add(data, data.adjoint()), 0.5)
+        data = _data.Dense(solution[:n ** 2].reshape((n, n)))
+        data = _data.mul(_data.add(data, data.conj()), 0.5)
         steady_state = Qobj(data, dims=self._sys_dims)
 
         solution = solution.reshape((self._n_ados, n, n))


Test patch:
diff --git a/qutip/tests/solver/heom/test_bofin_solvers.py b/qutip/tests/solver/heom/test_bofin_solvers.py
index b704b015d..ad6e5f229 100644
--- a/qutip/tests/solver/heom/test_bofin_solvers.py
+++ b/qutip/tests/solver/heom/test_bofin_solvers.py
@@ -8,8 +8,8 @@ from numpy.linalg import eigvalsh
 from scipy.integrate import quad
 
 from qutip import (
-    basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz, sigmay,
-    tensor, Qobj, QobjEvo, fidelity
+    basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz,
+    tensor, Qobj, QobjEvo
 )
 from qutip.core import data as _data
 from qutip.solver.heom.bofin_baths import (
@@ -741,24 +741,6 @@ class TestHEOMSolver:
         else:
             assert_raises_steady_state_time_dependent(hsolver)
 
-    def test_steady_state_bosonic_bath(
-        self, atol=1e-3
-    ):
-        H_sys = 0.25 * sigmaz() + 0.5 * sigmay()
-
-        bath = DrudeLorentzBath(sigmaz(), lam=0.025,
-                                gamma=0.05, T=1/0.95, Nk=2)
-        options = {"nsteps": 15000, "store_states": True}
-        hsolver = HEOMSolver(H_sys, bath, 5, options=options)
-
-        tlist = np.linspace(0, 500, 21)
-        rho0 = basis(2, 0) * basis(2, 0).dag()
-
-        result = hsolver.run(rho0, tlist)
-        rho_final, ado_state = hsolver.steady_state()
-        fid = fidelity(rho_final, result.states[-1])
-        np.testing.assert_allclose(fid, 1.0, atol=atol)
-
     @pytest.mark.parametrize(['terminator'], [
         pytest.param(True, id="terminator"),
         pytest.param(False, id="noterminator"),
diff --git a/qutip/tests/solver/test_floquet.py b/qutip/tests/solver/test_floquet.py
index 3c990d0c0..679a3cba8 100644
--- a/qutip/tests/solver/test_floquet.py
+++ b/qutip/tests/solver/test_floquet.py
@@ -51,7 +51,7 @@ class TestFloquet:
         states = sesolve(H, psi0, tlist).states
         for t, state in zip(tlist, states):
             from_floquet = floquet_basis.from_floquet_basis(floquet_psi0, t)
-            assert state.overlap(from_floquet) == pytest.approx(1., abs=8e-5)
+            assert state.overlap(from_floquet) == pytest.approx(1., abs=5e-5)
 
     def testFloquetUnitary(self):
         N = 10
diff --git a/qutip/tests/test_measurement.py b/qutip/tests/test_measurement.py
index 6bccd1a35..44b57ce28 100644
--- a/qutip/tests/test_measurement.py
+++ b/qutip/tests/test_measurement.py
@@ -3,8 +3,8 @@ import scipy
 import pytest
 from math import sqrt
 from qutip import (
-    Qobj, basis, ket2dm, sigmax, sigmay, sigmaz, identity, tensor,
-    rand_ket,
+    Qobj, basis, ket2dm, sigmax, sigmay, sigmaz, identity, num, tensor,
+    rand_ket
 )
 from qutip.measurement import (
     measure_povm, measurement_statistics_povm, measure_observable,
@@ -70,8 +70,6 @@ def _equivalent(left, right, tol=1e-8):
 
 
 @pytest.mark.parametrize(["op", "state", "pairs", "probabilities"], [
-    pytest.param(sigmaz(), basis(2, 0), SIGMAZ, [0, 1], id="sigmaz_ket"),
-    pytest.param(sigmaz(), basis(2, 0).proj(), SIGMAZ, [0, 1], id="sigmaz_dm"),
     pytest.param(sigmax(), basis(2, 0), SIGMAX, [0.5, 0.5], id="sigmax_ket"),
     pytest.param(sigmax(), basis(2, 0).proj(), SIGMAX, [0.5, 0.5],
                  id="sigmax_dm"),
@@ -82,19 +80,24 @@ def _equivalent(left, right, tol=1e-8):
 def test_measurement_statistics_observable(op, state, pairs, probabilities):
     """ measurement_statistics_observable: observables on basis states. """
 
-    evs, ess_or_projs, probs = measurement_statistics_observable(state, op)
-    np.testing.assert_almost_equal(evs, pairs.eigenvalues)
-    if state.isket:
-        ess = ess_or_projs
-        assert len(ess) == len(pairs.eigenstates)
-        for a, b in zip(ess, pairs.eigenstates):
-            assert _equivalent(a, b)
-    else:
-        projs = ess_or_projs
-        assert len(projs) == len(pairs.projectors)
-        for a, b in zip(projs, pairs.projectors):
-            assert _equivalent(a, b)
+    evs, projs, probs = measurement_statistics_observable(state, op)
+    assert len(projs) == len(probabilities)
     np.testing.assert_almost_equal(probs, probabilities)
+    for a, b in zip(projs, pairs.projectors):
+        assert _equivalent(a, b)
+
+
+def test_measurement_statistics_observable_degenerate():
+    """ measurement_statistics_observable: observables on basis states. """
+
+    state = basis(2, 1) & (basis(2, 0) + basis(2, 1)).unit()
+    op = sigmaz() & identity(2)
+    expected_projector = num(2) & identity(2)
+    evs, projs, probs = measurement_statistics_observable(state, op)
+    assert len(probs) == 1
+    np.testing.assert_almost_equal(probs, [1.])
+    np.testing.assert_almost_equal(evs, [-1.])
+    assert _equivalent(projs[0], expected_projector)
 
 
 @pytest.mark.parametrize(["ops", "state", "final_states", "probabilities"], [
diff --git a/qutip/tests/test_random.py b/qutip/tests/test_random.py
index 085e94044..8d6ff441e 100644
--- a/qutip/tests/test_random.py
+++ b/qutip/tests/test_random.py
@@ -222,7 +222,7 @@ def test_rand_super(dimensions, dtype, superrep):
     """
     random_qobj = rand_super(dimensions, dtype=dtype, superrep=superrep)
     assert random_qobj.issuper
-    with CoreOptions(atol=2e-9):
+    with CoreOptions(atol=1e-9):
         assert random_qobj.iscptp
     assert random_qobj.superrep == superrep
     _assert_metadata(random_qobj, dimensions, dtype, super=True)