Simon's Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import numpy as np\n",
    "import pyquil.api as api\n",
    "from pyquil.gates import *\n",
    "from pyquil.quil import Program\n",
    "from gaussian_elimination import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Problem Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The setup for Simon's problem consists of a given black-box operator that is a generalization from those given in the Deutsch and Deutsch-Jozsa problems, and maps $\\mathbf{f}: \\{0, 1\\}^n \\rightarrow \\{0, 1\\}^m$, such that<br>\n",
    "<br>\n",
    "$$ U_f : \\left\\vert \\mathbf{x} \\right\\rangle \\left\\vert \\mathbf{b} \\right\\rangle \\rightarrow \\left\\vert \\mathbf{x} \\right\\rangle \\left\\vert \\mathbf{b} \\oplus \\mathbf{f}(\\mathbf{x}) \\right\\rangle$$\n",
    "<br>\n",
    "where $\\mathbf{f}(\\mathbf{x}) \\in \\{0, 1\\}^m \\, \\, \\forall \\mathbf{x} \\in \\{ 0, 1 \\}^n$, $\\mathbf{b} \\in \\{0, 1\\}^m$, and the $\\oplus$ sign represents mod 2 addition on each of the components separately. The problem consists of finding \n",
    "$\\mathbf{s} \\in \\{0, 1\\}^n$ such that<br>\n",
    "<br>\n",
    "$$\\mathbf{f} (\\mathbf{x} \\oplus \\mathbf{s}) = \\mathbf{f} (\\mathbf{x})$$\n",
    "so that the function $\\mathbf{f}$ is periodic with period $\\mathbf{s}$.\n",
    "\n",
    "We solve by first preparing the state $\\left\\vert\\mathbf{x} \\right\\rangle \\left\\vert 0 \\right\\rangle$, applying the black-box to produce the state $\\left\\vert\\mathbf{x}\\right\\rangle \\left\\vert \\mathbf{f}(\\mathbf{x})\\right\\rangle$, then applying $H^{\\otimes n}$ to the first register $\\left\\vert\\mathbf{x}\\right\\rangle$, then measuring it and recording the value $\\mathbf{w}_i$, repeating these steps until $\\text{span}\\{\\mathbf{w}_i\\}$ equals $n-1$, at which point we solve the equation $\\mathbf{W}\\mathbf{s}^{T} = \\mathbf{0}^{T}$ via Gaussian elimination to obtain $\\mathbf{s}$ as the unique non-zero solution. To see _why_ this works, the reader is referred to \"An introduction to quantum computing\" by P. Kaye et al."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Implementation Notes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generalize the black-box operator from the Deutsch-Jozsa problem to construct the one required here\n",
    "$$U_f = \\sum_{\\mathbf{x}=0}^{2^{n} - 1} \\left\\vert \\mathbf{x} \\right\\rangle \\left\\langle \\mathbf{x} \\right\\vert \\otimes \\left[ I + f_{i} (\\mathbf{x}) \\left( X - I \\right) \\right]^{\\otimes_{i=m-1}^{i=0}}$$\n",
    "For example, if $m=2$, then\n",
    "$$ \\left[ I + f_{i} (\\mathbf{x}) \\left( X - I \\right) \\right]^{\\otimes_{i=m-1}^{i=0}} = \\left[ I + f_1(\\mathbf{x}) \\left( X - I \\right) \\right] \\otimes \\left[ I + f_0(\\mathbf{x}) \\left( X - I \\right)\\right]$$\n",
    "<br>\n",
    "and further if $n=3$, $\\mathbf{x} = 010$, and $\\mathbf{f}(\\mathbf{x}) = 10$, then\n",
    "$$ \\left[ I + f_{i} (\\mathbf{x}) \\left( X - I \\right) \\right]^{\\otimes_{i=m-1}^{i=0}} = \\left[ I + f_1(010) \\left( X - I\\right)\\right] \\otimes \\left[ I + f_0(010) \\left( X - I\\right)\\right] \\\\\n",
    "= \\left[ I + (1)(X-I)\\right] \\otimes \\left[ I + (0) (X-I)\\right] \\\\\n",
    "= X \\otimes I$$\n",
    "<br>\n",
    "The sampling of the $\\mathbf{w}_{i}$ is done in such a way as to ensure the reduced row-echelon form of the collective $\\mathbf{W}$ matrix (note that since we're working with mod 2 arithmetic, we automatically have reduced row-echelon, and not just row-echelon form). Back-substitution is modified to work with mod 2 arithmetic. The entire process is implemented in gaussian_elimination.py, and for an excellent discussion of the mathematical details involved, the reader is referred to Section 18.13 of \"<c|Q|c> : A Course in Quantum Computing (for the Community College)\", Vol. 1  by Michael Loceff."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simon's Algorithm using (n+m) qubits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubit_strings(n):\n",
    "    qubit_strings = []\n",
    "    for q in itertools.product(['0', '1'], repeat=n):\n",
    "        qubit_strings.append(''.join(q))\n",
    "    return qubit_strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def black_box_map(n, m, s):\n",
    "    \"\"\"\n",
    "    Black-box map f:{0,1}^n -> {0,1}^m, randomly taking values,\n",
    "    and periodic with period s\n",
    "    \"\"\"\n",
    "    # ensure s lives in {0,1}^n\n",
    "    if len(s) != n:\n",
    "        raise AssertionError(\"Length of period vector should equal n\")\n",
    "    # control qubits\n",
    "    cont_qubs = qubit_strings(n)\n",
    "    # target qubits\n",
    "    targ_qubs = qubit_strings(m)\n",
    "    \n",
    "    # initialize empty dictionary to store map values\n",
    "    d_blackbox = {}\n",
    "    # initialize counter over control qubits\n",
    "    i = 0\n",
    "    # randomly select values from {0,1}^m for the periodic function\n",
    "    while set(cont_qubs) - set(d_blackbox.keys()) != set():\n",
    "        # pick a random target\n",
    "        rand_targ = np.random.choice(targ_qubs)\n",
    "        # set the same value for x and x + s\n",
    "        d_blackbox[cont_qubs[i]] = rand_targ\n",
    "        d_blackbox[add_vec_mod2(cont_qubs[i], s)] = rand_targ\n",
    "        # avoid iterating over keys already assigned values\n",
    "        while cont_qubs[i] in d_blackbox.keys():\n",
    "            i = i + 1\n",
    "            if i >= n:\n",
    "                break\n",
    "        \n",
    "    return d_blackbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qubit_ket(qub_string):\n",
    "    \"\"\"\n",
    "    Form a basis ket out of n-bit string specified by the input 'qub_string', e.g.\n",
    "    '001' -> |001>\n",
    "    \"\"\"\n",
    "    e0 = np.array([[1], [0]])\n",
    "    e1 = np.array([[0], [1]])\n",
    "    d_qubstring = {'0': e0, '1': e1}\n",
    "\n",
    "    # initialize ket\n",
    "    ket = d_qubstring[qub_string[0]]\n",
    "    for i in range(1, len(qub_string)):\n",
    "        ket = np.kron(ket, d_qubstring[qub_string[i]])\n",
    "    \n",
    "    return ket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projection_op(qub_string):\n",
    "    \"\"\"\n",
    "    Creates a projection operator out of the basis element specified by 'qub_string', e.g.\n",
    "    '101' -> |101> <101|\n",
    "    \"\"\"\n",
    "    ket = qubit_ket(qub_string)\n",
    "    bra = np.transpose(ket)  # all entries real, so no complex conjugation necessary\n",
    "    proj = np.kron(ket, bra)\n",
    "    return proj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def black_box(n, m, s):\n",
    "    \"\"\"\n",
    "    Inputs:-\n",
    "    n: no. of control qubits\n",
    "    m: no. of target qubits\n",
    "    s: bit-string equal to the period of the black-box map\n",
    "    \n",
    "    Output:-\n",
    "    Unitary representation of the black-box operator\n",
    "    \"\"\"\n",
    "    d_bb = black_box_map(n, m, s)\n",
    "    # initialize unitary matrix\n",
    "    N = 2**(n+m)\n",
    "    unitary_rep = np.zeros(shape=(N, N))\n",
    "    # populate unitary matrix\n",
    "    for k, v in d_bb.items():\n",
    "        # initialize target qubit operator\n",
    "        targ_op = np.eye(2) + int(v[0])*(-np.eye(2) + np.array([[0, 1], [1, 0]]))\n",
    "        # fill out the rest of the target qubit operator\n",
    "        for i in range(1, m):\n",
    "            cont_op = np.eye(2) + int(v[i])*(-np.eye(2) + np.array([[0, 1], [1, 0]]))\n",
    "            targ_op = np.kron(targ_op, cont_op)\n",
    "        # complete the unitary operator for current control qubit-register\n",
    "        unitary_rep += np.kron(projection_op(k), targ_op)\n",
    "    \n",
    "    return unitary_rep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The period vector is found to be:  [1, 0, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "qvm = api.QVMConnection()\n",
    "# pick number of control qubits to be used\n",
    "n = 4\n",
    "# pick number of target qubits to be used\n",
    "m = 2\n",
    "# specify the period as an n bit-string\n",
    "s = '1011'\n",
    "# make sure s has the correct length\n",
    "if len(s) != n:\n",
    "    raise ValueError(\"s does not have correct bit-string length\")\n",
    "# make sure s is non-zero\n",
    "if s == '0' * n:\n",
    "    raise ValueError(\"s should not be zero vector\")\n",
    "# create the unitary black_box operator\n",
    "blackbox = black_box(n, m, s)\n",
    "# initialize the augmented matrix to be solved via Gaussian elimination\n",
    "W = []\n",
    "# initialize counter\n",
    "counter = 0\n",
    "# run main loop\n",
    "while rank(W) < n-1:\n",
    "    # initialize the program\n",
    "    p = Program()\n",
    "\n",
    "    # Define U_f\n",
    "    p.defgate(\"U_f\", blackbox)\n",
    "\n",
    "    # Prepare the initial state (1/sqrt[2])*(|0> + |1>)^(\\otimes n) \\otimes |0>^(\\otimes m)\n",
    "    for m_ in range(m):\n",
    "        p.inst(I(m_))\n",
    "    for n_ in range(m, n+m):\n",
    "        p.inst(H(n_))\n",
    "\n",
    "    # Apply U_f\n",
    "    p.inst((\"U_f\",) + tuple(range(n+m)[::-1]))\n",
    "\n",
    "    # Apply final H^(\\otimes n)\n",
    "    for n_ in range(m, n+m):\n",
    "        p.inst(H(n_))\n",
    "\n",
    "    # Final measurement\n",
    "    classical_regs = list(range(n))\n",
    "    for i, n_ in enumerate(list(range(m, n+m))[::-1]):\n",
    "        p.measure(n_, classical_regs[i])\n",
    "\n",
    "    measure_n_qubits = qvm.run(p, classical_regs)\n",
    "\n",
    "    # flatten out list\n",
    "    z = [item for sublist in measure_n_qubits for item in sublist]\n",
    "    z.append(0)\n",
    "    \n",
    "    # add (or not) the new sample z to W\n",
    "    W = new_sample(W, z)\n",
    "            \n",
    "    # increment counter\n",
    "    counter = counter + 1\n",
    "    del p\n",
    "    \n",
    "print (\"The period vector is found to be: \", solve_reduced_row_echelon_form(W))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}