Autodifferentiation support for PytorchBackend

src/qibo/backends/npmatrices.py

@@ -1,3 +1,5 @@
+import cmath
+import math
 from functools import cached_property
 
 from qibo.config import raise_error
@@ -15,9 +17,13 @@ def __init__(self, dtype):
     def _cast(self, x, dtype):
         return self.np.array(x, dtype=dtype)
 
+    # This method is used to cast the parameters of the gates to the right type for other backends
+    def _cast_parameter(self, x):
+        return x
+
     @cached_property
     def H(self):
-        return self._cast([[1, 1], [1, -1]], dtype=self.dtype) / self.np.sqrt(2)
+        return self._cast([[1, 1], [1, -1]], dtype=self.dtype) / math.sqrt(2)
 
     @cached_property
     def X(self):
@@ -50,19 +56,17 @@ def SDG(self):
     @cached_property
     def T(self):
         return self._cast(
-            [[1, 0], [0, self.np.exp(1j * self.np.pi / 4.0)]], dtype=self.dtype
+            [[1 + 0j, 0], [0, cmath.exp(1j * math.pi / 4.0)]], dtype=self.dtype
         )
 
     @cached_property
     def TDG(self):
         return self._cast(
-            [[1, 0], [0, self.np.exp(-1j * self.np.pi / 4.0)]], dtype=self.dtype
+            [[1 + 0j, 0], [0, cmath.exp(-1j * math.pi / 4.0)]], dtype=self.dtype
         )
 
     def I(self, n=2):
-        # dtype=complex is necessary for pytorch backend,
-        # _cast will take care of casting in the right dtype for all the backends
-        return self._cast(self.np.eye(n, dtype=complex), dtype=self.dtype)
+        return self._cast(self.np.eye(n), dtype=self.dtype)
 
     def Align(self, delay, n=2):
         return self._cast(self.I(n), dtype=self.dtype)
@@ -71,20 +75,25 @@ def M(self):  # pragma: no cover
         raise_error(NotImplementedError)
 
     def RX(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast([[cos, isin], [isin, cos]], dtype=self.dtype)
 
     def RY(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         sin = self.np.sin(theta / 2.0) + 0j
         return self._cast([[cos, -sin], [sin, cos]], dtype=self.dtype)
 
     def RZ(self, theta):
+        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         return self._cast([[self.np.conj(phase), 0], [0, phase]], dtype=self.dtype)
 
     def PRX(self, theta, phi):
+        theta = self._cast_parameter(theta)
+        phi = self._cast_parameter(phi)
         cos = self.np.cos(theta / 2)
         sin = self.np.sin(theta / 2)
         exponent1 = -1.0j * self.np.exp(-1.0j * phi)
@@ -95,29 +104,36 @@ def PRX(self, theta, phi):
         )
 
     def GPI(self, phi):
+        phi = self._cast_parameter(phi)
         phase = self.np.exp(1.0j * phi)
         return self._cast([[0, self.np.conj(phase)], [phase, 0]], dtype=self.dtype)
 
     def GPI2(self, phi):
+        phi = self._cast_parameter(phi)
         phase = self.np.exp(1.0j * phi)
         return self._cast(
             [[1, -1.0j * self.np.conj(phase)], [-1.0j * phase, 1]], dtype=self.dtype
-        ) / self.np.sqrt(2)
+        ) / math.sqrt(2)
 
     def U1(self, theta):
+        theta = self._cast_parameter(theta)
         phase = self.np.exp(1j * theta)
         return self._cast([[1, 0], [0, phase]], dtype=self.dtype)
 
     def U2(self, phi, lam):
+        phi = self._cast_parameter(phi)
+        lam = self._cast_parameter(lam)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
         eminus = self.np.exp(1j * (phi - lam) / 2.0)
         return self._cast(
-            [[self.np.conj(eplus), -self.np.conj(eminus)], [eminus, eplus]]
-            / self.np.sqrt(2),
+            [[self.np.conj(eplus), -self.np.conj(eminus)], [eminus, eplus]],
             dtype=self.dtype,
-        )
+        ) / math.sqrt(2)
 
     def U3(self, theta, phi, lam):
+        theta = self._cast_parameter(theta)
+        phi = self._cast_parameter(phi)
+        lam = self._cast_parameter(lam)
         cost = self.np.cos(theta / 2)
         sint = self.np.sin(theta / 2)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
@@ -131,8 +147,10 @@ def U3(self, theta, phi, lam):
         )
 
     def U1q(self, theta, phi):
+        theta = self._cast_parameter(theta)
+        phi = self._cast_parameter(phi)
         return self._cast(
-            self.U3(theta, phi - self.np.pi / 2, self.np.pi / 2 - phi), dtype=self.dtype
+            self.U3(theta, phi - math.pi / 2, math.pi / 2 - phi), dtype=self.dtype
         )
 
     @cached_property
@@ -161,7 +179,7 @@ def CZ(self):
 
     @cached_property
     def CSX(self):
-        a = (1 + 1j) / 2
+        a = self._cast_parameter((1 + 1j) / 2)
         b = self.np.conj(a)
         return self._cast(
             [
@@ -175,7 +193,7 @@ def CSX(self):
 
     @cached_property
     def CSXDG(self):
-        a = (1 - 1j) / 2
+        a = self._cast_parameter((1 - 1j) / 2)
         b = self.np.conj(a)
         return self._cast(
             [
@@ -188,6 +206,7 @@ def CSXDG(self):
         )
 
     def CRX(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         matrix = [
@@ -199,12 +218,14 @@ def CRX(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CRY(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         sin = self.np.sin(theta / 2.0) + 0j
         matrix = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, cos, -sin], [0, 0, sin, cos]]
         return self._cast(matrix, dtype=self.dtype)
 
     def CRZ(self, theta):
+        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         matrix = [
             [1, 0, 0, 0],
@@ -215,6 +236,7 @@ def CRZ(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU1(self, theta):
+        theta = self._cast_parameter(theta)
         phase = self.np.exp(1j * theta)
         matrix = [
             [1, 0, 0, 0],
@@ -225,8 +247,10 @@ def CU1(self, theta):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU2(self, phi, lam):
-        eplus = self.np.exp(1j * (phi + lam) / 2.0) / self.np.sqrt(2)
-        eminus = self.np.exp(1j * (phi - lam) / 2.0) / self.np.sqrt(2)
+        phi = self._cast_parameter(phi)
+        lam = self._cast_parameter(lam)
+        eplus = self.np.exp(1j * (phi + lam) / 2.0) / math.sqrt(2)
+        eminus = self.np.exp(1j * (phi - lam) / 2.0) / math.sqrt(2)
         matrix = [
             [1, 0, 0, 0],
             [0, 1, 0, 0],
@@ -236,6 +260,9 @@ def CU2(self, phi, lam):
         return self._cast(matrix, dtype=self.dtype)
 
     def CU3(self, theta, phi, lam):
+        theta = self._cast_parameter(theta)
+        phi = self._cast_parameter(phi)
+        lam = self._cast_parameter(lam)
         cost = self.np.cos(theta / 2)
         sint = self.np.sin(theta / 2)
         eplus = self.np.exp(1j * (phi + lam) / 2.0)
@@ -271,8 +298,8 @@ def SiSWAP(self):
         return self._cast(
             [
                 [1 + 0j, 0j, 0j, 0j],
-                [0j, 1 / self.np.sqrt(2) + 0j, 1j / self.np.sqrt(2), 0j],
-                [0j, 1j / self.np.sqrt(2), 1 / self.np.sqrt(2) + 0j, 0j],
+                [0j, 1 / math.sqrt(2) + 0j, 1j / math.sqrt(2), 0j],
+                [0j, 1j / math.sqrt(2), 1 / math.sqrt(2) + 0j, 0j],
                 [0j, 0j, 0j, 1 + 0j],
             ],
             dtype=self.dtype,
@@ -283,8 +310,8 @@ def SiSWAPDG(self):
         return self._cast(
             [
                 [1 + 0j, 0j, 0j, 0j],
-                [0j, 1 / self.np.sqrt(2) + 0j, -1j / self.np.sqrt(2), 0j],
-                [0j, -1j / self.np.sqrt(2), 1 / self.np.sqrt(2) + 0j, 0j],
+                [0j, 1 / math.sqrt(2) + 0j, -1j / math.sqrt(2), 0j],
+                [0j, -1j / math.sqrt(2), 1 / math.sqrt(2) + 0j, 0j],
                 [0j, 0j, 0j, 1 + 0j],
             ],
             dtype=self.dtype,
@@ -297,6 +324,8 @@ def FSWAP(self):
         )
 
     def fSim(self, theta, phi):
+        theta = self._cast_parameter(theta)
+        phi = self._cast_parameter(phi)
         cost = self.np.cos(theta) + 0j
         isint = -1j * self.np.sin(theta)
         phase = self.np.exp(-1j * phi)
@@ -312,12 +341,12 @@ def fSim(self, theta, phi):
 
     @cached_property
     def SYC(self):
-        cost = self.np.cos(self.np.pi / 2) + 0j
-        isint = -1j * self.np.sin(self.np.pi / 2)
-        phase = self.np.exp(-1j * self.np.pi / 6)
+        cost = math.cos(math.pi / 2) + 0j
+        isint = -1j * math.sin(math.pi / 2)
+        phase = cmath.exp(-1j * math.pi / 6)
         return self._cast(
             [
-                [1, 0, 0, 0],
+                [1 + 0j, 0, 0, 0],
                 [0, cost, isint, 0],
                 [0, isint, cost, 0],
                 [0, 0, 0, phase],
@@ -326,6 +355,7 @@ def SYC(self):
         )
 
     def GeneralizedfSim(self, u, phi):
+        phi = self._cast_parameter(phi)
         phase = self.np.exp(-1j * phi)
         return self._cast(
             [
@@ -338,6 +368,7 @@ def GeneralizedfSim(self, u, phi):
         )
 
     def RXX(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast(
@@ -351,6 +382,7 @@ def RXX(self, theta):
         )
 
     def RYY(self, theta):
+        theta = self._cast_parameter(theta)
         cos = self.np.cos(theta / 2.0) + 0j
         isin = -1j * self.np.sin(theta / 2.0)
         return self._cast(
@@ -364,6 +396,7 @@ def RYY(self, theta):
         )
 
     def RZZ(self, theta):
+        theta = self._cast_parameter(theta)
         phase = self.np.exp(0.5j * theta)
         return self._cast(
             [
@@ -376,6 +409,7 @@ def RZZ(self, theta):
         )
 
     def RZX(self, theta):
+        theta = self._cast_parameter(theta)
         cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j
         return self._cast(
             [
@@ -388,6 +422,7 @@ def RZX(self, theta):
         )
 
     def RXXYY(self, theta):
+        theta = self._cast_parameter(theta)
         cos, sin = self.np.cos(theta / 2) + 0j, self.np.sin(theta / 2) + 0j
         return self._cast(
             [
@@ -400,6 +435,11 @@ def RXXYY(self, theta):
         )
 
     def MS(self, phi0, phi1, theta):
+        phi0, phi1, theta = (
+            self._cast_parameter(phi0),
+            self._cast_parameter(phi1),
+            self._cast_parameter(theta),
+        )
         plus = self.np.exp(1.0j * (phi0 + phi1))
         minus = self.np.exp(1.0j * (phi0 - phi1))
         cos = self.np.cos(theta / 2) + 0j
@@ -415,6 +455,7 @@ def MS(self, phi0, phi1, theta):
         )
 
     def GIVENS(self, theta):
+        theta = self._cast_parameter(theta)
         return self._cast(
             [
                 [1, 0, 0, 0],
@@ -438,7 +479,7 @@ def ECR(self):
                 [-1j, 1 + 0j, 0j, 0j],
             ],
             dtype=self.dtype,
-        ) / self.np.sqrt(2)
+        ) / math.sqrt(2)
 
     @cached_property
     def TOFFOLI(self):
@@ -473,6 +514,7 @@ def CCZ(self):
         )
 
     def DEUTSCH(self, theta):
+        theta = self._cast_parameter(theta)
         sin = self.np.sin(theta) + 0j  # 0j necessary for right tensorflow dtype
         cos = self.np.cos(theta) + 0j
         return self._cast(

src/qibo/backends/numpy.py

@@ -156,7 +156,6 @@ def matrix_fused(self, fgate):
         return self.cast(matrix.toarray())
 
     def apply_gate(self, gate, state, nqubits):
-        state = self.cast(state)
         state = self.np.reshape(state, nqubits * (2,))
         matrix = gate.matrix(self)
         if gate.is_controlled_by:
@@ -718,31 +717,37 @@ def calculate_norm_density_matrix(self, state, order="nuc"):
         return self.np.linalg.norm(state, ord=order)
 
     def calculate_overlap(self, state1, state2):
-        return self.np.abs(self.np.sum(np.conj(self.cast(state1)) * self.cast(state2)))
+        return self.np.abs(
+            self.np.sum(self.np.conj(self.cast(state1)) * self.cast(state2))
+        )
 
     def calculate_overlap_density_matrix(self, state1, state2):
         return self.np.trace(
             self.np.matmul(self.np.conj(self.cast(state1)).T, self.cast(state2))
         )
 
-    def calculate_eigenvalues(self, matrix, k=6):
+    def calculate_eigenvalues(self, matrix, k=6, hermitian=True):
         if self.issparse(matrix):
             log.warning(
                 "Calculating sparse matrix eigenvectors because "
                 "sparse modules do not provide ``eigvals`` method."
             )
             return self.calculate_eigenvectors(matrix, k=k)[0]
-        return np.linalg.eigvalsh(matrix)
+        if hermitian:
+            return np.linalg.eigvalsh(matrix)
+        return np.linalg.eigvals(matrix)
 
-    def calculate_eigenvectors(self, matrix, k=6):
+    def calculate_eigenvectors(self, matrix, k=6, hermitian=True):
         if self.issparse(matrix):
             if k < matrix.shape[0]:
                 from scipy.sparse.linalg import eigsh
 
                 return eigsh(matrix, k=k, which="SA")
             else:  # pragma: no cover
                 matrix = self.to_numpy(matrix)
-        return np.linalg.eigh(matrix)
+        if hermitian:
+            return np.linalg.eigh(matrix)
+        return np.linalg.eig(matrix)
 
     def calculate_matrix_exp(self, a, matrix, eigenvectors=None, eigenvalues=None):
         if eigenvectors is None or self.issparse(matrix):