QComp - Coding Session #4: Oracle Synthesis¶
SOLUTIONS¶
This is the fourth 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 Coding Session #3, we saw a possible strategy to implement $e^{itA}$
- In this coding session, we will discuss how to implement an oracle such as the inversion oracle of HHL.
This session is only for practice.
First, make sure libraries are installed (you can uncomment those if needed)
# ! python -m pip install matplotlib sympy scipy
# ! python -m pip install qiskit qiskit-aer
Now, some libraries to load (nothing to modify here)
# Set-up for long lines
from IPython import get_ipython
from IPython.display import display, HTML
if get_ipython().__class__.__name__ == 'ZMQInteractiveShell':
print("Running inside Jupyter")
display(HTML("<style>pre { white-space: pre !important; }</style>"))
# Bunch of imports
import itertools, math, random
import sympy, scipy
import numpy as np
import matplotlib.pyplot as plt
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 linalg, optimize
%matplotlib inline
%config InlineBackend.figure_formats = ['svg']
Running inside Jupyter
Small library for pretty-printing (nothing to do)
def nat2bl(pad,n):
assert(n >= 0)
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):
s = [] + s # Make sure to not change 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))
0 - Our very own notion of circuits¶
In this lab session, we discuss reversible circuits: circuits built from X, NOT and Toffoli gates. As these are bijections on basis ket vectors, simulating a reversible circuit only require Boolean values: it should be doable efficiently.
Unfortunately, Qiskit does not come (afaik) with a simulation backend for such reversible circuits. We therefore have to build our very own.
For the sake of simplicity, we here define a very simple class for circuits.
from functools import reduce
class RevCirc:
def __init__(self,inputs):
self.gates = []
self.ancillas = set()
self.inputs = inputs
self.wires = set(inputs)
self.ctl = []
def x(self,i): self.gates.append({"target" : i, "ctl" : [] + self.ctl})
def cx(self,j,i): self.gates.append({"target" : i, "ctl" : [j] + self.ctl})
def ccx(self,j,k,i): self.gates.append({"target" : i, "ctl" : [j,k] + self.ctl})
def alloc(self):
w = max(list(self.wires) + [0]) + 1
self.ancillas.add(w)
self.wires.add(w)
return w
def allocn(self,n):
return [self.alloc() for _ in range(n)]
def append(self, circ, inputs):
assert(len(circ.inputs) == len(inputs))
ancillas = self.allocn(len(circ.ancillas))
renaming = dict(zip(circ.inputs + list(circ.ancillas), inputs + ancillas))
for g in circ.gates:
self.gates.append({"target" : renaming[g["target"]],
"ctl" : [renaming[c] for c in g["ctl"]] + self.ctl})
return renaming
def run(self,state):
assert(set(self.inputs) == set(state.keys()))
for a in self.ancillas: state[a] = 0
for g in self.gates: state[g["target"]] ^= reduce(lambda x, y: x & y, [1] + [state[i] for i in g["ctl"]])
return state
def draw(self):
height = max(self.wires)
res_string = ""
def isIn(wn, gn): return wn < max(gn["ctl"] + [gn["target"]]) and wn >= min(gn["ctl"] + [gn["target"]])
for wn in sorted(list(self.wires)):
if wn in self.ancillas: res_string += "a" + str(wn) + (" " * (len(str(height)) - len(str(wn)))) + " "
else: res_string += "i" + str(wn) + (" " * (len(str(height)) - len(str(wn)))) + " "
for gn in self.gates:
if gn["target"] == wn: res_string += "--X--"
elif wn in gn["ctl"]: res_string += "--.--"
elif isIn(wn,gn): res_string += "--|--"
else: res_string += "-----"
res_string += "\n " + (" " * (len(str(height))))
for gn in self.gates:
if isIn(wn,gn): res_string += " | "
else: res_string += " "
res_string += "\n"
print(res_string)
def getValue(valuation, register):
return [valuation[wire] for wire in register]
Reversible circuits are made of wires, here simply natural numbers. There are two types of wires: inputs wires, set at init time, and ancillas that can be added on the fly. Wires cannot be erased: they are all accessible at the end of the circuit through their ID number. The operations on circuits are
circ.x(q): apply an X gate onqcirc.cx(q1,q2): apply a CNOT onq2, the control isq1circ.ccx(q1,q2,q3)~ : apply a Toffoli gate onq3, the controls areq1andq2`q = circ.alloc(): allocate a new ancillaq(assumed initialized at 0)q = circ.allocn(n): allocatennew ancillas and store them in listq(they are assumed initialized at 0)circ.draw(): pretty-print the circuitvaluation = circ.run(input_dict): run the circuit with a dictionary of the form{q1 : 0, q2: 1, ...}whereq1,q2, etc are the input qubits. It returns a dictionaryvaluationstating the value of each wire of the circuit at the end of the computation.
The function getValue(valuation, register) is a wrapper to get the bitstring value of a register, given a valuation coming from runnning a circuit.
A small example of circuit is as follows.
def example():
circ = RevCirc([1,2])
circ.x(1)
q = circ.alloc()
circ.cx(1,q)
circ.ccx(1,2,q)
circ.draw()
input_val = dict(zip([1,2], [0,0])) # wires and 2 are set to 0
print("input valuation:" , input_val)
val = circ.run(input_val)
print("output valuation:", val)
print(f"The register {str([1,2,q])} is in state", getValue(val, [1,2,q]))
example()
i1 --X----.----.--
| |
i2 -------|----.--
| |
a3 -------X----X--
input valuation: {1: 0, 2: 0}
output valuation: {1: 1, 2: 0, 3: 1}
The register [1, 2, 3] is in state [1, 0, 1]
1 - A Ripple-Carry Adder¶
In this section, we will define the $V$-version of an adder: a circuit of shape
$$ V_{\text{adder}} : |{a}\rangle_n\otimes|{b}\rangle_n\otimes|{0\ldots0}\rangle_N\otimes|{0}\rangle_1\otimes|{0}\rangle_n \mapsto |{a}\rangle_n\otimes|{b}\rangle_n\otimes|{\text{garbage}}\rangle_N\otimes|{\text{carry}}\rangle_1\otimes|{a+b}\rangle_n $$
The notation $|{-}\rangle_k$ means that the register is of size $k$. In particular, our adder stores the last carry coming from the addition (if any) in an auxiliary register.
1.1 - Building-Blocks¶
Remember the building blocks HA (half-adder) and FA (full-adder) for an adder, as discussed in Section 5.3.9 of the lecture notes (and on the board during last lecture).
TODO Implement HA and FA with the following specification:
s,c = HA(circ, a, b)inputs a circuitcirc, two existing wiresaandbofcirc. It creates two ancillassandcand setssto $\texttt{a} \oplus\texttt{b}$ andcto $\texttt{a}\texttt{b}$.s,c = FA(circ, a, b, c)has a similar structure, but it setssandcto
$$ s = a\oplus b\oplus c \qquad c = ab\oplus ac\oplus bc $$
def HA(circ, a, b):
s = circ.alloc()
c = circ.alloc()
## TODO BEGIN
circ.cx(a,s)
circ.cx(b,s)
circ.ccx(a,b,c)
## TODO END
return (s,c)
def FA(circ, cin, a, b):
s = circ.alloc()
cout = circ.alloc()
## TODO BEGIN
circ.cx(cin,s)
circ.cx(a,s)
circ.cx(b,s)
circ.ccx(a,b,cout)
circ.ccx(cin,a,cout)
circ.ccx(cin,b,cout)
## TODO END
return (s,cout)
## Printing the circuit, to check
def test_HA_FA():
print("Half-adder")
circ = RevCirc([0,1])
HA(circ,0,1)
circ.draw()
print("Full-adder")
circ = RevCirc([0,1,2])
FA(circ,0,1,2)
circ.draw()
test_HA_FA()
Half-adder
i0 --.---------.--
| |
i1 --|----.----.--
| | |
a2 --X----X----|--
|
a3 ------------X--
Full-adder
i0 --.-------------------.----.--
| | |
i1 --|----.---------.----.----|--
| | | | |
i2 --|----|----.----.----|----.--
| | | | | |
a3 --X----X----X----|----|----|--
| | |
a4 -----------------X----X----X--
1.2 - Composing the Pieces Together¶
The complete adder (in $V$-form) is one half-adder followed by a chain of full-adder, each one sending its carry to the next full adder. The last carry is left in its own register.
We need to choose a binary encoding for natural numbers: let us set
the least significant bit on the right. This is the convention used in
nat2bl and bl2nat.
TODO Implement the $V$-form of the adder.
def adderV(circ, aa, bb): # assume aa and bb have same size
res = [] ## TO UPDATE
c = circ.alloc() ## TO UPDATE
## TODO BEGIN
s,c = HA(circ, aa[-1], bb[-1])
res.append(s)
for i in range(len(aa)-2,-1,-1):
s,c = FA(circ,c,aa[i],bb[i])
res.append(s)
res = list(reversed(res))
## TODO END
return res,c
Testing cell.
def testAdderV():
# Size of the bitstring to encode numbers
# Change to test
N = 4
# What to add together -- should fit on N bits
# Change to test
a = 10
b = 3
# Do not change the rest
ar = list(range(0,N))
br = list(range(N,2*N))
circ = RevCirc(ar + br)
s,c = adderV(circ,ar,br)
circ.draw()
ai = nat2bl(N,a)
bi = nat2bl(N,b)
print(f"{a} in binary is", ai)
print(f"{b} in binary is", bi)
d = dict(zip(ar + br, ai + bi))
output = circ.run(d)
si = getValue(output, s)
si_str = str(si)
s = bl2nat(si)
print(f"The sum of {a} and {b} is {a+b}")
print(f"The answer given by the circuit is {s}")
print("In binary form:", si_str)
print(f"The carry bit is {d[c]}")
testAdderV()
i0 ----------------------------------------------------------------------------------.---------.----.-------
| | |
i1 ----------------------------------------------------.---------.----.--------------|---------|----|-------
| | | | | |
i2 ----------------------.---------.----.--------------|---------|----|--------------|---------|----|-------
| | | | | | | | |
i3 --.---------.---------|---------|----|--------------|---------|----|--------------|---------|----|-------
| | | | | | | | | | |
i4 --|---------|---------|---------|----|--------------|---------|----|--------------|----.----.----|----.--
| | | | | | | | | | | | |
i5 --|---------|---------|---------|----|--------------|----.----.----|----.---------|----|----|----|----|--
| | | | | | | | | | | | | | |
i6 --|---------|---------|----.----.----|----.---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | | | | | |
i7 --|----.----.---------|----|----|----|----|---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | | | | | | |
a8 --|----|----|---------|----|----|----|----|---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | | | | | | |
a9 --X----X----|---------|----|----|----|----|---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | | | | |
a10 ------------X----.----|----|----|----.----.---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | | | | |
a11 -----------------X----X----X----|----|----|---------|----|----|----|----|---------|----|----|----|----|--
| | | | | | | | | | | | |
a12 --------------------------------X----X----X----.----|----|----|----.----.---------|----|----|----|----|--
| | | | | | | | | | |
a13 -----------------------------------------------X----X----X----|----|----|---------|----|----|----|----|--
| | | | | | | |
a14 --------------------------------------------------------------X----X----X----.----|----|----|----.----.--
| | | | | |
a15 -----------------------------------------------------------------------------X----X----X----|----|----|--
| | |
a16 --------------------------------------------------------------------------------------------X----X----X--
10 in binary is [1, 0, 1, 0]
3 in binary is [0, 0, 1, 1]
The sum of 10 and 3 is 13
The answer given by the circuit is 13
In binary form: [1, 1, 0, 1]
The carry bit is 0
2 - Multiplication¶
Consider $a$ and $b$ two natural numbers coded in binary on $5$ bits. The multiplication $a\cdot b$ of $a$ and $b$ fits on $10$ bits. The procedure is as you learned in elementary school.
Suppose that $a$ is $a_1a_2a_3a_4a_5$ and $b$ is $b_1b_2b_3b_4b_5$ in binary, least significant bit on the right. We can write $a\cdot b$ as follows:
$$ a \cdot b = b_1 \cdot 2^4 \cdot a + b_2 \cdot 2^3 \cdot a + b_3 \cdot 2^2 \cdot a + b_4 \cdot 2^1 \cdot a + b_5 \cdot 2^0 \cdot a $$
The multiplication of $a$ by $2^k$ corresponds to padding $a_1a_2a_3a_4a_5$ with $0$'s on the left. One simple possibility for implementing $a\cdot b$ is therefore to consider a series of registers of size $10$ as follows:
$$ \begin{array}{lllllll lllll} b_5 & {\cdot} & 0 & 0 & 0 & 0 & 0 & a_1 & a_2 & a_3 & a_4 & a_5 \\ b_4 & {\cdot} & 0 & 0 & 0 & 0 & a_1 & a_2 & a_3 & a_4 & a_5 & 0 \\ b_3 & {\cdot} & 0 & 0 & 0 & a_1 & a_2 & a_3 & a_4 & a_5 & 0 & 0 \\ b_2 & {\cdot} & 0 & 0 & a_1 & a_2 & a_3 & a_4 & a_5 & 0 & 0 & 0 \\ b_1 & {\cdot} & 0 & a_1 & a_2 & a_3 & a_4 & a_5 & 0 & 0 & 0 & 0 \\ \end{array} $$
The operation $b_k\cdot \texttt[register]$ is to be understood as: leave the corresponding register to $00000 00000$ when $b_k$ is $0$.
Once these registers are set up, it is just a matter of adding them
together with the function adderV you defined above.
TODO Implement the $V$-form of the multiplication of two registers of same size, as described in this section.
# Multiplication
def mult(circ, aa, bb):
# This is what I mean when I suggest to allocate auxiliary registers
aux = [circ.allocn(len(aa) + len(bb)) for i in range(len(bb))]
# You now have to fill them with the correct bitstring, add them
# together, and return the final register with the last carry
ss = [] # Placeholder for the final register, to redefine correctly
## TODO BEGIN
for i in range(len(bb)):
for j in range(len(aa)):
circ.ccx(aa[j], bb[-i-1], aux[i][j + len(bb)-i])
ss = aux[0]
for i in range(1,len(bb)):
ss,c = adderV(circ,ss,aux[i])
## TODO END
return ss # Return the final register
Testing cell.
def testMult():
# Size of registers (to play with)
N = 6
# Values to multiply (make sure they fit in registers of size N)
a = 33
b = 48
# Do not modify below
ar = list(range(0,N))
br = list(range(N,2*N))
circ = RevCirc(ar + br)
s = mult(circ,ar,br)
ai = nat2bl(N,a)
bi = nat2bl(N,b)
d = dict(zip(ar + br, ai + bi))
print("nber of gates:", len(circ.gates))
output = circ.run(d)
si = getValue(output, s)
print(f"First value is {a}, the register is", ai)
print(f"Second value is {b}, the register is", bi)
print(f"The circuit outputs {str(si)}, meaning {bl2nat(si)}")
print(f"The correct value should have been {a*b}, not counting the overflow")
testMult()
nber of gates: 381 First value is 33, the register is [1, 0, 0, 0, 0, 1] Second value is 48, the register is [1, 1, 0, 0, 0, 0] The circuit outputs [0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0], meaning 1584 The correct value should have been 1584, not counting the overflow
3 - Signed, fixed point representation of real numbers¶
The two functions you just wrote can be composed to generate polynomials. An interesting class of polynomials is the class of Taylor polynomials, for approximating functions such as $x \mapsto arcsin(\frac{1}{16x}$.
However, such functions works over real numbers... Going from bitstring representing natural numbers to bitstring representing real numbers can be done in two steps:
- Signed Representation. One consider the high bit as the bit representing the sign. So for instance, signed integers on $3$ bits are
$$ \begin{array}{lll} 0 0 0 & \simeq & 0\\ 0 0 1 & \simeq & 1\\ 0 1 0 & \simeq & 2\\ 0 1 1 & \simeq & 3\\ 1 0 0 & \simeq & -4\\ 1 0 1 & \simeq & -3\\ 1 1 0 & \simeq & -2\\ 1 1 1 & \simeq & -1\\ \end{array} $$
Using this convention, addition of signed integers can be done with the standard addition of unsigned integers you just defined. For multiplication, there is a slight twist but regular multiplication can also be reused almost directly.
- Fixed Point Representation: Some of the bits in the bitstring now represent fractional numbers. Said otherwise, the least significant bit does not represent the number $1$ but instead the number $\frac1{2^k}$ for some fixed $k$ (thus "fixed point").
On $3$ bits with a fractional part of $1$ bit, the possible bitstring are then
$$ \begin{array}{lll} 0 0 0 & \simeq & 0\\ 0 0 1 & \simeq & \frac12\\ 0 1 0 & \simeq & 1\\ 0 1 1 & \simeq & \frac32\\ 1 0 0 & \simeq & -2\\ 1 0 1 & \simeq & -\frac32\\ 1 1 0 & \simeq & -1\\ 1 1 1 & \simeq & -\frac12 \end{array} $$
On 6 bits with a fractional part of $3$ bits, we have for instance
$$ \begin{array}{lll} 0 ~~ 0 0 ~~ 0 0 0 & \simeq & 0\\ 0 ~~ 0 1 ~~ 0 0 1 & \simeq & 1 + \frac18\\ 0 ~~ 1 1 ~~ 0 0 0 & \simeq & 3\\ 1 ~~ 1 1 ~~ 1 1 1 & \simeq & -\frac18\\ 1 ~~ 1 1 ~~ 0 0 0 & \simeq & -1\\ 1 ~~ 1 0 ~~ 1 1 1 & \simeq & -1 - \frac18 \end{array} $$
3.1 - Lifting operations to Signed FP reals (nothing to do)¶
This section defines for you arithmetic operations on signed FP reals:
# Make a signed FP (Fixed Point) bitstring for real_nb
def makeSignedFP(n_frac, N, real_nb):
int_nb = int(real_nb * (2 ** n_frac))
if int_nb < 0:
binary = nat2bl(N, (2 ** N) + (int_nb))
else:
binary = nat2bl(N, int_nb)
return binary
for x in [0, 1+1/8, 3, -1/8, -1, -1-1/8]:
print("%s corresponds to %f" % (str(makeSignedFP(3, 6, x)), x))
[0, 0, 0, 0, 0, 0] corresponds to 0.000000 [0, 0, 1, 0, 0, 1] corresponds to 1.125000 [0, 1, 1, 0, 0, 0] corresponds to 3.000000 [1, 1, 1, 1, 1, 1] corresponds to -0.125000 [1, 1, 1, 0, 0, 0] corresponds to -1.000000 [1, 1, 0, 1, 1, 1] corresponds to -1.125000
# Opposite operation: get the real number from a bitstring
def readSignedFP(bitstring, n_frac):
si = bl2nat(bitstring)
sign = 1
if len(bitstring) > 0 and bitstring[0] == 1:
si = (2 ** (len(bitstring))) - si
sign = -1
return sign * si/(2 ** n_frac)
print(readSignedFP([1,0,1], 1))
-1.5
# getValue, but for real numbers in FP notation
def getSignedFP(reg, val, n_frac):
return readSignedFP(getValue(val,reg),n_frac)
# Initialize a register in the state corresponding to a real in FP notation.
def constFPsigned(circ, n_frac, N, real_nb):
binary = makeSignedFP(n_frac, N, real_nb)
s = circ.allocn(N)
for i in range(N):
if binary[i] == 1:
circ.x(s[i])
return s
# FP Signed Addition: nothing much to do beside forgetting the carry
def adderFPSigned(circ, n_frac, aa, bb):
ss,c = adderV(circ,aa,bb)
return ss
# FP Signed Multiplication
def multFPSigned(circ, n_frac, aa, bb):
na = len(aa)
nb = len(bb)
# signed extension
aa_ext = circ.allocn(nb) + aa
bb_ext = circ.allocn(na) + bb
for i in range(nb):
circ.cx(aa[0],aa_ext[i])
for i in range(na):
circ.cx(bb[0],bb_ext[i])
# Multiplication (omit the carry -- overflow if too large)
ss_ext = mult(circ, aa_ext, bb_ext)
n_int = len(aa)-n_frac
# Take relevant part
return ss_ext[-na-nb:][n_int:len(aa)+n_int]
In the following cell, we set up a circuit computing
$$ x \mapsto 2.5 * x - 1.25 $$
and run it with $x = 3.5$. Do you read the correct value ?
def exampleFPSigned():
# 2 bits for the fractional part: enough as we won't go beyond 1/4
# in precision
n_frac = 2
# 6 in total, leaving 3 for the integer part, enough for our
# computation
N = 6
# Input register
x = [i for i in range(N)]
# Set up a circuit
circ = RevCirc(x)
# Append all of the pieces, V-style
coef0 = constFPsigned(circ, n_frac, N, -1.25)
coef1 = constFPsigned(circ, n_frac, N, 2.5)
monom = multFPSigned(circ, n_frac, coef1, x)
result_reg = adderFPSigned(circ, n_frac, coef0, monom)
# Define an input
input = makeSignedFP(n_frac, N, 3.5)
# Input valuation
input_valuation = dict(zip(x, input))
# Run the circuit
valuation = circ.run(input_valuation)
# print result!
print("Resulting FP real: ", getSignedFP(result_reg, valuation, n_frac))
exampleFPSigned()
Resulting FP real: 7.5
3.2 - Polynomials¶
A polynomial is a function of the form
$$ x \mapsto 3\cdot x^2 - 2.5\cdot x + 1.125 $$
3.2.1 - Power function¶
One simple way to build a circuit for such a function is to define
$$ x \mapsto x^k $$
for a fixed $k$.
TODO Fill-in the following code.
aais a register containing a FP signed real numbernis a (regular, Python) positive integer
def powerFPSigned(circ, n_frac, aa, n):
assert(n > 0)
ss = circ.allocn(len(aa)) # Placeholder for the result register,
# possibly to change
## TODO BEGIN
for i in range(len(ss)):
circ.cx(aa[i],ss[i])
for i in range(n-1):
ss = multFPSigned(circ,n_frac, aa,ss)
## TODO END
return ss
Testing cell: you can play with N, n_frac, p and input_x.
Note: $2.5^3 = 15.625$, it should fit in N=8 and n_frac=3.
def testPowerFPSigned():
# Shape of the register
N = 8
n_frac = 3
# Power
p = 3
# value to power
input_x = 2.5
## Do not change below
# Input register
x = [i for i in range(N)]
# Set up a circuit
circ = RevCirc(x)
# Power of x
result_reg = powerFPSigned(circ, n_frac, x, p)
# Define an input
input_bl = makeSignedFP(n_frac, N, input_x)
# Input valuation
input_valuation = dict(zip(x, input_bl))
# Run the circuit
valuation = circ.run(input_valuation)
# print result!
print("Resulting FP real: ", getSignedFP(result_reg, valuation, n_frac))
testPowerFPSigned()
Resulting FP real: 15.625
3.2.2 - Implementing Polynomials¶
You are now ready to implement a function such as
$$ x \mapsto 3\cdot x^2 - 2.5\cdot x + 1.125 $$
The coefficients of the above polynomial are coded as
[3, -2.5, 1.125]
Said otherwise, the coeff of the monomials of highest degree is at the beginning.
TODO Fill-in the following code. The function returns a circuit realizing the polynomial.
def polyCirc(coeffs, n_frac, N):
inputs = [i for i in range(N)]
circ = RevCirc(inputs)
res = [] # placeholder for the register wires holding the
# answer. Change accordingly
## TODO BEGIN
exp = len(coeffs) - 1
res = constFPsigned(circ, n_frac, N, coeffs[-1])
for c in coeffs[:-1]:
bin_c = constFPsigned(circ, n_frac, N, c)
monom = powerFPSigned(circ, n_frac, inputs, exp)
monom = multFPSigned(circ, n_frac, bin_c, monom)
res = adderFPSigned(circ, n_frac, res, monom)
exp = exp - 1
## TODO END
return res, circ
Testing cell: you can play with N, n_frac, p and input_x.
def testPoly():
# Shape of registers
N = 8
n_frac = 3
# coeffs of the polynomial
coeffs = [3, -2.5, 1.125]
# input value
input_x = 1.5
# Do not change: setting up the circuit and running it
output, circ = polyCirc(coeffs, n_frac, N)
a = makeSignedFP(n_frac,N, input_x)
d = circ.run(dict(zip(circ.inputs, a)))
print("The circuit says: ", getSignedFP(output, d, n_frac))
exact_res = 0
exp = 0
for c in reversed(coeffs):
exact_res = exact_res + c * (input_x ** exp)
exp += 1
print("The exact result should be", exact_res)
testPoly()
The circuit says: 4.125 The exact result should be 4.125
4 - Implementing the Inversion Operator of HHL¶
We are ready to approximate the inversion operator ! Following what we already discussed, we shall focus on
$$ Inv : x \mapsto \arcsin\left(\frac1{16x}\right) $$
and we will use a Taylor expansion of the function.
4.1 - Taylor expansion of $Inv$ (nothing to do)¶
The function $Inv$ works on $[\frac1{16},1]$. We shall focus on the Taylor expansion centered $x=0.5$.
The library sympy is our friend here: we can compute the successive
derivatives of arcsin algebraically
# the function we care about
def Inv(x):
return np.arcsin(1/(16 * x))
# The degree of the polynomial
degree = 3
# Point at which computing the Taylor polynomial
at_x = 0.5
# The Taylor polynomial !
xvar = sympy.symbols('x')
# At degree 0
poly = Inv(at_x)
# Degrees 1, 2, ...
deriv = sympy.functions.elementary.trigonometric.asin(1/(16*xvar))
for k in range(1,degree+1):
deriv = sympy.diff(deriv, xvar)
poly += deriv.subs(xvar,at_x)/math.factorial(k) * (xvar - at_x) ** k
def Inv_approx(x):
return float(poly.subs(xvar,x))
print("The Taylor polynomial is")
print(" ", sympy.poly(poly))
# The list of coefficients using the convention we discussed when
# programming polynomials.
coeffs = sympy.poly(poly).all_coeffs()
print("The list of coefficients is")
print(" ", coeffs)
The Taylor polynomial is
Poly(-1.02669714579747*x**3 + 2.04799797342818*x**2 - 1.52995142941956*x + 0.506641195745486, x, domain='RR')
The list of coefficients is
[-1.02669714579747, 2.04799797342818, -1.52995142941956, 0.506641195745486]
4.2 - The Circuit for the Taylor Approximation¶
TODO
Generate a circuit for the Taylor approximation computed in the previous cell, using registers of size
N=16withn_frac=10bits for the fractional part. Note how the coefficients of the polynomial are stored incoeffs, so you can reusepolyCirc.Run it for several values between $0.1$ and $0.9$ (take, say, 50 points). Compare the number of bits of precision between
- the Taylor expansion and the function computed by the circuit
- The Taylor expansion and the "correct" $Inv$ function
You can use -numpy.log2(abs(exact_value - approx_value)) to compute
how many bits two values differ from.
- How does the error and the size of the circuit evolve with respect to
- the size of the register?
- the degree of the polynomial?
## TODO SOLUTION BEGIN
def runTaylorCircuit(N, n_frac):
output, circ = polyCirc(coeffs, n_frac, N)
xx = np.linspace(0.1, 0.9, num=500)
errors_circ_tayl = []
errors_tayl_exac = []
errors_circ_exac = []
print(" x | bits of precision | ")
print(" |---------------------------------------------------| ")
print(" | circ vs Taylor | Taylor vs exact | circ vs exact |")
counter_to_print = 0
for x in xx:
a = makeSignedFP(n_frac,N, x)
d = circ.run(dict(zip(circ.inputs, a)))
inv_circ = getSignedFP(output, d, n_frac)
inv_tayl = Inv_approx(x)
inv_exac = Inv(x)
error_circ_tayl = -np.log2(abs(inv_tayl - inv_circ))
error_tayl_exac = -np.log2(abs(inv_exac - inv_tayl))
error_circ_exac = -np.log2(abs(inv_exac - inv_circ))
if counter_to_print == 0: # print one every 10 values
print(" %8.5f | %8.5f | %8.5f | %8.5f |" %(x, error_circ_tayl, error_tayl_exac, error_circ_exac))
counter_to_print = 10
else:
counter_to_print -= 1
errors_circ_tayl.append(error_circ_tayl)
errors_tayl_exac.append(error_tayl_exac)
errors_circ_exac.append(error_circ_exac)
print("Size of circuit:", len(circ.gates))
plt.plot(xx, errors_circ_tayl, label="Prec Circ/Taylor (bits)")
plt.plot(xx, errors_tayl_exac, label="Prec Taylor/Exact (bits)")
plt.plot(xx, errors_circ_exac, label="Prec Circ/Exact (bits)")
plt.xlabel('x')
plt.ylabel("Bits of precision")
plt.title("Plot of precision (in bits)")
plt.legend()
plt.grid(True)
plt.show()
runTaylorCircuit(N=16, n_frac=10)
x | bits of precision |
|---------------------------------------------------|
| circ vs Taylor | Taylor vs exact | circ vs exact |
0.10000 | 14.21841 | 1.72723 | 1.72698 |
0.11764 | 9.14085 | 2.27366 | 2.26136 |
0.13527 | 8.13294 | 2.77919 | 2.74433 |
0.15291 | 8.66609 | 3.26303 | 3.22933 |
0.17054 | 8.73799 | 3.73639 | 3.69204 |
0.18818 | 9.68315 | 4.20698 | 4.17493 |
0.20581 | 9.95990 | 4.68077 | 4.64409 |
0.22345 | 8.56094 | 5.16292 | 5.03217 |
0.24108 | 8.82318 | 5.65828 | 5.50577 |
0.25872 | 10.64027 | 6.17176 | 6.10803 |
0.27635 | 8.50467 | 6.70870 | 6.34359 |
0.29399 | 9.18502 | 7.27517 | 6.93476 |
0.31162 | 9.48707 | 7.87843 | 7.46927 |
0.32926 | 8.52092 | 8.52748 | 7.52420 |
0.34689 | 9.14642 | 9.23399 | 8.18954 |
0.36453 | 8.21335 | 10.01358 | 7.84919 |
0.38216 | 10.55091 | 10.88822 | 9.70973 |
0.39980 | 8.69563 | 11.89033 | 8.54609 |
0.41743 | 9.17054 | 13.07106 | 9.07703 |
0.43507 | 9.34942 | 14.51856 | 9.30987 |
0.45271 | 9.14123 | 16.40579 | 9.13188 |
0.47034 | 8.69345 | 19.15467 | 8.69243 |
0.48798 | 8.66932 | 24.41893 | 8.66929 |
0.50561 | 10.14551 | 28.86912 | 10.14550 |
0.52325 | 9.31704 | 20.71684 | 9.31651 |
0.54088 | 9.33577 | 17.50759 | 9.33077 |
0.55852 | 9.14613 | 15.48494 | 9.12841 |
0.57615 | 8.83734 | 14.01027 | 8.79789 |
0.59379 | 10.24435 | 12.85230 | 10.02522 |
0.61142 | 8.62748 | 11.90074 | 8.48548 |
0.62906 | 11.18706 | 11.09436 | 10.13997 |
0.64669 | 8.25175 | 10.39558 | 7.95746 |
0.66433 | 9.60431 | 9.77971 | 8.68934 |
0.68196 | 8.81554 | 9.22964 | 8.00778 |
0.69960 | 8.86476 | 8.73305 | 7.79740 |
0.71723 | 9.92650 | 8.28074 | 7.88066 |
0.73487 | 9.14773 | 7.86569 | 7.36876 |
0.75251 | 9.49291 | 7.48243 | 7.16259 |
0.77014 | 8.48137 | 7.12657 | 6.65046 |
0.78778 | 8.92772 | 6.79459 | 6.49832 |
0.80541 | 8.93010 | 6.48358 | 6.24058 |
0.82305 | 9.13590 | 6.19114 | 6.01497 |
0.84068 | 9.70030 | 5.91524 | 5.81421 |
0.85832 | 9.62509 | 5.65419 | 5.56500 |
0.87595 | 9.01217 | 5.40651 | 5.29261 |
0.89359 | 8.82647 | 5.17095 | 5.06078 | Size of circuit: 77558
runTaylorCircuit(N=21, n_frac=15)
## SOLUTION END
x | bits of precision |
|---------------------------------------------------|
| circ vs Taylor | Taylor vs exact | circ vs exact |
0.10000 | 14.21841 | 1.72723 | 1.72698 |
0.11764 | 13.96645 | 2.27366 | 2.27323 |
0.13527 | 13.54993 | 2.77919 | 2.77836 |
0.15291 | 13.58456 | 3.26303 | 3.26190 |
0.17054 | 13.09484 | 3.73639 | 3.73419 |
0.18818 | 14.10515 | 4.20698 | 4.20547 |
0.20581 | 13.46294 | 4.68077 | 4.67750 |
0.22345 | 14.17933 | 5.16292 | 5.16014 |
0.24108 | 16.53569 | 5.65828 | 5.65751 |
0.25872 | 13.17993 | 6.17176 | 6.16060 |
0.27635 | 14.71631 | 6.70870 | 6.70311 |
0.29399 | 13.27893 | 7.27517 | 7.25286 |
0.31162 | 13.12723 | 7.87843 | 7.84097 |
0.32926 | 13.85802 | 8.52748 | 8.49207 |
0.34689 | 15.28076 | 9.23399 | 9.21233 |
0.36453 | 13.73443 | 10.01358 | 9.90811 |
0.38216 | 15.24658 | 10.88822 | 10.81954 |
0.39980 | 14.95444 | 11.89033 | 11.72739 |
0.41743 | 13.48163 | 13.07106 | 12.26179 |
0.43507 | 14.69715 | 14.51856 | 13.60509 |
0.45271 | 13.97734 | 16.40579 | 13.73152 |
0.47034 | 13.34381 | 19.15467 | 13.31834 |
0.48798 | 13.68561 | 24.41893 | 13.68476 |
0.50561 | 15.10473 | 28.86912 | 15.10463 |
0.52325 | 13.75283 | 20.71684 | 13.74132 |
0.54088 | 13.10813 | 17.50759 | 13.04134 |
0.55852 | 14.12412 | 15.48494 | 13.64970 |
0.57615 | 13.13665 | 14.01027 | 12.50832 |
0.59379 | 21.13253 | 12.85230 | 12.84766 |
0.61142 | 13.48641 | 11.90074 | 11.48588 |
0.62906 | 13.38918 | 11.09436 | 10.82677 |
0.64669 | 13.67558 | 10.39558 | 10.25421 |
0.66433 | 13.93079 | 9.77971 | 9.70071 |
0.68196 | 14.21004 | 9.22964 | 9.18465 |
0.69960 | 13.80500 | 8.73305 | 8.69078 |
0.71723 | 14.25804 | 8.28074 | 8.25802 |
0.73487 | 14.17550 | 7.86569 | 7.84762 |
0.75251 | 14.43661 | 7.48243 | 7.47084 |
0.77014 | 14.24656 | 7.12657 | 7.11624 |
0.78778 | 13.80590 | 6.79459 | 6.78345 |
0.80541 | 14.76455 | 6.48358 | 6.47895 |
0.82305 | 14.68218 | 6.19114 | 6.18713 |
0.84068 | 13.23938 | 5.91524 | 5.90627 |
0.85832 | 14.41875 | 5.65419 | 5.65087 |
0.87595 | 14.45162 | 5.40651 | 5.40378 |
0.89359 | 14.76091 | 5.17095 | 5.16907 |
Size of circuit: 134552
SOLUTION BEGIN
We can see that the precision of the circuit tops at 10 bits, which is exactly how many bits we've set up for the fractional part of our signed FP reals. If we were to change this value, the prevision would vary accordingly.
If we were to increase the degree of the polynomial, the precision of the approximation would increase, but the precision given by the circuit would top up to 10 bits, as this is the overall precision chosen for signed FP reals. Testing things out, if we choose Taylor expansion of degree 5, we reach go beyond 10 bits of precisions in the interval $[0.2,0.8]$. Within this interval, going beyond is useless if we do not increase the size of the fractional part of the registers of our signed FP reals.
The size of the circuit increases polynomially:
- The size of the adder is of order $O(N)$
- The size of the multiplier is of order $O(N^2)$
- The size of "power $d$" is of order $O(N^2\cdot d)$
- As currently coded, a polynomial of degree $d$ contains $d + d^2 + d^3 + \cdots + d^(d+1) = O(d^d)$ multiplications.
We can estimate how things grow as follows.
SOLUTION END
## SOLUTION BEGIN
def TaylorPoly(degree):
xvar = sympy.symbols('x')
at_x = 0.5
poly = Inv(at_x)
deriv = sympy.functions.elementary.trigonometric.asin(1/(16*xvar))
for k in range(1,degree+1):
deriv = sympy.diff(deriv, xvar)
poly += deriv.subs(xvar,at_x)/math.factorial(k) * (xvar - at_x) ** k
coeffs = sympy.poly(poly).all_coeffs()
return coeffs
def TaylorCircuitSize():
N = 10
n_frac = 4
MAX_DEGREE = 6
MAX_N = 10
# size[i][degree] is the size of the circuit for
# degree d and N=10+i, n_frac=4+i
sizes = [[] for i in range(MAX_DEGREE)]
for degree in range(1,MAX_DEGREE):
for n_inc in range(0,MAX_N):
_,circ = polyCirc(TaylorPoly(degree), n_frac+n_inc, N+n_inc)
sizes[degree].append(len(circ.gates))
for d in range(1,MAX_DEGREE):
xs = [N+i for i in range(MAX_N)]
coeffs = [int(c * 1000)/1000 for c in np.polyfit(xs, sizes[d], 2)]
function = np.poly1d(coeffs)
print(f"Polynomial approximation for degree {d}: ", sympy.poly(function(xvar)))
plt.plot(xs, sizes[d], label=f"degree {d}")
plt.xlabel('total size of register')
plt.ylabel("Circuit size")
plt.title("Size of circuit versus total size of register")
plt.legend()
plt.grid(True)
plt.show()
TaylorCircuitSize()
## SOLUTION END
Polynomial approximation for degree 1: Poly(52.045*x**2 - 21.166*x + 5.084, x, domain='RR') Polynomial approximation for degree 2: Poly(156.026*x**2 - 69.314*x + 2.865, x, domain='RR') Polynomial approximation for degree 3: Poly(311.833*x**2 - 139.99*x - 30.115, x, domain='RR') Polynomial approximation for degree 4: Poly(520.189*x**2 - 254.759*x + 48.524, x, domain='RR') Polynomial approximation for degree 5: Poly(779.992*x**2 - 381.362*x + 19.906, x, domain='RR')
