# Solving quantum equation using qiskit

I am struggling to solve the mathematical equations using qiskit in the below paper. https://arxiv.org/pdf/2205.00081 The mathematical equation I tried to solve is below:

I tried solving it using python as shown below (not sure whether it is correct or not).

import numpy as np
import pdb

n = 1 # one system qubit
a = 2  # two ancilla qubit

# Identity matrix for the n system qubits
I_n = np.eye(2**n)

zero_ket = np.array([[1], [0]])   # Column vector for |0>
zero_ket_extended = np.kron(zero_ket, zero_ket)   # |0> ⊗ |0> for two ancilla qubits

zero_bra = zero_ket.T      # <0|
zero_bra_extended = np.kron(zero_bra, zero_bra)   # <0| ⊗ <0|

ket_projector = np.kron(zero_ket_extended, I_n)   #  |0>^a ⊗ I_n
bra_projector = np.kron(zero_bra_extended, I_n)   #  <0|^a ⊗ I_n

# define a simple unitart matrix U on 2 qubits (e.g., CNOT followed by a Hadamard on the second qubit)
H = 1/np.sqrt(2)*np.array([[1, 1],
[1, -1]])
X = np.array([[0,1],
[1, 0]])
I = np.eye(2)

# Apply Hadamard to the first ancilla qubit, Identity to the second ancilla, and X to the system qubit
# Using the Kronecker product to construct the unitary U that applies H, I, and X respectively
U = np.kron(H, np.kron(I, X))

# Applying the block-encoding formula
A_tilde = bra_projector @ U @ ket_projector

# Result
print("Block-encoded operator A_tilde:")
print(A_tilde)


It would really great and kind help if someone can help me solving it using QISKITsimulator.

Edit 1

I think, my question is not clear. I am simplying question by first solving manually as shown below. I have also provided a circuit (not sure whether it is correct). Below is my python code

import numpy as np
# Adjust print options
np.set_printoptions(precision=3, suppress= True)

H = (1/np.sqrt(2))* np.array([[1, 1],[1, -1]])
I = np.array([[1,0],[0,1]])
U = np.kron(np.kron(I, H), H)
print(U)

[[ 0.5  0.5  0.5  0.5  0.   0.   0.   0. ]
[ 0.5 -0.5  0.5 -0.5  0.  -0.   0.  -0. ]
[ 0.5  0.5 -0.5 -0.5  0.   0.  -0.  -0. ]
[ 0.5 -0.5 -0.5  0.5  0.  -0.  -0.   0. ]
[ 0.   0.   0.   0.   0.5  0.5  0.5  0.5]
[ 0.  -0.   0.  -0.   0.5 -0.5  0.5 -0.5]
[ 0.   0.  -0.  -0.   0.5  0.5 -0.5 -0.5]
[ 0.  -0.  -0.   0.   0.5 -0.5 -0.5  0.5]]

zero_ket = np.array([[1], [0]])   #  |0>
I = np.array([[1,0],[0,1]])
U_a = np.kron(np.kron(zero_ket, I), I)
print(U_a)
[[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]]

zero_bra = zero_ket.T        # <0|
I = np.array([[1,0],[0,1]])
U_b = np.kron(np.kron(zero_bra, I), I)
print(U_b)
[[1 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0]
[0 0 0 1 0 0 0 0]]

A_tilde= U_b @ U @ U_a
print(A_tilde)
[[ 0.5  0.5  0.5  0.5]
[ 0.5 -0.5  0.5 -0.5]
[ 0.5  0.5 -0.5 -0.5]
[ 0.5 -0.5 -0.5  0.5]]


Please suggest me correction in my solution. Also, please help me how to solve this using qiskit statevector or qasm simulator?

## 3 Answers

The paper is about the FABLE block encoding. One way to implement this in qiskit is by first using the classiq SDK package to build the model on the functional level, synthesize to a quantum program, and then upload it to a qiskit circuit. (The codes below are in classiq 0.41.1).

## 1. Building the oracle

This is the quantum model you need to build

where the main interesting block is the oracle given in Eq. (6) in the paper: $$O_A |0\rangle_{\rm ind} |i\rangle_n |j\rangle_n = \left(a_{ij}|0\rangle_{\rm ind} + \sqrt{1-a_{ij}^2}|1\rangle_{\rm ind}\right)|i\rangle_n |j\rangle_n.$$ As described in the paper, this can be implemented as a series of multi-controlled rotations. We define a quantum function that gets the states $$i\rangle$$, $$|j\rangle$$, the indicator qubit, and the elements of the matrix as a flatten list

from classiq import *
from classiq.qmod.symbolic import acos
@qfunc
def my_encoding(entries: CArray[float], i: QNum, j:QNum, ind: QBit):
ij = QNum("ij",2*i.size,False,0)
bind([i,j], ij)
repeat(2**ij.size,
lambda k: control(ij==k,
lambda: RY(2*acos(entries[k]),ind)
)
)
bind(ij,[i, j])


This is how the oracle looks like for $$2\times 2$$ matrix:

## 2. Building the full block encoding

Now you can build the full algorithm. There are additional 2 blocks: a Hadamard transform, and a swap between $$i$$ and $$j$$. For the latter you need to define another function.

@qfunc
def swap_all(i: QArray, j: QArray):
repeat(i.len,
lambda k: SWAP(i[k],j[k])
)

@qfunc
def fable_block_encoding(entries: CArray[float], i: QNum, j: QNum, ind: QBit):
within_apply(lambda: hadamard_transform(j),
lambda: (my_encoding(entries, i, j, ind), swap_all(i,j))
)


## 3. Validation for a specific example

A validation is done by loading some quantum state $$|\vec{b}\rangle$$. Calling the block encoding for a specific matrix $$A$$, and some vector $$b$$ is done as follows:

import numpy as np

a_matrix = np.array([[1,2,5,7],[2,1,6,3],[-1,3,-4,5],[1,1,3,4]])
s_factor = np.linalg.norm(a_matrix, ord=2)
a_matrix = a_matrix/s_factor
b_vector = [7,1,2,3]
b_vector = (b_vector/np.linalg.norm(b_vector)).tolist()
entries = a_matrix.flatten().tolist()


where we normalized $$A$$ to have $$|A|\leq 1$$ (this is the $$\alpha$$ factor in the paper), and $$\vec{b}$$ must be normalized as we load it to a quantum state.

a. Now you build the model and synthesize it to a quantum program:

@qfunc
def main(i: Output[QArray], j: Output[QArray], ind: Output[QBit]):

prepare_amplitudes(b_vector,0.0,i)
allocate(i.len,j)
allocate(1,ind)
fable_block_encoding(entries, i, j, ind)

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


This is how the circuit looks like:

b. We expect to measure $$A\cdot \vec{b}$$ on the $$|i\rangle$$ register, conditioned on all the other variables being at state zero. Execution and post-selection is done as follows:

 from classiq.execution import (
ClassiqBackendPreferences,
ClassiqSimulatorBackendNames,
ExecutionPreferences,
set_quantum_program_execution_preferences
)

execution_prefs = ExecutionPreferences(
backend_preferences = ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR
)
)

qprog = set_quantum_program_execution_preferences(qprog, execution_prefs)

res = execute(qprog).result()
amps = np.zeros(len(b_vector)).astype(complex)
for sample in res[0].value.parsed_state_vector:
if sample.state["ind"] == sample.state["j"] ==0.0:
amps[int(sample.state["i"])] += sample.amplitude

global_phase = np.angle(res[0].value.parsed_state_vector[0].amplitude)

amps = np.real(amps/np.exp(1j*global_phase))
expected_result = a_matrix @ b_vector/len(b_vector)
print("measured:", amps)
print("expected:", expected_result)


This is what you find:

measured: [0.1028748 0.09258732 0.00771561 0.06686862] expected: [0.1028748 0.09258732 0.00771561 0.06686862]

## 4. Uploading to qiskit

Once you have the quantum program from classiq you can upload it to a QuantumCircuit object in qiskit as follows:

from qiskit import QuantumCircuit

program = QuantumProgram.parse_raw(qprog).transpiled_circuit.qasm
qc = QuantumCircuit.from_qasm_str(program)


Finally, the block encoding can be used to solve a system of linear equations, by a QSVT routine. See here for example.

• Wow thanks! can I use the *= Amplitude Encoding operator for Part 1? docs.classiq.io/latest/reference-manual/platform/qmod/… Commented May 29 at 7:44
• Thank you @tomergf. This looks like amazing code. Can you give this code in qiskit also?
– Manu
Commented May 30 at 4:12
• Yes @RonCohen, you can use the Amplitude Encoding assignment as well, in particular, the "Quantum Numeric Subscript Semantics", where you pass the entries variable as a lookup table (see details in the link you've sent). One drawback with this approach, compared to the one I've suggested, is that it cannot take into account the sparseness of the matrix. However, this approach is robust, and will work for any matrix. Commented Jun 4 at 11:16
• Thank you @Manu. Sorry, I do not know how to code this directly in qiskit. Again, you can get the qiskit version indirectly, by using the code above and uploading the qasm to a qiskit quantum circuit (step 4). I hope my answer and approach still answers your question and that you've found it useful and helpful. Commented Jun 4 at 11:26
• Hi @tomergf I am getting the error ClassiqAPIError: Call to API failed with code 401: Not authenticated while trying to run your code.
– Manu
Commented Jun 10 at 16:44

Using NumPy is probably the best way compute $$\tilde{A}$$. However, it is worth noting that the expression for $$\tilde{A}$$ shown below is just fancy bra-ket/matrix notation for a map that extracts from $$U$$ the first $$2^n \times 2^n$$ matrix block from it:

$$\tilde{A} = \left(\langle0|^{\otimes a} \otimes I_n \right) U \left(|0\rangle^{\otimes a} \otimes I_n \right).$$

So getting $$\tilde{A}$$ is as simple as indexing $$U$$ like this:

A_tilde = U[:2**n, :2**n]


Here's all the code you really need for the example you provided:

import numpy as np
from qiskit.quantum_info import Operator

U = Operator.from_label('HIX')
n = 1
N = 2**n

A_tilde = U.data[:N, :N]


Regarding Edit 1

Looking at the circuit diagram you drew, I am getting the impression that you are trying to treat $$\left(|0\rangle^{\otimes a} \otimes I_n \right)$$ and $$\left(|0\rangle^{\otimes a} \otimes I_n \right)$$ like unitary operators, but they are not. More formally, this calculation is a map where you are tracing out the zero-state over the last $$a$$ qubits of the $$m$$-qubit unitary $$U$$:

$$\text{Tr}_{0^{\otimes a}}[U] = \left(\langle0|^{\otimes a} \otimes I_n \right) U \left(|0\rangle^{\otimes a} \otimes I_n \right).$$

In other words, here $$\langle0| \left(\; \cdot \;\right) |0\rangle$$ is used to "remove" from $$U$$ the portion of it that acts on the subspace consisting of the last $$a$$ qubits, and $$I_n \left(\; \cdot \;\right) I_n$$ to "leave alone" the part of $$U$$ that acts on the first $$n$$ qubits.

If what you want is to implement a circuit in Qiskit to recover $$\tilde{A}$$, then why not just implement the circuit that actually applies $$\tilde{A}$$ on a state $$|\psi\rangle$$, which is shown in Fig 1 of the paper:

This circuit performs a measurement on the last $$a$$ qubits of the state, and post-selects over the all-zeros measurement. To recover $$\tilde{A}$$ you could then apply one by one the basis states over the subspace of $$|\psi\rangle$$, measure the top $$a$$ qubits, and add the output state in each of the columns of $$\tilde{A}$$. The problem is that you will not be able to get $$\tilde{A}$$ itself, but its elements up to an normalization factor $$\|\tilde{A} |\psi\rangle\|$$.

Here's an example that works OK. Imagine you want to apply a non-unitary matrix $$A$$ to $$|\psi\rangle$$ given by:

$$A = \begin{bmatrix} 1 & 1 \\ 0 & 0 \end{bmatrix}.$$

This is the set-to-zero operation, which results in state $$|0\rangle$$ independent if the input is $$|0\rangle$$ or $$|1\rangle$$. You can block-encode this unitary using a $$\text{SWAP}$$ gate followed by an $$H$$ gate in the ancillary qubit:

The unitary describing this circuit is:

\begin{aligned} U &= \text{SWAP} \, (H \otimes I) \\ \\ U &= \begin{bmatrix} \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & 0 & 0 \\ 0 & 0 & \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} \\ \frac{\sqrt{2}}{2} & - \frac{\sqrt{2}}{2} & 0 & 0 \\ 0 & 0 & \frac{\sqrt{2}}{2} & - \frac{\sqrt{2}}{2} \end{bmatrix} \end{aligned}

Which then makes $$A = \alpha \tilde{A}$$, with $$\alpha = \sqrt{2}$$.

The following code implements this circuit and extracts the output state for the two possible inputs:

from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator, partial_trace
from qiskit_aer import AerSimulator

# For loop to run circuit for |ψ⟩ = |0⟩ and then |ψ⟩ = |1⟩
for state in [0,1]:
qc = QuantumCircuit(2,1)

if state:
qc.x(0)

qc.swap(1,0)
qc.h(1)
qc.measure(1,0)
qc.save_statevector(pershot=True)
display(qc.draw())

simulator = AerSimulator()

# Implement while loop until ancilla result is '0'
anc = '1'
shots = 1
while(anc == '1'):
result = simulator.run(qc,shots=1,memory=True).result()
anc = result.get_memory()[0]

# If ancilla is '0', extract system statevector
if anc == '0':
ψ = result.data().get('statevector')[0]
ψ_s = partial_trace(ψ,[1]).to_statevector()
print(f'After {shots} shot(s), given the input state |{state}⟩, the output state is:')
display(ψ_s)

else:
shots += 1


Your output should look like this:

After 4 shot(s), given the input state |0⟩, the output state is:


$$|0\rangle$$

After 1 shot(s), given the input state |1⟩, the output state is:


$$|0\rangle$$

Note that I am using a while loop here, which for some matrices might be a bad choice. For instance, if $$A$$ was of the form:

$$A = \begin{bmatrix} 1 & 0 \\ 1 & 0 \end{bmatrix},$$

the expectation is that you should never see an output if the input is $$|1\rangle$$, which translates into never getting a $$|0\rangle$$ in the ancilla register (and therefore the while loop will get stuck).

Hope this clarifies your question.

• Thank you @diemilio for the response. I would be really great if you can explain me how you concluded that A_tilde=U[: 2**n, :2**n]? Also in your solution you have not explained about a in the mathematical equation?
– Manu
Commented May 27 at 5:31
• Hi @diemilio, I have edited my question to make it more clear. Please review.
– Manu
Commented May 27 at 7:17
• Great answer. Slicing is really the best way to go in this particular case. Commented May 27 at 8:49
• @Manu, I added a response to your edits. Commented May 28 at 10:36
• Thank you @diemilio for great answer. I still have some doubts. How is A = alpha *A_tilde with alpha= sqrt(2)? Also, when put a measurement gate in a circuit , then two lines and then 0, what exactly this means?
– Manu
Commented May 29 at 6:24

## TLDR

Yes, you are basically on the right track to calculate $$A$$. But there are some adjustments to make on how you model $$U$$. This is also where qiskit can be of help with the Operator API.

## How to model U...

### ... without qiskit

The way you define $$U$$ in your code is:

$$U = H \otimes \mathbb{I} \otimes X$$

which is not a CNOT followed by a Hadamard.

If you want to construct CNOT out of Paulis, the correct way to do so is by using

$$CX_{12} = \frac{1}{2}(1 + Z_1 + X_2 - Z_1 X_2)$$

where $$Z_1$$ is short for $$Z \otimes \mathbb{I}$$ and likewise for $$X_2$$.

In python, this would be

X = np.array([[0, 1], [1, 0]])
Z = np.array([[1, 0], [0, -1]])
CX = 1/2 * (np.eye(4) + np.kron(Z, I) + np.kron(I, X) - np.kron(Z, X))


Alternatively, you could also define CX using its numerical matrix representation.

If you want CNOT to be followed by H then you have to contract them. For example:

H_2 = np.kron(np.kron(I, H), I)
C_12 = np.kron(CX, I)
U = H_2 @ C_12


From your in-code comment, it is not clear to me on which of the 3 qubits you want your CNOT to act on (I just assumed qubits 1 and 2).

### ... with qiskit

With qiskit, it becomes a lot easier to model $$U$$, but you have to be careful with the ordering of the qubits, since qiskit counts the qubits from the right (in terms of bra- and ket-representation) which is also reflected in the matrix elements.

For example, if you want to produce the same $$U$$ as described above you could do the following:

circuit = QuantumCircuit(3)
circuit.cx(2, 1)
circuit.h(1)
u = Operator(circuit).data


which looks a lot simpler, but notice that I intentionally reversed the ordering of the qubits.

If you do not want to reverse the ordering of the qubits in qiskit you have to adjust your definitions of ket_projector and bra_projector accordingly.

Except of those adjustments your approach looks fine.

## Concerning Edit 1

In your edit 1, you defined $$U = \mathbb{I} \otimes H \otimes H$$. And your python calculations correctly yield

$$A = \langle 0| U |0\rangle = \langle 0| \mathbb{I} |0\rangle H \otimes H = H \otimes H = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & - 1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & - 1 & -1 & 1 \end{pmatrix}$$

As I said, the definition of $$U$$ can be effectively done with qiskit. In the case of $$U = \mathbb{I} \otimes H \otimes H$$, it is even a one-liner (like diemilio already pointed out in his answer):

U = Operator.from_label('IHH')

But as far as I know, qiskit cannot do the rest for you - the partial contraction with a bra and a ket vector. Though, in your special case of contracting with $$\langle 0|^{\otimes a}$$ and $$|0\rangle^{\otimes a}$$, diemilios one-line solution of calculating $$A$$ out of $$U$$ using slicing is very elegant and I do recommend it.

• Thank you @qubitzer for the response.
– Manu
Commented May 27 at 5:30
• Hi @qubitzer. I have edited my question to make it more clear. Please review,
– Manu
Commented May 27 at 7:16