TwoQubitDecompositionsIntoNativeGates.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize, NonlinearConstraint\n",
    "from scipy.linalg import fractional_matrix_power, expm\n",
    "from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitBasisDecomposer\n",
    "import qiskit as q\n",
    "\n",
    "CNOT_UNITARY = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])\n",
    "ISWAP_UNITARY = np.array([[1, 0, 0, 0], [0, 0, -1j, 0], [0, -1j, 0, 0], [0, 0, 0, 1]])\n",
    "CR90_UNITARY = np.array([[1, -1j, 0, 0], [-1j, 1, 0, 0], [0, 0, 1, 1j], [0, 0, 1j, 1]]) / np.sqrt(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Define Native Two-Qubit Gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NativeTwoQubitGate(object):\n",
    "    num_params = -1\n",
    "\n",
    "    def get_unitary(self, params=[]):\n",
    "        assert len(params) == self.num_params\n",
    "        assert self.num_params != -1, 'Subclass should set num_params'\n",
    "        assert self.num_params in [0, 1], 'Currently doesn\\'t handle multi-parameter native gates'\n",
    "\n",
    "        return self._get_unitary(params)\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        raise NotImplemented('Subclass should implement')\n",
    "        \n",
    "    def __str__(self):\n",
    "        return self.__class__.__name__\n",
    "\n",
    "    \n",
    "class NativeCNOT(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return CNOT_UNITARY\n",
    "\n",
    "\n",
    "class NativeiSWAP(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return ISWAP_UNITARY\n",
    "    \n",
    "\n",
    "class NativeRootiSWAP(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return fractional_matrix_power(ISWAP_UNITARY, 0.5)\n",
    "\n",
    "\n",
    "\n",
    "class NativeCR(NativeTwoQubitGate):\n",
    "    num_params = 1\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        theta = params[0]\n",
    "        return np.array([\n",
    "            [np.cos(theta/2), -1j * np.sin(theta/2), 0, 0],\n",
    "            [-1j * np.sin(theta/2), np.cos(theta/2), 0, 0],\n",
    "            [0, 0, np.cos(theta/2), 1j * np.sin(theta/2)],\n",
    "            [0, 0, 1j * np.sin(theta/2), np.cos(theta/2)]])\n",
    "\n",
    "\n",
    "class NativeCZ(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return np.diag([1,1,1,-1])\n",
    "    \n",
    "    \n",
    "class NativeCR90(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return CR90_UNITARY\n",
    "\n",
    "\n",
    "class NativeParametrizediSWAP(NativeTwoQubitGate):\n",
    "    num_params = 1\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        theta = params[0]\n",
    "        return np.array([\n",
    "            [1, 0, 0, 0],\n",
    "            [0, np.cos(theta), -1j*np.sin(theta), 0],\n",
    "            [0, -1j*np.sin(theta), np.cos(theta), 0],\n",
    "            [0, 0, 0, 1]])\n",
    "    \n",
    "    \n",
    "class NativeSWAPAlpha(NativeTwoQubitGate):\n",
    "    num_params = 1\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        alpha = params[0]\n",
    "        return np.array([\n",
    "            [1, 0, 0, 0],\n",
    "            [0, (1 + np.exp(1j*np.pi*alpha))/2, (1 - np.exp(1j*np.pi*alpha))/2, 0],\n",
    "            [0, (1 - np.exp(1j*np.pi*alpha))/2, (1 + np.exp(1j*np.pi*alpha))/2, 0],\n",
    "            [0, 0, 0, 1]\n",
    "        ])\n",
    "\n",
    "    \n",
    "class NativeBSWAP(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        return np.array([\n",
    "            [np.cos(np.pi/8), 0, 0, 1j * np.sin(np.pi/8)],\n",
    "            [0, np.cos(3*np.pi/8), 1j * np.sin(3*np.pi/8), 0],\n",
    "            [0, 1j * np.sin(3*np.pi/8), np.cos(3*np.pi/8), 0],\n",
    "            [1j * np.sin(np.pi/8), 0, 0, np.cos(np.pi/8)]\n",
    "        ])\n",
    "\n",
    "\n",
    "class NativeMAP(NativeTwoQubitGate):\n",
    "    num_params = 0\n",
    "\n",
    "    def _get_unitary(self, params=[]):\n",
    "        return expm(-1j * np.pi/4 * np.diag([1, -1, -1, 1]))\n",
    "\n",
    "    \n",
    "class NativeRIP(NativeTwoQubitGate):\n",
    "    pass\n",
    "\n",
    "\n",
    "class NativeMS(NativeTwoQubitGate):\n",
    "    num_params = 1\n",
    "    \n",
    "    def _get_unitary(self, params=[]):\n",
    "        theta = params[0]\n",
    "        return np.array([\n",
    "            [np.cos(np.pi*theta/2), 0, 0, -1j * np.sin(np.pi*theta/2)],\n",
    "            [0, np.cos(np.pi*theta/2), -1j*np.sin(np.pi*theta/2), 0],\n",
    "            [0, -1j*np.sin(np.pi*theta/2), np.cos(np.pi*theta/2), 0],\n",
    "            [-1j*np.sin(np.pi*theta/2), 0, 0, np.cos(np.pi*theta/2)]\n",
    "        ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Target Two-Qubit Operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TargetTwoQubitOperation(object):\n",
    "    def get_unitaries(self, params=None):\n",
    "        raise NotImplemented('Subclass should implement')\n",
    "        \n",
    "    def __str__(self):\n",
    "        return self.__class__.__name__\n",
    "\n",
    "\n",
    "class TargetCNOT(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self):\n",
    "        return [CNOT_UNITARY]\n",
    "    \n",
    "\n",
    "class TargetSWAP(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self):\n",
    "        return [np.array([[1, 0, 0, 0],\n",
    "                          [0, 0, 1, 0],\n",
    "                          [0, 1, 0 ,0],\n",
    "                          [0, 0, 0, 1]])]\n",
    "\n",
    "    \n",
    "class TargetZZInteraction(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self, params=None):\n",
    "        if params is None:\n",
    "            params = np.random.random(10) * 2 * np.pi\n",
    "        return [np.diag([1, np.exp(1j * param), np.exp(1j * param), 1]) for param in params]\n",
    "    \n",
    "    \n",
    "class TargetZZSWAP(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self, params=None):\n",
    "        if params is None:\n",
    "            params = np.random.random(10) * 2 * np.pi\n",
    "        unitaries = []\n",
    "        for param in params:\n",
    "            unitaries.append(np.array([[0, 0, 0, 1],\n",
    "                                       [0, np.exp(1j*param), 0, 0],\n",
    "                                       [0, 0, np.exp(1j*param), 0],\n",
    "                                       [1, 0, 0, 0]]))\n",
    "        return unitaries\n",
    "    \n",
    "\n",
    "class TargetFermionicSimulation(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self, params=None):\n",
    "        params = [np.random.random(2) * 2 * np.pi for _ in range(20)]\n",
    "        unitaries = []\n",
    "        for T, V in params:\n",
    "            unitaries.append(np.array([[1, 0, 0, 0],\n",
    "                                       [0, -1j*np.sin(T), np.cos(T), 0],\n",
    "                                       [0, np.cos(T), -1j*np.sin(T), 0],\n",
    "                                       [0, 0, 0, -np.exp(-1j*V)]]))\n",
    "        return unitaries\n",
    "    \n",
    "    \n",
    "class TargetFermionicFourierTransform(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self):\n",
    "        return [np.array([\n",
    "            [1, 0, 0, 0],\n",
    "            [0, 1/np.sqrt(2), 1/np.sqrt(2), 0],\n",
    "            [0, 1/np.sqrt(2), -1/np.sqrt(2), 0],\n",
    "            [0, 0, 0, -1]\n",
    "        ])]\n",
    "    \n",
    "\n",
    "class TargetBogoliubovTransform(TargetTwoQubitOperation):\n",
    "    def get_unitaries(self, params=None):\n",
    "        if params is None:\n",
    "            params = np.random.random(10) * 2 * np.pi\n",
    "\n",
    "        def _get_unitary(expo):  # adapted from Cirq\n",
    "            # --X--S--|iSWAP^expo|--S^1.5--X--\n",
    "            # --------|iSWAP^expo|------------\n",
    "            U = np.kron(np.array([[0, 1], [1, 0]]), np.eye(2))\n",
    "            U = np.kron(fractional_matrix_power(np.diag([1, -1]), 1.5), np.eye(2)) @ U\n",
    "            U = fractional_matrix_power(ISWAP_UNITARY, expo) @ U\n",
    "            U = np.kron(np.array(np.diag([1, -1])), np.eye(2)) @ U\n",
    "            U = np.kron(np.array([[0, 1], [1, 0]]), np.eye(2)) @ U\n",
    "            return U\n",
    "\n",
    "        return [_get_unitary(param) for param in params]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete Native Gates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "native_two_qubit_gates = [NativeCNOT(), NativeCR90(), NativeCZ(), NativeiSWAP(), NativeBSWAP(), NativeMAP(), NativeRootiSWAP()]\n",
    "target_two_qubit_operations = [TargetCNOT(), TargetSWAP(),\n",
    "                               TargetZZInteraction(), TargetZZSWAP(),\n",
    "                               TargetFermionicSimulation(),\n",
    "                               TargetFermionicFourierTransform(), TargetBogoliubovTransform()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With target operation of TargetCNOT, costs are:\n",
      "1 via NativeCNOT\n",
      "1 via NativeCR90\n",
      "1 via NativeCZ\n",
      "2 via NativeiSWAP\n",
      "2 via NativeBSWAP\n",
      "1 via NativeMAP\n",
      "1.0 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetSWAP, costs are:\n",
      "3 via NativeCNOT\n",
      "3 via NativeCR90\n",
      "3 via NativeCZ\n",
      "3 via NativeiSWAP\n",
      "3 via NativeBSWAP\n",
      "3 via NativeMAP\n",
      "1.5 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetZZInteraction, costs are:\n",
      "2 via NativeCNOT\n",
      "2 via NativeCR90\n",
      "2 via NativeCZ\n",
      "2 via NativeiSWAP\n",
      "2 via NativeBSWAP\n",
      "2 via NativeMAP\n",
      "1.0 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetZZSWAP, costs are:\n",
      "3 via NativeCNOT\n",
      "3 via NativeCR90\n",
      "3 via NativeCZ\n",
      "3 via NativeiSWAP\n",
      "3 via NativeBSWAP\n",
      "3 via NativeMAP\n",
      "1.5 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetFermionicSimulation, costs are:\n",
      "3 via NativeCNOT\n",
      "3 via NativeCR90\n",
      "3 via NativeCZ\n",
      "3 via NativeiSWAP\n",
      "3 via NativeBSWAP\n",
      "3 via NativeMAP\n",
      "1.5 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetFermionicFourierTransform, costs are:\n",
      "2 via NativeCNOT\n",
      "2 via NativeCR90\n",
      "2 via NativeCZ\n",
      "2 via NativeiSWAP\n",
      "2 via NativeBSWAP\n",
      "2 via NativeMAP\n",
      "1.0 via NativeRootiSWAP\n",
      "----------------------------\n",
      "With target operation of TargetBogoliubovTransform, costs are:\n",
      "2 via NativeCNOT\n",
      "2 via NativeCR90\n",
      "2 via NativeCZ\n",
      "2 via NativeiSWAP\n",
      "2 via NativeBSWAP\n",
      "2 via NativeMAP\n",
      "1.0 via NativeRootiSWAP\n",
      "----------------------------\n"
     ]
    }
   ],
   "source": [
    "for target_operation in target_two_qubit_operations:\n",
    "    print('With target operation of %s, costs are:' % target_operation)\n",
    "\n",
    "    for native_gate in native_two_qubit_gates:\n",
    "        decomposer = TwoQubitBasisDecomposer(q.extensions.UnitaryGate(native_gate.get_unitary()))\n",
    "        num_basis_gates = 0\n",
    "        for unitary in target_operation.get_unitaries():\n",
    "            num_basis_gates = max(num_basis_gates, decomposer.num_basis_gates(unitary))\n",
    "        if type(native_gate) == NativeRootiSWAP:\n",
    "            num_basis_gates *= 0.5  # sqrt(iSWAP) gate cost should be normalized to iSWAP cost\n",
    "        print(num_basis_gates, 'via', native_gate)\n",
    "    print('----------------------------')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continuous Native Gates:\n",
    "First define minimization problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_best_decomposition(target_two_qubit_operation, native_two_qubit_gate):\n",
    "    for k in [1, 2, 3]:\n",
    "        if has_k_step_decomposition(target_two_qubit_operation, native_two_qubit_gate, k):\n",
    "            return k\n",
    "    assert False, 'did not find any 3-step decomposition'\n",
    "            \n",
    "\n",
    "def has_k_step_decomposition(target_two_qubit_operation, native_two_qubit_gate, k):\n",
    "    for target_unitary in target_two_qubit_operation.get_unitaries():\n",
    "        num_failures = 0\n",
    "        while True:\n",
    "            result = k_step_minimization(target_unitary, native_two_qubit_gate, k)\n",
    "            if result.fun > -0.99:\n",
    "                num_failures += 1  # retry up to 10 times before giving up\n",
    "                if num_failures >= 10:\n",
    "                    return False\n",
    "            else:\n",
    "                break\n",
    "    return True\n",
    "\n",
    "\n",
    "def k_step_minimization(target_unitary, native_two_qubit_gate, k):\n",
    "    # there are k+1 single qubit layers and k two qubit layers\n",
    "    # yielding a total of 2(k+1) single qubit gates and k two qubit gates\n",
    "    # each single qubit gate has 3 parameters, so they contribute to 6k + 6 parameters\n",
    "    num_params = (6 * k + 6) + k * (native_two_qubit_gate.num_params)\n",
    "    params = [np.random.random() for _ in range(num_params)]\n",
    "    \n",
    "    from functools import partial\n",
    "    foo = partial(get_fidelity, target_unitary=target_unitary, native_two_qubit_gate=native_two_qubit_gate, k=k)\n",
    "    nonlinear_constraint= NonlinearConstraint(foo, lb=-np.inf, ub=-0.999);\n",
    "    return minimize(two_qubit_duration, params,\n",
    "                    args=(target_unitary, native_two_qubit_gate, k),\n",
    "                    method='trust-constr', constraints=[nonlinear_constraint],\n",
    "                    options={'maxiter': 1000})\n",
    "\n",
    "\n",
    "def get_fidelity(params, target_unitary, native_two_qubit_gate, k):\n",
    "    params = list(params)  # so that we can call .pop()\n",
    "    single_qubit_u = []  # list of 2(k+1) single qubit gates\n",
    "    for i in range(2*(k + 1)):\n",
    "        a, b, c = params[3*i], params[3*i + 1], params[3*i + 2]  # generic single-qubit unitary has 3 params\n",
    "        single_qubit_u.append(np.array([[np.exp(1j * a) * np.cos(b), np.exp(1j * c) * np.sin(b)],\n",
    "                                        [-np.exp(-1j * c) * np.sin(b), np.exp(-1j * a) * np.cos(b)]]))\n",
    "\n",
    "    two_qubit_u = []  # list of k two-qubit gates\n",
    "    if native_two_qubit_gate.num_params == 0:\n",
    "        for _ in range(k):\n",
    "            two_qubit_u.append(native_two_qubit_gate.get_unitary())\n",
    "    else:\n",
    "        for i in range(6*k+6, 6*k+9):\n",
    "            two_qubit_u.append(native_two_qubit_gate.get_unitary([params[i]]))\n",
    "\n",
    "    actual_unitary = np.eye(4)\n",
    "    for i in range(k):\n",
    "        actual_unitary = actual_unitary @ np.kron(single_qubit_u.pop(), single_qubit_u.pop())\n",
    "        actual_unitary = actual_unitary @ two_qubit_u.pop()\n",
    "    actual_unitary = actual_unitary @ np.kron(single_qubit_u.pop(), single_qubit_u.pop())\n",
    "\n",
    "    fidelity = np.abs(np.trace(actual_unitary @ np.linalg.inv(target_unitary))) / 4\n",
    "    return -fidelity  # return negative, since we will minimize\n",
    "\n",
    "\n",
    "def two_qubit_duration(params, target_unitary, native_two_qubit_gate, k):\n",
    "    two_qubit_params = params[-k:]\n",
    "    return sum(map(np.abs, two_qubit_params)) / (np.pi / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_step_minimization(TargetSWAP().get_unitaries()[0], NativeParametrizediSWAP(), 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(map(np.abs, [-0.77709967,0.78626331,1.56080778]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "native_two_qubit_gates = [NativeCNOT(), NativeiSWAP(), NativeCR(), NativeParametrizediSWAP()]\n",
    "target_two_qubit_operations = [TargetCNOT(), TargetSWAP(),\n",
    "                               TargetZZInteraction(), TargetZZSWAP(),\n",
    "                               TargetFermionicSimulation(),\n",
    "                               TargetFermionicFourierTransform(), TargetBogoliubovTransform()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "for target_two_qubit_operation in target_two_qubit_operations:\n",
    "    for native_two_qubit_gate in native_two_qubit_gates:\n",
    "        k = get_best_decomposition(target_two_qubit_operation, native_two_qubit_gate)\n",
    "        results[(target_two_qubit_operation, native_two_qubit_gate)] = k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_format =\"{:>20}\" * (len(native_two_qubit_gates) + 1)\n",
    "print(row_format.format(\"\", *map(lambda s: str(s)[6:], native_two_qubit_gates)))\n",
    "for target_two_qubit_operation in target_two_qubit_operations:\n",
    "    row = [results[target_two_qubit_operation, native_two_qubit_gate] for native_two_qubit_gate in native_two_qubit_gates]\n",
    "    print(row_format.format(str(target_two_qubit_operation)[6:], *row))"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}