# 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/… Apr 3, 2021 at 16:36
• @ryanhill1 Thank you. i tried this but could not get it to work unfortunately Apr 3, 2021 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 Apr 3, 2021 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$. Apr 3, 2021 at 20:33
• I just made these changes as seen in the above edit but I still do not get the correct result Apr 3, 2021 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]]) Apr 3, 2021 at 21:01
• Ok I tried it with the new swapcnot matrix but it still is not working for some reason Apr 3, 2021 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 Apr 3, 2021 at 19:39

Michele Amoretti's answer is essentially correct but there is a neater and less error-prone way to do it. There's no need to try to find the correct unitary matrix expressed with the correct qubit significance-ordering because the unitary matrices are all available in Qiskit's Gate classes as NumPy arrays. Taking them from here also saves you from manually typing the matrix out or copy-pasting it to have an explicit array in code (they can be printed anyway).

import numpy as np
import cmath
from qiskit.circuit.library import CXGate, IGate
from qiskit.aqua.utils import tensorproduct

cx_gate = CXGate().to_matrix()
id_gate = IGate().to_matrix()

layer1 = tensorproduct(id_gate, cx_gate)
layer2 = tensorproduct(cx_gate, id_gate)

qcirc_uni = np.matmul(layer2, layer1)

print(qcirc_uni.real)

array([[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.]])


Since this unitary matrix is real I used cmath's real method to display it more clearly.

The reason that order of the matrices in your tensor products in your first computation is incorrect is because of the qubit significance-ordering convention used in Qiskit: Why does Qiskit order its qubits the way it does? 1 Minute Qiskit. It's a really good idea to get on top of this.

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 Apr 3, 2021 at 19:39
• Hi @Jared, I added a piece of code to help do the whole thing without any simulators. I hope it helps. Apr 3, 2021 at 20:05
• Thank you very much, definitely helps! Apr 3, 2021 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! Apr 4, 2021 at 12:29