Skip to content

Commit

Permalink
Address review comments.
Browse files Browse the repository at this point in the history
  • Loading branch information
@fdmalone
fdmalone committed Oct 13, 2023
1 parent 0832f5a commit 388416a
Show file tree
Hide file tree
Showing 2 changed files with 112 additions and 87 deletions.
114 changes: 76 additions & 38 deletions qualtran/bloqs/chemistry/trotter.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# `Trotter Costs`\n",
"# Trotter Costs\n",
"\n",
"We want to estimate the cost of implementing time evolution of a wavefunction:\n",
"\n",
Expand All @@ -32,7 +32,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```First Quantized Grid Based Hamiltonian```\n",
"## First Quantized Grid Based Hamiltonian\n",
"\n",
"Ultimately we are interested in understanding the dynamics of real chemical systems. This requires us to study the ab-initio chemistry Hamiltonian which is given as (in atomic units)\n",
"\n",
Expand Down Expand Up @@ -63,6 +63,7 @@
"print(\"delta = \", delta)\n",
"x_scl = delta * x_int\n",
"assert (x_scl[-1] - x_scl[0]) - L < 1e-12\n",
"\n",
"print(f\"unscaled grid points = {x_int}\")\n",
"print(f\"scaled grid points = {x_scl}\")\n"
]
Expand Down Expand Up @@ -144,7 +145,13 @@
"$$\n",
"|\\psi(t)\\rangle \\approx \\mathrm{QFT} e^{-i\\delta t T} \\mathrm{QFT}^{\\dagger} e^{-i\\delta t U} e^{-i \\delta t V} |\\psi(0)\\rangle\n",
"$$\n",
"so that all the terms can be implemented via a gate which implements something of the form $e^{-i \\delta t \\phi_\\alpha }$ via phasing gates.\n",
"so that all the terms can be implemented via a gate which implements something of the form $e^{-i \\delta t \\phi_\\alpha }$ via phasing gates.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"We can represent our $\\eta$-electron wavefunction as \n",
"$$\n",
Expand All @@ -157,11 +164,23 @@
"\n",
"\\begin{align}\n",
"&\\sum_{r_1\\cdots r_\\eta} c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle \\\\\n",
" &\\rightarrow \\sum_{r_1\\cdots r_\\eta} c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle|V(r_1\\cdots r_\\eta)\\rangle \\hspace{10em} \\text{Compute pairwise potential in ancilla registers} \\\\\n",
" &\\rightarrow \\sum_{r_1\\cdots r_\\eta} e^{-i \\delta t V(r_1\\cdots r_\\eta)} c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle|V(r_1\\cdots r_\\eta)\\rangle \\hspace{5.5em} \\text{Phase the state with computed potential} \\\\ \n",
" &\\rightarrow \\sum_{r_1\\cdots r_\\eta} e^{-i \\delta t V(r_1\\cdots r_\\eta)} c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle|0\\cdots0\\rangle \\hspace{7.8em} \\text{Uncompute potential in ancilla register} \\\\ \n",
" \\rightarrow &\\sum_{r_1\\cdots r_\\eta} c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle \n",
" |V(r_1\\cdots r_\\eta)\\rangle\n",
" \\hspace{10.4em} \\text{Compute pairwise potential in ancilla registers} \\\\\n",
" \\rightarrow &\\sum_{r_1\\cdots r_\\eta} e^{-i \\delta t V(r_1\\cdots r_\\eta)}\n",
" c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle|V(r_1\\cdots r_\\eta)\\rangle\n",
" \\hspace{5.5em} \\text{Phase the state with computed potential} \\\\ \n",
" \\rightarrow &\\sum_{r_1\\cdots r_\\eta} e^{-i \\delta t V(r_1\\cdots r_\\eta)}\n",
" c(r_1, \\cdots, r_\\eta) |r_1\\cdots r_\\eta\\rangle|0\\cdots0\\rangle\n",
" \\hspace{7.8em} \\text{Uncompute potential in ancilla register} \\\\ \n",
"\\end{align}\n",
"in the above the ancilla register storing the value of the potential is of size $> 2n + 2$.\n",
"in the above the ancilla register storing the value of the potential is of size $> 2n + 2$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"To compute the potential we need to evaluate $\\frac{1}{r_{ij}} \\equiv \\frac{1}{|r_i - r_j|}$ which can be done in two steps:\n",
"\n",
Expand All @@ -187,7 +206,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```Function interpolation```\n",
"## Function interpolation\n",
"\n",
"The basic idea is to first fit a polynomial to $\\frac{1}{r_{ij}}$ to use as input for a Newton-Raphson step. In practice we do not want to try to fit $r_{ij}^{-1}$ over the entire range of values as this will lead to a poor quality fit. Instead we combine two ideas: piecewise interpolation and scaling. Piecewise interpolation allows us to boost the accuracy of our fit by splitting our domain into $k$ subintervals and fitting a polynomial in these smaller subintervals, which can then be pieced together on fly depending on the input value $\\bar{x}$. The second idea is to use the scaling properties of our function to only fit in a small subinterval, and appropriately scaling it if we need to know its value outside of this interval. That is, we know that $\\frac{1}{\\sqrt{x}} = \\frac{1}{\\sqrt{a}\\sqrt{x'}}$ for $x' \\in [1, 2]$, say.\n",
"\n",
Expand Down Expand Up @@ -229,8 +248,9 @@
"poly_fit_one = polynomial_approx_range_one(xs_one)\n",
"xs_two = np.linspace(1.5, 2.0, 10)\n",
"poly_fit_two = polynomial_approx_range_two(xs_two)\n",
"plt.plot(xs_one, np.abs(1.0 / xs_one**0.5 - poly_fit_one), marker='o', lw=0)\n",
"plt.plot(xs_two, np.abs(1.0 / xs_two**0.5 - poly_fit_two), marker='x', lw=0)\n",
"plt.plot(xs_one, np.abs(1.0 / xs_one**0.5 - poly_fit_one), marker='o', lw=0, label=\"range 1\")\n",
"plt.plot(xs_two, np.abs(1.0 / xs_two**0.5 - poly_fit_two), marker='x', lw=0, label=\"range 2\")\n",
"plt.legend()\n",
"plt.xlabel(\"$x$\")\n",
"plt.ylabel(\"Absolute Error\")\n",
"plt.yscale(\"log\")\n",
Expand All @@ -256,12 +276,36 @@
" yprime = 0.5 * y0 * (3 + delta - y0**2 * x) \n",
" return yprime\n",
"delta_range_one = 5.1642030908180720584e-9\n",
"delta_range_two = 3.6279794522852781448e-10\n",
"\n",
"nr_one = newton_raphson_step(xs_one, poly_fit_one, delta=delta_range_one)\n",
"nr_two = newton_raphson_step(xs_two, poly_fit_two, delta=delta_range_two)\n",
"plt.plot(xs_one, np.abs(1.0 / xs_one**0.5 - nr_one), marker='o', lw=0)\n",
"plt.plot(xs_two, np.abs(1.0 / xs_two**0.5 - nr_two), marker='x', lw=0)\n",
"delta_range_two = 3.6279794522852781448e-10"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nr_one = newton_raphson_step(xs_one, poly_fit_one, delta=delta_range_one)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nr_two = newton_raphson_step(xs_two, poly_fit_two, delta=delta_range_two)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.plot(xs_one, np.abs(1.0 / xs_one**0.5 - nr_one), marker='o', lw=0, label=\"range 1\")\n",
"plt.plot(xs_two, np.abs(1.0 / xs_two**0.5 - nr_two), marker='x', lw=0, label=\"range 2\")\n",
"plt.legend()\n",
"plt.xlabel(\"$x$\")\n",
"plt.ylabel(\"Absolute Error\")\n",
"plt.yscale(\"log\")\n",
Expand All @@ -274,7 +318,7 @@
"source": [
"Ok great! We have a very accurate approximation to our function using if we can load some polynomial coefficients and perform some arithmetic. Next we will describe how we can use variable spaced QROM to efficiently load only a small amount of data to enable this procedure. \n",
"\n",
"## ```Variable Spaced QROM and Fitting the Entire Domain```\n",
"## Variable Spaced QROM and Fitting the Entire Domain\n",
"\n",
"In this section we will discuss how we can use QROM along with function interpolation to approximate our function over the entire domain to high precision. But first let us convince ourselves that our two subintervals are sufficient and we can indeed fit our function by appropriately scaling our polynomial. For reasons that will become clear later, let us first group our allowed values of $r_{ij}^2$ in a logarithmic fashion as $0, 1, 2, 3, [4, 5], [6, 7], [8, 11], [12, 15], [16, 23], [24, 31], ...$ and so on. We can associate each of these integer intervals with our polynomial domains by noting that $[4, 5] = 2^2[1, 3/2]$, $[12, 15] = 2^3[3/2, 2]$, etc. Thus, we can evaluate our function for any value $r_{ij}^2$ after determining which \"bin\" it belongs to, determining the appropriate power of two, and scaling our polynomial appropriately. Let's test this below:"
]
Expand Down Expand Up @@ -375,7 +419,7 @@
"source": [
"Ok! You might be wondering what is special about the intervals we chose to bin our integer $r_{ij}^2$ values by, and how QROM is connected with this. We'll try to explain that now.\n",
"\n",
"### ```QROM```\n",
"### QROM\n",
"Recall that ultimately we need to load our polynomial coefficients onto the quantum computer before we can evaluate the inverse square root. Data loading can be achieved through QROM which implements: \n",
"\n",
"$$\n",
Expand All @@ -385,9 +429,9 @@
"That is, given a selection register $l$, QROM can load the binary representation of the $l$-th data element of $d$ into a target register of a given size.\n",
"In our case $|r_{ij}^2\\rangle$ is the selection register of size $2 n + 2$, and we want to load $b$-bit binary representations of our polynomial coefficients $\\{a_0, a_1, a_2, a_3\\}$ and $\\{b_0, b_1, b_2, b_3\\}$ for the two ranges.\n",
"\n",
"In principle we could load all of the (scaled) coefficients for each value of $r_{ij}$ in the register, but this would incur a cost that goes like $L-1$ Toffolis, where $L=2^{2n + 2}$ which would be unacceptably large. But we just saw that much of this data would be duplicated as we only need the coefficients for each subinterval, of which there are a logarithmic amount. Fortunately we can exploit the structure of the unary iteration circuit to exploit this redundancy and significantly reduce our costs and arrive at `variably space QROM`.\n",
"In principle we could load all of the (scaled) coefficients for each value of $r_{ij}$ in the register, but this would incur a cost that goes like $L-1$ Toffolis, where $L=2^{2n + 2}$ which would be unacceptably large. But we just saw that much of this data would be duplicated as we only need the coefficients for each subinterval, of which there are a logarithmic amount. Fortunately we can exploit the structure of the unary iteration circuit to exploit this redundancy and significantly reduce our costs and arrive at variably space QROM.\n",
"\n",
"To understand where the reduction in cost comes from, consider the `unary iteration` circuit which is used during QROM construction. There, depending on the particular binary representation of the selection register, different parts of the tree are traversed before writing the data to a register. As we only care about certain subintervals (which were conveniently chosen as powers of two), we can delete parts of our unary iteration circuit where we want to repeat the data for that range of integers. The number of allowed integers in each range is determined in the following way: for the starting integer of the range $l$: if $k$ is the most significant bit of $l$, then the number of integers in our range is $2^{k-2}$ (we only bin integers when $k \\ge 2$). This is identical to how we set up our binning procedure earlier. Below is some code to generate the ranges given the values of the most significant bits of the selection register. "
"To understand where the reduction in cost comes from, consider the unary iteration circuit which is used during QROM construction. There, depending on the particular binary representation of the selection register, different parts of the tree are traversed before writing the data to a register. As we only care about certain subintervals (which were conveniently chosen as powers of two), we can delete parts of our unary iteration circuit where we want to repeat the data for that range of integers. The number of allowed integers in each range is determined in the following way: for the starting integer of the range $l$: if $k$ is the most significant bit of $l$, then the number of integers in our range is $2^{k-2}$ (we only bin integers when $k \\ge 2$). This is identical to how we set up our binning procedure earlier. Below is some code to generate the ranges given the values of the most significant bits of the selection register. "
]
},
{
Expand Down Expand Up @@ -438,7 +482,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# ```Trotter Unitaries```\n",
"## Trotter Unitaries\n",
"\n",
"Ok, now that we have the setup straight in our heads lets build our bloqs and perform resource estimations!"
]
Expand All @@ -448,7 +492,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```Kinetic Energy Bloq```\n",
"## Kinetic Energy Bloq\n",
"Recall that the kinetic energy is diagonal in the momentum basis, which we are assuming our state is in following a ```QFT```. The basic algorithm is then\n",
"\n",
"1. For each electron $i$ compute $|\\mathbf{k}^2 \\rangle = |k_x^2 + k_y^2 + k_j^2\\rangle$ of size $2n + 2$.\n",
Expand Down Expand Up @@ -485,7 +529,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```Potential Energy Bloq```\n",
"## Potential Energy Bloq\n",
"\n",
"Here we consider the electron-electron interaction\n",
"\\begin{align}\n",
Expand Down Expand Up @@ -539,7 +583,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```Comparison to Costs in Paper```\n",
"## Comparison to Costs in Paper\n",
"\n",
"We're now in a position to compare our qualtran costs to those in the paper, which were rough estimates."
]
Expand Down Expand Up @@ -583,20 +627,14 @@
" cirq.ops.common_gates.CXPowGate,\n",
")\n",
"def generalize(bloq):\n",
" if isinstance(bloq, Split):\n",
" return None\n",
" if isinstance(bloq, CirqGateAsBloq):\n",
" if isinstance(bloq.gate, cirq_ft.algos.And) and (len(bloq.gate.cv) == 2):\n",
" return And(cv1=and_cv0, cv2=and_cv1, adjoint=bloq.gate.adjoint)\n",
" if \"QROM\" in bloq.pretty_name():\n",
" if isinstance(bloq.gate, cirq_ft.QROM):\n",
" bloq.gate.__class__.__repr__ = custom_repr\n",
" if isinstance(bloq.gate, cirq_cliffords_ignore):\n",
" return None\n",
" if isinstance(bloq, Join):\n",
" return None\n",
" if isinstance(bloq, Allocate):\n",
" return None\n",
" if isinstance(bloq, Free):\n",
" if isinstance(bloq, (Join, Split, Allocate, Free)):\n",
" return None\n",
" if isinstance(bloq, RotationBloq):\n",
" return attrs.evolve(bloq, angle=phi)\n",
Expand All @@ -609,8 +647,8 @@
"metadata": {},
"outputs": [],
"source": [
"from qualtran.bloqs.chemistry.trotter import PolynmomialEvaluation \n",
"poly_eval = PolynmomialEvaluation(14, 15, 24)\n",
"from qualtran.bloqs.chemistry.trotter import PolynmomialEvaluationInverseSquareRoot \n",
"poly_eval = PolynmomialEvaluationInverseSquareRoot(14, 15, 24)\n",
"graph, sigma = get_bloq_counts_graph(poly_eval, generalizer=generalize)\n",
"GraphvizCounts(graph).get_svg()"
]
Expand All @@ -631,8 +669,8 @@
"metadata": {},
"outputs": [],
"source": [
"from qualtran.bloqs.chemistry.trotter import NewtonRaphson\n",
"nr = NewtonRaphson(14, 15, 24)\n",
"from qualtran.bloqs.chemistry.trotter import NewtonRaphsonApproxInverseSquareRoot\n",
"nr = NewtonRaphsonApproxInverseSquareRoot(14, 15, 24)\n",
"graph, sigma = get_bloq_counts_graph(nr, generalizer=generalize)\n",
"GraphvizCounts(graph).get_svg()\n"
]
Expand Down Expand Up @@ -707,11 +745,11 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## ```Appendix```\n",
"## Appendix\n",
"\n",
"Note in practice we really use fixed point representation of floating numbers.\n",
"\n",
"### ```Fixed point and binary arithmetic```\n",
"### Fixed point and binary arithmetic\n",
"\n",
"To evaluate the polynomial and Newton-Raphson step we need fixed-point arithmetic (addition, multiplication, squaring and scaling). Recall that fixed-point real valued (between 0 and 1) numbers are approximated as (using a convention where the most significant bit is yielded first)\n",
"\n",
Expand Down
Loading

0 comments on commit 388416a

Please sign in to comment.