fix simdiag not returning orthonormal eigenvectors

doc/changes/2269.bugfix

@@ -0,0 +1 @@
+Fixed simdiag not returning orthonormal eigenvectors.

qutip/simdiag.py

@@ -22,7 +22,7 @@ def _degen(tol, vecs, ops, i=0):
                               / (1 - np.abs(dot)**2)**0.5)
 
     subspace = vecs.conj().T @ ops[i].full() @ vecs
-    eigvals, eigvecs = la.eig(subspace)
+    eigvals, eigvecs = la.eigh(subspace)
 
     perm = np.argsort(eigvals)
     eigvals = eigvals[perm]

qutip/tests/test_simdiag.py

@@ -73,6 +73,26 @@ def test_simdiag_degen_large():
                 atol=1e-13
             )
 
+def test_simdiag_orthonormal_eigenvectors():
+
+    # Special matrix that used to be problematic (see Issue #2268)
+    a = np.array([[1,  0, 1, -1, 0],
+                  [0,  4, 0,  0, 1],
+                  [1,  0, 4,  1, 0],
+                  [-1, 0, 1,  4, 0],
+                  [0,  1, 0,  0, 4]])
+
+    _, evecs = qutip.simdiag([qutip.Qobj(a), qutip.qeye(5)])
+    evecs = np.array([evec.full() for evec in evecs]).squeeze()
+
+    # Check that eigenvectors form an othonormal basis
+    # (<=> matrix of eigenvectors is unitary)
+    np.testing.assert_allclose(
+        evecs@evecs.conj().T,
+        np.eye(len(evecs)),
+        atol=1e-13
+    )
+
 
 def test_simdiag_no_input():
     with pytest.raises(ValueError):