Merge pull request #53 from qiboteam/gci_opt

Integrate optimizers into compiling XXZ code

compile.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+#SBATCH --job-name=boost
+#SBATCH --output=boost.log
+
+OPTIMIZATION_METHOD="cma"
+OPTIMIZATION_CONFIG="{ \"maxiter\": 2}"
+
+python3 compiling.py --backend numpy --path "./results/XXZ_5seeds/sgd_10q_6l_27/" \
+                     --epoch 500 --steps 2 --optimization_method $OPTIMIZATION_METHOD \
+                     --optimization_config "$OPTIMIZATION_CONFIG"

compiling.py

@@ -0,0 +1,248 @@
+import argparse
+import json
+import logging
+import pathlib
+import time
+from copy import deepcopy
+
+import numpy as np
+import qibo
+
+qibo.set_backend("numpy")
+from qibo import hamiltonians
+from qibo.backends import construct_backend
+from qibo.quantum_info.metrics import fidelity
+
+from boostvqe.ansatze import VQE, build_circuit
+from boostvqe.models.dbi.double_bracket_evolution_oracles import (
+    FrameShiftedEvolutionOracle,
+    MagneticFieldEvolutionOracle,
+    XXZ_EvolutionOracle,
+)
+from boostvqe.models.dbi.group_commutator_iteration_transpiler import (
+    DoubleBracketRotationType,
+    GroupCommutatorIterationWithEvolutionOracles,
+)
+from boostvqe.utils import (
+    OPTIMIZATION_FILE,
+    PARAMS_FILE,
+    build_circuit,
+    optimize_D,
+    select_recursion_step_gd_circuit,
+)
+
+logging.basicConfig(level=logging.INFO)
+
+
+def dump_config(config: dict, path):
+    config["path"] = config["path"]
+    config["db_rotation"] = config["db_rotation"].value
+    (path / "config.json").write_text(json.dumps(config, indent=4))
+
+
+def main(args):
+    """VQE training."""
+    path = pathlib.Path(args.path)
+    dump_path = (
+        path
+        / f"{args.db_rotation.name}_{args.optimization_method}_{args.epoch}e_{args.steps}s"
+    )
+    dump_path.mkdir(parents=True, exist_ok=True)
+
+    config = json.loads((path / OPTIMIZATION_FILE).read_text())
+    dump_config(deepcopy(vars(args)), path=dump_path)
+
+    if args.optimization_config is None:
+        opt_options = {}
+    else:
+        opt_options = json.loads(args.optimization_config)
+
+    try:
+        params = np.load(path / f"parameters/params_ite{args.epoch}.npy")
+    except FileNotFoundError:
+        params = np.array(
+            np.load(path / PARAMS_FILE, allow_pickle=True).tolist()[0][args.epoch]
+        )
+
+    nqubits = config["nqubits"]
+    nlayers = config["nlayers"]
+    vqe_backend = construct_backend(backend=config["backend"])
+    # TODO: remove delta hardcoded
+    hamiltonian = getattr(hamiltonians, config["hamiltonian"])(
+        nqubits=nqubits, delta=0.5, backend=vqe_backend
+    )
+    vqe = VQE(
+        build_circuit(
+            nqubits=nqubits,
+            nlayers=nlayers,
+        ),
+        hamiltonian=hamiltonian,
+    )
+    vqe.circuit.set_parameters(params)
+
+    base_oracle = XXZ_EvolutionOracle.from_nqubits(
+        nqubits=nqubits, delta=0.5, steps=args.steps, order=args.order
+    )
+    oracle = FrameShiftedEvolutionOracle.from_evolution_oracle(
+        before_circuit=vqe.circuit.invert(),
+        after_circuit=vqe.circuit,
+        base_evolution_oracle=base_oracle,
+    )
+
+    gci = GroupCommutatorIterationWithEvolutionOracles(
+        oracle,
+        args.db_rotation,
+    )
+
+    # TODO: remove hardcoded magnetic field
+    eo_d = MagneticFieldEvolutionOracle.from_b(
+        [4 - np.sin(x / 3) for x in range(nqubits)]
+    )
+    gci.eo_d = eo_d
+    print(
+        f"The gci mode is {gci.double_bracket_rotation_type} rotation with {eo_d.__class__.__name__} as the oracle.\n"
+    )
+    metadata = {}
+
+    this_report = report(vqe, hamiltonian, gci)
+    metadata[0] = this_report
+
+    for gci_step_nmb in range(args.steps):
+        logging.info(
+            "\n################################################################################\n"
+            + f"Optimizing GCI step {gci_step_nmb+1} with optimizer {args.optimization_method}"
+            + "\n################################################################################\n"
+        )
+        it = time.time()
+        if args.optimization_method == "sgd":
+            _, best_s, _, eo_d = select_recursion_step_gd_circuit(
+                gci,
+                mode_dbr_list=[args.db_rotation],
+                step_grid=np.linspace(1e-5, 2e-2, 30),
+                lr_range=(1e-3, 1),
+                nmb_gd_epochs=opt_options["gd_epochs"],
+                threshold=1e-4,
+                max_eval_gd=30,
+                please_be_visual=False,
+                save_path="gci_step",
+            )
+
+        else:
+            if gci_step_nmb == 0:
+                p0 = [0.01]
+                p0.extend([4 - np.sin(x / 3) for x in range(nqubits)])
+            else:
+                p0 = [best_s]
+                p0.extend(best_b)
+            optimized_params = optimize_D(
+                params=p0,
+                gci=gci,
+                method=args.optimization_method,
+                **opt_options,
+            )
+            best_s = optimized_params[0]
+            best_b = optimized_params[1:]
+            eo_d = MagneticFieldEvolutionOracle.from_b(best_b)
+
+        step_data = dict(
+            best_s=best_s,
+            eo_d_name=eo_d.__class__.__name__,
+            eo_d_params=eo_d.params,
+        )
+        logging.info(f"Total optimization time required: {time.time() - it} seconds")
+        gci.mode_double_bracket_rotation = args.db_rotation
+        gci.eo_d = eo_d
+        gci(best_s)
+
+        this_report = report(vqe, hamiltonian, gci)
+        print_report(this_report)
+        metadata[gci_step_nmb + 1] = this_report | step_data
+
+    (dump_path / "boosting_data.json").write_text(json.dumps(metadata, indent=4))
+
+
+def report(vqe, hamiltonian, gci):
+    energies = hamiltonian.eigenvalues()
+    ground_state_energy = float(energies[0])
+    vqe_energy = float(hamiltonian.expectation(vqe.circuit().state()))
+    gci_loss = float(gci.loss())
+    gap = float(energies[1] - energies[0])
+
+    return (
+        dict(
+            nqubits=hamiltonian.nqubits,
+            gci_loss=float(gci_loss),
+            vqe_energy=float(vqe_energy),
+            target_energy=ground_state_energy,
+            diff_vqe_target=vqe_energy - ground_state_energy,
+            diff_gci_target=gci_loss - ground_state_energy,
+            gap=gap,
+            diff_vqe_target_perc=abs(vqe_energy - ground_state_energy)
+            / abs(ground_state_energy)
+            * 100,
+            diff_gci_target_perc=abs(gci_loss - ground_state_energy)
+            / abs(ground_state_energy)
+            * 100,
+            fidelity_witness_vqe=1 - (vqe_energy - ground_state_energy) / gap,
+            fidelity_witness_gci=1 - (gci_loss - ground_state_energy) / gap,
+            fidelity_vqe=fidelity(vqe.circuit().state(), hamiltonian.ground_state()),
+            fidelity_gci=fidelity(
+                gci.get_composed_circuit()().state(), hamiltonian.ground_state()
+            ),
+        )
+        | gci.get_gate_count_dict()
+    )
+
+
+def print_report(report: dict):
+    print(
+        f"\
+    The target energy is {report['target_energy']}\n\
+    The VQE energy is {report['vqe_energy']} \n\
+    The DBQA energy is {report['gci_loss']}. \n\
+    The difference is for VQE is {report['diff_vqe_target']} \n\
+    and for the DBQA {report['diff_gci_target']} \n\
+    which can be compared to the spectral gap {report['gap']}.\n\
+    The relative difference is \n\
+        - for VQE {report['diff_vqe_target_perc']}% \n\
+        - for DBQA {report['diff_gci_target_perc']}%.\n\
+    The energetic fidelity witness of the ground state is: \n\
+        - for the VQE  {report['fidelity_witness_vqe']} \n\
+        - for DBQA {report['fidelity_witness_gci']}\n\
+    The true fidelity is \n\
+        - for the VQE  {report['fidelity_vqe']}\n\
+        - for DBQA {report['fidelity_gci']}\n\
+                    "
+    )
+    print(
+        f"The boosting circuit used {report['nmb_cnot']} CNOT gates coming from compiled XXZ evolution and {report['nmb_cz']} CZ gates from VQE.\n\
+For {report['nqubits']} qubits this gives n_CNOT/n_qubits = {report['nmb_cnot_relative']} and n_CZ/n_qubits = {report['nmb_cz_relative']}"
+    )
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Boosting VQE with DBI.")
+    parser.add_argument("--backend", default="qibojit", type=str, help="Qibo backend")
+    parser.add_argument("--path", type=str, help="Output folder")
+    parser.add_argument(
+        "--epoch", default=-1, type=int, help="VQE epoch where DBI will be applied."
+    )
+    parser.add_argument("--steps", default=2, type=int, help="DBI steps")
+    parser.add_argument("--order", default=2, type=int, help="Suzuki-Trotter order")
+    parser.add_argument(
+        "--db_rotation",
+        type=lambda arg: DoubleBracketRotationType[arg],
+        choices=DoubleBracketRotationType,
+        default="group_commutator_reduced",
+        help="DB rotation type.",
+    )
+    parser.add_argument(
+        "--optimization_method", default="sgd", type=str, help="Optimization method"
+    )
+    parser.add_argument(
+        "--optimization_config",
+        type=str,
+        help="Options to customize the optimizer training.",
+    )
+    args = parser.parse_args()
+    main(args)

extras/load_vqe.py

@@ -0,0 +1,85 @@
+import json
+import time
+from pathlib import Path
+
+import numpy as np
+import qibo
+from qibo import hamiltonians, set_backend
+from qibo.models.dbi.double_bracket import (
+    DoubleBracketGeneratorType,
+    DoubleBracketIteration,
+)
+
+from boostvqe.ansatze import VQE, build_circuit
+from boostvqe.utils import apply_dbi_steps, rotate_h_with_vqe
+
+qibo.set_backend("numpy")
+
+# set the path string which define the results
+path = "../results/vqe_data/with_params/10q3l/sgd_10q_3l_42/"
+
+# set the target epoch to which apply DBQA and the number of steps
+target_epoch = 2000
+dbi_steps = 1
+
+# upload system configuration and parameters for all the training
+with open(path + "optimization_results.json") as file:
+    config = json.load(file)
+
+losses = dict(np.load(path + "energies.npz"))["0"]
+params = np.load(path + f"parameters/params_ite{target_epoch}.npy")
+
+
+# build circuit, hamiltonian and VQE
+hamiltonian = hamiltonians.XXZ(nqubits=config["nqubits"], delta=0.5)
+circuit = build_circuit(config["nqubits"], config["nlayers"], "numpy")
+vqe = VQE(circuit, hamiltonian)
+zero_state = hamiltonian.backend.zero_state(config["nqubits"])
+target_energy = np.min(hamiltonian.eigenvalues())
+
+
+# set target parameters into the VQE
+vqe.circuit.set_parameters(params)
+vqe_state = vqe.circuit().state()
+
+ene1 = hamiltonian.expectation(vqe_state)
+
+# DBQA stuff
+t0 = time.time()
+print("Rotating with VQE")
+new_hamiltonian_matrix = rotate_h_with_vqe(hamiltonian=hamiltonian, vqe=vqe)
+new_hamiltonian = hamiltonians.Hamiltonian(
+    config["nqubits"], matrix=new_hamiltonian_matrix
+)
+print(time.time() - t0)
+dbi = DoubleBracketIteration(
+    hamiltonian=new_hamiltonian,
+    mode=DoubleBracketGeneratorType.single_commutator,
+)
+
+zero_state_t = np.transpose([zero_state])
+energy_h0 = float(dbi.h.expectation(np.array(zero_state_t)))
+fluctuations_h0 = float(dbi.h.energy_fluctuation(zero_state_t))
+
+print("Applying DBI steps")
+(
+    _,
+    dbi_energies,
+    dbi_fluctuations,
+    _,
+    _,
+    _,
+) = apply_dbi_steps(
+    dbi=dbi,
+    nsteps=dbi_steps,
+)
+print(time.time() - t0)
+print(
+    f"\nReached accuracy before DBI at iter {target_epoch}: {np.abs(target_energy - ene1)}"
+)
+print(
+    f"Reached accuracy after DBI at iter {target_epoch}: {np.abs(target_energy - dbi_energies[-1])}"
+)
+print(
+    f"Reached accuracy in the end of VQE long training: {np.abs(target_energy - losses[-1])}"
+)

extras/polynomial_s.py

@@ -0,0 +1,30 @@
+import math
+import numpy as np
+
+
+
+def energy_expectation_polynomial_expansion_smin(
+    dbi_object, d: np.array, n: int = 3, state=0
+):
+    """Return the first root of the Taylor expansion coefficients of energy expectation of `dbi_object` with respect to double bracket rotation duration `s`."""
+    # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H
+    Gamma_list = dbi_object.generate_gamma_list(n + 1, d)
+    # coefficients
+    coef = np.empty(n)
+    state_cast = dbi_object.backend.cast(state)
+    state_dag = dbi_object.backend.cast(state.conj().T)
+    exp_list = np.array([1 / math.factorial(k) for k in range(n + 1)])
+    for i in range(n):
+        coef[i] = np.real(
+            exp_list[i] * state_dag @ Gamma_list[i+1] @ state_cast
+        )
+    coef = list(reversed(coef))
+
+    roots = np.roots(coef)
+    real_positive_roots = [
+        np.real(root) for root in roots if np.imag(root) < 1e-3 and np.real(root) > 0
+    ]
+    # solution exists, return minimum s
+    if len(real_positive_roots) > 0:
+        losses = [dbi_object.loss(step=root, d=d) for root in real_positive_roots]
+    return real_positive_roots[losses.index(min(losses))]

notebooks/compiling_XXZ_test.py

@@ -0,0 +1,15 @@
+from qibo import hamiltonians
+import numpy as np
+from boostvqe.compiling_XXZ import *
+
+t = 0.01
+steps = 3
+delta=0.5
+nqubits=6
+
+h_xxz = hamiltonians.XXZ(nqubits=nqubits, delta = delta)
+u = h_xxz.exp(t)
+circ = nqubit_XXZ_decomposition(nqubits=nqubits,t=t,delta=delta,steps=steps)
+v = circ.unitary()
+print(np.linalg.norm(u-np.exp(nqubits*steps*1j*np.pi/4)*v))
+print(circ.draw())