Merge pull request #423 from ionq/bugfix-kraus-decomp-with-degeneracies

Fixes bug in DenseOperator.kraus_operators when Choi matrix has degen…

pygsti/modelmembers/operations/denseop.py

@@ -416,11 +416,11 @@ def kraus_operators(self):
         #CHECK 1 (to unit test?) REMOVE
         #tmp_std = _bt.change_basis(superop_mx, self._basis, 'std')
         #B = _bt.basis_matrices('std', superop_mx.shape[0])
-        #check_superop = sum([ choi_mx[i,j] * _np.kron(B[i], B[j].T) for i in range(d*d) for j in range(d*d)])
+        #check_superop = sum([ choi_mx[i,j] * _np.kron(B[i], B[j].conjugate()) for i in range(d*d) for j in range(d*d)])
         #assert(_np.allclose(check_superop, tmp_std))
 
-        evals, evecs = _np.linalg.eig(choi_mx)
-        #assert(_np.allclose(evecs @ _np.diag(evals) @ (evecs.conjugate().T), choi_mx))
+        evals, evecs = _np.linalg.eigh(choi_mx)
+        assert(_np.allclose(evecs @ _np.diag(evals) @ (evecs.conjugate().T), choi_mx))
         TOL = 1e-7  # consider lowering this tolerance as it leads to errors of this order in the Kraus decomp
         if any([ev <= -TOL for ev in evals]):
             raise ValueError("Cannot compute Kraus decomposition of non-positive-definite superoperator!")

test/unit/modelmembers/test_kraus_interface.py

@@ -4,7 +4,7 @@
 import numpy as np
 from pygsti.modelpacks import smq1Q_XYI
 from pygsti.baseobjs import QubitSpace, Basis
-from pygsti.modelmembers.operations import StochasticNoiseOp
+from pygsti.modelmembers.operations import StochasticNoiseOp, DepolarizeOp
 from pygsti.circuits import Circuit
 from pygsti.models import create_explicit_model
 from pygsti.modelmembers.operations.composedop import ComposedOp
@@ -79,6 +79,18 @@ def test_dense_op(self):
         kkdag = [kop @ kop.conjugate().T for kop in op.kraus_operators]
         assert(np.allclose(sum(kkdag), np.identity(2)))
 
+    def test_kraus_ops(self):
+        # test that kraus operators for a depolarization op can recover that depolarization op
+        op = DepolarizeOp(QubitSpace(1), initial_rate=0.1)
+        kraus_ops = op.kraus_operators
+        test = FullArbitraryOp.from_kraus_operators(kraus_ops)
+        assert np.allclose(op.to_dense(), test.to_dense())
+
+        op = FullArbitraryOp(op.to_dense(), 'pp')
+        kraus_ops = op.kraus_operators
+        test = FullArbitraryOp.from_kraus_operators(kraus_ops)
+        assert np.allclose(op.to_dense(), test.to_dense())
+
     def test_stochastic_errorgen_equivalence_single(self):
         #Check that StochasticOp and 'S'-type elementary errorgen give the same op
         B = Basis.cast('PP', 4)