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

Question

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.]]

Answer 1

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.]]

Answer 2

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} $

Answer 3

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.

Answer 4

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)

Answer 5

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)