VQE with LabOne Q and Qiskit¶

The following example demonstrates a simple end-to-end implementation of the variational quantum eigensolver (VQE) with LabOne Q.

Qiskit is used to generate the quantum circuits to prepare a parameterized state ansatz and to perform the transformations needed to compute the Hamiltonian expectation value.

The circuits are then exported as OpenQASM programs and subsequently transformed into a LabOne Q experiment. The compiled experiment is then used in an example implementation of a VQE objective function that updates the parameters of the state ansatz, executes the experiment, and returns the expectation value. Finally, this objective function is then minimized with the help of an external optimizer.

Note, there are many ways to improve the efficiency of the VQE algorithm. The comparatively simple implementation demonstrated in this example serves to illustrate basic concepts and might serve as a starting point for users to optimize its performance.

Imports¶

In [1]:

# Numpy:
import numpy as np
from laboneq._utils import id_generator
from laboneq.contrib.example_helpers.generate_device_setup import (
    generate_device_setup_qubits,
)
from laboneq.openqasm3.gate_store import GateStore

# LabOne Q:
from laboneq.simple import *

# Qiskit:
from qiskit import QuantumCircuit, qasm3
from qiskit.circuit.library import EfficientSU2
from qiskit.quantum_info import SparsePauliOp

# scipy optimizer
from scipy.optimize import minimize

# Enable/disable emulation mode
use_emulation = True

# Numpy: import numpy as np from laboneq._utils import id_generator from laboneq.contrib.example_helpers.generate_device_setup import ( generate_device_setup_qubits, ) from laboneq.openqasm3.gate_store import GateStore # LabOne Q: from laboneq.simple import * # Qiskit: from qiskit import QuantumCircuit, qasm3 from qiskit.circuit.library import EfficientSU2 from qiskit.quantum_info import SparsePauliOp # scipy optimizer from scipy.optimize import minimize # Enable/disable emulation mode use_emulation = True

Hamiltonian (Qiskit)¶

The variational quantum eigensolver (VQE) is an algorithm to approximate the ground state $\psi_{0}$ with eigenvalue $\epsilon_{0}$ of a problem described by a Hamiltonian $\hat{H}$. VQE achieves this by minimizing the expectation value of $\hat{H}$

$\epsilon_{\theta} = \langle\psi_{\theta}|\hat{H}|\psi_{\theta}\rangle$

with respect to a trial ansatz $\psi_{\theta}$ parameterized by $\theta$. Because of the variational principle $\epsilon_{\theta}$ is always an upper bound for $\epsilon_{0}$. A successful minimization of $\epsilon_{\theta}$ will thus provide an increasingly accurate approximation to the sought solution.

VQE requires the repeated computation of expectation values. This can be done efficiently for Hamiltonians which have the form of a (sparse) series,

$\hat{H} = \sum\limits_{i}c_{i}\ \hat{P}_{i}$

of multi-qubit Pauli operators $\hat{P}_{i} \in \{\hat{I}, \hat{X}, \hat{Y}, \hat{Z}\}^{\otimes n}$ with coefficients $c_{i}$.

Quantum computing software packages like Qiskit have dedicated modules that construct Hamiltonians in this form for various quantum computational applications. For this example, we use Qiskit's SparsePauliOp class to define such a Hamiltonian explicitly.

In [2]:

ham = SparsePauliOp.from_list(
    [("YZ", 0.3980), ("ZI", -0.3980), ("ZZ", -1.0113), ("XX", 0.1810)]
)

ham = SparsePauliOp.from_list( [("YZ", 0.3980), ("ZI", -0.3980), ("ZZ", -1.0113), ("XX", 0.1810)] )

Parameterized wavefunction ansatz (Qiskit)¶

In this section, we construct a parameterized circuit that prepares the trial state $\psi_{\theta}$. Qiskit's EfficientSU2 class allows us to define an ansatz for this state. For two qubits this results in four layers. Each layer applies two parameterized SU2 gates (here Rx and Ry) to each qubit, respectively, followed by a cnot gate to create an entangled state.

Note, that we use the Rx,Ry pair of SU2 gates over the default Ry, Rz pair since optimizing over virtual z-gate angles would require replacing phase increments in compiled LabOne Q experiments, which at the time of writing is still an experimental feature.

The resulting circuit is appended to the initial state reset.

In [3]:

trial_state = EfficientSU2(ham.num_qubits, su2_gates=["rx", "ry"], insert_barriers=True)
trial_state.barrier()

trial_state.decompose().draw(style="iqp")

trial_state = EfficientSU2(ham.num_qubits, su2_gates=["rx", "ry"], insert_barriers=True) trial_state.barrier() trial_state.decompose().draw(style="iqp")

Out[3]:

     ┌──────────┐┌──────────┐ ░       ░ ┌──────────┐┌──────────┐ ░       ░ »
q_0: ┤ Rx(θ[0]) ├┤ Ry(θ[2]) ├─░───■───░─┤ Rx(θ[4]) ├┤ Ry(θ[6]) ├─░───■───░─»
     ├──────────┤├──────────┤ ░ ┌─┴─┐ ░ ├──────────┤├──────────┤ ░ ┌─┴─┐ ░ »
q_1: ┤ Rx(θ[1]) ├┤ Ry(θ[3]) ├─░─┤ X ├─░─┤ Rx(θ[5]) ├┤ Ry(θ[7]) ├─░─┤ X ├─░─»
     └──────────┘└──────────┘ ░ └───┘ ░ └──────────┘└──────────┘ ░ └───┘ ░ »
«     ┌──────────┐┌───────────┐ ░       ░ ┌───────────┐┌───────────┐ ░ 
«q_0: ┤ Rx(θ[8]) ├┤ Ry(θ[10]) ├─░───■───░─┤ Rx(θ[12]) ├┤ Ry(θ[14]) ├─░─
«     ├──────────┤├───────────┤ ░ ┌─┴─┐ ░ ├───────────┤├───────────┤ ░ 
«q_1: ┤ Rx(θ[9]) ├┤ Ry(θ[11]) ├─░─┤ X ├─░─┤ Rx(θ[13]) ├┤ Ry(θ[15]) ├─░─
«     └──────────┘└───────────┘ ░ └───┘ ░ └───────────┘└───────────┘ ░ 

The above trial state $|\psi_{\theta}\rangle$ is parameterized by $\theta=(\theta_{0}, \theta_{1}, \dots, \theta_{15})$.

Transformation to Pauli eigenbasis (Qiskit)¶

After preparing the parameterized trial state $|\psi_{\theta}\rangle$ we need to evaluate the Hamiltonian expectation value for this state. We note the following aspects:

• Expectation values are computed for each Pauli term $\hat{P}_{i}$ individually. For the sake of simplicity, we use circuits and ultimately separate LabOne Q experiments. We note here, that optimizing the approach discussed here is
• We can only obtain information from a quantum circuit by measuring states in the computational basis $\left\{|0\rangle, |1\rangle\right\}$.
• The computational basis of a qubit corresponds to the eigenfunctions of the Pauli $\hat{Z}$ operator acting on that qubit. All occurrences of $\hat{X}$ or $\hat{Y}$ in a Pauli term are therefore transformed into $\hat{Z}$ operators by applying a transformation $\hat{U}$ to the trial state $|\psi_{\theta}\rangle$.

In [4]:

vqe_terms = []
for p in ham.paulis:
    P_i = p.to_label()
    u_transform = QuantumCircuit(ham.num_qubits, ham.num_qubits)

    for q in range(ham.num_qubits):
        if P_i[q] == "X":
            u_transform.h(q)

        elif P_i[q] == "Y":
            u_transform.sdg(q)
            u_transform.h(q)

    # u_transform.barrier()

    # for q in range(ham.num_qubits):
    #     u_transform.measure(q, q)

    vqe_terms.append(trial_state & u_transform)

vqe_terms = [] for p in ham.paulis: P_i = p.to_label() u_transform = QuantumCircuit(ham.num_qubits, ham.num_qubits) for q in range(ham.num_qubits): if P_i[q] == "X": u_transform.h(q) elif P_i[q] == "Y": u_transform.sdg(q) u_transform.h(q) # u_transform.barrier() # for q in range(ham.num_qubits): # u_transform.measure(q, q) vqe_terms.append(trial_state & u_transform)

We can now inspect the final circuits of the VQE algorithm.

In [5]:

vqe_terms[0].decompose().draw(style="iqp")
# vqe_terms[1].decompose().draw(style="iqp")
# vqe_terms[2].decompose().draw(style="iqp")
# vqe_terms[3].decompose().draw(style="iqp")

vqe_terms[0].decompose().draw(style="iqp") # vqe_terms[1].decompose().draw(style="iqp") # vqe_terms[2].decompose().draw(style="iqp") # vqe_terms[3].decompose().draw(style="iqp")

Out[5]:

     ┌──────────┐┌──────────┐ ░       ░ ┌──────────┐┌──────────┐ ░       ░ »
q_0: ┤ Rx(θ[0]) ├┤ Ry(θ[2]) ├─░───■───░─┤ Rx(θ[4]) ├┤ Ry(θ[6]) ├─░───■───░─»
     ├──────────┤├──────────┤ ░ ┌─┴─┐ ░ ├──────────┤├──────────┤ ░ ┌─┴─┐ ░ »
q_1: ┤ Rx(θ[1]) ├┤ Ry(θ[3]) ├─░─┤ X ├─░─┤ Rx(θ[5]) ├┤ Ry(θ[7]) ├─░─┤ X ├─░─»
     └──────────┘└──────────┘ ░ └───┘ ░ └──────────┘└──────────┘ ░ └───┘ ░ »
c: 2/══════════════════════════════════════════════════════════════════════»
                                                                           »
«     ┌──────────┐┌───────────┐ ░       ░ ┌───────────┐┌───────────┐ ░ »
«q_0: ┤ Rx(θ[8]) ├┤ Ry(θ[10]) ├─░───■───░─┤ Rx(θ[12]) ├┤ Ry(θ[14]) ├─░─»
«     ├──────────┤├───────────┤ ░ ┌─┴─┐ ░ ├───────────┤├───────────┤ ░ »
«q_1: ┤ Rx(θ[9]) ├┤ Ry(θ[11]) ├─░─┤ X ├─░─┤ Rx(θ[13]) ├┤ Ry(θ[15]) ├─░─»
«     └──────────┘└───────────┘ ░ └───┘ ░ └───────────┘└───────────┘ ░ »
«c: 2/═════════════════════════════════════════════════════════════════»
«                                                                      »
«     ┌──────────┐┌─────────┐
«q_0: ┤ U1(-π/2) ├┤ U2(0,π) ├
«     └──────────┘└─────────┘
«q_1: ───────────────────────
«                            
«c: 2/═══════════════════════
«                            

Conversion to OpenQASM code (Qiskit)¶

We can use Qiskit's qasm3 functionality to turn the circuits defined above into qasm programs.

In [6]:

qasm_programs = [qasm3.dumps(c.decompose()) for c in vqe_terms]

print(qasm_programs[0])

qasm_programs = [qasm3.dumps(c.decompose()) for c in vqe_terms] print(qasm_programs[0])

OPENQASM 3.0;
include "stdgates.inc";
input float[64] _θ_0_;
input float[64] _θ_1_;
input float[64] _θ_2_;
input float[64] _θ_3_;
input float[64] _θ_4_;
input float[64] _θ_5_;
input float[64] _θ_6_;
input float[64] _θ_7_;
input float[64] _θ_8_;
input float[64] _θ_9_;
input float[64] _θ_10_;
input float[64] _θ_11_;
input float[64] _θ_12_;
input float[64] _θ_13_;
input float[64] _θ_14_;
input float[64] _θ_15_;
bit[2] c;
qubit[2] q;
rx(_θ_0_) q[0];
ry(_θ_2_) q[0];
rx(_θ_1_) q[1];
ry(_θ_3_) q[1];
barrier q[0], q[1];
cx q[0], q[1];
barrier q[0], q[1];
rx(_θ_4_) q[0];
ry(_θ_6_) q[0];
rx(_θ_5_) q[1];
ry(_θ_7_) q[1];
barrier q[0], q[1];
cx q[0], q[1];
barrier q[0], q[1];
rx(_θ_8_) q[0];
ry(_θ_10_) q[0];
rx(_θ_9_) q[1];
ry(_θ_11_) q[1];
barrier q[0], q[1];
cx q[0], q[1];
barrier q[0], q[1];
rx(_θ_12_) q[0];
ry(_θ_14_) q[0];
rx(_θ_13_) q[1];
ry(_θ_15_) q[1];
barrier q[0], q[1];
u1(-pi/2) q[0];
u2(0, pi) q[0];

Note, that the parameters now appear as input variables _θ_0_ to _θ_15_.

Parameter Initialization¶

When later assigning values to these parameters, we will need to provide matching parameter names.

Also parameter updates to the objective function are achieved by replacing the pulses containing these parameters. For this we need to be able to identify such parameter dependent pulses in the compiled experiments.

To achieve this, we define the ScalarParameter class that holds the parameter value as well as its name.

In [7]:

class ScalarParameter:
    def __init__(self, val, identifier):
        self.val = val
        self.identifier = identifier

class ScalarParameter: def __init__(self, val, identifier): self.val = val self.identifier = identifier

We initialize the parameters with random values in the interval $[0, 1)$.

In [8]:

np.random.seed(seed=42)
theta = [
    ScalarParameter(np.random.random(), f"_θ_{i}_")
    for i in range(trial_state.num_parameters)
]

np.random.seed(seed=42) theta = [ ScalarParameter(np.random.random(), f"_θ_{i}_") for i in range(trial_state.num_parameters) ]

Gate Definitions (LabOne Q)¶

The construction of circuits and qasm programs made use of abstract gate definitions. To run the VQE algorithm on actual hardware, we need to define the gates $U_1$, $U_2$, $CX$, $R_x$, and $R_y$ as well as reset and measurement operations in terms of pulse sequences in LabOne Q.

Note that the definitions below are for demonstration purposes. Actual gate and pulse implementations are typically hardware-specific might even vary for individual qubits.

In [9]:

def u1(qubit: Qubit):
    """Return a parameterized U1 gate for the specified qubit.

    The gate is a function that takes the angle to rotate and
    returns a LabOne Q section that performs the rotation.
    """

    def u1_gate(angle):
        """Equivalent to Rz(theta) as U1(theta) differs from it only by a global phase which is ignored.

        Theta is in radians - implements a virtual z-gate.
        """
        gate = Section(uid=id_generator(f"p_{qubit.uid}_u1_{int(180 * angle / np.pi)}"))
        gate.play(
            signal=qubit.signals["drive"],
            pulse=None,
            increment_oscillator_phase=angle,
        )
        return gate

    return u1_gate

def u1(qubit: Qubit): """Return a parameterized U1 gate for the specified qubit. The gate is a function that takes the angle to rotate and returns a LabOne Q section that performs the rotation. """ def u1_gate(angle): """Equivalent to Rz(theta) as U1(theta) differs from it only by a global phase which is ignored. Theta is in radians - implements a virtual z-gate. """ gate = Section(uid=id_generator(f"p_{qubit.uid}_u1_{int(180 * angle / np.pi)}")) gate.play( signal=qubit.signals["drive"], pulse=None, increment_oscillator_phase=angle, ) return gate return u1_gate

In [10]:

def u2(qubit: Qubit):
    """Return a parameterized U2 gate for the specified qubit.

    The gate is a function that takes two angles to rotate and
    returns a LabOne Q section that performs the rotation.
    """

    def u2_gate(phi, lam):
        """U2(ϕ,λ) = RZ(ϕ).RY(π/2).RZ(λ) - following definition in qiskit.circuit.library.U2Gate.

        Variables phi and lambda are in radians and serve to implement
        the two virtual z-gates.
        """
        pulse = pulse_library.drag(
            uid=f"{qubit.uid}_ry_in_u2",
            length=qubit.parameters.user_defined["pulse_length"],
            amplitude=qubit.parameters.user_defined["amplitude_pi"] / 2,
        )
        gate = Section(
            uid=id_generator(
                f"p_{qubit.uid}_u2_{int(180 * phi / np.pi)}_{int(180 * lam / np.pi)}"
            )
        )
        gate.play(
            signal=qubit.signals["drive"],
            pulse=None,
            increment_oscillator_phase=lam,
        )
        gate.play(
            signal=qubit.signals["drive"],
            pulse=pulse,
            phase=np.pi / 2,
        )
        gate.play(
            signal=qubit.signals["drive"],
            pulse=None,
            increment_oscillator_phase=phi,
        )
        return gate

    return u2_gate

def u2(qubit: Qubit): """Return a parameterized U2 gate for the specified qubit. The gate is a function that takes two angles to rotate and returns a LabOne Q section that performs the rotation. """ def u2_gate(phi, lam): """U2(ϕ,λ) = RZ(ϕ).RY(π/2).RZ(λ) - following definition in qiskit.circuit.library.U2Gate. Variables phi and lambda are in radians and serve to implement the two virtual z-gates. """ pulse = pulse_library.drag( uid=f"{qubit.uid}_ry_in_u2", length=qubit.parameters.user_defined["pulse_length"], amplitude=qubit.parameters.user_defined["amplitude_pi"] / 2, ) gate = Section( uid=id_generator( f"p_{qubit.uid}_u2_{int(180 * phi / np.pi)}_{int(180 * lam / np.pi)}" ) ) gate.play( signal=qubit.signals["drive"], pulse=None, increment_oscillator_phase=lam, ) gate.play( signal=qubit.signals["drive"], pulse=pulse, phase=np.pi / 2, ) gate.play( signal=qubit.signals["drive"], pulse=None, increment_oscillator_phase=phi, ) return gate return u2_gate

In [11]:

def cx(control: Qubit, target: Qubit):
    """Return a controlled X gate for the specified control and target qubits.

    The CX gate function takes no arguments and returns a LabOne Q section that performs
    the controlled X gate.
    """

    def cx_gate():
        cx_id = f"cx_{control.uid}_{target.uid}"

        gate = Section(uid=id_generator(cx_id))

        # define X pulses for target and control
        x180_pulse_control = pulse_library.drag(
            uid=f"{control.uid}_x180",
            length=control.parameters.user_defined["pulse_length"],
            amplitude=control.parameters.user_defined["amplitude_pi"],
        )
        x180_pulse_target = pulse_library.drag(
            uid=f"{target.uid}_x180",
            length=target.parameters.user_defined["pulse_length"],
            amplitude=target.parameters.user_defined["amplitude_pi"],
        )

        # define cancellation pulses for target and control
        cancellation_control_n = pulse_library.gaussian_square(uid="CR-")
        cancellation_control_p = pulse_library.gaussian_square(uid="CR+")
        cancellation_target_p = pulse_library.gaussian_square(uid="q1+")
        cancellation_target_n = pulse_library.gaussian_square(uid="q1-")

        # play X pulses on both target and control
        x180_both = Section(uid=id_generator(f"{cx_id}_x_both"))
        x180_both.play(signal=control.signals["drive"], pulse=x180_pulse_control)
        x180_both.play(signal=target.signals["drive"], pulse=x180_pulse_target)
        gate.add(x180_both)

        # First cross-resonance component
        cancellation_p = Section(
            uid=id_generator(f"{cx_id}_canc_p"), play_after=x180_both
        )
        cancellation_p.play(signal=target.signals["drive"], pulse=cancellation_target_p)
        cancellation_p.play(
            signal=control.signals["flux"], pulse=cancellation_control_n
        )
        gate.add(cancellation_p)

        # play X pulse on control
        x180_control = Section(
            uid=id_generator(f"{cx_id}_x_q0"), play_after=cancellation_p
        )
        x180_control.play(signal=control.signals["drive"], pulse=x180_pulse_control)
        gate.add(x180_control)

        # Second cross-resonance component
        cancellation_n = Section(
            uid=id_generator(f"cx_{cx_id}_canc_n"), play_after=x180_control
        )
        cancellation_n.play(signal=target.signals["drive"], pulse=cancellation_target_n)
        cancellation_n.play(
            signal=control.signals["flux"], pulse=cancellation_control_p
        )
        gate.add(cancellation_n)

        return gate

    return cx_gate

def cx(control: Qubit, target: Qubit): """Return a controlled X gate for the specified control and target qubits. The CX gate function takes no arguments and returns a LabOne Q section that performs the controlled X gate. """ def cx_gate(): cx_id = f"cx_{control.uid}_{target.uid}" gate = Section(uid=id_generator(cx_id)) # define X pulses for target and control x180_pulse_control = pulse_library.drag( uid=f"{control.uid}_x180", length=control.parameters.user_defined["pulse_length"], amplitude=control.parameters.user_defined["amplitude_pi"], ) x180_pulse_target = pulse_library.drag( uid=f"{target.uid}_x180", length=target.parameters.user_defined["pulse_length"], amplitude=target.parameters.user_defined["amplitude_pi"], ) # define cancellation pulses for target and control cancellation_control_n = pulse_library.gaussian_square(uid="CR-") cancellation_control_p = pulse_library.gaussian_square(uid="CR+") cancellation_target_p = pulse_library.gaussian_square(uid="q1+") cancellation_target_n = pulse_library.gaussian_square(uid="q1-") # play X pulses on both target and control x180_both = Section(uid=id_generator(f"{cx_id}_x_both")) x180_both.play(signal=control.signals["drive"], pulse=x180_pulse_control) x180_both.play(signal=target.signals["drive"], pulse=x180_pulse_target) gate.add(x180_both) # First cross-resonance component cancellation_p = Section( uid=id_generator(f"{cx_id}_canc_p"), play_after=x180_both ) cancellation_p.play(signal=target.signals["drive"], pulse=cancellation_target_p) cancellation_p.play( signal=control.signals["flux"], pulse=cancellation_control_n ) gate.add(cancellation_p) # play X pulse on control x180_control = Section( uid=id_generator(f"{cx_id}_x_q0"), play_after=cancellation_p ) x180_control.play(signal=control.signals["drive"], pulse=x180_pulse_control) gate.add(x180_control) # Second cross-resonance component cancellation_n = Section( uid=id_generator(f"cx_{cx_id}_canc_n"), play_after=x180_control ) cancellation_n.play(signal=target.signals["drive"], pulse=cancellation_target_n) cancellation_n.play( signal=control.signals["flux"], pulse=cancellation_control_p ) gate.add(cancellation_n) return gate return cx_gate

In [12]:

def measurement(qubit: Qubit):
    """Return a measurement operation of the specified qubit.

    The operation is a function that takes the measurement handle (a string)
    and returns a LabOne Q section that performs the measurement.
    """

    def measurement_gate(handle: str):
        """Perform a measurement.

        Handle is the name of where to store the measurement result. E.g. "meas[0]".
        """
        measure_pulse = pulse_library.gaussian_square(
            uid=f"{qubit.uid}_readout_pulse",
            length=qubit.parameters.user_defined["readout_length"],
            amplitude=qubit.parameters.user_defined["readout_amplitude"],
        )
        integration_kernel = pulse_library.const(
            uid=f"{qubit.uid}_integration_kernel",
            length=qubit.parameters.user_defined["readout_length"],
        )
        gate = Section(uid=id_generator(f"meas_{qubit.uid}_{handle}"))
        gate.reserve(signal=qubit.signals["drive"])
        gate.play(signal=qubit.signals["measure"], pulse=measure_pulse)
        gate.acquire(
            signal=qubit.signals["acquire"],
            handle=handle,
            kernel=integration_kernel,
        )
        return gate

    return measurement_gate

def measurement(qubit: Qubit): """Return a measurement operation of the specified qubit. The operation is a function that takes the measurement handle (a string) and returns a LabOne Q section that performs the measurement. """ def measurement_gate(handle: str): """Perform a measurement. Handle is the name of where to store the measurement result. E.g. "meas[0]". """ measure_pulse = pulse_library.gaussian_square( uid=f"{qubit.uid}_readout_pulse", length=qubit.parameters.user_defined["readout_length"], amplitude=qubit.parameters.user_defined["readout_amplitude"], ) integration_kernel = pulse_library.const( uid=f"{qubit.uid}_integration_kernel", length=qubit.parameters.user_defined["readout_length"], ) gate = Section(uid=id_generator(f"meas_{qubit.uid}_{handle}")) gate.reserve(signal=qubit.signals["drive"]) gate.play(signal=qubit.signals["measure"], pulse=measure_pulse) gate.acquire( signal=qubit.signals["acquire"], handle=handle, kernel=integration_kernel, ) return gate return measurement_gate

In [13]:

def reset(qubit: Qubit):
    """Reset the specified qubit to the ground state with the supplied reset pulse.

    The reset gate function takes no arguments and returns a LabOne Q section that performs
    the reset.
    """

    def reset_gate():
        reset_pulse = pulse_library.drag(
            uid=f"{qubit.uid}_reset",
            length=qubit.parameters.user_defined["pulse_length"],
            amplitude=qubit.parameters.user_defined["amplitude_pi"],
        )
        # Reset Section
        reset = Section(uid=f"{qubit.uid}_reset")
        # qubit state readout
        readout = measurement(qubit)(f"{qubit.uid}_reset")
        # delay after measurement
        readout.delay(
            signal=qubit.signals["acquire"],
            time=qubit.parameters.user_defined["reset_delay_length"],
        )
        # real-time feedback, fetching measurement data identified
        # by handle locally from the QA unit of the SHFQC
        match_case = Match(
            uid=f"{qubit.uid}_feedback",
            handle=f"{qubit.uid}_reset",
            play_after=readout,
        )
        # measurement result 0 - ground state
        case_0 = Case(uid=f"{qubit.uid}_0_Case", state=0)
        case_0.play(signal=qubit.signals["drive"], pulse=reset_pulse, amplitude=0.01)
        # measurement result 1 - excited state
        case_1 = Case(uid=f"{qubit.uid}_1_Case", state=1)
        # play x180 pulse
        case_1.play(signal=qubit.signals["drive"], pulse=reset_pulse)
        match_case.add(case_0)
        match_case.add(case_1)

        reset.add(readout)
        reset.add(match_case)
        return reset

    return reset_gate

def reset(qubit: Qubit): """Reset the specified qubit to the ground state with the supplied reset pulse. The reset gate function takes no arguments and returns a LabOne Q section that performs the reset. """ def reset_gate(): reset_pulse = pulse_library.drag( uid=f"{qubit.uid}_reset", length=qubit.parameters.user_defined["pulse_length"], amplitude=qubit.parameters.user_defined["amplitude_pi"], ) # Reset Section reset = Section(uid=f"{qubit.uid}_reset") # qubit state readout readout = measurement(qubit)(f"{qubit.uid}_reset") # delay after measurement readout.delay( signal=qubit.signals["acquire"], time=qubit.parameters.user_defined["reset_delay_length"], ) # real-time feedback, fetching measurement data identified # by handle locally from the QA unit of the SHFQC match_case = Match( uid=f"{qubit.uid}_feedback", handle=f"{qubit.uid}_reset", play_after=readout, ) # measurement result 0 - ground state case_0 = Case(uid=f"{qubit.uid}_0_Case", state=0) case_0.play(signal=qubit.signals["drive"], pulse=reset_pulse, amplitude=0.01) # measurement result 1 - excited state case_1 = Case(uid=f"{qubit.uid}_1_Case", state=1) # play x180 pulse case_1.play(signal=qubit.signals["drive"], pulse=reset_pulse) match_case.add(case_0) match_case.add(case_1) reset.add(readout) reset.add(match_case) return reset return reset_gate

The remaining Rx and Ry gates depend on the parameters θ.

The following function generates a drag pulse from an amplitude parameter.

In [14]:

def pulse_from_parameter(theta):
    val = theta.val % 2.0  # full rotations
    identifier = f"pulse{theta.identifier}"
    return pulse_library.drag(
        uid=identifier, amplitude=val if val <= 1.0 else val - 2.0
    )

def pulse_from_parameter(theta): val = theta.val % 2.0 # full rotations identifier = f"pulse{theta.identifier}" return pulse_library.drag( uid=identifier, amplitude=val if val <= 1.0 else val - 2.0 )

In [15]:

def rx(qubit: Qubit):
    """Return a parameterized Rx gate for the specified qubit.

    The gate is a function that takes the angle to rotate and
    returns a LabOne Q section that performs the rotation.
    """

    def rx_gate(angle: ScalarParameter):
        """Rx(theta).

        Theta is in units of pi - pulse amplitude
        is adjusted according to the chosen angle.
        """
        gate = Section(uid=id_generator(f"p_{qubit.uid}_rx_{int(180 * angle.val)}"))
        gate.play(
            signal=qubit.signals["drive"],
            pulse=pulse_from_parameter(angle),
            amplitude=qubit.parameters.user_defined["amplitude_pi"],
        )
        return gate

    return rx_gate

def rx(qubit: Qubit): """Return a parameterized Rx gate for the specified qubit. The gate is a function that takes the angle to rotate and returns a LabOne Q section that performs the rotation. """ def rx_gate(angle: ScalarParameter): """Rx(theta). Theta is in units of pi - pulse amplitude is adjusted according to the chosen angle. """ gate = Section(uid=id_generator(f"p_{qubit.uid}_rx_{int(180 * angle.val)}")) gate.play( signal=qubit.signals["drive"], pulse=pulse_from_parameter(angle), amplitude=qubit.parameters.user_defined["amplitude_pi"], ) return gate return rx_gate

In [16]:

def ry(qubit: Qubit):
    """Return a parameterized Ry gate for the specified qubit.

    The gate is a function that takes the angle to rotate and
    returns a LabOne Q section that performs the rotation.
    """

    def ry_gate(angle: ScalarParameter):
        """Ry(theta).

        Theta is in units of pi - pulse amplitude
        is adjusted according to the chosen angle.
        """
        gate = Section(uid=id_generator(f"p_{qubit.uid}_ry_{int(180 * angle.val)}"))

        gate.play(
            signal=qubit.signals["drive"],
            pulse=pulse_from_parameter(angle),
            amplitude=qubit.parameters.user_defined["amplitude_pi"],
            phase=np.pi / 2,
        )
        return gate

    return ry_gate

def ry(qubit: Qubit): """Return a parameterized Ry gate for the specified qubit. The gate is a function that takes the angle to rotate and returns a LabOne Q section that performs the rotation. """ def ry_gate(angle: ScalarParameter): """Ry(theta). Theta is in units of pi - pulse amplitude is adjusted according to the chosen angle. """ gate = Section(uid=id_generator(f"p_{qubit.uid}_ry_{int(180 * angle.val)}")) gate.play( signal=qubit.signals["drive"], pulse=pulse_from_parameter(angle), amplitude=qubit.parameters.user_defined["amplitude_pi"], phase=np.pi / 2, ) return gate return ry_gate

Device setup, qubits, and gate store (LabOne Q)¶

Before we can proceed with LabOne Q, we need to define the necessary qubit objects (and their respective parameters). Here we use a helper function to generate the device_setup and qubit objects from a pre-defined parameter set.

In [17]:

# specify the number of qubits you want to use
number_of_qubits = 2

# generate the device setup and the qubit objects using a helper function
device_setup, qubits = generate_device_setup_qubits(
    number_qubits=number_of_qubits,
    pqsc=[{"serial": "DEV10001"}],
    hdawg=[{"serial": "DEV8001", "zsync": 0, "number_of_channels": 8, "options": None}],
    shfqc=[
        {
            "serial": "DEV12001",
            "zsync": 1,
            "number_of_channels": 6,
            "readout_multiplex": 6,
            "options": None,
        }
    ],
    include_flux_lines=True,
    server_host="localhost",
    setup_name=f"my_{number_of_qubits}_tuneable_qubit_setup",
)

q0, q1 = qubits[:2]

# specify the number of qubits you want to use number_of_qubits = 2 # generate the device setup and the qubit objects using a helper function device_setup, qubits = generate_device_setup_qubits( number_qubits=number_of_qubits, pqsc=[{"serial": "DEV10001"}], hdawg=[{"serial": "DEV8001", "zsync": 0, "number_of_channels": 8, "options": None}], shfqc=[ { "serial": "DEV12001", "zsync": 1, "number_of_channels": 6, "readout_multiplex": 6, "options": None, } ], include_flux_lines=True, server_host="localhost", setup_name=f"my_{number_of_qubits}_tuneable_qubit_setup", ) q0, q1 = qubits[:2]

The gate store defines the mapping between logical operations (i.e. those that appear in OpenQASM statements) and the physical operations (i.e. functions that define LabOne Q sections to play) above.

The qubit names provided in this mapping therefore need to match the qubit names of the OpenQASM program.

In [18]:

# Initialize
gate_store = GateStore()

# map LabOne Q qubits to qasm qubits
qubit_map = {"q[0]": device_setup.qubits[0], "q[1]": device_setup.qubits[1]}

# Single qubit gates:
for q in qubit_map:
    gate_store.register_gate_section("rx", (q,), rx(qubit_map[q]))
    gate_store.register_gate_section("ry", (q,), ry(qubit_map[q]))
    gate_store.register_gate_section("u1", (q,), u1(qubit_map[q]))
    gate_store.register_gate_section("u2", (q,), u2(qubit_map[q]))
    gate_store.register_gate_section("measure", (q,), measurement(qubit_map[q]))
    gate_store.register_gate_section("reset", (q,), reset(qubit_map[q]))

# Two qubit gates:
gate_store.register_gate_section("cx", ("q[0]", "q[1]"), cx(q0, q1))

# Initialize gate_store = GateStore() # map LabOne Q qubits to qasm qubits qubit_map = {"q[0]": device_setup.qubits[0], "q[1]": device_setup.qubits[1]} # Single qubit gates: for q in qubit_map: gate_store.register_gate_section("rx", (q,), rx(qubit_map[q])) gate_store.register_gate_section("ry", (q,), ry(qubit_map[q])) gate_store.register_gate_section("u1", (q,), u1(qubit_map[q])) gate_store.register_gate_section("u2", (q,), u2(qubit_map[q])) gate_store.register_gate_section("measure", (q,), measurement(qubit_map[q])) gate_store.register_gate_section("reset", (q,), reset(qubit_map[q])) # Two qubit gates: gate_store.register_gate_section("cx", ("q[0]", "q[1]"), cx(q0, q1))

/tmp/ipykernel_8101/3121741620.py:2: FutureWarning: GateStore is deprecated. Use `OpenQASMTranspiler()` and `QuantumOperations()` instead.
  gate_store = GateStore()

LabOne Q Experiment (LabOne Q)¶

We are now ready to use LabOne Q's exp_from_qasm_list function to convert the OpenQASM programs into a single Experiment instance and note the following

• the initial parameter values are assigned to the inputs argument
• averaging_mode and acquisition_type are set to allow to use the state indices for the computation of the expectation value
• we set the repetition time to the approximative length of the VQE sequence manually
• the do_reset option is used to initialize all qubits to the ground state before each VQE sequence
• the add_measurement option allows to add measurements of all qubit after each VQE sequence

In [19]:

exp = exp_from_qasm_list(
    qasm_programs,
    qubits=qubit_map,
    gate_store=gate_store,
    inputs={t.identifier: t for t in theta},
    count=1024,
    averaging_mode=AveragingMode.SINGLE_SHOT,
    acquisition_type=AcquisitionType.DISCRIMINATION,
    batch_execution_mode="pipeline",
    repetition_time=2.2e-6,
    do_reset=True,
    add_measurement=True,
)

exp = exp_from_qasm_list( qasm_programs, qubits=qubit_map, gate_store=gate_store, inputs={t.identifier: t for t in theta}, count=1024, averaging_mode=AveragingMode.SINGLE_SHOT, acquisition_type=AcquisitionType.DISCRIMINATION, batch_execution_mode="pipeline", repetition_time=2.2e-6, do_reset=True, add_measurement=True, )

/tmp/ipykernel_8101/2518688964.py:1: FutureWarning: `exp_from_qasm_list()` is deprecated. Use `OpenQASMTranspiler.batch_experiment()` instead.
  exp = exp_from_qasm_list(

A Session instance is created and used to compile all experiments.

In [20]:

session = Session(device_setup=device_setup)
session.connect(do_emulation=use_emulation)

cexp = session.compile(exp)

session = Session(device_setup=device_setup) session.connect(do_emulation=use_emulation) cexp = session.compile(exp)

[2024.12.19 17:10:12.381] INFO    Logging initialized from [Default inline config in laboneq.laboneq_logging] logdir is /builds/qccs/laboneq-applications/docs/sources/how-to-guides/sources/04_qasm/laboneq_output/log

[2024.12.19 17:10:12.383] INFO    VERSION: laboneq 2.43.0

[2024.12.19 17:10:12.384] INFO    Connecting to data server at localhost:8004

[2024.12.19 17:10:12.387] INFO    Connected to Zurich Instruments LabOne Data Server version 24.10 at localhost:8004

[2024.12.19 17:10:12.387] WARNING HDAWG:dev8001: Include the device options 'HDAWG8/MF/ME/SKW/PC' in the device setup ('options' field of the 'instruments' list in the device setup descriptor, 'device_options' argument when constructing instrument objects to be added to 'DeviceSetup' instances). This will become a strict requirement in the future.

[2024.12.19 17:10:12.389] WARNING SHFQC/QA:dev12001: Include the device options 'SHFQC/QC6CH' in the device setup ('options' field of the 'instruments' list in the device setup descriptor, 'device_options' argument when constructing instrument objects to be added to 'DeviceSetup' instances). This will become a strict requirement in the future.

[2024.12.19 17:10:12.390] INFO    Configuring the device setup

[2024.12.19 17:10:12.429] INFO    The device setup is configured

[2024.12.19 17:10:12.449] INFO    Resolved modulation type of oscillator 'q0_readout_acquire_osc' on signal '/logical_signal_groups/q0/acquire_line' to SOFTWARE

[2024.12.19 17:10:12.450] INFO    Resolved modulation type of oscillator 'q0_drive_ge_osc' on signal '/logical_signal_groups/q0/drive_line' to HARDWARE

[2024.12.19 17:10:12.450] INFO    Resolved modulation type of oscillator 'q1_readout_acquire_osc' on signal '/logical_signal_groups/q1/acquire_line' to SOFTWARE

[2024.12.19 17:10:12.451] INFO    Resolved modulation type of oscillator 'q1_drive_ge_osc' on signal '/logical_signal_groups/q1/drive_line' to HARDWARE

[2024.12.19 17:10:12.452] INFO    Starting LabOne Q Compiler run...

[2024.12.19 17:10:12.462] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.462] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.462] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns

[2024.12.19 17:10:12.462] INFO    Schedule completed. [0.006 s]

[2024.12.19 17:10:12.491] INFO    Code generation completed for all AWGs. [0.028 s]

[2024.12.19 17:10:12.492] INFO    Completed compilation step 1 of 4. [0.035 s]

[2024.12.19 17:10:12.496] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.496] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.496] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.496] INFO      [2 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.497] INFO    Schedule completed. [0.004 s]

[2024.12.19 17:10:12.522] INFO    Code generation completed for all AWGs. [0.025 s]

[2024.12.19 17:10:12.524] INFO    Completed compilation step 2 of 4. [0.031 s]

[2024.12.19 17:10:12.528] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.528] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.528] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.528] INFO      [4 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.528] INFO    Schedule completed. [0.004 s]

[2024.12.19 17:10:12.554] INFO    Code generation completed for all AWGs. [0.025 s]

[2024.12.19 17:10:12.556] INFO    Completed compilation step 3 of 4. [0.032 s]

[2024.12.19 17:10:12.560] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.560] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.560] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.560] INFO      [6 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.561] INFO    Schedule completed. [0.004 s]

[2024.12.19 17:10:12.588] INFO    Code generation completed for all AWGs. [0.027 s]

[2024.12.19 17:10:12.589] INFO    Completed compilation step 4 of 4. [0.033 s]

[2024.12.19 17:10:12.600] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.600] INFO      Step   Device       AWG   SeqC LOC   CT entries   Waveforms   Samples  

[2024.12.19 17:10:12.601] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.601] INFO      0      hdawg_0        0         23            3           2       832  

[2024.12.19 17:10:12.602] INFO             shfqc_0        0          9            0           2     16000  

[2024.12.19 17:10:12.602] INFO             shfqc_0_sg     0         62           13          10      5696  

[2024.12.19 17:10:12.603] INFO             shfqc_0_sg     1         58           12          10      6464  

[2024.12.19 17:10:12.603] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.603] INFO      1      hdawg_0        0         23            3           2       832  

[2024.12.19 17:10:12.604] INFO             shfqc_0        0          9            0           2     16000  

[2024.12.19 17:10:12.604] INFO             shfqc_0_sg     0         48            8           6      4032  

[2024.12.19 17:10:12.604] INFO             shfqc_0_sg     1         48            9           7      4832  

[2024.12.19 17:10:12.605] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.605] INFO      2      hdawg_0        0         23            3           2       832  

[2024.12.19 17:10:12.605] INFO             shfqc_0        0          9            0           2     16000  

[2024.12.19 17:10:12.605] INFO             shfqc_0_sg     0         48            8           6      4032  

[2024.12.19 17:10:12.606] INFO             shfqc_0_sg     1         48            9           7      4832  

[2024.12.19 17:10:12.607] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.607] INFO      3      hdawg_0        0         23            3           2       832  

[2024.12.19 17:10:12.607] INFO             shfqc_0        0          9            0           2     16000  

[2024.12.19 17:10:12.608] INFO             shfqc_0_sg     0         62           13          10      5696  

[2024.12.19 17:10:12.608] INFO             shfqc_0_sg     1         59           13          10      6848  

[2024.12.19 17:10:12.608] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.609] INFO             TOTAL                   415           77                 52064  

[2024.12.19 17:10:12.609] INFO     ─────────────────────────────────────────────────────────────────────── 

[2024.12.19 17:10:12.625] INFO    Finished LabOne Q Compiler run.

We can use the pulse sheet viewer to inspect the compiled experiment. Note, that the VQE experiments are relatively long and we need to raise the number of events to see a full sequence of all VQE terms.

In [21]:

show_pulse_sheet("vqe", cexp, max_events_to_publish=3000)

show_pulse_sheet("vqe", cexp, max_events_to_publish=3000)

[2024.12.19 17:10:12.754] INFO    Recompiling the experiment due to missing extra information in the compiled experiment. Compile with `OUTPUT_EXTRAS=True` and `MAX_EVENTS_TO_PUBLISH=3000` to bypass this step with a small impact on the compilation time.

[2024.12.19 17:10:12.771] INFO    Resolved modulation type of oscillator 'q0_readout_acquire_osc' on signal '/logical_signal_groups/q0/acquire_line' to SOFTWARE

[2024.12.19 17:10:12.772] INFO    Resolved modulation type of oscillator 'q0_drive_ge_osc' on signal '/logical_signal_groups/q0/drive_line' to HARDWARE

[2024.12.19 17:10:12.773] INFO    Resolved modulation type of oscillator 'q1_readout_acquire_osc' on signal '/logical_signal_groups/q1/acquire_line' to SOFTWARE

[2024.12.19 17:10:12.773] INFO    Resolved modulation type of oscillator 'q1_drive_ge_osc' on signal '/logical_signal_groups/q1/drive_line' to HARDWARE

[2024.12.19 17:10:12.774] INFO    Starting LabOne Q Compiler run...

[2024.12.19 17:10:12.784] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.784] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.784] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns

[2024.12.19 17:10:12.796] INFO    Schedule completed. [0.017 s]

[2024.12.19 17:10:12.823] INFO    Code generation completed for all AWGs. [0.027 s]

[2024.12.19 17:10:12.824] INFO    Completed compilation step 1 of 4. [0.045 s]

[2024.12.19 17:10:12.828] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.828] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.828] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.828] INFO      [2 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.839] INFO    Schedule completed. [0.015 s]

[2024.12.19 17:10:12.864] INFO    Code generation completed for all AWGs. [0.025 s]

[2024.12.19 17:10:12.866] INFO    Completed compilation step 2 of 4. [0.042 s]

[2024.12.19 17:10:12.870] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.870] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.870] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.870] INFO      [4 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.881] INFO    Schedule completed. [0.015 s]

[2024.12.19 17:10:12.906] INFO    Code generation completed for all AWGs. [0.024 s]

[2024.12.19 17:10:12.908] INFO    Completed compilation step 3 of 4. [0.042 s]

[2024.12.19 17:10:12.913] INFO    Due to feedback latency constraints, the timing of the following match sections with corresponding handles were changed:
[2024.12.19 17:10:12.913] INFO      - 'q0_feedback' with handle 'q0_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.913] INFO      - 'q1_feedback' with handle 'q1_reset', delayed by 272.00 ns
[2024.12.19 17:10:12.913] INFO      [6 similar messages suppressed, see log file for full set of messages.]

[2024.12.19 17:10:12.924] INFO    Schedule completed. [0.015 s]

[2024.12.19 17:10:12.950] INFO    Code generation completed for all AWGs. [0.026 s]

[2024.12.19 17:10:12.952] INFO    Completed compilation step 4 of 4. [0.043 s]

[2024.12.19 17:10:12.955] INFO    Finished LabOne Q Compiler run.

[2024.12.19 17:10:12.999] INFO    Writing html file to /builds/qccs/laboneq-applications/docs/sources/how-to-guides/sources/04_qasm/vqe_2024-12-19-17-10-12.html

Out[21]:

vqe_2024-12-19-17-10-12.html

VQE objective function¶

We have now established the parts needed to compute the expectation value $\epsilon_{\theta} = \langle\psi_{\theta}|\hat{H}|\psi_{\theta}\rangle$. For this we first define the following auxiliary functions:

At present, the emulation mode does not provide the the qubit state indices. To make the acquired results compatible, the function format_emulated_results defines the measured qubit states as random integer vectors with entries $\{0, 1\}$.

In [22]:

def format_emulated_results(results, num_qubits):
    for i in range(num_qubits):
        results.acquired_results[f"q{i}"].data = np.array(
            np.random.randint(
                0, 2, size=results.acquired_results[f"q{i}"].data.shape
            )
        )
    return results

def format_emulated_results(results, num_qubits): for i in range(num_qubits): results.acquired_results[f"q{i}"].data = np.array( np.random.randint( 0, 2, size=results.acquired_results[f"q{i}"].data.shape ) ) return results

The function expectation_value_vqe computes $\langle\psi_{\theta}|\hat{P}_{i}|\psi_{\theta}\rangle$.

• The acquired results (qubit state indices) are collected into a single array with axes spanning 1. qubits, 2. single shots, and 3. VQE terms.
• The expectation value of each multi-qubit Pauli operator is either $1$ or $-1$ depending on whether an even or odd number of contributing qubits is found in state 1. Note, that this corresponds to the parity of the sum of the state indices. The parity is computed over all qubits that contribute to the expectation value, i.e. onto which a Pauli operator acts.
• The resulting parities are then averaged over the number of single shots, converted into expectation values contributions and summed up to the final $\epsilon_{\theta}$.

In [23]:

def expectation_value_vqe(acquired_results, ham):
    # sort measurement data
    state_indices = np.array(
        [acquired_results[f"q{q}"].data for q in range(ham.num_qubits)]
    ).transpose((2, 1, 0))
    # expectation value
    exp_val = 0.0
    for i, (pauli, c) in enumerate(ham.to_list()):
        exp_val += (
            c.real
            - c.real
            * np.mean(
                np.sum(state_indices[i], axis=-1, where=[p for p in pauli if p != "I"])
                % 2,
                axis=-1,
            )
            * 2.0
        )
    return exp_val

def expectation_value_vqe(acquired_results, ham): # sort measurement data state_indices = np.array( [acquired_results[f"q{q}"].data for q in range(ham.num_qubits)] ).transpose((2, 1, 0)) # expectation value exp_val = 0.0 for i, (pauli, c) in enumerate(ham.to_list()): exp_val += ( c.real - c.real * np.mean( np.sum(state_indices[i], axis=-1, where=[p for p in pauli if p != "I"]) % 2, axis=-1, ) * 2.0 ) return exp_val

Finally, we collect all the above functionality for parameter updates, experiment execution, results emulation, and VQE-postprocessing into a single objective function that returns $\langle\psi_{\theta}|\hat{P}_{i}|\psi_{\theta}\rangle$ in terms of a vector of parameter values.

In [24]:

def objective_function_vqe(ham, theta, cexp, session):
    def hamiltonian_expectation_value(x):
        if use_emulation:  # ensure consistent results at every pass
            np.random.seed(seed=10101)

        # Parameter Update
        for xx, t in zip(x, theta):
            t.val = xx
            new_pulse = pulse_from_parameter(t)
            cexp.replace_pulse(new_pulse.uid, new_pulse)

        # Execute Experiment
        res = session.run(cexp)

        # Format Results
        if use_emulation:
            res = format_emulated_results(res, ham.num_qubits)

        # Obtain expectation value
        y = expectation_value_vqe(res.acquired_results, ham)

        # Log iteration result and return
        print(f"VQE expectation value {y:20.4e}")
        return y

    return hamiltonian_expectation_value

def objective_function_vqe(ham, theta, cexp, session): def hamiltonian_expectation_value(x): if use_emulation: # ensure consistent results at every pass np.random.seed(seed=10101) # Parameter Update for xx, t in zip(x, theta): t.val = xx new_pulse = pulse_from_parameter(t) cexp.replace_pulse(new_pulse.uid, new_pulse) # Execute Experiment res = session.run(cexp) # Format Results if use_emulation: res = format_emulated_results(res, ham.num_qubits) # Obtain expectation value y = expectation_value_vqe(res.acquired_results, ham) # Log iteration result and return print(f"VQE expectation value {y:20.4e}") return y return hamiltonian_expectation_value

Optimization¶

The objective function constructed in the previous section can now directly be used for minimization. We use the "constrained optimization by linear approximation" (COBYLA) method from the scipy library for this. However, any gradient-less optimization method should be suitable for this step.

Note, that we restrict the number of evaluation steps to 1 in emulation mode as we would otherwise optimize a constant function.

In [25]:

res = minimize(
    fun=objective_function_vqe(ham, theta, cexp, session),
    x0=[t.val for t in theta],
    method="COBYLA",
    options={"maxiter": 1} if use_emulation else {},
)

res = minimize( fun=objective_function_vqe(ham, theta, cexp, session), x0=[t.val for t in theta], method="COBYLA", options={"maxiter": 1} if use_emulation else {}, )

[2024.12.19 17:10:13.039] INFO    Starting near-time execution...

[2024.12.19 17:10:13.097] INFO    Finished near-time execution.

VQE expectation value           3.5789e-03

Finally, we can retrieve $\min\limits_{\theta}\{\epsilon_{\theta}\}$ as VQE approximation of the ground state energy of $\hat{H}$

In [26]:

res.fun

res.fun

Out[26]:

0.0035789062500000968

We can also inspect the full results object returned by the minimizer for the metainformation of the optimization.

In [27]:

res

res

Out[27]:

 message: Maximum number of function evaluations has been exceeded.
 success: False
  status: 2
     fun: 0.0035789062500000968
       x: [ 3.745e-01  9.507e-01 ...  1.818e-01  1.834e-01]
    nfev: 1
   maxcv: 0.0