Linear Combination of Unitaries and Qubitization¶
$$ \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} $$
In this tutorial, we'll review block encodings using two techniques: linear combination of unitaries (LCU) and qubitization. We'll also demonstrate WBA Qubricks that implement these techniques.
For a more in-depth review of these techniques, please refer to "Lecture Notes on Quantum Algorithms for Scientific Computation" by Lin Lin (arXiv:2201.08309 ⧉).
from psiqworkbench import Qubits, QPU
from workbench_algorithms import PrepareNaive, SelectNaive, SawtoothSelect, LCU, Qubitization
from workbench_algorithms.utils.paulimask import PauliSum, PauliMask, pauli_sum_to_numpy
import numpy as np
Block encoding¶
Block encoding is a method for encoding a non-unitary matrix $A$ as a "block" inside a larger unitary matrix $U_A$. Block encoding allows us to apply the operation $A$ to the target quantum register probabilistically by applying the unitary matrix $U_A$ to a larger quantum register composed of the target register and some auxiliary qubits.
Generally speaking, the matrix of $U_A$ will look something like this:
$$ U_A = \begin{pmatrix} A & * \\ * & * \end{pmatrix} $$
When we want to apply the non-unitary operator $A$ to the state $\ket{\psi}$, we will instead apply the unitary operator $U_A$ to a larger state, $\ket{\psi}$ padded with auxiliary qubits. The resulting state will be
$$ U_A \left( \ket{0}^{\otimes M}\ket{\psi} \right) = \ket{0}^{\otimes M} A \ket{\psi} + \ket{\perp} $$
Here $\ket{\perp}$ is a quantum state that is not normalized and which is completely orthogonal to our desired state $\ket{0}^{\otimes M} A \ket{\psi}$.
Linear combination of unitaries¶
In this tutorial, we will look at a specific method of block encoding called linear combination of unitaries (LCU). This protocol consists of two main unitary subroutines:
- $\text{PREPARE}$ is a unitary operation that, when acting on the $\ket{0}$ state, prepares an auxiliary qubit register in a state that encodes arbitrary real-valued coefficients as the amplitudes of the computational basis states.
- $\text{SELECT}$ is used to apply unitary operators to the target register, controlled on the state of the auxiliary register.
As we go through the protocol, we will see how these two subroutines work in tandem to implement LCU block encoding.
Problem setup¶
As the name implies, for the LCU method we write $A$ as a linear combination of several unitaries with real, positive-valued coefficients:
$$ A = \sum_{i=0}^{K-1} \alpha_i U_i $$
Here $\alpha_i \in \mathbb{R}$, $\alpha_i > 0$, and all $U_i$ are unitary. It is common to choose $U_i$ from the set of Pauli matrices $I, X, Y, Z$.
Let's consider a single-qubit matrix $A$ represented as the following linear combination of unitaries:
$$A = I + 2X + 3Z = \begin{bmatrix} 4 & 2 \\ 2 & -2 \end{bmatrix}$$
In Workbench, we can represent this matrix as a PauliSum object:
prep_probs = [1, 2, 3]
A = PauliSum(
[prep_probs[0], PauliMask(0, 0)],
[prep_probs[1], PauliMask(1, 0)],
[prep_probs[2], PauliMask(0, 1)]
)
A_matrix = pauli_sum_to_numpy(A)
print(A_matrix)
# Confirm that the matrix is not unitary
assert not np.allclose(A_matrix.T.conj() @ A_matrix, np.eye(2))
[[ 4.+0.j 2.+0.j] [ 2.+0.j -2.+0.j]]
$\text{PREPARE}$¶
As mentioned above, $\text{PREPARE}$ is a routine that prepares a quantum state that encodes a list of arbitrary real numbers as the amplitudes of the basis states.
In the case of LCU, we want to prepare the auxiliary register in a state where the amplitudes encode the coefficients $\alpha_i$ (the coefficients of the unitaries in the linear combination that represents our operator $A$).
For our toy example, $A = I + 2X + 3Z$, so $\alpha_0 = 1, \alpha_1 = 2, \alpha_2 = 3$. If we tried to prepare a state in which amplitudes were square roots of these coefficients, it would look as follows:
$$ \sqrt{\alpha_0} \ket{0} + \sqrt{\alpha_1} \ket{1} + \sqrt{\alpha_2} \ket{2} $$
You can see that this state is unphysical, since the sum of these coefficients $\alpha_0 + \alpha_1 + \alpha_2 = 6$ is greater than $1$.
In order to correct this, we normalize the state we want to prepare (and the operator $A$ itself) by introducing a new set of coefficients $\alpha_i^\prime = \frac{\alpha_i}{\alpha}$, where the normalization coefficient $\alpha = \sum_{i=0}^{K-1} \alpha_i$. This means that $\sum_{i=0}^{K-1} \alpha_i^\prime = 1$, and we can use the normalized coefficients as the probabilities associated with basis states.
As a result, we will actually be applying this normalized version of $A$, which we will call $A^\prime$. The implication here is that the probability at which we apply $A$ onto our target system (the success probability) will be proportional to $\frac{1}{\alpha^2}$.
For our toy example, the normalization would look as follows:
$$\alpha = 1 + 2 + 3 = 6$$
$$A^\prime = \tfrac16 (I + 2X + 3Z)$$
$$\alpha_0^\prime = \tfrac16, \alpha_1^\prime = \tfrac26 = \tfrac13, \alpha_2^\prime = \tfrac36 = \tfrac12$$
Therefore $\text{PREPARE}$ should prepare the auxiliary register in the following state:
$$\text{PREPARE}\ket{0} = \sqrt{\alpha_0^\prime} \ket{0} + \sqrt{\alpha_1\prime} \ket{1} + \sqrt{\alpha_2\prime} \ket{2}$$ $$ = \sqrt{\tfrac16}\ket{0} + \sqrt{\tfrac13}\ket{1} + \sqrt{\tfrac12}\ket{2}$$ $$ = 0.4082\ket{0} + 0.5774\ket{1} + 0.7071\ket{2}$$
Workbench Algorithms offers multiple Qubricks that implement the $\text{PREPARE}$ routine. The following code snippet shows how to use one of them, PrepareNaive, to prepare the required state.
You can see that we pass to the Qubrick constructor the coefficients $\alpha_i$ of the original LCU decomposition, and it takes care of the normalization under the hood. The amplitudes of the prepared state are exactly the $\sqrt{\alpha_i^\prime}$ that we need.
# Initialize the QPU and registers
qpu = QPU(num_qubits=3)
aux_reg = Qubits(2, "aux_reg", qpu)
# Initialize the Qubrick
prep = PrepareNaive(prep_probs)
# Prepare the auxiliary register
prep.compute(aux_reg)
qpu.print_state_vector()
qpu.draw(show_qubricks=True)
|aux_reg|?> |0|.> 0.408248+0.000000j |1|.> 0.577350+0.000000j |2|.> 0.707107+0.000000j
$\text{SELECT}$¶
The next step is the $\text{SELECT}$ routine, used to apply unitary operators to the target register, controlled on the state of the auxiliary register. This is where we bring in our target register, $\ket{\psi}$, and apply the unitary operators comprising $A$ to $\ket{\psi}$ controlled on the auxiliary register that we prepared using $\text{PREPARE}$.
Together, these two routines act to prepare the following state:
$$\ket{0}_{aux} \ket{\psi}_{target} \xrightarrow{\text{PREPARE}} \sum_{i=0}^{K-1} \sqrt{\alpha_i^\prime}\ket{i}_{aux} \ket{\psi}_{target} $$ $$\xrightarrow{\text{SELECT}} \sum_{i=0}^{K-1} \sqrt{\alpha_i^\prime}\ket{i}_{aux} U_i\ket{\psi}_{target} $$
In this state, each of the unitary operators $U_i$ is applied to the target register if the auxiliary register is in the associated basis state $\ket{i}$. The coefficient for this term in the superposition is $\sqrt{\alpha_i^\prime}$, the amplitude of that basis state in the auxiliary register. This means that we apply the unitary operator to the state $\ket{\psi}$ with an amplitude proportional to the square root of the coefficient of that term in the definition of $A$.
For our toy example, after applying the $\text{SELECT}$ routine, we expect to get something the following state:
$$0.4082\ket{0} I \ket{\psi} + 0.5774\ket{1} X \ket{\psi} + 0.7071\ket{2} Z \ket{\psi}$$
Assuming we begin with $\ket{\psi} = \ket{0}$, this should give us the explicit state:
$$0.4082\ket{0} I \ket{0} + 0.5774\ket{1} X \ket{0} + 0.7071\ket{2} Z \ket{0}$$ $$= \left( 0.4082\ket{0} + 0.7071\ket{2} \right) \ket{0} + 0.5774\ket{1} \ket{1}$$
The following code snippet continues the previous example by allocating the target register and using SelectNaive Qubrick to apply the $\text{SELECT}$ routine.
You can see that the resulting state matches the state we expect.
# Initialize the target registers
psi_reg = Qubits(1, 'psi_reg', qpu)
# Initialize the Qubricks
select = SelectNaive()
select.compute(aux_reg, psi_reg, A)
qpu.print_state_vector()
qpu.draw(show_qubricks=True)
|aux_reg|psi_reg> |0|0> 0.408248+0.000000j |1|1> 0.577350+0.000000j |2|0> 0.707107+0.000000j
Un-$\text{PREPARE}$¶
The last step of the LCU protocol is to uncompute preparation of the auxiliary register using adjoint of the $\text{PREPARE}$ routine we used on the first step. This will bring the system to the state we're aiming for with block encoding:
$$\ket{0}^{\otimes M} A \ket{\psi} + \ket{\perp}$$
In other words, if the auxiliary register is in the $\ket{0}$ state, the target register is in the (normalized) state $A \ket{\psi}$.
The following code snippet adds uncomputation to our ongoing example. You can see that, indeed, if we only consider the terms for which aux_reg is in the $\ket{0}$ state, the main register psi_reg ends up in the state that is exactly the (normalized) first column of the matrix $A$ (the amplitude of $\ket{0}$ is twice the amplitude of $\ket{1}$).
prep.uncompute()
qpu.print_state_vector()
qpu.draw(show_qubricks=True)
|aux_reg|psi_reg> |0|0> 0.666667+0.000000j |0|1> 0.333333+0.000000j |1|0> -0.235702+0.000000j |1|1> 0.235702+0.000000j |2|0> 0.333333+0.000000j |2|1> -0.333333-0.000000j |3|0> 0.235702+0.000000j |3|1> -0.235702+0.000000j
Now, we can calculate the success probability of applying our non-unitary operator $A$ as the probability to measure the auxiliary register in the $\ket{0}$ state, and compare it with the expected success probability
$$ p(\ket{aux} = \ket{0}) = \bra{\psi}{A^\prime}^\dagger A^\prime\ket{\psi} $$
The following code snippet shows how to get both probabilities and confirm that they are equal.
print(f"Success probability: {aux_reg.peek_read_probability(0) * 100:.2f}%")
A_prime_matrix = pauli_sum_to_numpy(A) / A.norm()
psi_init = np.array([1, 0])
psi_final = (A_prime_matrix @ psi_init)
overlap = np.dot(psi_final.T.conj(), psi_final)
print(f"Expected success probability: {overlap.real * 100:.2f}%")
Success probability: 55.56% Expected success probability: 55.56%
LCU Qubrick¶
Workbench Algorithms offers a handy Qubrick called LCU which does all the heavy lifting we went through in the earlier sequence of examples.
The LCU Qubrick implements block encoding of an operator $A$, using arbitrary $\text{PREPARE}$ and $\text{SELECT}$ routines as building blocks.
The following example shows how to use it to reproduce the result of the manually constructed LCU protocol using SawtoothSelect Qubrick as the $\text{SELECT}$ routine instead of SelectNaive. Now both the target register and the auxiliary register are passed to the compute method of the LCU Qubrick instead of the $\text{PREPARE}$ and $\text{SELECT}$ Qubricks individually, as well as the operator data A. You can see that the resulting state is the same as in our first example.
# Instantiate QPU and Qubits objects
qpu = QPU(num_qubits=4)
aux_reg = Qubits(2, 'aux_reg', qpu)
psi_reg = Qubits(1, 'psi_reg', qpu)
prep = PrepareNaive(prep_probs)
select = SawtoothSelect()
block_encoding = LCU(prep, select)
block_encoding.compute(psi_reg, aux_reg, data=A)
qpu.print_state_vector()
qpu.draw(show_qubricks=True)
|aux_reg|psi_reg|?> |0|0|.> 0.666667+0.000000j |0|1|.> 0.333333+0.000000j |1|0|.> -0.235702+0.000000j |1|1|.> 0.235702+0.000000j |2|0|.> 0.333333+0.000000j |2|1|.> -0.333333-0.000000j |3|0|.> 0.235702+0.000000j |3|1|.> -0.235702-0.000000j
Validating block encoding via getting the unitary matrix¶
Another way to check that our code actually implemented block encoding is via the use of UnitaryMatrixFilter to get the unitary matrix of the transformation done by the code. The following example shows how to use this filter in the code and to analyze its return.
UnitaryMatrixFiltercan be used to get the entire unitary implemented by the quantum circuit and then parse and inspect it to validate the circuit construction. However, since this requires constructing and storing the full unitary explicitly, this operation comes at a large computational cost and should only be used for debugging/validating small systems.
qpu = QPU(num_qubits=3, filters=['>>unitary>>'])
# Allocate auxiliary register after target register,
# so that aux_reg=0 will correspond to top left block of the matrix
# as returned by Workbench tools (little-endian encoding)
psi_reg = Qubits(1, 'psi_reg', qpu)
aux_reg = Qubits(2, 'aux_reg', qpu)
prep = PrepareNaive(prep_probs)
select = SelectNaive()
block_encoding = LCU(prep, select)
block_encoding.compute(psi_reg, aux_reg, A)
# Get the unitary matrix
ufilter = qpu.get_filter_by_name('>>unitary>>')
# We know the matrix is real-valued, so print only the real parts
full_matrix = ufilter.get().real
print(full_matrix.round(3))
print(A_prime_matrix.real.round(3))
# Check that the upper left block of the full matrix matches A'
upper_left_block = full_matrix[:2, :2]
assert np.allclose(upper_left_block, A_prime_matrix)
[[ 0.667 0.333 -0.236 0.236 0.333 -0.333 0.236 -0.236] [ 0.333 -0.333 0.236 -0.236 -0.333 -0.667 -0.236 0.236] [-0.236 0.236 0.833 0.167 0.236 -0.236 0.167 -0.167] [ 0.236 -0.236 0.167 0.833 -0.236 0.236 -0.167 0.167] [ 0.333 -0.333 0.236 -0.236 0.667 0.333 -0.236 0.236] [-0.333 -0.667 -0.236 0.236 0.333 -0.333 0.236 -0.236] [ 0.236 -0.236 0.167 -0.167 -0.236 0.236 0.833 0.167] [-0.236 0.236 -0.167 0.167 0.236 -0.236 0.167 0.833]] [[ 0.667 0.333] [ 0.333 -0.333]]
Qubitization¶
Qubitization is another form of block encoding. In this protocol, instead of applying the operator $A$ to our target state probabilistically, we apply Chebyshev polynomials of our operator $A$.
We won't go into details of the algorithm here and refer to Section 7.4 of "Lecture Notes on Quantum Algorithms for Scientific Computation" by Lin Lin (arXiv:2201.08309 ⧉).
Instead, we'll demonstrate the Qubitization Qubrick which allows us to construct a qubitization iterate easily.
The following example shows how to use Qubitization Qubrick. This Qubrick relies on another block encoding Qubrick as a building block; in our case, we'll use LCU Qubrick for this. Otherwise, using Qubitization Qubrick is similar to using LCU Qubrick directly.
qpu = QPU(num_qubits=3, filters=['>>buffer>>', '>>unitary>>'])
psi_reg = Qubits(1, 'psi_reg', qpu)
aux_reg = Qubits(2, 'aux_reg', qpu)
prep = PrepareNaive(prep_probs)
select = SelectNaive()
block_encoding = LCU(prep, select)
qubitization = Qubitization(block_encoding)
qubitization.compute(psi_reg, aux_reg, A)
# Get the unitary matrix
ufilter = qpu.get_filter_by_name('>>unitary>>')
# We know the matrix is real-valued, so print only the real parts
full_matrix = ufilter.get().real
print(full_matrix.round(3))
print(A_prime_matrix.real.round(3))
qpu.draw(show_qubricks=True)
# Check that the upper left block of the full matrix matches A'
# (up to a global phase of -1)
upper_left_block = full_matrix[:2, :2]
assert np.allclose(upper_left_block, -A_prime_matrix)
[[-0.667 -0.333 0.236 -0.236 -0.333 0.333 -0.236 0.236] [-0.333 0.333 -0.236 0.236 0.333 0.667 0.236 -0.236] [-0.236 0.236 0.833 0.167 0.236 -0.236 0.167 -0.167] [ 0.236 -0.236 0.167 0.833 -0.236 0.236 -0.167 0.167] [ 0.333 -0.333 0.236 -0.236 0.667 0.333 -0.236 0.236] [-0.333 -0.667 -0.236 0.236 0.333 -0.333 0.236 -0.236] [ 0.236 -0.236 0.167 -0.167 -0.236 0.236 0.833 0.167] [-0.236 0.236 -0.167 0.167 0.236 -0.236 0.167 0.833]] [[ 0.667 0.333] [ 0.333 -0.333]]
You'll notice that the matrix of qubitization iterate doesn't match that produced by LCU, though the operator $A$ is embedded in both (in the case of qubitization, up to a global phase). LCU block encoding implements a reflection, whereas qubitization iterate implements a rotation.
In practice, this means that successive applications of LCU block encodings will cancel each other out, whereas back-to-back qubitization iterates will perform non-trivial transformations of the block-encoded matrix.
Summary¶
In this tutorial, we explored block encodings and Workbench Algorithms Qubricks that implement them.
LCUQubrick implements a general block encoding using linear combination of unitaries.QubitizationQubrick implements qubitization iterate based on a block encoding Qubrick (such asLCU).