add DensityMatrix.to_statevector() for pure states (#4433)

* add DensityMatrix.to_statevector()

* lint

* updates

* review

* nit

* docs

* ..

Co-authored-by: Christopher J. Wood <cjwood@us.ibm.com>

qiskit/quantum_info/states/densitymatrix.py

@@ -585,6 +585,38 @@ class DensityMatrix(QuantumState):
         vec._append_instruction(obj, qargs=qargs)
         return vec
 
+    def to_statevector(self, atol=None, rtol=None):
+        """Return a statevector from a pure density matrix.
+
+        Args:
+            atol (float): Absolute tolerance for checking operation validity.
+            rtol (float): Relative tolerance for checking operation validity.
+
+        Returns:
+            Statevector: The pure density matrix's corresponding statevector.
+                Corresponds to the eigenvector of the only non-zero eigenvalue.
+
+        Raises:
+            QiskitError: if the state is not pure.
+        """
+        if atol is None:
+            atol = self.atol
+        if rtol is None:
+            rtol = self.rtol
+
+        if not is_hermitian_matrix(self._data, atol=atol, rtol=rtol):
+            raise QiskitError("Not a valid density matrix (non-hermitian).")
+
+        evals, evecs = np.linalg.eig(self._data)
+
+        nonzero_evals = evals[abs(evals) > atol]
+        if len(nonzero_evals) != 1 or not np.isclose(nonzero_evals[0], 1,
+                                                     atol=atol, rtol=rtol):
+            raise QiskitError("Density matrix is not a pure state")
+
+        psi = evecs[:, np.argmax(evals)]  # eigenvectors returned in columns.
+        return Statevector(psi)
+
     def to_counts(self):
         """Returns the density matrix as a counts dict of probabilities.
 

test/python/quantum_info/states/test_densitymatrix.py

@@ -359,6 +359,22 @@ class TestDensityMatrix(QiskitTestCase):
                     target[key] = 2 * j + i + 1
             self.assertDictAlmostEqual(target, vec.to_dict())
 
+    def test_densitymatrix_to_statevector_pure(self):
+        """Test converting a pure density matrix to statevector."""
+        state = 1/np.sqrt(2) * (np.array([1, 0, 0, 0, 0, 0, 0, 1]))
+        psi = Statevector(state)
+        rho = DensityMatrix(psi)
+        phi = rho.to_statevector()
+        self.assertTrue(psi.equiv(phi))
+
+    def test_densitymatrix_to_statevector_mixed(self):
+        """Test converting a pure density matrix to statevector."""
+        state_1 = 1/np.sqrt(2) * (np.array([1, 0, 0, 0, 0, 0, 0, 1]))
+        state_2 = 1/np.sqrt(2) * (np.array([0, 0, 0, 0, 0, 0, 1, 1]))
+        psi = 0.5 * (Statevector(state_1) + Statevector(state_2))
+        rho = DensityMatrix(psi)
+        self.assertRaises(QiskitError, rho.to_statevector)
+
     def test_probabilities_product(self):
         """Test probabilities method for product state"""