estimation prob workflow

Estimation problem workflow.ipynb

@@ -0,0 +1,315 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Estimation problem workflow\n",
+    "\n",
+    "For a multivariate random variable $X$ and a real-valued function $f$ the goal is to compute \n",
+    "\n",
+    "$$\n",
+    "    \\mathbb{E}[f(X)], \\text{ where } f: \\mathbb{R}^n \\mapsto I \\subset \\mathbb R\n",
+    "$$\n",
+    "\n",
+    "Here, this is achieved using quantum amplitude estimation (QAE), based on a $\\mathcal{A}$ operator in the format\n",
+    "\n",
+    "$$\n",
+    "    \\mathcal{A}|0\\rangle = \\sqrt{1 - a}|\\Psi_0\\rangle + \\sqrt{a}|\\Psi_1\\rangle\n",
+    "      = \\sum_{\\hat x = 0}^{2^n - 1} \\sqrt{p_{\\phi(\\hat x)} (1 - \\hat f(\\phi(\\hat x)))} |x\\rangle |0\\rangle\n",
+    "      + \\sum_{\\hat x = 0}^{2^n - 1} \\sqrt{p_{\\phi(\\hat x)} \\hat f(\\phi(\\hat x))} |x\\rangle |1\\rangle\n",
+    "$$\n",
+    "\n",
+    "To revert the scalings from $f$ to $\\hat f$ and possible post-processing from the function approximation, the result of QAE, $\\tilde a$, is mapped to the result via a post-processing function $\\eta$\n",
+    "\n",
+    "$$\n",
+    "\\mathbb{E}[f(X)] = \\eta(a) \\approx \\eta(\\tilde a)\n",
+    "$$\n",
+    "\n",
+    "**Notes:**\n",
+    "\n",
+    "An amplitude function $F$ of a function $\\hat f$ is a mapping\n",
+    "\n",
+    "$$\n",
+    "F: |x\\rangle|0\\rangle \\mapsto \\sqrt{1 - \\hat f(x)}|x\\rangle|0\\rangle + \\sqrt{\\hat f(x)}|x\\rangle|1\\rangle\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$$\n",
+    "\\hat f: \\mathbb{N}_0 \\rightarrow [0, 1].\n",
+    "$$\n",
+    "\n",
+    "If we have a function $f$ that does not act on these spaces we must normalize the image and apply an affine transformation from the input domain to $\\mathbb{N}_0$. Let the function $f$ be defined on $2^n$ equidistant points in $[a, b]$:\n",
+    "\n",
+    "$$\n",
+    "f: [a, b] \\rightarrow [c, d]\n",
+    "$$\n",
+    "\n",
+    "then the rescaled version is\n",
+    "\n",
+    "$$\n",
+    "\\hat f(x) = \\frac{f(\\phi(x)) - c}{d - c}\n",
+    "$$\n",
+    "\n",
+    "where $\\phi(x) = a + (b - a) x / 2^n$ is the affine transformation.\n",
+    "\n",
+    "The normalized function $\\hat f$ can be implemented in different manners. In general this requires a post-processing step, $\\zeta$:\n",
+    "\n",
+    "$$\n",
+    "\\hat f \\approx \\zeta(T\\hat f)\n",
+    "$$\n",
+    "\n",
+    "where $T\\hat f$ is the amplitude we can map onto the qubits.\n",
+    "One example is a Taylor approximation of $\\sin^2$ which can be mapped to qubits using RY gates.\n",
+    "\n",
+    "Examples:\n",
+    "\n",
+    "* quadratic terms can be approximated with\n",
+    "    $$\n",
+    "    ax^2 \\approx \\sin^2(c\\sqrt{a}x) / c^2\n",
+    "    $$\n",
+    "    where the sine part is estimated with amplitude estimation and the rescaling is $\\zeta(x) = x/c^2$\n",
+    "    \n",
+    "* linear terms can be approximated as\n",
+    "    $$\n",
+    "    ax \\approx \\zeta\\left( \\sin^2\\left(\\frac{\\pi}{4} + \\frac{\\pi c}{2} \\left(f(x) - \\frac{1}{2}\\right)\\right) \\right) \n",
+    "    $$\n",
+    "    with\n",
+    "    $$\n",
+    "    \\zeta(x) = \\frac{2}{\\pi c}\\left(x - \\frac{1}{2}\\right) + \\frac{1}{2}\n",
+    "    $$\n",
+    "\n",
+    "Let $\\tilde a$ be the output of amplitude estimation, then in general we have to apply a post-processing in the form of\n",
+    "\n",
+    "$$\n",
+    "c + (d - c) \\zeta(\\tilde a)\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Defining the random variable\n",
+    "\n",
+    "The distribution of the random variable is loaded from the circuit library. The function will be evaluated on the x-values ``[0, ..., 2 ** num_qubits - 1]``. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit.circuit.library import NormalDistribution\n",
+    "X = NormalDistribution(num_qubits=3, mu=1, sigma=0.5, bounds=(0, 2))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Defining the objective function"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Workflow #1\n",
+    "\n",
+    "Similar to the oracle compiler: take a Python function and synthesize the circuit automatically."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def classical_f(x):\n",
+    "    if x > 1:\n",
+    "        return x ** 2\n",
+    "    return x\n",
+    "\n",
+    "\n",
+    "from qiskit.circuit.library import AmplitudeFunction\n",
+    "f = AmplitudeFunction(classical_f, domain=(0, 2))  # domain must be the same as the bounds of the probability dist\n",
+    "post_processing = f.post_processing  # c + (d - c) zeta"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Workflow #2\n",
+    "\n",
+    "Similar to the `UnivariatePiecewiseLinearObjective`, specify the function and construct the circuit from this information. Includes rescalings."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit.circuit.library import LinearAmplitudeFunction\n",
+    "f = LinearAmplitudeFunction(num_state_qubits=3, \n",
+    "                            slope=[-1, 1],\n",
+    "                            offset=[1, 0],\n",
+    "                            breakpoints=[0, 1],\n",
+    "                            domain=(0, 2))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Workflow #3\n",
+    "\n",
+    "Manual."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit.circuit.library import PiecewiseLinearPauliRotations\n",
+    "\n",
+    "slopes = np.array([-1, 1])\n",
+    "offsets = np.array([1, 0])\n",
+    "breakpoints = np.array([0, 1])\n",
+    "domain = (0, 2)\n",
+    "num_state_qubits = 3\n",
+    "N = 2 ** num_state_qubits\n",
+    "\n",
+    "# apply rescaling to the right domain \n",
+    "phi_inv = lambda x: (x - domain[0]) / (domain[1] - domain[0]) * N\n",
+    "breakpoints = phi_inv(breakpoints)\n",
+    "slopes = (domain[1] - domain[0]) / (N - 1) * slopes\n",
+    "# offsets remain the same\n",
+    "\n",
+    "# apply normalization of the image space\n",
+    "offsets = (offsets - fmin) / (fmax - fmin)\n",
+    "slopes = slopes / (fmax - fmin)\n",
+    "\n",
+    "# apply taylor \n",
+    "offsets *= c * np.pi / 2\n",
+    "slopes *= c * np.pi / 2\n",
+    "\n",
+    "f = PiecewiseLinearPauliRotations(num_state_qubits, \n",
+    "                                  breakpoints=breakpoints,\n",
+    "                                  slopes=2 * slopes,  \n",
+    "                                  offsets=2 * offsets)  "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Conclusion\n",
+    "\n",
+    "Option #2 makes sense for now. This means porting most of the `UnivariatePiecewiseLinearObjective` to the circuit library."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Defining the A operator\n",
+    "\n",
+    "### Workflow #1\n",
+    "\n",
+    "Stack together the function circuit and probability distribution circuit."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit.circuit import QuantumCircuit\n",
+    "\n",
+    "A = QuantumCircuit(f.num_qubits)\n",
+    "A.compose(X, inplace=True)\n",
+    "A.compose(f, inplace=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Workflow #2\n",
+    "\n",
+    "Can introduce a class for that but I don't really think the work it does justifies a class, see rather the section about encapsulating the entire workflow."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit.circuit.library import EstimationProblem  # or aqua algorithms\n",
+    "\n",
+    "A = EstimationProblem(f, X)  # but here f is not the circuit but a function! and we do all the synthesisa"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Run QAE\n",
+    "\n",
+    "QAE in the basic form takes the $\\mathcal{A}$ operator and a post processing function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import Aer\n",
+    "from qiskit.aqua.algorithms import AmplitudeEstimation\n",
+    "\n",
+    "ae = AmplitudeEstimation(num_eval_qubits=4, state_in=A, post_processing=f.post_processing)\n",
+    "result = ae.run(Aer.get_backend('qasm_simulator'))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Encapsulating the workflow\n",
+    "\n",
+    "Add a class that is fully equivalent to `UnivariatePiecewiseLineaObjective`? Essentially this means accepting the probabilit"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.7.7 64-bit ('py377': conda)",
+   "language": "python",
+   "name": "python37764bitpy377condada92b85e600c45a79ff7c5f35e6f13ab"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}