Error propagation analysis for Fixed point arithmetic using fxpmath#
In this notebook, we use library fxpmath to analyze the propagation of errors due to fixed point arithmetic for algorithms described in Compilation of Fault-Tolerant Quantum Heuristics for Combinatorial Optimization.
Qualtran should provide tools to enable such analysis for algorithms that use fixed point arithmetic Bloqs as subroutines.
from typing import Optional
from fxpmath import Fxp
import numpy as np
# Initial setup - Helpers
def fxp(x: float, d: int, *, d_word: Optional[int]=None) -> Fxp:
"""Creates an unsigned fixed-point representation of `x` with `d` bits of precision after decimal."""
assert 0 <= x < 1
d_word = d if d_word is None else d_word
# Dtype Format: f'fxp-{signed/unsigned}{total_bitsize}/{frac_bitsize}'
dtype = f'fxp-u{d_word}/{d}'
# `op_sizing` and `const_op_sizing` controls the behavior of fixed point type when performing
# arithmetic with other fixed point types. Setting the value to `same` ensures that no addition
# bits are used during arithmetic.
# shifting='trunc' sets the behavior of fixed point type to truncate shifted bits when doing a
# bitwise left/right shift operation.
return Fxp(x, dtype=dtype, op_sizing='same', const_op_sizing='same', shifting='trunc')
def assert_allclose(x: Fxp, y: float, eps: float):
np.testing.assert_allclose(x.get_val(), y, atol=eps)
Appendix D4: Multiplying an integer to a real number#
Goal: Given quantum registers A and B with real number \(\kappa\) (\(0 \leq \kappa \lt 1\)) and a \(d_{B}\)-bit integer \(\lambda\), our goal is to compute an approximation \(\widetilde{\gamma}\) of \(\gamma = \kappa * \lambda\) s.t. \(|\widetilde{\gamma} - \gamma| < \epsilon\)
We analyze the error due to fixed point arithmetic for algorithm described in Appendix D4 of https://arxiv.org/abs/2007.07391
# Multiplying a real numbers with an d_B-bit integer.
def get_bitsize_for_fxp_mul_with_integer(eps: float, d_B: int):
return d_B + int(np.ceil(np.log2(d_B / eps))) # Equation D7
def mul_with_int_via_repeated_add(a: float, b: int, d_A: int, d_B: int):
"""Multiplicaiton via repeated additions algorithm described in Appendix D5"""
a_fxp = fxp(a, d=d_A, d_word=d_A+d_B-1) # The paper proposes to use `d_A-1` bits of `a` but we instead need to use `d_A` bits.
res = fxp(0, d=d_A-d_B, d_word=d_A)
for i in range(d_B):
b_i = (b >> i) & 1
a_i = (a_fxp << i).like(res)
res += a_i * b_i
return res
def test_multiplication_with_integer_for_eps(eps: float, d_B: int, d_A: int):
rng = np.random.default_rng(int(eps * d_B * 1e9))
try:
for _ in range(100):
a, = rng.random(1)
b = rng.integers(0, 1 << d_B)
res = mul_with_int_via_repeated_add(a, b, d_A, d_B)
assert_allclose(res, a * b, eps)
print(f'Success! {eps=}, {d_A=}, {d_B=}')
except AssertionError:
print(f'Failed! {eps=}, {d_A=}, {d_B=}')
for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-9]:
for d_B in [5, 8, 11]:
d_A = get_bitsize_for_fxp_mul_with_integer(eps, d_B)
test_multiplication_with_integer_for_eps(eps, d_B, d_A)
Success! eps=0.001, d_A=18, d_B=5
Success! eps=0.001, d_A=21, d_B=8
Success! eps=0.001, d_A=25, d_B=11
Success! eps=0.0001, d_A=21, d_B=5
Success! eps=0.0001, d_A=25, d_B=8
Success! eps=0.0001, d_A=28, d_B=11
Success! eps=1e-05, d_A=24, d_B=5
Success! eps=1e-05, d_A=28, d_B=8
Success! eps=1e-05, d_A=32, d_B=11
Success! eps=1e-06, d_A=28, d_B=5
Success! eps=1e-06, d_A=31, d_B=8
Success! eps=1e-06, d_A=35, d_B=11
Success! eps=1e-07, d_A=31, d_B=5
Success! eps=1e-07, d_A=35, d_B=8
Success! eps=1e-07, d_A=38, d_B=11
Success! eps=1e-09, d_A=38, d_B=5
Success! eps=1e-09, d_A=41, d_B=8
Success! eps=1e-09, d_A=45, d_B=11
Appendix D5: Multiplying two different real numbers#
Goal: Given quantum registers A and B with real numbers \(\kappa\) (\(0 \leq \kappa \lt 1\)) and \(\lambda\) (\(0 \leq \lambda \lt 1\)), our goal is to compute an approximation \(\widetilde{\gamma}\) of \(\gamma = \kappa * \lambda\) s.t. \(|\widetilde{\gamma} - \gamma| < \epsilon\)
We analyze the error due to fixed point arithmetic for algorithm described in Appendix D5 of https://arxiv.org/abs/2007.07391
# Multiplying two real numbers
def get_bitsize_for_fxp_mul(eps: float):
return int(np.ceil(1 + np.log2(1/eps)) + np.log2(1 + np.log2(1/eps))) # Equation D17
def mul_via_repeated_add(a: float, b: float, d: int):
"""Multiplicaiton via repeated additions algorithm described in Appendix D5"""
a_fxp, b_fxp = fxp(a, d=d), fxp(b, d=d)
res = fxp(0, d=d)
for i, b_bin in enumerate(b_fxp.bin()[:-1]):
a_fxp >>= 1
if int(b_bin):
res += a_fxp
return res
def test_multiplication_for_eps(eps: float, d: int):
rng = np.random.default_rng(int(eps * 1e9))
try:
for _ in range(100):
a, b = rng.random(2)
res = mul_via_repeated_add(a, b, d)
assert_allclose(res, a * b, eps)
print(f'Success! {d=}, {eps=}')
except AssertionError:
print(f'Failed! {d=}, {eps=}')
for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
d = get_bitsize_for_fxp_mul(eps)
test_multiplication_for_eps(eps, d - 1) # Works for d-1 as well, bounds in the paper are probably loose?
Success! d=13, eps=0.001
Success! d=17, eps=0.0001
Success! d=21, eps=1e-05
Success! d=24, eps=1e-06
Success! d=28, eps=1e-07
Appendix D6: Squaring a real number#
Goal: Given a quantum register A with a real number \(\kappa\) (\(0 \leq \kappa \lt 1\)), our goal is to compute an approximation \(\widetilde{\gamma}\) of \(\gamma = \kappa^2\) s.t. \(|\widetilde{\gamma} - \gamma| < \epsilon\)
We analyze the error due to fixed point arithmetic for algorithm described in Appendix D6 of https://arxiv.org/abs/2007.07391
# Squaring a real number
def get_bitsize_for_fxp_square(eps: float):
return int(np.ceil(np.log2(1/eps) + np.log2(11/3 + np.log2(1/eps)))) # Equation D36
def square_via_repeated_add(a: float, d: int):
a = fxp(a, d)
res = fxp(0, d=d)
a_bin = [int(x) for x in a.bin()]
one = fxp(0.5, d=d)
# Equation D23 & D29
for n in range(d//2, d-1):
res += (a >> n) * a_bin[n]
mask = fxp(0, d=d)
for n in range((d - 1) // 2):
res += (((a & mask) >> n) + ((one >> (2*n+1)) * a_bin[n])) * a_bin[n]
mask |= (one >> n)
if not d & 1: # Equation D29
n = d // 2 - 1
res += ((a & mask) >> n) * a_bin[n]
return res
def test_square_for_eps(eps: float, d: int):
rng = np.random.default_rng(int(eps * 1e9))
try:
for _ in range(100):
a, = rng.random(1)
res = square_via_repeated_add(a, d)
assert_allclose(res, a**2, eps)
print(f'Success! {d=}, {eps=}')
except AssertionError:
print(f'Failed! {d=}, {eps=}')
for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
d = get_bitsize_for_fxp_square(eps) # Need a +1 to make it work for smaller `eps`. Bound in the paper is too tight?
test_square_for_eps(eps, d)
Failed! d=14, eps=0.001
Failed! d=18, eps=0.0001
Failed! d=21, eps=1e-05
Success! d=25, eps=1e-06
Success! d=29, eps=1e-07