# How do I get the Unitary matrix of a circuit without using the 'unitary_simulator'?

I am using jupyter notebook and qiskit. I have a simple quantum circuit and I want to know how to get the unitary matrix of the circuit without using 'get_unitary' from the Aer unitary_simulator. i.e.:By just using matrix manipulation, how do I get the unitary matrix of the circuit below by just using numpy and normal matrix properties?

The result should equal this:

[[1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

[0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j]

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j]

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j]

[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j]

[0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]

[0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]]

# I tried doing the following but it did not result in the correct Unitary matrix

swapcnot = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])

layer1 = np.kron( swapcnot,np.eye(2) )

layer2 = np.kron( np.eye(2),swapcnot )

print( np.matmul(layer2,layer1) )

# The result:

[[1. 0. 0. 0. 0. 0. 0. 0.]

[0. 0. 0. 0. 0. 0. 0. 1.]

[0. 0. 0. 0. 0. 0. 1. 0.]

[0. 1. 0. 0. 0. 0. 0. 0.]

[0. 0. 0. 0. 1. 0. 0. 0.]

[0. 0. 0. 1. 0. 0. 0. 0.]

[0. 0. 1. 0. 0. 0. 0. 0.]

[0. 0. 0. 0. 0. 1. 0. 0.]]

• See this question, which isn't technical but outlines the procedure you're looking for: quantumcomputing.stackexchange.com/questions/15651/… – rjh324 Apr 3 at 16:36
• @ryanhill1 Thank you. i tried this but could not get it to work unfortunately – Jared Apr 3 at 19:33

For the first layer of your circuit, compute the tensor product between the unitary matrix of the (swapped) CNOT gate and the identity matrix (using numpy's kron()). Do a similar operation for the second layer. You will obtain two 8x8 matrices. Then multiply them using numpy's matmul().

Here you have the working code:

import numpy as np

swapcnot = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])

layer1 = np.kron(np.eye(2),swapcnot )

layer2 = np.kron( swapcnot, np.eye(2) )

print( np.matmul(layer2,layer1) )

Output:

[[1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 1. 0.]
[0. 0. 0. 0. 0. 1. 0. 0.]
[0. 0. 0. 0. 1. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 1.]
[0. 0. 1. 0. 0. 0. 0. 0.]
[0. 1. 0. 0. 0. 0. 0. 0.]]
• Thank you. I tried this but as you can see in the edit of my question I get a 8*8 Identity matrix – Jared Apr 3 at 19:35
• The error in your code is that you do the same tensor product for both layers, but the layers are not equal. Try with layer2 = np.kron( np.eye(2),cnot ). Morever, you have to invert the order of the layers in matmul(), as you are looking for the matrix $U_2U_1$. – Michele Amoretti Apr 3 at 20:33
• I just made these changes as seen in the above edit but I still do not get the correct result – Jared Apr 3 at 20:43
• Ok, you are using a cnot with the control and target that are exchanged, so you have to use the following matrix: swapcnot = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]) – Michele Amoretti Apr 3 at 21:01
• Ok I tried it with the new swapcnot matrix but it still is not working for some reason – Jared Apr 3 at 21:13

It seems you need the Operator class, from the quantum_info module:

from qiskit import QuantumCircuit

circuit = QuantumCircuit(3)
circuit.cx(0, 1)
circuit.cx(1, 2)
circuit.draw('latex')

from qiskit.quantum_info import Operator
from qiskit.visualization.array import array_to_latex

array_to_latex(Operator(circuit))

$$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix}$$

• Thank you for replying! Unfortunately I am trying to figure it out by not using any simulators, only algebra and normal matrix manipulation by using things such as tensor product – Jared Apr 3 at 19:39

The following script does the work, but it is based on the statevector_simulator. I do not know if it is part of what you want to avoid using.

import qiskit
from qiskit import QuantumCircuit, Aer
import numpy as np

#We need the statevector simulatro in this case
simulator = Aer.get_backend('statevector_simulator')

#Any circuit you want to analyse
qc = QuantumCircuit(3)
qc.cx(0,1)
qc.cx(1,2)

#Setting up unitary matrix
for reg in qc.qregs:
nbits=len(reg)
H = np.zeros((2**nbits, 2**nbits), dtype=complex)

#Control quantum circuit
qc_test = QuantumCircuit(nbits)
#testing all initialize vector and filling up H
for i in range(2**nbits):
init=list(range(2**nbits))
for j in range(2 ** nbits):
init[j]=0
init[i]=1
qc_test.initialize(init, range(nbits))
qc_last = qc_test + qc
result = qiskit.execute(qc_last, simulator).result()
statevector = result.get_statevector()
for j in range(2 ** nbits):
H[j,i]=statevector[j]
#Printing desired output
print(H)

In your case the output is :

[[1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 1.+0.j]
[0.+0.j 0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]
[0.+0.j 1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j]]

If you really need to avoid using any simulator you can try using the following : (I have not coded the permutation part, which is the complex part I would say)

import qiskit
from qiskit import QuantumCircuit, Aer
import numpy as np

#Any circuit you want to analyse
qc = QuantumCircuit(3)
qc.cx(1,0)
qc.cx(1,2)

for reg in qc.qregs:
nbits=len(reg)
H = np.identity(2**nbits, dtype=complex)
Htmp = np.identity(2**nbits, dtype=complex)
#This part will loop through all operators and construct the corresponding matrix
for instr, qargs, cargs in qc._data:
#This allows you to get the corresponding matrix of the local operator
loc_mat = instr.to_matrix()
comp = int(2**(np.log2(np.shape(H)[0]) - np.log2(np.shape(loc_mat)[0])))
#construct a complementary identity operator to tensor product with the local operator
Hkro=np.identity(comp, dtype=complex)
Htmp = np.kron(Hkro, loc_mat)
# should permut H first by rearranging the circuit by using the index of
# the used qregister, this part I did not have the courage to do
# (that's why people have coded it so well in unitary_simulator
# so why bother if you understand what is happening behind the scene)
for j in range(len(qargs)):
i_pos=qargs[j].index
for k in range(len(qargs)):
j_pos=qargs[k].index

H=np.matmul(H, Htmp)
#should back permut
print(qargs)
print(H)
• Thank you so much, i really appreciate it! Unfortunately I am trying to figure it out by not using any simulators, only algebra and normal matrix manipulation by using things such as tensor product – Jared Apr 3 at 19:39
• Hi @Jared, I added a piece of code to help do the whole thing without any simulators. I hope it helps. – PilouPili Apr 3 at 20:05
• Thank you very much, definitely helps! – Jared Apr 3 at 20:13

I found that if you multiply layer1 by layer2 (not layer2 @ layer1), the answer is consistent with the Qiskit one. So I tested through IBM experience and found that the state is reverse in Qiskit, which is a bit weird for me (that is if you apply a not gate to q0 in q0=0, q1=0, q2=0, a typical |000> state, it will become |001> but not |100>). Follow this idea you would find that the kronecker product in your first layer should be $$\ \mathbb{I}\ \otimes\$$SWAP-CNOT, exactly the reverse of your numpy computation. I'm also noob in Quantum Computing but it probably is the answer. (See Qiskit document https://qiskit.org/textbook/ch-gates/multiple-qubits-entangled-states.html to get more sense on how it works)

• Wow ok that makes a lot of sense now, was very confusing but that cleared it up. Thanks! – Jared Apr 4 at 12:29