2D Hidden Linear Function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hidden Linear Function Problem\n",
    "\n",
    "In this notebook we consider a problem from paper [1] and build a quantum cirquit, which solves it, in Cirq. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The problem\n",
    "\n",
    "Consider $A \\in \\mathbb{F}_2^{n \\times n}$ - upper-triangular binary matrix, $b \\in \\mathbb{F}_2^n$ - binary vector.\n",
    "\n",
    "Define a function $q : \\mathbb{F}_2^n \\to \\mathbb{Z}_4 $:\n",
    "\n",
    "$$q(x) = (2 x^T A x + b^T x) ~\\text{mod}~ 4, $$ \n",
    "\n",
    "Also define\n",
    "\n",
    "$$\\mathcal{L}_q = \\Big\\{x \\in  \\mathbb{F}_2^n : q(x \\oplus y) = (q(x) + q(y)) ~\\text{mod}~ 4 ~~ \\forall y \\in \\mathbb{F}_2^n \\Big\\}.$$\n",
    "\n",
    "Turns out that restriction of $q$ on $\\mathcal{L}_q$ is a linear function, i.e. there exists such $z \\in \\mathbb{F}_2^n$, that\n",
    "\n",
    "$$q(x) = 2 z^T x \\forall x \\in \\mathcal{L}_q.$$\n",
    "\n",
    "Our task is, given $A$ and $b$, to find $z$. There may be multiple answers - we need to find any such answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation and bruteforce solution\n",
    "\n",
    "For small values of $n$ we can solve this problem with a trivial bruteforce solution. First, we need to build $\\mathcal{L}_q$ by checking for all $2^n$ binary vectors, which of them belongs to it(by definition). Then we need to try all possible $z \\in \\mathbb{F}_2^n$, and for each of them and for each $x \\in \\mathcal{L}_q$ check whether $q(x) = 2 z^T x$.\n",
    "\n",
    "Below we implement a class which represents instance of a problem and solves it with a bruteforce solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cirq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HiddenLinearFunctionProblem:\n",
    "    def __init__(self, A, b):\n",
    "        self.n = A.shape[0]\n",
    "        assert A.shape == (self.n, self.n)\n",
    "        assert b.shape == (self.n, )\n",
    "        for i in range(self.n):\n",
    "            for j in range(i+1):\n",
    "                assert A[i][j] == 0, 'A[i][j] can be 1 only if i<j'\n",
    "        \n",
    "        self.A = A\n",
    "        self.b = b\n",
    "        \n",
    "    def q(self, x):\n",
    "        assert x.shape == (self.n, )\n",
    "        return (2 * (x @ self.A @ x) + (self.b @ x)) % 4\n",
    "        \n",
    "    def bruteforce_solve(self):\n",
    "        all_vectors = [np.array([(m>>i)%2 for i in range(self.n)]) for m in range(2**self.n)]\n",
    "\n",
    "        def vector_in_L(x):\n",
    "            for y in all_vectors:\n",
    "                if self.q( (x + y)%2 ) != (self.q(x) + self.q(y))%4:\n",
    "                    return False\n",
    "            return True\n",
    "\n",
    "        self.L = [x for x in all_vectors if vector_in_L(x)]\n",
    "\n",
    "        self.all_zs = [z for z in all_vectors if self.is_z(z)]\n",
    "    \n",
    "    # Whether given vector z is solution to this problem.\n",
    "    def is_z(self, z):\n",
    "        assert z.shape == (self.n, )\n",
    "        assert self.L is not None\n",
    "        for x in self.L:\n",
    "            if self.q(x) != 2 * ((z @ x) % 2):\n",
    "                return False\n",
    "        return True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For testing, we need to generate an instance of a problem. We can generate random $A$ and $b$. However, for some $A$ and $b$ problem is trivial - that is, $\\mathcal{L}_q = \\{0\\}$ and therefore any $z$ is a solution. In fact, product of $|\\mathcal{L}_q|$ and number of solutions is always equal to $2^n$ (see prrof in [1]), so we want a problem with large $\\mathcal{L}_q$.\n",
    "\n",
    "Code below can be used to generate random problem with given size of $\\mathcal{L}_q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 256\n"
     ]
    }
   ],
   "source": [
    "def random_problem(n, seed=None):\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed) \n",
    "    A = np.random.randint(0, 2, size=(n,n))\n",
    "    for i in range(n):\n",
    "        for j in range(i+1):\n",
    "            A[i][j] = 0\n",
    "    b = np.random.randint(0, 2, size=n)\n",
    "    problem = HiddenLinearFunctionProblem(A, b)\n",
    "    return problem\n",
    "        \n",
    "def find_interesting_problem(n, min_L_size):\n",
    "    for _ in range(1000):\n",
    "        problem = random_problem(n)\n",
    "        problem.bruteforce_solve()\n",
    "        if len(problem.L) >= min_L_size and not np.max(problem.A) == 0:\n",
    "            return problem\n",
    "    return None\n",
    "\n",
    "problem = find_interesting_problem(10, 4)\n",
    "print(len(problem.L), len(problem.all_zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We found a problem with $n=10$ and $|\\mathcal{L}_q|=16$, so only 64 of 1024 possible vectors are solutions. So, chance of randomly guessing a solution is $\\frac{1}{16}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16 64\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[0, 1, 1, 0, 0, 1, 0, 0, 1, 1],\n",
    "              [0, 0, 0, 1, 1, 1, 1, 1, 1, 1],\n",
    "              [0, 0, 0, 0, 0, 0, 1, 1, 0, 1],\n",
    "              [0, 0, 0, 0, 0, 0, 0, 0, 1, 0],\n",
    "              [0, 0, 0, 0, 0, 1, 0, 0, 0, 1],\n",
    "              [0, 0, 0, 0, 0, 0, 1, 1, 0, 0],\n",
    "              [0, 0, 0, 0, 0, 0, 0, 0, 0, 1],\n",
    "              [0, 0, 0, 0, 0, 0, 0, 0, 1, 1],\n",
    "              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
    "              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])\n",
    "b = np.array([0, 0, 0, 0, 1, 1, 1, 0, 0, 1])\n",
    "problem_10_64 = HiddenLinearFunctionProblem(A, b)\n",
    "problem_10_64.bruteforce_solve()\n",
    "print(len(problem_10_64.L), len(problem_10_64.all_zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution with a quantum circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in [1], given problem can be solved by a quantum circuit, which implements operator $H ^ {\\otimes n} U_q H ^ {\\otimes n}$, where\n",
    "\n",
    "$$U_q = \\prod_{1 < i < j < n} CZ_{ij}^{A_{ij}} \\cdot \\bigotimes_{j=1}^{n} S_j^{b_j} .$$\n",
    "\n",
    "We need to apply this operator to $| 0^n \\rangle$ and measure the result - result is guaranteed to be one of the solutions. Moreover, we can get any solution with equal probability.\n",
    "\n",
    "Let's implement code which would generate such circuit and simulate it.\n",
    "\n",
    "Note that: \n",
    "\n",
    "* We use Cirq S gate, whose matrix is $\\left(\\begin{smallmatrix}1 & 0\\\\0 & i\\end{smallmatrix}\\right)$. In the paper [1] matrix of S gate is defined as $\\left(\\begin{smallmatrix}1 & 0\\\\0 & -i\\end{smallmatrix}\\right)$. But for this problem it doesn't matter.\n",
    "\n",
    "* We reorder CZ gates in such a way so they take less moments. This is a problem of minimal [edge coloring](https://en.wikipedia.org/wiki/Edge_coloring), and we solve it here with a simple greedy algorithm. We can do that because CZ gates commute (because their matrices are diagonal).\n",
    "\n",
    "* All gates are Clifford gates, so we can use Clifford simulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          ┌───┐   ┌──┐   ┌───┐   ┌───┐   ┌───┐\n",
      "0: ───H────@───────@──────@───────@───────@──────────────────────H───M───\n",
      "           │       │      │       │       │\n",
      "1: ───H────@───────┼@─────┼@──────┼@──────┼@─────@───@───@───────H───M───\n",
      "                   ││     ││      ││      ││     │   │   │\n",
      "2: ───H────@───────@┼─────┼┼@─────┼┼@─────┼┼─────┼───┼───┼───────H───M───\n",
      "           │        │     │││     │││     ││     │   │   │\n",
      "3: ───H────┼@───────@─────┼┼┼─────┼┼┼─────┼┼─────┼───┼───┼───────H───M───\n",
      "           ││             │││     │││     ││     │   │   │\n",
      "4: ───H────┼┼@─────@──────┼@┼─────┼┼┼─────┼┼─────┼───┼───┼───S───H───M───\n",
      "           │││     │      │ │     │││     ││     │   │   │\n",
      "5: ───H────┼┼@─────┼@─────@─┼─────┼@┼─────┼┼@────┼───┼───┼───S───H───M───\n",
      "           ││      ││       │     │ │     │││    │   │   │\n",
      "6: ───H────@┼──────┼@─────@─┼─────┼─┼─────┼@┼────┼───┼───┼───S───H───M───\n",
      "            │      │      │ │     │ │     │ │    │   │   │\n",
      "7: ───H────@┼──────┼@─────┼─@─────┼─┼─────┼─@────@───┼───┼───────H───M───\n",
      "           ││      ││     │       │ │     │          │   │\n",
      "8: ───H────┼@──────┼@─────┼───────@─┼─────┼──────────@───┼───────H───M───\n",
      "           │       │      │         │     │              │\n",
      "9: ───H────@───────@──────@─────────@─────@──────────────@───S───H───M───\n",
      "          └───┘   └──┘   └───┘   └───┘   └───┘\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1, 0, 1, 1, 1, 0, 0, 1, 0, 1], dtype=uint8)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Given adjacency matrix A, returns list of lists of edges, such as \n",
    "# edges in each list do not have common vertex.\n",
    "def graph_coloring(A):\n",
    "    A = np.copy(A)\n",
    "    n = A.shape[0]\n",
    "    ans = []\n",
    "    while np.max(A) != 0:\n",
    "        edges_group = []\n",
    "        used = np.zeros(n, dtype=np.bool)\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                if A[i][j] == 1 and not used[i] and not used[j]:\n",
    "                    edges_group.append((i, j))\n",
    "                    A[i][j] = 0\n",
    "                    used[i] = used[j] = True\n",
    "        ans.append(edges_group)\n",
    "    return ans\n",
    "    \n",
    "def generate_circuit_for_problem(problem):\n",
    "    qubits = cirq.LineQubit.range(problem.n)\n",
    "    circuit = cirq.Circuit()\n",
    "    \n",
    "    # Hadamard gates at the beginning.\n",
    "    circuit += cirq.Moment([cirq.H(qubits[i]) for i in range(problem.n)])\n",
    "    \n",
    "    # Controlled-Z gates encoding the matrix A.\n",
    "    for layer in graph_coloring(problem.A):\n",
    "        for i, j in layer:\n",
    "            circuit += cirq.CZ(qubits[i], qubits[j])\n",
    "        \n",
    "    # S gates encoding the vector b.\n",
    "    S_moment = cirq.Moment()\n",
    "    for i in range(problem.n):\n",
    "        if problem.b[i] == 1:\n",
    "            S_moment += cirq.S.on(qubits[i])\n",
    "    circuit += S_moment\n",
    "            \n",
    "    # Hadamard gates at the end.\n",
    "    circuit += cirq.Moment([cirq.H(qubits[i]) for i in range(problem.n)])\n",
    "    \n",
    "    # Measurements.\n",
    "    circuit += cirq.Moment([cirq.measure(qubits[i], key=str(i)) for i in range(problem.n)]) \n",
    "    \n",
    "    return circuit\n",
    "\n",
    "def solve_problem(problem, print_circuit=False):\n",
    "    circuit = generate_circuit_for_problem(problem)\n",
    "        \n",
    "    if print_circuit:\n",
    "        print(circuit)\n",
    "    \n",
    "    sim = cirq.CliffordSimulator()\n",
    "    result = sim.simulate(circuit)\n",
    "    z = np.array([result.measurements[str(i)][0] for i in range(problem.n)])\n",
    "    return z\n",
    "\n",
    "solve_problem(problem_10_64, print_circuit=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's test this algorithm. Let's solve it with a circuit 100 times and each time check that measurement result is indeed an answer to the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OK\n"
     ]
    }
   ],
   "source": [
    "def test_problem(problem):\n",
    "    problem.bruteforce_solve()\n",
    "    for _ in range(100):\n",
    "        z = solve_problem(problem)\n",
    "        assert problem.is_z(z)\n",
    "    \n",
    "test_problem(problem_10_64)\n",
    "print('OK')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat that for 10 other problems with $n=8$ and chance of random guessing at most $\\frac{1}{4}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OK\n"
     ]
    }
   ],
   "source": [
    "for _ in range(10):\n",
    "    test_problem(find_interesting_problem(8, 4))\n",
    "print('OK')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's run our algorithm on a problem with $n=200$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 9.26 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1,\n",
       "       0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1,\n",
       "       1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0,\n",
       "       0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0,\n",
       "       1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0,\n",
       "       1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,\n",
       "       1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1,\n",
       "       1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1,\n",
       "       1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0,\n",
       "       1, 1], dtype=uint8)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "problem = random_problem(200, seed=0)\n",
    "solve_problem(problem, print_circuit=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why is this problem interesting?\n",
    "\n",
    "### 1. It's a problem without an oracle\n",
    "\n",
    "This problem is similar to a problem solved by [Bernstein–Vazirani algorithm](https://en.wikipedia.org/wiki/Bernstein%E2%80%93Vazirani_algorithm). It also finds coefficients of unknown linear function. But in Bernstein-Vazirani algorithm this function is represented by an oracle. In this problem, the linear function is \"hidden\" in inputs $A$ and $b$.\n",
    "\n",
    "### 2. Quantum circuits have advantage over classical when solving this problem\n",
    "\n",
    "According to [Gottesman–Knill theorem](https://en.wikipedia.org/wiki/Gottesman%E2%80%93Knill_theorem), this problem can be solved in polynomial time on classical computer, because it can be solved by simulating Clifford circuit. So, it might look like quantum comuters aren't better than classical ones in solving this problem.\n",
    "\n",
    "However, if we apply certain restrictions on matrix $A$, the circuit will have fixed depth (i.e. number of Moments). Namely, if the matrix $A$ is an adjacency matrix of a \"grid\" graph (whose edges can be colored in 4 colors), all CZ gates will fit in 4 moments, and overall we will have only 8 moments - and this doesn't depend on $n$.\n",
    "\n",
    "But for classical circuits it can be proven (see [1]) that even if we restrict matrix $A$ in the same way, the depth of classical circuit (with gates of bounded fan-in) must grow as $n$ grows (in fact, it grows as $\\log(n)$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "\n",
    "[1] [Quantum advantage with shallow circuits](https://arxiv.org/pdf/1704.00690.pdf) by Sergey Bravyi, David Gosset and Robert König."
   ]
  }
 ],
 "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}