Implementing and running new quantum hortogonal ansatz (with RBS gates)

compiling.py

@@ -12,7 +12,7 @@
 from qibo.backends import construct_backend
 from qibo.quantum_info.metrics import fidelity
 
-from boostvqe.ansatze import VQE, build_circuit
+from boostvqe import ansatze
 from boostvqe.models.dbi import double_bracket_evolution_oracles
 from boostvqe.models.dbi.double_bracket_evolution_oracles import (
     FrameShiftedEvolutionOracle,
@@ -27,7 +27,6 @@
 from boostvqe.utils import (
     OPTIMIZATION_FILE,
     PARAMS_FILE,
-    build_circuit,
     optimize_D,
     select_recursion_step_gd_circuit,
 )
@@ -74,11 +73,12 @@ def main(args):
     hamiltonian = getattr(hamiltonians, config["hamiltonian"])(
         nqubits=nqubits, delta=0.5, backend=vqe_backend
     )
-    vqe = VQE(
-        build_circuit(
-            nqubits=nqubits,
-            nlayers=nlayers,
-        ),
+
+    # construct circuit from parsed ansatz name
+    circ = getattr(ansatze, config["ansatz"])(config["nqubits"], config["nlayers"])
+
+    vqe = ansatze.VQE(
+        circuit=circ,
         hamiltonian=hamiltonian,
     )
     vqe.circuit.set_parameters(params)

main.py

@@ -16,7 +16,7 @@
 )
 
 # boostvqe's
-from boostvqe.ansatze import build_circuit
+from boostvqe import ansatze
 from boostvqe.plotscripts import plot_gradients, plot_loss
 from boostvqe.training_utils import Model, vqe_loss
 from boostvqe.utils import (
@@ -59,10 +59,12 @@ def main(args):
     path = pathlib.Path(create_folder(generate_path(args)))
     ham = getattr(Model, args.hamiltonian)(args.nqubits)
     target_energy = np.real(np.min(np.asarray(ham.eigenvalues())))
-    circ0 = build_circuit(
-        nqubits=args.nqubits,
-        nlayers=args.nlayers,
-    )
+
+    # construct circuit from parsed ansatz name
+    circ0 = getattr(ansatze, args.ansatz)(args.nqubits, args.nlayers)
+
+    logging.info(circ0.draw())
+
     circ = circ0.copy(deep=True)
     backend = ham.backend
     zero_state = backend.zero_state(args.nqubits)
@@ -326,5 +328,10 @@ def main(args):
         type=int,
         help="number of shots",
     )
+    parser.add_argument(
+        "--ansatz",
+        type=str,
+        help="Parametric quantum circuit ansatz. It can be hw_preserving or hdw_efficient",
+    )
     args = parser.parse_args()
     main(args)

single_commutator_boosting.py

@@ -2,9 +2,12 @@
 import time
 import argparse
 import logging
+import copy
 import pathlib
 
 import numpy as np
+from scipy.optimize import minimize
+
 import qibo
 from qibo import hamiltonians, set_backend
 from qibo.backends import construct_backend
@@ -14,18 +17,39 @@
     DoubleBracketIteration,
 )
 
-from boostvqe.ansatze import VQE, build_circuit
+
+from boostvqe import ansatze
+from boostvqe.training_utils import Model
+
 from boostvqe.utils import (
     OPTIMIZATION_FILE,
     PARAMS_FILE,
-    build_circuit_RBS, 
     apply_dbi_steps, 
     rotate_h_with_vqe,
 )
 
 logging.basicConfig(level=logging.INFO)
 qibo.set_backend("numpy")
 
+def cost_function(s, dbi, zero_state):
+    """Compute the final energy we get applying step `s` to DBI."""
+    original_h = copy.deepcopy(dbi.h)
+    dbi(step=s, d=dbi.diagonal_h_matrix)
+    cost = dbi.h.expectation(np.asarray(zero_state))
+    dbi.h = original_h
+    return cost
+
+def optimize_s(dbi, zero_state):
+    """Optimize DBI stepsize using Powell minimizer."""
+    opt_result = minimize(
+        cost_function, 
+        x0=[0.01], 
+        args=(dbi, zero_state),
+        method="Powell",
+        options={"maxiter": 200}
+    )
+    return opt_result
+
 def main(args):
     path = pathlib.Path(args.path)
     dump_path = (
@@ -48,20 +72,19 @@ def main(args):
     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_RBS(
-            nqubits=nqubits,
-            nlayers=nlayers,
-        ),
+    hamiltonian = getattr(Model, config["hamiltonian"])(config["nqubits"])
+
+    # construct circuit from parsed ansatz name
+    circ = getattr(ansatze, config["ansatz"])(config["nqubits"], config["nlayers"])
+
+    vqe = ansatze.VQE(
+        circuit=circ,
         hamiltonian=hamiltonian,
     )
     print(vqe.circuit.draw())
     vqe.circuit.set_parameters(params)
 
-    zero_state = hamiltonian.backend.zero_state(config["nqubits"])
+    zero_state = np.asarray(hamiltonian.backend.zero_state(config["nqubits"]))
     target_energy = np.min(np.array(hamiltonian.eigenvalues()))
 
     # set target parameters into the VQE
@@ -76,6 +99,7 @@ def main(args):
         nqubits, matrix=new_hamiltonian_matrix
     )
 
+
     dbi = DoubleBracketIteration(
         hamiltonian=new_hamiltonian,
         mode=DoubleBracketGeneratorType.single_commutator,
@@ -85,19 +109,20 @@ def main(args):
     energy_h0 = float(dbi.h.expectation(np.array(zero_state_t)))
     fluctuations_h0 = float(dbi.h.energy_fluctuation(zero_state_t))
 
+    energies = []
 
-    dbi_results = apply_dbi_steps(
-        dbi=dbi,
-        nsteps=args.steps,
-        #optimize_step=True,
-        stepsize=0.01,
-    )
-
-    dbi_energies = dbi_results[1]
+    for s in range(args.steps):
+        logging.info(f"Optimizing step {s+1}")
+        s = optimize_s(dbi, zero_state).x
+        dbi(step=s, d=dbi.diagonal_h_matrix)
+        this_energy = dbi.h.expectation(zero_state)
+        logging.info(f"Best found energy: {this_energy}")
+        energies.append(this_energy)
+
     dict_results = {
         "VQE energy": float(ene1),
         "Ground state": float(target_energy),
-        "dbi_energies": dbi_energies
+        "dbi_energies": energies
     }
 
     (dump_path / "boosting_data.json").write_text(json.dumps(dict_results, indent=4))

src/boostvqe/ansatze.py

@@ -5,8 +5,18 @@
 
 from boostvqe.training_utils import vqe_loss
 
+def connect_qubits(circuit, jumpsize=1, start_from=0):
+    def get_circular_index(n, index):
+        circular_index = index % n
+        return circular_index
+    for q in range(start_from, circuit.nqubits, jumpsize+1):
+        ctrl_index = q
+        targ_index = get_circular_index(circuit.nqubits, q+jumpsize)
+        circuit.add(gates.RBS(q0=ctrl_index, q1=targ_index, theta=0.))
+    return circuit
+
 
-def build_circuit(nqubits, nlayers):
+def hdw_efficient(nqubits, nlayers):
     """Build qibo's aavqe example circuit."""
 
     circuit = Circuit(nqubits)
@@ -21,7 +31,42 @@ def build_circuit(nqubits, nlayers):
     circuit.add(gates.RY(q, theta=0) for q in range(nqubits))
 
     return circuit
+
+
+def hw_preserving(nqubits, nlayers=1):
+
+    if nqubits%2 != 0:
+        raise_error(
+            ValueError,
+            "To use this ansatz please be sure number of qubits is even."
+    )
+    c = Circuit(nqubits)
 
+    for q in range(int(nqubits/2)):
+        c.add(gates.X(q))
+
+    for _ in range(int(nlayers)):
+        c = connect_qubits(c, jumpsize=1, start_from=0)
+        c = connect_qubits(c, jumpsize=1, start_from=1)
+        c = connect_qubits(c, jumpsize=2, start_from=0)
+        c = connect_qubits(c, jumpsize=2, start_from=1)
+        c = connect_qubits(c, jumpsize=2, start_from=3)
+
+    return c
+
+def su2_preserving(nqubits, nlayers):
+    """SU2 invariant circuit."""
+    c = Circuit(nqubits)
+    for _ in range(nlayers):
+        for q in range(1, nqubits, 2):
+            c.add(gates.RXX(q0=q, q1=(q+1)%nqubits, theta=0.))
+            c.add(gates.RYY(q0=q, q1=(q+1)%nqubits, theta=0.))
+            c.add(gates.RZZ(q0=q, q1=(q+1)%nqubits, theta=0.))
+        for q in range(0, nqubits, 2):
+            c.add(gates.RXX(q0=q, q1=(q+1)%nqubits, theta=0.))
+            c.add(gates.RYY(q0=q, q1=(q+1)%nqubits, theta=0.))
+            c.add(gates.RZZ(q0=q, q1=(q+1)%nqubits, theta=0.))
+    return c
 
 def compute_gradients(parameters, circuit, hamiltonian):
     """

src/boostvqe/training_utils.py

@@ -20,19 +20,22 @@
 
 DEFAULT_DELTA = 0.5
 """Default `delta` value of XXZ Hamiltonian"""
-DEFAULT_DELTAS = [0.5, 0.5]
+DEFAULT_DELTAS = [0.0, 2.0]
+TLFIM_h = [1.0, 2.0]
+J1J2_h = [1.0, 0.2]
 
 
 class Model(Enum):
     XXZ = lambda nqubits: hamiltonians.XXZ(
         nqubits=nqubits, delta=DEFAULT_DELTA, dense=False
     )
-    XYZ = lambda nqubits: XYZ(nqubits=nqubits, delta=DEFAULT_DELTAS, dense=False)
+    XYZ = lambda nqubits: XYZ(nqubits=nqubits, delta=[0.5, 0.5], dense=False)
     TFIM = lambda nqubits: hamiltonians.TFIM(nqubits=nqubits, h=nqubits, dense=False)
-    TLFIM = lambda nqubits: TLFIM(nqubits=nqubits, h=[nqubits, nqubits], dense=False)
+    TLFIM = lambda nqubits: TLFIM(nqubits=nqubits, h=TLFIM_h, dense=False)
+    J1J2 = lambda nqubits: J1J2(nqubits=nqubits, h=J1J2_h, dense=False)
 
 
-def TLFIM(nqubits, h=[0.5, 0.5], dense=True, backend=None):
+def TLFIM(nqubits, h=TLFIM_h, dense=True, backend=None):
     """Transverse and longitudinal field Ising model with periodic boundary conditions.
 
     .. math::
@@ -45,7 +48,7 @@ def TLFIM(nqubits, h=[0.5, 0.5], dense=True, backend=None):
             :class:`qibo.core.hamiltonians.Hamiltonian`, otherwise it creates
             a :class:`qibo.core.hamiltonians.SymbolicHamiltonian`.
     """
-    h = [nqubits, nqubits]
+    print(h)
     if nqubits < 2:
         raise_error(ValueError, "Number of qubits must be larger than one.")
     if dense:
@@ -55,6 +58,7 @@ def TLFIM(nqubits, h=[0.5, 0.5], dense=True, backend=None):
             if h[m] != 0:
                 condition = lambda i, j: i == j % nqubits
                 ham += h[m] * _build_spin_model(nqubits, matrix, condition)
+        ham *= -1
         return Hamiltonian(nqubits, ham, backend=backend)
 
     matrix = -(
@@ -69,6 +73,38 @@ def TLFIM(nqubits, h=[0.5, 0.5], dense=True, backend=None):
     return ham
 
 
+def J1J2(nqubits, h=J1J2_h, dense=True, backend=None):
+    """Heisenberg J1-J2 model."""
+    logging.info("Building the J1-J2 Heisenberg model.")
+    if nqubits < 3:
+        raise_error(ValueError, "Number of qubits must be larger than two.")
+    if dense:
+        condition_1 = lambda i, j: i in {j % nqubits, (j + 1) % nqubits}
+        hx = _build_spin_model(nqubits, matrices.X, condition_1)
+        hy = _build_spin_model(nqubits, matrices.Y, condition_1)
+        hz = _build_spin_model(nqubits, matrices.Z, condition_1)
+        condition_2 = lambda i, j: i in {j % nqubits, (j + 2) % nqubits}
+        hx2 = _build_spin_model(nqubits, matrices.X, condition_2)
+        hy2 = _build_spin_model(nqubits, matrices.Y, condition_2)
+        hz2 = _build_spin_model(nqubits, matrices.Z, condition_2)
+        matrix = h[0] * (hx + hy + hz) + h[1] * (hx2 + hy2 + hz2)
+        return Hamiltonian(nqubits, matrix, backend=backend)
+
+    hx = multikron([matrices.X, matrices.X])
+    hy = multikron([matrices.Y, matrices.Y])
+    hz = multikron([matrices.Z, matrices.Z])
+    matrix_1 = h[0] * (hx + hy + hz)
+    matrix_2 = h[1] * (hx + hy + hz)
+    terms = [HamiltonianTerm(matrix_1, i, i + 1) for i in range(nqubits - 1)]
+    terms.extend([HamiltonianTerm(matrix_2, i, i + 2) for i in range(nqubits - 2)])
+    terms.append(HamiltonianTerm(matrix_1, nqubits - 1, 0))
+    terms.append(HamiltonianTerm(matrix_2, nqubits - 2, 0))
+    terms.append(HamiltonianTerm(matrix_2, nqubits - 1, 1))
+    ham = SymbolicHamiltonian(backend=backend)
+    ham.terms = terms
+    return ham
+
+
 def XYZ(nqubits, deltas=[0.5, 0.5], dense=True, backend=None):
     """XYZ model with periodic boundary conditions.