GateOps usage example
Simple gate-level operations
This notebook demonstrates some of the simple wrapping provided to examine the values of cascaded gate operations. Support is provided to manipulate the pybind11 wrapped matrices (DCM), or mapping them directly to numpy types, allowing use of the full scipy ecosystem.
%matplotlib notebook
from matplotlib import rc
rc('text', usetex=True)
import matplotlib.pyplot as plt
# Load the wrqpped simulator as before
from PyQNLPSimulator import PyQNLPSimulator as p
# Load the wrapped 2x2 matrix type; dcm => double complex matrix
from PyQNLPSimulator import DCMatrix as dcm
# Load the QNLP python packages
import QNLP as q
# Useful for examining various product types
from itertools import product
# Load numpy and matrix-based functions from scipy
import numpy as np
from scipy.linalg import expm, sinm, cosm
# Define a matrix using lexicographical ordering
ll=[0,1,1,0]
m = dcm(ll)
Using the matrix m we can perform a variety of oeprations using the underlying C++ wrapped methods.
Note: We do not restrict ourselves to these being unitary.
print(m*m)
print(m + m)
print(2*m - 5*m)
{{1.000+0.000 0.000+0.000 }{0.000+0.000 1.000+0.000 }{
{{0.000+0.000 2.000+0.000 }{2.000+0.000 0.000+0.000 }{
{{0.000+0.000 -3.000+0.000 }{-3.000+0.000 0.000+0.000 }{
Additionally, to avail of the numpy/scipy ecosystem, it makes sense to have conversions between the numpy datatypes and the above matrices. This allows using to perform a variety of useful operations.
# Create numpy variant
nm = m.as_numpy()
Below we demonstrate comparisons between both types, and the operations available to each. We can load the matrices directly using the ‘DCM’ type, from the PyQNLPSimulator class ‘GateOps’, or access them from an initialised simulator object.
from PyQNLPSimulator import GateOps
sim = p(3, False)
gops = GateOps(sim)
#Direct from simulator
X = sim.getGateX()
Y = sim.getGateY()
Z = sim.getGateZ()
H = sim.getGateH()
#GateOps access
Xgo = gops.X
#Numpy variant
Xnp = np.matrix(X.as_numpy())
Hnp_go = gops.Hnp
print(H * X)
print(Hnp_go * Xnp)
{{0.707+0.000 0.707+0.000 }{-0.707+0.000 0.707+0.000 }{
[[ 0.70710678+0.j 0.70710678+0.j]
[-0.70710678+0.j 0.70710678+0.j]]
Using the numpy variants, we can explicitly pass our matrices to the matrix functions loaded from scipy.linalg. In this case, we calculate $\exp(iH)$, where $H$ is the Hadamard gate.
expm(1j*H.as_numpy())
array([[ 5.40302306e-01+0.59500984j, 0.00000000e+00+0.59500984j],
[-5.55111512e-17+0.59500984j, 5.40302306e-01-0.59500984j]])
We can thus define rotation matrices from the above gate definitions using the $R_a(\theta) = \exp\left(-\frac{i\theta a}{2}\right)$ definition. The below listed methods are all available within the GateOps class
def RX(theta):
return np.matrix(expm(-0.5*1j*theta*X.as_numpy()))
def RY(theta):
return np.matrix(expm(-0.5*1j*theta*Y.as_numpy()))
def RZ(theta):
return np.matrix(expm(-0.5*1j*theta*Z.as_numpy()))
def RHam(theta, Ham):
return np.matrix(expm(-0.5*1j*theta*Ham.as_numpy()))
Looking at the application of this matrix on a simple vector can be performed as follows:
for i in range(10):
data = (RY(np.pi*(i/10))*np.array([[1],[0]]))
print(data)
[[1.+0.j]
[0.+0.j]]
[[0.98768834+0.j]
[0.15643447+0.j]]
[[0.95105652+0.j]
[0.30901699+0.j]]
[[0.89100652+0.j]
[0.4539905 +0.j]]
[[0.80901699+0.j]
[0.58778525+0.j]]
[[0.70710678+0.j]
[0.70710678+0.j]]
[[0.58778525+0.j]
[0.80901699+0.j]]
[[0.4539905 +0.j]
[0.89100652+0.j]]
[[0.30901699+0.j]
[0.95105652+0.j]]
[[0.15643447+0.j]
[0.98768834+0.j]]
We can visualise this rotation as follows:
arr = np.array([[1],[0]])
for i in range(10):
arr = RY(np.pi/10)*arr
arr
matrix([[-1.94289029e-16+0.j],
[ 1.00000000e+00+0.j]])
bin(296)
'0b100101000'
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm, colors
fig = plt.figure(figsize=plt.figaspect(1))
ax = fig.add_subplot(111, projection='3d')
# Bloch sphere
r = 0.99
phi, theta = np.mgrid[0.0:np.pi:100j, 0.0:2.0*np.pi:100j]
x0 = r*np.sin(phi)*np.cos(theta)
y0 = r*np.sin(phi)*np.sin(theta)
z0 = r*np.cos(phi)
ax.plot_surface(x0, y0, z0, color='violet', linewidth=1.0, alpha=0.55)
data = np.array([[0],[1]])
arrow0 = ax.quiver(0,0,0, np.real(data[0]), 0, np.real(data[1]), length=1.0, color='black')
for i in range(10):
data = (RY(np.pi*(i/10))*np.array([[1],[0]]))
arrow = ax.quiver(0,0,0, np.real(data[0]), 0, np.real(data[1]), length=1.0)
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
ax.set_zlim([-1,1])
#ax.set_aspect("equal")
plt.tight_layout()
ax.grid(False)
plt.pause(0.05)
plt.draw()
To examine the product of these gates, as would be offered by the quantum circuit model, we can create tensor products as:
We can, as demonstrated previously, raise these to exponents. Care should be taken though, as Kronecker products grow the memory footprint greatly.
s1 = np.kron(gops.Inp, gops.Znp)
s2 = 0.5*(gops.Znp + gops.Inp)
s3 = np.kron(Xnp, np.kron(s2,s1))
s3
matrix([[ 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, -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, 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, -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, 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],
[ 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, 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, 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,
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, 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, 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, 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, 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,
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, 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, 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, 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, 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],
[ 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, 0.+0.j, -0.+0.j, 0.+0.j, -0.+0.j,
0.+0.j, -0.+0.j]])
where the matrix s3 is given by $\mathbf{I}\otimes(\vert 0 \rangle\langle 0 \vert )\otimes(\mathbf{I}\otimes\sigma_z)$.
aa = expm(-0.5*1j*np.pi*s3)
fig = plt.figure()
plt.subplot(121)
plt.imshow(abs(aa)**2)
plt.colorbar()
plt.subplot(122)
plt.imshow(np.angle(aa))
plt.colorbar()
plt.show()
Next, we can demonstrate the usage of the UnitaryFinder class. As an example, say we are given an arbitrary 2x2 unitary matrix, and a set of gates implemenetd on our simulator/hardware platform. We cannot directly implement the given gate, but would like to decompose it (e.g. compile it) to the available gate-set. Using the matrix form given by createS(i), we attempt to find the optimal gates to create this matrix.
def createS(i):
a = np.sqrt((i-1)/i)
b = 1/np.sqrt(i)
return dcm([a,b,-b,a])
np.matrix(createS(4).as_numpy()).flatten().tolist()[0]
[(0.8660254037844386+0j), (0.5+0j), (-0.5+0j), (0.8660254037844386+0j)]
from PyQNLPSimulator import UnitaryFinder
# Creating the UnitaryFinder object, we give the S matrix as an argument
uf = UnitaryFinder(sim, np.matrix(createS(4).as_numpy()).flatten().tolist()[0])
Next, we generate all possible combinations of gates in the given gate-set (taken from GateOps) up to a given depth.
gate_op_combo = uf.genOps(uf.gates, depth=3)
print(gate_op_combo[0:4])
[('X', 'Y', 'X'), ('X', 'Y', 'Z'), ('X', 'Y', 'H'), ('X', 'Y', 'RX')]
Lastly, we iterate through this gate combination list, and aim to find the gate combinations that best match our target matrix, to within a given tolerance.
for i in gate_op_combo:
res, vals = uf.findOps(i)
if res:
print(i, vals)
('X', 'Z', 'RY') [2.09207924]
('X', 'H', 'RY') [0.52239173]
('X', 'RY', 'X') [1.04575906]
('X', 'RY', 'Z') [-2.09207924]
('X', 'RY', 'H') [-0.52239173]
('Y', 'RY', 'Y') [-1.04575906]
('Z', 'H', 'RY') [-2.61570235]
('Z', 'RY', 'Z') [1.04575906]
('Z', 'RY', 'H') [2.61570235]
('Z', 'RY', 'RZPh') [3.13371396 1.04713375]
('Z', 'RZPh', 'RY') [-1.04713398 3.13371395]
('H', 'X', 'RY') [-2.61570235]
('H', 'Z', 'RY') [0.52239173]
('H', 'RY', 'X') [2.61570235]
('H', 'RY', 'Z') [-0.52239173]
('H', 'RY', 'H') [1.04575906]
('H', 'RY', 'RZPh') [ 3.13596229 -0.52362299]
('H', 'RZPh', 'RY') [0.52362323 3.135996 ]
('RX', 'RY', 'RX') [-5.18605916e-04 -1.04613486e+00 -5.18605924e-04]
('RX', 'RY', 'RZ') [-1.85780375e-09 -1.04561088e+00 -2.24742458e-09]
('RX', 'RY', 'RZPh') [-1.38746315e-03 -1.04758429e+00 9.40171898e-04]
('RX', 'RZ', 'RY') [-1.04566595e+00 4.46520933e-04 4.46537155e-04]
('RX', 'RZPh', 'RY') [-1.04711238e+00 -8.72097633e-04 -1.39162598e-03]
('RY', 'X', 'Z') [2.09207924]
('RY', 'X', 'H') [0.52239173]
('RY', 'Z', 'H') [-2.61570235]
('RY', 'Z', 'RZPh') [ 3.13371396 -1.04713369]
('RY', 'H', 'X') [-2.61570235]
('RY', 'H', 'Z') [0.52239173]
('RY', 'H', 'RZPh') [3.13596229 0.52362299]
('RY', 'RX', 'RY') [-5.22950526e-01 -9.62041915e-05 -5.22950524e-01]
('RY', 'RX', 'RZ') [ 4.46528063e-04 4.46538670e-04 -1.04566596e+00]
('RY', 'RX', 'RZPh') [ 8.53317080e-04 -8.54122439e-05 -1.04720800e+00]
('RY', 'RZ', 'RX') [-6.64169503e-05 -7.95999326e-05 -1.04765239e+00]
('RY', 'RZ', 'RY') [-5.22950591e-01 -9.61878893e-05 -5.22950592e-01]
('RY', 'RZ', 'RZPh') [-0.00109832 0.00257097 -1.04706923]
('RY', 'RZPh', 'Z') [ 3.13371396 -1.04713369]
('RY', 'RZPh', 'H') [ 3.13535753 -2.61797294]
('RY', 'RZPh', 'RX') [ 5.13018218e-04 -5.40212364e-04 -1.04685560e+00]
('RY', 'RZPh', 'RY') [-0.52363945 0.00106341 -0.52363945]
('RY', 'RZPh', 'RZ') [ 3.71069550e-04 -1.19844950e-03 -1.04711682e+00]
('RZ', 'RX', 'RY') [-1.04765239e+00 -6.64180325e-05 -7.96009587e-05]
('RZ', 'RY', 'RX') [ 0.00393974 -1.04720165 0.00341443]
('RZ', 'RY', 'RZ') [-5.18605855e-04 -1.04613486e+00 -5.18605828e-04]
('RZ', 'RY', 'RZPh') [ 5.36333318e-04 -1.04571135e+00 4.26673895e-04]
('RZ', 'RZPh', 'RY') [-1.04612097e+00 5.12574535e-04 -2.50861681e-04]
('RZPh', 'Z', 'RY') [-1.04713398 3.13371395]
('RZPh', 'H', 'RY') [-2.61763872 3.13477955]
('RZPh', 'RX', 'RY') [-1.04559432e+00 8.78686271e-06 7.95029099e-05]
('RZPh', 'RY', 'Z') [1.04713398 3.13371395]
('RZPh', 'RY', 'H') [2.61763874 3.13477954]
('RZPh', 'RY', 'RX') [-1.19556857e-04 -1.04509814e+00 3.74598119e-05]
('RZPh', 'RY', 'RZ') [-8.54517929e-06 -1.04564993e+00 5.48402952e-06]
('RZPh', 'RY', 'RZPh') [-3.66001552e-04 -1.04537300e+00 -3.66001226e-04]
('RZPh', 'RZ', 'RY') [-1.04722635e+00 -5.07168483e-04 -4.17892558e-04]
The values returned above list the gate combinations (applied right to left) and the parameters (angles for rotation matrices) that best-approximated the given unitary.
For our example above of S4, some simple arithmetic shows that $S4 = \sigma_x RY(\cos^{-1}(1 / 2))\sigma_x$. A further investigation gives us a relationship for $S_i = \sigma_x RY(\cos^{-1}((i-2) / 2)))\sigma_x$ as one of many possible combinations, which can be further reduced to $S_i = RY(-\cos^{-1}( (i-2) / 2 )))$.