Stavros' suggestion

src/qibojit/backends/gpu.py

@@ -1,7 +1,7 @@
 from typing import Union
 
 import numpy as np
-from qibo.backends.numpy import NumpyBackend
+from qibo.backends.numpy import NumpyBackend, _calculate_negative_power_singular_matrix
 from qibo.config import log, raise_error
 
 from qibojit.backends.cpu import NumbaBackend
@@ -550,12 +550,8 @@ def calculate_matrix_power(
  # negative powers of singular matrices via SVD
  determinant = self.cp.linalg.det(matrix)
  if abs(determinant) < precision_singularity:
- U, S, Vh = self.cp.linalg.svd(matrix)
- S_inv = self.cp.where(
- self.cp.abs(S) < precision_singularity, 0.0, S**power
- )
- return (
- self.cp.linalg.inv(Vh) @ self.cp.diag(S_inv) @ self.cp.linalg.inv(U)
+ return _calculate_negative_power_singular_matrix(
+ matrix, power, precision_singularity, self.cp, self
  )
 
  copied = self.to_numpy(matrix)