diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py
index ccdeea4d8..620f07105 100644
--- a/python/damask/mechanics.py
+++ b/python/damask/mechanics.py
@@ -337,7 +337,7 @@ def _polar_decomposition(T: _np.ndarray,
     if isinstance(requested, str): requested = [requested]
 
     u, _, vh = _np.linalg.svd(T)
-    R = _np.einsum('...ij,...jk',u,vh)
+    R = u @ vh
 
     output = []
     if 'R' in requested:
diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py
index 825d88463..a53472c8b 100644
--- a/python/tests/test_mechanics.py
+++ b/python/tests/test_mechanics.py
@@ -158,7 +158,7 @@ class TestMechanics:
     @pytest.mark.parametrize('side',[('left','V'),('right','U')])
     def test_polar_decomposition(self,side):
         F     = np.random.rand(self.n,3,3)
-        F     = np.einsum('...ij,...jk',F,F) # positive determinant
+        F     = F @ F # positive determinant
         F_vec = np.reshape(F,(self.n//10,10,3,3))
         p = mechanics._polar_decomposition(F_vec,side[1])
         for p_,F_ in zip(np.reshape(p,F.shape),F):