Repository: qutip/qutip (python)

Issue:
Add `einsum` for `Qobj`.
### Problem Description

Having a version of `einsum` working on `Qobj`'s dimensions would be useful (qutip/qutip-qip#225).
A dense implementation should be easy enough to develop using `qutip.core.dimension.to_tensor_rep`.


### Proposed Solution

Use `to_tensor_rep`, `from_tensor_rep` with `np.einsum` to create qutip version of einsum.

### Alternate Solutions

If would be even better to develop it without conversion to numpy array to allow support for sparse, jax, cupy, etc. matrices.
However this would probably not be an easy task doable in the scope of `good first issue`.

### Additional Context

_No response_

Pull request description:
Add einsum
**Description**
Added a einsum, a numpy.einsum (Evaluates the Einstein summation convention on the operands.) version for Qobj dimensions.

**Related issues or PRs**
fixes #2340

** Question: Not sure if I should update the documentation and how. Just in case I did not modify the docs yet.

Patch:
diff --git a/doc/changes/2545.feature b/doc/changes/2545.feature
new file mode 100644
index 000000000..3514cff82
--- /dev/null
+++ b/doc/changes/2545.feature
@@ -0,0 +1,1 @@
+Implements a numpy.einsum version for Qobj dimensions (Evaluates the Einstein summation convention on the operands.)
\ No newline at end of file
diff --git a/qutip/core/dimensions.py b/qutip/core/dimensions.py
index f00afc1ff..56fe35396 100644
--- a/qutip/core/dimensions.py
+++ b/qutip/core/dimensions.py
@@ -352,6 +352,34 @@ def _frozen(*args, **kwargs):
     raise RuntimeError("Dimension cannot be modified.")
 
 
+def einsum(subscripts, *operands):
+    """
+    Implementation of numpy.einsum for Qobj.
+    Evaluates the Einstein summation convention on the operands.
+    Parameters
+    ----------
+    subscripts: str
+        Specifies the subscripts for summation as comma
+        separated list of subscript labels.
+    operands: list of array_like
+        These are the arrays for the operation.
+
+    Returns
+    -------
+    Qobj (numpy.complex128)
+        Result of einsum as Qobj (numpy.complex128 if result is scalar)
+    """
+    operands_array = [to_tensor_rep(op) for op in operands]
+    result = np.einsum(subscripts, *operands_array)
+    if result.shape == ():
+        return result
+    dims = [
+        [d for d in result.shape[:result.ndim // 2]],
+        [d for d in result.shape[result.ndim // 2:]]
+    ]
+    return from_tensor_rep(result, dims)
+
+
 class MetaSpace(type):
     def __call__(cls, *args: SpaceLike, rep: str = None) -> "Space":
         """


Test patch:
diff --git a/qutip/tests/core/test_dimensions.py b/qutip/tests/core/test_dimensions.py
index 7be9d2760..96b2ed4cd 100644
--- a/qutip/tests/core/test_dimensions.py
+++ b/qutip/tests/core/test_dimensions.py
@@ -6,7 +6,7 @@ import qutip
 from qutip.core.dimensions import (
     flatten, unflatten, enumerate_flat, deep_remove, deep_map,
     dims_idxs_to_tensor_idxs, dims_to_tensor_shape, dims_to_tensor_perm,
-    collapse_dims_super, collapse_dims_oper, Dimensions
+    collapse_dims_super, collapse_dims_oper, einsum, Dimensions
 )
 
 
@@ -189,3 +189,17 @@ def test_dims_comparison():
     assert Dimensions([[1], [2]])[0] != Dimensions([[1], [2]])[1]
     assert not Dimensions([[1], [2]])[1] != Dimensions([[1], [2]])[1]
     assert not Dimensions([[1], [2]])[0] != Dimensions([[1], [2]])[0]
+
+
+@pytest.mark.parametrize(["subscripts", "operands", "expected"], [
+    pytest.param("ii", [qutip.sigmaz()], 0),
+    pytest.param("ij", [qutip.sigmax()], qutip.sigmax()),
+    pytest.param("ij->ji", [qutip.sigmay()], qutip.sigmay().trans()),
+    pytest.param("ij,ji", [qutip.sigmaz(), qutip.sigmaz()], 2),
+    pytest.param("ijij", [qutip.tensor(qutip.thermal_dm(2,1), qutip.thermal_dm(2,1))], 1),
+    pytest.param("ikjl,jm->ikml", [qutip.tensor(qutip.sigmaz(), qutip.sigmaz()), 
+                             qutip.sigmaz()], qutip.tensor(qutip.qeye(2), qutip.sigmaz())),
+    pytest.param("ijkl->kjil", [qutip.tensor(qutip.sigmam(), qutip.sigmaz())], qutip.tensor(qutip.sigmap(), qutip.sigmaz()))
+])
+def test_einsum(subscripts, operands, expected):
+    assert einsum(subscripts, *operands) == expected
\ No newline at end of file