Dyson Series Expansion¶
1. Dynamic Pictures Overview¶
In this tutorial, we consider the Hamiltonian of the form $H(t)= H_0 + V(t)$, where the operator $V(t)$ denotes a transition (often called "perturbation") between eigenstates of the $H_0$ operator.
The Quantum Mechanics theory uses three approaches in describing an evolution of the Hamiltonian $H(t)$ over time $t$:
- The Schrodinger Picture: the evolution is reflected as a change in the state vector $\psi(t)$ while the time-independent operator $H$ is fixed.
- The Heisenberg Picture: the evolution of the system is observed into a quantum operator. In contrast to the Schrodinger picture, in the Heisenberg picture the state $\psi$ does not change with time, rather the quantum operator $H(t)$ changes.
- The Interaction Picture: the evolution of a quantum system is represented as a hybrid approach of the Shrodinger and Heisenberg pictures. In other words, some portion of the quantum dynamics is carried by the quantum state $\psi(t)$ and some - by the quantum operator $H(t)$.
The Table 1 gives the summaries of the differences among the three pictures.
Quantum Mechanics Pictures Comparison¶
from IPython.display import Image
Image(filename='../images/qp_comparison.png', width=800, height=600)
Table 1. The evolution comparison in Quantum Mechanics pictures. Source: Wikipedia: Interaction picture ⧉.
In this tutorial, we will explore the Hamiltonian simulation in the Interaction picture using the Dyson Series Qubrick available in WBA.
2. Hamiltonian Simulation in the Interaction Picture¶
We start with the Schrödinger's Picture and slowly transition to the Interaction Picture by introducing appropriate "picture-transition operators". One can represent the Hamiltonian as the sum of an unperturbed term $H_0$ and a perturbation term $V(t)$:
$$ H_{S}(t) = H_{0, S} + V_{S}(t). \tag{1} $$
Note that $ H_{0, S}$ has no explicit time dependence. At the start time $t=0$, this operator coincides with the one in the Interaction picture $H_{0, S} =H_{0, I}$. Similarly, at $t=0$ the dynamic operator in the Heisenberg picture coincides with the operator in the Schrödinger picture $V_H(0)= V_S(0)$. For simplicity, we drop the $H$, $S$ and $I$ subscripts for operators as soon as it becomes apparent on what picture we are referring to.
We use the conversion operator $U_{conv}(t) = \exp(-i H_{0,S}t)$ to transition between a quantum state in the Schrödinger Picture and the Interaction Picture:
$$ |ψ_I(t)⟩ = U_{conv}^{\dagger}(t)|ψ_S(t)⟩ = \exp(-i H_{0,S}t)|ψ_S(t)⟩. \tag{2} $$
We then can construct a Schrödinger-like equation of motion in the Interaction Picture:
$$ i \frac{d}{dt}|ψ_I(t)⟩ = U_{conv}^{\dagger}(t)H_{S}(t) U_{conv}(t)|ψ_I(t)⟩. \tag{3} $$
Now, we are interested to find such an operator $U_{D}(t)$ that describes how a quantum state $|ψ_I(t)⟩$ evolves from the time $t=0$ to some time $t$. More generally, we are interested in solving the equation of the form
$$ i\frac{d}{dt}U_{D}(t) = W(t) U_{D}(t) \tag{4} $$
where $W(t) = U_{conv}^{\dagger}(t)V U_{conv}(t) = \exp(-i H_{0}t)V \exp(i H_{0}t)$. The solution to this equation is given by the Dyson Operator.
In general, we can split the Hamiltonian $H(t)$ as we pleased. However, to accommodate the easiness of mapping the Hamiltonian to a quantum circuit, we split $H(t)$ into a diagonal operator (fast-forwardable operator) and the off-diagonal operator which can be represented through a block-encoding routine.
3. Dyson Series¶
Formally, the solution to the equation $(4)$ (we solve for the $U_{D}$ operator) we introduced in the last section is given by the Dyson Series expansion:
$$ U_{D}(t) = 1 + \sum_{n=1}^{\infty}(-i)^n \int_{0 < t_1}\int_{t_1 < t_2} \ldots \int_{t_n < t} dt_n \ldots dt_1 W(t_n) \ldots W(t_1). \tag{5} $$
Note the time ordering inside the integrals: the rightmost operator $W(t_1)$ acts first, followed by the second rightmost operator $W(t_2)$ and so on. In literature such ordering is usually referred to as the time-ordering operator denoted by $\mathcal{T}$ and is used to compactly represent the Dyson Series Expansion:
$$ U_{D}(t)= \mathcal{T} \left ( \exp (-i \int_{0}^{t}W(t)dt) \right). \tag{6} $$
But before we look at an example, there are two crucial parameters we need to introduce: the truncation parameter $K$ and the time-discretization parameter $M$.
The truncation parameter $K$ arises due to the fact that computationally we replace an infinite sum with its Taylor-truncated version by ignoring the tems of a higher order. The discretization parameter $M$ appears due to the replacement of each of the integrals in $U_{D}$ by finite sums with $M$ terms. In practice, instead of the exact solution, we solve for the $(K, M)$-truncated-(time)-discretized version of the Dyson Series expansion:
$$ D_{(K, M)}(t) = I + \sum_{k=1}^{K}{\frac{(-it)^k}{M^k k!}} \left( \sum^{M-1}_{m_k=k-1} \ldots \sum^{m_3-1}_{m_2=1} \sum^{m_2-1}_{m_1=0} W \left(\frac{m_kt}{M} \right) \ldots W\left(\frac{m_2 t}{M} \right) W \left(\frac{m_1t}{M} \right) \right) \tag{7} $$
where $I$ is the identity operator. In the Workbench Algorithms we have DysonSeriesLCU qubrick implemented that can handle for us the summation and integration.
Note that the $D_{(K, M)}(t)$ approximates the Dyson Operator $U_D(t)$, and the larger the truncation $K$ and the discretization $M$ parameters are, the better the approximation is.
4. Dyson Series with Unary Truncation Register in WBA¶
%load_ext autoreload
%autoreload 2
import numpy as np
from psiqworkbench.qubricks import NaiveAdd
from workbench_algorithms.utils import PauliSum, PauliMask, pauli_sum_to_numpy
from workbench_algorithms.utils.antisymmetrization_utils import total_sort_num_comps
from workbench_algorithms.subroutines.dyson import (
DysonSeriesPrepareUnary,
DysonSeriesPrepareBinary,
DysonSeriesSelectUnary,
DysonSeriesSelectBinary,
DysonSeriesLCU,
FlagCollisionsUnary,
FlagCollisionsBinary
)
from workbench_algorithms.utils.dyson_utils import (
DysonSeriesUnaryIndex,
DysonSeriesBinaryIndex,
DysonSeriesData,
DysonSelectStrategyNaive,
dyson_target_prob,
num_qubits_DysonSeriesExpansionIndex,
norm_beta,
dyson_series_numpy,
)
from psiqworkbench import QPU, Qubits, Qubrick, Units
from psiqworkbench.utils.numpy_utils import reverse_numpy_op
from workbench_algorithms import PrepareClockState
from workbench_algorithms.subroutines.naive_state_prep import PrepareNaive
from workbench_algorithms.subroutines.select import SelectNaive
from workbench_algorithms.subroutines.block_encoding import LCU
import random
Although the expression above may look intimidating at first glance, it is easily implementable as $Prep \cdot Sel \cdot Prep ^{\dagger}$ routine, or you can also think about it in terms of the LCU routine, except that we need to add an extra logic to keep track of the truncation and time-discretization parameters. To do so, we allocate two quantum registers: one of the size $K$ which represents the outer summation, and the second register of the size $K\log_2(M)$ to represent the time-discretization at each of the $k$-th term of the outer sum. We have compactly incorporated both registers as a CompositeRegister WB's data structure.
On the truncation register, we encode the coefficients of the Taylor series expansion $p_k := \frac{t^k}{\beta k!}$ where $\beta$ is a normalization factor defined as $\beta = \sum_{k=0}^{K} \frac{t^k}{k!}$: $$ PREP |0⟩^{\otimes n} = \sum_{k=0}^n \sqrt{p_k}|1⟩^{\otimes k}|0⟩^{\otimes n - k} \tag{8} $$ which you can think of as a unary encoding of numbers with amplitudes such that $\sum_{k} p_k =1$. For example, the probability distribution over numbers $0, 1, 2, 3$ is
$$ \sqrt{p_0}|000⟩ + \sqrt{p_1}|100⟩ + \sqrt{p_2}|110⟩ + \sqrt{p_3}|111⟩ $$
where $|.⟩$ is an integer in unary representation, for instance, $111_{1} = 3_{10}, 110_{1} = 2_{10}$.
The dyson_target_fun function calculates the probabilities. The PrepareClockState Qubrick is a subroutine that encode the coefficients of the expansion to a quantum register.
We then entangle the truncation register with the discretization register using the CNOT gates conditioned to the truncation register such that for eack $k < K$ we multiply $k$ many $W \left(t \right)$ operators. We then sort the discretization register in increasing order of the $m_k$ parameter utilizing BitonicSort().
We represent the off-diagonal operator $V$ of the Hamiltonian we want to simulate as a block-encoding routine $BE(V)$ that acts on the system and an ancillary register $|0⟩_{\text{anc}}$
$$ HAM(t) = U^{\dagger}(t) \cdot V \cdot U(t) \approx ⟨0|_{\text{anc}} U_{conv}^{\dagger}(t) \cdot BE(V) \cdot U_{conv}(t) |0⟩_{\text{anc}} = ⟨0|_{\text{anc}} \exp(i H_{0}t)\cdot BE(V) \cdot \exp(-i H_{0}t) |0⟩_{\text{anc}}. \tag{9} $$
When we project the ancillary register onto the zero subspace, we recover the action of the Hamiltonian $HAM(t)$ .
Note. We also use a dataclass DysonSeriesData to keep all the required parameters such as truncation, discretization, time, time-independent Hamiltonian, and boolean collision for a Dyson Series in one place.
Hamiltonian Simulation Example¶
Now, it is time to define the Hamiltonian of the form $H(t)= H_0 + V(t)$ and simulate it quantumly!
For the scope of this tutorial, we use a toy example of the LCU of OFF_DIAG_OP = X + Y as the block-encoding part (off-diagonal operator) of the Hamiltonian which is wrapped up in the PauliSum and the FastForwardableOperator representing the time-dependent portion.
Note that the fast-forwarding subroutine, FastForwardableOperator, is represented in our example as a phase-kickback routine but does not have to be defined as a phase kickback in general, and constructs the time-dependent unitary $U_{conv}(t)=\exp{(-iH_0t)}$.
# off-diagonal operator
OFF_DIAG_OP = PauliSum(
[1, PauliMask(1, 0)],
[1, PauliMask(1, 1)]
)
# fast-forwarable operator
class FastForwardableOperator(Qubrick):
r"""Time-dependent Hamiltonian (fast-forwardable) ```H_0```.
Define: ```H_0 = exp(- i * t * theta (I - Z) / 2)```.
```\ket{time}\ket{psi} -> exp(- i * t * theta (1 - (-1)^(psi)) / 2)\ket{time}\ket{psi} ```.
Note that time is expected to be represented as an integer of ```b``` bits,
so that ```t = time / 2^b```.
"""
def __init__(self, norm, theta=1.0, **kwargs):
super().__init__(name="[PauliHam]", **kwargs)
self.norm = norm # discritization param in practice
self.theta = theta
self.ham = PauliSum(
[-0.5, PauliMask(0, 1)],
[0.5, PauliMask(0, 0)]
)
def _compute(self, psi: Qubits, time: Qubits, ctrl=0):
"""Apply Phase Kickback.
Args:
psi (Qubits): a system register.
ctrl (Qubits): a control register.
"""
time_size = time.num_qubits
# norm_factor = 2 ** time_size
for indx in range(time_size):
theta = -1 * self.theta * (2 ** indx) / self.norm
time[indx].phase(theta=theta * Units.radian, cond=psi | ctrl)
Now we set up all the parameters as $t= 0.25, K=2, M=3$. Note that this is some approximation of the target unitary, which gets closer to the target as $(K, M)$ gets larger.
time, trunc_bound, disc_bound = 0.25, 2, 3
num_system = max(OFF_DIAG_OP.width(), 1)
be_size = int(np.ceil(np.log2(len(OFF_DIAG_OP))))
num_comps = total_sort_num_comps(trunc_bound)
num_dyson_idx_reg = num_qubits_DysonSeriesExpansionIndex(
trunc_bound, disc_bound, num_system, be_size, reg_type="unary"
)
num_anc = int(np.ceil(np.log2(disc_bound)))
total_num_qubits = num_dyson_idx_reg + num_comps + num_anc
onenorm = OFF_DIAG_OP.norm()
be_matrix = pauli_sum_to_numpy(OFF_DIAG_OP)
diag_op = FastForwardableOperator(disc_bound, time / onenorm)
diag_matrix = pauli_sum_to_numpy(diag_op.ham)
# construct the Dyson Series circuit
qc = QPU()
qc.reset(total_num_qubits)
psi = Qubits(num_system, "psi_reg", qc)
dyson_reg = DysonSeriesUnaryIndex.initialize(
qc, "DysonReg", trunc_bound, disc_bound, psi, be_size
)
flag_t_collisions = False
data = DysonSeriesData(disc_bound, trunc_bound, flag_t_collisions, OFF_DIAG_OP, time)
trunc_reg_probs = dyson_target_prob(time, trunc_bound)
flag_qbk = FlagCollisionsUnary()
prepare = PrepareClockState()
dyson_prep = DysonSeriesPrepareUnary(prepare)
lcu_prepare = PrepareNaive(OFF_DIAG_OP.get_coefficients())
lcu_select = SelectNaive()
lcu = LCU(lcu_prepare, lcu_select)
dyson_select = DysonSeriesSelectUnary(diag_op, lcu, flag_qbk)
dyson_op = DysonSeriesLCU(dyson_prep, dyson_select)
dyson_op.compute(dyson_reg, ctrl=0, data=data)
qc.draw()
Collisions¶
When we approximate the intergrals with the discretized sum, the (time)-collisions in the time-discretization registers can appear. To keep record of the collisions occurence, we can use the functionality of a FlagCollisionsUnary Qubrick for the unary Dyson Series and FlagCollisionsBinary Qubrick for the binary Dyson Series. Note that we do not keep track in which of the time-discretization registers a collision occured, we only flag if at least one such an event occured.
As a small example explaining the difference between including and excuding the time collisions, consider the terms for the truncation parameter $k=3$, $k \leq K$ and the discretization parameter $m=4$, $m \leq M$. The expansion of the Dyson series without time collisions contain the terms such as
$$ W\left(\frac{2t}{M}\right)W \left( \frac{1t}{M} \right)W\left(\frac{0t}{M}\right) \tag{10} $$
while the Dyson Series expansion with time collisions has terms such as
$$ W\left(\frac{2t}{M}\right)W \left( \frac{2t}{M} \right)W\left(\frac{1t}{M}\right) + W\left(\frac{3t}{M}\right)W \left( \frac{3t}{M} \right)W\left(\frac{2t}{M}\right) \tag{11}. $$
# record at least one time collision in the discretization register
flag_t_collisions = True
if flag_t_collisions:
total_num_qubits += num_anc * 2 + 1
qc = QPU()
qc.reset(total_num_qubits)
psi = Qubits(num_system, "psi_reg", qc)
dyson_reg = DysonSeriesUnaryIndex.initialize(
qc, "DysonReg", trunc_bound, disc_bound, psi, be_size
)
data = DysonSeriesData(disc_bound, trunc_bound, flag_t_collisions, OFF_DIAG_OP, time)
trunc_reg_probs = dyson_target_prob(time, trunc_bound)
# allocation of a flag qubrick is not nessessary and can also happen inside one of the Dyson Select routines
flag_qbk = FlagCollisionsUnary()
prepare = PrepareClockState()
dyson_prep = DysonSeriesPrepareUnary(prepare)
lcu_prepare = PrepareNaive(OFF_DIAG_OP.get_coefficients())
lcu_select = SelectNaive()
lcu = LCU(lcu_prepare, lcu_select)
dyson_select = DysonSeriesSelectUnary(diag_op, lcu, flag_qbk)
dyson_op = DysonSeriesLCU(dyson_prep, dyson_select)
dyson_op.compute(dyson_reg, ctrl=0, data=data)
qc.draw()
5. Dyson Series Quantum Circuit Verification: Dyson Operator Norm and State Vector Simulation¶
But how can we verify that our circuit construction for the Dyson Series is correct? In this section we look at the Dyson Series Qubrick using the block encoding concept. We can calculate the norm of the Dyson Series operator classically, and compare classical results with measurement from the quantum circuit that represents the Dyson Series Qubrick.
If $D_{(K, M)}(t)$ represent the unitary approximation of the Dyson Series, then we can also think about it as as a block-encoding subroutine:
$$ BE(D_{(K, M)}(t))|0⟩^{\otimes n}|ψ⟩ = |0⟩^{\otimes n} D_{(K, M)}(t) |ψ⟩ + |⊥⟩ \tag{12} $$
where $|0⟩^{\otimes n}$ represents the zero space of all registers required to construct the Dyson Series Qubrick (truncation, time-disretization, block-encoded ancillaes, etc.). If we collect and measure out all these registers in the zero space with some probability $p_D$, it means that with the probability $p_D$ the Dyson operator $D_{(K, M)}(t)$ has been successfully applied to the system state $|ψ⟩$.
To verify the correctness of the Dyson Series Qubrick implementation, we can do the following steps:
- Classically, we construct the operator $D_{(K, M)}(t)$ using the
numpyPython library withing theget_expected_parametersfunction. The functiondyson_series_numpyconstructs the matrix representation of operator $D_{(K, M)}(t) = I + \sum_{k=1}^{K}{\frac{(-it)^k}{M^k k!}} \left( \sum^{M-1}_{m_k=k-1} \ldots \sum^{m_3-1}_{m_2=1} \sum^{m_2-1}_{m_1=0} W \left(\frac{m_kt}{M} \right) \ldots W\left(\frac{m_2 t}{M} \right) W \left(\frac{m_1t}{M} \right) \right)$ for the given parameters $K$, $M$, $t$. - We find the eigenvectors of the $D_{(K, M)}(t)$ operator and push any eigenvector to the quantum circuit (this is the $|ψ⟩$ or system state). As an alternative approach, we could also push some random normalized state.
- We compare the classical result of $ \lVert D_{(K, M)}(t) \cdot |ψ⟩ \rVert_2^2 $ with the probability $p_D$ of measuring all ancilae registers in the zero space.
- We also verify that the classically computed final state $D_{(K, M)}(t) |ψ⟩$ equals to the quantum state we obtain from the quantum circuit using
qc.pull_state().
Note: For a more in-depth explanation of the block-encoding concept check out the notebook on LCU and Qubitization.
def get_expected_parameters(trunc_param, t_disc_param, t0, onenorm, diag_op, be_op, num_qubits, flag_collision):
"""A helper function to calculate the expected state after applying a matrix representation
of the (K, M)-Dyson Series.
Args:
K (int): truncation bound.
M (int): time discretization bound.
t0 (int, float): time.
onenorm (float): L1 norm of the time-independent Hamiltonian.
diag_op (): the time-dependent Hamiltonian
be_op (): the time-independent Hamiltonian
num_qubits (int): the expected number of qubits in the circuit.
flag_collision (boolean): flag/not flag the time collisions.
Return:
expeced_prob (float): expected probability of applying the (K, M)-Dyson operattor to a system state.
initial_state ((np.array): system state the (K, M)-Dyson operattor is applied to.
final_expected_state (np.array): expected normalized state after applying (K, M)-Dyson operator.
"""
# construct the matrix representation of D(K, M) operator
be_op_matrix = pauli_sum_to_numpy(be_op)
final_op = dyson_series_numpy(trunc_param, t_disc_param, t0 / onenorm, diag_op, be_op_matrix, flag_collision)
final_norm = norm_beta(t0, trunc_param)
op = reverse_numpy_op(final_op) / final_norm
# generate an initial state randomly
seed = 42
random.seed(seed)
num_system = max(be_op.width(), 1)
initial_state_np = [random.random() for _ in range(1 << num_system)]
initial_state_norm = initial_state_np / np.linalg.norm(initial_state_np)
initial_state = np.kron(np.eye(1, 2 ** (num_qubits - num_system)), initial_state_norm)[0]
final_expected_state = np.kron(np.eye(1, 2 ** (num_qubits - num_system)), (op @ initial_state_norm))[0]
expeced_prob = np.linalg.norm(op @ initial_state_norm) ** 2
final_expected_state = final_expected_state / np.linalg.norm(final_expected_state)
return expeced_prob, initial_state, final_expected_state
# classically compute the norm and the final state
expeced_prob, initial_state, expected_state = get_expected_parameters(
trunc_bound, disc_bound, time, onenorm, diag_matrix, OFF_DIAG_OP, total_num_qubits, flag_t_collisions
)
# construct a quantum circuit representing the Dyson series
qc = QPU()
qc.reset(total_num_qubits)
psi = Qubits(num_system, "psi_reg", qc)
dyson_reg = DysonSeriesUnaryIndex.initialize(
qc, "DysonReg", trunc_bound, disc_bound, psi, be_size
)
qc.push_state(initial_state) # push random state as a quantum state
# construct Dyson Series circuit
data = DysonSeriesData(disc_bound, trunc_bound, True, OFF_DIAG_OP, time)
trunc_reg_probs = dyson_target_prob(time, trunc_bound)
flag_qbk = FlagCollisionsUnary()
prepare = PrepareClockState()
dyson_prep = DysonSeriesPrepareUnary(prepare)
lcu_prepare = PrepareNaive(OFF_DIAG_OP.get_coefficients())
lcu_select = SelectNaive()
lcu = LCU(lcu_prepare, lcu_select)
dyson_select = DysonSeriesSelectUnary(diag_op, lcu, flag_qbk)
dyson_op = DysonSeriesLCU(dyson_prep, dyson_select)
dyson_op.compute(dyson_reg, ctrl=0, data=data)
# collect all registers (truncation, discretization, block-encoding ancillae, flag)
# postselect, and get the system state back
flag = flag_qbk.get_result_qreg()
all_registers = 0
all_registers |= dyson_reg.trunc_reg
for reg in dyson_reg.be_anc:
all_registers |= reg
for reg in dyson_reg.discret_reg:
all_registers |= reg
all_registers |= flag
# total_prob must be equal to ||D(K,M) |psi> || ** 2
actual_prob = all_registers.peek_read_probability(0)
# compare probabilities from quantum circuit and classical computed
print(f"Expected probability is {expeced_prob}")
print(f"Actual probability is {actual_prob}")
# compare quantum states from quantum circuit and classical computed
# postselect (partial measurement) to get the system register
all_registers.postselect(0)
actual_state = qc.pull_state()[:1 << num_system]
expected_state = expected_state[: 1 << num_system]
print(f"Expected final state is {expected_state}")
print(f"Actual final state is {actual_state}")
Expected probability is 0.6155727155960966 Actual probability is 0.6155727155960967 Expected final state is [0.97862678-0.00407587j 0.1677214 -0.11892244j] Actual final state is [0.97862678-0.00407587j 0.1677214 -0.11892244j]
Success 🙌! The quantum circuit for the Dyson Series produces all expected results.
6. Dyson Series with Binary Truncation Register in WBA¶
The Dyson Series with the binary version of the truncation register utilizes similar to the Dyson Series with a unary register structure of the $Prep \cdot Sel \cdot Prep ^{\dagger}$, except some modifications in $Prep$ and $Sel$ subroutines themselves:
- The coefficients of the Taylor series expansion $p_k := \frac{t^k}{\beta k!}$ are represented in a binary encoding. We use here available in the WBA
PrepareNaivestate preparation subroutine to encode these coefficients into a quantum register. Since we use a binary version of a truncation register, its size is $\left\lceil \log_2{(K + 1)} \right\rceil + 1$, in contrast to $K$ for the truncation register size in the unary version of the Dyson series. The extra qubit is required to use a compression gadget subroutine for marking the zeroe subspace of the block-encoding of the product of operators $BE(V_n \cdots V_m)$, where $n \in [K, 1]$ and $m \in [K-1, 0]$. For more details, please refer to Appendix B of the research paper Hamiltonian Simulation in the Interaction Picture ⧉ by G.H. Low and N. Wiebe. - Each of the block-encoding of an off-diagonal operators $BE(V_i)$ uses the same zero subspace due to the Compression Gadget subroutine. In the unary version of the Dyson Series, each of the $BE(V_i)$ uses separate qubit(s). The logic of the Compression Gadget subroutine correctly identifies the block encoding of the product of the operators $BE(V_l \cdot V_{l-1} \ldots V_1)$. It also allows to reduce the qubits overhead in favour of extra quantum gates.
In this section we repeat the experiemnt with the norm and the system register for the binary Dyson Series. In fact, we are going to construct the binary Dyson Series for the same truncation, time-discretization, time parameters and the Hamiltonian and compare the results.
flag_t_collisions = True
num_system = max(OFF_DIAG_OP.width(), 1)
be_size = int(np.ceil(np.log2(len(OFF_DIAG_OP))))
num_comps = total_sort_num_comps(trunc_bound)
num_dyson_idx_reg = num_qubits_DysonSeriesExpansionIndex(
trunc_bound, disc_bound, num_system, be_size, reg_type="binary"
)
num_anc = int(np.ceil(np.log2(disc_bound)))
total_num_qubits = num_dyson_idx_reg + num_comps + num_anc + 4
if flag_t_collisions:
total_num_qubits += num_anc * 2 + 1
# get DysonBinary norm classically
expected_prob, initial_state, expected_state = get_expected_parameters(
trunc_bound, disc_bound, time, onenorm, diag_matrix, OFF_DIAG_OP, total_num_qubits, flag_t_collisions
)
# Construct DysonSeriesBinary quantum circuit
qc = QPU()
qc.reset(total_num_qubits)
psi = Qubits(num_system, "psi_reg", qc)
dyson_reg = DysonSeriesBinaryIndex.initialize(
qc, "DysonReg", trunc_bound, disc_bound, psi, be_size
)
qc.push_state(initial_state) # push one of the eigenvectors as a quantum state
data = DysonSeriesData(disc_bound, trunc_bound, flag_t_collisions, OFF_DIAG_OP, time)
lcu_prepare = PrepareNaive(OFF_DIAG_OP.get_coefficients())
lcu_select = SelectNaive()
lcu = LCU(lcu_prepare, lcu_select)
flag_qbk_binary = FlagCollisionsBinary()
naive_strategy = DysonSelectStrategyNaive(trunc_bound, diag_op, OFF_DIAG_OP, lcu)
dyson_select = DysonSeriesSelectBinary(naive_strategy, flag_qbk_binary)
prepare = PrepareNaive(data.probabilities())
dyson_prep = DysonSeriesPrepareBinary(prepare)
dyson_op = DysonSeriesLCU(dyson_prep, dyson_select)
dyson_op.compute(dyson_reg, ctrl=0, data=data)
flag_reg = flag_qbk_binary.get_result_qreg() if flag_t_collisions else 0
all_ancillae = 0
all_ancillae |= dyson_reg.trunc_reg
for reg in dyson_reg.be_anc:
all_ancillae |= reg
for reg in dyson_reg.discret_reg:
all_ancillae |= reg
all_ancillae = all_ancillae | flag_reg
all_ancillae |= dyson_op.get_result_qreg()
total_prob = all_ancillae.peek_read_probability(0)
actual_probs = all_ancillae.peek_read_probability(0)
# compare probabilities from quantum circuit and classical computed
print(f"Expected probability is {expected_prob}")
print(f"Actual probability is {actual_probs}")
# compare quantum states from quantum circuit and classical computed
all_ancillae.postselect(0)
actual_state = qc.pull_state()[:1 << num_system]
expected_state = expected_state[: 1 << num_system]
print(f"Expected final state is {expected_state}")
print(f"Actual final state is {actual_state}")
qc.draw(show_qubricks=False)
Expected probability is 0.6155727155960966 Actual probability is 0.6155727155960967 Expected final state is [0.97862678-0.00407587j 0.1677214 -0.11892244j] Actual final state is [0.97862678-0.00407587j 0.1677214 -0.11892244j]
We can observe that we are getting the same norm and the final state vector as in unary version of the Dyson Series circuit which concludes the equivalnce of both qubricks.
7. Quantum Circuits Comparison¶
We can do a high level comparison of the quantum circuit structures for both versions of the Dyson operators.
# This code is to geerate the high-level circuit comparison
def generate_unary_cir():
qc = QPU()
qc.reset(10)
psi = Qubits(2, "psi_reg", qc)
be_anc = Qubits(2, "BE_ancillae", qc)
trunc_reg = Qubits(2, "trunc_reg", qc)
discret_reg = Qubits(2, "t_discr_reg", qc)
flag_collisions = Qubits(2, "flag_collisions", qc)
qc.label("DysonPrepareUnary")
trunc_reg.box_open("UnaryPrepare")
qc.nop(repeat=3)
trunc_reg.box_close("UnaryPrepare")
discret_reg.box_open("Entangle", cond=trunc_reg)
qc.nop(repeat=2)
discret_reg.box_open("Entangle", cond=trunc_reg)
discret_reg.box_open("BitonicSort")
qc.nop(repeat=2)
discret_reg.box_close("BitonicSort")
qc.label()
qc.label("DysonSelectUnary")
trunc_reg.box_open("Phases")
qc.nop()
trunc_reg.box_close("Phases")
(psi | be_anc).box("V(t)", cond=~trunc_reg|discret_reg)
(psi | be_anc).box("V(t)", cond=trunc_reg[0] | ~trunc_reg[1]|discret_reg)
flag_collisions.box_open("FlagCollisions", cond=trunc_reg|~discret_reg)
qc.nop(repeat=2)
flag_collisions.box_close("FlagCollisions", cond=trunc_reg|~discret_reg)
qc.label()
qc.label("DysonPrepareUnaryUncompute")
discret_reg.box_open("BitonicSortUncomp")
qc.nop(repeat=4)
discret_reg.box_close("BitonicSortUncomp")
discret_reg.box_open("EntangleUncomp", cond=trunc_reg)
qc.nop(repeat=3)
discret_reg.box_open("EntangleUncomp", cond=trunc_reg)
trunc_reg.box_open("UnaryPrepareUncomp")
qc.nop(repeat=5)
trunc_reg.box_close("UnaryPrepareUncomp")
qc.label()
qc.draw("unary_demo.svg")
def generate_binary_cir():
qc = QPU()
qc.reset(12)
psi = Qubits(2, "psi_reg", qc)
be_anc = Qubits(1, "BE_ancillae", qc)
trunc_reg = Qubits(2, "trunc_reg", qc)
discret_reg = Qubits(2, "t_discr_reg", qc)
flag_collisions = Qubits(2, "flag_collisions", qc)
CG_counter = Qubits(2, "CompGadget_reg", qc)
qc.label("DysonPrepareBinary")
trunc_reg.box_open("BinaryPrepare")
qc.nop(repeat=3)
trunc_reg.box_close("BinaryPrepare")
discret_reg.box_open("BitonicSort")
qc.nop(repeat=2)
discret_reg.box_close("BitonicSort")
qc.label()
qc.label("DysonSelectBinary")
trunc_reg.box_open("Phases")
qc.nop()
trunc_reg.box_close("Phases")
flag_collisions.box_open("FlagCollisions", cond=trunc_reg|~discret_reg)
qc.nop(repeat=2)
flag_collisions.box_close("FlagCollisions", cond=trunc_reg|~discret_reg)
(psi | be_anc).box("V(t)", cond=~trunc_reg|discret_reg)
qc.box_open("CompressionGadget")
CG_counter.box_open("CG_Flag", cond=~be_anc)
qc.nop()
CG_counter.box_close("CG_Flag", cond=~be_anc)
trunc_reg.box_open("CG_Adder", cond=~be_anc)
qc.nop()
trunc_reg.box_close("CG_Adder", cond=~be_anc)
qc.box_close("CompressionGadget")
(psi | be_anc).box("V(t)", cond=trunc_reg[0] | ~trunc_reg[1]|discret_reg)
qc.box_open("CompressionGadget")
CG_counter.box_open("CG_Flag", cond=~be_anc)
qc.nop()
CG_counter.box_close("CG_Flag", cond=~be_anc)
trunc_reg.box_open("CG_Adder", cond=~be_anc)
qc.nop()
trunc_reg.box_close("CG_Adder", cond=~be_anc)
trunc_reg.box_open("CG Adder")
qc.nop()
trunc_reg.box_close("CG Adder")
qc.box_close("CompressionGadget")
qc.label()
qc.label("DysonPrepareBinaryUncompute")
discret_reg.box_open("BitonicSortUncomp")
qc.nop(repeat=4)
discret_reg.box_close("BitonicSortUncomp")
trunc_reg.box_open("BinaryPrepareUncomp")
qc.nop(repeat=5)
trunc_reg.box_close("BinaryPrepareUncomp")
qc.label()
qc.draw("binary_demo.svg")
print("\nDysonSeriesUnary Quantum circuit: \n")
generate_unary_cir()
print("\nDysonSeriesBinary Quantum circuit:\n")
generate_binary_cir()
DysonSeriesUnary Quantum circuit:
DysonSeriesBinary Quantum circuit:
