Source code for qibo.models.dbi.double_bracket

from copy import copy
from enum import Enum, auto
from typing import Optional

import numpy as np
import optuna

from qibo.config import raise_error
from qibo.hamiltonians import Hamiltonian
from qibo.models.dbi.utils import *
from qibo.models.dbi.utils_scheduling import (
    grid_search_step,
    hyperopt_step,
    polynomial_step,
    simulated_annealing_step,
)
from qibo.quantum_info.linalg_operations import commutator, matrix_exponentiation


[docs]class DoubleBracketGeneratorType(Enum): """Define DBF evolution.""" canonical = auto() """Use canonical commutator.""" single_commutator = auto() """Use single commutator.""" group_commutator = auto() """Use group commutator approximation""" group_commutator_third_order = auto() """Implements Eq. (8) of Ref. [1], i.e. .. math:: e^{\\frac{\\sqrt{5}-1}{2}isH} \\, e^{\\frac{\\sqrt{5}-1}{2}isD} \\, e^{-isH} \\, e^{isD} \\, e^{\\frac{3-\\sqrt{5}}{2}isH} \\, e^{isD} \\approx e^{-s^2[H,D]} + O(s^4) \\, . :math:`s` must be taken as :math:`\\sqrt{s}` to approximate the flow using the commutator. """
class DoubleBracketCostFunction(str, Enum): """Define the DBI cost function.""" off_diagonal_norm = "off_diagonal_norm" """Use off-diagonal norm as cost function.""" least_squares = "least_squares" """Use least squares as cost function.""" energy_fluctuation = "energy_fluctuation" """Use energy fluctuation as cost function.""" class DoubleBracketScheduling(Enum): """Define the DBI scheduling strategies.""" hyperopt = hyperopt_step """Use optuna package to hyperoptimize the DBI step.""" grid_search = grid_search_step """Use greedy grid search.""" polynomial_approximation = polynomial_step """Use polynomial expansion (analytical) of the loss function.""" simulated_annealing = simulated_annealing_step """Use simulated annealing algorithm"""
[docs]class DoubleBracketIteration: """ Class implementing the Double Bracket iteration algorithm. For more details, see Ref. [1]. Args: hamiltonian (:class:`qibo.hamiltonians.Hamiltonian`): starting Hamiltonian. mode (:class:`qibo.models.dbi.double_bracket.DoubleBracketGeneratorType`): type of generator of the evolution. scheduling (:class:`qibo.models.dbi.double_bracket.DoubleBracketScheduling`): type of scheduling strategy. cost (:class:`qibo.models.dbi.double_bracket.DoubleBracketCostFunction`): type of cost function. ref_state (ndarray): reference state for computing the energy fluctuation. Example: .. testcode:: from qibo.models.dbi.double_bracket import DoubleBracketIteration, DoubleBracketGeneratorType from qibo.quantum_info import random_hermitian from qibo.hamiltonians import Hamiltonian nqubits = 4 h0 = random_hermitian(2**nqubits, seed=2) dbf = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0)) # diagonalized matrix dbf.h References: 1. M. Gluza, *Double-bracket quantum algorithms for diagonalization*. `arXiv:2206.11772 [quant-ph] <https://arxiv.org/abs/2206.11772>`_. """ def __init__( self, hamiltonian, mode=DoubleBracketGeneratorType.canonical, scheduling=DoubleBracketScheduling.grid_search, cost=DoubleBracketCostFunction.off_diagonal_norm, ref_state=None, ): self.h = hamiltonian self.h0 = copy(self.h) self.mode = mode self.scheduling = scheduling self.cost = cost self.ref_state = ref_state def __call__(self, step: float, mode=None, d=None): """We use the following convention: .. math:: H^{'} = U^{\\dagger} \\, H \\, U \\, , where :math:`U=e^{-s\\,W}`, and :math:`W =[D, H]` (or depending on ``mode`` an approximation, see `eval_dbr_unitary`). If :math:`s > 0`, then, for :math:`D = \\Delta(H)`, the GWW DBR will give a :math:`\\sigma`-decrease. References: 1. M. Gluza, *Double-bracket quantum algorithms for diagonalization*. `arXiv:2206.11772 [quant-ph] <https://arxiv.org/abs/2206.11772>`_.""" operator = self.eval_dbr_unitary(step, mode, d) operator_dagger = self.backend.cast( np.array(np.matrix(self.backend.to_numpy(operator)).getH()) ) self.h.matrix = operator_dagger @ self.h.matrix @ operator return operator
[docs] def eval_dbr_unitary( self, step: float, mode=None, d=None, ): """In :meth:`qibo.models.dbi.double_bracket.DoubleBracketIteration.__call__`, we are working in the following convention: .. math:: H^{'} = U^{\\dagger} \\, H \\, U \\, , where :math:`U = e^{-s\\,W}`, and :math:`W = [D, H]` (or an approximation of that by a group commutator). That is convenient because if we switch from the DBI in the Heisenberg picture for the Hamiltonian, we get that the transformation of the state is :math:`|\\psi'\\rangle = U \\, |\\psi\\rangle`, so that .. math:: \\langle H\\rangle_{\\psi'} = \\langle H' \\rangle_\\psi \\, , i.e. when writing the unitary acting on the state dagger notation is avoided). The group commutator must approximate :math:`U = e^{-s\\, [D,H]}`. This is achieved by setting :math:`r = \\sqrt{s}` so that .. math:: V = e^{-i\\,r\\,H} \\, e^{i\\,r\\,D} \\, e^{i\\,r\\,H} \\, e^{-i\\,r\\,D} because .. math:: e^{-i\\,r\\,H} \\, D \\, e^{i\\,r\\,H} = D + i\\,r\\,[D, H] +O(r^2) so .. math:: V \\approx \\exp\\left(i\\,r\\,D + i^2 \\, r^2 \\, [D, H] + O(r^2) -i\\,r\\,D\\right) \\approx U \\, . See the Appendix in Ref. [1] for the complete derivation. References: 1. M. Gluza, *Double-bracket quantum algorithms for diagonalization*. `arXiv:2206.11772 [quant-ph] <https://arxiv.org/abs/2206.11772>`_. """ if mode is None: mode = self.mode if mode is DoubleBracketGeneratorType.canonical: operator = matrix_exponentiation( -1.0j * step, self.commutator(self.diagonal_h_matrix, self.h.matrix), backend=self.backend, ) elif mode is DoubleBracketGeneratorType.single_commutator: if d is None: d = self.diagonal_h_matrix operator = matrix_exponentiation( -1.0j * step, self.commutator(self.backend.cast(d), self.h.matrix), backend=self.backend, ) elif mode is DoubleBracketGeneratorType.group_commutator: if d is None: d = self.diagonal_h_matrix operator = ( self.h.exp(step)
[docs] @ matrix_exponentiation(-step, d, backend=self.backend) @ self.h.exp(-step) @ matrix_exponentiation(step, d, backend=self.backend) ) elif mode is DoubleBracketGeneratorType.group_commutator_third_order: if d is None: d = self.diagonal_h_matrix operator = ( self.h.exp(-step * (np.sqrt(5) - 1) / 2) @ matrix_exponentiation( -step * (np.sqrt(5) - 1) / 2, d, backend=self.backend ) @ self.h.exp(step) @ matrix_exponentiation( step * (np.sqrt(5) + 1) / 2, d, backend=self.backend ) @ self.h.exp(-step * (3 - np.sqrt(5)) / 2) @ matrix_exponentiation(-step, d, backend=self.backend) ) operator = ( matrix_exponentiation(step, d, backend=self.backend) @ self.h.exp(step * (3 - np.sqrt(5)) / 2) @ matrix_exponentiation( -step * (np.sqrt(5) + 1) / 2, d, backend=self.backend ) @ self.h.exp(-step) @ matrix_exponentiation( step * (np.sqrt(5) - 1) / 2, d, backend=self.backend ) @ self.h.exp(step * (np.sqrt(5) - 1) / 2) ) else: raise_error( NotImplementedError, f"The mode {mode} is not supported" ) # pragma: no cover return operator
@staticmethod def commutator(a, b): """Compute commutator between two arrays.""" return commutator(a, b)
@property def diagonal_h_matrix(self): """Diagonal H matrix.""" return self.backend.cast(np.diag(np.diag(self.backend.to_numpy(self.h.matrix)))) @property def off_diag_h(self): """Off-diagonal H matrix.""" return self.h.matrix - self.diagonal_h_matrix @property def off_diagonal_norm(self): """Hilbert Schmidt norm of off-diagonal part of H matrix, namely :math:`\\text{Tr}(\\sqrt{A^{\\dagger} A})`.""" off_diag_h_dag = self.backend.cast( np.matrix(self.backend.to_numpy(self.off_diag_h)).getH() ) return np.sqrt( np.real(np.trace(self.backend.to_numpy(off_diag_h_dag @ self.off_diag_h))) ) @property def backend(self): """Get Hamiltonian's backend.""" return self.h0.backend @property def nqubits(self): """Number of qubits.""" return self.h.nqubits
[docs] def least_squares(self, d: np.array): """Least squares cost function.""" d = self.backend.to_numpy(d) return np.real( 0.5 * np.linalg.norm(d) ** 2 - np.trace(self.backend.to_numpy(self.h.matrix) @ d) )
[docs] def choose_step( self, d: Optional[np.array] = None, scheduling: Optional[DoubleBracketScheduling] = None, **kwargs, ): """Calculate the optimal step using respective the `scheduling` methods.""" if scheduling is None: scheduling = self.scheduling step = scheduling(self, d=d, **kwargs) # TODO: write test for this case if ( step is None and scheduling is DoubleBracketScheduling.polynomial_approximation ): # pragma: no cover kwargs["n"] = kwargs.get("n", 3) kwargs["n"] += 1 # if n==n_max, return None step = scheduling(self, d=d, **kwargs) # if for a given polynomial order n, no solution is found, we increase the order of the polynomial by 1 return step
[docs] def loss(self, step: float, d: np.array = None, look_ahead: int = 1): """ Compute loss function distance between `look_ahead` steps. Args: step (float): iteration step. d (np.array): diagonal operator, use canonical by default. look_ahead (int): number of iteration steps to compute the loss function; """ # copy initial hamiltonian h_copy = copy(self.h) for _ in range(look_ahead): self.__call__(mode=self.mode, step=step, d=d) # loss values depending on the cost function if self.cost is DoubleBracketCostFunction.off_diagonal_norm: loss = self.off_diagonal_norm elif self.cost is DoubleBracketCostFunction.least_squares: loss = self.least_squares(d) elif self.cost == DoubleBracketCostFunction.energy_fluctuation: loss = self.energy_fluctuation(self.ref_state) # set back the initial configuration self.h = h_copy return loss
[docs] def energy_fluctuation(self, state): """ Evaluate energy fluctuation. .. math:: \\Xi(\\mu) = \\sqrt{\\langle\\mu|\\hat{H}^2|\\mu\\rangle - \\langle\\mu|\\hat{H}|\\mu\\rangle^2} \\, for a given state :math:`|\\mu\\rangle`. Args: state (np.ndarray): quantum state to be used to compute the energy fluctuation with H. """ return self.h.energy_fluctuation(state)
[docs] def sigma(self, h: np.array): """Returns the off-diagonal restriction of matrix `h`.""" return self.backend.cast(h) - self.backend.cast( np.diag(np.diag(self.backend.to_numpy(h))) )
[docs] def generate_gamma_list(self, n: int, d: np.array): r"""Computes the n-nested Gamma functions, where $\Gamma_k=[W,...,[W,[W,H]]...]$, where we take k nested commutators with $W = [D, H]$""" W = self.commutator(self.backend.cast(d), self.sigma(self.h.matrix)) gamma_list = [self.h.matrix] for _ in range(n - 1): gamma_list.append(self.commutator(W, gamma_list[-1])) return gamma_list
def cost_expansion(self, d, n): d = self.backend.cast(d) if self.cost is DoubleBracketCostFunction.off_diagonal_norm: coef = off_diagonal_norm_polynomial_expansion_coef(self, d, n) elif self.cost is DoubleBracketCostFunction.least_squares: coef = least_squares_polynomial_expansion_coef( self, d, n, backend=self.backend ) elif self.cost is DoubleBracketCostFunction.energy_fluctuation: coef = energy_fluctuation_polynomial_expansion_coef( self, d, n, self.ref_state ) else: # pragma: no cover raise ValueError(f"Cost function {self.cost} not recognized.") return coef