Skip to content

The Qmod Workshop - Part 3: Execution Flows

View on GitHub

This is the third part of the Qmod workshop, covering exercises 11 and 12. Make sure to go through Part 1 and 2 before continuing with this notebook.

Exercise 11 - Execution with Parameters

In this exercise, we will modify the state-preparation model from exercise 10 to accept a rotation angle as an execution parameter (see solution to exercise 10 at the end of QMOD_Workshop_Part_2.ipynb). See more under Quantum Entry Point.

  1. Add a classical real-type parameter rotation_angle to function main.
  2. Pass rotation_angle as a parameter to the controlled RY instead of using pi/3 directly.
  3. Sample the outputs of the quantum program using ExecutionSession.sample() for rotation angle values pi/3 and make sure the results are consistent with the ones in exercise 10, with the hard-coded rotation angle value.
  4. Do the same in a loop for all values in the \(i\frac{\pi}{4}\) for i values 0..7
from classiq import *
from classiq.qmod.symbolic import pi


@qfunc
def main(rotation_angle: CReal, res: Output[QArray[QBit]]) -> None:
    allocate(1, res)
    # Your model code here


model = create_model(main)
qprog = synthesize(model)

# Your execution code here

Exercise 12 - VQE

The Variational Quantum Eigensolver is an algorithm that finds the minimal eigenvalue of a matrix by executing a parametric circuit (also referred to as an ansatz), estimating the expected value of the matrix for the state the circuit creates (from the distribution received by the execution), and using a classical optimizer to select the next set of parameters for the circuit, until reaching convergence (or exceeding a set amount of maximum iterations).

The estimation of the expectation value is done on Pauli based matrices, so any matrix we want to perform this operation on, need to be decomposed into a sum of Pauli terms.

In this exercise, we will create a simple VQE algorithm that estimates the minimal eigenvalue of a 2x2 matrix. Fill in the gaps in the following snippet to find the minimal eigenvalue and it corresponding eigenstate for

[[1, -1], [-1, 0]] = 1/2*I + 1/2*Z - X

Hint for HamiltonianHAMILTONIAN = QConstant(...)
Hint for cmain res = vqe(
 hamiltonian=...,
 maximize=...,
 initial_point=[],
 optimizer=Optimizer.COBYLA,
 max_iteration=1000,
 tolerance=0.001,
 step_size=0,
 skip_compute_variance=False,
 alpha_cvar=1.0,
)
save({"result": res})
from typing import List

from classiq import *

# Your code here:


@qfunc
def main(q: Output[QBit], angles: CArray[CReal, 3]) -> None:
    allocate(1, q)
    U(angles[0], angles[1], angles[2], 0, q)


@cfunc
def cmain() -> None:
    # Your code here
    pass


qmod = create_model(main, classical_execution_function=cmain)
qprog = synthesize(qmod)
show(qprog)
res = execute(qprog)

# TODO: Upon completion uncomment the code below
# vqe_result = res.result()[0].value
# print(vqe_result.energy, vqe_result.optimal_parameters, vqe_result.eigenstate)

Note: The U gate is a general rotation matrix on a single qubit, so the given model creates an ansatz that spans all of the space for a single qubit, and thus gives us a full search space for this specific problem.

Bonus: Exercise 13 - Quantum Counting with Iterative Quantum Amplitude Estimation

Quantum Counting algorithm is an algorithm for efficiently estimating the number of valid solutions to a search problem, based on the amplitude estimation algorithm. It demonstrates a quadratic improvement in regard to a classical algorithm with black-box oracle access to the function \(f\).

More precisely, the counting problem is, given a boolean function \(f :\{0, 1\}^n\rightarrow\{0,1\}\), estimate the number of inputs \(x\) to \(f\) such that \(f(x)=1\).

Let's demonstrate how to estimate the Counting problem using a specific variant of the Amplitude Estimation algorithm - the Iterative Quantum Amplitude Estimation (IQAE).

The IQAE does not rely on the Quantum Phase Estimation algorithm, but purely on applications of the grover operator:

\[Q\equiv - A S_0 A^{\dagger} S_{\psi_1},\]

Hence reducing the required amount of qubits and gates of the circuit, at the expense of additional multiplicative factor poly-logarithmic in the error \(\epsilon\).

For that we need a state preparation with an indicator qubit point the valid solution among the states:

from classiq import QArray, QBit, QNum, bind, hadamard_transform, qfunc

A_SIZE = 2
B_SIZE = 2
DOMAIN_SIZE = A_SIZE + B_SIZE


@qfunc
def arith_equation(a: QNum, b: QNum, res: QBit):
    res ^= a + b <= 2


@qfunc
def iqae_state_preparation(a: QNum, b: QNum, res: QBit):
    reg = QArray("reg")
    bind([a, b, res], reg)
    hadamard_transform(reg[0:DOMAIN_SIZE])
    bind(reg, [a, b, res])
    arith_equation(a, b, res)

Then, the quantum circuit that is needed to the iterative QAE scheme, which needs to apply powers of the grover operator:

from classiq import CInt, QCallable, grover_operator, power


@qfunc
def my_iqae_algorithm(
    k: CInt,
    oracle_operand: QCallable[QArray[QBit]],
    sp_operand: QCallable[QArray[QBit]],
    x: QArray[QBit],
):
    sp_operand(x)
    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))

Now, after we already know how to attach execution parameters, let's build a parametric main, and a cmain assigning values to the parameters, in order to actually run hybrid execution algorithm. So the parametric main should get as an input a parameter named k:

from classiq import Output, Z, allocate


@qfunc
def main(
    k: CInt,
    ind_reg: Output[QBit],
) -> None:
    full_reg = QArray("full_reg")
    allocate(DOMAIN_SIZE + 1, full_reg)
    my_iqae_algorithm(
        k=k,
        oracle_operand=lambda x: Z(x[x.len - 1]),
        sp_operand=lambda x: iqae_state_preparation(
            x[0:A_SIZE], x[A_SIZE : x.len - 1], x[x.len - 1]
        ),
        x=full_reg,
    )
    state_reg = QArray("state_reg", length=DOMAIN_SIZE)
    bind(full_reg, [state_reg, ind_reg])

And the cmain will attach the relevant k, behind the scenes, all you need is: * Call iqae in cmain, with the required epsilon=1 / ((2**DOMAIN_SIZE_QCONST) * 2) and alpha=0.01. * Save its result into a variable named iqae_res using save.

Fill in the cmain:

from classiq import QConstant, bind, cfunc, iqae, save

DOMAIN_SIZE_QCONST = QConstant("DOMAIN_SIZE_QCONST", int, DOMAIN_SIZE)


@cfunc
def cmain():
    # your code
    pass

Now let's run it and examine how many results apply the condition a+b<=2 (there should be 6 of them: (a,b)=(0,0),(0,1),(1,0),(1,1),(2,0),(0,2)).

from classiq import Constraints, create_model, show, synthesize

constraints = Constraints(optimization_parameter="width")
qmod_iqae = create_model(
    entry_point=main,
    constraints=constraints,
    classical_execution_function=cmain,
)

qprog_iqae = synthesize(qmod_iqae)
show(qprog_iqae)

Uncomment this cell and run it after completing the cmain function.

# from classiq import execute

# res = execute(qprog_iqae).result()

# iqae_res = res[0].value
# print(
#     f"IQAE result: {iqae_res.estimation}, confidence interval: {iqae_res.confidence_interval}"
# )

# print(
#     f"Number of solutions: {(2**DOMAIN_SIZE) * iqae_res.estimation}, accuracy: "
#     f"{(2**DOMAIN_SIZE)*(iqae_res.confidence_interval[1]-iqae_res.confidence_interval[0])}"
# )

Solutions

Exercise 11

# Solution to exercise 11:


from math import pi

from classiq import *
from classiq.execution import ExecutionSession


@qfunc
def main(rotation_angle: CReal, res: Output[QArray[QBit]]) -> None:
    x: QArray[QBit] = QArray("x")
    allocate(3, x)
    hadamard_transform(x)

    ls_bit = QBit("ls_bit")
    ms_bits = QNum("ms_bits", 2, False, 0)
    bind(x, [ls_bit, ms_bits])

    control(ms_bits == 1, lambda: RY(rotation_angle, ls_bit))

    bind([ls_bit, ms_bits], res)


model = create_model(main)
qprog = synthesize(model)
show(qprog)

es = ExecutionSession(qprog)

ed = es.sample({"rotation_angle": pi / 3})
print(ed.parsed_counts)

for i in range(8):
    ed = es.sample({"rotation_angle": i * pi / 4})
    print(ed.parsed_counts)

Exercise 12

# Solution to exercise 12:


from typing import List

from classiq import *

HAMILTONIAN = QConstant(
    "HAMILTONIAN",
    List[PauliTerm],
    [PauliTerm([Pauli.I], 0.5), PauliTerm([Pauli.Z], 0.5), PauliTerm([Pauli.X], -1.0)],
)


@qfunc
def main(q: Output[QBit], angles: CArray[CReal, 3]) -> None:
    allocate(1, q)
    U(angles[0], angles[1], angles[2], 0, q)


@cfunc
def cmain() -> None:
    res = vqe(
        HAMILTONIAN,
        False,
        [],
        optimizer=Optimizer.COBYLA,
        max_iteration=1000,
        tolerance=0.001,
        step_size=0,
        skip_compute_variance=False,
        alpha_cvar=1.0,
    )
    save({"result": res})


qmod = create_model(main, classical_execution_function=cmain)
qprog = synthesize(qmod)
show(qprog)
res = execute(qprog)
vqe_result = res.result()[0].value
print(vqe_result.energy, vqe_result.optimal_parameters, vqe_result.eigenstate)

Exercise 13

from classiq import QArray, QBit, QCallable, QNum, bind, hadamard_transform, qfunc

A_SIZE = 2
B_SIZE = 2
DOMAIN_SIZE = A_SIZE + B_SIZE


@qfunc
def arith_equation(a: QNum, b: QNum, res: QBit):
    res ^= a + b <= 2


@qfunc
def iqae_state_preparation(a: QNum, b: QNum, res: QBit):
    reg = QArray("reg")
    bind([a, b, res], reg)
    hadamard_transform(reg[0:DOMAIN_SIZE])
    bind(reg, [a, b, res])
    arith_equation(a, b, res)


from classiq import CInt, grover_operator, power


@qfunc
def my_iqae_algorithm(
    k: CInt,
    oracle_operand: QCallable[QArray[QBit]],
    sp_operand: QCallable[QArray[QBit]],
    x: QArray[QBit],
):
    sp_operand(x)
    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))


from classiq import Output, Z, allocate


@qfunc
def main(
    k: CInt,
    ind_reg: Output[QBit],
) -> None:
    full_reg = QArray("full_reg")
    allocate(DOMAIN_SIZE + 1, full_reg)
    my_iqae_algorithm(
        k=k,
        oracle_operand=lambda x: Z(x[x.len - 1]),
        sp_operand=lambda x: iqae_state_preparation(
            x[0:A_SIZE], x[A_SIZE : x.len - 1], x[x.len - 1]
        ),
        x=full_reg,
    )
    state_reg = QArray("state_reg", length=DOMAIN_SIZE)
    bind(full_reg, [state_reg, ind_reg])


from classiq import QConstant, bind, cfunc, iqae, save

DOMAIN_SIZE_QCONST = QConstant("DOMAIN_SIZE_QCONST", int, DOMAIN_SIZE)


@cfunc
def cmain():
    iqae_res = iqae(epsilon=1 / ((2**DOMAIN_SIZE_QCONST) * 2), alpha=0.01)
    save({"iqae_res": iqae_res})


from classiq import Constraints, create_model, show, synthesize

constraints = Constraints(optimization_parameter="width")
qmod_iqae = create_model(
    main,
    constraints=constraints,
    classical_execution_function=cmain,
)

qprog_iqae = synthesize(qmod_iqae)
show(qprog_iqae)

from classiq import execute

res = execute(qprog_iqae).result()

iqae_res = res[0].value
print(
    f"IQAE result: {iqae_res.estimation}, confidence interval: {iqae_res.confidence_interval}"
)
print(
    f"Number of solutions: {(2**DOMAIN_SIZE) * iqae_res.estimation}, accuracy: "
    f"{(2**DOMAIN_SIZE)*(iqae_res.confidence_interval[1]-iqae_res.confidence_interval[0])}"
)
for i, iteration in enumerate(iqae_res.iterations_data):
    print(
        f"iteration_id: {i}, num grover iterations: {iteration.grover_iterations}, counts: {iteration.sample_results.counts}"
    )

Exercise 13

from classiq import QArray, QBit, QCallable, QNum, bind, hadamard_transform, qfunc

A_SIZE = 2
B_SIZE = 2
DOMAIN_SIZE = A_SIZE + B_SIZE


@qfunc
def arith_equation(a: QNum, b: QNum, res: QBit):
    res ^= a + b <= 2


@qfunc
def iqae_state_preparation(a: QNum, b: QNum, res: QBit):
    reg = QArray("reg")
    bind([a, b, res], reg)
    hadamard_transform(reg[0:DOMAIN_SIZE])
    bind(reg, [a, b, res])
    arith_equation(a, b, res)


from classiq import CInt, grover_operator, power


@qfunc
def my_iqae_algorithm(
    k: CInt,
    oracle_operand: QCallable[QArray[QBit]],
    sp_operand: QCallable[QArray[QBit]],
    x: QArray[QBit],
):
    sp_operand(x)
    power(k, lambda: grover_operator(oracle_operand, sp_operand, x))


from classiq import Output, Z, allocate


@qfunc
def main(
    k: CInt,
    ind_reg: Output[QBit],
) -> None:
    full_reg = QArray("full_reg")
    allocate(DOMAIN_SIZE + 1, full_reg)
    my_iqae_algorithm(
        k=k,
        oracle_operand=lambda x: Z(x[x.len - 1]),
        sp_operand=lambda x: iqae_state_preparation(
            x[0:A_SIZE], x[A_SIZE : x.len - 1], x[x.len - 1]
        ),
        x=full_reg,
    )
    state_reg = QArray("state_reg", length=DOMAIN_SIZE)
    bind(full_reg, [state_reg, ind_reg])


from classiq import QConstant, bind, cfunc, iqae, save

DOMAIN_SIZE_QCONST = QConstant("DOMAIN_SIZE_QCONST", int, DOMAIN_SIZE)


@cfunc
def cmain():
    iqae_res = iqae(epsilon=1 / ((2**DOMAIN_SIZE_QCONST) * 2), alpha=0.01)
    save({"iqae_res": iqae_res})


from classiq import Constraints, create_model, show, synthesize

constraints = Constraints(optimization_parameter="width")
qmod_iqae = create_model(
    main,
    constraints=constraints,
    classical_execution_function=cmain,
)

qprog_iqae = synthesize(qmod_iqae)
show(qprog_iqae)

from classiq import execute

res = execute(qprog_iqae).result()

iqae_res = res[0].value
print(
    f"IQAE result: {iqae_res.estimation}, confidence interval: {iqae_res.confidence_interval}"
)
print(
    f"Number of solutions: {(2**DOMAIN_SIZE) * iqae_res.estimation}, accuracy: "
    f"{(2**DOMAIN_SIZE)*(iqae_res.confidence_interval[1]-iqae_res.confidence_interval[0])}"
)
for i, iteration in enumerate(iqae_res.iterations_data):
    print(
        f"iteration_id: {i}, num grover iterations: {iteration.grover_iterations}, counts: {iteration.sample_results.counts}"
    )