{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Benchmark Pauli operator projection\"\n",
        "description: \"Benchmark Pauli operator projection for the latest version of Sample-based quantum diagonalization (SQD)\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "55900f6c-2fb6-4d2c-8745-29bad8b66d9f",
      "metadata": {},
      "source": [
        "# Benchmark Pauli operator projection\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5caa5b4",
      "metadata": {},
      "source": [
        "### Pauli string action on a computational basis state\n",
        "\n",
        "The action of a Pauli string on a computational basis state is rather trivial, and is a single computational basis state itself. This is a direct consequence of the structure of Pauli matrices, which only have a single nonzero element on each of their rows. Consequently, their action on a qubit is:\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_x |0 \\rangle = |1 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_x |1 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_y |0 \\rangle = i|1 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_y |1 \\rangle = -i|0 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_z |0 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_z |1 \\rangle = -|1 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "I |0 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "I |1 \\rangle = |1 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "Each bit on the bitstring labeling the computational basis will be labeled by $x \\in \\{0, 1 \\}$. In order to keep the implementation at light as possible, we will\n",
        "represent the bitstrings with `bool` variables: $0\\rightarrow \\textrm{False}$ and\n",
        "$1\\rightarrow \\textrm{True}$.\n",
        "\n",
        "To represent the action of each Pauli operator in a computational basis state,\n",
        "we will assign three variables to it: `diag`, `sign`, `imag`.\n",
        "\n",
        "* `diag` labels whether the operator is diagonal:\n",
        "\n",
        "  * $\\textrm{diag}(I) = \\textrm{True}$\n",
        "  * $\\textrm{diag}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{diag}(\\sigma_y) = \\textrm{False}$\n",
        "  * $\\textrm{diag}(\\sigma_z) = \\textrm{True}$\n",
        "\n",
        "* `sign` Identifies if there is a sign change in the matrix element connected\n",
        "  to either 0 or 1:\n",
        "\n",
        "  * $\\textrm{sign}(I) = \\textrm{False}$\n",
        "  * $\\textrm{sign}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{sign}(\\sigma_y) = \\textrm{True}$\n",
        "  * $\\textrm{sign}(\\sigma_z) = \\textrm{True}$\n",
        "\n",
        "* `imag` Identifies if there is a complex component to the matrix element:\n",
        "\n",
        "  * $\\textrm{imag}(I) = \\textrm{False}$\n",
        "  * $\\textrm{imag}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{imag}(\\sigma_y) = \\textrm{True}$\n",
        "  * $\\textrm{imag}(\\sigma_z) = \\textrm{False}$\n",
        "\n",
        "We label an arbitrary Pauli operator as $\\sigma \\in \\{ I, \\sigma_x, \\sigma_y\n",
        "\\sigma_z\\}$. The action of the Pauli operator on a computational basis state\n",
        "can then be represented by the logic operation:\n",
        "\n",
        "$$\n",
        "\\sigma |x \\rangle = |x == \\textrm{diag}(\\sigma) \\rangle (-1)^{x\\textrm{ and sign}(\\sigma)}\n",
        "(i)^{\\textrm{imag}(\\sigma)}.\n",
        "$$\n",
        "\n",
        "The same is straightforwardly generalized to arbitrary number of qubits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a8fd7e11",
      "metadata": {},
      "source": [
        "Let's check that this works:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "dcb15308",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "-------------------\n",
            "I\n",
            "|False> -->  |False>    ME:(1+0j)\n",
            "|True> -->  |True>    ME:(1+0j)\n",
            "-------------------\n",
            "SX\n",
            "|False> -->  |True>    ME:(1+0j)\n",
            "|True> -->  |False>    ME:(1+0j)\n",
            "-------------------\n",
            "SZ\n",
            "|False> -->  |False>    ME:(1+0j)\n",
            "|True> -->  |True>    ME:(-1+0j)\n",
            "-------------------\n",
            "SY\n",
            "|False> -->  |True>    ME:1j\n",
            "|True> -->  |False>    ME:(-0-1j)\n"
          ]
        }
      ],
      "source": [
        "def connected_element_and_amplitude_bool(\n",
        "    x: bool, diag: bool, sign: bool, imag: bool\n",
        ") -> tuple[bool, complex]:\n",
        "    \"\"\"\n",
        "    Finds the connected element to computational basis state |x> under\n",
        "    the action of the Pauli operator represented by (diag, sign, imag).\n",
        "\n",
        "    Args:\n",
        "        x: Value of the bit, either True or False.\n",
        "        diag: Whether the Pauli operator is diagonal (I, Z)\n",
        "        sigma: Whether the Pauli operator's rows differ in sign (Y, Z)\n",
        "        imag: Whether the Pauli operator is purely imaginary (Y)\n",
        "\n",
        "    Returns:\n",
        "        A length-2 tuple:\n",
        "            - The connected element to x, either False or True\n",
        "            - The matrix element\n",
        "    \"\"\"\n",
        "    return x == diag, (-1) ** (x and sign) * (1j) ** (imag)\n",
        "\n",
        "\n",
        "sigma_indices = [0, 1, 2, 3]\n",
        "sigma_string = [\"I\", \"SX\", \"SZ\", \"SY\"]\n",
        "sigma_diag = [True, False, True, False]\n",
        "sigma_sign = [False, False, True, True]\n",
        "sigma_imag = [False, False, False, True]\n",
        "qubit_values = [False, True]\n",
        "\n",
        "for xi in sigma_indices:\n",
        "    print(\"-------------------\")\n",
        "    print(sigma_string[xi])\n",
        "    for x in qubit_values:\n",
        "        x_p, matrix_element = connected_element_and_amplitude_bool(\n",
        "            x, sigma_diag[xi], sigma_sign[xi], sigma_imag[xi]\n",
        "        )\n",
        "        print(\n",
        "            \"|\"\n",
        "            + str(x)\n",
        "            + \"> -->  |\"\n",
        "            + str(x_p)\n",
        "            + \">    ME:\"\n",
        "            + str(matrix_element)\n",
        "        )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f40a3463",
      "metadata": {},
      "source": [
        "We generate some large number of bitstrings (50 M) for a 40-qubit system:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "74c16c91-cc5c-46ce-aff8-17d0e71ac50f",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Total number of unique bitstrings: 49998839\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from qiskit_addon_sqd.qubit import sort_and_remove_duplicates\n",
        "\n",
        "rand_seed = 22\n",
        "np.random.seed(rand_seed)\n",
        "\n",
        "# Generate some random bitstrings for testing\n",
        "\n",
        "\n",
        "def random_bitstrings(n_samples, n_qubits):\n",
        "    return (\n",
        "        np.round(np.random.rand(n_samples, n_qubits))\n",
        "        .astype(\"int\")\n",
        "        .astype(\"bool\")\n",
        "    )\n",
        "\n",
        "\n",
        "n_qubits = 40\n",
        "bts_matrix = random_bitstrings(50_000_000, n_qubits)\n",
        "\n",
        "# We need to sort the bitstrings and only keep the unique ones\n",
        "# NOTE: It is essential for the projection code to have the bitstrings sorted!\n",
        "bts_matrix = sort_and_remove_duplicates(bts_matrix).astype(\"bool\")\n",
        "\n",
        "# Final subspace dimension after getting rid of duplicated bitstrings\n",
        "d = bts_matrix.shape[0]\n",
        "\n",
        "print(\"Total number of unique bitstrings: \" + str(d))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "184bf287",
      "metadata": {},
      "source": [
        "### Benchmark SQD Pauli projection functions\n",
        "\n",
        "The Pauli string under consideration is $\\sigma_z \\otimes ... \\otimes \\sigma_z$.\n",
        "\n",
        "Different subspace dimensions are considered by slicing the matrix of bitstrings. We time the subspace projection for the different subspace sizes.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "8fe182bc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 0 took 0.201246s\n",
            "Iteration 1 took 0.348222s\n",
            "Iteration 2 took 0.576333s\n",
            "Iteration 3 took 0.78356s\n",
            "Iteration 4 took 1.016162s\n",
            "Iteration 5 took 1.305325s\n",
            "Iteration 6 took 1.392751s\n",
            "Iteration 7 took 1.632433s\n",
            "Iteration 8 took 1.826521s\n",
            "Iteration 9 took 2.02903s\n",
            "Iteration 10 took 2.297458s\n",
            "Iteration 11 took 2.588042s\n",
            "Iteration 12 took 2.738746s\n",
            "Iteration 13 took 2.906144s\n",
            "Iteration 14 took 3.148833s\n",
            "Iteration 15 took 3.323253s\n",
            "Iteration 16 took 3.664171s\n",
            "Iteration 17 took 3.680663s\n",
            "Iteration 18 took 4.008313s\n",
            "Iteration 19 took 4.173532s\n"
          ]
        }
      ],
      "source": [
        "import time\n",
        "\n",
        "from qiskit.quantum_info import Pauli\n",
        "from qiskit_addon_sqd.qubit import matrix_elements_from_pauli\n",
        "\n",
        "pauli = Pauli(\"Z\" * n_qubits)\n",
        "\n",
        "# Different subspace sizes to test\n",
        "d_list = np.linspace(d / 1000, d, 20).astype(\"int\")\n",
        "\n",
        "# To store the walltime\n",
        "time_array = np.zeros(20)\n",
        "\n",
        "for i in range(20):\n",
        "    int_bts_matrix = bts_matrix[: d_list[i], :]\n",
        "    time_1 = time.time()\n",
        "    _ = matrix_elements_from_pauli(int_bts_matrix, pauli)\n",
        "    time_array[i] = time.time() - time_1\n",
        "    print(f\"Iteration {i} took {round(time_array[i], 6)}s\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "9abb110f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/qiskit-addon-sqd/guides/benchmark-pauli-projection/extracted-outputs/9abb110f-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Data for energies plot\n",
        "x1 = d_list\n",
        "y1 = time_array\n",
        "\n",
        "# Plot energies\n",
        "plt.title(\"Runtime vs subspace dimension 40 qubits\")\n",
        "plt.xlabel(\"Subspace dimension (millions)\")\n",
        "plt.ylabel(\"Wall time [s]\")\n",
        "plt.xticks([1e7, 2e7, 3e7, 4e7, 5e7], [str(i) for i in [10, 20, 30, 40, 50]])\n",
        "plt.plot(x1, y1, marker=\".\", markersize=20)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8c486bc7",
      "metadata": {},
      "source": [
        "We now do the same for 60 qubits:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "359ed3f3",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Total number of unique bitstrings: 50000000\n"
          ]
        }
      ],
      "source": [
        "n_qubits = 60\n",
        "bts_matrix = random_bitstrings(50_000_000, n_qubits)\n",
        "\n",
        "# We need to sort the bitstrings and just keep the unique ones\n",
        "bts_matrix = sort_and_remove_duplicates(bts_matrix).astype(\"bool\")\n",
        "\n",
        "# Final subspace dimension after getting rid of duplicated bitstrings\n",
        "d = bts_matrix.shape[0]\n",
        "\n",
        "print(\"Total number of unique bitstrings: \" + str(d))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "7bb0d8d8",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 0 took 0.236567s\n",
            "Iteration 1 took 0.424116s\n",
            "Iteration 2 took 0.673399s\n",
            "Iteration 3 took 0.905164s\n",
            "Iteration 4 took 1.168936s\n",
            "Iteration 5 took 1.454204s\n",
            "Iteration 6 took 1.74778s\n",
            "Iteration 7 took 1.920795s\n",
            "Iteration 8 took 2.259994s\n",
            "Iteration 9 took 2.550674s\n",
            "Iteration 10 took 2.681287s\n",
            "Iteration 11 took 3.04411s\n",
            "Iteration 12 took 3.293262s\n",
            "Iteration 13 took 3.471247s\n",
            "Iteration 14 took 3.726639s\n",
            "Iteration 15 took 4.072854s\n",
            "Iteration 16 took 4.221037s\n",
            "Iteration 17 took 4.498535s\n",
            "Iteration 18 took 4.741108s\n",
            "Iteration 19 took 5.159038s\n"
          ]
        }
      ],
      "source": [
        "pauli = Pauli(\"Z\" * n_qubits)\n",
        "\n",
        "# Different subspace sizes to test\n",
        "d_list = np.linspace(d / 1000, d, 20).astype(\"int\")\n",
        "\n",
        "# It is better to do this once\n",
        "row_array = np.arange(d)\n",
        "\n",
        "# To store the walltime\n",
        "time_array = np.zeros(20)\n",
        "\n",
        "for i in range(20):\n",
        "    int_bts_matrix = bts_matrix[: d_list[i], :]\n",
        "    int_row_array = row_array[: d_list[i]]\n",
        "    time_1 = time.time()\n",
        "    _ = matrix_elements_from_pauli(int_bts_matrix, pauli)\n",
        "    time_array[i] = time.time() - time_1\n",
        "    print(f\"Iteration {i} took {round(time_array[i], 6)}s\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "6b961c81",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/qiskit-addon-sqd/guides/benchmark-pauli-projection/extracted-outputs/6b961c81-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Data for energies plot\n",
        "x1 = d_list\n",
        "y1 = time_array\n",
        "\n",
        "fig, axs = plt.subplots(1, 1, figsize=(6, 6))\n",
        "\n",
        "# Plot energies\n",
        "axs.plot(x1, y1, marker=\".\", markersize=20)\n",
        "axs.set_title(\"Runtime vs subspace dimension 60 qubits\")\n",
        "axs.set_xlabel(\"Subspace dimension (millions)\")\n",
        "plt.xticks([1e7, 2e7, 3e7, 4e7, 5e7], [str(i) for i in [10, 20, 30, 40, 50]])\n",
        "axs.set_ylabel(\"Wall time [s]\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "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"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}