QComp - Coding Session #3: Trotterization¶
SOLUTIONS¶
This is the third part of a series of lab sessions to implement and analyze the HHL algorithm.
- In Coding Session #1, we implemented QPE
- In Coding Session #2, we implemented HHL
- In this coding session, we will analyze different strategy for implementing the gate $e^{iA}$, and its powers.
This coding session is for practice, there is nothing to hand out.
First, make sure libraries are installed (you can uncomment those if needed)
# ! python -m pip install matplotlib sympy
# ! python -m pip install qiskit qiskit-aer
Now, some libraries to load (nothing to modify here)
import itertools
import random
import numpy as np
import matplotlib.pyplot as plt
from scipy import linalg
from math import pi, gcd
from qiskit import *
from qiskit.circuit import *
from qiskit.circuit.library import *
from qiskit.compiler import transpile
from qiskit.quantum_info.operators import Operator
from qiskit_aer import AerSimulator, StatevectorSimulator, UnitarySimulator
from scipy import optimize
%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
Small library for pretty-printing (nothing to do)
def nat2bl(pad,n):
if n == 0: r = [0 for i in range(pad)]
elif n % 2 == 1: r = nat2bl(pad-1,(n-1)//2); r.append(1)
else: r = nat2bl(pad-1,n//2); r.append(0)
return r
def bl2nat(s):
if len(s) == 0: return 0
else: a = s.pop(); return (a + 2*bl2nat(s))
def bl2bs(l):
if len(l) == 0: return ""
else: a = l.pop(); return (bl2bs(l) + str(a))
def nat2bs(pad,i): return bl2bs(nat2bl(pad,i))
def bs2bl(s):
l = []
for i in range(len(s)): l.append(int(s[i]))
return l
def bs2nat(s): return bl2nat(bs2bl(s))
def processOneState(st): # Length = power of 2
s = list(st)
if len(s) == 2: return {'0' : s[0], '1' : s[1]}
else:
a0 = processOneState(s[:len(s)//2])
a1 = processOneState(s[len(s)//2:])
r = {}
for k in a0: r['0' + k] = a0[k]
for k in a1: r['1' + k] = a1[k]
return r
def printOneState(d, keys=None): # get a dict as per processStates output
isThereMore = True
if keys == None:
keys = d.keys()
isThereMore = False
for k in sorted(keys):
if k not in d:
print(f"Error: {k} not valid. Maybe the size of the bitstring is invalid?")
isThereMore = False
break
im = d[k].imag
re = d[k].real
if abs(im) >= 0.001 or abs(re) >= 0.001:
print("% .5f + % .5fj |%s>" % (re,im,k))
if isThereMore: print(" ... (there might be more hidden terms)")
def printFinalRes(result, keys=None):
printOneState(processOneState(list(np.asarray(result))), keys)
def runStateVector(qc, keys=None):
simulator = StatevectorSimulator(statevector_parallel_threshold=6)
job = simulator.run(qc.decompose(reps=6), memory=True)
job_result = job.result()
result = job_result.results[0].to_dict()['data']['statevector']
printFinalRes(result, keys)
def runStateVectorSeveralTimes(qc, howmany):
qc.save_statevector(label = 'collect', pershot = True)
simulator = StatevectorSimulator()
job = simulator.run(qc.decompose(reps=6), memory=True, shots=howmany)
result = job.result()
memory = result.data(0)['memory']
collect = result.data(0)['collect']
r = {}
for i in range(len(collect)):
r[str(collect[i])] = (0, collect[i])
for i in range(len(collect)):
n, v = r[str(collect[i])]
r[str(collect[i])] = (n+1, v)
for k in r:
i, v = r[k]
print(f"With {i} occurences:")
printFinalRes(v)
def plotDistrib(d):
sorted_items = sorted(d.items())
keys = [k for k, _ in sorted_items]
values = [v for _, v in sorted_items]
plt.figure()
plt.bar(keys, values)
plt.xticks(rotation=90)
plt.show()
def getSample(qc,howmany):
simulator = AerSimulator()
job = simulator.run(qc.decompose(reps=6), shots=howmany)
res = dict(job.result().get_counts(qc))
return res
def plotSample(qc,howmany):
d = getSample(qc,howmany)
ld = len(list(d.keys())[0])
for i in range(2 ** ld):
s = nat2bs(ld,i)
if s not in d: d[s] = 0
plotDistrib(d)
def flipEndian(matrix):
lm = len(matrix)
n = int(np.log2(lm))
if lm != 2 ** n:
raise Exception("dimension error")
res = [[0 for i in range(lm)] for j in range(lm)]
for i in range(lm):
ii = bs2nat(nat2bs(n,i)[::-1])
for j in range(lm):
jj = bs2nat(nat2bs(n,j)[::-1])
res[i][j] = matrix[ii][jj]
return np.array(res)
q = QuantumRegister(2, name="q")
qc = QuantumCircuit(q)
qc.cx(q[0],q[1])
matrix_cnot = np.array(Operator.from_circuit(qc))
matrix_cnot
array([[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, 1.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j]])
You might notice that we asked for a CNOT gate... Qiskit is by default placing the least significant bit on the left. I wrote for you a small wrapper function flipping its position, when needed.
flipEndian(matrix_cnot)
array([[1.+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, 1.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j]])
0.2 - Estimating the distance between two matrices.¶
Let us add a very small perturbation to the circuit qc.
qc.rx(0.1, q[0])
matrix_cnot_perturb = np.array(Operator.from_circuit(qc))
matrix_cnot_perturb
array([[0.99875026+0.j , 0. +0.j ,
0. +0.j , 0. -0.04997917j],
[0. -0.04997917j, 0. +0.j ,
0. +0.j , 0.99875026+0.j ],
[0. +0.j , 0. -0.04997917j,
0.99875026+0.j , 0. +0.j ],
[0. +0.j , 0.99875026+0.j ,
0. -0.04997917j, 0. +0.j ]])
The library scipy contains a method to estimate the spectral norm of a
matrix. We can then estimate the error between matrix_cnot and
matrix_cnot_perturb as follows.
linalg.interpolative.estimate_spectral_norm(matrix_cnot - matrix_cnot_perturb)
0.049994791829424665
This matrix norm is the one that we shall be using in the rest of this lab session.
0.3 - Circuit decomposition¶
QisKit comes with an automated process to decompose circuits into elementary gates. It is using a completely generic procedure, akin to what can be found in Section 4.3 of the lecture notes.
We show here how to use it to decompose the matrix that we call $A_2$:
A2 = np.array([[1,2,0,0],
[2,1,2,0],
[0,2,1,2],
[0,0,2,1]])
t = 9/16
U = UnitaryGate(Operator(linalg.expm(1j*t*A2)), label="U")
q = QuantumRegister(2, name="q")
qc = QuantumCircuit(q)
qc.append(U,q)
# We ask for a decomposition with CNOTs, Rx and Rz gates.
qc_elementary = transpile(qc, basis_gates=['cx', 'rx', 'rz'])
print(qc_elementary.draw())
d = dict(qc_elementary.count_ops())
print("each gate appears as follows:", d)
r = 0
for i in d:
r += d[i]
print("total number of gates:", r)
global phase: 5.2749
┌─────────┐ ┌─────────┐ ┌─────────────┐ ┌─────────────┐»
q_0: ┤ Rz(π/2) ├───┤ Rx(π/2) ├──┤ Rz(-1.2236) ├──■──┤ Rx(-2.2626) ├»
├─────────┴┐┌─┴─────────┴─┐└─────────────┘┌─┴─┐└┬────────────┤»
q_1: ┤ Rz(-π/2) ├┤ Rx(-1.7229) ├───────────────┤ X ├─┤ Rz(2.0329) ├»
└──────────┘└─────────────┘ └───┘ └────────────┘»
« ┌────────────┐ ┌─────────┐ ┌──────────┐
«q_0: ┤ Rz(1.2236) ├────────────────■────┤ Rx(π/2) ├──┤ Rz(-π/2) ├
« ├────────────┤┌────────────┐┌─┴─┐┌─┴─────────┴─┐├─────────┬┘
«q_1: ┤ Rx(1.6561) ├┤ Rz(-3.004) ├┤ X ├┤ Rx(0.15207) ├┤ Rz(π/2) ├─
« └────────────┘└────────────┘└───┘└─────────────┘└─────────┘
each gate appears as follows: {'rz': 8, 'rx': 6, 'cx': 2}
total number of gates: 16
B = np.array([[1,3],[3,-1]])
t = 1/3
exact = linalg.expm(1j * t * B)
d = {}
def trotterB(p): # p is the number of iterations
q = QuantumRegister(1, name="q")
qc = QuantumCircuit(q)
for i in range(p):
qc.rx(-6*t/p, q)
qc.rz(-2*t/p, q)
return qc
approx = np.array(Operator.from_circuit(trotterB(100)))
error = linalg.interpolative.estimate_spectral_norm(exact - approx)
print("error:", error)
error: 0.0027494676292157167
1.2 - Behavior of the Error w.r.t. Iterations Count¶
In class, we discussed the fact that the error between $e^{itM}$ and its Trotter implementation $U^M_{t,n}$ using $n$ iterations hs the following asymptotic behavior:
$$ |\!| e^{itM} - U^{M}_{t,n} |\!| \simeq O\left(\frac{t^2}{n}\right) $$
In the next cell, plot the function $n \mapsto 1/f(n)$, where $f(n)$ is $$ f(n) = |\!| e^{i\frac1{3}B} - U^B_{\frac1{3},n} |\!| $$
You can use $n=10, 20, \ldots 200$.
How many iterations would be needed to reach an error of $10^{-8}$ ?
Hint The function np.polyfit([x1,x2,x3],[y1,y2,y3],1) gives the coefficients of the best linear approximation going through the points (x1,y1), (x2,y2) and (x3,y3). It works for any number of points.
exact = linalg.expm(1j * t * B)
d = {}
for p in range(10,200,10):
approx = np.array(Operator.from_circuit(trotterB(p)))
d[p] = linalg.interpolative.estimate_spectral_norm(exact - approx)
sorted_items = sorted(d.items())
keys = [k for k, _ in sorted_items]
values = [1/v for _, v in sorted_items]
[a,b] = np.polyfit(keys,values,1)
print(f"linear approximation : {a}*x + {b}")
print(f"We get an error of 10^-4 with {np.ceil(10000/a)} iterations")
print(f"We get an error of 10^-8 with {np.ceil((10**8)/a)} iterations")
plt.plot(keys, values)
plt.show()
linear approximation : 3.6371607404495445*x + -0.011373713282339836 We get an error of 10^-4 with 2750.0 iterations We get an error of 10^-8 with 27493974.0 iterations
2 - More involved matrices¶
2.1 - The easy case¶
Consider the matrix $C$, of size $4\times4$:
$$ C = \left( \begin{matrix} 1 & 2 & 0 & 2\\ 2 & 1 & 2 & 0\\ 0 & 2 & 1 & 2\\ 2 & 0 & 2 & 1 \end{matrix} \right) $$
TODO
- Generate a decomposition as a linear combination of Pauli strings.
- Using this linear combination, write a decomposition for $e^{i\frac12 C}$ as an exact circuit.
- Verify that it work by computing the error with the exact matrix.
Hint
Global phase can be added to a circuit by simply doing qc.global_phase = [something]
SOLUTION
$$ C = I\otimes I + 2\cdot I\otimes X + 2\cdot X\otimes X $$
They commute, so
$$ e^{i\frac12(I\otimes I + 2\cdot I\otimes X + 2\cdot X\otimes X)} = e^{i\frac12I\otimes I}e^{2i\frac12I\otimes X}e^{2i\frac12X\otimes X} $$
C = np.array([[1,2,0,2],
[2,1,2,0],
[0,2,1,2],
[2,0,2,1]])
q = QuantumRegister(2, name="q")
qc = QuantumCircuit(q)
t = 1/2
qc.global_phase = t
qc.rx(-4*t, q[0])
qc.h(q)
qc.cx(q[0],q[1])
qc.rz(-4*t, q[1])
qc.cx(q[0],q[1])
qc.h(q)
print(qc.draw())
approx = np.array(Operator.from_circuit(qc))
exact = linalg.expm(1j * t * C)
diff = linalg.interpolative.estimate_spectral_norm(exact - approx)
print("Error is", diff)
global phase: 0.5
┌────────┐┌───┐ ┌───┐
q_0: ┤ Rx(-2) ├┤ H ├──■──────────────■──┤ H ├
└─┬───┬──┘└───┘┌─┴─┐┌────────┐┌─┴─┐├───┤
q_1: ──┤ H ├────────┤ X ├┤ Rz(-2) ├┤ X ├┤ H ├
└───┘ └───┘└────────┘└───┘└───┘
Error is 5.96282874802339e-16
2.2 - A small change¶
Let us remove the extremal $2$ to reach the matrix
$$ A_2 = \left( \begin{matrix} 1 & 2 & 0 & 0\\ 2 & 1 & 2 & 0\\ 0 & 2 & 1 & 2\\ 0 & 0 & 2 & 1 \end{matrix} \right) $$
TODO
- Write $A_2$ as a linear decomposition of Pauli strings. You will need the fact that
$$ Y\otimes Y = \left( \begin{matrix} 0 & 0 & 0 &-1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \end{matrix} \right) $$
- Use Trotterization to approximate $e^{i\frac12A_2}$
- How many iteration steps are needed for reaching a precision of $10^{-3}$?
SOLUTION
$$ A_2 = I\otimes I + 2\cdot I \otimes X + Y\otimes Y + X \otimes X $$
A2 = np.array([[1,2,0,0],
[2,1,2,0],
[0,2,1,2],
[0,0,2,1]])
t = 1/2
def IX(qc,q,N):
qc.rx(-4*t/N, q[0])
def XX(qc,q,N):
qc.h(q)
qc.cx(q[0],q[1])
qc.rz(-2*t/N, q[1])
qc.cx(q[0],q[1])
qc.h(q)
def YY(qc,q,N):
qc.sdg(q)
qc.h(q)
qc.cx(q[0],q[1])
qc.rz(-2*t/N, q[1])
qc.cx(q[0],q[1])
qc.h(q)
qc.s(q)
def trotterA2(N):
q = QuantumRegister(2, name="q")
qc = QuantumCircuit(q)
qc.global_phase = t
for _ in range(N):
IX(qc,q,N)
XX(qc,q,N)
YY(qc,q,N)
return qc
qc = trotterA2(403)
approx = np.array(Operator.from_circuit(qc))
exact = linalg.expm(1j * t * A2)
diff = linalg.interpolative.estimate_spectral_norm(exact - approx)
print("Error is", diff)
Error is 0.000997899523928073
3 - Going larger¶
We shall now consider a family of matrices of the shape of $A_2$: the matrices $A_n$ are defined inductively following the following pattern.
$A_1$ of size $2\times2$: $$ \left( \begin{matrix} 1 & 2\\ 2 & 1 \end{matrix} \right) $$
$A_2$ of size $4\times4$:
$$ \left( \begin{matrix} 1 & 2 & 0 & 0\\ 2 & 1 & 2 & 0\\ 0 & 2 & 1 & 2\\ 0 & 0 & 2 & 1 \end{matrix} \right) $$
$A_3$ is of size $8\times8$:
$$ \left( \begin{matrix} 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0\\ 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0\\ 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0\\ 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0\\ 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2\\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 \end{matrix} \right) $$
$A_4$ of size $16\times16$:
$$ \left( \begin{matrix} 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1 & 2\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 1\\ \end{matrix} \right) $$
3.1 - Coding the family $A_n$¶
TODO Fill in the following code so that A(n) is the matrix $A_n$.
def An(n):
a = 1
b = 2
N = 2 ** n
A = [[0 for i in range(N)] for j in range(N)]
for i in range(N):
A[i][i] = a
if i > 0:
A[i-1][i] = b
A[i][i-1] = b
return np.array(A)
print(An(2))
print(An(3))
print(An(4))
[[1 2 0 0] [2 1 2 0] [0 2 1 2] [0 0 2 1]] [[1 2 0 0 0 0 0 0] [2 1 2 0 0 0 0 0] [0 2 1 2 0 0 0 0] [0 0 2 1 2 0 0 0] [0 0 0 2 1 2 0 0] [0 0 0 0 2 1 2 0] [0 0 0 0 0 2 1 2] [0 0 0 0 0 0 2 1]] [[1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0] [2 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0] [0 2 1 2 0 0 0 0 0 0 0 0 0 0 0 0] [0 0 2 1 2 0 0 0 0 0 0 0 0 0 0 0] [0 0 0 2 1 2 0 0 0 0 0 0 0 0 0 0] [0 0 0 0 2 1 2 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 2 1 2 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 2 1 2 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 2 1 2 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 2 1 2 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 2 1 2 0 0 0 0] [0 0 0 0 0 0 0 0 0 0 2 1 2 0 0 0] [0 0 0 0 0 0 0 0 0 0 0 2 1 2 0 0] [0 0 0 0 0 0 0 0 0 0 0 0 2 1 2 0] [0 0 0 0 0 0 0 0 0 0 0 0 0 2 1 2] [0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 1]]
3.2 - Pauli decompositions of $A_n$¶
We want to generate a decomposition of $A_n$ as a sum of Pauli strings. We could look for an algebraic formulation by deriving an inductive formula, but I propose instead to reformulate the problem into a linear system of equations. If this is very crude and could be done better, it somehow works for matrices over a handful of qubits, enough for this lab session.
The following code shows one way to do it (nothing to do until the TODO below).
pauli = ["i", "x", "y", "z"]
pmat = {"i" : np.matrix([[1,0],[0,1]]),
"x" : np.matrix([[0,1],[1,0]]),
"y" : np.matrix([[0,-1j],[1j,0]]),
"z" : np.matrix([[1,0],[0,-1]])}
def tens(a,b):
na = len(a)
nb = len(b)
c = [[0 for i in range(na*nb)] for j in range(na*nb)]
for ia in range(na):
for ja in range(na):
for ib in range(nb):
for jb in range(nb):
c[ia*nb + ib][ja*nb + jb] = a[ia,ja]*b[ib,jb]
return np.matrix(c)
def ntens(a):
res = np.matrix([[1]])
for k in range(len(a)): res = tens(res, a[k])
return res
def paulistring_m(a):
return ntens([pmat[p] for p in a])
print("Y tensor Y")
print(paulistring_m(["y","y"]))
print("Z tensor X")
print(paulistring_m(["z","x"]))
print("I tensor Z tensor X")
print(paulistring_m(["i", "z","x"]))
Y tensor Y [[ 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 1.-0.j 0.+0.j 0.-0.j] [-1.+0.j 0.+0.j 0.+0.j 0.+0.j]] Z tensor X [[ 0 1 0 0] [ 1 0 0 0] [ 0 0 0 -1] [ 0 0 -1 0]] I tensor Z tensor X [[ 0 1 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 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 0 0 0 0 -1 0]]
def paulistrings(n):
res = list(itertools.product(pauli, repeat=n))
return [list(item) for item in res]
print("All Pauli strings of size 3:")
print(paulistrings(3))
All Pauli strings of size 3: [['i', 'i', 'i'], ['i', 'i', 'x'], ['i', 'i', 'y'], ['i', 'i', 'z'], ['i', 'x', 'i'], ['i', 'x', 'x'], ['i', 'x', 'y'], ['i', 'x', 'z'], ['i', 'y', 'i'], ['i', 'y', 'x'], ['i', 'y', 'y'], ['i', 'y', 'z'], ['i', 'z', 'i'], ['i', 'z', 'x'], ['i', 'z', 'y'], ['i', 'z', 'z'], ['x', 'i', 'i'], ['x', 'i', 'x'], ['x', 'i', 'y'], ['x', 'i', 'z'], ['x', 'x', 'i'], ['x', 'x', 'x'], ['x', 'x', 'y'], ['x', 'x', 'z'], ['x', 'y', 'i'], ['x', 'y', 'x'], ['x', 'y', 'y'], ['x', 'y', 'z'], ['x', 'z', 'i'], ['x', 'z', 'x'], ['x', 'z', 'y'], ['x', 'z', 'z'], ['y', 'i', 'i'], ['y', 'i', 'x'], ['y', 'i', 'y'], ['y', 'i', 'z'], ['y', 'x', 'i'], ['y', 'x', 'x'], ['y', 'x', 'y'], ['y', 'x', 'z'], ['y', 'y', 'i'], ['y', 'y', 'x'], ['y', 'y', 'y'], ['y', 'y', 'z'], ['y', 'z', 'i'], ['y', 'z', 'x'], ['y', 'z', 'y'], ['y', 'z', 'z'], ['z', 'i', 'i'], ['z', 'i', 'x'], ['z', 'i', 'y'], ['z', 'i', 'z'], ['z', 'x', 'i'], ['z', 'x', 'x'], ['z', 'x', 'y'], ['z', 'x', 'z'], ['z', 'y', 'i'], ['z', 'y', 'x'], ['z', 'y', 'y'], ['z', 'y', 'z'], ['z', 'z', 'i'], ['z', 'z', 'x'], ['z', 'z', 'y'], ['z', 'z', 'z']]
import sympy as sym
#pauli = ["i", "x", "y"] # for An, we don't need Z... Add it back if you work with another matrix shape.
def decomposeAsPauliStrings(A):
n = int(np.log2(len(A)))
ps = paulistrings(n)
pvar = [sym.Symbol("".join(ps[k])) for k in range(len(ps))]
res = 0
for k in range(len(ps)):
res = res + pvar[k] * paulistring_m(ps[k])
eqs = [item for sublist in (res - A).tolist() for item in sublist]
solx = sym.solve(eqs, pvar)
sol = {}
for k in solx:
if solx[k] != 0.0:
sol[k] = solx[k]
return sol
# Let us decompose C
print(decomposeAsPauliStrings(C))
{ii: 1.00000000000000, ix: 2.00000000000000, xx: 2.00000000000000}
TODO
- For $n=2,3,4,5$, count how many Pauli strings are in the decomposition of $A_n$.
- What is your guess on the behavior of the number of terms in the decomposition for $n$, in general?
- Can you give a bound on the number of gates in the Trotterization of $A_n$, parameterized by $n$ and
p, the number of steps in the iteration? - If you had enough time and memory and asked Qiskit to give you an exact decomposition for $A_n$ in elementary gates, can you estimate the size of the circuit as a function of $n$?
- Consider 1000 iterations: For which $n$ would Trotterization win over exact synthesis?
for k in range(2,5):
print(f"Decomposition of A{k}:")
d = decomposeAsPauliStrings(An(k))
print(f"There are {len(d)} terms, with coefficients as follows:")
print(d)
print()
Decomposition of A2:
There are 4 terms, with coefficients as follows:
{ii: 1.00000000000000, ix: 2.00000000000000, xx: 1.00000000000000, yy: 1.00000000000000}
Decomposition of A3:
There are 8 terms, with coefficients as follows:
{iii: 1.00000000000000, iix: 2.00000000000000, ixx: 1.00000000000000, iyy: 1.00000000000000, xxx: 0.500000000000000, xyy: -0.500000000000000, yxy: 0.500000000000000, yyx: 0.500000000000000}
Decomposition of A4:
There are 16 terms, with coefficients as follows:
{iiii: 1.00000000000000, iiix: 2.00000000000000, iixx: 1.00000000000000, iiyy: 1.00000000000000, ixxx: 0.500000000000000, ixyy: -0.500000000000000, iyxy: 0.500000000000000, iyyx: 0.500000000000000, xxxx: 0.250000000000000, xxyy: -0.250000000000000, xyxy: -0.250000000000000, xyyx: -0.250000000000000, yxxy: 0.250000000000000, yxyx: 0.250000000000000, yyxx: 0.250000000000000, yyyy: -0.250000000000000}
The number of Pauli strings in the decomposition of $A_n$ is then literally $2^n$.
A Pauli string of size $n$ yields a subcircuit containing $2n$ CNOT gates, and at most $4n$ one-qubit gates. The circuit size of one iteration is therefore bounded by $6n2^n$.
On the other hand, Qiskit generates circuits of size:
#
keys = []
values = []
for n in range(2,7):
U = UnitaryGate(Operator(linalg.expm(1j*t*An(n))), label="U")
q = QuantumRegister(n, name="q")
qc = QuantumCircuit(q)
qc.append(U,q)
qc_elementary = transpile(qc, basis_gates=['cx', 'rx', 'rz'])
d = dict(qc_elementary.count_ops())
r = 0
for i in d:
r += d[i]
print(f"total number of gates for A{n}: {r}")
keys.append(n)
values.append(np.log2(r))
[a,b] = np.polyfit(keys,values,1)
print(f"linear approximation : {a}*x + {b}")
print(f"We get an error of 10^-4 with {np.ceil(10000/a)} iterations")
plt.plot(keys, values)
plt.show()
total number of gates for A2: 16 total number of gates for A3: 73 total number of gates for A4: 413
total number of gates for A5: 1724
total number of gates for A6: 7160 linear approximation : 2.217320724451232*x + -0.3818608054968873 We get an error of 10^-4 with 4510.0 iterations
The number of gates in Qiskit decomposition for $A(n)$ is therefore a bit above $4^n$.
For a fixed number $p$ of iteration, $4^n$ eventually becomes larger than $6np2^n$.
print("Example for p=1000")
for n in range(2,20):
tr = 6000*n*(2**n)
qi = 4 ** n
if tr < qi:
print(f"For A({n}), Trotter wins with {tr} gates, while Qiskit {qi} gates")
else:
print(f"For A({n}), Trotter lags behind with {tr} gates, while Qiskit {qi} gates")
Example for p=1000 For A(2), Trotter lags behind with 48000 gates, while Qiskit 16 gates For A(3), Trotter lags behind with 144000 gates, while Qiskit 64 gates For A(4), Trotter lags behind with 384000 gates, while Qiskit 256 gates For A(5), Trotter lags behind with 960000 gates, while Qiskit 1024 gates For A(6), Trotter lags behind with 2304000 gates, while Qiskit 4096 gates For A(7), Trotter lags behind with 5376000 gates, while Qiskit 16384 gates For A(8), Trotter lags behind with 12288000 gates, while Qiskit 65536 gates For A(9), Trotter lags behind with 27648000 gates, while Qiskit 262144 gates For A(10), Trotter lags behind with 61440000 gates, while Qiskit 1048576 gates For A(11), Trotter lags behind with 135168000 gates, while Qiskit 4194304 gates For A(12), Trotter lags behind with 294912000 gates, while Qiskit 16777216 gates For A(13), Trotter lags behind with 638976000 gates, while Qiskit 67108864 gates For A(14), Trotter lags behind with 1376256000 gates, while Qiskit 268435456 gates For A(15), Trotter lags behind with 2949120000 gates, while Qiskit 1073741824 gates For A(16), Trotter lags behind with 6291456000 gates, while Qiskit 4294967296 gates For A(17), Trotter wins with 13369344000 gates, while Qiskit 17179869184 gates For A(18), Trotter wins with 28311552000 gates, while Qiskit 68719476736 gates For A(19), Trotter wins with 59768832000 gates, while Qiskit 274877906944 gates