Revert "added full diagonal gradients and a tutorial file"

This reverts commit 93eb20a.

src/qibo/models/dbi/double_bracket.py

@@ -82,6 +82,7 @@ def __init__(
         self.mode = mode
         self.scheduling = scheduling
         self.cost = cost
+        self.cost_str = cost.name
         self.state = state
 
     def __call__(
@@ -99,7 +100,7 @@ def __call__(
             if d is None:
                 d = self.diagonal_h_matrix
             operator = self.backend.calculate_matrix_exp(
-                1.0j*step,
+                1.0j * step,
                 self.commutator(d, self.h.matrix),
             )
         elif mode is DoubleBracketGeneratorType.group_commutator:
@@ -114,7 +115,6 @@ def __call__(
         operator_dagger = self.backend.cast(
             np.matrix(self.backend.to_numpy(operator)).getH()
         )
-
         self.h.matrix = operator @ self.h.matrix @ operator_dagger
 
     @staticmethod
@@ -140,17 +140,20 @@ def off_diagonal_norm(self):
         return np.sqrt(
             np.real(np.trace(self.backend.to_numpy(off_diag_h_dag @ self.off_diag_h)))
         )
+    @property
+    def least_squares(self,d: np.array):
+        """Least squares cost function."""
+        H = self.backend.cast(
+            np.matrix(self.backend.to_numpy(self.h)).getH()
+        )
+        D = d
+        return -(np.linalg.trace(H@D)-0.5(np.linalg.norm(H)**2+np.linalg.norm(D)**2))
 
     @property
     def backend(self):
         """Get Hamiltonian's backend."""
         return self.h0.backend
 
-    def least_squares(self, D: np.array):
-        """Least squares cost function."""
-        H = self.h.matrix
-        return -np.real(np.trace(H@D)-0.5*(np.linalg.norm(H)**2+np.linalg.norm(D)**2))
-
     def choose_step(
         self,
         d: Optional[np.array] = None,
@@ -189,7 +192,7 @@ def loss(self, step: float, d: np.array = None, look_ahead: int = 1):
         if self.cost == DoubleBracketCostFunction.off_diagonal_norm:
             loss = self.off_diagonal_norm
         elif self.cost == DoubleBracketCostFunction.least_squares:
-            loss = self.least_squares(d)
+            loss = self.least_squares(d=d)
         else:
             loss = self.energy_fluctuation(self.state)
 
@@ -210,9 +213,7 @@ def energy_fluctuation(self, state):
         Args:
             state (np.ndarray): quantum state to be used to compute the energy fluctuation with H.
         """
-        state_vector = np.zeros(len(self.h.matrix))
-        state_vector[state] = 1.0
-        return np.real(self.h.energy_fluctuation(state_vector))
+        return self.h.energy_fluctuation(state)
 
     def sigma(self, h: np.array):
         return h - self.backend.cast(np.diag(np.diag(self.backend.to_numpy(h))))

src/qibo/models/dbi/utils_scheduling.py

@@ -1,26 +1,21 @@
 import math
 from functools import partial
 from typing import Optional
-from copy import deepcopy
+
 import hyperopt
 import numpy as np
 
-
 error = 1e-3
 
-def commutator(A, B):
-    """Compute commutator between two arrays."""
-    return A@B-B@A
-
 def variance(A, state):
     """Calculates the variance of a matrix A with respect to a state: Var($A$) = $\langle\mu|A^2|\mu\rangle-\langle\mu|A|\mu\rangle^2$"""
     B = A@A
     return B[state,state]-A[state,state]**2
 
 def covariance(A, B, state):
     """Calculates the covariance of two matrices A and B with respect to a state: Cov($A,B$) = $\langle\mu|AB|\mu\rangle-\langle\mu|A|\mu\rangle\langle\mu|B|\mu\rangle$"""
-    C = A@B+B@A
-    return C[state,state]-2*A[state,state]*B[state,state]
+    C = A@B
+    return C[state,state]-A[state,state]*B[state,state]
 
 def grid_search_step(
     dbi_object,
@@ -49,7 +44,6 @@ def grid_search_step(
         d = dbi_object.diagonal_h_matrix
 
     loss_list = [dbi_object.loss(step, d=d) for step in space]
-
     idx_max_loss = np.argmin(loss_list)
     return space[idx_max_loss]
 
@@ -89,14 +83,12 @@ def hyperopt_step(
         d = dbi_object.diagonal_h_matrix
 
     space = space("step", step_min, step_max)
-
-
     best = hyperopt.fmin(
-    fn=partial(dbi_object.loss, d=d, look_ahead=look_ahead),
-    space=space,
-    algo=optimizer.suggest,
-    max_evals=max_evals,
-    verbose=verbose,
+        fn=partial(dbi_object.loss, d=d, look_ahead=look_ahead),
+        space=space,
+        algo=optimizer.suggest,
+        max_evals=max_evals,
+        verbose=verbose,
     )
     return best["step"]
 
@@ -120,7 +112,6 @@ def polynomial_step(
     """
     if cost is None:
         cost = dbi_object.cost.name
-
     if d is None:
         d = dbi_object.diagonal_h_matrix
 
@@ -144,11 +135,7 @@ def polynomial_step(
     ]
     # solution exists, return minimum s
     if len(real_positive_roots) > 0:
-        sol = min(real_positive_roots)
-        for s in real_positive_roots:
-            if dbi_object.loss(s, d) < dbi_object.loss(sol, d):
-                sol = s
-        return sol
+        return min(real_positive_roots)
     # solution does not exist, return None
     else:
         return None
@@ -180,81 +167,25 @@ def least_squares_polynomial_expansion_coef(dbi_object, d, n):
     if d is None:
         d = dbi_object.diagonal_h_matrix
     # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H
-    Gamma_list = dbi_object.generate_Gamma_list(n+1, d)
+    W = dbi_object.commutator(d, dbi_object.sigma(dbi_object.h.matrix))
+    Gamma_list = dbi_object.generate_Gamma_list(n, d)
     exp_list = np.array([1 / math.factorial(k) for k in range(n + 1)])
     # coefficients
     coef = np.empty(n)
     for i in range(n):
-        coef[i] = np.real(exp_list[i]*np.trace(d@Gamma_list[i+1]))
-    coef = list(reversed(coef))
+        coef[i] = exp_list[i]*np.trace(d@Gamma_list[i+1])
+
     return coef
 
 #TODO: add a general expansion formula not stopping at 3rd order
 def energy_fluctuation_polynomial_expansion_coef(dbi_object, d, n, state):
     if d is None:
         d = dbi_object.diagonal_h_matrix
     # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H
-    Gamma_list = dbi_object.generate_Gamma_list(n+1, d)
+    Gamma_list = dbi_object.generate_Gamma_list(n, d)
     # coefficients
     coef = np.empty(3)
-    coef[0] = np.real(2*covariance(Gamma_list[0], Gamma_list[1],state))
-    coef[1] = np.real(2*variance(Gamma_list[1],state))
-    coef[2] = np.real(covariance(Gamma_list[0], Gamma_list[3],state)+3*covariance(Gamma_list[1], Gamma_list[2],state))
-    coef = list(reversed(coef))
+    coef[0] = 2*covariance(Gamma_list[0], Gamma_list[1],state)
+    coef[1] = 2*variance(Gamma_list[1],state)
+    coef[2] = covariance(Gamma_list[0], Gamma_list[3],state)+3*covariance(Gamma_list[1], Gamma_list[2],state)
     return coef
-
-def dGamma_diDiagonal(dbi_object, d, H, n, i,dGamma, Gamma_list):
-    # Derivative of gamma with respect to diagonal elements of D (full-diagonal ansatz)
-    A = np.zeros(d.shape)
-    A[i,i] = 1
-    B = commutator(commutator(A,H),Gamma_list[n-1])
-    W = commutator(d,H)
-    return B + commutator(W,dGamma[-1])
-
-def dpolynomial_diDiagonal(dbi_object, d,H,i):
-    # Derivative of polynomial approximation of potential function with respect to diagonal elements of d (full-diagonal ansatz)
-    # Formula can be expanded easily to any order, with n=3 corresponding to cubic approximation
-    derivative = 0
-    s = polynomial_step(dbi_object, n=3, d=d)
-    A = np.zeros(d.shape)
-    Gamma_list = dbi_object.generate_Gamma_list(4, d)
-    A[i,i] = 1
-    dGamma = [commutator(A,H)]
-    derivative += np.real(np.trace(Gamma_list[0]@A)+np.trace(dGamma[0]@d+Gamma_list[1]@A)*s)
-    for n in range(2,4):
-        dGamma.append(dGamma_diDiagonal(dbi_object,d,H,n,i,dGamma,Gamma_list))
-        derivative += np.real(np.trace(dGamma[-1]@d + Gamma_list[n]@A)*s**n/math.factorial(n))
-
-    return derivative
-
-def gradientDiagonal(dbi_object,d,H):
-    # Gradient of potential function with respect to diagonal elements of D (full-diagonal ansatz)
-    grad = np.zeros(len(d))
-    for i in range(len(d)):
-        derivative = dpolynomial_diDiagonal(dbi_object,d,H,i)
-        grad[i] = d[i,i]-derivative
-    return grad
-
-def gradient_ascent(dbi_object, d, step, iterations):
-    H = dbi_object.h.matrix
-    loss = np.zeros(iterations+1)
-    grad = np.zeros((iterations,len(d)))
-    dbi_new = deepcopy(dbi_object)
-    s = polynomial_step(dbi_object, n = 3, d=d)
-    dbi_new(s,d=d)
-    loss[0] = dbi_new(d)
-    diagonals = np.empty((len(d),iterations+1))
-    diagonals[:,0] = np.diag(d)
-
-    for i in range(iterations):
-        dbi_new = deepcopy(dbi_object)
-        grad[i,:] = gradientDiagonal(dbi_object, d, H)
-        for j in range(len(d)):
-            d[j,j] = d[j,j] - step*grad[i,j] 
-        s = polynomial_step(dbi_object, n = 3, d=d)
-        dbi_new(s,d=d)
-        loss[i+1] = dbi_new.least_squares(d)
-        diagonals[:,i+1] = np.diag(d)
-
-
-    return d,loss,grad,diagonals