{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "509f7bd9-b597-4d23-b3af-0a76a7b4d33d",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Diagonalização quântica de Krylov de hamiltonianos de rede\"\n",
        "description: \"Implemente o Algoritmo de Diagonalização Quântica de Krylov (KQD) no contexto dos padrões Qiskit.\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore prefactors */}\n",
        "\n",
        "<span id=\"krylov-quantum-diagonalization-of-lattice-hamiltonians\" />\n",
        "\n",
        "# Diagonalização quântica de Krylov de hamiltonianos de rede\n",
        "\n",
        "*Estimativa de uso: 20 minutos em um Heron r2 (OBSERVAÇÃO: essa é apenas uma estimativa. Seu tempo de execução pode variar)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "921c7b04-5b5d-4cfc-aba0-15a61334e619",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## Segundo plano\n",
        "\n",
        "Este tutorial demonstra como implementar o Algoritmo de Diagonalização Quântica de Krylov (KQD) no contexto dos padrões Qiskit. Primeiro, você aprenderá sobre a teoria por trás do algoritmo e, em seguida, verá uma demonstração de sua execução em uma QPU.\n",
        "\n",
        "Em todas as disciplinas, estamos interessados em aprender as propriedades do estado fundamental dos sistemas quânticos. Os exemplos incluem a compreensão da natureza fundamental das partículas e forças, a previsão e a compreensão do comportamento de materiais complexos e a compreensão das interações e reações bioquímicas. Devido ao crescimento exponencial do espaço de Hilbert e à correlação que surge em sistemas emaranhados, o algoritmo clássico tem dificuldade para resolver esse problema em sistemas quânticos de tamanho cada vez maior. Em uma extremidade do espectro está a abordagem existente que aproveita o foco do hardware quântico em métodos quânticos variacionais (por exemplo, [eigensolver quântico variacional](/docs/tutorials/spin-chain-vqe) ). Essas técnicas enfrentam desafios com os dispositivos atuais devido ao alto número de chamadas de função necessárias no processo de otimização, o que adiciona uma grande sobrecarga de recursos quando são introduzidas técnicas avançadas de atenuação de erros, limitando assim sua eficácia a sistemas pequenos. No outro extremo do espectro, há métodos quânticos tolerantes a falhas com garantias de desempenho (por exemplo, [estimativa de fase quântica](https://arxiv.org/abs/quant-ph/0604193) ), que exigem circuitos profundos que só podem ser executados em um dispositivo tolerante a falhas. Por esses motivos, apresentamos aqui um algoritmo quântico baseado em métodos de subespaço (conforme descrito neste [artigo de revisão](https://arxiv.org/abs/2312.00178) ), o algoritmo de diagonalização quântica de Krylov (KQD). Esse algoritmo tem bom desempenho em larga escala [\\[1\\]](#references) no hardware quântico existente, compartilha [garantias de desempenho](https://arxiv.org/abs/2110.07492) semelhantes às da estimativa de fase, é compatível com técnicas avançadas de atenuação de erros e pode fornecer resultados que são classicamente inacessíveis.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5f698c82-95ca-4dc4-b9fa-d6e741e2c02c",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Requisitos\n",
        "\n",
        "Antes de iniciar este tutorial, verifique se você tem os seguintes itens instalados:\n",
        "\n",
        "* Qiskit SDK v2.0 ou posterior, com suporte [para visualização](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 ou mais tarde ( `pip install qiskit-ibm-runtime` )\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c44956e4-47ab-4b0f-9d6d-553080110062",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Instalação\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ce7d5adc-3ef2-4654-b865-14d5141ce41a",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import scipy as sp\n",
        "import matplotlib.pylab as plt\n",
        "from typing import Union, List\n",
        "import itertools as it\n",
        "import copy\n",
        "from sympy import Matrix\n",
        "import warnings\n",
        "\n",
        "warnings.filterwarnings(\"ignore\")\n",
        "\n",
        "from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState\n",
        "from qiskit.circuit import Parameter, IfElseOp\n",
        "from qiskit import QuantumCircuit, QuantumRegister\n",
        "from qiskit.circuit.library import PauliEvolutionGate\n",
        "from qiskit.synthesis import LieTrotter\n",
        "from qiskit.transpiler import Target, CouplingMap\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "\n",
        "\n",
        "from qiskit_ibm_runtime import (\n",
        "    QiskitRuntimeService,\n",
        "    EstimatorV2 as Estimator,\n",
        ")\n",
        "\n",
        "\n",
        "def solve_regularized_gen_eig(\n",
        "    h: np.ndarray,\n",
        "    s: np.ndarray,\n",
        "    threshold: float,\n",
        "    k: int = 1,\n",
        "    return_dimn: bool = False,\n",
        ") -> Union[float, List[float]]:\n",
        "    \"\"\"\n",
        "    Method for solving the generalized eigenvalue problem with regularization\n",
        "\n",
        "    Args:\n",
        "        h (numpy.ndarray):\n",
        "            The effective representation of the matrix in the Krylov subspace\n",
        "        s (numpy.ndarray):\n",
        "            The matrix of overlaps between vectors of the Krylov subspace\n",
        "        threshold (float):\n",
        "            Cut-off value for the eigenvalue of s\n",
        "        k (int):\n",
        "            Number of eigenvalues to return\n",
        "        return_dimn (bool):\n",
        "            Whether to return the size of the regularized subspace\n",
        "\n",
        "    Returns:\n",
        "        lowest k-eigenvalue(s) that are the solution of the\n",
        "        regularized generalized eigenvalue problem\n",
        "\n",
        "\n",
        "    \"\"\"\n",
        "    s_vals, s_vecs = sp.linalg.eigh(s)\n",
        "    s_vecs = s_vecs.T\n",
        "    good_vecs = np.array(\n",
        "        [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]\n",
        "    )\n",
        "    h_reg = good_vecs.conj() @ h @ good_vecs.T\n",
        "    s_reg = good_vecs.conj() @ s @ good_vecs.T\n",
        "    if k == 1:\n",
        "        if return_dimn:\n",
        "            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)\n",
        "        else:\n",
        "            return sp.linalg.eigh(h_reg, s_reg)[0][0]\n",
        "    else:\n",
        "        if return_dimn:\n",
        "            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)\n",
        "        else:\n",
        "            return sp.linalg.eigh(h_reg, s_reg)[0][:k]\n",
        "\n",
        "\n",
        "def single_particle_gs(H_op, n_qubits):\n",
        "    \"\"\"\n",
        "    Find the ground state of the single particle(excitation) sector\n",
        "    \"\"\"\n",
        "    H_x = []\n",
        "    for p, coeff in H_op.to_list():\n",
        "        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))\n",
        "\n",
        "    H_z = []\n",
        "    for p, coeff in H_op.to_list():\n",
        "        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))\n",
        "\n",
        "    H_c = H_op.coeffs\n",
        "\n",
        "    print(\"n_sys_qubits\", n_qubits)\n",
        "\n",
        "    n_exc = 1\n",
        "    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))\n",
        "    print(\"n_exc\", n_exc, \", subspace dimension\", sub_dimn)\n",
        "\n",
        "    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)\n",
        "\n",
        "    # list all of the possible sets of n_exc indices of 1s in\n",
        "    # n_exc-particle states\n",
        "    sparse_vecs = [\n",
        "        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)\n",
        "    ]\n",
        "\n",
        "    m = 0\n",
        "    for i, i_set in enumerate(sparse_vecs):\n",
        "        for j, j_set in enumerate(sparse_vecs):\n",
        "            m += 1\n",
        "\n",
        "            if len(i_set.symmetric_difference(j_set)) <= 2:\n",
        "                for p_x, p_z, coeff in zip(H_x, H_z, H_c):\n",
        "                    if i_set.symmetric_difference(j_set) == p_x:\n",
        "                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (\n",
        "                            (-1) ** len(i_set.intersection(p_z))\n",
        "                        )\n",
        "                    else:\n",
        "                        sgn = 0\n",
        "\n",
        "                    few_particle_H[i, j] += sgn * coeff\n",
        "\n",
        "    gs_en = min(np.linalg.eigvalsh(few_particle_H))\n",
        "    print(\"single particle ground state energy: \", gs_en)\n",
        "    return gs_en"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "76a1c616-a749-48a2-b5fd-8beff760760b",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "## Passo 1: Mapear entradas clássicas para um problema quântico\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a166e2a1-3003-4799-9c15-59f5e78b6ea8",
      "metadata": {},
      "source": [
        "<span id=\"the-krylov-space\" />\n",
        "\n",
        "### O espaço de Krylov\n",
        "\n",
        "O espaço de Krylov $\\mathcal{K}^r$ de ordem $r$ é o espaço abrangido pelos vetores obtidos pela multiplicação das potências mais altas de uma matriz $A$, até $r-1$, com um vetor de referência $\\vert v \\rangle$.\n",
        "\n",
        "$$\n",
        "\\mathcal{K}^r = \\left\\{ \\vert v \\rangle, A \\vert v \\rangle, A^2 \\vert v \\rangle, ..., A^{r-1} \\vert v \\rangle \\right\\}\n",
        "$$\n",
        "\n",
        "Se a matriz $A$ for o Hamiltoniano $H$, nos referiremos ao espaço correspondente como o espaço de Krylov de potência $\\mathcal{K}_P$. No caso em que $A$ é o operador de evolução temporal gerado pelo Hamiltoniano $U=e^{-iHt}$, vamos nos referir ao espaço como o espaço de Krylov unitário $\\mathcal{K}_U$. O subespaço de Krylov de potência que usamos classicamente não pode ser gerado diretamente em um computador quântico, pois $H$ não é um operador unitário. Em vez disso, podemos usar o operador de evolução temporal $U = e^{-iHt}$, que pode ser mostrado como [garantia de convergência](https://arxiv.org/abs/2110.07492) semelhante ao método de potência. As potências de $U$ tornam-se, então, etapas de tempo diferentes $U^k = e^{-iH(kt)}$.\n",
        "\n",
        "$$\n",
        "\\mathcal{K}_U^r = \\left\\{ \\vert \\psi \\rangle, U \\vert \\psi \\rangle, U^2 \\vert \\psi \\rangle, ..., U^{r-1} \\vert \\psi \\rangle \\right\\}\n",
        "$$\n",
        "\n",
        "Consulte o Apêndice para obter uma derivação detalhada de como o espaço unitário de Krylov permite representar com precisão os estados próprios de baixa energia.\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "5573ca7e-ab16-4488-88b3-d8d1eba9e20c",
      "metadata": {},
      "source": [
        "<span id=\"krylov-quantum-diagonalization-algorithm\" />\n",
        "\n",
        "### Algoritmo de diagonalização quântica de Krylov\n",
        "\n",
        "Dado um Hamiltoniano $H$ que desejamos diagonalizar, primeiro consideramos o espaço de Krylov unitário correspondente $\\mathcal{K}_U$. O objetivo é encontrar uma representação compacta do Hamiltoniano em $\\mathcal{K}_U$, que chamaremos de $\\tilde{H}$. Os elementos da matriz de $\\tilde{H}$, a projeção do Hamiltoniano no espaço de Krylov, podem ser calculados pelos seguintes valores de expectativa\n",
        "\n",
        "$$\n",
        "\\tilde{H}_{mn} = \\langle \\psi_m \\vert H \\vert \\psi_n \\rangle =\n",
        "$$\n",
        "\n",
        "$$\n",
        "= \\langle \\psi \\vert  e^{i H t_m}   H e^{-i H t_n} \\vert \\psi \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "= \\langle \\psi \\vert  e^{i H m dt}   H e^{-i H n dt} \\vert \\psi \\rangle\n",
        "$$\n",
        "\n",
        "Onde $\\vert \\psi_n \\rangle = e^{-i H t_n} \\vert \\psi \\rangle$ são os vetores do espaço unitário de Krylov e $t_n = n dt$ são os múltiplos da etapa de tempo $dt$ escolhida. Em um computador quântico, o cálculo de cada elemento da matriz pode ser feito com qualquer algoritmo que permita obter a sobreposição entre os estados quânticos. Este tutorial se concentra no teste de Hadamard. Como o $\\mathcal{K}_U$ tem dimensão $r$, o Hamiltoniano projetado no subespaço terá dimensões $r \\times r$. Com $r$ suficientemente pequeno (geralmente $r<<100$ é suficiente para obter a convergência das estimativas das energias próprias), podemos então diagonalizar facilmente o Hamiltoniano projetado $\\tilde{H}$. Entretanto, não podemos diagonalizar diretamente $\\tilde{H}$ devido à não ortogonalidade dos vetores do espaço de Krylov. Teremos que medir suas sobreposições e construir uma matriz $\\tilde{S}$\n",
        "\n",
        "$$\n",
        "\\tilde{S}_{mn} = \\langle \\psi_m \\vert \\psi_n \\rangle\n",
        "$$\n",
        "\n",
        "Isso nos permite resolver o problema do valor próprio em um espaço não ortogonal (também chamado de problema do valor próprio generalizado)\n",
        "\n",
        "$$\n",
        "\\tilde{H} \\ \\vec{c} = E \\ \\tilde{S} \\ \\vec{c}\n",
        "$$\n",
        "\n",
        "Pode-se, então, obter estimativas dos autovalores e dos estados próprios de $H$ observando os de $\\tilde{H}$. Por exemplo, a estimativa da energia do estado fundamental é obtida tomando-se o menor autovalor de $c$ e o estado fundamental do vetor próprio correspondente $\\vec{c}$. Os coeficientes em $\\vec{c}$ determinam a contribuição dos diferentes vetores que abrangem $\\mathcal{K}_U$.\n",
        "\n",
        "![fig1.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/fc662b76-8ad7-4a6c-8c49-5f08c125aee8.avif)\n",
        "\n",
        "A figura mostra uma representação do circuito do teste de Hadamard modificado, um método usado para calcular a sobreposição entre diferentes estados quânticos. Para cada elemento da matriz $\\tilde{H}_{i,j}$, é realizado um teste de Hadamard entre os estados $\\vert \\psi_i \\rangle$, $\\vert \\psi_j \\rangle$. Isso é destacado na figura pelo esquema de cores dos elementos da matriz e das operações $\\text{Prep} \\; \\psi_i$, $\\text{Prep} \\; \\psi_j$ correspondentes. Assim, é necessário um conjunto de testes de Hadamard para todas as combinações possíveis de vetores do espaço de Krylov para calcular todos os elementos da matriz do Hamiltoniano projetado $\\tilde{H}$. O fio superior no circuito de teste de Hadamard é um qubit ancilla que é medido na base X ou Y, e seu valor de expectativa determina o valor da sobreposição entre os estados. O fio inferior representa todos os qubits do sistema Hamiltoniano. A operação $\\text{Prep} \\; \\psi_i$ prepara o qubit do sistema no estado $\\vert \\psi_i \\rangle$ controlado pelo estado do qubit ancilla (da mesma forma para $\\text{Prep} \\; \\psi_j$ ) e a operação $P$ representa a decomposição de Pauli do sistema Hamiltoniano $H = \\sum_i P_i$. Uma derivação mais detalhada das operações calculadas pelo teste de Hadamard é apresentada abaixo.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1a6d7a4f-6c1f-4069-93d1-f5b670645d7f",
      "metadata": {},
      "source": [
        "<span id=\"define-hamiltonian\" />\n",
        "\n",
        "#### Defina Hamiltoniano\n",
        "\n",
        "Vamos considerar o Hamiltoniano de Heisenberg para $N$ qubits em uma cadeia linear: $H= \\sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "82163249-bafd-4bc6-9741-20a4077971a7",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]\n"
          ]
        }
      ],
      "source": [
        "# Define problem Hamiltonian.\n",
        "n_qubits = 30\n",
        "J = 1  # coupling strength for ZZ interaction\n",
        "\n",
        "# Define the Hamiltonian:\n",
        "H_int = [[\"I\"] * n_qubits for _ in range(3 * (n_qubits - 1))]\n",
        "for i in range(n_qubits - 1):\n",
        "    H_int[i][i] = \"Z\"\n",
        "    H_int[i][i + 1] = \"Z\"\n",
        "for i in range(n_qubits - 1):\n",
        "    H_int[n_qubits - 1 + i][i] = \"X\"\n",
        "    H_int[n_qubits - 1 + i][i + 1] = \"X\"\n",
        "for i in range(n_qubits - 1):\n",
        "    H_int[2 * (n_qubits - 1) + i][i] = \"Y\"\n",
        "    H_int[2 * (n_qubits - 1) + i][i + 1] = \"Y\"\n",
        "H_int = [\"\".join(term) for term in H_int]\n",
        "H_tot = [(term, J) if term.count(\"Z\") == 2 else (term, 1) for term in H_int]\n",
        "\n",
        "# Get operator\n",
        "H_op = SparsePauliOp.from_list(H_tot)\n",
        "print(H_tot)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f26578ba-24ad-42f3-baad-c348d3c05699",
      "metadata": {},
      "source": [
        "<span id=\"set-parameters-for-the-algorithm\" />\n",
        "\n",
        "#### Definir parâmetros para o algoritmo\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8d376c0e-fb7a-4a41-838a-ae041d5c9afa",
      "metadata": {},
      "source": [
        "Escolhemos heuristicamente um valor para a etapa de tempo `dt` (com base nos limites superiores da norma Hamiltoniana). A Ref [\\[2\\]](#references) mostrou que um intervalo de tempo suficientemente pequeno é $\\pi/\\vert \\vert H \\vert \\vert$, e que é preferível, até certo ponto, subestimar esse valor em vez de superestimá-lo, pois a superestimação pode permitir que as contribuições de estados de alta energia corrompam até mesmo o estado ideal no espaço de Krylov. Por outro lado, escolher $dt$ para ser muito pequeno leva a um condicionamento pior do subespaço de Krylov, já que os vetores da base de Krylov diferem menos de um intervalo de tempo para outro.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "963dd2e9-f4ea-456d-b4ed-e3d05a114082",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "np.float64(0.10833078115826875)"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Get Hamiltonian restricted to single-particle states\n",
        "single_particle_H = np.zeros((n_qubits, n_qubits))\n",
        "for i in range(n_qubits):\n",
        "    for j in range(i + 1):\n",
        "        for p, coeff in H_op.to_list():\n",
        "            p_x = Pauli(p).x\n",
        "            p_z = Pauli(p).z\n",
        "            if all(\n",
        "                p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)\n",
        "            ):\n",
        "                sgn = (\n",
        "                    (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))\n",
        "                ) * ((-1) ** p_z[i])\n",
        "            else:\n",
        "                sgn = 0\n",
        "            single_particle_H[i, j] += sgn * coeff\n",
        "for i in range(n_qubits):\n",
        "    for j in range(i + 1, n_qubits):\n",
        "        single_particle_H[i, j] = np.conj(single_particle_H[j, i])\n",
        "\n",
        "# Set dt according to spectral norm\n",
        "dt = np.pi / np.linalg.norm(single_particle_H, ord=2)\n",
        "dt"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d2bcc2b9-ca7c-4147-8e6b-0f15208297ce",
      "metadata": {},
      "source": [
        "E definir outros parâmetros do algoritmo. Para fins deste tutorial, vamos nos limitar a usar um espaço de Krylov com apenas cinco dimensões, o que é bastante limitado.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "ddde76b8-446f-4cd5-bc4b-8f00e2ec726c",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Set parameters for quantum Krylov algorithm\n",
        "krylov_dim = 5  # size of Krylov subspace\n",
        "num_trotter_steps = 6\n",
        "dt_circ = dt / num_trotter_steps"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3c0dcc08-6963-4e14-a50c-21c5dd7f8ade",
      "metadata": {},
      "source": [
        "<span id=\"state-preparation\" />\n",
        "\n",
        "#### Preparação do estado\n",
        "\n",
        "Escolha um estado de referência $\\vert \\psi \\rangle$ que tenha alguma sobreposição com o estado fundamental. Para esse Hamiltoniano, usamos o estado a com uma excitação no qubit do meio $\\vert 00..010...00 \\rangle$ como nosso estado de referência.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "70161afe-8ace-4642-894a-cd21ed77a3b9",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/70161afe-8ace-4642-894a-cd21ed77a3b9-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "qc_state_prep = QuantumCircuit(n_qubits)\n",
        "qc_state_prep.x(int(n_qubits / 2) + 1)\n",
        "qc_state_prep.draw(\"mpl\", scale=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1786de6b-476f-47ea-9c20-3464d3cfbbdb",
      "metadata": {},
      "source": [
        "<span id=\"time-evolution\" />\n",
        "\n",
        "#### Evolução temporal\n",
        "\n",
        "Podemos realizar o operador de evolução temporal gerado por um determinado Hamiltoniano: $U=e^{-iHt}$ por meio da [aproximação de Lie-Trotter](/docs/api/qiskit/qiskit.synthesis.LieTrotter).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "a23b9e9c-5dc8-447f-8c73-b4fa01630c8f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<qiskit.circuit.instructionset.InstructionSet at 0x11eef9be0>"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "t = Parameter(\"t\")\n",
        "\n",
        "## Create the time-evo op circuit\n",
        "evol_gate = PauliEvolutionGate(\n",
        "    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)\n",
        ")\n",
        "\n",
        "qr = QuantumRegister(n_qubits)\n",
        "qc_evol = QuantumCircuit(qr)\n",
        "qc_evol.append(evol_gate, qargs=qr)"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "41c2cc43-51a2-42c3-88ee-a3b7b30eabb4",
      "metadata": {},
      "source": [
        "<span id=\"hadamard-test\" />\n",
        "\n",
        "#### teste de Hadamard\n",
        "\n",
        "![fig2.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/c5263851-6067-4ca2-8e0c-a835631cdc7f.avif)\n",
        "\n",
        "$$\n",
        "\\begin{equation*}\n",
        "    |0\\rangle|0\\rangle^N \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle + |1\\rangle \\Big)|0\\rangle^N \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle|0\\rangle^N+|1\\rangle |\\psi_i\\rangle\\Big) \\quad\\longrightarrow\\quad \\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |0\\rangle^N+|1\\rangle P |\\psi_i\\rangle\\Big) \\quad\\longrightarrow\\quad\\frac{1}{\\sqrt{2}}\\Big(|0\\rangle |\\psi_j\\rangle+|1\\rangle P|\\psi_i\\rangle\\Big)\n",
        "\\end{equation*}\n",
        "$$\n",
        "\n",
        "Onde $P$ é um dos termos na decomposição do Hamiltoniano $H=\\sum P$ e $\\text{Prep} \\; \\psi_i$, $\\text{Prep} \\; \\psi_j$ são operações controladas que preparam $|\\psi_i\\rangle$, $|\\psi_j\\rangle$ vetores do espaço unitário de Krylov, com $|\\psi_k\\rangle = e^{-i H k dt } \\vert \\psi \\rangle = e^{-i H k dt } U_{\\psi} \\vert 0 \\rangle^N$. Para medir $X$, primeiro aplique $H$...\n",
        "\n",
        "$$\n",
        "\\begin{equation*}\n",
        "    \\longrightarrow\\quad\\frac{1}{2}|0\\rangle\\Big( |\\psi_j\\rangle + P|\\psi_i\\rangle\\Big) + \\frac{1}{2}|1\\rangle\\Big(|\\psi_j\\rangle - P|\\psi_i\\rangle\\Big)\n",
        "\\end{equation*}\n",
        "$$\n",
        "\n",
        "... e depois medir:\n",
        "\n",
        "$$\n",
        "\\begin{equation*}\n",
        "\\begin{split}\n",
        "    \\Rightarrow\\quad\\langle X\\rangle &= \\frac{1}{4}\\Bigg(\\Big\\|| \\psi_j\\rangle + P|\\psi_i\\rangle \\Big\\|^2-\\Big\\||\\psi_j\\rangle - P|\\psi_i\\rangle\\Big\\|^2\\Bigg) \\\\\n",
        "    &= \\text{Re}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n",
        "\\end{split}\n",
        "\\end{equation*}\n",
        "$$\n",
        "\n",
        "A partir da identidade $|a + b\\|^2 = \\langle a + b | a + b \\rangle = \\|a\\|^2 + \\|b\\|^2 + 2\\text{Re}\\langle a | b \\rangle$. Da mesma forma, a medição de $Y$ resulta em\n",
        "\n",
        "$$\n",
        "\\begin{equation*}\n",
        "    \\langle Y\\rangle = \\text{Im}\\Big[\\langle\\psi_j| P|\\psi_i\\rangle\\Big].\n",
        "\\end{equation*}\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "7c1efca7-7db9-43a9-bcb7-053ba274d6f6",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Circuit for calculating the real part of the overlap in S via Hadamard test\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/7c1efca7-7db9-43a9-bcb7-053ba274d6f6-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "## Create the time-evo op circuit\n",
        "evol_gate = PauliEvolutionGate(\n",
        "    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n",
        ")\n",
        "\n",
        "## Create the time-evo op dagger circuit\n",
        "evol_gate_d = PauliEvolutionGate(\n",
        "    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)\n",
        ")\n",
        "evol_gate_d = evol_gate_d.inverse()\n",
        "\n",
        "# Put pieces together\n",
        "qc_reg = QuantumRegister(n_qubits)\n",
        "qc_temp = QuantumCircuit(qc_reg)\n",
        "qc_temp.compose(qc_state_prep, inplace=True)\n",
        "for _ in range(num_trotter_steps):\n",
        "    qc_temp.append(evol_gate, qargs=qc_reg)\n",
        "for _ in range(num_trotter_steps):\n",
        "    qc_temp.append(evol_gate_d, qargs=qc_reg)\n",
        "qc_temp.compose(qc_state_prep.inverse(), inplace=True)\n",
        "\n",
        "# Create controlled version of the circuit\n",
        "controlled_U = qc_temp.to_gate().control(1)\n",
        "\n",
        "# Create hadamard test circuit for real part\n",
        "qr = QuantumRegister(n_qubits + 1)\n",
        "qc_real = QuantumCircuit(qr)\n",
        "qc_real.h(0)\n",
        "qc_real.append(controlled_U, list(range(n_qubits + 1)))\n",
        "qc_real.h(0)\n",
        "\n",
        "print(\n",
        "    \"Circuit for calculating the real part of the overlap in S via Hadamard test\"\n",
        ")\n",
        "qc_real.draw(\"mpl\", fold=-1, scale=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c72fd5e5-a586-498e-bbc1-0f48ecea49d5",
      "metadata": {},
      "source": [
        "O circuito de teste Hadamard pode ser um circuito profundo quando decomposto em portas nativas (o que aumentará ainda mais se levarmos em conta a topologia do dispositivo)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "06f22c05-6e1c-4540-9dac-884b3400a50e",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Number of layers of 2Q operations 112753\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    \"Number of layers of 2Q operations\",\n",
        "    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),\n",
        ")"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "8e15903d-d4a8-4fbe-8f9e-e960e686629e",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "## Etapa 2: Otimizar o problema para execução em hardware quântico\n",
        "\n",
        "<span id=\"efficient-hadamard-test\" />\n",
        "\n",
        "### Teste de Hadamard eficiente\n",
        "\n",
        "Podemos otimizar os circuitos profundos para o teste de Hadamard que obtivemos introduzindo algumas aproximações e confiando em algumas suposições sobre o modelo Hamiltoniano. Por exemplo, considere o seguinte circuito para o teste de Hadamard:\n",
        "\n",
        "![fig3.png](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-subspace-diagonalization/35b13797-5a46-486c-b50e-97c205cc9747.avif)\n",
        "\n",
        "Suponha que possamos calcular classicamente $E_0$, o valor próprio de $|0\\rangle^N$ sob o Hamiltoniano $H$. Isso é satisfeito quando o Hamiltoniano preserva a simetria U(1). Embora isso possa parecer uma suposição forte, há muitos casos em que é seguro supor que há um estado de vácuo (nesse caso, ele é mapeado para o estado $|0\\rangle^N$ ) que não é afetado pela ação do Hamiltoniano. Isso é verdadeiro, por exemplo, para Hamiltonianos químicos que descrevem moléculas estáveis (em que o número de elétrons é conservado).\n",
        "Considerando que a porta $\\text{Prep} \\; \\psi$ prepara o estado de referência desejado $\\ket{psi} = \\text{Prep} \\; \\psi \\ket{0} = e^{-i H 0 dt} U_{\\psi} \\ket{0}$, por exemplo, preparar o estado HF para a química $\\text{Prep} \\; \\psi$ seria um produto de NOTs de um único qubit, portanto, controlado- $\\text{Prep} \\; \\psi$ é apenas um produto de CNOTs.\n",
        "Então, o circuito acima implementa o seguinte estado antes da medição:\n",
        "\n",
        "$$\n",
        "\\begin{equation}\n",
        "\\begin{split}\n",
        "    \\ket{0} \\ket{0}^N\\xrightarrow{H}&\\frac{1}{\\sqrt{2}}\n",
        "    \\left(\n",
        "    \\ket{0}\\ket{0}^N+ \\ket{1} \\ket{0}^N\n",
        "    \\right)\\\\\n",
        "    \\xrightarrow{\\text{1-ctrl-init}}&\\frac{1}{\\sqrt{2}}\\left(|0\\rangle|0\\rangle^N+|1\\rangle|\\psi\\rangle\\right)\\\\\n",
        "    \\xrightarrow{U}&\\frac{1}{\\sqrt{2}}\\left(e^{i\\phi}\\ket{0}\\ket{0}^N+\\ket{1} U\\ket{\\psi}\\right)\\\\\n",
        "    \\xrightarrow{\\text{0-ctrl-init}}&\\frac{1}{\\sqrt{2}}\n",
        "    \\left(\n",
        "    e^{i\\phi}\\ket{0} \\ket{\\psi}\n",
        "    +\\ket{1} U\\ket{\\psi}\n",
        "    \\right)\\\\\n",
        "    =&\\frac{1}{2}\n",
        "    \\left(\n",
        "    \\ket{+}\\left(e^{i\\phi}\\ket{\\psi}+U\\ket{\\psi}\\right)\n",
        "    +\\ket{-}\\left(e^{i\\phi}\\ket{\\psi}-U\\ket{\\psi}\\right)\n",
        "    \\right)\\\\\n",
        "    =&\\frac{1}{2}\n",
        "    \\left(\n",
        "    \\ket{+i}\\left(e^{i\\phi}\\ket{\\psi}-iU\\ket{\\psi}\\right)\n",
        "    +\\ket{-i}\\left(e^{i\\phi}\\ket{\\psi}+iU\\ket{\\psi}\\right)\n",
        "    \\right)\n",
        "\\end{split}\n",
        "\\end{equation}\n",
        "$$\n",
        "\n",
        "em que usamos o deslocamento de fase simulável clássico $ U\\ket{0}^N = e^{i\\phi}\\ket{0}^N$ na terceira linha. Portanto, os valores de expectativa são obtidos como\n",
        "\n",
        "$$\n",
        "\\begin{equation}\n",
        "\\begin{split}\n",
        "    \\langle X\\otimes P\\rangle&=\\frac{1}{4}\n",
        "    \\Big(\n",
        "    \\left(e^{-i\\phi}\\bra{\\psi}+\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}+U\\ket{\\psi}\\right)\n",
        "    \\\\\n",
        "    &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi}-\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}-U\\ket{\\psi}\\right)\n",
        "    \\Big)\\\\\n",
        "    &=\\text{Re}\\left[e^{-i\\phi}\\bra{\\psi}PU\\ket{\\psi}\\right],\n",
        "\\end{split}\n",
        "\\end{equation}\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\begin{equation}\n",
        "\\begin{split}\n",
        "    \\langle Y\\otimes P\\rangle&=\\frac{1}{4}\n",
        "    \\Big(\n",
        "    \\left(e^{-i\\phi}\\bra{\\psi}+i\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}-iU\\ket{\\psi}\\right)\n",
        "    \\\\\n",
        "    &\\qquad-\\left(e^{-i\\phi}\\bra{\\psi}-i\\bra{\\psi}U^\\dagger\\right)P\\left(e^{i\\phi}\\ket{\\psi}+iU\\ket{\\psi}\\right)\n",
        "    \\Big)\\\\\n",
        "    &=\\text{Im}\\left[e^{-i\\phi}\\bra{\\psi}PU\\ket{\\psi}\\right].\n",
        "\\end{split}\n",
        "\\end{equation}\n",
        "$$\n",
        "\n",
        "Usando essas premissas, conseguimos escrever os valores de expectativa dos operadores de interesse com menos operações controladas. Na verdade, só precisamos implementar a preparação do estado controlado $\\text{Prep} \\; \\psi$ e não as evoluções de tempo controladas. Reenquadrar nosso cálculo como acima nos permitirá reduzir bastante a profundidade dos circuitos resultantes.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2e29da5f-b5f5-4b9e-96fc-3fa7ab698398",
      "metadata": {},
      "source": [
        "<span id=\"decompose-time-evolution-operator-with-trotter-decomposition\" />\n",
        "\n",
        "### Decomponha o operador de evolução temporal com a decomposição de Trotter\n",
        "\n",
        "Em vez de implementar exatamente o operador de evolução temporal, podemos usar a decomposição de Trotter para implementar uma aproximação dele. Repetir várias vezes uma determinada ordem de decomposição de Trotter nos proporciona uma redução adicional do erro introduzido pela aproximação. A seguir, criamos diretamente a implementação do Trotter da maneira mais eficiente para o gráfico de interação do Hamiltoniano que estamos considerando (somente interações de vizinhos mais próximos). Na prática, inserimos as rotações de Pauli $R_{xx}$, $R_{yy}$, $R_{zz}$ com um ângulo parametrizado $t$ que corresponde à implementação aproximada de $e^{-i (XX + YY + ZZ) t}$. Dada a diferença na definição das rotações de Pauli e a evolução temporal que estamos tentando implementar, teremos de usar o parâmetro $2*dt$ para obter uma evolução temporal de $dt$. Além disso, invertemos a ordem das operações para um número ímpar de repetições das etapas de Trotter, o que é funcionalmente equivalente, mas permite sintetizar operações adjacentes em uma única $SU(2)$ unitária. Isso proporciona um circuito muito mais raso do que o obtido com o uso da funcionalidade genérica do `PauliEvolutionGate()` .\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "267716dc-fa23-41bd-abe4-6d4e0499a0f4",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/267716dc-fa23-41bd-abe4-6d4e0499a0f4-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "t = Parameter(\"t\")\n",
        "\n",
        "# Create instruction for rotation about XX+YY-ZZ:\n",
        "Rxyz_circ = QuantumCircuit(2)\n",
        "Rxyz_circ.rxx(t, 0, 1)\n",
        "Rxyz_circ.ryy(t, 0, 1)\n",
        "Rxyz_circ.rzz(t, 0, 1)\n",
        "Rxyz_instr = Rxyz_circ.to_instruction(label=\"RXX+YY+ZZ\")\n",
        "\n",
        "interaction_list = [\n",
        "    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],\n",
        "    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],\n",
        "]  # linear chain\n",
        "\n",
        "qr = QuantumRegister(n_qubits)\n",
        "trotter_step_circ = QuantumCircuit(qr)\n",
        "for i, color in enumerate(interaction_list):\n",
        "    for interaction in color:\n",
        "        trotter_step_circ.append(Rxyz_instr, interaction)\n",
        "    if i < len(interaction_list) - 1:\n",
        "        trotter_step_circ.barrier()\n",
        "reverse_trotter_step_circ = trotter_step_circ.reverse_ops()\n",
        "\n",
        "qc_evol = QuantumCircuit(qr)\n",
        "for step in range(num_trotter_steps):\n",
        "    if step % 2 == 0:\n",
        "        qc_evol = qc_evol.compose(trotter_step_circ)\n",
        "    else:\n",
        "        qc_evol = qc_evol.compose(reverse_trotter_step_circ)\n",
        "\n",
        "qc_evol.decompose().draw(\"mpl\", fold=-1, scale=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0a9c3a4d-a678-41e4-b9dd-995fe34341fb",
      "metadata": {},
      "source": [
        "<span id=\"use-an-optimized-circuit-for-state-preparation\" />\n",
        "\n",
        "### Use um circuito otimizado para a preparação do estado\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "70411715-eed3-4cf5-961d-06a6f1e04efc",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/70411715-eed3-4cf5-961d-06a6f1e04efc-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "control = 0\n",
        "excitation = int(n_qubits / 2) + 1\n",
        "controlled_state_prep = QuantumCircuit(n_qubits + 1)\n",
        "controlled_state_prep.cx(control, excitation)\n",
        "controlled_state_prep.draw(\"mpl\", fold=-1, scale=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2421ff11-d97a-488a-a382-a652df8c94d6",
      "metadata": {},
      "source": [
        "<span id=\"template-circuits-for-calculating-matrix-elements-of-$tilde{s}$-and-$tilde{h}$-via-hadamard-test\" />\n",
        "\n",
        "### Circuitos modelo para o cálculo de elementos matriciais de $\\tilde{S}$ e $\\tilde{H}$ através do teste de Hadamard\n",
        "\n",
        "A única diferença entre os circuitos usados no teste de Hadamard será a fase no operador de evolução temporal e os observáveis medidos. Portanto, podemos preparar um circuito modelo que represente o circuito genérico para o teste Hadamard, com espaços reservados para as portas que dependem do operador de evolução temporal.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "2f102112-4ddc-41ea-999c-db5863bc77ac",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Parameters for the template circuits\n",
        "parameters = []\n",
        "for idx in range(1, krylov_dim):\n",
        "    parameters.append(2 * dt_circ * (idx))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "33ec7c29-904e-4445-a654-405214349a4d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/33ec7c29-904e-4445-a654-405214349a4d-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Create modified hadamard test circuit\n",
        "qr = QuantumRegister(n_qubits + 1)\n",
        "qc = QuantumCircuit(qr)\n",
        "qc.h(0)\n",
        "qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)\n",
        "qc.barrier()\n",
        "qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)\n",
        "qc.barrier()\n",
        "qc.x(0)\n",
        "qc.compose(\n",
        "    controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True\n",
        ")\n",
        "qc.x(0)\n",
        "\n",
        "qc.decompose().draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "157356c9-06bb-411c-87ef-cd1d2f73be6f",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The optimized circuit has 2Q gates depth:  74\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    \"The optimized circuit has 2Q gates depth: \",\n",
        "    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5a59bbc4-72c7-4a55-9a9e-57bb7b469021",
      "metadata": {},
      "source": [
        "Reduzimos consideravelmente a profundidade do teste de Hadamard com uma combinação de aproximação de Trotter e unidades não controladas\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ca78bfe3-684f-4d3a-b8c2-3d4e55c9ec30",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "## Passo 3: Execute usando Qiskit primitives\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a082f023-d752-40e4-b693-c4d0da5ba102",
      "metadata": {},
      "source": [
        "Instanciar o backend e definir os parâmetros de tempo de execução\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0d90e4df-e262-4852-a811-ff6a1d3232ae",
      "metadata": {},
      "outputs": [],
      "source": [
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(operational=True, simulator=False)\n",
        "if (\n",
        "    \"if_else\" not in backend.target.operation_names\n",
        "):  # Needed as \"op_name\" could be \"if_else\"\n",
        "    backend.target.add_instruction(IfElseOp, name=\"if_else\")\n",
        "print(backend.name)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b3f7227-e6ef-4a11-93d6-d3696054b9cc",
      "metadata": {},
      "source": [
        "<span id=\"transpiling-to-a-qpu\" />\n",
        "\n",
        "### Transpiling para um QPU\n",
        "\n",
        "Primeiro, vamos selecionar subconjuntos do mapa de acoplamento com qubits de “bom” desempenho (onde “bom” é bastante arbitrário aqui, queremos principalmente evitar qubits de desempenho realmente ruim) e criar um novo alvo para transpilagem\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "e123fda1-6454-4893-9612-2b8591e8cfb9",
      "metadata": {},
      "outputs": [],
      "source": [
        "target = backend.target\n",
        "cmap = target.build_coupling_map(filter_idle_qubits=True)\n",
        "cmap_list = list(cmap.get_edges())\n",
        "\n",
        "cust_cmap_list = copy.deepcopy(cmap_list)\n",
        "for q in range(target.num_qubits):\n",
        "    meas_err = target[\"measure\"][(q,)].error\n",
        "    t2 = target.qubit_properties[q].t2 * 1e6\n",
        "    if meas_err > 0.02 or t2 < 100:\n",
        "        for q_pair in cmap_list:\n",
        "            if q in q_pair:\n",
        "                try:\n",
        "                    cust_cmap_list.remove(q_pair)\n",
        "                except:\n",
        "                    continue\n",
        "\n",
        "for q in cmap_list:\n",
        "    op_name = list(target.operation_names_for_qargs(q))[0]\n",
        "    twoq_gate_err = target[f\"{op_name}\"][q].error\n",
        "    if twoq_gate_err > 0.005:\n",
        "        for q_pair in cmap_list:\n",
        "            if q == q_pair:\n",
        "                try:\n",
        "                    cust_cmap_list.remove(q)\n",
        "                except:\n",
        "                    continue\n",
        "\n",
        "\n",
        "cust_cmap = CouplingMap(cust_cmap_list)\n",
        "cust_target = Target.from_configuration(\n",
        "    basis_gates=backend.configuration().basis_gates,\n",
        "    coupling_map=cust_cmap,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "96f79e35-a265-4de4-b98c-f9c84f58f0d5",
      "metadata": {},
      "source": [
        "Em seguida, transpile o circuito virtual para o melhor layout físico nesse novo destino\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "id": "62109c95-79fe-4076-b591-9bd920fd51f4",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "depth 52\n",
            "num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])\n",
            "physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]\n"
          ]
        }
      ],
      "source": [
        "basis_gates = list(target.operation_names)\n",
        "pm = generate_preset_pass_manager(\n",
        "    optimization_level=3,\n",
        "    target=cust_target,\n",
        "    basis_gates=basis_gates,\n",
        ")\n",
        "\n",
        "qc_trans = pm.run(qc)\n",
        "\n",
        "print(\"depth\", qc_trans.depth(lambda x: x[0].num_qubits == 2))\n",
        "print(\"num 2q ops\", qc_trans.count_ops())\n",
        "print(\n",
        "    \"physical qubits\",\n",
        "    sorted(\n",
        "        [\n",
        "            idx\n",
        "            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()\n",
        "            if qb._register.name != \"ancilla\"\n",
        "        ]\n",
        "    ),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9b9b0170-b96a-47b9-b5ab-5112d02338a3",
      "metadata": {},
      "source": [
        "<span id=\"create-pubs-for-execution-with-estimator\" />\n",
        "\n",
        "### Criar PUBs para execução com o Estimador\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "bcfd6693-d1fb-44e4-9b06-c70e6766e877",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Define observables to measure for S\n",
        "observable_S_real = \"I\" * (n_qubits) + \"X\"\n",
        "observable_S_imag = \"I\" * (n_qubits) + \"Y\"\n",
        "\n",
        "observable_op_real = SparsePauliOp(\n",
        "    observable_S_real\n",
        ")  # define a sparse pauli operator for the observable\n",
        "observable_op_imag = SparsePauliOp(observable_S_imag)\n",
        "\n",
        "layout = qc_trans.layout  # get layout of transpiled circuit\n",
        "observable_op_real = observable_op_real.apply_layout(\n",
        "    layout\n",
        ")  # apply physical layout to the observable\n",
        "observable_op_imag = observable_op_imag.apply_layout(layout)\n",
        "observable_S_real = (\n",
        "    observable_op_real.paulis.to_labels()\n",
        ")  # get the label of the physical observable\n",
        "observable_S_imag = observable_op_imag.paulis.to_labels()\n",
        "\n",
        "observables_S = [[observable_S_real], [observable_S_imag]]\n",
        "\n",
        "\n",
        "# Define observables to measure for H\n",
        "# Hamiltonian terms to measure\n",
        "observable_list = []\n",
        "for pauli, coeff in zip(H_op.paulis, H_op.coeffs):\n",
        "    # print(pauli)\n",
        "    observable_H_real = pauli[::-1].to_label() + \"X\"\n",
        "    observable_H_imag = pauli[::-1].to_label() + \"Y\"\n",
        "    observable_list.append([observable_H_real])\n",
        "    observable_list.append([observable_H_imag])\n",
        "\n",
        "layout = qc_trans.layout\n",
        "\n",
        "observable_trans_list = []\n",
        "for observable in observable_list:\n",
        "    observable_op = SparsePauliOp(observable)\n",
        "    observable_op = observable_op.apply_layout(layout)\n",
        "    observable_trans_list.append([observable_op.paulis.to_labels()])\n",
        "\n",
        "observables_H = observable_trans_list\n",
        "\n",
        "\n",
        "# Define a sweep over parameter values\n",
        "params = np.vstack(parameters).T\n",
        "\n",
        "\n",
        "# Estimate the expectation value for all combinations of\n",
        "# observables and parameter values, where the pub result will have\n",
        "# shape (# observables, # parameter values).\n",
        "pub = (qc_trans, observables_S + observables_H, params)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8bd27e84-d673-4326-9004-f3be5a65aed3",
      "metadata": {},
      "source": [
        "<span id=\"run-circuits\" />\n",
        "\n",
        "### Circuitos de corrida\n",
        "\n",
        "Os circuitos para $t=0$ são calculáveis de forma clássica\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "768cd443-2ed8-4ba2-ba9d-98491b36fa54",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "(25+0j)\n"
          ]
        }
      ],
      "source": [
        "qc_cliff = qc.assign_parameters({t: 0})\n",
        "\n",
        "\n",
        "# Get expectation values from experiment\n",
        "S_expval_real = StabilizerState(qc_cliff).expectation_value(\n",
        "    Pauli(\"I\" * (n_qubits) + \"X\")\n",
        ")\n",
        "S_expval_imag = StabilizerState(qc_cliff).expectation_value(\n",
        "    Pauli(\"I\" * (n_qubits) + \"Y\")\n",
        ")\n",
        "\n",
        "# Get expectation values\n",
        "S_expval = S_expval_real + 1j * S_expval_imag\n",
        "\n",
        "H_expval = 0\n",
        "for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n",
        "    # Get expectation values from experiment\n",
        "    expval_real = StabilizerState(qc_cliff).expectation_value(\n",
        "        Pauli(pauli[::-1].to_label() + \"X\")\n",
        "    )\n",
        "    expval_imag = StabilizerState(qc_cliff).expectation_value(\n",
        "        Pauli(pauli[::-1].to_label() + \"Y\")\n",
        "    )\n",
        "    expval = expval_real + 1j * expval_imag\n",
        "\n",
        "    # Fill-in matrix elements\n",
        "    H_expval += coeff * expval\n",
        "\n",
        "\n",
        "print(H_expval)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "794edf4c-d539-4864-865d-52a88a450b5c",
      "metadata": {},
      "source": [
        "Executar circuitos para $S$ e $\\tilde{H}$ com o Estimator\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "9059dc9e-9203-4572-92c4-76e9c85dfed9",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Experiment options\n",
        "num_randomizations = 300\n",
        "num_randomizations_learning = 30\n",
        "shots_per_randomization = 100\n",
        "noise_factors = [1, 1.2, 1.4]\n",
        "learning_pair_depths = [0, 4, 24, 48]\n",
        "\n",
        "\n",
        "experimental_opts = {}\n",
        "experimental_opts[\"resilience\"] = {\n",
        "    \"measure_mitigation\": True,\n",
        "    \"measure_noise_learning\": {\n",
        "        \"num_randomizations\": num_randomizations_learning,\n",
        "        \"shots_per_randomization\": shots_per_randomization,\n",
        "    },\n",
        "    \"zne_mitigation\": True,\n",
        "    \"zne\": {\"noise_factors\": noise_factors},\n",
        "    \"layer_noise_learning\": {\n",
        "        \"max_layers_to_learn\": 10,\n",
        "        \"layer_pair_depths\": learning_pair_depths,\n",
        "        \"shots_per_randomization\": shots_per_randomization,\n",
        "        \"num_randomizations\": num_randomizations_learning,\n",
        "    },\n",
        "    \"zne\": {\n",
        "        \"amplifier\": \"pea\",\n",
        "        \"extrapolated_noise_factors\": [0] + noise_factors,\n",
        "    },\n",
        "}\n",
        "experimental_opts[\"twirling\"] = {\n",
        "    \"num_randomizations\": num_randomizations,\n",
        "    \"shots_per_randomization\": shots_per_randomization,\n",
        "    \"strategy\": \"all\",\n",
        "}\n",
        "\n",
        "estimator = Estimator(mode=backend, options=experimental_opts)\n",
        "\n",
        "\n",
        "job = estimator.run([pub])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9524e609-b457-42e9-9de7-8dbd4464ac74",
      "metadata": {},
      "source": [
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "## Etapa 4: Pós-processamento e retorno do resultado no formato clássico desejado\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 46,
      "id": "30d72032-02e4-40e7-af28-cadda2f7f3bd",
      "metadata": {},
      "outputs": [],
      "source": [
        "results = job.result()[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8149d431-f069-4567-a0be-0cc4ed8d523b",
      "metadata": {},
      "source": [
        "<span id=\"calculate-effective-hamiltonian-and-overlap-matrices\" />\n",
        "\n",
        "### Calcule as matrizes Hamiltoniana Eficaz e de Sobreposição\n",
        "\n",
        "Primeiro, calcule a fase acumulada pelo estado $\\vert 0 \\rangle$ durante a evolução do tempo sem controle\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 47,
      "id": "34ba8b05-d52c-4b0f-8d28-acaf7a5ddffc",
      "metadata": {},
      "outputs": [],
      "source": [
        "prefactors = [\n",
        "    np.exp(-1j * sum([c for p, c in H_op.to_list() if \"Z\" in p]) * i * dt)\n",
        "    for i in range(1, krylov_dim)\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "36875c8b-36c8-464b-9f0e-7a79e17d5fef",
      "metadata": {},
      "source": [
        "Quando tivermos os resultados das execuções do circuito, poderemos pós-processar os dados para calcular os elementos da matriz de $S$\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 48,
      "id": "16dd2534-8d5f-40f7-a938-5d519475fd8d",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Assemble S, the overlap matrix of dimension D:\n",
        "S_first_row = np.zeros(krylov_dim, dtype=complex)\n",
        "S_first_row[0] = 1 + 0j\n",
        "\n",
        "# Add in ancilla-only measurements:\n",
        "for i in range(krylov_dim - 1):\n",
        "    # Get expectation values from experiment\n",
        "    expval_real = results.data.evs[0][0][\n",
        "        i\n",
        "    ]  # automatic extrapolated evs if ZNE is used\n",
        "    expval_imag = results.data.evs[1][0][\n",
        "        i\n",
        "    ]  # automatic extrapolated evs if ZNE is used\n",
        "\n",
        "    # Get expectation values\n",
        "    expval = expval_real + 1j * expval_imag\n",
        "    S_first_row[i + 1] += prefactors[i] * expval\n",
        "\n",
        "S_first_row_list = S_first_row.tolist()  # for saving purposes\n",
        "\n",
        "\n",
        "S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n",
        "\n",
        "# Distribute entries from first row across matrix:\n",
        "for i, j in it.product(range(krylov_dim), repeat=2):\n",
        "    if i >= j:\n",
        "        S_circ[j, i] = S_first_row[i - j]\n",
        "    else:\n",
        "        S_circ[j, i] = np.conj(S_first_row[j - i])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 49,
      "id": "01c6563d-87ca-4bd6-b487-dcdaece2d8c2",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$$\n",
              "\\displaystyle \\left[\\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\\end{matrix}\\right]\n",
              "$$"
            ],
            "text/plain": [
              "Matrix([\n",
              "[                                     1.0, -0.723052998582984 - 0.345085413575966*I,  0.467051960502366 + 0.516197865254034*I, -0.180546747798251 - 0.492624093654174*I, 0.0012070853532697 + 0.312052218182462*I],\n",
              "[-0.723052998582984 + 0.345085413575966*I,                                      1.0, -0.723052998582984 - 0.345085413575966*I,  0.467051960502366 + 0.516197865254034*I, -0.180546747798251 - 0.492624093654174*I],\n",
              "[ 0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I,                                      1.0, -0.723052998582984 - 0.345085413575966*I,  0.467051960502366 + 0.516197865254034*I],\n",
              "[-0.180546747798251 + 0.492624093654174*I,  0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I,                                      1.0, -0.723052998582984 - 0.345085413575966*I],\n",
              "[0.0012070853532697 - 0.312052218182462*I, -0.180546747798251 + 0.492624093654174*I,  0.467051960502366 - 0.516197865254034*I, -0.723052998582984 + 0.345085413575966*I,                                      1.0]])"
            ]
          },
          "execution_count": 49,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Matrix(S_circ)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a34c72c0-2984-47fb-8579-089ea795ef39",
      "metadata": {},
      "source": [
        "E os elementos da matriz de $\\tilde{H}$\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 50,
      "id": "9cde2419-8a29-4d7a-b5b1-9ef2d4c126d2",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Assemble S, the overlap matrix of dimension D:\n",
        "H_first_row = np.zeros(krylov_dim, dtype=complex)\n",
        "H_first_row[0] = H_expval\n",
        "\n",
        "for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):\n",
        "    # Add in ancilla-only measurements:\n",
        "    for i in range(krylov_dim - 1):\n",
        "        # Get expectation values from experiment\n",
        "        expval_real = results.data.evs[2 + 2 * obs_idx][0][\n",
        "            i\n",
        "        ]  # automatic extrapolated evs if ZNE is used\n",
        "        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][\n",
        "            i\n",
        "        ]  # automatic extrapolated evs if ZNE is used\n",
        "\n",
        "        # Get expectation values\n",
        "        expval = expval_real + 1j * expval_imag\n",
        "        H_first_row[i + 1] += prefactors[i] * coeff * expval\n",
        "\n",
        "H_first_row_list = H_first_row.tolist()\n",
        "\n",
        "H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)\n",
        "\n",
        "# Distribute entries from first row across matrix:\n",
        "for i, j in it.product(range(krylov_dim), repeat=2):\n",
        "    if i >= j:\n",
        "        H_eff_circ[j, i] = H_first_row[i - j]\n",
        "    else:\n",
        "        H_eff_circ[j, i] = np.conj(H_first_row[j - i])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 51,
      "id": "5800d892-5554-4a92-bfe0-3750ef0a3ed7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$$\n",
              "\\displaystyle \\left[\\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\\end{matrix}\\right]\n",
              "$$"
            ],
            "text/plain": [
              "Matrix([\n",
              "[                                  25.0, -14.2437089383409 - 6.50486277982165*I,   10.2857217968584 + 9.0431912203186*I, -5.15587257589417 - 8.88280836036843*I,   1.98818301405581 + 5.8897614762563*I],\n",
              "[-14.2437089383409 + 6.50486277982165*I,                                   25.0, -14.2437089383409 - 6.50486277982165*I,   10.2857217968584 + 9.0431912203186*I, -5.15587257589417 - 8.88280836036843*I],\n",
              "[  10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I,                                   25.0, -14.2437089383409 - 6.50486277982165*I,   10.2857217968584 + 9.0431912203186*I],\n",
              "[-5.15587257589417 + 8.88280836036843*I,   10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I,                                   25.0, -14.2437089383409 - 6.50486277982165*I],\n",
              "[  1.98818301405581 - 5.8897614762563*I, -5.15587257589417 + 8.88280836036843*I,   10.2857217968584 - 9.0431912203186*I, -14.2437089383409 + 6.50486277982165*I,                                   25.0]])"
            ]
          },
          "execution_count": 51,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Matrix(H_eff_circ)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "896553fa-3c28-44cc-89ad-5031d2f8d35f",
      "metadata": {},
      "source": [
        "Por fim, podemos resolver o problema do valor próprio generalizado para $\\tilde{H}$ :\n",
        "\n",
        "$\\tilde{H} \\vec{c} = c S \\vec{c}$\n",
        "\n",
        "e obter uma estimativa da energia do estado fundamental $c_{min}$\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 58,
      "id": "8b997d15-ee25-40eb-80e8-1f054a298ee9",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The estimated ground state energy is:  25.0\n",
            "The estimated ground state energy is:  22.572154819954875\n",
            "The estimated ground state energy is:  21.691509219286587\n",
            "The estimated ground state energy is:  21.23882298756386\n",
            "The estimated ground state energy is:  20.965499325470294\n"
          ]
        }
      ],
      "source": [
        "gnd_en_circ_est_list = []\n",
        "for d in range(1, krylov_dim + 1):\n",
        "    # Solve generalized eigenvalue problem for different size of the Krylov space\n",
        "    gnd_en_circ_est = solve_regularized_gen_eig(\n",
        "        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1\n",
        "    )\n",
        "    gnd_en_circ_est_list.append(gnd_en_circ_est)\n",
        "    print(\"The estimated ground state energy is: \", gnd_en_circ_est)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b39943f8-b75f-4bd4-b246-833ea4f9b710",
      "metadata": {},
      "source": [
        "Para um setor de partícula única, podemos calcular com eficiência o estado fundamental desse setor do Hamiltoniano de forma clássica\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 59,
      "id": "fcfe07e5-99e4-4276-a4d6-6f1c7e28c5f2",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "n_sys_qubits 30\n",
            "n_exc 1 , subspace dimension 31\n",
            "single particle ground state energy:  21.021912418526906\n"
          ]
        }
      ],
      "source": [
        "gs_en = single_particle_gs(H_op, n_qubits)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 60,
      "id": "4bc52594-0376-497f-8a61-0949415a1fe0",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/krylov-quantum-diagonalization/extracted-outputs/4bc52594-0376-497f-8a61-0949415a1fe0-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(\n",
        "    range(1, krylov_dim + 1),\n",
        "    gnd_en_circ_est_list,\n",
        "    color=\"blue\",\n",
        "    linestyle=\"-.\",\n",
        "    label=\"KQD estimate\",\n",
        ")\n",
        "plt.plot(\n",
        "    range(1, krylov_dim + 1),\n",
        "    [gs_en] * krylov_dim,\n",
        "    color=\"red\",\n",
        "    linestyle=\"-\",\n",
        "    label=\"exact\",\n",
        ")\n",
        "plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))\n",
        "plt.legend()\n",
        "plt.xlabel(\"Krylov space dimension\")\n",
        "plt.ylabel(\"Energy\")\n",
        "plt.title(\n",
        "    \"Estimating Ground state energy with Krylov Quantum Diagonalization\"\n",
        ")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2ef2fac9-b4b9-4457-a551-c98999a88710",
      "metadata": {},
      "source": [
        "<span id=\"appendix-krylov-subspace-from-real-time-evolutions\" />\n",
        "\n",
        "## Apêndice: Subespaço de Krylov a partir de evoluções em tempo real\n",
        "\n",
        "O espaço unitário de Krylov é definido como\n",
        "\n",
        "$$\n",
        "\\mathcal{K}_U(H, |\\psi\\rangle) = \\text{span}\\left\\{ |\\psi\\rangle,  e^{-iH\\,dt} |\\psi\\rangle, \\dots, e^{-irH\\,dt} |\\psi\\rangle \\right\\}\n",
        "$$\n",
        "\n",
        "para algum intervalo de tempo $dt$ que determinaremos posteriormente. Suponha temporariamente que $r$ seja par: em seguida, defina $d=r/2$. Observe que, quando projetamos o Hamiltoniano no espaço de Krylov acima, ele é indistinguível do espaço de Krylov\n",
        "\n",
        "$$\n",
        "\\mathcal{K}_U(H, |\\psi\\rangle) = \\text{span}\\left\\{ e^{i\\,d\\,H\\,dt}|\\psi\\rangle,  e^{i(d-1)H\\,dt} |\\psi\\rangle, \\dots, e^{-i(d-1)H\\,dt} |\\psi\\rangle, e^{-i\\,d\\,H\\,dt} |\\psi\\rangle \\right\\},\n",
        "$$\n",
        "\n",
        "ou seja, onde todas as evoluções temporais são deslocadas para trás em $d$ passos de tempo.\n",
        "O motivo pelo qual é indistinguível é que os elementos da matriz\n",
        "\n",
        "$$\n",
        "\\tilde{H}_{j,k} = \\langle\\psi|e^{i\\,j\\,H\\,dt}He^{-i\\,k\\,H\\,dt}|\\psi\\rangle=\\langle\\psi|He^{i(j-k)H\\,dt}|\\psi\\rangle\n",
        "$$\n",
        "\n",
        "são invariantes sob mudanças gerais do tempo de evolução, uma vez que as evoluções temporais são compatíveis com o Hamiltoniano. Para $r$ ímpares, podemos usar a análise para $r-1$.\n",
        "\n",
        "Queremos mostrar que, em algum lugar desse espaço de Krylov, é garantido que haja um estado de baixa energia. Fazemos isso por meio do seguinte resultado, que é derivado do Teorema 3.1 em [\\[3\\]](#references) :\n",
        "\n",
        "**Afirmação 1:** existe uma função $f$ de modo que, para as energias $E$ na faixa espectral do Hamiltoniano (ou seja, entre a energia do estado fundamental e a energia máxima)...\n",
        "\n",
        "1. $f(E_0)=1$\n",
        "2. $|f(E)|\\le2\\left(1 + \\delta\\right)^{-d}$ para todos os valores de $E$ que se encontram a uma distância de $\\ge\\delta$ de $E_0$, ou seja, é exponencialmente suprimido\n",
        "3. $f(E)$ é uma combinação linear de $e^{ijE\\,dt}$ para $j=-d,-d+1,...,d-1,d$\n",
        "\n",
        "Apresentamos uma prova abaixo, mas ela pode ser ignorada com segurança, a menos que se queira entender o argumento completo e rigoroso. Por enquanto, vamos nos concentrar nas implicações da afirmação acima. Pela propriedade 3 acima, podemos ver que o espaço de Krylov deslocado acima contém o estado $f(H)|\\psi\\rangle$. Esse é o nosso estado de baixa energia. Para entender o motivo, escreva $|\\psi\\rangle$ na base de energia própria:\n",
        "\n",
        "$$\n",
        "|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_k|E_k\\rangle,\n",
        "$$\n",
        "\n",
        "em que $|E_k\\rangle$ é o k-ésimo estado próprio de energia e $\\gamma_k$ é sua amplitude no estado inicial $|\\psi\\rangle$. Expresso em termos disso, $f(H)|\\psi\\rangle$ é dado por\n",
        "\n",
        "$$\n",
        "f(H)|\\psi\\rangle = \\sum_{k=0}^{N}\\gamma_kf(E_k)|E_k\\rangle,\n",
        "$$\n",
        "\n",
        "usando o fato de que podemos substituir $H$ por $E_k$ quando ele atua no estado próprio $|E_k\\rangle$. O erro de energia desse estado é, portanto\n",
        "\n",
        "$$\n",
        "\\text{energy error} = \\frac{\\langle\\psi|f(H)(H-E_0)f(H)|\\psi\\rangle}{\\langle\\psi|f(H)^2|\\psi\\rangle}\n",
        "$$\n",
        "\n",
        "$$\n",
        "= \\frac{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2}.\n",
        "$$\n",
        "\n",
        "Para transformar isso em um limite superior que seja mais fácil de entender, primeiro separamos a soma no numerador em termos com $E_k-E_0\\le\\delta$ e termos com $E_k-E_0>\\delta$ :\n",
        "\n",
        "$$\n",
        "\\text{energy error} = \\frac{\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} + \\frac{\\sum_{E_k> E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2}.\n",
        "$$\n",
        "\n",
        "Podemos limitar o primeiro termo por $\\delta$,\n",
        "\n",
        "$$\n",
        "\\frac{\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0)}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} < \\frac{\\delta\\sum_{E_k\\le E_0+\\delta}|\\gamma_k|^2f(E_k)^2}{\\sum_{k=0}^{N}|\\gamma_k|^2f(E_k)^2} \\le \\delta,\n",
        "$$\n",
        "\n",
        "em que a primeira etapa segue porque $E_k-E_0\\le\\delta$ para cada $E_k$ na soma, e a segunda etapa segue porque a soma no numerador é um subconjunto da soma no denominador. Para o segundo termo, primeiro reduzimos o denominador em $|\\gamma_0|^2$, já que $f(E_0)^2=1$ : somando tudo, obtemos\n",
        "\n",
        "$$\n",
        "\\text{energy error} \\le \\delta + \\frac{1}{|\\gamma_0|^2}\\sum_{E_k>E_0+\\delta}|\\gamma_k|^2f(E_k)^2(E_k-E_0).\n",
        "$$\n",
        "\n",
        "Para simplificar o que resta, observe que, para todos esses $E_k$, pela definição de $f$, sabemos que $f(E_k)^2 \\le 4\\left(1 + \\delta\\right)^{-2d}$. Além disso, o limite superior de $E_k-E_0<2\\|H\\|$ e o limite superior de $\\sum_{E_k>E_0+\\delta}|\\gamma_k|^2<1$ fornecem\n",
        "\n",
        "$$\n",
        "\\text{energy error} \\le \\delta + \\frac{8}{|\\gamma_0|^2}\\|H\\|\\left(1 + \\delta\\right)^{-2d}.\n",
        "$$\n",
        "\n",
        "Isso é válido para qualquer $\\delta>0$, portanto, se definirmos $\\delta$ igual ao nosso erro de meta, então o limite de erro acima converge exponencialmente para esse valor com a dimensão de Krylov $2d=r$. Observe também que, se $\\delta<E_1-E_0$, o termo $\\delta$ desaparece totalmente no limite acima.\n",
        "\n",
        "Para completar o argumento, primeiro observamos que o resultado acima é apenas o erro de energia do estado específico $f(H)|\\psi\\rangle$, em vez do erro de energia do estado de menor energia no espaço de Krylov. No entanto, pelo princípio variacional (Rayleigh-Ritz), o erro de energia do estado de menor energia no espaço de Krylov é limitado superiormente pelo erro de energia de qualquer estado no espaço de Krylov, de modo que o acima exposto também é um limite superior do erro de energia do estado de menor energia, ou seja, o resultado do algoritmo de diagonalização quântica de Krylov.\n",
        "\n",
        "Uma análise semelhante à anterior pode ser realizada, levando em conta, adicionalmente, o ruído e o procedimento de limiar discutido no bloco de notas. Consulte [\\[2\\]](#references) e [\\[4\\]](#references) para ver essa análise.\n",
        "\n",
        "<span id=\"appendix-proof-of-claim-1\" />\n",
        "\n",
        "## Apêndice: prova da alegação 1\n",
        "\n",
        "O seguinte é derivado principalmente de [\\[3\\]](#references), Teorema 3.1:. Seja $0 < a < b$ e $\\Pi^*_d$ o espaço de polinômios residuais (polinômios cujo valor em 0 é 1) de grau no máximo $d$. A solução para\n",
        "\n",
        "$$\n",
        "\\beta(a, b, d) = \\min_{p \\in \\Pi^*_d} \\max_{x \\in [a, b]} |p(x)| \\quad\n",
        "$$\n",
        "\n",
        "é\n",
        "\n",
        "$$\n",
        "p^*(x) = \\frac{T_d\\left(\\frac{b + a - 2x}{b - a}\\right)}{T_d\\left(\\frac{b + a}{b - a}\\right)}, \\quad\n",
        "$$\n",
        "\n",
        "e o valor mínimo correspondente é\n",
        "\n",
        "$$\n",
        "\\beta(a, b, d) = T_d^{-1}\\left(\\frac{b + a}{b - a}\\right).\n",
        "$$\n",
        "\n",
        "Queremos converter isso em uma função que possa ser expressa naturalmente em termos de exponenciais complexos, porque essas são as evoluções de tempo reais que geram o espaço de Krylov quântico.\n",
        "Para isso, é conveniente introduzir a seguinte transformação de energias dentro do intervalo espectral do Hamiltoniano para números no intervalo $[0,1]$ : definir\n",
        "\n",
        "$$\n",
        "g(E) = \\frac{1-\\cos\\big((E-E_0)dt\\big)}{2},\n",
        "$$\n",
        "\n",
        "onde $dt$ é um intervalo de tempo tal que $-\\pi < E_0dt < E_\\text{max}dt < \\pi$. Observe que $g(E_0)=0$ e $g(E)$ crescem à medida que $E$ se afasta de $E_0$.\n",
        "\n",
        "Agora, usando o polinômio $p^*(x)$ com os parâmetros a, b, d definidos como $a = g(E_0 + \\delta)$, $b = 1$ e d = int( r/2 ), definimos a função:\n",
        "\n",
        "$$\n",
        "f(E) = p^* \\left( g(E) \\right) = \\frac{T_d\\left(1 + 2\\frac{\\cos\\big((E-E_0)dt\\big) - \\cos\\big(\\delta\\,dt\\big)}{1 +\\cos\\big(\\delta\\,dt\\big)}\\right)}{T_d\\left(1 + 2\\frac{1-\\cos\\big(\\delta\\,dt\\big)}{1 + \\cos\\big(\\delta\\,dt\\big)}\\right)}\n",
        "$$\n",
        "\n",
        "em que $E_0$ é a energia do estado fundamental. Podemos ver, inserindo $\\cos(x)=\\frac{e^{ix}+e^{-ix}}{2}$, que $f(E)$ é um polinômio trigonométrico de grau $d$, ou seja, uma combinação linear de $e^{ijE\\,dt}$ para $j=-d,-d+1,...,d-1,d$. Além disso, a partir da definição de $p^*(x)$ acima, temos que $f(E_0)=p(0)=1$ e para qualquer $E$ na faixa espectral tal que $\\vert E-E_0 \\vert > \\delta$ temos\n",
        "\n",
        "$$\n",
        "|f(E)| \\le \\beta(a, b, d) = T_d^{-1}\\left(1 + 2\\frac{1-\\cos\\big(\\delta\\,dt\\big)}{1 + \\cos\\big(\\delta\\,dt\\big)}\\right)\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\leq 2\\left(1 + \\delta\\right)^{-d} = 2\\left(1 + \\delta\\right)^{-\\lfloor k/2\\rfloor}.\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ccd6c21d-ed1e-4894-9cbb-8d3a9d91afe2",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## Referências\n",
        "\n",
        "\\[1] N. Yoshioka, M. Amico, W. Kirby et al. \"Diagonalization of large many-body Hamiltonians on a quantum processor\" (Diagonalização de Hamiltonianos de muitos corpos grandes em um processador quântico). [arXiv:2407.14431](https://arxiv.org/abs/2407.14431)\n",
        "\n",
        "\\[2] Ethan N. Epperly, Lin Lin e Yuji Nakatsukasa. \"Uma teoria da diagonalização do subespaço quântico\". SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).\n",
        "\n",
        "\\[3] Å. Björck. \"Métodos numéricos em cálculos de matrizes\". Textos em Matemática Aplicada. Springer International Publishing. (2014).\n",
        "\n",
        "\\[4] William Kirby. \"Análise de algoritmos de Krylov quânticos com erros\". Quantum 8, 1457 (2024).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "74948cc7-041f-412c-ab16-2554bf164061",
      "metadata": {},
      "source": [
        "<span id=\"tutorial-survey\" />\n",
        "\n",
        "## Pesquisa tutorial\n",
        "\n",
        "Responda a esta breve pesquisa para fornecer feedback sobre este tutorial. Suas percepções nos ajudarão a melhorar nossas ofertas de conteúdo e a experiência do usuário.\n",
        "\n",
        "[Link para a pesquisa](https://your.feedback.ibm.com/jfe/form/SV_82nennpKIjjD8rQ)\n",
        "\n"
      ]
    },
    {
      "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"
    },
    "hours": 1.5,
    "qpuSeconds": 1200
  },
  "nbformat": 4,
  "nbformat_minor": 4
}