Using QSVT for fixed-point amplitude amplification
This demo will show how to use the QSVT framework for search problems. Specifically, we will implement fixed-point amplitude amplification (FPAA). In FPAA, we do not know in advance the concentration of solutions to the search problem, want to sample a solution with high probability. In contrast, for the original grover search algorithm, too much iterations might 'overshoot'.
The demo is based on the paper Grand unification of quantum algorithms.
Given \(|s\rangle\) the initial state and \(|t\rangle\) the 'good' states, we get an effective block encoding of a 1-dimensional matrix \(A=|t\rangle\langle s|\).
Given that \(a = \langle s|t\rangle\gt0\), we want to amplify \(a\). The signal operator \(U\) here will be \(I\) (and also \(\dagger{U}\)). Now we implement 2 projector-rotations - one in '\(|s\rangle\)' space and one in '\(|t\rangle\)' space, each one around the given state, giving phase to the specific state.
Defining the QSVT circuit for the problem
We start with the general qsvt framework definition. It accepts a unitary which block-encode a matrix together with projector-controlled-phase functions which rotate the state around each of the subspaces where the matrix is encoded.
It applies the qsvt_step
multiple times, iterating over the rotation angles provided which encode the polynomial transformation.
Notice - The last step is quite tricky and depend on the specific transformation we wish to perform. Here the code is suitable for the FP Amplitude Amplification case. Also - we wrap the auxilliary qubit with \(H\) gates.
from classiq import *
from classiq.qmod.symbolic import floor, logical_and
@qfunc
def my_qsvt_step(
phase1: CReal,
phase2: CReal,
proj_ctrl_phase_1: QCallable[CReal, QArray[QBit], QBit],
proj_ctrl_phase_2: QCallable[CReal, QArray[QBit], QBit],
u: QCallable[QArray[QBit]],
state: QArray[QBit],
aux: QBit,
):
u(state)
proj_ctrl_phase_2(phase1, state, aux)
invert(lambda: u(state))
proj_ctrl_phase_1(phase2, state, aux)
@qfunc
def my_qsvt(
phase_seq: CArray[CReal],
proj_ctrl_phase_1: QCallable[CReal, QArray[QBit], QBit],
proj_ctrl_phase_2: QCallable[CReal, QArray[QBit], QBit],
u: QCallable[QArray[QBit]],
state: QArray[QBit],
aux: QBit,
) -> None:
H(aux)
proj_ctrl_phase_1(phase_seq[0], state, aux)
repeat(
count=floor((phase_seq.len - 1) / 2),
iteration=lambda index: my_qsvt_step(
phase_seq[2 * index + 1],
phase_seq[2 * index + 2],
proj_ctrl_phase_1,
proj_ctrl_phase_2,
u,
state,
aux,
),
)
if_(
condition=phase_seq.len % 2 == 1,
then=lambda: IDENTITY(state),
else_=lambda: (
u(state),
proj_ctrl_phase_2(
phase_seq[phase_seq.len - 1],
state,
aux,
),
),
)
H(aux)
Initial state z-rotation - rotation around \(|s\rangle\)
def initial_state_rot(phase: CReal, state: QArray[QBit], aux: QBit):
hadamard_transform(state)
apply_to_all(X, state)
control(ctrl=state, stmt_block=lambda: X(aux))
RZ(phase, aux)
control(ctrl=state, stmt_block=lambda: X(aux))
apply_to_all(X, state)
hadamard_transform(state)
Good states z-rotation - rotation around \(|t\rangle\)
def target_state_rot(
phase: CReal,
arith_oracle: QCallable[QArray[QBit], QBit],
state: QArray[QBit],
aux: QBit,
):
arith_oracle(state, aux)
RZ(phase, aux)
arith_oracle(state, aux)
Defining the arithmetic Oracle
Implement the following equation:
(a + b) == 3 and (c - a) == 2
with a, b, c
in sizes 2, 1, 3.
class OracleVars(QStruct):
a: QNum[2, False, 0]
b: QNum[1, False, 0]
c: QNum[3, False, 0]
REG_SIZE = 6
@qfunc
def arith_equation(state: OracleVars, res: QBit):
res ^= logical_and((state.a + state.b) == 3, (state.c - state.a) == 2)
Wrapping everything together for the FPAA case
In the FPAA case, the provided unitary is just the Identity matrix! In addition, we provide both projector-controlled phase functions - the initial and target state rotations.
@qfunc
def qsvt_fpaa(
phase_seq: CArray[CReal],
arith_oracle: QCallable[QArray[QBit], QBit],
state: QArray[QBit],
aux: Output[QBit],
) -> None:
allocate(1, aux)
my_qsvt(
phase_seq,
lambda phase, state, aux: initial_state_rot(phase, state, aux),
lambda phase, state, aux: target_state_rot(phase, arith_oracle, state, aux),
lambda state: IDENTITY(target=state),
state,
aux,
)
get the phase sequence for the sign function
Now we will use the package pyqsp
in order to get the phases of the rotation sequence.
Get directly the coef of the sign function, based on the erfc approximation.
# The following code assumes pyqsp version 0.1.6
# !pip install pyqsp==0.1.6
import pyqsp
DEGREE = 25
X_BASIS = True
pg = pyqsp.poly.PolySign()
pcoefs, scale = pg.generate(
degree=DEGREE, delta=10, ensure_bounded=True, return_scale=True
)
[pyqsp.poly.PolySign] degree=25, delta=10
[PolyTaylorSeries] max [0.89999962] is at [1.]: normalizing
[PolyTaylorSeries] average error = 0.09627807622454267 in the domain [-1, 1] using degree 25
import numpy as np
from pyqsp.angle_sequence import Polynomial, QuantumSignalProcessingPhases
poly = Polynomial(pcoefs)
measurement = "x"
ang_seq = QuantumSignalProcessingPhases(
poly, signal_operator="Wx", method="laurent", measurement=measurement
)
pyqsp.response.PlotQSPResponse(ang_seq, signal_operator="Wx", measurement=measurement)
As \(R(a)=-i*e^{i\frac{\pi}{4}Z}W(a)e^{i\frac{\pi}{4}Z}\) and we have odd number of rotations, we get an \(i\) phase to our polynomial, so we get \(Im(P(a))\) instead of the real part. So we will get the result in the \(|1\rangle\) state in the ancilla. However, we can fix it by adding \(\pi/2\) phase to the last or first rotation.
#### change the R(x) to W(x), as the phases are in the W(x) conventions
phases = np.array(ang_seq)
phases[1:-1] = phases[1:-1] - np.pi / 2
phases[0] = phases[0] - np.pi / 4
phases[-1] = phases[-1] + (2 * (len(phases) - 1) - 1) * np.pi / 4
phases = (
-2 * phases
) # verify conventions. minus is due to exp(-i*phi*z) in qsvt in comparison to qsp
Create the full QSVT Model
@qfunc
def main(state: Output[OracleVars], aux: Output[QBit]):
allocate(state.size, state)
hadamard_transform(state)
qsvt_fpaa(
phase_seq=list(phases),
arith_oracle=arith_equation,
state=state,
aux=aux,
)
Synthesis and execution on a simulator
We will use Classiq's synthesis engine to translate the model to a quantum circuit, and execute on Classiq's simulator
from classiq.execution import (
ClassiqBackendPreferences,
ClassiqSimulatorBackendNames,
ExecutionPreferences,
)
constraints = Constraints(max_width=12)
NUM_SHOTS = 1000
# we will want to execute this qmod on a state-vector simulator:
execution_preferences = ExecutionPreferences(
num_shots=NUM_SHOTS,
backend_preferences=ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR
),
)
qmod = create_model(
main,
constraints,
execution_preferences,
out_file="qsvt_fixed_point_amplitude_amplification",
)
qprog = synthesize(qmod)
show(qprog)
Execute the circuit:
result = execute(qprog).result_value()
def equation(a, b, c):
return ((a + b) == 3) and ((c - a) == 2)
measured_good_shots = 0
for r in result.parsed_counts:
a, b, c, aux = (
r.state["state"]["a"],
r.state["state"]["b"],
r.state["state"]["c"],
r.state["aux"],
)
if equation(a, b, c) and (aux == 0):
print(
f"a: {a}, b: {b}, c: {c}, aux: {aux}, equation_result: {equation(a, b, c)}, counts={r.shots}"
)
measured_good_shots += r.shots
print("Measured good shots:", measured_good_shots)
a: 2, b: 1, c: 4, aux: 0.0, equation_result: True, counts=419
a: 3, b: 0, c: 5, aux: 0.0, equation_result: True, counts=394
Measured good shots: 813
What do we expect?
We need to subtitue the amplitude of \(|s\rangle\langle t|\) in \(P(x)\):
p_good_shot = poly(np.sqrt(2 / 2**6)) ** 2
print("Expected good shots:", NUM_SHOTS * p_good_shot)
Expected good shots: 832.9890148378344
Indeed, we got the expected result according to the polynomial we created with the QSVT sequence.
import scipy
assert np.isclose(
measured_good_shots,
NUM_SHOTS * p_good_shot,
atol=5 * scipy.stats.binom.std(NUM_SHOTS, p_good_shot),
)
References
[1]: Martyn JM, Rossi ZM, Tan AK, Chuang IL. Grand unification of quantum algorithms. PRX Quantum. 2021 Dec 3;2(4):040203.