From 3cd9338f7acc876786150305633c7d3f45f1d8b1 Mon Sep 17 00:00:00 2001 From: Sam-XiaoyueLi Date: Fri, 23 Feb 2024 17:49:47 +0800 Subject: [PATCH] Fix error: missing argument d in gradient_onsite_Z --- examples/dbi/dbi_strategies_compare.ipynb | 8 + examples/dbi/dbi_strategy_Pauli-Z 2.ipynb | 514 ++++++++++++++++++ .../dbi/dbi_strategy_magnetic_field.ipynb | 112 +++- src/qibo/models/dbi/utils.py | 21 +- tests/test_models_dbi_utils 2.py | 50 ++ 5 files changed, 687 insertions(+), 18 deletions(-) create mode 100644 examples/dbi/dbi_strategy_Pauli-Z 2.ipynb create mode 100644 tests/test_models_dbi_utils 2.py diff --git a/examples/dbi/dbi_strategies_compare.ipynb b/examples/dbi/dbi_strategies_compare.ipynb index 0afd3a9841..0bb5642abf 100644 --- a/examples/dbi/dbi_strategies_compare.ipynb +++ b/examples/dbi/dbi_strategies_compare.ipynb @@ -194,6 +194,7 @@ "steps_gradient_plot= [0]\n", "for _ in range(NSTEPS):\n", " step, d_coef, d = gradient_descent_onsite_Z(dbi_gradient, d_coef, d, onsite_Z_ops=onsite_Z_ops, max_evals=50, n_taylor=5)\n", + " dbi_gradient(d=d,step=step)\n", " off_diagonal_norm_history_gradient.append(dbi_gradient.off_diagonal_norm)\n", " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with d_coef {d_coef}, loss {dbi_gradient.off_diagonal_norm}\")\n", " steps_gradient_plot.append(steps_gradient_plot[-1]+step)" @@ -228,6 +229,13 @@ "plt.xlabel('Iteration')\n", "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/dbi/dbi_strategy_Pauli-Z 2.ipynb b/examples/dbi/dbi_strategy_Pauli-Z 2.ipynb new file mode 100644 index 0000000000..5f563a172a --- /dev/null +++ b/examples/dbi/dbi_strategy_Pauli-Z 2.ipynb @@ -0,0 +1,514 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Double-Bracket Iteration Strategy: Pauli-Z products\n", + "\n", + "In this example, we demonstrate the usage of a DBI strategy, where the diagonal operators for double bracket iterations are variationally chosen from all possible local Pauli-Z operators." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Initial setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "!python -m pip install hyperopt # required to optimize the DBF step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from copy import copy, deepcopy\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from qibo import hamiltonians, set_backend\n", + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketIteration\n", + "from qibo.models.dbi.utils import *" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below are some useful functions to visualize the diagonalization process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def visualize_matrix(matrix, title=\"\"):\n", + " \"\"\"Visualize absolute values of a matrix in a heatmap form.\"\"\"\n", + " fig, ax = plt.subplots(figsize=(5, 5))\n", + " ax.set_title(title)\n", + " try:\n", + " im = ax.imshow(np.absolute(matrix), cmap=\"inferno\")\n", + " except TypeError:\n", + " im = ax.imshow(np.absolute(matrix.get()), cmap=\"inferno\")\n", + " fig.colorbar(im, ax=ax)\n", + "\n", + "\n", + "def visualize_drift(h0, h):\n", + " \"\"\"Visualize drift of the evolved hamiltonian w.r.t. h0.\"\"\"\n", + " fig, ax = plt.subplots(figsize=(5, 5))\n", + " ax.set_title(r\"Drift: $|\\hat{H}_0 - \\hat{H}_{1}|$\")\n", + " try:\n", + " im = ax.imshow(np.absolute(h0 - h), cmap=\"inferno\")\n", + " except TypeError:\n", + " im = ax.imshow(np.absolute((h0 - h).get()), cmap=\"inferno\")\n", + "\n", + " fig.colorbar(im, ax=ax)\n", + "\n", + "\n", + "def plot_histories(loss_histories: list, steps: list, labels: list = None):\n", + " \"\"\"Plot off-diagonal norm histories over a sequential evolution.\"\"\"\n", + " plt.figure(figsize=(5, 5 * 6 / 8))\n", + " if len(steps) == 1:\n", + " # fixed_step\n", + " x_axis = [i * steps[0] for i in range(len(loss_histories))]\n", + " else:\n", + " x_axis = [sum(steps[:k]) for k in range(1, len(steps) + 1)]\n", + " plt.plot(x_axis, loss_histories, \"-o\")\n", + "\n", + " x_labels_rounded = [round(x, 2) for x in x_axis]\n", + " x_labels_rounded = [0] + x_labels_rounded[0:5] + [max(x_labels_rounded)]\n", + " x_labels_rounded.pop(3)\n", + " plt.xticks(x_labels_rounded)\n", + "\n", + " y_labels_rounded = [round(y, 1) for y in loss_histories]\n", + " y_labels_rounded = y_labels_rounded[0:5] + [min(y_labels_rounded)]\n", + " plt.yticks(y_labels_rounded)\n", + "\n", + " if labels is not None:\n", + " labels_copy = copy(labels)\n", + " labels_copy.insert(0, \"Initial\")\n", + " for i, label in enumerate(labels_copy):\n", + " plt.text(x_axis[i], loss_histories[i], label)\n", + "\n", + " plt.grid()\n", + " plt.xlabel(r\"Flow duration $s$\")\n", + " plt.title(\"Loss function histories\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example: TFIM\n", + "\n", + "As an example, we consider the Transverse Field Ising Model (TFIM):\n", + "$$ H_{\\rm TFIM} = - \\sum_{i=1}^{N}\\bigl( Z_i Z_{i+1} + h X_i \\bigr),$$\n", + "which is already implemented in `Qibo`. For this tutorial we set $N=5$ and $h=3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# set the qibo backend (we suggest qibojit if N >= 20)\n", + "# alternatives: tensorflow (not optimized), numpy (when CPU not supported by jit)\n", + "set_backend(\"qibojit\", \"numba\")\n", + "\n", + "# hamiltonian parameters\n", + "nqubits = 5\n", + "h = 3\n", + "\n", + "# define the hamiltonian\n", + "H_TFIM = hamiltonians.TFIM(nqubits=nqubits, h=h)\n", + "\n", + "# initialize class\n", + "# Note: use deepcopy to prevent h being edited\n", + "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator)\n", + "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(H_TFIM.matrix)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generate local Pauli-Z operators" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "generate_local_Z = generate_Z_operators(nqubits)\n", + "Z_ops = list(generate_local_Z.values())\n", + "Z_names = list(generate_local_Z.keys())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Iteration from a list of operators\n", + "The idea of this strategy is to chose the Z operator that reduces the off-diagonal norm of the hamiltonian most efficiently. Given a list of operators (np.array), the function `select_best_dbr_generator_and_run` searches for the maximum decrease in off-diagonal norm for each operator and runs one double bracket rotation using the optimal operator from the list.\n", + "\n", + "Note that the hyperopt settings can be set as positional arguments." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NSTEPS = 15\n", + "max_evals = 100\n", + "step_max = 1\n", + "Z_optimal = []\n", + "# add in initial values for plotting\n", + "off_diagonal_norm_history = [dbi.off_diagonal_norm]\n", + "steps = [0]\n", + "scheduling = DoubleBracketScheduling.use_hyperopt\n", + "for _ in range(NSTEPS):\n", + " dbi, idx, step, flip_sign = select_best_dbr_generator(dbi, Z_ops, scheduling=scheduling, compare_canonical=False, max_evals=max_evals, step_max=step_max)\n", + " off_diagonal_norm_history.append(dbi.off_diagonal_norm)\n", + " steps.append(steps[-1]+step)\n", + " if flip_sign < 0:\n", + " Z_optimal.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with operator {Z_optimal[-1]}, loss {dbi.off_diagonal_norm}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_histories(off_diagonal_norm_history, steps, Z_optimal)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is worth noting that due to the nature of `hyperopt`, the iterations may be unstable and multiple runs may be required for the optimal result (alternatively, we can perform a grid search on the optimal step). Hence, it is sometimes needed to adjust its parameters including the following:\n", + "\n", + "- step_min\n", + "- step_max\n", + "- max_evals" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compare with canonical\n", + "\n", + "We compare the effectiveness at diagonalzation between the Pauli-Z operators and the canonical generator:\n", + "\n", + "$$ d = [H,\\sigma(H)]$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# set the qibo backend (we suggest qibojit if N >= 20)\n", + "# alternatives: tensorflow (not optimized), numpy (when CPU not supported by jit)\n", + "set_backend(\"qibojit\", \"numba\")\n", + "\n", + "\n", + "# initialize class|\n", + "# Note: use deepcopy to prevent h being edited\n", + "dbi_canonical = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical)\n", + "print(\"Initial off diagonal norm\", dbi_canonical.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "off_diagonal_norm_history_canonical = [dbi_canonical.off_diagonal_norm]\n", + "steps_canonical = [0]\n", + "steps_canonical_plot = [0]\n", + "for s in range(NSTEPS):\n", + " # same settings as iteration from list\n", + " step = dbi_canonical.hyperopt_step(\n", + " step_min = 1e-5,\n", + " step_max = 1,\n", + " space = hp.uniform,\n", + " optimizer = tpe,\n", + " )\n", + " dbi_canonical(step=step)\n", + " print(f\"New optimized step at iteration {s+1}/{NSTEPS}: {step}, loss {dbi_canonical.off_diagonal_norm}\")\n", + " off_diagonal_norm_history_canonical.append(dbi_canonical.off_diagonal_norm)\n", + " steps_canonical.append(step)\n", + " steps_canonical_plot.append(steps_canonical_plot[-1]+step)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "# plt.plot(steps, off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "# plt.plot(steps_canonical, off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.plot(off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "plt.plot(off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.xlabel(\"Iterations\")\n", + "plt.ylabel(\"Norm off-diagonal restriction\")\n", + "plt.title(\"Compare Variational Pauli-Z with Canonical\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(off_diagonal_norm_history)\n", + "print(off_diagonal_norm_history_canonical)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, we make 2 observations:\n", + "\n", + "1. The canonical strategy has a steeper decrease at the beginning than Pauli-Z operators.\n", + "2. However, the canonical strategy is also prone to getting stuck at a local minimum and hence resultting in a lesser degree of diagonalization.\n", + "\n", + "Therefore, we explore the possibility of mixing the two strategies by including the canonical generator in the list." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mixed strategy: optimal at each step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi_mixed = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator)\n", + "print(\"Initial off diagonal norm\", dbi_mixed.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi_eval = deepcopy(dbi_mixed)\n", + "dbi_eval.mode = DoubleBracketGeneratorType.canonical\n", + "if step is None:\n", + " step = dbi_eval.hyperopt_step(\n", + " step_max=step_max,\n", + " space=hp.uniform,\n", + " optimizer=tpe,\n", + " max_evals=max_evals,\n", + " )\n", + "dbi_eval(step=step)\n", + "print('canonical norm', dbi_eval.off_diagonal_norm, 'step', step)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Z_optimal_mixed = []\n", + "# add in initial values for plotting\n", + "off_diagonal_norm_history_mixed = [dbi_mixed.off_diagonal_norm]\n", + "steps = [0]\n", + "for _ in range(NSTEPS):\n", + " dbi_mixed, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed, Z_ops, scheduling=scheduling, compare_canonical=True, max_evals=max_evals, step_max=step_max)\n", + " off_diagonal_norm_history_mixed.append(dbi_mixed.off_diagonal_norm)\n", + " steps.append(steps[-1]+step)\n", + " if idx == len(Z_ops):\n", + " Z_optimal_mixed.append('Canonical')\n", + " elif flip_sign < 0:\n", + " Z_optimal_mixed.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal_mixed.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with operator {Z_optimal_mixed[-1]}, loss {dbi_mixed.off_diagonal_norm}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "# plt.plot(steps, off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "# plt.plot(steps_canonical, off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.plot(off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "plt.plot(off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.plot(off_diagonal_norm_history_mixed, label=\"Mixed\")\n", + "plt.xlabel(\"Iterations\")\n", + "plt.ylabel(\"Norm off-diagonal restriction\")\n", + "plt.title(\"Compare Variational Pauli-Z with Canonical\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After a few tests, we realize that the mixed strategy does not always outperform just using Pauli-Z operators. This could be caused by 2 reasons: \n", + "\n", + "1. Unstability of hyperopt\n", + "2. Tendency of canonical operator to get stuck at a near local minimum" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mixed strategy: initial canonical\n", + "\n", + "Since the canonical double bracket iteration performs better at the initial steps, we attempt to combine the two strategies: iterate a few steps using the canonical bracket before switching to the variational Z-operators." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi_mixed_can= DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical)\n", + "print(\"Initial off diagonal norm\", dbi_mixed_can.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Run the initial iterations using canonical iterations\n", + "off_diagonal_norm_history_mixed_can = [dbi_mixed_can.off_diagonal_norm]\n", + "steps_mixed_can = [0]\n", + "cannonical_NSTEPS = 2\n", + "for i in range(cannonical_NSTEPS):\n", + " step = steps_canonical[i+1]\n", + " dbi_mixed_can(step=step)\n", + " off_diagonal_norm_history_mixed_can.append(dbi_mixed_can.off_diagonal_norm)\n", + " steps_mixed_can.append(step)\n", + " \n", + "print(\"After 2 steps, off diagonal norm:\", dbi_mixed_can.off_diagonal_norm)\n", + "print(\"By comparison, the Pauli-Z:\", off_diagonal_norm_history[2])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Continue the remaining steps with Pauli-Z operators\n", + "Z_optimal_mixed_can = [\"Cannonical\" for _ in range(cannonical_NSTEPS)]\n", + "remaining_NSTEPS = NSTEPS - cannonical_NSTEPS\n", + "dbi_mixed_can.mode = DoubleBracketGeneratorType.single_commutator\n", + "for _ in range(remaining_NSTEPS):\n", + " dbi_mixed_can, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed_can, Z_ops, scheduling=scheduling, compare_canonical=False, max_evals=max_evals, step_max=step_max)\n", + " off_diagonal_norm_history_mixed_can.append(dbi_mixed_can.off_diagonal_norm)\n", + " steps_mixed_can.append(step)\n", + " if idx == len(Z_ops):\n", + " Z_optimal_mixed.append('Canonical')\n", + " elif flip_sign < 0:\n", + " Z_optimal_mixed.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal_mixed.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{remaining_NSTEPS}: {step} with operator {Z_optimal_mixed_can[-1]}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "# plt.plot(steps, off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "# plt.plot(steps_canonical, off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.plot(off_diagonal_norm_history, label=\"Pauli-Z\")\n", + "plt.plot(off_diagonal_norm_history_canonical, label=\"Canonical\")\n", + "plt.plot(off_diagonal_norm_history_mixed, label=\"Mixed: optimal steps\")\n", + "plt.plot(off_diagonal_norm_history_mixed_can, label=\"Mixed: initial canonical\")\n", + "plt.xlabel(\"Iterations\")\n", + "plt.ylabel(\"Norm off-diagonal restriction\")\n", + "plt.title(\"Compare Variational Pauli-Z with Canonical\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example also shows that the canonical generator is more likely to drive the model into a local minimum than variationally assigned diagonal operator, and that it is hard to get it unstuck even with the Pauli-Z operators." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/dbi/dbi_strategy_magnetic_field.ipynb b/examples/dbi/dbi_strategy_magnetic_field.ipynb index d324efdb38..cb65b76541 100644 --- a/examples/dbi/dbi_strategy_magnetic_field.ipynb +++ b/examples/dbi/dbi_strategy_magnetic_field.ipynb @@ -67,7 +67,10 @@ "# initialize dbi object\n", "nqubits = 5\n", "h0 = random_hermitian(2**nqubits, seed=2)\n", - "dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0))\n", + "scheduling = DoubleBracketScheduling.use_hyperopt\n", + "mode = DoubleBracketGeneratorType.single_commutator\n", + "n_taylor = 5\n", + "dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), scheduling=scheduling, mode=mode)\n", "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)\n", "visualize_matrix(dbi.h.matrix, title=f'Random hamiltonian with L={nqubits}')" ] @@ -82,7 +85,7 @@ "onsite_Z_ops = generate_onsite_Z_ops(nqubits)\n", "d_coef = onsite_Z_decomposition(dbi.h.matrix, onsite_Z_ops)\n", "d = sum([d_coef[i] * onsite_Z_ops[i] for i in range(nqubits)])\n", - "grad, s = gradient_onsite_Z(dbi,d,3, onsite_Z_ops)\n", + "grad, s = gradient_onsite_Z(dbi,d,5, onsite_Z_ops)\n", "print('The initial D coefficients:', d_coef)\n", "print('Gradient:', grad)\n", "print('s:', s)" @@ -95,11 +98,12 @@ "outputs": [], "source": [ "iters = 30\n", - "off_diagonal_norm_tot = [dbi.off_diagonal_norm]\n", + "off_diagonal_norm = [dbi.off_diagonal_norm]\n", "s_step = [0]\n", "for i in range(iters):\n", " s, d_coef, d = gradient_descent_onsite_Z(dbi, d_coef, d, onsite_Z_ops=onsite_Z_ops, max_evals=100)\n", - " off_diagonal_norm_tot.append(dbi.off_diagonal_norm)\n", + " dbi(step=s, d=d)\n", + " off_diagonal_norm.append(dbi.off_diagonal_norm)\n", " s_step.append(s)" ] }, @@ -110,7 +114,7 @@ "outputs": [], "source": [ "plt.title(str(nqubits) + ' spins magnetic field diagonalization')\n", - "plt.plot(off_diagonal_norm_tot)\n", + "plt.plot(off_diagonal_norm)\n", "plt.xlabel('Iteration')\n", "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" ] @@ -151,7 +155,7 @@ "# backend\n", "set_backend(\"qibojit\", \"numba\")\n", "# initialize dbi object\n", - "dbi_TFIM = DoubleBracketIteration(deepcopy(H_TFIM))" + "dbi_TFIM = DoubleBracketIteration(deepcopy(H_TFIM), scheduling=scheduling, mode=mode)" ] }, { @@ -164,7 +168,7 @@ "onsite_Z_ops = generate_onsite_Z_ops(nqubits)\n", "d_coef = onsite_Z_decomposition(dbi_TFIM.h.matrix, onsite_Z_ops)\n", "d = sum([d_coef[i] * onsite_Z_ops[i] for i in range(nqubits)])\n", - "grad, s = gradient_onsite_Z(dbi_TFIM,d,3, onsite_Z_ops)\n", + "grad, s = gradient_onsite_Z(dbi_TFIM,d,5, onsite_Z_ops)\n", "print('Initial off-diagonal norm:', dbi.off_diagonal_norm)\n", "print('The initial D coefficients:', d_coef)\n", "print('Gradient:', grad)\n", @@ -177,13 +181,14 @@ "metadata": {}, "outputs": [], "source": [ - "NSTEPS = 30\n", - "off_diagonal_norm_tot = [dbi_TFIM.off_diagonal_norm]\n", - "s_step = [0]\n", + "NSTEPS = 15\n", + "off_diagonal_norm_delta = [dbi_TFIM.off_diagonal_norm]\n", + "s_step_delta = [0]\n", "for _ in range(NSTEPS):\n", - " s, d_coef, d = gradient_descent_onsite_Z(dbi_TFIM, d_coef, d, onsite_Z_ops=onsite_Z_ops, max_evals=100)\n", - " off_diagonal_norm_tot.append(dbi_TFIM.off_diagonal_norm)\n", - " s_step.append(s)\n", + " s, d_coef, d = gradient_descent_onsite_Z(dbi_TFIM, d_coef, d, onsite_Z_ops=onsite_Z_ops, max_evals=100, n_taylor=5, use_ds=True)\n", + " dbi_TFIM(step=s, d=d)\n", + " off_diagonal_norm_delta.append(dbi_TFIM.off_diagonal_norm)\n", + " s_step_delta.append(s)\n", " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {s} with d_coef {d_coef}, loss {dbi_TFIM.off_diagonal_norm}\")" ] }, @@ -194,7 +199,7 @@ "outputs": [], "source": [ "plt.title(str(nqubits) + ' spins TFIM magnetic field diagonalization')\n", - "plt.plot(off_diagonal_norm_tot)\n", + "plt.plot(off_diagonal_norm_delta)\n", "plt.xlabel('Iteration')\n", "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" ] @@ -216,6 +221,85 @@ "## Different initial `d`\n", "Next, we show the effect of different choices of the initial direction of the gradient descent method." ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "H = H_TFIM.matrix\n", + "L = int(np.log2(H.shape[0]))\n", + "N = np.diag(np.linspace(np.min(np.diag(H)),np.max(np.diag(H)),2**L))\n", + "d_coef = onsite_Z_decomposition(N, onsite_Z_ops)\n", + "print(d_coef)\n", + "d = sum([d_coef[i] * onsite_Z_ops[i] for i in range(nqubits)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "visualize_matrix(H, 'Initial hamiltonian')\n", + "visualize_matrix(N, 'Min-max diagonal matrix')\n", + "visualize_matrix(d, 'Min-max projection onsite-Z')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we see that the min-max diagonal operator can be correctly decomposed into onsite-Z operators. Then we generate the diagonalization curve and compare with other initializations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# backend\n", + "set_backend(\"qibojit\", \"numba\")\n", + "# initialize dbi object\n", + "dbi_TFIM_MMH = DoubleBracketIteration(deepcopy(H_TFIM), scheduling=scheduling, mode=mode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NSTEPS = 15\n", + "off_diagonal_norm_MMH = [dbi_TFIM_MMH.off_diagonal_norm]\n", + "s_step_MMH = [0]\n", + "# d = np.diag(np.linspace(np.min(np.diag(dbi_TFIM_MMH.h.matrix)),np.max(np.diag(dbi_TFIM_MMH.h.matrix)),2**nqubits))\n", + "# d_coef = onsite_Z_decomposition(d, onsite_Z_ops)\n", + "for _ in range(NSTEPS):\n", + " d = np.diag(np.linspace(np.min(np.diag(dbi_TFIM_MMH.h.matrix)),np.max(np.diag(dbi_TFIM_MMH.h.matrix)),2**nqubits))\n", + " d_coef = onsite_Z_decomposition(d, onsite_Z_ops)\n", + " s, d_coef, d = gradient_descent_onsite_Z(dbi_TFIM_MMH, d_coef, d, onsite_Z_ops=onsite_Z_ops, max_evals=100)\n", + " dbi_TFIM_MMH(d=d, step=s)\n", + " off_diagonal_norm_MMH.append(dbi_TFIM_MMH.off_diagonal_norm)\n", + " s_step_MMH.append(s)\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {s} with d_coef {d_coef}, loss {dbi_TFIM_MMH.off_diagonal_norm}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' spins TFIM magnetic field diagonalization')\n", + "plt.plot(off_diagonal_norm_MMH, label='MMH')\n", + "plt.plot(off_diagonal_norm_delta, label='delta')\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')\n", + "plt.legend()" + ] } ], "metadata": { diff --git a/src/qibo/models/dbi/utils.py b/src/qibo/models/dbi/utils.py index b0a89c116a..f42ab174d3 100644 --- a/src/qibo/models/dbi/utils.py +++ b/src/qibo/models/dbi/utils.py @@ -225,7 +225,11 @@ def derivative_product(k1, k2): def gradient_onsite_Z( - dbi_object: DoubleBracketIteration, d: np.array, n_taylor=3, onsite_Z_ops=None + dbi_object: DoubleBracketIteration, + d: np.array, + n_taylor=3, + onsite_Z_ops=None, + use_ds=False, ): """Calculate the gradient of loss function with respect to onsite Pauli-Z coefficients""" # n is the highest order for calculating s @@ -236,12 +240,15 @@ def gradient_onsite_Z( onsite_Z_ops = generate_onsite_Z_ops(nqubits) grad = np.zeros(nqubits) s, coef = dbi_object.polynomial_step( + d=d, n=n_taylor, backup_scheduling=DoubleBracketScheduling.use_polynomial_approximation, ) a, b, c = coef[len(coef) - 3 :] for i in range(nqubits): da, db, dc, ds = ds_di_onsite_Z(dbi_object, d, i, [a, b, c], onsite_Z_ops) + if use_ds is True: + ds = 0 grad[i] = ( s**3 / 3 * da + s**2 / 2 * db @@ -281,7 +288,6 @@ def gradient_descent_onsite_Z( d: np.array = None, n_taylor: int = 3, onsite_Z_ops=None, - grad_tol: float = 1e-2, lr_min: float = 1e-5, lr_max: float = 1, max_evals: int = 100, @@ -289,6 +295,7 @@ def gradient_descent_onsite_Z( optimizer: callable = None, look_ahead: int = 1, verbose: bool = False, + use_ds: bool = True, ): nqubits = int(np.log2(dbi_object.h.matrix.shape[0])) if onsite_Z_ops is None: @@ -296,7 +303,7 @@ def gradient_descent_onsite_Z( if d is None: d = sum([d_coef[i] * onsite_Z_ops[i] for i in range(nqubits)]) grad, s = gradient_onsite_Z( - dbi_object, d, n_taylor=n_taylor, onsite_Z_ops=onsite_Z_ops + dbi_object, d, n_taylor=n_taylor, onsite_Z_ops=onsite_Z_ops, use_ds=use_ds ) # optimize gradient descent step with hyperopt if space is None: @@ -320,5 +327,11 @@ def func_loss_to_lr(lr): d_coef = [d_coef[j] - grad[j] * lr for j in range(nqubits)] d = sum([d_coef[i] * onsite_Z_ops[i] for i in range(nqubits)]) - dbi_object(step=s, d=d) return s, d_coef, d + + +def diagonal_min_max(matrix: np.array): + L = int(np.log2(matrix.shape[0])) + D = np.linspace(np.min(np.diag(matrix)), np.max(np.diag(matrix)), 2**L) + D = np.diag(D) + return D diff --git a/tests/test_models_dbi_utils 2.py b/tests/test_models_dbi_utils 2.py new file mode 100644 index 0000000000..cd9f74e9de --- /dev/null +++ b/tests/test_models_dbi_utils 2.py @@ -0,0 +1,50 @@ +""""Testing utils for DoubleBracketIteration model""" + +import numpy as np +import pytest + +from qibo import set_backend +from qibo.hamiltonians import Hamiltonian +from qibo.models.dbi.double_bracket import ( + DoubleBracketGeneratorType, + DoubleBracketIteration, +) +from qibo.models.dbi.utils import * +from qibo.quantum_info import random_hermitian + +NSTEPS = 5 +"""Number of steps for evolution.""" + + +@pytest.mark.parametrize("nqubits", [3, 4, 5]) +def test_generate_Z_operators(backend, nqubits): + h0 = random_hermitian(2**nqubits) + dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0)) + generate_Z = generate_Z_operators(nqubits) + Z_ops = list(generate_Z.values()) + + delta_h0 = dbi.diagonal_h_matrix + dephasing_channel = (sum([Z_op @ h0 @ Z_op for Z_op in Z_ops]) + h0) / 2**nqubits + norm_diff = np.linalg.norm(delta_h0 - dephasing_channel) + + assert norm_diff < 1e-3 + + +@pytest.mark.parametrize("nqubits", [3, 4, 5]) +@pytest.mark.parametrize("step", [0.1, None]) +def test_select_best_dbr_generator(backend, nqubits, step): + h0 = random_hermitian(2**nqubits, seed=1, backend=backend) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + ) + generate_Z = generate_Z_operators(nqubits) + Z_ops = list(generate_Z.values()) + initial_off_diagonal_norm = dbi.off_diagonal_norm + + for _ in range(NSTEPS): + dbi, idx, step_optimize, flip = select_best_dbr_generator( + dbi, Z_ops, step=step, compare_canonical=True + ) + + assert initial_off_diagonal_norm > dbi.off_diagonal_norm