diff --git a/main.py b/main.py index 2d94851..45ee58d 100644 --- a/main.py +++ b/main.py @@ -2,7 +2,7 @@ import json import logging import pathlib -from typing import Optional +from functools import partial import numpy as np @@ -18,12 +18,13 @@ # boostvqe's from boostvqe.ansatze import build_circuit from boostvqe.plotscripts import plot_gradients, plot_loss -from boostvqe.shotnoise import loss_shots +from boostvqe.training_utils import vqe_loss from boostvqe.utils import ( DBI_D_MATRIX, DBI_ENERGIES, DBI_FLUCTUATIONS, DBI_STEPS, + DELTA, FLUCTUATION_FILE, GRADS_FILE, HAMILTONIAN_FILE, @@ -36,7 +37,6 @@ results_dump, rotate_h_with_vqe, train_vqe, - vqe_loss, ) DEFAULT_DELTA = 0.5 @@ -49,37 +49,32 @@ def main(args): """VQE training.""" # set backend and number of classical threads - accuracy = args.accuracy - - if accuracy == 0.0: - accuracy = None - if args.platform is not None: qibo.set_backend(backend=args.backend, platform=args.platform) else: qibo.set_backend(backend=args.backend) args.platform = GlobalBackend().platform - qibo.set_threads(args.nthreads) + if args.optimizer_options is None: + opt_options = {} + else: + opt_options = json.loads(args.optimizer_options) # setup the results folder logging.info("Set VQE") path = pathlib.Path(create_folder(generate_path(args))) ham = getattr(hamiltonians, args.hamiltonian)(nqubits=args.nqubits) - target_energy = float(min(ham.eigenvalues())) - circ0 = build_circuit(nqubits=args.nqubits, nlayers=args.nlayers) + target_energy = np.real(np.min(np.asarray(ham.eigenvalues()))) + circ0 = build_circuit( + nqubits=args.nqubits, + nlayers=args.nlayers, + ) circ = circ0.copy(deep=True) backend = ham.backend zero_state = backend.zero_state(args.nqubits) - # build hamiltonian and variational quantum circuit - if args.shot_train: - loss = lambda params, circ, _ham: loss_shots( - params, circ, _ham, delta=DEFAULT_DELTA, nshots=args.nshots - ) - else: - loss = vqe_loss + loss = partial(vqe_loss, delta=DELTA, nshots=args.nshots) # fix numpy seed to ensure replicability of the experiment np.random.seed(int(args.seed)) @@ -90,7 +85,7 @@ def main(args): # dbi lists boost_energies, boost_fluctuations_dbi, boost_steps, boost_d_matrix = {}, {}, {}, {} # hamiltonian history - hamiltonians_history = [] + fun_eval, hamiltonians_history = [], [] hamiltonians_history.append(ham.matrix) new_hamiltonian = ham args.nboost += 1 @@ -115,13 +110,17 @@ def main(args): niterations=args.boost_frequency, nmessage=1, loss=loss, - accuracy=accuracy, + training_options=opt_options, ) # append results to global lists params_history[b] = np.array(partial_params_history) loss_history[b] = np.array(partial_loss_history) grads_history[b] = np.array(partial_grads_history) fluctuations[b] = np.array(partial_fluctuations) + # this works with scipy.optimize.minimize only + if args.optimizer not in ["sgd", "cma"]: + fun_eval.append(int(partial_results[2].nfev)) + # build new hamiltonian using trained VQE if b != args.nboost - 1: new_hamiltonian_matrix = rotate_h_with_vqe(hamiltonian=ham, vqe=vqe) @@ -176,16 +175,34 @@ def main(args): initial_parameters = np.zeros(len(initial_parameters)) circ.set_parameters(initial_parameters) + # reduce the learning rate after DBI has been applied + if "learning_rate" in opt_options: + opt_options["learning_rate"] *= args.decay_rate_lr + + best_loss = min(np.min(array) for array in loss_history.values()) + + opt_results = partial_results[2] # save final results output_dict = vars(args) output_dict.update( { "true_ground_energy": target_energy, - "accuracy": args.accuracy, + "feval": list(fun_eval), "energy": float(vqe.hamiltonian.expectation(zero_state)), "fluctuations": float(vqe.hamiltonian.energy_fluctuation(zero_state)), + "reached_accuracy": float(np.abs(target_energy - best_loss)), } ) + # this works only with scipy.optimize.minimize + if args.optimizer not in ["sgd", "cma"]: + output_dict.update( + { + "best_loss": float(opt_results.fun), + "success": bool(opt_results.success), + "message": opt_results.message, + "feval": list(fun_eval), + } + ) np.savez( path / LOSS_FILE, **{json.dumps(key): np.array(value) for key, value in loss_history.items()}, @@ -243,11 +260,19 @@ def main(args): parser.add_argument( "--optimizer", default="Powell", type=str, help="Optimizer used by VQE" ) + parser.add_argument( + "--optimizer_options", + type=str, + help="Options to customize the optimizer training", + ) parser.add_argument( "--tol", default=TOL, type=float, help="Absolute precision to stop VQE training" ) parser.add_argument( - "--accuracy", default=0.0, type=float, help="Threshold accuracy" + "--decay_rate_lr", + default=1.0, + type=float, + help="Decay factor of the learning rate if sgd is used", ) parser.add_argument( "--nqubits", default=6, type=int, help="Number of qubits for Hamiltonian / VQE" @@ -303,15 +328,9 @@ def main(args): default=SEED, help="Random seed", ) - parser.add_argument( - "--shot_train", - action=argparse.BooleanOptionalAction, - help="If True the Hamiltonian expactation value is evaluate with the shots, otherwise with the state vector", - ) parser.add_argument( "--nshots", type=int, - default=10000, help="number of shots", ) args = parser.parse_args() diff --git a/run.sh b/run.sh index c30cc8c..294120a 100755 --- a/run.sh +++ b/run.sh @@ -1,20 +1,25 @@ #!/bin/bash -#SBATCH --job-name=dbi8q1l -#SBATCH --output=dbi_8q1l.log +#SBATCH --job-name=adamch +#SBATCH --output=bp_regime.log -NQUBITS=8 +NQUBITS=5 NLAYERS=1 -DBI_STEPS=2 -NBOOST=1 -NSHOTS=10000 -OPTIMIZER="BFGS" -TOL=0.00001 -BOOST_FREQUENCY=10 -ACC=0.1 -python main.py --nqubits $NQUBITS --nlayers $NLAYERS --optimizer $OPTIMIZER \ - --output_folder results/pure_vqe --backend numpy --tol $TOL \ - --dbi_step $DBI_STEPS --seed 42 \ - --boost_frequency $BOOST_FREQUENCY --accuracy $ACC --nboost $NBOOST +DBI_STEPS=0 +NBOOST=0 +BOOST_FREQUENCY=100 + +NSHOTS=1000 +SEED=42 - +OPTIMIZER="sgd" +BACKEND="tensorflow" +OPTIMIZER_OPTIONS="{ \"optimizer\": \"Adam\", \"learning_rate\": 0.1, \"nmessage\": 1, \"nepochs\": $BOOST_FREQUENCY }" +DECAY_RATE_LR=0.05 + +python main.py --nqubits $NQUBITS --nlayers $NLAYERS --optimizer $OPTIMIZER \ + --output_folder results/debugging_decay --backend $BACKEND \ + --dbi_step $DBI_STEPS --seed $SEED \ + --boost_frequency $BOOST_FREQUENCY --nboost $NBOOST \ + --optimizer_options "$OPTIMIZER_OPTIONS" \ + --decay_rate_lr $DECAY_RATE_LR diff --git a/run_sgd.sh b/run_sgd.sh new file mode 100755 index 0000000..5374984 --- /dev/null +++ b/run_sgd.sh @@ -0,0 +1,24 @@ +#!/bin/bash +#SBATCH --job-name=boostvqe +#SBATCH --output=boostvqe.log + +NQUBITS=4 +NLAYERS=2 + +DBI_STEPS=2 +NBOOST=2 +BOOST_FREQUENCY=10 + +NSHOTS=10000 +TOL=1e-8 +ACC=0.5 + +OPTIMIZER="sgd" +BACKEND="tensorflow" +OPTIMIZER_OPTIONS="{ \"optimizer\": \"Adagrad\", \"learning_rate\": 0.1, \"nmessage\": 1, \"nepochs\": $BOOST_FREQUENCY }" + +python main.py --nqubits $NQUBITS --nlayers $NLAYERS --optimizer $OPTIMIZER \ + --output_folder results/debugging --backend $BACKEND --tol $TOL \ + --dbi_step $DBI_STEPS --seed 42 \ + --boost_frequency $BOOST_FREQUENCY --nboost $NBOOST \ + --optimizer_options "$OPTIMIZER_OPTIONS" diff --git a/run_sgd_exact.sh b/run_sgd_exact.sh new file mode 100755 index 0000000..c50c447 --- /dev/null +++ b/run_sgd_exact.sh @@ -0,0 +1,23 @@ +#!/bin/bash +#SBATCH --job-name=boostvqe +#SBATCH --output=boostvqe.log + +NQUBITS=7 +NLAYERS=3 + +DBI_STEPS=0 +NBOOST=0 +BOOST_FREQUENCY=1000 + +NSHOTS=1000 +TOL=1e-8 + +OPTIMIZER="sgd" +BACKEND="tensorflow" +OPTIMIZER_OPTIONS="{ \"optimizer\": \"Adagrad\", \"learning_rate\": 0.05, \"nmessage\": 1, \"nepochs\": $BOOST_FREQUENCY }" + +python main.py --nqubits $NQUBITS --nlayers $NLAYERS --optimizer $OPTIMIZER \ + --output_folder results/sgd_exact --backend $BACKEND --tol $TOL \ + --dbi_step $DBI_STEPS --seed 42 \ + --boost_frequency $BOOST_FREQUENCY --nboost $NBOOST \ + --optimizer_options "$OPTIMIZER_OPTIONS" diff --git a/run_sgd_hybrid.sh b/run_sgd_hybrid.sh new file mode 100755 index 0000000..f0881e2 --- /dev/null +++ b/run_sgd_hybrid.sh @@ -0,0 +1,23 @@ +#!/bin/bash +#SBATCH --job-name=boostvqe +#SBATCH --output=boostvqe.log + +NQUBITS=7 +NLAYERS=3 + +DBI_STEPS=2 +NBOOST=1 +BOOST_FREQUENCY=100 + +NSHOTS=1000 +TOL=1e-8 + +OPTIMIZER="sgd" +BACKEND="tensorflow" +OPTIMIZER_OPTIONS="{ \"optimizer\": \"Adagrad\", \"learning_rate\": 0.05, \"nmessage\": 1, \"nepochs\": $BOOST_FREQUENCY }" + +python main.py --nqubits $NQUBITS --nlayers $NLAYERS --optimizer $OPTIMIZER \ + --output_folder results/sgd_exact_hybrid --backend $BACKEND --tol $TOL \ + --dbi_step $DBI_STEPS --seed 42 \ + --boost_frequency $BOOST_FREQUENCY --nboost $NBOOST \ + --optimizer_options "$OPTIMIZER_OPTIONS" diff --git a/src/boostvqe/ansatze.py b/src/boostvqe/ansatze.py index 237e50f..95b5d99 100644 --- a/src/boostvqe/ansatze.py +++ b/src/boostvqe/ansatze.py @@ -1,7 +1,10 @@ -from qibo import gates +import qibo +from qibo import gates, get_backend from qibo.backends import construct_backend from qibo.models import Circuit +from boostvqe.training_utils import vqe_loss + def build_circuit(nqubits, nlayers): """Build qibo's aavqe example circuit.""" @@ -38,3 +41,96 @@ def compute_gradients(parameters, circuit, hamiltonian): ) return hamiltonian.backend.cast(tape.gradient(expectation, parameters)) + + +class VQE: + """This class implements the variational quantum eigensolver algorithm.""" + + from qibo import optimizers + + def __init__(self, circuit, hamiltonian): + """Initialize circuit ansatz and hamiltonian.""" + self.circuit = circuit + self.hamiltonian = hamiltonian + self.backend = hamiltonian.backend + + def minimize( + self, + initial_state, + method="Powell", + loss_func=None, + jac=None, + hess=None, + hessp=None, + bounds=None, + constraints=(), + tol=None, + callback=None, + options=None, + compile=False, + processes=None, + ): + """Search for parameters which minimizes the hamiltonian expectation. + + Args: + initial_state (array): a initial guess for the parameters of the + variational circuit. + method (str): the desired minimization method. + See :meth:`qibo.optimizers.optimize` for available optimization + methods. + loss (callable): loss function, the default one is :func:`qibo.models.utils.vqe_loss`. + jac (dict): Method for computing the gradient vector for scipy optimizers. + hess (dict): Method for computing the hessian matrix for scipy optimizers. + hessp (callable): Hessian of objective function times an arbitrary + vector for scipy optimizers. + bounds (sequence or Bounds): Bounds on variables for scipy optimizers. + constraints (dict): Constraints definition for scipy optimizers. + tol (float): Tolerance of termination for scipy optimizers. + callback (callable): Called after each iteration for scipy optimizers. + options (dict): a dictionary with options for the different optimizers. + compile (bool): whether the TensorFlow graph should be compiled. + processes (int): number of processes when using the paralle BFGS method. + + Return: + The final expectation value. + The corresponding best parameters. + The optimization result object. For scipy methods it returns + the ``OptimizeResult``, for ``'cma'`` the ``CMAEvolutionStrategy.result``, + and for ``'sgd'`` the options used during the optimization. + """ + if loss_func is None: + loss_func = vqe_loss + if compile: + loss = self.hamiltonian.backend.compile(loss_func) + else: + loss = loss_func + + if method == "cma": + dtype = self.hamiltonian.backend.np.float64 + loss = lambda p, c, h: dtype(loss_func(p, c, h)) + elif method != "sgd": + loss = lambda p, c, h: self.hamiltonian.backend.to_numpy(loss_func(p, c, h)) + + result, parameters, extra = self.optimizers.optimize( + loss, + initial_state, + args=(self.circuit, self.hamiltonian), + method=method, + jac=jac, + hess=hess, + hessp=hessp, + bounds=bounds, + constraints=constraints, + tol=tol, + callback=callback, + options=options, + compile=compile, + processes=processes, + backend=self.hamiltonian.backend, + ) + self.circuit.set_parameters(parameters) + return result, parameters, extra + + def energy_fluctuation(self, state): + """Compute Energy Fluctuation.""" + return self.hamiltonian.energy_fluctuation(state) diff --git a/src/boostvqe/plotscripts.py b/src/boostvqe/plotscripts.py index 0f6d1f0..d035c64 100644 --- a/src/boostvqe/plotscripts.py +++ b/src/boostvqe/plotscripts.py @@ -190,29 +190,114 @@ def plot_gradients( plt.savefig(f"{path}/grads_{title}.pdf", bbox_inches="tight") -def plot_nruns_result( +def plot_loss_nruns( path, training_specs, title="", save=True, width=0.5, ): - - losses = [] + """ + Plot loss with confidence belt. + """ + + losses_dbi = [] + losses_vqe = [] - for f in os.listdir(path): + for i, f in enumerate(os.listdir(path)): + this_path = path + "/" + f + "/" + if i == 0: + with open(this_path + OPTIMIZATION_FILE, 'r') as file: + config = json.load(file) + target_energy = config["true_ground_energy"] + # accumulating dictionaries with results for each boost if training_specs in f: - losses.append(np.load(path + "/" + f)) - - losses = np.array(losses) - means = np.mean(losses, axis=0) - stds = np.std(losses, axis=0) + losses_vqe.append(dict(np.load(this_path + f"{LOSS_FILE + '.npz'}"))) + losses_dbi.append(dict(np.load(this_path + f"{DBI_ENERGIES + '.npz'}"))) + loss_vqe, dbi_energies, stds_vqe, stds_dbi = {}, {}, {}, {} + plt.figure(figsize=(10 * width, 10 * width * 6 / 8)) - plt.plot(means, color=BLUE) - plt.fill_between(means - stds, means + stds, alpha=0.3, color=BLUE) plt.title(title) - plt.xlabel("Iteration") - plt.ylabel(r"$\langle L \rangle$") + + for i in range(config["nboost"]): + this_vqe_losses, this_dbi_losses = [], [] + for d in range(len(loss_vqe)): + this_vqe_losses.append(losses_vqe[d][str(i)]) + this_dbi_losses.append(losses_dbi[d][str(i)]) + + loss_vqe.update({str(i): np.mean(np.asarray(this_vqe_losses), axis=0)}) + dbi_energies.update({str(i): np.mean(np.asarray(this_dbi_losses), axis=0)}) + stds_vqe.update({str(i): np.std(np.asarray(this_vqe_losses), axis=0)}) + stds_dbi.update({str(i): np.std(np.asarray(this_dbi_losses), axis=0)}) + + for i in range(config["nboost"]): + start = ( + 0 + if str(i - 1) not in loss_vqe + else sum( + len(loss_vqe[str(j)]) + len(dbi_energies[str(j)]) + for j in range(config["nboost"]) + if j < i + ) + - 2 * i + ) + plt.plot( + np.arange(start, len(loss_vqe[str(i)]) + start), + loss_vqe[str(i)], + color=BLUE, + lw=1.5, + label="VQE", + ) + plt.plot( + np.arange( + len(loss_vqe[str(i)]) + start - 1, + len(dbi_energies[str(i)]) + len(loss_vqe[str(i)]) + start - 1, + ), + dbi_energies[str(i)], + color=RED, + lw=1.5, + label="DBI", + ) + plt.fill_between( + np.arange(start, len(loss_vqe[str(i)]) + start), + loss_vqe[str(i)] - stds_vqe[str(i)], + loss_vqe[str(i)] + stds_vqe[str(i)], + color=BLUE, + alpha=0.4, + ) + plt.fill_between( + np.arange( + len(loss_vqe[str(i)]) + start - 1, + len(dbi_energies[str(i)]) + len(loss_vqe[str(i)]) + start - 1, + ), + dbi_energies[str(i)] - stds_vqe[str(i)], + dbi_energies[str(i)] + stds_dbi[str(i)], + color=RED, + alpha=0.4, + ) + + max_length = ( + sum(len(l) for l in list(dbi_energies.values())) + + sum(len(l) for l in list(loss_vqe.values())) + - 2 * config["nboost"] + + 1 + ) + plt.hlines( + target_energy, + 0, + max_length, + color="black", + lw=1, + label="Target energy", + ls="-", + ) + plt.xlabel("Iterations") + plt.ylabel("Loss") + handles, labels = plt.gca().get_legend_handles_labels() + by_label = dict(zip(labels, handles)) + plt.legend(by_label.values(), by_label.keys()) if save: - plt.savefig("runs.png") \ No newline at end of file + plt.savefig(f"{path}/loss_{title}.pdf", bbox_inches="tight") + + diff --git a/src/boostvqe/shotnoise.py b/src/boostvqe/shotnoise.py deleted file mode 100644 index bb2d493..0000000 --- a/src/boostvqe/shotnoise.py +++ /dev/null @@ -1,34 +0,0 @@ -from qibo import gates, hamiltonians -from qibo.symbols import Z - -NITER = 20 - - -def loss_shots( - params, - circ, - _ham, # We need it to uniform all the loss functions - delta, - nshots, -): - """Train the VQE with the shots, this function is specific to the XXZ Hamiltonian""" - - # Diagonal Hamiltonian to evaluate the expactation value - hamiltonian = sum(Z(i) * Z(i + 1) for i in range(circ.nqubits - 1)) - hamiltonian += Z(0) * Z(circ.nqubits - 1) - hamiltonian = hamiltonians.SymbolicHamiltonian(hamiltonian) - circ.set_parameters(params) - coefficients = [1, 1, delta] - mgates = ["X", "Y", "Z"] - expectation_value = 0 - for i, mgate in enumerate(mgates): - circ1 = circ.copy(deep=True) - if mgate != "Z": - circ1.add(gates.M(*range(circ1.nqubits), basis=getattr(gates, mgate))) - else: - circ1.add(gates.M(*range(circ1.nqubits))) - result = circ1(nshots=nshots) - expectation_value += coefficients[i] * hamiltonian.expectation_from_samples( - result.frequencies() - ) - return expectation_value diff --git a/src/boostvqe/training_utils.py b/src/boostvqe/training_utils.py new file mode 100644 index 0000000..3f67782 --- /dev/null +++ b/src/boostvqe/training_utils.py @@ -0,0 +1,146 @@ +import os + +os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3" + +import numpy as np +import tensorflow as tf + +tf.get_logger().setLevel("ERROR") + +from qibo import gates, hamiltonians +from qibo.backends import NumpyBackend, TensorflowBackend +from qibo.config import raise_error +from qibo.hamiltonians import AbstractHamiltonian +from qibo.symbols import Z + + +def vqe_loss(params, circuit, hamiltonian, nshots=None, delta=0.5): + """ + Evaluate the hamiltonian expectation values of the + circuit final state. + + TODO: fix the following statement. + IMPORTANT: this works only for Heisemberg hamiltonians XXZ. + """ + circ = circuit.copy(deep=True) + circ.set_parameters(params) + + if isinstance(hamiltonian.backend, TensorflowBackend) and nshots is not None: + expectation_value = _exp_with_tf(circ, hamiltonian, nshots, delta) + elif nshots is None: + expectation_value = _exact(circ, hamiltonian) + else: + expectation_value = _with_shots(circ, hamiltonian, nshots, delta) + return expectation_value + + +def _exact(circ, hamiltonian): + """Helper function to compute expectation value of Heisemberg hamiltonian.""" + expectation_value = hamiltonian.expectation( + hamiltonian.backend.execute_circuit(circuit=circ).state() + ) + return expectation_value + + +def _with_shots(circ, hamiltonian, nshots, delta=0.5, exec_backend=None): + """Helper function to compute XXZ expectation value from frequencies.""" + + # we may prefer run this on a different backend (e.g. with TF and PSR) + if exec_backend is None: + exec_backend = hamiltonian.backend + + hamiltonian = sum(Z(i) * Z(i + 1) for i in range(circ.nqubits - 1)) + hamiltonian += Z(0) * Z(circ.nqubits - 1) + hamiltonian = hamiltonians.SymbolicHamiltonian(hamiltonian) + coefficients = [1, 1, delta] + mgates = ["X", "Y", "Z"] + expectation_value = 0 + for i, mgate in enumerate(mgates): + circ1 = circ.copy(deep=True) + if mgate != "Z": + circ1.add(gates.M(*range(circ1.nqubits), basis=getattr(gates, mgate))) + else: + circ1.add(gates.M(*range(circ1.nqubits))) + + expval_contribution = exec_backend.execute_circuit( + circuit=circ1, nshots=nshots + ).expectation_from_samples(hamiltonian) + expectation_value += coefficients[i] * expval_contribution + return expectation_value + + +def _exp_with_tf(circuit, hamiltonian, nshots=None, delta=0.5): + params = circuit.get_parameters() + nparams = len(circuit.get_parameters()) + + @tf.custom_gradient + def _expectation(params): + def grad(upstream): + gradients = [] + for p in range(nparams): + gradients.append( + upstream + * parameter_shift( + circuit=circuit, + hamiltonian=hamiltonian, + parameter_index=p, + nshots=nshots, + delta=delta, + exec_backend=NumpyBackend(), + ) + ) + return gradients + + if nshots is None: + expectation_value = _exact(circuit, hamiltonian) + else: + expectation_value = _with_shots(circuit, hamiltonian, nshots, delta) + return expectation_value, grad + + return _expectation(params) + + +def parameter_shift( + hamiltonian, + circuit, + parameter_index, + exec_backend, + nshots=None, + delta=0.5, +): + """Parameter Shift Rule.""" + + if parameter_index > len(circuit.get_parameters()): + raise_error(ValueError, """This index is out of bounds.""") + + if not isinstance(hamiltonian, AbstractHamiltonian): + raise_error( + TypeError, + "hamiltonian must be a qibo.hamiltonians.Hamiltonian or qibo.hamiltonians.SymbolicHamiltonian object", + ) + + gate = circuit.associate_gates_with_parameters()[parameter_index] + generator_eigenval = gate.generator_eigenvalue() + + s = np.pi / (4 * generator_eigenval) + + original = np.asarray(circuit.get_parameters()).copy() + shifted = original.copy() + + shifted[parameter_index] += s + circuit.set_parameters(shifted) + + if nshots is None: + forward = _exact(circ=circuit, hamiltonian=hamiltonian) + shifted[parameter_index] -= 2 * s + circuit.set_parameters(shifted) + backward = _exact(circ=circuit, hamiltonian=hamiltonian) + + else: + forward = _with_shots(circuit, hamiltonian, nshots, delta, exec_backend) + shifted[parameter_index] -= 2 * s + circuit.set_parameters(shifted) + backward = _with_shots(circuit, hamiltonian, nshots, delta, exec_backend) + + circuit.set_parameters(original) + return float(generator_eigenval * (forward - backward)) diff --git a/src/boostvqe/utils.py b/src/boostvqe/utils.py index 737d2b4..1100f44 100644 --- a/src/boostvqe/utils.py +++ b/src/boostvqe/utils.py @@ -1,12 +1,11 @@ import json import logging -import warnings from pathlib import Path import numpy as np -from qibo.models.variational import VQE +from qibo import get_backend -from boostvqe.ansatze import compute_gradients +from boostvqe.ansatze import VQE, compute_gradients OPTIMIZATION_FILE = "optimization_results.json" PARAMS_FILE = "parameters_history.npy" @@ -17,6 +16,7 @@ GRADS_FILE = "gradients" HAMILTONIAN_FILE = "hamiltonian_matrix.npz" SEED = 42 +DELTA = 0.5 TOL = 1e-10 DBI_ENERGIES = "dbi_energies" DBI_FLUCTUATIONS = "dbi_fluctuations" @@ -65,22 +65,9 @@ def callback_energy_fluctuations(params, circuit, hamiltonian): circ.set_parameters(params) result = hamiltonian.backend.execute_circuit(circ) final_state = result.state() - return hamiltonian.energy_fluctuation(final_state) -def vqe_loss(params, circuit, hamiltonian): - """ - Evaluate the hamiltonian expectation values of the - circuit final state. - """ - circ = circuit.copy(deep=True) - circ.set_parameters(params) - result = hamiltonian.backend.execute_circuit(circ) - final_state = result.state() - return hamiltonian.expectation(final_state) - - def train_vqe( circ, ham, @@ -90,16 +77,21 @@ def train_vqe( loss, niterations=None, nmessage=1, - accuracy=None, + training_options=None, ): """Helper function which trains the VQE according to `circ` and `ham`.""" - params_history, loss_list, fluctuations, results_history, grads_history = ( - [], + params_history, loss_list, fluctuations, grads_history = ( [], [], [], [], ) + + if training_options is None: + options = {} + else: + options = training_options + circ.set_parameters(initial_parameters) vqe = VQE( @@ -107,8 +99,6 @@ def train_vqe( hamiltonian=ham, ) - target_energy = float(np.min(ham.eigenvalues())) - def callbacks( params, vqe=vqe, @@ -138,24 +128,18 @@ def callbacks( if niterations is not None and iteration_count % nmessage == 0: logging.info(f"Optimization iteration {iteration_count}/{niterations}") - - if iteration_count >= niterations: - logging.info("Maximum number of iterations reached.") - raise StopIteration("Target accuracy reached.") - - if accuracy is not None: - if np.abs(np.min(loss_list) - target_energy) <= accuracy: - logging.info("Target accuracy reached.") - raise StopIteration("Target accuracy reached.") + logging.info(f"Loss {energy:.5}") callbacks(initial_parameters) logging.info("Minimize the energy") + results = vqe.minimize( initial_parameters, method=optimizer, callback=callbacks, tol=tol, loss_func=loss, + options=options, ) return ( @@ -204,4 +188,6 @@ def apply_dbi_steps(dbi, nsteps, stepsize=0.01, optimize_step=False): energies.append(dbi.h.expectation(zero_state)) fluctuations.append(dbi.energy_fluctuation(zero_state)) hamiltonians.append(dbi.h.matrix) + + logging.info(f"DBI energies: {energies}") return hamiltonians, energies, fluctuations, steps, d_matrix, operators