From 2bd6331eca5755d3348f8deba06f74e530f218fb Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Thu, 21 Nov 2024 15:41:16 +0100 Subject: [PATCH] Implementation of the dynamic MPF feature (#33) * Implementation of the dynamic MPF feature Co-authored-by: Kate Marshall Co-authored-by: Max Rossmannek Co-authored-by: Niall Robertson Co-authored-by: Alberto Baiardi * Fix issue with empty attributes section * docs: some minor fixes based on offline review Co-authored-by: Kate Marshall * docs: update reference after its peer-review publication * docs: update module overview in landing pages * docs: attempt to fix alert rendering * Revert "docs: attempt to fix alert rendering" This reverts commit 32091cbc7b92a7c70f2be5e7d3a62be91f2bab4e. * fix: close div in README --------- Co-authored-by: Kate Marshall Co-authored-by: Kate Marshall Co-authored-by: Niall Robertson Co-authored-by: Alberto Baiardi Co-authored-by: Eric Arellano <14852634+Eric-Arellano@users.noreply.github.com> --- README.md | 22 +- .../class_without_inherited_members.rst | 34 ++ docs/apidocs/index.rst | 3 + ...iskit_addon_mpf.backends.quimb_circuit.rst | 8 + ...qiskit_addon_mpf.backends.quimb_layers.rst | 8 + .../qiskit_addon_mpf.backends.quimb_tebd.rst | 8 + docs/apidocs/qiskit_addon_mpf.backends.rst | 8 + ...qiskit_addon_mpf.backends.tenpy_layers.rst | 8 + .../qiskit_addon_mpf.backends.tenpy_tebd.rst | 8 + docs/apidocs/qiskit_addon_mpf.costs.rst | 8 + docs/apidocs/qiskit_addon_mpf.dynamic.rst | 8 + docs/apidocs/qiskit_addon_mpf.static.rst | 18 - docs/conf.py | 2 + docs/explanations/mpf_stability.ipynb | 2 +- docs/how_tos/choose_trotter_steps.ipynb | 10 +- docs/how_tos/using_approximate_model.ipynb | 18 +- docs/index.rst | 25 +- docs/install.rst | 2 +- docs/tutorials/01_getting_started.ipynb | 26 +- pyproject.toml | 18 +- qiskit_addon_mpf/backends/__init__.py | 64 +++ qiskit_addon_mpf/backends/interface.py | 117 ++++++ .../backends/quimb_circuit/__init__.py | 110 ++++++ .../backends/quimb_circuit/evolver.py | 92 +++++ .../backends/quimb_circuit/state.py | 78 ++++ .../backends/quimb_layers/__init__.py | 230 +++++++++++ .../backends/quimb_layers/evolver.py | 75 ++++ .../backends/quimb_layers/model.py | 150 +++++++ .../backends/quimb_tebd/__init__.py | 139 +++++++ .../backends/quimb_tebd/evolver.py | 160 ++++++++ qiskit_addon_mpf/backends/quimb_tebd/state.py | 234 +++++++++++ .../backends/tenpy_layers/__init__.py | 245 ++++++++++++ .../backends/tenpy_layers/evolver.py | 171 ++++++++ .../backends/tenpy_layers/model.py | 162 ++++++++ .../backends/tenpy_tebd/__init__.py | 163 ++++++++ .../backends/tenpy_tebd/evolver.py | 152 +++++++ qiskit_addon_mpf/backends/tenpy_tebd/state.py | 136 +++++++ qiskit_addon_mpf/costs/__init__.py | 44 +++ qiskit_addon_mpf/{static => costs}/exact.py | 21 +- qiskit_addon_mpf/costs/frobenius.py | 126 ++++++ qiskit_addon_mpf/costs/lse.py | 73 ++++ .../sum_of_squares.py} | 23 +- qiskit_addon_mpf/dynamic.py | 372 ++++++++++++++++++ qiskit_addon_mpf/{static/lse.py => static.py} | 73 ++-- .../notes/dynamic-mpf-cf18e0b4d2bdeaff.yaml | 19 + test/{static => backends}/__init__.py | 0 .../backends/quimb_circuit}/__init__.py | 14 - test/backends/quimb_circuit/test_e2e.py | 86 ++++ test/backends/quimb_circuit/test_state.py | 36 ++ test/backends/quimb_layers/__init__.py | 11 + test/backends/quimb_layers/test_e2e.py | 183 +++++++++ test/backends/quimb_layers/test_model.py | 89 +++++ test/backends/quimb_tebd/__init__.py | 11 + test/backends/quimb_tebd/test_e2e.py | 124 ++++++ test/backends/quimb_tebd/test_state.py | 29 ++ test/backends/tenpy_layers/__init__.py | 11 + test/backends/tenpy_layers/test_e2e.py | 181 +++++++++ test/backends/tenpy_layers/test_model.py | 90 +++++ test/backends/tenpy_tebd/__init__.py | 11 + test/backends/tenpy_tebd/test_e2e.py | 120 ++++++ test/backends/tenpy_tebd/test_state.py | 32 ++ test/costs/__init__.py | 11 + test/costs/test_exact.py | 135 +++++++ test/costs/test_lse.py | 38 ++ test/costs/test_sum_of_squares.py | 75 ++++ test/static/test_approximate.py | 79 ---- test/static/test_exact.py | 137 ------- test/{static/test_lse.py => test_static.py} | 31 +- 68 files changed, 4636 insertions(+), 371 deletions(-) create mode 100644 docs/_templates/autosummary/class_without_inherited_members.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.quimb_circuit.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.quimb_layers.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.quimb_tebd.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.tenpy_layers.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.backends.tenpy_tebd.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.costs.rst create mode 100644 docs/apidocs/qiskit_addon_mpf.dynamic.rst create mode 100644 qiskit_addon_mpf/backends/__init__.py create mode 100644 qiskit_addon_mpf/backends/interface.py create mode 100644 qiskit_addon_mpf/backends/quimb_circuit/__init__.py create mode 100644 qiskit_addon_mpf/backends/quimb_circuit/evolver.py create mode 100644 qiskit_addon_mpf/backends/quimb_circuit/state.py create mode 100644 qiskit_addon_mpf/backends/quimb_layers/__init__.py create mode 100644 qiskit_addon_mpf/backends/quimb_layers/evolver.py create mode 100644 qiskit_addon_mpf/backends/quimb_layers/model.py create mode 100644 qiskit_addon_mpf/backends/quimb_tebd/__init__.py create mode 100644 qiskit_addon_mpf/backends/quimb_tebd/evolver.py create mode 100644 qiskit_addon_mpf/backends/quimb_tebd/state.py create mode 100644 qiskit_addon_mpf/backends/tenpy_layers/__init__.py create mode 100644 qiskit_addon_mpf/backends/tenpy_layers/evolver.py create mode 100644 qiskit_addon_mpf/backends/tenpy_layers/model.py create mode 100644 qiskit_addon_mpf/backends/tenpy_tebd/__init__.py create mode 100644 qiskit_addon_mpf/backends/tenpy_tebd/evolver.py create mode 100644 qiskit_addon_mpf/backends/tenpy_tebd/state.py create mode 100644 qiskit_addon_mpf/costs/__init__.py rename qiskit_addon_mpf/{static => costs}/exact.py (74%) create mode 100644 qiskit_addon_mpf/costs/frobenius.py create mode 100644 qiskit_addon_mpf/costs/lse.py rename qiskit_addon_mpf/{static/approximate.py => costs/sum_of_squares.py} (75%) create mode 100644 qiskit_addon_mpf/dynamic.py rename qiskit_addon_mpf/{static/lse.py => static.py} (63%) create mode 100644 releasenotes/notes/dynamic-mpf-cf18e0b4d2bdeaff.yaml rename test/{static => backends}/__init__.py (100%) rename {qiskit_addon_mpf/static => test/backends/quimb_circuit}/__init__.py (60%) create mode 100644 test/backends/quimb_circuit/test_e2e.py create mode 100644 test/backends/quimb_circuit/test_state.py create mode 100644 test/backends/quimb_layers/__init__.py create mode 100644 test/backends/quimb_layers/test_e2e.py create mode 100644 test/backends/quimb_layers/test_model.py create mode 100644 test/backends/quimb_tebd/__init__.py create mode 100644 test/backends/quimb_tebd/test_e2e.py create mode 100644 test/backends/quimb_tebd/test_state.py create mode 100644 test/backends/tenpy_layers/__init__.py create mode 100644 test/backends/tenpy_layers/test_e2e.py create mode 100644 test/backends/tenpy_layers/test_model.py create mode 100644 test/backends/tenpy_tebd/__init__.py create mode 100644 test/backends/tenpy_tebd/test_e2e.py create mode 100644 test/backends/tenpy_tebd/test_state.py create mode 100644 test/costs/__init__.py create mode 100644 test/costs/test_exact.py create mode 100644 test/costs/test_lse.py create mode 100644 test/costs/test_sum_of_squares.py delete mode 100644 test/static/test_approximate.py delete mode 100644 test/static/test_exact.py rename test/{static/test_lse.py => test_static.py} (73%) diff --git a/README.md b/README.md index fc0a82f..886a6b6 100644 --- a/README.md +++ b/README.md @@ -12,6 +12,7 @@ [![Downloads](https://img.shields.io/pypi/dm/qiskit-addon-mpf.svg?label=Downloads)](https://pypi.org/project/qiskit-addon-mpf/) [![Tests](https://github.com/Qiskit/qiskit-addon-mpf/actions/workflows/test_latest_versions.yml/badge.svg)](https://github.com/Qiskit/qiskit-addon-mpf/actions/workflows/test_latest_versions.yml) [![Coverage](https://coveralls.io/repos/github/Qiskit/qiskit-addon-mpf/badge.svg?branch=main)](https://coveralls.io/github/Qiskit/qiskit-addon-mpf?branch=main) + # Qiskit addon: multi-product formulas (MPF) @@ -34,9 +35,10 @@ This package contains the Qiskit addon for multi-product formulas (MPFs). These can be used to reduce the Trotter error of Hamiltonian dynamics. -This package currently contains the following submodules: +This package currently contains the following main entry points for users: - `qiskit_addon_mpf.static` for working with static MPFs [1-2](#references) +- `qiskit_addon_mpf.dynamic` for working with dynamic MPFs [2-3](#references) ---------------------------------------------------------------------------------------------------- @@ -56,6 +58,21 @@ pip install 'qiskit-addon-mpf' For more installation information refer to these [installation instructions](docs/install.rst). +#### Optional dependencies + +The `qiskit-addon-mpf` package has a number of optional dependencies which enable certain features. +The dynamic MPF feature (see [2-3](#references)) is one such example. +You can install the related optional dependencies like so: + +```bash +pip install 'qiskit-addon-mpf[dynamic]' +``` + +> [!IMPORTANT] +> The optional dependency [TeNPy](https://github.com/tenpy/tenpy) was previously offered under a GPLv3 license. +> As of the release of [v1.0.4](https://github.com/tenpy/tenpy/releases/tag/v1.0.4) on October 2nd, 2024, it has been offered under the Apache v2 license. +> The license of this package is only compatible with Apache-licensed versions of TeNPy. + ---------------------------------------------------------------------------------------------------- ### Deprecation Policy @@ -85,7 +102,8 @@ We use [GitHub issues](https://github.com/Qiskit/qiskit-addon-mpf/issues/new/cho ### References 1. A. Carrera Vazquez, D. J. Egger, D. Ochsner, and S. Wörner, [Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation](https://quantum-journal.org/papers/q-2023-07-25-1067/), Quantum 7, 1067 (2023). -2. S. Zhuk, N. Robertson, and S. Bravyi, [Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation](https://arxiv.org/abs/2306.12569v2), arXiv:2306.12569 [quant-ph]. +2. S. Zhuk, N. Robertson, and S. Bravyi, [Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309), Phys. Rev. Research 6, 033309 (2024). +3. N. Robertson, et al. [Tensor Network enhanced Dynamic Multiproduct Formulas](https://arxiv.org/abs/2407.17405v2), arXiv:2407.17405v2 [quant-ph]. ---------------------------------------------------------------------------------------------------- diff --git a/docs/_templates/autosummary/class_without_inherited_members.rst b/docs/_templates/autosummary/class_without_inherited_members.rst new file mode 100644 index 0000000..66fe739 --- /dev/null +++ b/docs/_templates/autosummary/class_without_inherited_members.rst @@ -0,0 +1,34 @@ +{# + We show all the class's methods and attributes on the same page. By default, we document + all methods, including those defined by parent classes. +-#} + +{{ objname | escape | underline }} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :no-members: + :no-inherited-members: + :no-special-members: + :show-inheritance: + +{% block attributes_summary %} + {% set wanted_attributes = (attributes | reject('in', inherited_members) | list) %} + {% if wanted_attributes %} + .. rubric:: Attributes + {% for item in wanted_attributes %} + .. autoattribute:: {{ item }} + {%- endfor %} + {% endif %} +{% endblock -%} + +{% block methods_summary %} + {% set wanted_methods = (methods | reject('==', '__init__') | reject('in', inherited_members) | list) %} + {% if wanted_methods %} + .. rubric:: Methods + {% for item in wanted_methods %} + .. automethod:: {{ item }} + {%- endfor %} + {% endif %} +{% endblock %} diff --git a/docs/apidocs/index.rst b/docs/apidocs/index.rst index 7985618..8c4e52b 100644 --- a/docs/apidocs/index.rst +++ b/docs/apidocs/index.rst @@ -6,3 +6,6 @@ :maxdepth: 1 qiskit_addon_mpf.static + qiskit_addon_mpf.dynamic + qiskit_addon_mpf.costs + qiskit_addon_mpf.backends diff --git a/docs/apidocs/qiskit_addon_mpf.backends.quimb_circuit.rst b/docs/apidocs/qiskit_addon_mpf.backends.quimb_circuit.rst new file mode 100644 index 0000000..595a2d4 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.quimb_circuit.rst @@ -0,0 +1,8 @@ +============================================================================ +Quimb circuit-based backend (:mod:`qiskit_addon_mpf.backends.quimb_circuit`) +============================================================================ + +.. automodule:: qiskit_addon_mpf.backends.quimb_circuit + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.backends.quimb_layers.rst b/docs/apidocs/qiskit_addon_mpf.backends.quimb_layers.rst new file mode 100644 index 0000000..8af3a4e --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.quimb_layers.rst @@ -0,0 +1,8 @@ +========================================================================= +Quimb layer-based backend (:mod:`qiskit_addon_mpf.backends.quimb_layers`) +========================================================================= + +.. automodule:: qiskit_addon_mpf.backends.quimb_layers + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.backends.quimb_tebd.rst b/docs/apidocs/qiskit_addon_mpf.backends.quimb_tebd.rst new file mode 100644 index 0000000..cdd2e60 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.quimb_tebd.rst @@ -0,0 +1,8 @@ +================================================================ +Quimb TEBD backend (:mod:`qiskit_addon_mpf.backends.quimb_tebd`) +================================================================ + +.. automodule:: qiskit_addon_mpf.backends.quimb_tebd + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.backends.rst b/docs/apidocs/qiskit_addon_mpf.backends.rst new file mode 100644 index 0000000..beb3b59 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.rst @@ -0,0 +1,8 @@ +=========================================== +Backends (:mod:`qiskit_addon_mpf.backends`) +=========================================== + +.. automodule:: qiskit_addon_mpf.backends + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.backends.tenpy_layers.rst b/docs/apidocs/qiskit_addon_mpf.backends.tenpy_layers.rst new file mode 100644 index 0000000..b05ae5a --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.tenpy_layers.rst @@ -0,0 +1,8 @@ +========================================================================= +TeNPy layer-based backend (:mod:`qiskit_addon_mpf.backends.tenpy_layers`) +========================================================================= + +.. automodule:: qiskit_addon_mpf.backends.tenpy_layers + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.backends.tenpy_tebd.rst b/docs/apidocs/qiskit_addon_mpf.backends.tenpy_tebd.rst new file mode 100644 index 0000000..e57fdb5 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.backends.tenpy_tebd.rst @@ -0,0 +1,8 @@ +================================================================ +TeNPy TEBD backend (:mod:`qiskit_addon_mpf.backends.tenpy_tebd`) +================================================================ + +.. automodule:: qiskit_addon_mpf.backends.tenpy_tebd + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.costs.rst b/docs/apidocs/qiskit_addon_mpf.costs.rst new file mode 100644 index 0000000..7ad82f9 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.costs.rst @@ -0,0 +1,8 @@ +============================================== +Cost Functions (:mod:`qiskit_addon_mpf.costs`) +============================================== + +.. automodule:: qiskit_addon_mpf.costs + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.dynamic.rst b/docs/apidocs/qiskit_addon_mpf.dynamic.rst new file mode 100644 index 0000000..5880f67 --- /dev/null +++ b/docs/apidocs/qiskit_addon_mpf.dynamic.rst @@ -0,0 +1,8 @@ +============================================== +Dynamic MPFs (:mod:`qiskit_addon_mpf.dynamic`) +============================================== + +.. automodule:: qiskit_addon_mpf.dynamic + :no-members: + :no-inherited-members: + :no-special-members: diff --git a/docs/apidocs/qiskit_addon_mpf.static.rst b/docs/apidocs/qiskit_addon_mpf.static.rst index 70cf5ff..7a8264d 100644 --- a/docs/apidocs/qiskit_addon_mpf.static.rst +++ b/docs/apidocs/qiskit_addon_mpf.static.rst @@ -6,21 +6,3 @@ Static MPFs (:mod:`qiskit_addon_mpf.static`) :no-members: :no-inherited-members: :no-special-members: - -.. currentmodule:: qiskit_addon_mpf.static - -Linear system of equations utilities ------------------------------------- - -.. autoclass:: LSE -.. autofunction:: setup_lse - -Exact static MPF coefficients ------------------------------ - -.. autofunction:: setup_exact_model - -Approximate static MPF coefficients ------------------------------------ - -.. autofunction:: setup_approximate_model diff --git a/docs/conf.py b/docs/conf.py index 858e344..cdd40eb 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -107,6 +107,8 @@ "qiskit": ("https://docs.quantum.ibm.com/api/qiskit/", None), "rustworkx": ("https://www.rustworkx.org/", None), "qiskit_addon_utils": ("https://qiskit.github.io/qiskit-addon-utils/", None), + "quimb": ("https://quimb.readthedocs.io/en/latest/", None), + "tenpy": ("https://tenpy.readthedocs.io/en/latest/", None), } plot_working_directory = "." diff --git a/docs/explanations/mpf_stability.ipynb b/docs/explanations/mpf_stability.ipynb index a2010a0..bb9d266 100644 --- a/docs/explanations/mpf_stability.ipynb +++ b/docs/explanations/mpf_stability.ipynb @@ -17,7 +17,7 @@ "Over in the guide [How to choose the Trotter steps for an MPF](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) we introduced a heuristic that limits the **smallest** $k_j$ value ($k_{\\text{min}}$) based on the total evolution time, $t$, that we would like to reach.\n", "In particular, we state that $t/k_{\\text{min}} \\lt 1$ must be satisfied (noting that $t/k_{\\text{min}} \\leq 1$ empirically seems to work fine, too).\n", "\n", - "On this page, we will analyze the behavior of MPFs is not guaranteed to work well, when this constraint gets violated." + "On this page, we will analyze when the behavior of MPFs is not guaranteed to work well, such as when this constraint gets violated." ] }, { diff --git a/docs/how_tos/choose_trotter_steps.ipynb b/docs/how_tos/choose_trotter_steps.ipynb index 47eefbf..51218fb 100644 --- a/docs/how_tos/choose_trotter_steps.ipynb +++ b/docs/how_tos/choose_trotter_steps.ipynb @@ -76,7 +76,7 @@ "\n", "where $\\chi$ is the order of our individual product formulas.\n", "\n", - "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569" + "[Zhuk et al., 2023]: https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309" ] }, { @@ -93,7 +93,7 @@ "source": [ "### Preparing our LSE\n", "\n", - "In a first step, we prepare the LSE, $Ax=b$, leveraging a feature of the [setup_lse](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_lse.html) function that allows us to use a [cvxpy.Parameter](https://www.cvxpy.org/api_reference/cvxpy.expressions.html#parameter) object for our $k_j$ values." + "In a first step, we prepare the LSE, $Ax=b$, leveraging a feature of the [setup_static_lse](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_static_lse.html) function that allows us to use a [cvxpy.Parameter](https://www.cvxpy.org/api_reference/cvxpy.expressions.html#parameter) object for our $k_j$ values." ] }, { @@ -133,10 +133,10 @@ "outputs": [], "source": [ "import cvxpy as cp\n", - "from qiskit_addon_mpf.static import setup_lse\n", + "from qiskit_addon_mpf.static import setup_static_lse\n", "\n", "ks = cp.Parameter(3, integer=True)\n", - "lse = setup_lse(ks, order=order, symmetric=symmetric)" + "lse = setup_static_lse(ks, order=order, symmetric=symmetric)" ] }, { @@ -208,7 +208,7 @@ "source": [ "Finally, we sort our $k_j$ by the re-weighted L1-norm presented by [Zhuk et al., 2023]:\n", "\n", - "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569" + "[Zhuk et al., 2023]: https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309" ] }, { diff --git a/docs/how_tos/using_approximate_model.ipynb b/docs/how_tos/using_approximate_model.ipynb index 776e486..c85e4d3 100644 --- a/docs/how_tos/using_approximate_model.ipynb +++ b/docs/how_tos/using_approximate_model.ipynb @@ -14,7 +14,7 @@ "source": [ "# How to use the approximate model\n", "\n", - "You have already seen a first example of the [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) function in the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial.\n", + "You have already seen a first example of the [setup_sum_of_squares_problem](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.setup_sum_of_squares_problem.html) function in the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial.\n", "In this guide, we are going to revisit that function in more detail." ] }, @@ -92,7 +92,7 @@ "source": [ "## Choosing $k_j$\n", "\n", - "In this guide, we are going to focus on the tweaks that we can apply to the [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) and its solver that are returned by the [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html). Therefore, we are not concerned with choosing any specific $k_j$ values and can simply reuse the ones from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial." + "In this guide, we are going to focus on the tweaks that we can apply to the [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) and its solver that are returned by the [setup_sum_of_squares_problem](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.setup_sum_of_squares_problem.html). Therefore, we are not concerned with choosing any specific $k_j$ values and can simply reuse the ones from the [Getting started](https://qiskit.github.io/qiskit-addon-mpf/tutorials/01_getting_started.html) tutorial." ] }, { @@ -135,9 +135,9 @@ }, "outputs": [], "source": [ - "from qiskit_addon_mpf.static import setup_lse\n", + "from qiskit_addon_mpf.static import setup_static_lse\n", "\n", - "lse = setup_lse(trotter_steps, order=2, symmetric=True)" + "lse = setup_static_lse(trotter_steps, order=2, symmetric=True)" ] }, { @@ -222,14 +222,14 @@ "\n", "Now, we will construct the approximate model.\n", "The idea of using approximate solutions $\\tilde{x}$ in order to ensure a smaller L1-norm was proposed by [Zhuk et al., 2023].\n", - "The model is already explained in detail in the API documentation of [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) so we refrain from repeating it here.\n", + "The model is already explained in detail in the API documentation of [setup_sum_of_squares_problem](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.setup_sum_of_squares_problem.html) so we refrain from repeating it here.\n", "Suffice to say, that the model consists of two constraints which ensure that the coefficients, $\\tilde{x}$, sum to $1$ and it enforces their L1-norm to be smaller than the user-provided `max_l1_norm` value.\n", "The optimization objective is to minimize the deviation of $A\\tilde{x}=b$.\n", "\n", "We can see all of this in the output below.\n", - "Here, we are setting `max_l1_norm=3` because we already know that the L1-norm of the exact solution is approximately $3.4$ and we want to find a $\\tilde{x}$ that has a smaller L1-norm (since otherwise this exercise would be rather pointless).\n", + "Here, we are setting `max_l1_norm=3` because we already know that the L1-norm of the exact solution is approximately $3.4$ and we want to find a $\\tilde{x}$ that has a smaller L1-norm.\n", "\n", - "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569" + "[Zhuk et al., 2023]: https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309" ] }, { @@ -255,9 +255,9 @@ } ], "source": [ - "from qiskit_addon_mpf.static import setup_approximate_model\n", + "from qiskit_addon_mpf.costs import setup_sum_of_squares_problem\n", "\n", - "model, coeff_var = setup_approximate_model(lse, max_l1_norm=3)\n", + "model, coeff_var = setup_sum_of_squares_problem(lse, max_l1_norm=3)\n", "print(model)" ] }, diff --git a/docs/index.rst b/docs/index.rst index 7cfc54d..f95573c 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,9 +7,10 @@ Qiskit addon: multi-product formulas (MPF) This package contains the Qiskit addon for multi-product formulas (MPFs). These can be used to reduce the Trotter error of Hamiltonian dynamics. -This package currently contains the following submodules: +This package currently contains the following main entry points for users: - ``qiskit_addon_mpf.static`` for working with static MPFs [`1-2 <#references>`_]. +- ``qiskit_addon_mpf.dynamic`` for working with dynamic MPFs [`2-3 <#references>`_]. Documentation ------------- @@ -28,6 +29,25 @@ We encourage installing this package via ``pip``, when possible: For more installation information refer to the `installation instructions `_ in the documentation. +Optional dependencies ++++++++++++++++++++++ + +The ``qiskit-addon-mpf`` package has a number of optional dependencies which enable certain features. +The dynamic MPF feature (see [`2-3 <#references>`_]) is one such example. +You can install the related optional dependencies like so: + +.. code-block:: bash + + pip install 'qiskit-addon-mpf[dynamic]' + +.. caution:: + The optional dependency `TeNPy `_ was previously offered under a + GPLv3 license. + As of the release of `v1.0.4 `_ on October + 2nd, 2024, it has been offered under the Apache v2 license. + The license of this package is only compatible with Apache-licensed versions of TeNPy. + + Deprecation Policy ------------------ @@ -54,7 +74,8 @@ References ---------- 1. A. Carrera Vazquez, D. J. Egger, D. Ochsner, and S. Wörner, `Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation `_, Quantum 7, 1067 (2023). -2. S. Zhuk, N. Robertson, and S. Bravyi, `Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation `_, arXiv:2306.12569 [quant-ph]. +2. S. Zhuk, N. Robertson, and S. Bravyi, `Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation `_, Phys. Rev. Research 6, 033309 (2024). +3. N. Robertson, et al. `Tensor Network enhanced Dynamic Multiproduct Formulas `_, arXiv:2407.17405v2 [quant-ph]. License ------- diff --git a/docs/install.rst b/docs/install.rst index 6239e69..1096e71 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -56,7 +56,7 @@ If so, the first step is to clone the ``qiskit-addon-mpf`` repository. git clone git@github.com:Qiskit/qiskit-addon-mpf.git -Next, upgrade pip and enter the repository. +Next, upgrade ``pip`` and enter the repository. .. code:: sh diff --git a/docs/tutorials/01_getting_started.ipynb b/docs/tutorials/01_getting_started.ipynb index ae2eb69..cb42503 100644 --- a/docs/tutorials/01_getting_started.ipynb +++ b/docs/tutorials/01_getting_started.ipynb @@ -241,7 +241,7 @@ "\\mu(t) = \\sum_{j} x_j \\rho^{k_j}_{j}\\left(\\frac{t}{k_j}\\right) + \\text{some remaining Trotter error} \\, ,\n", "$$\n", "\n", - "where $x_j$ are our weighting coefficients, $\\rho^{k_j}_j$ is the density matrix corresponding to the pure state obtained by evolving the initial state with the product formula, $S^{k_j}$, involving $k_j$ Trotter steps, and $j$ indexes the number of PFs that make up the MPF.\n", + "where $x_j$ are our weighting coefficients, $\\rho^{k_j}_j$ is the density matrix corresponding to the pure state obtained by evolving the initial state with the product formula, $S^{k_j}$, involving $k_j$ Trotter steps, and $j$ indexes the number of product formulas that make up the MPF.\n", "\n", "The key here is that the remaining Trotter error is smaller than the Trotter error that one would obtain by simply using the largest $k_j$ value!\n", "\n", @@ -269,7 +269,7 @@ "\n", "Determining the static MPF coefficients for a given set of $k_j$ values amounts to solving a linear system of equations:\n", "$Ax=b$, where $x$ are our coefficients of interest, $A$ is a matrix depending on $k_j$ and the type of PF we use ($S$), and $b$ is a vector of constraints.\n", - "For brevity, we are not going to go into more detail here and instead refer you to the documentation of [LSE](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.LSE.html).\n", + "For brevity, we are not going to go into more detail here and instead refer you to the documentation of [LSE](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.LSE.html).\n", "\n", "We can find a solution for $x$ analytically as $x = A^{-1}b$, see e.g. [Carrera Vazquez et al., 2023] or [Zhuk et al., 2023].\n", "However, this exact solution can be _\"ill-conditioned\"_ resulting in very large L1-norms of our coefficients, $x$, which can lead to bad performance of the MPF.\n", @@ -278,7 +278,7 @@ "In the following, you will learn how to do all of this.\n", "\n", "[Carrera Vazquez et al., 2023]: https://quantum-journal.org/papers/q-2023-07-25-1067/\n", - "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569" + "[Zhuk et al., 2023]: https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309" ] }, { @@ -343,7 +343,7 @@ "Here, we are going to use a second-order Suzuki-Trotter formula yielding `order=2` and we will set `symmetric=True`.\n", "\n", "[Carrera Vazquez et al., 2023]: https://quantum-journal.org/papers/q-2023-07-25-1067/\n", - "[Zhuk et al., 2023]: https://arxiv.org/abs/2306.12569" + "[Zhuk et al., 2023]: https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309" ] }, { @@ -369,9 +369,9 @@ } ], "source": [ - "from qiskit_addon_mpf.static import setup_lse\n", + "from qiskit_addon_mpf.static import setup_static_lse\n", "\n", - "lse = setup_lse(trotter_steps, order=2, symmetric=True)\n", + "lse = setup_static_lse(trotter_steps, order=2, symmetric=True)\n", "print(lse)" ] }, @@ -433,7 +433,7 @@ "source": [ "#### Optimizing for $x$ using an exact model\n", "\n", - "Alternatively to computing $x=A^{-1}b$, you can also use [setup_exact_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) to construct a [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance which uses the LSE as constraints and whose optimal solution will yield $x$.\n", + "Alternatively to computing $x=A^{-1}b$, you can also use [setup_exact_problem](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.setup_exact_problem.html) to construct a [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance which uses the LSE as constraints and whose optimal solution will yield $x$.\n", "\n", "In the next section, it will be clear why this interface exists." ] @@ -459,9 +459,9 @@ } ], "source": [ - "from qiskit_addon_mpf.static import setup_exact_model\n", + "from qiskit_addon_mpf.costs import setup_exact_problem\n", "\n", - "model_exact, coeffs_exact = setup_exact_model(lse)\n", + "model_exact, coeffs_exact = setup_exact_problem(lse)\n", "model_exact.solve()\n", "print(coeffs_exact.value)" ] @@ -524,7 +524,7 @@ "It may happen that the L1 norm for the chosen set of $k_j$ values is deemed too high.\n", "If that is the case and you cannot choose a different set of $k_j$ values, you can use an approximate solution to the LSE instead of an exact one.\n", "\n", - "To do so, simply use [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) to construct a different [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance which constrains the L1-norm to a chosen threshold while minimizing the difference of $Ax$ and $b$." + "To do so, simply use [setup_sum_of_squares_problem](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.costs.setup_sum_of_squares_problem.html) to construct a different [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance which constrains the L1-norm to a chosen threshold while minimizing the difference of $Ax$ and $b$." ] }, { @@ -550,9 +550,9 @@ } ], "source": [ - "from qiskit_addon_mpf.static import setup_approximate_model\n", + "from qiskit_addon_mpf.costs import setup_sum_of_squares_problem\n", "\n", - "model_approx, coeffs_approx = setup_approximate_model(lse, max_l1_norm=3.0)\n", + "model_approx, coeffs_approx = setup_sum_of_squares_problem(lse, max_l1_norm=3.0)\n", "model_approx.solve()\n", "print(coeffs_approx.value)\n", "print(np.linalg.norm(coeffs_approx.value, ord=1))" @@ -589,7 +589,7 @@ "### 1c: Setting up the Trotter circuits\n", "\n", "At this point, we have found our expansion coefficients, $x$, and all that is left to do is to generate the Trotterized quantum circuits.\n", - "Once again, the [qiskit_addon_utils.problem_generators](https://qiskit.github.io/qiskit-addon-utils/stubs/qiskit_addon_utils.problem_generators.html) module comes to the rescue with a handy function do just that:" + "Once again, the [qiskit_addon_utils.problem_generators](https://qiskit.github.io/qiskit-addon-utils/stubs/qiskit_addon_utils.problem_generators.html) module comes to the rescue to do just that:" ] }, { diff --git a/pyproject.toml b/pyproject.toml index 4d6c7cc..b5eb56f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,19 +27,31 @@ requires-python = ">=3.9" dependencies = [ "numpy>=1.23", - "cvxpy>=1.5", + "cvxpy>=1.6", "qiskit>=1.2, <2", ] [project.optional-dependencies] +quimb = [ + "quimb>=1.8.3, <2", + "qiskit-quimb", +] +tenpy = [ + "physics-tenpy>=1.0.4, <2", +] +dynamic = [ + "qiskit-addon-mpf[quimb,tenpy]", +] dev = [ "qiskit-addon-mpf[test,doctest,nbtest,lint,docs]", "tox>=4.4.3", ] basetest = [ - "ddt>=1.7", "pytest>=8.0", "pytest-cov>=5.0", + "pytest-subtests>=0.13", + "qiskit-addon-mpf[dynamic]", + "qiskit-addon-utils", ] test = [ "qiskit-addon-mpf[basetest]", @@ -67,7 +79,6 @@ lint = [ notebook-dependencies = [ "qiskit-addon-mpf", "rustworkx[graphviz]>=0.15", - "qiskit-addon-utils", "matplotlib", "ipywidgets", "pylatexenc", @@ -114,7 +125,6 @@ py-version = "3.9" disable = ["all"] enable = [ "reimported", - "no-self-use", "no-else-raise", "redefined-argument-from-local", "redefined-builtin", diff --git a/qiskit_addon_mpf/backends/__init__.py b/qiskit_addon_mpf/backends/__init__.py new file mode 100644 index 0000000..b0bb92f --- /dev/null +++ b/qiskit_addon_mpf/backends/__init__.py @@ -0,0 +1,64 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Optional backends for the :class:`.DynamicMPF` algorithm. + +.. currentmodule:: qiskit_addon_mpf.backends + +Availability +------------ + +Whether a certain backend can be used depends on the availability of the underlying tensor network +library. This can easily be asserted at runtime using the following indicators: + +.. autoclass:: HAS_QUIMB + +.. autoclass:: HAS_TENPY + +Backends +-------- + +Depending on the availability (see above), the following backends are available: + +.. autosummary:: + :toctree: + + quimb_tebd + quimb_layers + quimb_circuit + tenpy_tebd + tenpy_layers + +Interface +--------- + +The interface implemented by any one of these optional backends is made up of the following classes: + +.. autoclass:: Evolver + +.. autoclass:: State +""" + +from qiskit.utils import LazyImportTester as _LazyImportTester + +from .interface import Evolver, State + +HAS_QUIMB = _LazyImportTester("quimb", install="pip install qiskit-addon-mpf[quimb]") +"""Indicates whether the optional :external:mod:`quimb` dependency is installed.""" + +HAS_TENPY = _LazyImportTester("tenpy", install="pip install qiskit-addon-mpf[tenpy]") +"""Indicates whether the optional :external:mod:`tenpy` dependency is installed.""" + +__all__ = [ + "Evolver", + "State", +] diff --git a/qiskit_addon_mpf/backends/interface.py b/qiskit_addon_mpf/backends/interface.py new file mode 100644 index 0000000..1d17f49 --- /dev/null +++ b/qiskit_addon_mpf/backends/interface.py @@ -0,0 +1,117 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""The interface for :class:`.DynamicMPF` backends.""" + +from __future__ import annotations + +from abc import ABC, abstractmethod +from typing import Any + + +class Evolver(ABC): + """The interface for the time-evolution algorithms used within :class:`.DynamicMPF`. + + This time-evolution interface is used by the :attr:`.DynamicMPF.lhs` and :attr:`.DynamicMPF.rhs` + and should time-evolve a :class:`.State` object under its hood. The exact mechanism of the + algorithm is described in more detail in :class:`.DynamicMPF`, :class:`.State`, and + :func:`.setup_dynamic_lse`. + """ + + @abstractmethod + def step(self) -> None: + """Perform a single time step of this time-evolution algorithm. + + This should act on the internally referenced :class:`.State` object (for which no name is + prescribed by this interface). Whether this time-evolution algorithm instance should evolve + the :class:`.State` from the left- or right-hand side, depends on the value of + :attr:`.conjugate`. + """ + + @property + @abstractmethod + def evolved_time(self) -> float: + """Returns the current evolution time.""" + + @property + @abstractmethod + def conjugate(self) -> bool: + """Returns whether this time-evolver instance acts on the right-hand side.""" + + @conjugate.setter + @abstractmethod + def conjugate(self, conjugate: bool) -> None: ... # pragma: no cover + + +class State(ABC): + """The interface for the :attr:`.DynamicMPF.evolution_state`. + + This time-evolution state is shared between the LHS and RHS :class:`.Evolver` instances of the + :class:`.DynamicMPF` instance. In most cases where a concrete backend implementing this + interface is based on tensor networks, this state will be a matrix product operator (MPO). + This is because most time-evolution algorithms would normally evolve a matrix product state + (MPS) as shown pictorially below, where time evolution blocks (``U#``) are successively applied + to a 1-dimensional MPS (``S#``). Here, the tensor network grows towards the right as time goes + on. + + .. code-block:: + + MPS Evolution + + S0┄┄┲━━━━┱┄┄┄┄┄┄┄┄┲━━━━┱┄ + │ ┃ U1 ┃ ┃ U5 ┃ + S1┄┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄ + │ ┃ U3 ┃ + S2┄┄┲━━━━┱┄┺━━━━┹┄┲━━━━┱┄ ... + │ ┃ U2 ┃ ┃ U6 ┃ + S3┄┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄ + │ ┃ U4 ┃ + S4┄┄┄┄┄┄┄┄┄┺━━━━┹┄┄┄┄┄┄┄┄ + + However, in our case, we want two time-evolution engines to share a single state. In order to + achieve that, we can have one of them evolve the state from the right (just as before, ``U#``), + but have the second one evolve the state from the left (``V#``). This requires the state to also + have bonds going of in that direction, rendering it a 2-dimensional MPO (``M#``) rather than the + 1-dimensional MPS from before. + + .. code-block:: + + MPO Evolution + + ┄┲━━━━┱┄┄┄┄┄┄┄┄┲━━━━┱┄┄M0┄┄┲━━━━┱┄┄┄┄┄┄┄┄┲━━━━┱┄ + ┃ V5 ┃ ┃ V1 ┃ │ ┃ U1 ┃ ┃ U5 ┃ + ┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄┄M1┄┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄ + ┃ V3 ┃ │ ┃ U3 ┃ + ... ┄┲━━━━┱┄┺━━━━┹┄┲━━━━┱┄┄M2┄┄┲━━━━┱┄┺━━━━┹┄┲━━━━┱┄ ... + ┃ V6 ┃ ┃ V2 ┃ │ ┃ U2 ┃ ┃ U6 ┃ + ┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄┄M3┄┄┺━━━━┹┄┲━━━━┱┄┺━━━━┹┄ + ┃ V4 ┃ │ ┃ U4 ┃ + ┄┄┄┄┄┄┄┄┺━━━━┹┄┄┄┄┄┄┄┄┄M4┄┄┄┄┄┄┄┄┄┺━━━━┹┄┄┄┄┄┄┄┄ + """ + + @abstractmethod + def overlap(self, initial_state: Any) -> complex: + """Compute the overlap of this state with the provided initial state. + + .. warning:: + A concrete implementation of this method should raise a :class:`TypeError` if the + provided ``initial_state`` object is not supported by the implementing backend. + + Args: + initial_state: the initial state with which to compute the overlap. + + Raises: + TypeError: if the provided initial state has an incompatible type. + + Returns: + The overlap of this state with the provided one. + """ diff --git a/qiskit_addon_mpf/backends/quimb_circuit/__init__.py b/qiskit_addon_mpf/backends/quimb_circuit/__init__.py new file mode 100644 index 0000000..233b78b --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_circuit/__init__.py @@ -0,0 +1,110 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A circuit-based time-evolution backend using :external:mod:`quimb`. + +.. currentmodule:: qiskit_addon_mpf.backends.quimb_circuit + +.. warning:: + This backend is only available if the optional dependencies have been installed: + + .. code-block:: + + pip install "qiskit-addon-mpf[quimb]" + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + :template: autosummary/class_without_inherited_members.rst + + CircuitEvolver + CircuitState + + +Underlying method +----------------- + +Quimb boasts direct support for the simulation of quantum circuits in the form of its tensor-network +based :external:class:`quimb.tensor.Circuit` representation. We can leverage this, to bypass any +explicit time-evolution algorithm and instead directly encode the time-evolution in a +:external:class:`~qiskit.circuit.QuantumCircuit` and use :external:mod:`quimb` to compute the +overlap between two such circuits. For more information, check out their guide on +:external:std:doc:`tensor-circuit`. + + +Code example +------------ + +This section shows a simple example to get you started with using this backend. The example shows +how to create the three factory functions required for the :func:`.setup_dynamic_lse`. + +The :class:`.IdentityStateFactory` protocol is already fulfilled by the +:class:`~.quimb_circuit.CircuitState` constructor, rendering the ``identity_factory`` argument +trivial: + +>>> from qiskit_addon_mpf.backends.quimb_circuit import CircuitState +>>> identity_factory = CircuitState + +The setup of the :class:`~.quimb_circuit.CircuitEvolver` is slightly more involved. It requires a +**parametrized** :external:class:`~qiskit.circuit.QuantumCircuit` object as its input where the +:class:`~qiskit.circuit.Parameter` should take the place of the Trotter methods time step (``dt``). + +To show how such a parametrized Trotter circuit template is constructed, we reuse the same +Hamiltonian and second-order Suzuki-Trotter formula as in :mod:`.quimb_layers`. + +>>> from qiskit.quantum_info import SparsePauliOp +>>> hamil = SparsePauliOp.from_sparse_list( +... [("ZZ", (i, i+1), 1.0) for i in range(0, 9, 2)] + +... [("Z", (i,), 0.5) for i in range(10)] + +... [("ZZ", (i, i+1), 1.0) for i in range(1, 9, 2)] + +... [("X", (i,), 0.25) for i in range(10)], +... num_qubits=10, +... ) + +But this time, we specify a :class:`~qiskit.circuit.Parameter` as the ``time`` argument when +constructing the actual circuits. + +>>> from functools import partial +>>> from qiskit.circuit import Parameter +>>> from qiskit.synthesis import SuzukiTrotter +>>> from qiskit_addon_mpf.backends.quimb_circuit import CircuitEvolver +>>> from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit +>>> dt = Parameter("dt") +>>> suzuki_2 = generate_time_evolution_circuit(hamil, time=dt, synthesis=SuzukiTrotter(order=2)) +>>> approx_evolver_factory = partial(CircuitEvolver, circuit=suzuki_2) + +.. caution:: + It is **necessary** that the name of the :class:`~qiskit.circuit.Parameter` is ``dt``! + +We can choose a higher order Trotter formula for the ``exact_evolver_factory``. But note, that we +must once again use a parametrized circuit, even if we immediately bind its parameter when +constructing the ``partial`` function. + +>>> suzuki_4 = generate_time_evolution_circuit(hamil, time=dt, synthesis=SuzukiTrotter(order=4)) +>>> exact_evolver_factory = partial(CircuitEvolver, circuit=suzuki_4, dt=0.05) + +These factory functions may now be used to run the :func:`.setup_dynamic_lse`. Refer to its +documentation for more details on that. +""" + +# ruff: noqa: E402 +from .. import HAS_QUIMB + +HAS_QUIMB.require_now(__name__) + +from .evolver import CircuitEvolver +from .state import CircuitState + +__all__ = [ + "CircuitEvolver", + "CircuitState", +] diff --git a/qiskit_addon_mpf/backends/quimb_circuit/evolver.py b/qiskit_addon_mpf/backends/quimb_circuit/evolver.py new file mode 100644 index 0000000..dec35d2 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_circuit/evolver.py @@ -0,0 +1,92 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A time-evolution engine based on quantum circuits.""" + +from __future__ import annotations + +from qiskit.circuit import QuantumCircuit +from qiskit_quimb import quimb_circuit + +from .. import Evolver +from .state import CircuitState + + +class CircuitEvolver(Evolver): + """A time-evolution engine based on quantum circuits. + + This algorithm performs time-evolution by means of successively applying a quantum circuit + corresponding to a single Trotter step to its internal state. More specifically, it builds out a + tensor network in the :class:`~.quimb_circuit.CircuitState`. As required by the + :class:`.DynamicMPF` algorithm, it tracks a left- and right-hand side of the time-evolution for + computing the overlap of two circuits. Depending on :attr:`conjugate`, an instance of this + engine will apply the quantum gates of its template circuit to the corresponding side (see + :mod:`.quimb_circuit` for more details). + """ + + def __init__(self, evolution_state: CircuitState, circuit: QuantumCircuit, dt: float) -> None: + """Initialize a :class:`CircuitEvolver` instance. + + Args: + evolution_state: a reference to the time-evolution state. + circuit: the template circuit encoding the time-evolution of a single Trotter step. This + circuit **must** be parametrized (see :external:class:`~qiskit.circuit.Parameter` in + place of the Trotter methods time step. This parameter must be named ``dt``. + dt: the time step that will be used and later bound to the + :external:class:`~qiskit.circuit.Parameter` of the ``circuit`` object. + """ + self.evolution_state = evolution_state + """The time-evolution state (see also :attr:`.DynamicMPF.evolution_state`).""" + self.evolved_time = 0 + self.dt = dt + self.circuit = quimb_circuit(circuit.assign_parameters({"dt": dt}, inplace=False)) + """The parameterized :external:class:`~qiskit.circuit.QuantumCircuit` describing the Trotter + step.""" + self._conjugate = False + + @property + def conjugate(self) -> bool: + """Returns whether this time-evolver instance acts on the right-hand side.""" + return self._conjugate + + @conjugate.setter + def conjugate(self, conjugate: bool) -> None: + self._conjugate = conjugate + + @property + def evolved_time(self) -> float: + """Returns the current evolution time.""" + return self._evolved_time + + @evolved_time.setter + def evolved_time(self, evolved_time: float) -> None: + self._evolved_time = evolved_time + + def step(self) -> None: + """Perform a single time step of TEBD. + + This will apply the gates of the :attr:`circuit` to the :attr:`evolution_state`. If + :attr:`conjugate` is ``True``, it applies to :attr:`.CircuitState.lhs`, otherwise to + :attr:`.CircuitState.rhs`. + """ + self.evolved_time += self.dt + + if self.conjugate: + if self.evolution_state.lhs is None: + self.evolution_state.lhs = self.circuit.copy() + else: + self.evolution_state.lhs.apply_gates(self.circuit.gates) + else: + if self.evolution_state.rhs is None: + self.evolution_state.rhs = self.circuit.copy() + else: + self.evolution_state.rhs.apply_gates(self.circuit.gates) diff --git a/qiskit_addon_mpf/backends/quimb_circuit/state.py b/qiskit_addon_mpf/backends/quimb_circuit/state.py new file mode 100644 index 0000000..870b7c7 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_circuit/state.py @@ -0,0 +1,78 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A circuit-based MPO-like time-evolution state based on quimb.""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +from qiskit.circuit import QuantumCircuit +from qiskit_quimb import quimb_circuit +from quimb.tensor import Circuit, TensorNetwork + +from .. import State + + +class CircuitState(State): + """An MPO-like representation of a time-evolution state based on quantum circuits. + + This time-evolution state can be evolved on its left- and right-hand side as required by the + :class:`.DynamicMPF` algorithm. + """ + + def __init__(self) -> None: + """Initialize a :class:`CircuitState` instance.""" + self.lhs: Circuit | None = None + """The left-hand side circuit in form of a tensor network.""" + self.rhs: Circuit | None = None + """The right-hand side circuit in form of a tensor network.""" + + def overlap(self, initial_state: Any) -> complex: + """Compute the overlap of this state with the provided initial state. + + .. warning:: + This implementation only supports instances of + :external:class:`qiskit.circuit.QuantumCircuit` for ``initial_state``. + + Args: + initial_state: the initial state with which to compute the overlap. + + Raises: + TypeError: if the provided initial state has an incompatible type. + + Returns: + The overlap of this state with the provided one. + """ + if not isinstance(initial_state, QuantumCircuit): + raise TypeError( + "CircuitState.overlap is only implemented for qiskit.QuantumCircuit! " + "But not for states of type '%s'", + type(initial_state), + ) + + if self.lhs is None or self.rhs is None: + raise RuntimeError("You must evolve the state before an overlap can be computed!") + + lhs = quimb_circuit(initial_state) + lhs.apply_gates(self.lhs.gates) + + rhs = quimb_circuit(initial_state) + rhs.apply_gates(self.rhs.gates) + + # TODO: find a good way to inject arguments into .contract() below + # For example, specifcying backend="jax" would allow us to run this on a GPU (if available + # and installed properly). + ovlp = TensorNetwork((lhs.psi.H, rhs.psi)).contract() + + return float(np.abs(ovlp) ** 2) diff --git a/qiskit_addon_mpf/backends/quimb_layers/__init__.py b/qiskit_addon_mpf/backends/quimb_layers/__init__.py new file mode 100644 index 0000000..a2aa193 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_layers/__init__.py @@ -0,0 +1,230 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A layer-wise time-evolution backend using :external:mod:`quimb`. + +.. currentmodule:: qiskit_addon_mpf.backends.quimb_layers + +.. warning:: + This backend is only available if the optional dependencies have been installed: + + .. code-block:: + + pip install "qiskit-addon-mpf[quimb]" + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + :template: autosummary/class_without_inherited_members.rst + + LayerwiseEvolver + LayerModel + + +Underlying method +----------------- + +This module provides a time-evolution backend similar to the TEBD-based one provided by the +:mod:`~qiskit_addon_mpf.backends.quimb_tebd` module. The main difference is that this module gives +the user full flexibility for defining their product formulas, thereby not limiting them to the +options built into the :external:mod:`quimb` library. + +At its core, the algorithm provided by this module is still a TEBD [1] algorithm. However, rather +than enforcing the alternating updates to the even and odd bonds of the time-evolution state (see +also :meth:`.quimb_tebd.TEBDEvolver.sweep`) this implementation outsources the responsibility of +updating bonds in alternating fashion to the definition of multiple time-evolution **layers**. + +This is best explained with an example. Let us assume, we have some generic Hamiltonian acting on a +1-dimensional chain of sites. + +.. hint:: + Below we are very deliberate about the order of the Hamiltonian's Pauli terms because this + directly impacts the structure of the time-evolution circuit later on. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from qiskit.quantum_info import SparsePauliOp + >>> hamil = SparsePauliOp.from_sparse_list( + ... [("ZZ", (i, i+1), 1.0) for i in range(1, 9, 2)] + + ... [("Z", (i,), 0.5) for i in range(10)] + + ... [("ZZ", (i, i+1), 1.0) for i in range(0, 9, 2)], + ... num_qubits=10, + ... ) + +Let us now inspect the time-evolution circuit of this Hamiltonian using a second-order +Suzuki-Trotter formula. + +.. plot:: + :context: close-figs + :include-source: + + >>> from qiskit.synthesis import SuzukiTrotter + >>> from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit + >>> circuit = generate_time_evolution_circuit(hamil, time=1.0, synthesis=SuzukiTrotter(order=2)) + >>> circuit.draw("mpl") # doctest: +ELLIPSIS +
+ +In the circuit above, we can clearly identify its layer-wise structure. We can emphasize this +further, by splitting the circuit into multiple layers as shown below (we recombine the ``layers`` +into a single circuit with barriers between them to ease the visualization). + +.. plot:: + :context: close-figs + :include-source: + + >>> from qiskit_addon_utils.slicing import combine_slices, slice_by_gate_types + >>> layers = slice_by_gate_types(circuit) + >>> combine_slices(layers, include_barriers=True).draw("mpl") # doctest: +ELLIPSIS +
+ +.. hint:: + The asymmetry of the central layers is a result of the implementation of Qiskit's + :external:class:`~qiskit.synthesis.SuzukiTrotter` formula. In its second-order form, it combines + the two half-time evolutions of the final term in the Hamiltonian into a single one of twice the + length. We could transpile this circuit to collapse all such subequent gates in the central two + layers (just like the last one), but for the sake of simplicity of this example, we will not do + that here. + +It is not possible to instruct Quimb's TEBD algorithm to simulate the exact structure of the circuit +shown above. The reason for that is a limitation in its interface, as it only accepts the full +Hamiltonian to be provided which is then time-evolved using pre-defined Trotter formulas. However, +in doing so it does not treat the order of the Pauli terms in a Hamiltonian with any significance +(like we do here). + +If one wants to compute the dynamic MPF coefficients of a time-evolution employing a product formula +structure other than the ones implemented in Quimb (like the example above), then one can use the +time-evolution algorithm provided by this module. Rather than taking a single monolithic Hamiltonian +whose time-evolution is to be modeled, the :class:`~.quimb_layers.LayerwiseEvolver` accepts a list +of :class:`~.quimb_layers.LayerModel` objects, each one describing an individual layer of the +product formula. This gives the user full flexibility in defining the Trotter decomposition down to +the most granular level. + +However, caution must be applied to ensure that the property of TEBD to update even and odd bonds in +an alternating manner is still guaranteed. Luckily, for quantum circuits consisting of at most +two-qubit gates, this property is satisfied by construction. + + +Code example +------------ + +In this section, we build up on the example above and show how to take a custom Trotter formula and +use it to construct a :class:`~.quimb_layers.LayerwiseEvolver` which can be used to replace the +:class:`.quimb_tebd.TEBDEvolver` in the workflow described in :mod:`.quimb_tebd`. + +.. hint:: + The overall workflow of using this module is the same as of the :mod:`.quimb_tebd` module, so be + sure to read those instructions as well. + +Simply put, we must convert each one of the circuit ``layers`` (see above) into a +:class:`~.quimb_layers.LayerModel` instance. For this purpose, we can use its +:meth:`~.quimb_layers.LayerModel.from_quantum_circuit` method. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from qiskit_addon_mpf.backends.quimb_layers import LayerModel + >>> model0 = LayerModel.from_quantum_circuit(layers[0]) + >>> layer_models = [model0] + +In the code above you can see how simple the conversion is for layers which contain only two-qubit +gates acting on mutually exclusive qubits (which layers of depth 1 guarantee). + +However, we must be more careful with layers including single-qubit gates. The reason for that is +that the TEBD algorithm underlying the :class:`~.quimb_layers.LayerwiseEvolver` must update even and +odd bonds in an alternating manner. And because single-qubit gates are not applied on a site, but +instead are split in half and applied to the bonds on either side, a layer of single-qubit gates +acting on all qubits would break this assumption. + +To circumvent this problem, we can take any layer consisting of only single-qubit gates, and apply +twice (once on the even and once on the odd bonds) while scaling each one with a factor of ``0.5``. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> model1a = LayerModel.from_quantum_circuit(layers[1], keep_only_odd=True, scaling_factor=0.5) + >>> model1b = LayerModel.from_quantum_circuit(layers[1], keep_only_odd=False, scaling_factor=0.5) + >>> layer_models.extend([model1a, model1b]) + +Now that we know how to treat layers consisting of two-qubit and single-qubit gates, we can +transform the remaining layers. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> for layer in layers[2:]: + ... num_qubits = len(layer.data[0].qubits) + ... if num_qubits == 2: + ... layer_models.append(LayerModel.from_quantum_circuit(layer)) + ... else: + ... layer_models.append( + ... LayerModel.from_quantum_circuit(layer, keep_only_odd=True, scaling_factor=0.5) + ... ) + ... layer_models.append( + ... LayerModel.from_quantum_circuit(layer, keep_only_odd=False, scaling_factor=0.5) + ... ) + >>> assert len(layer_models) == 8 + +In the end, we have 8 :class:`~.quimb_layers.LayerModel`'s, one for each of the 4 two-qubit layers, +and two for each of the 2 single-qubit layers. + +Finally, we can define our :class:`.ApproxEvolverFactory` protocol to be used within the +:func:`.setup_dynamic_lse` function. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from functools import partial + >>> from qiskit_addon_mpf.backends.quimb_layers import LayerwiseEvolver + >>> approx_evolver_factory = partial( + ... LayerwiseEvolver, + ... layers=layer_models, + ... split_opts={"max_bond": 10, "cutoff": 1e-5}, + ... ) + +.. caution:: + It should be noted, that in this workflow we have not yet fixed the time step used by the Trotter + formula. We have also only set up a single repetition of the Trotter formula as the rest will be + done by the internal :class:`.DynamicMPF` algorithm, executed during :func:`.setup_dynamic_lse`. + +Of course, you could also use this to specify a :class:`.ExactEvolverFactory`. But you can also +mix-and-match a :class:`.quimb_layers.LayerwiseEvolver` with a :class:`.quimb_tebd.TEBDEvolver`. + + +Resources +--------- + +[1]: https://en.wikipedia.org/wiki/Time-evolving_block_decimation +""" + +# ruff: noqa: E402 +from .. import HAS_QUIMB + +HAS_QUIMB.require_now(__name__) + +from .evolver import LayerwiseEvolver +from .model import LayerModel + +__all__ = [ + "LayerModel", + "LayerwiseEvolver", +] diff --git a/qiskit_addon_mpf/backends/quimb_layers/evolver.py b/qiskit_addon_mpf/backends/quimb_layers/evolver.py new file mode 100644 index 0000000..479dd51 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_layers/evolver.py @@ -0,0 +1,75 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A layer-wise time-evolver using quimb.""" + +from __future__ import annotations + +from .. import quimb_tebd +from .model import LayerModel + + +class LayerwiseEvolver(quimb_tebd.TEBDEvolver): + """A special case of the :class:`~.quimb_tebd.TEBDEvolver` based on layer-wise evolution models. + + As also explained in :mod:`.quimb_layers`, this implementation extracts the alternating even/odd + bond updates implemented inside of the original :external:class:`quimb.tensor.TEBD` to become + the end users responsibility. It does so, by replacing the single Hamiltonian provided to the + :class:`~.quimb_tebd.TEBDEvolver` instance with a sequence of :class:`~.quimb_layers.LayerModel` + instances. Every single instance of these encodes a single **layer** of interactions. These + should enforce the alternating updates of even and odd bonds of the underlying tensor network. + + The motivation for this more complicated interface is that is provides a lot more flexbility and + enables users to define custom Trotter product formulas rather than being limited to the ones + implemented by ``quimb`` directly. + """ + + def __init__( + self, evolution_state: quimb_tebd.MPOState, layers: list[LayerModel], *args, **kwargs + ) -> None: + """Initialize a :class:`LayerwiseEvolver` instance. + + Args: + evolution_state: forwarded to :class:`~.quimb_tebd.TEBDEvolver`. Please refer to its + documentation for more details. + layers: the list of models describing single layers of interactions. See above as well + as the explanations provided in :mod:`.quimb_layers`. + args: any further positional arguments will be forwarded to the + :class:`~.quimb_tebd.TEBDEvolver` constructor. + kwargs: any further keyword arguments will be forwarded to the + :class:`~.quimb_tebd.TEBDEvolver` constructor. + """ + super().__init__(evolution_state, layers[0], *args, **kwargs) + self.layers = layers + """The layers of interactions used to implement the time-evolution.""" + + def step(self) -> None: + # pylint: disable=attribute-defined-outside-init + # NOTE: Somehow, pylint does not properly pick up on the attributes defined in the external + # base classes. + """Perform a single time step of TEBD. + + This will iterate over the :attr:`layers` and apply their interaction to the internal state. + """ + dt = self._dt + + for layer in self.layers: + self.H = layer + for i in range(self.L): + sites = (i, (i + 1) % self.L) + gate = self._get_gate_from_ham(1.0, sites) + if gate is None: + continue + self._pt.gate_split_(gate, sites, conj=self.conjugate, **self.split_opts) + + self.t += dt + self._err += float("NaN") diff --git a/qiskit_addon_mpf/backends/quimb_layers/model.py b/qiskit_addon_mpf/backends/quimb_layers/model.py new file mode 100644 index 0000000..468fe67 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_layers/model.py @@ -0,0 +1,150 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A quimb-based model for describing a single layer of interactions.""" + +from __future__ import annotations + +from typing import cast + +import numpy as np +from qiskit.circuit import QuantumCircuit +from quimb.gen.operators import pauli +from quimb.tensor import LocalHam1D + + +class LayerModel(LocalHam1D): + """A model for representing a layer of time-evolution interactions. + + Essentially, this class is a simple wrapper of :external:class:`quimb.tensor.LocalHam1D`. Its + main purpose is to provide a simple interface for constructing a Quimb-compatible Hamiltonian + from Qiskit objects. + """ + + def __init__(self, L, H2, H1=None, cyclic=False, keep_only_odd=None) -> None: + """Initialize a :class:`LayerModel` instance. + + Most of the arguments below are simply forwarded to + :external:class:`quimb.tensor.LocalHam1D` so check out its documentation for more details. + + Args: + L: the number of qubits. + H2: the two-site interactions. + H1: the optional on-site interactions. + cyclic: whether to apply periodic boundary conditions. + keep_only_odd: whether to keep only odd bond interactions. For more details see + :attr:`keep_only_odd`. + """ + super().__init__(L, H2, H1, cyclic) + self.keep_only_odd = keep_only_odd + """Whether to keep only interactions on bonds with odd indices.""" + + def get_gate_expm(self, where: tuple[int, int], x: float) -> np.ndarray | None: + """Get the local term at the sites ``where``, matrix exponentiated by ``x``. + + If ``where`` applies to an even bond index and :attr:`keep_only_odd` is ``True``, this + method will return ``None``. + + Args: + where: the pair of site indices of the local term to get. This identifies the bond + index. + x: the value with which to matrix exponentiate the interaction term. + + Returns: + The interaction in terms of an array or ``None`` depending on :attr:`keep_only_odd` (see + above). + """ + if self.keep_only_odd is not None and where[0] % 2 - self.keep_only_odd: + return None + try: + return cast(np.ndarray, self._expm_cached(self.get_gate(where), x)) + except KeyError: + return None + + @classmethod + def from_quantum_circuit( + cls, + circuit: QuantumCircuit, + *, + scaling_factor: float = 1.0, + keep_only_odd: bool | None = None, + **kwargs, + ) -> LayerModel: + """Construct a :class:`LayerModel` from a :external:class:`~qiskit.circuit.QuantumCircuit`. + + You can see an example of this function in action in the docs of :mod:`quimb_layers`. + + Args: + circuit: the quantum circuit to parse. + scaling_factor: a factor with which to scale the term strengths. This can be used to + apply (for example) a time step scaling factor. It may also be used (e.g.) to split + onsite terms into two layers (even and odd) with $0.5$ of the strength, each. + keep_only_odd: the value to use for :attr:`keep_only_odd`. + kwargs: any additional keyword arguments to pass to the :class:`LayerModel` constructor. + + Returns: + A new LayerModel instance. + + Raises: + NotImplementedError: if an unsupported quantum gate is encountered. + """ + H2: dict[tuple[int, int] | None, np.ndarray] = {} + H1: dict[int, np.ndarray] = {} + paulis_cache: dict[str, np.ndarray] = {} + + for instruction in circuit.data: + op = instruction.operation + sites = tuple(circuit.find_bit(qubit)[0] for qubit in instruction.qubits) + + # NOTE: the hard-coded scaling factors below account for the Pauli matrix conversion + if op.name == "rzz": + term = paulis_cache.get(op.name, None) + if term is None: + paulis_cache[op.name] = pauli("Z") & pauli("Z") + term = paulis_cache[op.name] + if sites in H2: + H2[sites] += 0.5 * scaling_factor * op.params[0] * term + else: + H2[sites] = 0.5 * scaling_factor * op.params[0] * term + elif op.name == "xx_plus_yy": + term = paulis_cache.get(op.name, None) + if term is None: + paulis_cache["rxx"] = pauli("X") & pauli("X") + paulis_cache["ryy"] = pauli("Y") & pauli("Y") + paulis_cache[op.name] = paulis_cache["rxx"] + paulis_cache["ryy"] + term = paulis_cache[op.name] + if sites in H2: + H2[sites] += 0.25 * scaling_factor * op.params[0] * term + else: + H2[sites] = 0.25 * scaling_factor * op.params[0] * term + elif op.name == "rz": + term = paulis_cache.get(op.name, None) + if term is None: + paulis_cache[op.name] = pauli("Z") + term = paulis_cache[op.name] + if sites[0] in H1: + H1[sites[0]] += 2.0 * scaling_factor * op.params[0] * term + else: + H1[sites[0]] = 2.0 * scaling_factor * op.params[0] * term + else: + raise NotImplementedError(f"Cannot handle gate of type {op.name}") + + if len(H2) == 0: + H2[None] = np.zeros((4, 4)) + + return cls( + circuit.num_qubits, + H2, + H1, + keep_only_odd=keep_only_odd, + **kwargs, + ) diff --git a/qiskit_addon_mpf/backends/quimb_tebd/__init__.py b/qiskit_addon_mpf/backends/quimb_tebd/__init__.py new file mode 100644 index 0000000..df14a4c --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_tebd/__init__.py @@ -0,0 +1,139 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A :external:mod:`quimb`-based TEBD backend. + +.. currentmodule:: qiskit_addon_mpf.backends.quimb_tebd + +.. warning:: + This backend is only available if the optional dependencies have been installed: + + .. code-block:: + + pip install "qiskit-addon-mpf[quimb]" + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + :template: autosummary/class_without_inherited_members.rst + + TEBDEvolver + MPOState + +Underlying method +----------------- + +This module provides a time-evolution backend for computing dynamic MPF coefficients based on the +time-evolving block decimation (TEBD) algorithm [1] implemented in the :external:mod:`quimb` tensor +network library. + +The classes provided by this module serve two purposes: + +1. Connecting :external:mod:`quimb`'s implementation to the interface set out by + :mod:`qiskit_addon_mpf.backends`. +2. Extending :external:mod:`quimb`'s TEBD implementation to handle an internal MPO (rather than + MPS) state (see also :class:`.State` for more details). + + +In the simplest sense, this module provides a straight-forward extension of the TEBD algorithm to +evolve an internal MPO state. +As such, if you wish to use this backend for your dynamic MPF algorithm, you must encode the +Hamiltonian that you wish to time-evolve, in a :external:mod:`quimb`-native form. To be more +concrete, the :class:`~qiskit_addon_mpf.backends.quimb_tebd.TEBDEvolver` class (which is a subclass +of :external:class:`quimb.tensor.TEBD`) works with a Hamiltonian in the form of a +:external:class:`quimb.tensor.LocalHam1D`. Quimb provides a number of convenience methods for +constructing such Hamiltonians in its :external:mod:`quimb.tensor.tensor_builder` module. +If none of those fulfill your needs, you can consider using the +:class:`~qiskit_addon_mpf.backends.quimb_layers.LayerModel` class which implements some conversion +methods from Qiskit-native objects. + +Code example +------------ + +This section shows a simple example to get you started with using this backend. The example shows +how to create the three factory functions required for the :func:`.setup_dynamic_lse`. + +First, we create the ``identity_factory`` which has to match the :class:`.IdentityStateFactory` +protocol. We do so simply by using the :external:func:`quimb.tensor.MPO_identity` function and +wrapping the resulting :external:class:`quimb.tensor.MatrixProductOperator` with our custom +:class:`~.quimb_tebd.MPOState` interface. + +>>> from qiskit_addon_mpf.backends.quimb_tebd import MPOState +>>> from quimb.tensor import MPO_identity +>>> num_qubits = 10 +>>> identity_factory = lambda: MPOState(MPO_identity(num_qubits)) + +Next, before being able to define the :class:`.ExactEvolverFactory` and +:class:`.ApproxEvolverFactory` protocols, we must define the Hamiltonian which we would like to +time-evolve. Here, we simply choose one of :external:mod:`quimb`'s convenience methods. + +>>> from quimb.tensor import ham_1d_heis +>>> hamil = ham_1d_heis(num_qubits, 0.8, 0.3, cyclic=False) + +We can now construct the exact and approximate time-evolution instance factories. To do so, we can +simply use :func:`functools.partial` to bind the pre-defined values of the +:class:`~qiskit_addon_mpf.backends.quimb_tebd.TEBDEvolver` initializer, reducing it to the correct +interface as expected by the :class:`.ExactEvolverFactory` and :class:`.ApproxEvolverFactory` +protocols, respectively. + +>>> from functools import partial +>>> from qiskit_addon_mpf.backends.quimb_tebd import TEBDEvolver +>>> exact_evolver_factory = partial( +... TEBDEvolver, +... H=hamil, +... dt=0.05, +... order=4, +... ) + +Notice, how we have fixed the ``dt`` value to a small time step and have used a higher-order +Suzuki-Trotter decomposition to mimic the exact time evolution above. + +Below, we do not fix the ``dt`` value and use only a second-order Suzuki-Trotter formula for the +approximate time evolution. Additionally, we also specify some truncation settings. + +>>> approx_evolver_factory = partial( +... TEBDEvolver, +... H=hamil, +... order=2, +... split_opts={"max_bond": 10, "cutoff": 1e-5}, +... ) + +Of course, you are not limited to the examples shown here, and we encourage you to play around with +the other settings provided by the :external:class:`quimb.tensor.TEBD` implementation. + +Limitations +----------- + +Finally, we point out a few known limitations on what kind of Hamiltonians can be treated by this +backend: + +* all interactions must be 1-dimensional +* the interactions must be acylic + +Resources +--------- + +[1]: https://en.wikipedia.org/wiki/Time-evolving_block_decimation +""" + +# ruff: noqa: E402 +from .. import HAS_QUIMB + +HAS_QUIMB.require_now(__name__) + +from .evolver import TEBDEvolver +from .state import MPOState + +__all__ = [ + "TEBDEvolver", + "MPOState", +] diff --git a/qiskit_addon_mpf/backends/quimb_tebd/evolver.py b/qiskit_addon_mpf/backends/quimb_tebd/evolver.py new file mode 100644 index 0000000..7e7385e --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_tebd/evolver.py @@ -0,0 +1,160 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A quimb-based TEBD algorithm for evolving an internal MPO.""" + +from __future__ import annotations + +from typing import Literal, cast + +from quimb.tensor import TEBD + +from .. import Evolver +from .state import MPOState + + +class TEBDEvolver(TEBD, Evolver): + """A TEBD algorithm for evolving an internal MPO. + + As discussed in more detail in :mod:`~qiskit_addon_mpf.backends.quimb_tebd`, this extension of + ``quimb``'s existing :external:class:`quimb.tensor.TEBD` implementation time-evolves an internal + matrix product operator (MPO) rather than the conventional matrix product state (MPS). + + More concretely, the internal object is expected to be an :class:`~.quimb_tebd.MPOState`. + + .. warning:: + The API of this class is actually much larger than shown here, because it inherits additional + methods from the :external:class:`quimb.tensor.TEBD` base class. However, we do not duplicate + that API here. + """ + + def __init__(self, evolution_state: MPOState, *args, order: int = 2, **kwargs) -> None: + """Initialize a :class:`TEBDEvolver` instance. + + Args: + evolution_state: a reference to the time-evolution state. This overwrites the ``p0`` + argument of the underlying :external:class:`quimb.tensor.TEBD` class. + + .. warning:: + In contrast to the default behavior, this state will **NOT** be canonicalized. + Instead, it is taken as is and is kept **by reference** (i.e. no copy is + created). This ensures that the same object can be shared between two instances + of this class, as required by the :class:`.DynamicMPF` algorithm. + + args: any further positional arguments will be forwarded to the + :external:class:`quimb.tensor.TEBD` constructor. + order: the order of the builtin Suzuki-Trotter formula to use during time evolution. + This will be the value forwarded to the :external:meth:`quimb.tensor.TEBD.step` + method. + kwargs: any further keyword arguments will be forwarded to the + :external:class:`quimb.tensor.TEBD` constructor. + """ + super().__init__(evolution_state, *args, **kwargs) + # WARNING: we must forcefully overwrite self._pt to ensure that p0 is kept by reference! + self._pt = evolution_state + self._order = order + self._conjugate = False + + @property + def evolved_time(self) -> float: + """Returns the current evolution time.""" + return cast(float, self.t) + + @property + def conjugate(self) -> bool: + """Returns whether this time-evolver instance acts on the right-hand side.""" + return self._conjugate + + @conjugate.setter + def conjugate(self, conjugate: bool) -> None: + self._conjugate = conjugate + + def step(self) -> None: + """Perform a single time step of TEBD. + + This essentially calls :external:meth:`quimb.tensor.TEBD.step` and forwards the value of + the ``order`` attribute that was provided upon construction. + """ + TEBD.step(self, order=self._order) + + def sweep( + self, + direction: Literal["left", "right"], + dt_frac: float, + dt: float | None = None, + queue: bool = False, + ) -> None: + """Perform a single sweep of the TEBD algorithm [1]. + + The TEBD algorithm updates the even and odd bonds of the underlying tensor network in + alternating fashion. In the implementation of the :external:class:`quimb.tensor.TEBD` base + class, this is realized in the form of alternating "sweeps" in left and right directions + over the internal state. + + We are overwriting the behavior of this method in this subclass, in order to call the + specialized :meth:`~.quimb_tebd.MPOState.gate_split` method. + + Args: + direction: the direction of the sweep. This must be either of the literal strings, + ``"left"`` or ``"right"``. + dt_frac: what fraction of the internal time step (``dt``) to time-evolve for. This is + how any builtin Suzuki-Trotter formula specifies its splitting. + dt: an optional value to overwrite the internal time step. + queue: setting this to ``True`` will raise a :class:`NotImplementedError`. + + Raises: + NotImplementedError: if ``queue=True``. + NotImplementedError: if :external:attr:`~quimb.tensor.TEBD.cyclic` is ``True``. + NotImplementedError: if :external:attr:`~quimb.tensor.TEBD.imag` is ``True``. + RuntimeError: if an invalid ``direction`` is provided. + + References: + [1]: https://en.wikipedia.org/wiki/Time-evolving_block_decimation + """ + if queue: + raise NotImplementedError( # pragma: no cover + "The TEBDEvolver does not support queued operations!" + ) + + if self.cyclic: + raise NotImplementedError( # pragma: no cover + "The TEBDEvolver does not support PBCs!" + ) + + if self.imag: + raise NotImplementedError( # pragma: no cover + "The TEBDEvolver does not support imaginary time!" + ) + + # if custom dt set, scale the dt fraction + if dt is not None: + dt_frac *= dt / self._dt # pragma: no cover + + final_site_ind = self.L - 1 + if direction == "right": + for i in range(0, final_site_ind, 2): + sites = (i, (i + 1) % self.L) + gate = self._get_gate_from_ham(dt_frac, sites) + self._pt.gate_split_(gate, sites, conj=self.conjugate, **self.split_opts) + + elif direction == "left": + for i in range(1, final_site_ind, 2): + sites = (i, (i + 1) % self.L) + gate = self._get_gate_from_ham(dt_frac, sites) + self._pt.gate_split_(gate, sites, conj=self.conjugate, **self.split_opts) + + else: + # NOTE: it should not be possible to reach this but we do a sanity check to ensure that + # no faulty direction was provided to this algorithm for an unknown reason + raise RuntimeError( # pragma: no cover + f"Expected the direction to be 'left' or 'right', not {direction}!" + ) diff --git a/qiskit_addon_mpf/backends/quimb_tebd/state.py b/qiskit_addon_mpf/backends/quimb_tebd/state.py new file mode 100644 index 0000000..49f27d2 --- /dev/null +++ b/qiskit_addon_mpf/backends/quimb_tebd/state.py @@ -0,0 +1,234 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A quimb-based MPO state.""" + +from __future__ import annotations + +from functools import partialmethod +from itertools import zip_longest +from typing import Any, cast + +import numpy as np +from autoray import do +from quimb.tensor import ( + MatrixProductOperator, + MatrixProductState, + Tensor, + expec_TN_1D, +) +from quimb.tensor.tensor_core import ( + group_inds, + rand_uuid, + tensor_contract, + tensor_split, +) + +from .. import State + + +class MPOState(MatrixProductOperator, State): + """An MPO enforcing the Vidal gauge. + + This specialization of quimb's existing :external:class:`quimb.tensor.MatrixProductOperator` + enforces the Vidal gauge throughout its existence. This ensures a stable behavior of the + :class:`.DynamicMPF` algorithm when using the + :class:`~qiskit_addon_mpf.backends.quimb_tebd.TEBDEvolver`. + """ + + def __init__(self, *args, **kwargs) -> None: + """Initialize a :class:`MPOState` instance. + + .. hint:: + All arguments (positional and keyword) are simply forwarded to the + :external:class:`quimb.tensor.MatrixProductOperator` constructor. Additionally, the + :attr:`vidal_singular_values` attribute gets initialized to a list of empty lists of + length equal to the number of sites in this MPO. + + Args: + args: all positional arguments will be forwarded to the + :external:class:`quimb.tensor.MatrixProductState` constructor. + kwargs: all keyword arguments will be forwarded to the + :external:class:`quimb.tensor.MatrixProductState` constructor. + """ + super().__init__(*args, **kwargs) + + self.vidal_singular_values: list[list[float]] = [[]] * self._L + """A nested list of singular values. The outer list is of equal length as this MPO itself + (:external:attr:`quimb.tensor.TensorNetwork1D.L`). Every item is another list of all the + singular values for determining the Vidal gauge at that site. + """ + + # TODO: extend the capabilities of this function to work for gates not acting on 2 qubits. + def gate_split( + self, + gate: np.ndarray, + where: tuple[int, int], + inplace: bool = False, + conj: bool = False, + **split_opts, + ) -> MPOState: + """Apply a two-site gate and contract it back into the MPO. + + The basic principle of this method is the same as that of + :external:meth:`quimb.tensor.MatrixProductState.gate_split`. However, the implementation + ensures that the Vidal gauge is conserved. + + Args: + gate: the gate to be applied to the MPO. Its shape should be either ``(d**2, d**2)`` for + a physical dimension of ``d``, or a reshaped version thereof. + where: the indices of the sites where the gate should be applied. + inplace: whether to perform the gate application in-place or return a new + :class:`MPOState` with the gate applied to it. + conj: whether the gate should be applied to the lower (``conj=False``, the default, + :external:meth:`~quimb.tensor.tensor_arbgeom.TensorNetworkGenOperator.lower_ind`) + or upper (``conj=True``, + :external:meth:`~quimb.tensor.tensor_arbgeom.TensorNetworkGenOperator.upper_ind`) + indices of the underlying MPO. + + .. note:: + This is essentially how the LHS and RHS of the :class:`.DynamicMPF` are + differentiated, by passing their :attr:`.Evolver.conjugate` property to this + argument. + split_opts: additional keyword arguments that will be forwarded to the + :external:func:`quimb.tensor.tensor_split` function. These can be used to affect the + truncation of the tensor before it gets contracted back into the MPO. + + Returns: + The :class:`MPOState` with the ``gate`` applied and contracted. + """ + t_n = self if inplace else self.copy() + + # By default, the indices of the MPO to which we apply the gate are the "upper indices" + # (commonly labeled "k" in quimb). + inds = tuple(map(self.lower_ind, where)) + if conj: + # However, when we are dealing with the conjugated side, we apply the gate to the "lower + # indices" (commonly labeled "b" in quimb). + inds = tuple(map(self.upper_ind, where)) + # Obviously we must also conjugate the actual gate being applied. + gate = gate.conj() + + # Now that we know the indices, we can extract the left and right tensors on which this gate + # acts. + t_l, t_r = t_n._inds_get(*inds) + # And we can extract the indices of the left and right legs. The middle item of the tuple + # contains the indices of any legs shared between the tensors. In this case, we know that + # there is exactly one such index but we do not need it later and, thus, ignore it here. + left_inds, (_,), right_inds = group_inds(t_l, t_r) + + # Next, we must prepare the gate by converting it to a Tensor. + if gate.ndim != 2 * len(inds): # pragma: no branch + # The gate was supplied as a matrix, so we must factorize it. + dims = [t_n.ind_size(ix) for ix in inds] + gate = do("reshape", gate, dims * 2) + + # We must generate new indices to join the old physical sites to the new gate. + bnds = [rand_uuid() for _ in range(len(inds))] + # And we keep track of mapping inds to bnds. + reindex_map = dict(zip(inds, bnds)) + + # Now we actually create the tensor of our gate. + gate_tensor = Tensor(gate, inds=(*inds, *bnds), left_inds=bnds) + # And contract it with the left and right tensors of our MPO. + C = tensor_contract(t_l.reindex_(reindex_map), t_r.reindex_(reindex_map), gate_tensor) + # At this point we have built and contracted a network which looks as follows: + # + # | | + # --tL--tR-- + # | | + # G_gate + # | | + + # Now we create a copy of C because we want to scale the left tensor with the singular + # values ensuring our Vidal gauge. We call this new tensor `theta`: + # + # | | + # --S--tL--tR-- + # | | + # G_gate + # | | + theta = C.copy() + theta.modify( + data=np.asarray( + [ + (s or 1.0) * t + for s, t in zip_longest(self.vidal_singular_values[where[0]], theta.data) + ] + ) + ) + + # Next, we split the contracted tensor network using SVD. Since we want to obtain the + # singular values, we must set absorb=None. + # NOTE: we can ignore U because we will rely solely in V. + _, S, V = tensor_split( + theta, + left_inds=left_inds, + right_inds=right_inds, + get="tensors", + **split_opts, + absorb=None, # NOTE: this takes precedence over what might be specified in split_opts + ) + # We now store the renormalized singular values on the right index. + renorm = np.linalg.norm(S.data) + t_n.vidal_singular_values[where[1]] = S.data / renorm + + # And since we have renormalized the singular values, the same must be done for U, which we + # compute by contracting C with the conjugated V tensor. + U = tensor_contract(C, V.conj()) + U.modify(data=U.data / renorm) + + # And finally, before being able to update our left and right tensors, we must reverse the + # reindexing which was applied earlier. + rev_reindex_map = {v: k for k, v in reindex_map.items()} + t_l.reindex_(rev_reindex_map).modify(data=U.transpose_like_(t_l).data) + t_r.reindex_(rev_reindex_map).modify(data=V.transpose_like_(t_r).data) + + return t_n + + gate_split_ = partialmethod(gate_split, inplace=True) + """The ``inplace=True`` variant of :meth:`gate_split`.""" + + def overlap(self, initial_state: Any) -> complex: + """Compute the overlap of this state with the provided initial state. + + .. warning:: + This implementation only supports instances of + :external:class:`quimb.tensor.MatrixProductState` for ``initial_state``. + + Args: + initial_state: the initial state with which to compute the overlap. + + Raises: + TypeError: if the provided initial state has an incompatible type. + + Returns: + The overlap of this state with the provided one. + """ + if not isinstance(initial_state, MatrixProductState): + raise TypeError( + "MPOState.overlap is only implemented for quimb.tensor.MatrixProductState states! " + "But not for states of type '%s'", + type(initial_state), + ) + + for k in range(self._L): + C = self[k] + C.modify(data=np.sqrt(2.0) * C.data) + + ret = expec_TN_1D(initial_state, self, initial_state) + + for k in range(self._L): + C = self[k] + C.modify(data=(1.0 / np.sqrt(2.0)) * C.data) + + return cast(complex, ret) diff --git a/qiskit_addon_mpf/backends/tenpy_layers/__init__.py b/qiskit_addon_mpf/backends/tenpy_layers/__init__.py new file mode 100644 index 0000000..385a98f --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_layers/__init__.py @@ -0,0 +1,245 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A layer-wise time-evolution backend using :external:mod:`tenpy`. + +.. currentmodule:: qiskit_addon_mpf.backends.tenpy_layers + +.. caution:: + The optional dependency `TeNPy `_ was previously offered under a + GPLv3 license. + As of the release of `v1.0.4 `_ on October + 2nd, 2024, it has been offered under the Apache v2 license. + The license of this package is only compatible with Apache-licensed versions of TeNPy. + +.. warning:: + This backend is only available if the optional dependencies have been installed: + + .. code-block:: + + pip install "qiskit-addon-mpf[tenpy]" + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + :template: autosummary/class_without_inherited_members.rst + + LayerwiseEvolver + LayerModel + +Underlying method +----------------- + +This module provides a time-evolution backend similar to the TEBD-based one provided by the +:mod:`~qiskit_addon_mpf.backends.tenpy_tebd` module. The main difference is that this module gives +the user full flexibility for defining their product formulas, thereby not limiting them to the +options built into the :external:mod:`tenpy` library. + +At its core, the algorithm provided by this module is still a TEBD [1] algorithm. However, rather +than enforcing the alternating updates to the even and odd bonds of the time-evolution state (see +also :external:meth:`~tenpy.algorithms.tebd.TEBDEngine.evolve_step`) this implementation outsources the +responsibility of updating bonds in alternating fashion to the definition of multiple time-evolution +**layers**. + +This is best explained with an example. Let us assume, we have some generic Hamiltonian acting on a +1-dimensional chain of sites. + +.. hint:: + Below we are very deliberate about the order of the Hamiltonian's Pauli terms because this + directly impacts the structure of the time-evolution circuit later on. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from qiskit.quantum_info import SparsePauliOp + >>> hamil = SparsePauliOp.from_sparse_list( + ... [("ZZ", (i, i+1), 1.0) for i in range(1, 9, 2)] + + ... [("Z", (i,), 0.5) for i in range(10)] + + ... [("ZZ", (i, i+1), 1.0) for i in range(0, 9, 2)], + ... num_qubits=10, + ... ) + +Let us now inspect the time-evolution circuit of this Hamiltonian using a second-order +Suzuki-Trotter formula. + +.. plot:: + :context: close-figs + :include-source: + + >>> from qiskit.synthesis import SuzukiTrotter + >>> from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit + >>> circuit = generate_time_evolution_circuit(hamil, time=1.0, synthesis=SuzukiTrotter(order=2)) + >>> circuit.draw("mpl") # doctest: +ELLIPSIS +
+ +In the circuit above, we can clearly identify its layer-wise structure. We can emphasize this +further, by splitting the circuit into multiple layers as shown below (we recombine the ``layers`` +into a single circuit with barriers between them to ease the visualization). + +.. plot:: + :context: close-figs + :include-source: + + >>> from qiskit_addon_utils.slicing import combine_slices, slice_by_gate_types + >>> layers = slice_by_gate_types(circuit) + >>> combine_slices(layers, include_barriers=True).draw("mpl") # doctest: +ELLIPSIS +
+ +.. hint:: + The asymmetry of the central layers is a result of the implementation of Qiskit's + :external:class:`~qiskit.synthesis.SuzukiTrotter` formula. In its second-order form, it combines + the two half-time evolutions of the final term in the Hamiltonian into a single one of twice the + length. We could transpile this circuit to collapse all such subequent gates in the central two + layers (just like the last one), but for the sake of simplicity of this example, we will not do + that here. + +It is not possible to instruct TeNPy's TEBD algorithm to simulate the exact structure of the circuit +shown above. The reason for that is a limitation in its interface, as it only accepts the full +Hamiltonian to be provided which is then time-evolved using pre-defined Trotter formulas. However, +in doing so it does not treat the order of the Pauli terms in a Hamiltonian with any significance +(like we do here). + +If one wants to compute the dynamic MPF coefficients of a time-evolution employing a product formula +structure other than the ones implemented in TeNPy (like the example above), then one can use the +time-evolution algorithm provided by this module. Rather than taking a single monolithic Hamiltonian +whose time-evolution is to be modeled, the :class:`~.tenpy_layers.LayerwiseEvolver` accepts a list +of :class:`~.tenpy_layers.LayerModel` objects, each one describing an individual layer of the +product formula. This gives the user full flexibility in defining the Trotter decomposition down to +the most granular level. + +However, caution must be applied to ensure that the property of TEBD to update even and odd bonds in +an alternating manner is still guaranteed. Luckily, for quantum circuits consisting of at most +two-qubit gates, this property is satisfied by construction. + + +Code example +------------ + +In this section, we build up on the example above and show how to take a custom Trotter formula and +use it to construct a :class:`~.tenpy_layers.LayerwiseEvolver` which can be used to replace the +:class:`.tenpy_tebd.TEBDEvolver` in the workflow described in :mod:`.tenpy_tebd`. + +.. hint:: + The overall workflow of using this module is the same as of the :mod:`.tenpy_tebd` module, so be + sure to read those instructions as well. + +Simply put, we must convert each one of the circuit ``layers`` (see above) into a +:class:`~.tenpy_layers.LayerModel` instance. For this purpose, we can use its +:meth:`~.tenpy_layers.LayerModel.from_quantum_circuit` method. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from qiskit_addon_mpf.backends.tenpy_layers import LayerModel + >>> model0 = LayerModel.from_quantum_circuit(layers[0]) + >>> layer_models = [model0] + +In the code above you can see how simple the conversion is for layers which contain only two-qubit +gates acting on mutually exclusive qubits (which layers of depth 1 guarantee). + +However, we must be more careful with layers including single-qubit gates. The reason for that is +that the TEBD algorithm underlying the :class:`~.tenpy_layers.LayerwiseEvolver` must update even and +odd bonds in an alternating manner. And because single-qubit gates are not applied on a site, but +instead are split in half and applied to the bonds on either side, a layer of single-qubit gates +acting on all qubits would break this assumption. + +To circumvent this problem, we can take any layer consisting of only single-qubit gates, and apply +it twice (once on the even and once on the odd bonds) while scaling each one with a factor of +``0.5``. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> model1a = LayerModel.from_quantum_circuit(layers[1], keep_only_odd=True, scaling_factor=0.5) + >>> model1b = LayerModel.from_quantum_circuit(layers[1], keep_only_odd=False, scaling_factor=0.5) + >>> layer_models.extend([model1a, model1b]) + +Now that we know how to treat layers consisting of two-qubit and single-qubit gates, we can +transform the remaining layers. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> for layer in layers[2:]: + ... num_qubits = len(layer.data[0].qubits) + ... if num_qubits == 2: + ... layer_models.append(LayerModel.from_quantum_circuit(layer)) + ... else: + ... layer_models.append( + ... LayerModel.from_quantum_circuit(layer, keep_only_odd=True, scaling_factor=0.5) + ... ) + ... layer_models.append( + ... LayerModel.from_quantum_circuit(layer, keep_only_odd=False, scaling_factor=0.5) + ... ) + >>> assert len(layer_models) == 8 + +In the end, we have 8 :class:`~.tenpy_layers.LayerModel`'s, one for each of the 4 two-qubit layers, +and two for each of the 2 single-qubit layers. + +Finally, we can define our :class:`.ApproxEvolverFactory` protocol to be used within the +:func:`.setup_dynamic_lse` function. + +.. plot:: + :context: + :nofigs: + :include-source: + + >>> from functools import partial + >>> from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver + >>> approx_evolver_factory = partial( + ... LayerwiseEvolver, + ... layers=layer_models, + ... options={ + ... "preserve_norm": False, + ... "trunc_params": { + ... "chi_max": 10, + ... "svd_min": 1e-5, + ... "trunc_cut": None, + ... }, + ... }, + ... ) + +.. caution:: + It should be noted, that in this workflow we have not yet fixed the time step used by the Trotter + formula. We have also only set up a single repetition of the Trotter formula as the rest will be + done by the internal :class:`.DynamicMPF` algorithm, executed during :func:`.setup_dynamic_lse`. + +Of course, you could also use this to specify a :class:`.ExactEvolverFactory`. But you can also +mix-and-match a :class:`.tenpy_layers.LayerwiseEvolver` with a :class:`.tenpy_tebd.TEBDEvolver`. + + +Resources +--------- + +[1]: https://en.wikipedia.org/wiki/Time-evolving_block_decimation +""" + +# ruff: noqa: E402 +from .. import HAS_TENPY + +HAS_TENPY.require_now(__name__) + +from .evolver import LayerwiseEvolver +from .model import LayerModel + +__all__ = [ + "LayerModel", + "LayerwiseEvolver", +] diff --git a/qiskit_addon_mpf/backends/tenpy_layers/evolver.py b/qiskit_addon_mpf/backends/tenpy_layers/evolver.py new file mode 100644 index 0000000..eaa0930 --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_layers/evolver.py @@ -0,0 +1,171 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A layer-wise time-evolver using TeNPy.""" + +from __future__ import annotations + +from tenpy.linalg import TruncationError + +from .. import tenpy_tebd +from .model import LayerModel + + +class LayerwiseEvolver(tenpy_tebd.TEBDEvolver): + """A special case of the :class:`~.tenpy_tebd.TEBDEvolver` based on layer-wise evolution models. + + As also explained in :mod:`.tenpy_layers`, this implementation extracts the alternating even/odd + bond updates implemented inside of the original + :external:class:`~tenpy.algorithms.tebd.TEBDEngine` to become the end users responsibility. It + does so, by replacing the single Hamiltonian provided to the :class:`~.tenpy_tebd.TEBDEvolver` + instance with a sequence of :class:`~.tenpy_layers.LayerModel` instances. Every single instance + of these encodes a single **layer** of interactions. These should enforce the alternating + updates of even and odd bonds of the underlying tensor network. + + The motivation for this more complicated interface is that is provides a lot more flexbility and + enables users to define custom Trotter product formulas rather than being limited to the ones + implemented by TeNPy directly. + """ + + def __init__( + self, + evolution_state: tenpy_tebd.MPOState, + layers: list[LayerModel], + *args, + **kwargs, + ) -> None: + """Initialize a :class:`LayerwiseEvolver` instance. + + Args: + evolution_state: forwarded to :class:`~.tenpy_tebd.TEBDEvolver`. Please refer to its + documentation for more details. + layers: the list of models describing single layers of interactions. See above as well + as the explanations provided in :mod:`.tenpy_layers`. + args: any further positional arguments will be forwarded to the + :class:`~.tenpy_tebd.TEBDEvolver` constructor. + kwargs: any further keyword arguments will be forwarded to the + :class:`~.tenpy_tebd.TEBDEvolver` constructor. + """ + super().__init__(evolution_state, layers[0], *args, **kwargs) + self.layers = layers + """The layers of interactions used to implement the time-evolution.""" + + def evolve(self, N_steps: int, dt: float) -> TruncationError: + # pylint: disable=attribute-defined-outside-init + # NOTE: Somehow, pylint does not properly pick up on the attributes defined in the external + # base classes. + """Perform a single time step of TEBD. + + Args: + N_steps: should always be ``1`` in this case. See + :external:class:`~tenpy.algorithms.tebd.TEBDEngine` for more details. + dt: the time-step to use. + + Returns: + The truncation error. + """ + if dt is not None: # pragma: no branch + assert dt == self._U_param["delta_t"] + + trunc_err = TruncationError() + for U_idx_dt in range(len(self.layers)): + Us = self._U[U_idx_dt] + for i_bond in range(self.psi.L): + if Us[i_bond] is None: + # NOTE: in the original TeNPy implementation this handles finite vs. infinite + # boundary conditions + # NOTE: here, we leverage the same principle to automatically skip even and odd + # bond indices based on the LayerModel.calc_H_bond output + continue + self._update_index = (U_idx_dt, i_bond) + update_err = self.update_bond(i_bond, Us[i_bond]) + trunc_err += update_err + self._update_index = None # type: ignore[assignment] + self.evolved_time = self.evolved_time + N_steps * self._U_param["tau"] + self.trunc_err: TruncationError = self.trunc_err + trunc_err # not += : make a copy! + # (this is done to avoid problems of users storing self.trunc_err after each `evolve`) + return trunc_err + + @staticmethod + def suzuki_trotter_time_steps(order: int) -> list[float]: # pylint: disable=unused-argument + """Returns an empty list. + + .. note:: + This method is undefined for this subclass but we cannot raise an error upon calling + it because of the internal algorithm flow. Instead, the Trotter decomposition in this + class is encoded directly into the :attr:`layers`. + + Args: + order: is being ignored. + + Returns: + An empty list. + """ + return [] + + @staticmethod + def suzuki_trotter_decomposition( + order: int, # pylint: disable=unused-argument + N_steps: int, # pylint: disable=unused-argument + ) -> list[tuple[int, int]]: + """Returns an empty list. + + .. note:: + This method is undefined for this subclass but we cannot raise an error upon calling + it because of the internal algorithm flow. Instead, the Trotter decomposition in this + class is encoded directly into the :attr:`layers`. + + Args: + order: is being ignored. + N_steps: is being ignored. + + Returns: + An empty list. + """ + return [] # pragma: no cover + + def calc_U( + self, + order: int, + delta_t: float, + type_evo: str = "real", + E_offset: list[float] | None = None, + ) -> None: + # pylint: disable=attribute-defined-outside-init + # NOTE: Somehow, pylint does not properly pick up on the attributes defined in the external + # base classes. + """Calculates the local bond updates. + + This adapts :external:meth:`~tenpy.algorithms.tebd.TEBDEngine.calc_U` to work with the + layer-wise implementation. + + Args: + order: this is being ignored. + delta_t: the time-step to use. + type_evo: the type of time-evolution. Imaginary time-evolution is not supported at this + time. + E_offset: a constant energy offset to be applied. + """ + super().calc_U(order, delta_t, type_evo=type_evo, E_offset=E_offset) + + # NOTE: since suzuki_trotter_time_steps returns an empty list, we did not yet compute + # self._U in the super-call above and do so manually, here + L = self.psi.L + self._U = [] + for layer in self.layers: + self.model = layer + U_bond = [ + self._calc_U_bond(i_bond, delta_t, type_evo, E_offset) + for i_bond in range(L) + # NOTE: this iteration over L implies a 1D chain! + ] + self._U.append(U_bond) diff --git a/qiskit_addon_mpf/backends/tenpy_layers/model.py b/qiskit_addon_mpf/backends/tenpy_layers/model.py new file mode 100644 index 0000000..dc0e24e --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_layers/model.py @@ -0,0 +1,162 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A TeNPy-based model for describing a single layer of interactions.""" + +from __future__ import annotations + +from collections import defaultdict +from typing import cast + +from qiskit.circuit import QuantumCircuit +from tenpy.models import CouplingMPOModel, NearestNeighborModel +from tenpy.networks import Site, SpinHalfSite +from tenpy.tools.params import Config + + +class LayerModel(CouplingMPOModel, NearestNeighborModel): + """A model for representing a layer of time-evolution interactions. + + Essentially, this class is a simple wrapper of + :external:class:`tenpy.models.model.CouplingMPOModel` and + :external:class:`tenpy.models.model.NearestNeighborModel`. Its main purpose is to provide a + simple interface for constructing a TeNPy-compatible Hamiltonian from Qiskit objects. + """ + + def init_sites(self, model_params: Config) -> Site: + """Initializes the sites of this Hamiltonian. + + See :external:meth:`~tenpy.models.model.CouplingMPOModel.init_sites` for more details. + + .. caution:: + Currently, this enforces ``Sz`` conservation on all sites. + + Args: + model_params: the model parameters. + + Returns: + The site to be used internally. + """ + # WARN: we use our own default to NOT sort charges (contrary to TeNPy default: `True`) + sort_charge = model_params.get("sort_charge", False, bool) + # TODO: currently, this default Sz conservation is coupled to the LegCharge definition in + # MPOState.initialize_from_lattice. Not conserving Sz errors with 'Different ChargeInfo'. + conserve = model_params.get("conserve", "Sz", str) + return SpinHalfSite(conserve=conserve, sort_charge=sort_charge) + + def init_terms(self, model_params: Config) -> None: + """Initializes the terms of this Hamiltonian. + + See :external:meth:`~tenpy.models.model.CouplingMPOModel.init_terms` for more details. + + Args: + model_params: the model parameters. + """ + coupling_terms = model_params.get("coupling_terms", {}) + onsite_terms = model_params.get("onsite_terms", {}) + + for category, terms in coupling_terms.items(): + for term in terms: + self.add_coupling_term(*term, category=category) + + for category, terms in onsite_terms.items(): + for term in terms: + self.add_onsite_term(*term, category=category) + + def calc_H_bond(self, tol_zero: float = 1e-15) -> list: + """Calculate the interaction Hamiltonian based on the coupling and onsite terms. + + Essentially, this class overwrites + :external:meth:`~tenpy.models.model.CouplingModel.calc_H_bond` and takes care of removing + even or odd bond interaction Hamiltonians depending on the value of ``keep_only_odd`` (see + :mod:`.tenpy_layers` for more details). + + Args: + tol_zero: the threshold for values considered to be zero. + + Returns: + The list of interaction Hamiltonians for all bonds. + """ + H_bond = super().calc_H_bond(tol_zero=tol_zero) + + keep_only_odd = self.options.get("keep_only_odd", None, bool) + if keep_only_odd is None: + # return H_bond unchanged + return cast(list, H_bond) + + # NOTE: if `keep_only_odd` was specified, that means we explicitly overwrite those H_bond + # values with `None` which we do not want to keep. In the case of (for example) coupling + # layers, this should have no effect since those bonds were `None` to begin with. However, + # in the case of onsite layers, this will remove half of the bonds ensuring that we split + # the bond updates into even and odd parts (as required by the TEBD algorithm). + for i in range(0 if keep_only_odd else 1, self.lat.N_sites, 2): + H_bond[i] = None + + return cast(list, H_bond) + + @classmethod + def from_quantum_circuit( + cls, + circuit: QuantumCircuit, + *, + scaling_factor: float = 1.0, + **kwargs, + ) -> LayerModel: + """Construct a :class:`LayerModel` from a :external:class:`~qiskit.circuit.QuantumCircuit`. + + You can see an example of this function in action in the docs of :mod:`tenpy_layers`. + + Args: + circuit: the quantum circuit to parse. + scaling_factor: a factor with which to scale the term strengths. This can be used to + apply (for example) a time step scaling factor. It may also be used (e.g.) to split + onsite terms into two layers (even and odd) with $0.5$ of the strength, each. + kwargs: any additional keyword arguments to pass to the :class:`LayerModel` constructor. + + Returns: + A new LayerModel instance. + + Raises: + NotImplementedError: if an unsupported quantum gate is encountered. + """ + coupling_terms = defaultdict(list) + onsite_terms = defaultdict(list) + + for instruction in circuit.data: + op = instruction.operation + sites = [circuit.find_bit(qubit)[0] for qubit in instruction.qubits] + + # NOTE: the hard-coded scaling factors below account for the Pauli matrix conversion + if op.name == "rzz": + coupling_terms["Sz_i Sz_j"].append( + (2.0 * scaling_factor * op.params[0], *sites, "Sz", "Sz", "Id") + ) + elif op.name == "xx_plus_yy": + coupling_terms["Sp_i Sm_j"].append( + (0.5 * scaling_factor * op.params[0], *sites, "Sp", "Sm", "Id") + ) + coupling_terms["Sp_i Sm_j"].append( + (0.5 * scaling_factor * op.params[0], *sites, "Sm", "Sp", "Id") + ) + elif op.name == "rz": + onsite_terms["Sz"].append((2.0 * scaling_factor * op.params[0], *sites, "Sz")) + else: + raise NotImplementedError(f"Cannot handle gate of type {op.name}") + + return cls( + { + "L": circuit.num_qubits, + "coupling_terms": coupling_terms, + "onsite_terms": onsite_terms, + **kwargs, + } + ) diff --git a/qiskit_addon_mpf/backends/tenpy_tebd/__init__.py b/qiskit_addon_mpf/backends/tenpy_tebd/__init__.py new file mode 100644 index 0000000..571bcf7 --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_tebd/__init__.py @@ -0,0 +1,163 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A :external:mod:`tenpy`-based TEBD backend. + +.. currentmodule:: qiskit_addon_mpf.backends.tenpy_tebd + +.. caution:: + The optional dependency `TeNPy `_ was previously offered under a + GPLv3 license. + As of the release of `v1.0.4 `_ on October + 2nd, 2024, it has been offered under the Apache v2 license. + The license of this package is only compatible with Apache-licensed versions of TeNPy. + +.. warning:: + This backend is only available if the optional dependencies have been installed: + + .. code-block:: + + pip install "qiskit-addon-mpf[tenpy]" + +.. autosummary:: + :toctree: ../stubs/ + :nosignatures: + :template: autosummary/class_without_inherited_members.rst + + TEBDEvolver + MPOState + MPS_neel_state + +Underlying method +----------------- + +This module provides a time-evolution backend for computing dynamic MPF coefficients based on the +time-evolving block decimation (TEBD) algorithm [1] implemented in the :external:mod:`tenpy` tensor +network library. + +The classes provided by this module serve two purposes: + +1. Connecting :external:mod:`tenpy`'s implementation to the interface set out by + :mod:`qiskit_addon_mpf.backends`. +2. Extending :external:mod:`tenpy`'s TEBD implementation to handle an internal MPO (rather than + MPS) state (see also :class:`.State` for more details). + + +In the simplest sense, this module provides a straight-forward extension of the TEBD algorithm to +evolve an internal MPO state. +As such, if you wish to use this backend for your dynamic MPF algorithm, you must encode the +Hamiltonian that you wish to time-evolve, in a :external:mod:`tenpy`-native form. To be more +concrete, the :class:`~qiskit_addon_mpf.backends.tenpy_tebd.TEBDEvolver` class (which is a subclass +of :external:class:`tenpy.algorithms.tebd.TEBDEngine`) works with a Hamiltonian in the form of a +:external:class:`~tenpy.models.model.Model`. TeNPy provides a number of convenience methods for +constructing such Hamiltonians in its :external:mod:`tenpy.models` module. +If none of those fulfill your needs, you can consider using the +:class:`~qiskit_addon_mpf.backends.tenpy_layers.LayerModel` class which implements some conversion +methods from Qiskit-native objects. + +Code example +------------ + +This section shows a simple example to get you started with using this backend. The example shows +how to create the three factory functions required for the :func:`.setup_dynamic_lse`. + +First of all, we define the Hamiltonian which we would like to time-evolve. Here, we simply choose +one of :external:mod:`tenpy`'s convenience methods. + +>>> from tenpy.models import XXZChain2 +>>> hamil = XXZChain2( +... { +... "L": 10, +... "Jz": 0.8, +... "Jxx": 0.7, +... "hz": 0.3, +... "bc_MPS": "finite", +... "sort_charge": False, +... } +... ) + +Next, we can create the ``identity_factory`` which has to match the :class:`.IdentityStateFactory` +protocol. We do so by using the :func:`~.tenpy_tebd.MPOState.initialize_from_lattice` convenience +method which takes the lattice underlying the Hamiltonian which we just defined as its only input. + +>>> from functools import partial +>>> from qiskit_addon_mpf.backends.tenpy_tebd import MPOState +>>> identity_factory = partial(MPOState.initialize_from_lattice, hamil.lat), + +We can now construct the :class:`.ExactEvolverFactory` and :class:`.ApproxEvolverFactory` +time-evolution instance factories. To do so, we can simply bind the pre-defined values of the +:class:`~qiskit_addon_mpf.backends.tenpy_tebd.TEBDEvolver` initializer, reducing it to the correct +interface as expected by the respective function protocols. + +>>> from qiskit_addon_mpf.backends.tenpy_tebd import TEBDEvolver +>>> exact_evolver_factory = partial( +... TEBDEvolver, +... model=hamil, +... dt=0.05, +... options={ +... "order": 4, +... "preserve_norm": False, +... }, +... ) + +Notice, how we have fixed the ``dt`` value to a small time step and have used a higher-order +Suzuki-Trotter decomposition to mimic the exact time-evolution above. + +Below, we do not fix the ``dt`` value and use only a second-order Suzuki-Trotter formula for the +approximate time-evolution. Additionally, we also specify some truncation settings. + +>>> approx_evolver_factory = partial( +... TEBDEvolver, +... model=hamil, +... options={ +... "order": 2, +... "preserve_norm": False, +... "trunc_params": { +... "chi_max": 10, +... "svd_min": 1e-5, +... "trunc_cut": None, +... }, +... }, +... ) + +Of course, you are not limited to the examples shown here, and we encourage you to play around with +the other settings provided by TeNPy's :external:class:`~tenpy.algorithms.tebd.TEBDEngine` +implementation. + +Limitations +----------- + +Finally, we point out a few known limitations on what kind of Hamiltonians can be treated by this +backend: + +* all interactions must be 1-dimensional +* the interactions must use finite boundary conditions + +Resources +--------- + +[1]: https://en.wikipedia.org/wiki/Time-evolving_block_decimation +""" + +# ruff: noqa: E402 +from .. import HAS_TENPY + +HAS_TENPY.require_now(__name__) + +from .evolver import TEBDEvolver +from .state import MPOState, MPS_neel_state + +__all__ = [ + "MPS_neel_state", + "MPOState", + "TEBDEvolver", +] diff --git a/qiskit_addon_mpf/backends/tenpy_tebd/evolver.py b/qiskit_addon_mpf/backends/tenpy_tebd/evolver.py new file mode 100644 index 0000000..40fa5f1 --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_tebd/evolver.py @@ -0,0 +1,152 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A tenpy-based TEBD algorithm for evolving an internal MPO.""" + +from __future__ import annotations + +import logging +from typing import cast + +from tenpy import TEBDEngine, svd_theta +from tenpy.linalg import np_conserved as npc + +from .. import Evolver + +LOGGER = logging.getLogger(__name__) + + +class TEBDEvolver(TEBDEngine, Evolver): + """A TEBD algorithm for evolving an internal MPO. + + As discussed in more detail in :mod:`~qiskit_addon_mpf.backends.tenpy_tebd`, this extension of + TeNPy's existing :external:class:`~tenpy.algorithms.tebd.TEBDEngine` implementation time-evolves + an internal matrix product operator (MPO) rather than the conventional matrix product state + (MPS). + + More concretely, the internal object is expected to be an :class:`~.tenpy_tebd.MPOState`. + + .. warning:: + The API of this class is actually much larger than shown here, because it inherits additional + methods from the :external:class:`~tenpy.algorithms.tebd.TEBDEngine` base class. However, we + do not duplicate that API here. + """ + + def __init__(self, *args, dt: float = 0.1, **kwargs) -> None: + """Initialize a :class:`TEBDEvolver` instance. + + Args: + args: any positional arguments will be forwarded to the + :external:class:`~tenpy.algorithms.tebd.TEBDEngine` constructor. + dt: the time step to be used by this time-evolution instance. + kwargs: any further keyword arguments will be forwarded to the + :external:class:`~tenpy.algorithms.tebd.TEBDEngine` constructor. + """ + super().__init__(*args, **kwargs) + self.dt = dt + """The time step to be used by this time-evolution instance.""" + self._conjugate = False + + @property + def evolved_time(self) -> float: + """Returns the current evolution time.""" + return self._evolved_time + + @evolved_time.setter + def evolved_time(self, evolved_time: float) -> None: + self._evolved_time = evolved_time + + @property + def conjugate(self) -> bool: + """Returns whether this time-evolver instance acts on the right-hand side.""" + return self._conjugate + + @conjugate.setter + def conjugate(self, conjugate: bool) -> None: + self._conjugate = conjugate + + def step(self) -> None: + """Perform a single time step of TEBD. + + This essentially calls :external:meth:`~tenpy.algorithms.tebd.TEBDEngine.run_evolution` and + forwards the value of :attr:`dt` that was provided upon construction. + """ + self.run_evolution(1, self.dt) + + def update_bond(self, i: int, U_bond: npc.Array) -> float: + """Update the specified bond. + + Overwrites the original (MPS-based) implementation to support an MPO as the shared state. + + Args: + i: the bond index. + U_bond: the bond to update. + + Returns: + The truncation error. + """ + i0, i1 = i - 1, i + LOGGER.debug("Update sites (%d, %d)", i0, i1) + # Construct the theta matrix + C0 = self.psi.get_W(i0) + C1 = self.psi.get_W(i1) + + C = npc.tensordot(C0, C1, axes=(["wR"], ["wL"])) + new_labels = ["wL", "p0", "p0*", "p1", "p1*", "wR"] + C.iset_leg_labels(new_labels) + + if self.conjugate: + C = npc.tensordot(U_bond.conj(), C, axes=(["p0", "p1"], ["p0*", "p1*"])) # apply U + else: + C = npc.tensordot(U_bond, C, axes=(["p0*", "p1*"], ["p0", "p1"])) # apply U + C.itranspose(["wL", "p0", "p0*", "p1", "p1*", "wR"]) + theta = C.scale_axis(self.psi.Ss[i0], "wL") + # now theta is the same as if we had done + # theta = self.psi.get_theta(i0, n=2) + # theta = npc.tensordot(U_bond, theta, axes=(['p0*', 'p1*'], ['p0', 'p1'])) # apply U + # but also have C which is the same except the missing "S" on the left + # so we don't have to apply inverses of S (see below) + + # theta = theta.combine_legs([("wL", "p0", "p0*"), ("p1", "p1*", "wR")], qconj=[+1, -1]) + theta = theta.combine_legs([[0, 1, 2], [3, 4, 5]], qconj=[+1, -1]) + # Perform the SVD and truncate the wavefunction + U, S, V, trunc_err, renormalize = svd_theta( + theta, self.trunc_params, [None, None], inner_labels=["wR", "wL"] + ) + + # Split tensor and update matrices + B_R = V.split_legs(1).ireplace_labels(["p1", "p1*"], ["p", "p*"]) + + # In general, we want to do the following: + # U = U.iscale_axis(S, 'vR') + # B_L = U.split_legs(0).iscale_axis(self.psi.get_SL(i0)**-1, 'vL') + # B_L = B_L.ireplace_label('p0', 'p') + # i.e. with SL = self.psi.get_SL(i0), we have ``B_L = SL**-1 U S`` + # + # However, the inverse of SL is problematic, as it might contain very small singular + # values. Instead, we use ``C == SL**-1 theta == SL**-1 U S V``, + # such that we obtain ``B_L = SL**-1 U S = SL**-1 U S V V^dagger = C V^dagger`` + # here, C is the same as theta, but without the `S` on the very left + # (Note: this requires no inverse if the MPS is initially in 'B' canonical form) + B_L = npc.tensordot( + C.combine_legs(("p1", "p1*", "wR"), pipes=theta.legs[1]), + V.conj(), + axes=["(p1.p1*.wR)", "(p1*.p1.wR*)"], + ) + B_L.ireplace_labels(["wL*", "p0", "p0*"], ["wR", "p", "p*"]) + B_L /= renormalize # re-normalize to = 1 + + self.psi.Ss[i1] = S + self.psi.set_W(i0, B_L) + self.psi.set_W(i1, B_R) + self._trunc_err_bonds[i] = self._trunc_err_bonds[i] + trunc_err + return cast(float, trunc_err) diff --git a/qiskit_addon_mpf/backends/tenpy_tebd/state.py b/qiskit_addon_mpf/backends/tenpy_tebd/state.py new file mode 100644 index 0000000..16a02be --- /dev/null +++ b/qiskit_addon_mpf/backends/tenpy_tebd/state.py @@ -0,0 +1,136 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""A TeNPy-based MPO state.""" + +from __future__ import annotations + +from typing import Any, cast + +import numpy as np +from tenpy.linalg import np_conserved as npc +from tenpy.models import Lattice +from tenpy.networks import MPO, MPS + +from .. import State + + +class MPOState(MPO, State): + """A mediator class to make TeNPy's MPO match the :class:`.State` interface. + + This class simply ensures that a :external:class:`tenpy.networks.mpo.MPO` object can work as a + :class:`.State` instance. + """ + + @classmethod + def initialize_from_lattice(cls, lat: Lattice) -> MPOState: + """Construct an identity :class:`MPOState` instance matching the provided lattice shape. + + Given a lattice, this method constructs a new MPO identity matching the shape of the + lattice. + + Args: + lat: the lattice describing the MPO sites. + + Returns: + An identity MPO. + """ + # creates a 4d array filled with zeros - shape 1x2x2x1 + B = np.zeros([1, 2, 2, 1], dtype=float) + # sets element values of B array to 1 + # creates a tensor that represents identity for MPO + B[0, 0, 0, 0] = 1 + B[0, 1, 1, 0] = 1 + + labels = ["wL", "p", "p*", "wR"] + + # creates a list of tensor leg charge objects encoding charges + conjugations for tensor + # legs (i.e. dimensions) + leg_charge = [ + # e.g. charge information for tensor leg / dimension [1] and label ["2*Sz"] + # creates a LegCharge object from the flattened list of charges + # one for each of four legs or dimensions on B + npc.LegCharge.from_qflat(npc.ChargeInfo([1], ["2*Sz"]), [1], qconj=1), + npc.LegCharge.from_qflat(npc.ChargeInfo([1], ["2*Sz"]), [1, -1], qconj=1), + npc.LegCharge.from_qflat(npc.ChargeInfo([1], ["2*Sz"]), [1, -1], qconj=-1), + npc.LegCharge.from_qflat(npc.ChargeInfo([1], ["2*Sz"]), [1], qconj=-1), + ] + + B_array = npc.Array.from_ndarray(B, legcharges=leg_charge, labels=labels) + + # FIXME: the following is supposed to allow `conserve=None` for our `SpinHalfSite` but it + # does not seem to work, even in the `conserve="Sz"` case. + # + # B_array = npc.Array.from_ndarray_trivial(B, labels=labels) + + num_sites = lat.N_sites + # initialize the MPO psi with the wavepacket and an identity operator + psi = cls.from_wavepacket(lat.mps_sites(), [1.0] * num_sites, "Id") + + # set the wavefunction at each site in psi to B_array + for k in range(num_sites): + psi.set_W(k, B_array) + + # srt the bond strengths of psi to a list of lists with all elements as 1.0 + psi.Ss = [[1.0]] * num_sites + # psi is now an MPO representing the identity operator + # psi consists of an identical B_array at each site + # psi is a product of local identity operators since the bond dimensions are all 1 + return cast(MPOState, psi) + + def overlap(self, initial_state: Any) -> complex: + """Compute the overlap of this state with the provided initial state. + + .. warning:: + This implementation only supports instances of + :external:class:`tenpy.networks.mps.MPS` for ``initial_state``. + + Args: + initial_state: the initial state with which to compute the overlap. + + Raises: + TypeError: if the provided initial state has an incompatible type. + + Returns: + The overlap of this state with the provided one. + """ + if not isinstance(initial_state, MPS): + raise TypeError( + "MPOState.overlap is only implemented for tenpy.networks.mps.MPS states! " + "But not for states of type '%s'", + type(initial_state), + ) + + for k in range(self.L): + self.set_W(k, np.sqrt(2.0) * self.get_W(k)) + + overlap = self.expectation_value(initial_state) + + for k in range(self.L): + self.set_W(k, (1.0 / np.sqrt(2.0)) * self.get_W(k)) + + return cast(complex, overlap) + + +def MPS_neel_state(lat: Lattice) -> MPS: + """Constructs the Néel state as an MPS. + + Args: + lat: the lattice describing the MPS sites. + + Returns: + A Néel state as an MPS. + """ + num_sites = lat.N_sites + product_state = ["up", "down"] * (num_sites // 2) + (num_sites % 2) * ["up"] + initial_state = MPS.from_product_state(lat.mps_sites(), product_state, bc=lat.bc_MPS) + return initial_state diff --git a/qiskit_addon_mpf/costs/__init__.py b/qiskit_addon_mpf/costs/__init__.py new file mode 100644 index 0000000..4d5ab9b --- /dev/null +++ b/qiskit_addon_mpf/costs/__init__.py @@ -0,0 +1,44 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Cost functions for MPF coefficients. + +.. currentmodule:: qiskit_addon_mpf.costs + +This module provides a number of optimization problem generator functions, each implementing a +different cost function as the problem's target objective. All of the functions provided by this +module take a linear system of equations (:class:`.LSE`) encoding the parameters of the optimization +problem as their first argument. + +.. autoclass:: LSE + +Optimization problem constructors +--------------------------------- + +.. autofunction:: setup_exact_problem + +.. autofunction:: setup_sum_of_squares_problem + +.. autofunction:: setup_frobenius_problem +""" + +from .exact import setup_exact_problem +from .frobenius import setup_frobenius_problem +from .lse import LSE +from .sum_of_squares import setup_sum_of_squares_problem + +__all__ = [ + "LSE", + "setup_exact_problem", + "setup_sum_of_squares_problem", + "setup_frobenius_problem", +] diff --git a/qiskit_addon_mpf/static/exact.py b/qiskit_addon_mpf/costs/exact.py similarity index 74% rename from qiskit_addon_mpf/static/exact.py rename to qiskit_addon_mpf/costs/exact.py index 126e64d..f2aa305 100644 --- a/qiskit_addon_mpf/static/exact.py +++ b/qiskit_addon_mpf/costs/exact.py @@ -10,7 +10,7 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Exact static MPF coefficients.""" +"""Exact MPF coefficients.""" from __future__ import annotations @@ -19,17 +19,17 @@ from .lse import LSE -def setup_exact_model(lse: LSE) -> tuple[cp.Problem, cp.Variable]: - r"""Construct a :external:class:`cvxpy.Problem` for finding exact static MPF coefficients. +def setup_exact_problem(lse: LSE) -> tuple[cp.Problem, cp.Variable]: + r"""Construct a :external:class:`cvxpy.Problem` for finding the exact MPF coefficients. .. note:: The coefficients found via this optimization problem will be identical to the analytical ones obtained from the :meth:`.LSE.solve` method. This additional interface exists to highlight - the parallel to the :func:`.setup_approximate_model` interface. It also serves educational - purposes for how-to approach optimization problems targeting MPF coefficients. + the parallel to the other cost functions provided by this module. It also serves educational + purposes for how to approach optimization problems targeting MPF coefficients. - The optimization problem constructed by this class is defined as follows: + The optimization problem constructed by this function is defined as follows: - the cost function minimizes the L1-norm (:external:class:`~cvxpy.atoms.norm1.norm1`) of the variables (:attr:`.LSE.x`) @@ -41,16 +41,17 @@ def setup_exact_model(lse: LSE) -> tuple[cp.Problem, cp.Variable]: Here is an example: - >>> from qiskit_addon_mpf.static import setup_lse, setup_exact_model - >>> lse = setup_lse([1,2,3], order=2, symmetric=True) - >>> problem, coeffs = setup_exact_model(lse) + >>> from qiskit_addon_mpf.costs import setup_exact_problem + >>> from qiskit_addon_mpf.static import setup_static_lse + >>> lse = setup_static_lse([1,2,3], order=2, symmetric=True) + >>> problem, coeffs = setup_exact_problem(lse) >>> print(problem) minimize norm1(x) subject to Sum([1. 1. 1.] @ x, None, False) == 1.0 Sum([1. 0.25 0.11111111] @ x, None, False) == 0.0 Sum([1. 0.0625 0.01234568] @ x, None, False) == 0.0 - You can then solve the problem and extract the expansion coefficients like so: + You can then solve the problem and access the expansion coefficients like so: >>> final_cost = problem.solve() >>> print(coeffs.value) # doctest: +FLOAT_CMP diff --git a/qiskit_addon_mpf/costs/frobenius.py b/qiskit_addon_mpf/costs/frobenius.py new file mode 100644 index 0000000..4fe23bb --- /dev/null +++ b/qiskit_addon_mpf/costs/frobenius.py @@ -0,0 +1,126 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Frobenius norm MPF coefficients.""" + +from __future__ import annotations + +import cvxpy as cp + +from .lse import LSE + + +def setup_frobenius_problem( + lse: LSE, *, max_l1_norm: float = 10.0 +) -> tuple[cp.Problem, cp.Variable]: + r"""Construct a :external:class:`cvxpy.Problem` for finding approximate MPF coefficients. + + The optimization problem constructed by this function is defined as follows: + + - the cost function minimizes the following quadratic expression: + + .. math:: + 1 + x^T A x - 2 x^T b + + As shown in [1] and [2], this expression arises from the Frobenius norm of the error between + an exact time evolution state and a dynamic MPF. As such, taking the :class:`.LSE` constructed + by :func:`.setup_dynamic_lse` and plugging it into this function will yield Eq. (20) of [1] + (which is identical to Eq. (2) of [2]), which we repeat below + + .. math:: + 1 + \sum_{i,j} A_{ij}(t) x_i(t) x_j(t) - 2 \sum_i b_i(t) x_i(t) \, , + + where $A$ and $b$ of our :class:`.LSE` correspond to the Gram matrix ($M$ in [1] and [2]) and + the overlap vector ($L$ in [1] and [2]), respectively. Additionally, we use $x(t)$ to denote + the MPF variables (or coefficients) rather than $c$ in [1] and [2]. + + - two constraints are set: + + 1. the variables must sum to 1: :math:`\sum_i x_i == 1` + 2. the L1-norm (:external:class:`~cvxpy.atoms.norm1.norm1`) of the variables is bounded by + ``max_l1_norm`` + + Below is an example which uses the ``lse`` object constructed in the example for + :func:`.setup_dynamic_lse`. + + .. testsetup:: + >>> from functools import partial + >>> from qiskit_addon_mpf.backends.quimb_tebd import MPOState, TEBDEvolver + >>> from qiskit_addon_mpf.dynamic import setup_dynamic_lse + >>> from quimb.tensor import ham_1d_heis, MPO_identity, MPS_neel_state + >>> trotter_steps = [3, 4] + >>> time = 0.9 + >>> num_qubits = 10 + >>> initial_state = MPS_neel_state(num_qubits) + >>> hamil = ham_1d_heis(num_qubits, 0.8, 0.3, cyclic=False) + >>> identity_factory = lambda: MPOState(MPO_identity(num_qubits)) + >>> exact_evolver_factory = partial( + ... TEBDEvolver, + ... H=hamil, + ... dt=0.05, + ... order=4, + ... split_opts={"max_bond": 10, "cutoff": 1e-5}, + ... ) + >>> approx_evolver_factory = partial( + ... TEBDEvolver, + ... H=hamil, + ... order=2, + ... split_opts={"max_bond": 10, "cutoff": 1e-5}, + ... ) + >>> lse = setup_dynamic_lse( + ... trotter_steps, + ... time, + ... identity_factory, + ... exact_evolver_factory, + ... approx_evolver_factory, + ... initial_state, + ... ) + + .. doctest:: + >>> from qiskit_addon_mpf.costs import setup_frobenius_problem + >>> problem, coeffs = setup_frobenius_problem(lse, max_l1_norm=3.0) + >>> print(problem) # doctest: +FLOAT_CMP + minimize 1.0 + QuadForm(x, [[1.00 1.00] + [1.00 1.00]]) + -[2.00003171 1.99997911] @ x + subject to Sum(x, None, False) == 1.0 + norm1(x) <= 3.0 + + You can then solve the problem and access the expansion coefficients like so: + + .. testsetup:: + >>> import sys, pytest + >>> if not sys.platform.startswith("linux"): + ... pytest.skip("This doctest only converges to numerically identical values on Linux") + + .. doctest:: + >>> final_cost = problem.solve() + >>> print(coeffs.value) # doctest: +FLOAT_CMP + [0.50596416 0.49403584] + + Args: + lse: the linear system of equations from which to build the model. + max_l1_norm: the upper limit to use for the constrain of the L1-norm of the variables. + + Returns: + The optimization problem and coefficients variable. + + References: + [1]: S. Zhuk et al., Phys. Rev. Research 6, 033309 (2024). + https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309 + [2]: N. Robertson et al., arXiv:2407.17405v2 (2024). + https://arxiv.org/abs/2407.17405v2 + """ + coeffs = lse.x + cost = 1.0 + cp.quad_form(coeffs, lse.A) - 2.0 * lse.b.T @ coeffs + constraints = [cp.sum(coeffs) == 1, cp.norm1(coeffs) <= max_l1_norm] + problem = cp.Problem(cp.Minimize(cost), constraints) + return problem, coeffs diff --git a/qiskit_addon_mpf/costs/lse.py b/qiskit_addon_mpf/costs/lse.py new file mode 100644 index 0000000..0a3bc50 --- /dev/null +++ b/qiskit_addon_mpf/costs/lse.py @@ -0,0 +1,73 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Linear system of equations supporting cvxpy variable and parameter objects.""" + +from __future__ import annotations + +from typing import NamedTuple, cast + +import cvxpy as cp +import numpy as np + + +class LSE(NamedTuple): + """A :class:`.namedtuple` representing a linear system of equations. + + .. math:: + A x = b + """ + + A: np.ndarray + """The left hand side of the LSE.""" + + b: np.ndarray + """The right hand side of the LSE.""" + + @property + def x(self) -> cp.Variable: + """Returns the $x$ :external:class:`~cvxpy.expressions.variable.Variable`.""" + return cp.Variable(shape=len(self.b), name="x") + + def solve(self) -> np.ndarray: + r"""Return the solution to this LSE: :math:`x=A^{-1}b`. + + Returns: + The solution to this LSE. + + Raises: + ValueError: if this LSE is parameterized with unassigned values. + ValueError: if this LSE does not include a row ensuring that :math:`\sum_i x_i == 1` + which is a requirement for valid MPF coefficients. + """ + if self.A.ndim == 1: + # self.A is a vector of cp.Expression objects + mat_a = np.array([row.value for row in self.A]) + if any(row is None for row in mat_a): + raise ValueError( + "This LSE contains unassigned parameter values! Assign a value to them first " + "before trying to solve this LSE again." + ) + else: + mat_a = self.A + + vec_b = self.b + ones = [all(row == 1) for row in mat_a] + if not any(ones) or not np.isclose(vec_b[np.where(ones)], 1.0): + raise ValueError( + "This LSE does not enforce the sum of all coefficients to be equal to 1 which is " + "required for valid MPF coefficients. To find valid coefficients for this LSE use " + "one of the non-exact cost functions provided in this module and find its optimal " + "solution." + ) + + return cast(np.ndarray, np.linalg.solve(mat_a, vec_b)) diff --git a/qiskit_addon_mpf/static/approximate.py b/qiskit_addon_mpf/costs/sum_of_squares.py similarity index 75% rename from qiskit_addon_mpf/static/approximate.py rename to qiskit_addon_mpf/costs/sum_of_squares.py index 44a58fa..41ab145 100644 --- a/qiskit_addon_mpf/static/approximate.py +++ b/qiskit_addon_mpf/costs/sum_of_squares.py @@ -10,7 +10,7 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Approximate static MPF coefficients.""" +"""Sum-of-squares MPF coefficients.""" from __future__ import annotations @@ -19,12 +19,12 @@ from .lse import LSE -def setup_approximate_model( +def setup_sum_of_squares_problem( lse: LSE, *, max_l1_norm: float = 10.0 ) -> tuple[cp.Problem, cp.Variable]: - r"""Construct a :external:class:`cvxpy.Problem` for finding approximate static MPF coefficients. + r"""Construct a :external:class:`cvxpy.Problem` for finding approximate MPF coefficients. - The optimization problem constructed by this class is defined as follows: + The optimization problem constructed by this function is defined as follows: - the cost function minimizes the sum of squares (:external:func:`~cvxpy.atoms.sum_squares.sum_squares`) of the distances to an exact solution @@ -41,17 +41,18 @@ def setup_approximate_model( Here is an example: - >>> from qiskit_addon_mpf.static import setup_lse, setup_approximate_model - >>> lse = setup_lse([1,2,3], order=2, symmetric=True) - >>> problem, coeffs = setup_approximate_model(lse, max_l1_norm=3.0) - >>> print(problem) + >>> from qiskit_addon_mpf.costs import setup_sum_of_squares_problem + >>> from qiskit_addon_mpf.static import setup_static_lse + >>> lse = setup_static_lse([1,2,3], order=2, symmetric=True) + >>> problem, coeffs = setup_sum_of_squares_problem(lse, max_l1_norm=3.0) + >>> print(problem) # doctest: +FLOAT_CMP minimize quad_over_lin(Vstack([1. 1. 1.] @ x + -1.0, [1. 0.25 0.11111111] @ x + -0.0, [1. 0.0625 0.01234568] @ x + -0.0), 1.0) subject to Sum(x, None, False) == 1.0 norm1(x) <= 3.0 - You can then solve the problem and extract the expansion coefficients like so: + You can then solve the problem and access the expansion coefficients like so: >>> final_cost = problem.solve() >>> print(coeffs.value) # doctest: +FLOAT_CMP @@ -65,8 +66,8 @@ def setup_approximate_model( The optimization problem and coefficients variable. References: - [1]: S. Zhuk et al., arXiv:2306.12569 (2023). - https://arxiv.org/abs/2306.12569 + [1]: S. Zhuk et al., Phys. Rev. Research 6, 033309 (2024). + https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309 """ coeffs = lse.x # NOTE: the following list comprehension is required to support parameterized LSE objects diff --git a/qiskit_addon_mpf/dynamic.py b/qiskit_addon_mpf/dynamic.py new file mode 100644 index 0000000..f62e2ae --- /dev/null +++ b/qiskit_addon_mpf/dynamic.py @@ -0,0 +1,372 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Dynamic MPF coefficients. + +.. currentmodule:: qiskit_addon_mpf.dynamic + +This module provides the generator function for the linear system of equations (:class:`.LSE`) for +computing dynamic (that is, time-dependent) MPF coefficients. + +.. autofunction:: setup_dynamic_lse + +Factory Protocols +----------------- + +The following protocols define the function signatures for the various object factory arguments. + +.. autoclass:: IdentityStateFactory + +.. autoclass:: ExactEvolverFactory + +.. autoclass:: ApproxEvolverFactory + +Core algorithm +-------------- + +.. autoclass:: DynamicMPF +""" + +from __future__ import annotations + +import logging +from typing import Any, Protocol + +import numpy as np + +from .backends.interface import Evolver, State +from .costs import LSE + +LOGGER = logging.getLogger(__name__) + + +class DynamicMPF: + """The dynamic MPF algorithm. + + Instantiated with a LHS and RHS :class:`.Evolver` this algorithm will + :meth:`~.DynamicMPF.evolve` a shared :class:`.State` up to a target evolution time. + Afterwards, the :meth:`.DynamicMPF.overlap` of the time-evolved :class:`.State` with some + initial state can be computed. See :func:`.setup_dynamic_lse` for a more detailed explanation on + how this is used to compute the elements :math:`M_{ij}` and :math:`L_i` making up the + :class:`.LSE` of the dynamic MPF coefficients. + + References: + [1]: S. Zhuk et al., Phys. Rev. Research 6, 033309 (2024). + https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309 + [2]: N. Robertson et al., arXiv:2407.17405 (2024). + https://arxiv.org/abs/2307.17405 + """ + + def __init__(self, evolution_state: State, lhs: Evolver, rhs: Evolver) -> None: + """Construct a :class:`.DynamicMPF` instance. + + Args: + evolution_state: the state to be shared by the LHS and RHS time-evolution engines. + lhs: the LHS time-evolution engine. + rhs: the RHS time-evolution engine. + """ + self.evolution_state = evolution_state + """The state shared between the LHS and RHS time-evolution engines.""" + self.lhs = lhs + """The LHS time-evolution engine.""" + self.rhs = rhs + """The RHS time-evolution engine.""" + + def evolve(self, time: float) -> None: + """Evolve the dynamic MPF algorithm up to the provided time. + + This actually runs the dynamic MPF algorithm by time-evolving + :attr:`.DynamicMPF.evolution_state` up to the specified time using the LHS and RHS + :class:`.Evolver` instances. + + Args: + time: the total target evolution time. + + Raises: + RuntimeError: if the LHS and RHS evolved times are not equal at the end. + """ + while np.round(self.lhs.evolved_time, 8) < time: + while self.rhs.evolved_time < self.lhs.evolved_time: + self.rhs.step() + LOGGER.info("Stepped RHS to %s", self.rhs.evolved_time) + self.lhs.step() + LOGGER.info("Stepped LHS to %s", self.lhs.evolved_time) + # we must ensure that the RHS can catch up with the LHS + while np.round(self.rhs.evolved_time, 8) < time: + self.rhs.step() + LOGGER.info("Stepped RHS to %s", self.rhs.evolved_time) + + LOGGER.info("Final times are %s and %s", self.lhs.evolved_time, self.rhs.evolved_time) + + if not np.isclose(self.lhs.evolved_time, self.rhs.evolved_time): + # NOTE: it should not be possible to reach this but we do a sanity check nonetheless to + # ensure that the user knows if something went really wrong + raise RuntimeError( # pragma: no cover + "The evolved times of the LHS and RHS are not equal: ", + self.lhs.evolved_time, + self.rhs.evolved_time, + ) + + def overlap(self, initial_state: Any) -> complex: + """Compute the overlap of :attr:`.DynamicMPF.evolution_state` with the provided state. + + .. warning:: + The type of the provided ``initial_state`` will depend on the chosen backend used for the + :class:`.State` and :class:`.Evolver` instances provided to this :class:`.DynamicMPF` + instance. In other words, a backend may only support a specific type of ``initial_state`` + objects for this overlap computation. See also the explanations of the ``initial_state`` + argument to the :func:`.setup_dynamic_lse` for more details. + + Args: + initial_state: the initial state with which to compute the overlap. + + Raises: + TypeError: if the provided initial state has an incompatible type. + + Returns: + The overlap of :attr:`.DynamicMPF.evolution_state` with the provided one. + """ + return self.evolution_state.overlap(initial_state) + + +class IdentityStateFactory(Protocol): + r"""The factory function protocol for constructing an identity :class:`.State` instance. + + As explained in more detail in :func:`.setup_dynamic_lse`, this factory function is called to + initialize the :attr:`.DynamicMPF.evolution_state` with an identity or empty state. This + function should not take any arguments and return a :class:`.State` instance. + """ + + def __call__(self) -> State: # noqa: D102 + ... # pragma: no cover + + +class ExactEvolverFactory(Protocol): + r"""The factory function protocol for constructing an exact :class:`.Evolver` instance. + + As explained in more detail in :func:`.setup_dynamic_lse`, this factory function is called to + initialize the :attr:`.DynamicMPF.lhs` instances of :class:`.Evolver` which produce the exact + time-evolution state, :math:`\rho(t)`, when computing the elements :math:`L_i`. + """ + + def __call__(self, evolution_state: State, /) -> Evolver: # noqa: D102 + ... # pragma: no cover + + +class ApproxEvolverFactory(Protocol): + r"""The factory function protocol for constructing an approximate :class:`.Evolver` instance. + + As explained in more detail in :func:`.setup_dynamic_lse`, this factory function is called to + initialize either the :attr:`.DynamicMPF.rhs` instances of :class:`.Evolver` when computing the + elements :math:`L_i` or both sides (:attr:`.DynamicMPF.lhs` and :attr:`.DynamicMPF.rhs`) when + computing elements :math:`M_{ij}`. Since these approximate time evolution states depend on the + Trotter step (:math:`\rho_{k_i}(t)`), this function requires the time step of the time evolution + to be provided as a keyword argument called ``dt``. + """ + + def __call__(self, evolution_state: State, /, *, dt: float = 1.0) -> Evolver: # noqa: D102 + ... # pragma: no cover + + +def setup_dynamic_lse( + trotter_steps: list[int], + time: float, + identity_factory: IdentityStateFactory, + exact_evolver_factory: ExactEvolverFactory, + approx_evolver_factory: ApproxEvolverFactory, + initial_state: Any, +) -> LSE: + r"""Return the linear system of equations for computing dynamic MPF coefficients. + + This function uses the :class:`.DynamicMPF` algorithm to compute the components of the Gram + matrix (:attr:`.LSE.A`, :math:`M` in [1] and [2]) and the overlap vector (:attr:`.LSE.b`, + :math:`L` in [1] and [2]) for the provided time-evolution parameters. + + The elements of the Gram matrix, :math:`M_{ij}`, and overlap vector, :math:`L_i`, are defined as + + .. math:: + M_{ij} &= \text{Tr}(\rho_{k_i}(t)\rho_{k_j}(t)) \, , \\ + L_i &= \text{Tr}(\rho(t)\rho_{k_i}(t)) \, , + + where :math:`\rho(t)` is the exact time-evolution state at time :math:`t` and + :math:`\rho_{k_i}(t)` is the time-evolution state approximated using :math:`k_i` Trotter steps. + + Computing the dynamic (that is, time-dependent) MPF coefficients from :math:`M` and :math:`L` + amounts to finding a solution to the :class:`.LSE` (similarly to how the :mod:`.static` MPF + coefficients are computed) while enforcing the constraint that all coefficients must sum to 1 + (:math:`\sum_i x_i = 1`), which is not enforced as part of this LSE (unlike in the static case). + Optimization problems which include this additional constraint are documented in the + :mod:`.costs` module. The one suggested by [1] and [2] is the + :meth:`~qiskit_addon_mpf.costs.setup_frobenius_problem`. + + Evaluating every element :math:`M_{ij}` and :math:`L_i` requires computing the overlap between + two time-evolution states. The :class:`.DynamicMPF` algorithm does so by means of tensor network + calculations, provided by one of the optional dependencies. The available backends are listed + and explained in more detail in the :mod:`.backends` module. + + Below, we provide an example using the :mod:`~.qiskit_addon_mpf.backends.quimb_tebd` backend. + We briefly explain each element. + + First, we initialize a simple Heisenberg Hamiltonian which we would like to time-evolve. Since + we are using a time-evolver based on :external:mod:`quimb`, we also initialize the Hamiltonian + using that library. + + >>> from quimb.tensor import ham_1d_heis + >>> num_qubits = 10 + >>> hamil = ham_1d_heis(num_qubits, 0.8, 0.3, cyclic=False) + + Next, we define the number of Trotter steps to make up our MPF, the target evolution time as + well as the initial state (:math:`\psi_{in}` in [1] and :math:`\psi_0` in [2], resp.) with + respect to which we compute the overlap between the time-evolution states. Here, we simply use + the Néel state which we also construct using :external:mod:`quimb`: + + >>> trotter_steps = [3, 4] + >>> time = 0.9 + + >>> from quimb.tensor import MPS_neel_state + >>> initial_state = MPS_neel_state(num_qubits) + + Since we must run the full :class:`.DynamicMPF` algorithm for computing every element of + :math:`M_{ij}` and :math:`L_i`, we must provide factory methods for initializing the input + arguments of the :class:`.DynamicMPF` instances. To this end, we must provide three functions. + To construct these, we will use the :func:`functools.partial` function. + + >>> from functools import partial + + First, we need a function to initialize an empty time-evolution state (see also + :attr:`.DynamicMPF.evolution_state` for more details). This constructor function may not take + any positional or keyword arguments and must return a :class:`.State` object. + + >>> from qiskit_addon_mpf.backends.quimb_tebd import MPOState + >>> from quimb.tensor import MPO_identity + >>> identity_factory = lambda: MPOState(MPO_identity(num_qubits)) + + The second and third function must construct the left- and right-hand side time-evolution + engines (see also :attr:`.DynamicMPF.lhs` and :attr:`.DynamicMPF.rhs` for more details). These + functions should follow the :class:`.ExactEvolverFactory` and :class:`.ApproxEvolverFactory` + protocols, respectively. + + The :class:`.ExactEvolverFactory` function should take a :class:`.State` object as its only + positional argument and should return a :class:`.Evolver` object, which will be used for + computing the LHS of the :math:`L_i` elements (i.e. it should produce the exact time-evolution + state, :math:`\rho(t)`). + + Here, we approximate the exact time-evolved state with a fourth-order Suzuki-Trotter formula + using a small time step of 0.05. We also specify some :external:mod:`quimb`-specific truncation + options to bound the maximum bond dimension of the underlying tensor network as well as the + minimum singular values of the split tensor network bonds. + + >>> from qiskit_addon_mpf.backends.quimb_tebd import TEBDEvolver + >>> exact_evolver_factory = partial( + ... TEBDEvolver, + ... H=hamil, + ... dt=0.05, + ... order=4, + ... split_opts={"max_bond": 10, "cutoff": 1e-5}, + ... ) + + The :class:`.ApproxEvolverFactory` function should also take a :class:`.State` object as its + only positional argument and additionally a keyword argument called ``dt`` to specify the time + step of the time-evolution. It should also return a :class:`.Evolver` object which produces the + approximate time-evolution states, :math:`\rho_{k_i}(t)`, where :math:`k_i` is determined by the + chosen time step, ``dt``. As such, these instances will be used for computing the RHS of the + :math:`L_i` as well as both sides of the :math:`M_{ij}` elements. + + Here, we use a second-order Suzuki-Trotter formula with the same truncation settings as before. + + >>> approx_evolver_factory = partial( + ... TEBDEvolver, + ... H=hamil, + ... order=2, + ... split_opts={"max_bond": 10, "cutoff": 1e-5}, + ... ) + + Finally, we can initialize and run the :func:`.setup_dynamic_lse` function to obtain the + :class:`.LSE` described at the top. + + >>> from qiskit_addon_mpf.dynamic import setup_dynamic_lse + >>> lse = setup_dynamic_lse( + ... trotter_steps, + ... time, + ... identity_factory, + ... exact_evolver_factory, + ... approx_evolver_factory, + ... initial_state, + ... ) + >>> print(lse.A) # doctest: +FLOAT_CMP + [[1. 0.99998513] + [0.99998513 1. ]] + >>> print(lse.b) # doctest: +FLOAT_CMP + [1.00001585 0.99998955] + + Args: + trotter_steps: the sequence of trotter steps to be used. + time: the total target evolution time. + identity_factory: a function to generate an empty :class:`.State` object. + exact_evolver_factory: a function to initialize the :class:`.Evolver` instance which + produces the exact time-evolution state, :math:`\rho(t)`. + approx_evolver_factory: a function to initialize the :class:`.Evolver` instance which + produces the approximate time-evolution state, :math:`\rho_{k_i}(t)`, for different + values of :math:`k_i` depending on the provided time step, ``dt``. + initial_state: the initial state (:math:`\psi_{in}` or :math:`\psi_0`) with respect to which + to compute the elements :math:`M_{ij}` of :attr:`.LSE.A` and :math:`L_i` of + :attr:`.LSE.b`. The type of this object must match the tensor network backend chosen for + the previous arguments. + + Returns: + The :class:`.LSE` to find the dynamic MPF coefficients as described above. + + References: + [1]: S. Zhuk et al., Phys. Rev. Research 6, 033309 (2024). + https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309 + [2]: N. Robertson et al., arXiv:2407.17405v2 (2024). + https://arxiv.org/abs/2407.17405v2 + """ + Bs = np.zeros(len(trotter_steps), dtype=np.complex128) + for idx, k in enumerate(trotter_steps): + dt = time / k + + common_state = identity_factory() + lhs = exact_evolver_factory(common_state) + lhs.conjugate = True + rhs = approx_evolver_factory(common_state, dt=dt) + + algo = DynamicMPF(common_state, lhs, rhs) + algo.evolve(time) + + overlap = algo.overlap(initial_state) + + Bs[idx] = overlap + + Bs = np.power(np.abs(Bs), 2) + + As = np.eye(len(trotter_steps), dtype=np.complex128) + for idx1, idx2 in zip(*np.triu_indices_from(As, k=1)): + dta = time / trotter_steps[idx1] + dtb = time / trotter_steps[idx2] + + common_state = identity_factory() + lhs = approx_evolver_factory(common_state, dt=dta) + lhs.conjugate = True + rhs = approx_evolver_factory(common_state, dt=dtb) + + algo = DynamicMPF(common_state, lhs, rhs) + algo.evolve(time) + + overlap = algo.overlap(initial_state) + + As[idx1][idx2] = overlap + As[idx2][idx1] = overlap + + As = np.power(np.abs(As), 2) + + return LSE(As, Bs) diff --git a/qiskit_addon_mpf/static/lse.py b/qiskit_addon_mpf/static.py similarity index 63% rename from qiskit_addon_mpf/static/lse.py rename to qiskit_addon_mpf/static.py index 8ec0a37..186bff6 100644 --- a/qiskit_addon_mpf/static/lse.py +++ b/qiskit_addon_mpf/static.py @@ -10,63 +10,31 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Linear system of equations for static MPF coefficients.""" +"""Static MPF coefficients. -from __future__ import annotations - -from typing import NamedTuple, cast - -import cvxpy as cp -import numpy as np - - -class LSE(NamedTuple): - """A :class:`.namedtuple` representing a linear system of equations. - - .. math:: - A x = b - """ +.. currentmodule:: qiskit_addon_mpf.static - A: np.ndarray - """The left hand side of the LSE.""" +This module provides the generator function for the linear system of equations (:class:`.LSE`) for +computing static (that is, time-independent) MPF coefficients. - b: np.ndarray - """The right hand side of the LSE.""" +.. autofunction:: setup_static_lse +""" - @property - def x(self) -> cp.Variable: - """Returns the $x$ :external:class:`~cvxpy.expressions.variable.Variable`.""" - return cp.Variable(shape=len(self.b), name="x") - - def solve(self) -> np.ndarray: - """Return the solution to this LSE: :math:`x=A^{-1}b`. +from __future__ import annotations - Returns: - The solution to this LSE. +import cvxpy as cp +import numpy as np - Raises: - ValueError: if this LSE is parameterized with unassigned values. - """ - if self.A.ndim == 1: - # self.A is a vector of cp.Expression objects - mat_a = np.array([row.value for row in self.A]) - if any(row is None for row in mat_a): - raise ValueError( - "This LSE contains unassigned parameter values! Assign a value to them first " - "before trying to solve this LSE again." - ) - else: - mat_a = self.A - return cast(np.ndarray, np.linalg.solve(mat_a, self.b)) +from .costs import LSE -def setup_lse( +def setup_static_lse( trotter_steps: list[int] | cp.Parameter, *, order: int = 1, symmetric: bool = False, ) -> LSE: - r"""Return the linear system of equations for computing the static MPF coefficients. + r"""Return the linear system of equations for computing static MPF coefficients. This function constructs the following linear system of equations: @@ -87,8 +55,8 @@ def setup_lse( Here is an example: - >>> from qiskit_addon_mpf.static import setup_lse - >>> lse = setup_lse([1,2,3], order=2, symmetric=True) + >>> from qiskit_addon_mpf.static import setup_static_lse + >>> lse = setup_static_lse([1,2,3], order=2, symmetric=True) >>> print(lse.A) [[1. 1. 1. ] [1. 0.25 0.11111111] @@ -106,8 +74,19 @@ def setup_lse( symmetric: whether the individual product formulas making up the MPF are symmetric. For example, the Lie-Trotter formula is `not` symmetric, while Suzuki-Trotter `is`. + .. note:: + Making use of this value is equivalent to the static MPF coefficient description + provided by [1]. In contrast, [2] disregards the symmetry of the individual product + formulas, effectively always setting ``symmetric=False``. + Returns: - An :class:`.LSE`. + The :class:`.LSE` to find the static MPF coefficients as described above. + + References: + [1]: A. Carrera Vazquez et al., Quantum 7, 1067 (2023). + https://quantum-journal.org/papers/q-2023-07-25-1067/ + [2]: S. Zhuk et al., Phys. Rev. Research 6, 033309 (2024). + https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309 """ symmetric_factor = 2 if symmetric else 1 diff --git a/releasenotes/notes/dynamic-mpf-cf18e0b4d2bdeaff.yaml b/releasenotes/notes/dynamic-mpf-cf18e0b4d2bdeaff.yaml new file mode 100644 index 0000000..21a8429 --- /dev/null +++ b/releasenotes/notes/dynamic-mpf-cf18e0b4d2bdeaff.yaml @@ -0,0 +1,19 @@ +--- +features: + - | + Adds the ability to compute dynamic (i.e. time-dependent) MPF coefficients. + For more details, refer to :mod:`qiskit_addon_mpf.dynamic`. +upgrade: + - | + The code for the static MPF coefficients has been moved. + It functions the same as before, but you have to update your imports and + function names as summarized in the table below: + + =============================================================== =================================================================== + Old New + =============================================================== =================================================================== + ``from qiskit_addon_mpf.static import setup_lse`` ``from qiskit_addon_mpf.static import setup_static_lse`` + ``from qiskit_addon_mpf.static import LSE`` ``from qiskit_addon_mpf.costs import LSE`` + ``from qiskit_addon_mpf.static import setup_exact_model`` ``from qiskit_addon_mpf.costs import setup_exact_problem`` + ``from qiskit_addon_mpf.static import setup_approximate_model`` ``from qiskit_addon_mpf.costs import setup_sum_of_squares_problem`` + =============================================================== =================================================================== diff --git a/test/static/__init__.py b/test/backends/__init__.py similarity index 100% rename from test/static/__init__.py rename to test/backends/__init__.py diff --git a/qiskit_addon_mpf/static/__init__.py b/test/backends/quimb_circuit/__init__.py similarity index 60% rename from qiskit_addon_mpf/static/__init__.py rename to test/backends/quimb_circuit/__init__.py index adfa652..f1a7934 100644 --- a/qiskit_addon_mpf/static/__init__.py +++ b/test/backends/quimb_circuit/__init__.py @@ -9,17 +9,3 @@ # Any modifications or derivative works of this code must retain this # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. - -# Reminder: update the RST file in docs/apidocs when adding new interfaces. -"""Static MPFs.""" - -from .approximate import setup_approximate_model -from .exact import setup_exact_model -from .lse import LSE, setup_lse - -__all__ = [ - "LSE", - "setup_lse", - "setup_exact_model", - "setup_approximate_model", -] diff --git a/test/backends/quimb_circuit/test_e2e.py b/test/backends/quimb_circuit/test_e2e.py new file mode 100644 index 0000000..4313901 --- /dev/null +++ b/test/backends/quimb_circuit/test_e2e.py @@ -0,0 +1,86 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from functools import partial + +import numpy as np +import pytest +from qiskit.circuit import Parameter, QuantumCircuit +from qiskit.quantum_info import SparsePauliOp +from qiskit.synthesis import SuzukiTrotter +from qiskit_addon_mpf.backends import HAS_QUIMB +from qiskit_addon_mpf.costs import setup_frobenius_problem +from qiskit_addon_mpf.dynamic import setup_dynamic_lse +from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_circuit import CircuitEvolver, CircuitState + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestEndToEnd: + @pytest.mark.parametrize( + ["time", "expected_A", "expected_b", "expected_coeffs"], + [ + ( + 0.5, + [[1.0, 0.99961572], [0.99961572, 1.0]], + [0.99939447, 0.99804667], + [2.25365838, -1.25365812], + ), + ], + ) + def test_end_to_end(self, time, expected_A, expected_b, expected_coeffs): + np.random.seed(0) + + # constants + L = 6 + W = 0.5 + epsilon = 0.5 + + J = np.random.rand(L - 1) + W * np.ones(L - 1) + # ZZ couplings + Jz = 1.0 + # XX and YY couplings + Jxx = epsilon + hz = 0.000000001 * np.array([(-1) ** i for i in range(L)]) + + hamil = SparsePauliOp.from_sparse_list( + [("Z", [k], hz[k]) for k in range(L)] + + [("ZZ", [k, k + 1], J[k] * Jz) for k in range(L - 1)] + + [("YY", [k, k + 1], J[k] * Jxx) for k in range(L - 1)] + + [("XX", [k, k + 1], J[k] * Jxx) for k in range(L - 1)], + num_qubits=L, + ) + + dt = Parameter("dt") + suz_4 = generate_time_evolution_circuit(hamil, synthesis=SuzukiTrotter(order=4), time=dt) + suz_2 = generate_time_evolution_circuit(hamil, synthesis=SuzukiTrotter(order=2), time=dt) + + initial_state = QuantumCircuit(L) + for k in range(1, L, 2): + initial_state.x(k) + + model = setup_dynamic_lse( + [4, 3], + time, + CircuitState, + partial(CircuitEvolver, circuit=suz_4, dt=0.05), + partial(CircuitEvolver, circuit=suz_2), + initial_state, + ) + np.testing.assert_allclose(model.b, expected_b, rtol=1e-4) + np.testing.assert_allclose(model.A, expected_A, rtol=1e-4) + + prob, coeffs = setup_frobenius_problem(model) + prob.solve() + np.testing.assert_allclose(coeffs.value, expected_coeffs, rtol=1e-4) diff --git a/test/backends/quimb_circuit/test_state.py b/test/backends/quimb_circuit/test_state.py new file mode 100644 index 0000000..7c5a1f9 --- /dev/null +++ b/test/backends/quimb_circuit/test_state.py @@ -0,0 +1,36 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit_addon_mpf.backends import HAS_QUIMB + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_circuit import CircuitState + from quimb.tensor import MPS_neel_state + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestCircuitState: + def test_empty_state_handling(self): + """Test the handling of a non-evolved state.""" + state = CircuitState() + circ = QuantumCircuit(5) + with pytest.raises(RuntimeError): + state.overlap(circ) + + def test_unsupported_state(self): + """Test the handling of a non-supported state object.""" + state = CircuitState() + neel = MPS_neel_state(5) + with pytest.raises(TypeError): + state.overlap(neel) diff --git a/test/backends/quimb_layers/__init__.py b/test/backends/quimb_layers/__init__.py new file mode 100644 index 0000000..f1a7934 --- /dev/null +++ b/test/backends/quimb_layers/__init__.py @@ -0,0 +1,11 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/backends/quimb_layers/test_e2e.py b/test/backends/quimb_layers/test_e2e.py new file mode 100644 index 0000000..a1e438c --- /dev/null +++ b/test/backends/quimb_layers/test_e2e.py @@ -0,0 +1,183 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from functools import partial + +import numpy as np +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import XXPlusYYGate +from qiskit_addon_mpf.backends import HAS_QUIMB +from qiskit_addon_mpf.costs import setup_frobenius_problem +from qiskit_addon_mpf.dynamic import setup_dynamic_lse + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_layers import LayerModel, LayerwiseEvolver + from qiskit_addon_mpf.backends.quimb_tebd import MPOState, TEBDEvolver + from quimb.tensor import MPO_identity, MPS_neel_state, SpinHam1D + + +def gen_ext_field_layer(n, hz): + qc = QuantumCircuit(n) + for q in range(n): + qc.rz(-hz[q], q) + return qc + + +def trotter_step(qc, q, Jxx, Jz): + qc.rzz(Jz, q, q + 1) + qc.append(XXPlusYYGate(2.0 * Jxx), [q, q + 1]) + + +def gen_odd_coupling_layer(n, Jxx, Jz, J): + qc = QuantumCircuit(n) + for q in range(0, n - 1, 2): + trotter_step(qc, q, J[q] * Jxx, J[q] * Jz) + return qc + + +def gen_even_coupling_layer(n, Jxx, Jz, J): + qc = QuantumCircuit(n) + for q in range(1, n - 1, 2): + trotter_step(qc, q, J[q] * Jxx, J[q] * Jz) + return qc + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestEndToEnd: + @pytest.mark.parametrize( + ["time", "expected_A", "expected_b", "expected_coeffs"], + [ + ( + 0.5, + [[1.0, 0.99975555], [0.99975555, 1.0]], + [0.99958611, 0.99843783], + [2.84866034, -1.84866008], + ), + ( + 1.0, + [[1.0, 0.99184952], [0.99184952, 1.0]], + [0.98289207, 0.96107058], + [1.83866282, -0.83866282], + ), + ( + 1.5, + [[1.0, 0.95394736], [0.95394736, 1.0]], + [0.92251831, 0.76729061], + [2.18532883, -1.18532883], + ), + ], + ) + def test_end_to_end_custom_suzuki(self, time, expected_A, expected_b, expected_coeffs): + np.random.seed(0) + + # constants + L = 10 + W = 0.5 + epsilon = 0.5 + + J = np.random.rand(L - 1) + W * np.ones(L - 1) + # ZZ couplings + Jz = 1.0 + # XX and YY couplings + Jxx = epsilon + + # base coupling + # external field + hz = 0.000000001 * np.array([(-1) ** i for i in range(L)]) + + # Initialize the builder for a spin 1/2 chain + builder = SpinHam1D(S=1 / 2) + + # Add XX and YY couplings for neighboring sites + for i in range(L - 1): + builder[i, i + 1] += 2.0 * Jxx * J[i], "-", "+" + builder[i, i + 1] += 2.0 * Jxx * J[i], "+", "-" + + # Add ZZ couplings for neighboring sites + for i in range(L - 1): + builder[i, i + 1] += 4.0 * Jz * J[i], "Z", "Z" + + # Add the external Z-field (hz) to each site + for i in range(L): + builder[i] += -2.0 * hz[i], "Z" + + # Build the local Hamiltonian + exact_model = builder.build_local_ham(L) + + # NOTE: below we are building each layer at a time, but we could also have built a single + # Trotter circuit and sliced it using `qiskit_addon_utils.slicing`. + odd_coupling_layer = LayerModel.from_quantum_circuit( + gen_odd_coupling_layer(L, Jxx, Jz, J), + cyclic=False, + ) + even_coupling_layer = LayerModel.from_quantum_circuit( + gen_even_coupling_layer(L, Jxx, Jz, 2.0 * J), # factor 2 because its the central layer + cyclic=False, + ) + odd_onsite_layer = LayerModel.from_quantum_circuit( + gen_ext_field_layer(L, hz), + keep_only_odd=True, + scaling_factor=0.5, + cyclic=False, + ) + even_onsite_layer = LayerModel.from_quantum_circuit( + gen_ext_field_layer(L, hz), + keep_only_odd=False, + scaling_factor=0.5, + cyclic=False, + ) + # Our layers combine to form a second-order Suzuki-Trotter formula as follows: + layers = [ + odd_coupling_layer, + odd_onsite_layer, + even_onsite_layer, + even_coupling_layer, + even_onsite_layer, + odd_onsite_layer, + odd_coupling_layer, + ] + + split_opts = { + "max_bond": 10, + "cutoff": 1e-5, + "cutoff_mode": "rel", + "method": "svd", + "renorm": False, + } + + initial_state = MPS_neel_state(L) + + model = setup_dynamic_lse( + [4, 3], + time, + lambda: MPOState(MPO_identity(L)), + partial( + TEBDEvolver, + H=exact_model, + dt=0.05, + order=4, + split_opts=split_opts, + ), + partial( + LayerwiseEvolver, + layers=layers, + split_opts=split_opts, + ), + initial_state, + ) + np.testing.assert_allclose(model.b, expected_b, rtol=1e-4) + np.testing.assert_allclose(model.A, expected_A, rtol=1e-4) + + prob, coeffs = setup_frobenius_problem(model) + prob.solve() + np.testing.assert_allclose(coeffs.value, expected_coeffs, rtol=1e-4) diff --git a/test/backends/quimb_layers/test_model.py b/test/backends/quimb_layers/test_model.py new file mode 100644 index 0000000..563cb55 --- /dev/null +++ b/test/backends/quimb_layers/test_model.py @@ -0,0 +1,89 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + + +import numpy as np +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import XXPlusYYGate +from qiskit_addon_mpf.backends import HAS_QUIMB + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_layers import LayerModel + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestLayerModel: + def test_from_quantum_circuit(self): + L = 6 + + qc = QuantumCircuit(L) + for i in range(0, L - 1, 2): + qc.rzz(1.0, i, i + 1) + for i in range(L): + qc.rz(1.0, i) + for i in range(1, L - 1, 2): + qc.append(XXPlusYYGate(1.0), [i, i + 1]) + + qc = qc.repeat(2).decompose() + + model = LayerModel.from_quantum_circuit(qc) + expected_terms = { + (0, 1): np.array( + [ + [7.0, 0.0, 0.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, -3.0, 0.0], + [0.0, 0.0, 0.0, -5.0], + ] + ), + (2, 3): np.array( + [ + [5.0, 0.0, 0.0, 0.0], + [0.0, -1.0, 0.0, 0.0], + [0.0, 0.0, -1.0, 0.0], + [0.0, 0.0, 0.0, -3.0], + ] + ), + (1, 2): np.array( + [ + [4.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -4.0], + ] + ), + (4, 5): np.array( + [ + [7.0, 0.0, 0.0, 0.0], + [0.0, -3.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, -5.0], + ] + ), + (3, 4): np.array( + [ + [4.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -4.0], + ] + ), + } + for sites in model.terms: + np.testing.assert_allclose(model.terms[sites], expected_terms[sites]) + + def test_handling_unsupportedown_gate(self): + qc = QuantumCircuit(1) + qc.rx(1.0, 0) + with pytest.raises(NotImplementedError): + LayerModel.from_quantum_circuit(qc) diff --git a/test/backends/quimb_tebd/__init__.py b/test/backends/quimb_tebd/__init__.py new file mode 100644 index 0000000..f1a7934 --- /dev/null +++ b/test/backends/quimb_tebd/__init__.py @@ -0,0 +1,11 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/backends/quimb_tebd/test_e2e.py b/test/backends/quimb_tebd/test_e2e.py new file mode 100644 index 0000000..99a50d6 --- /dev/null +++ b/test/backends/quimb_tebd/test_e2e.py @@ -0,0 +1,124 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from functools import partial + +import numpy as np +import pytest +from qiskit_addon_mpf.backends import HAS_QUIMB +from qiskit_addon_mpf.costs import setup_frobenius_problem +from qiskit_addon_mpf.dynamic import setup_dynamic_lse + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_tebd import MPOState, TEBDEvolver + from quimb.tensor import MPO_identity, MPS_neel_state, SpinHam1D + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestEndToEnd: + @pytest.mark.parametrize( + ["time", "expected_A", "expected_b", "expected_coeffs"], + [ + ( + 0.5, + [[1.0, 0.9997562], [0.9997562, 1.0]], + [0.99952645, 0.99854528], + [2.51220065, -1.51220037], + ), + ( + 1.0, + [[1.0, 0.99189288], [0.99189288, 1.0]], + [0.98461028, 0.96466791], + [1.72992884, -0.72992884], + ), + ( + 1.5, + [[1.0, 0.95352741], [0.95352741, 1.0]], + [0.91618286, 0.79874277], + [1.76354156, -0.76354156], + ), + ], + ) + def test_end_to_end_builtin_suzuki(self, time, expected_A, expected_b, expected_coeffs): + np.random.seed(0) + + # constants + L = 10 + W = 0.5 + epsilon = 0.5 + + J = np.random.rand(L - 1) + W * np.ones(L - 1) + # ZZ couplings + Jz = 1.0 + # XX and YY couplings + Jxx = epsilon + + # base coupling + # external field + hz = 0.000000001 * np.array([(-1) ** i for i in range(L)]) + + # Initialize the builder for a spin 1/2 chain + builder = SpinHam1D(S=1 / 2) + + # Add XX and YY couplings for neighboring sites + for i in range(L - 1): + builder[i, i + 1] += 2.0 * Jxx * J[i], "-", "+" + builder[i, i + 1] += 2.0 * Jxx * J[i], "+", "-" + + # Add ZZ couplings for neighboring sites + for i in range(L - 1): + builder[i, i + 1] += 4.0 * Jz * J[i], "Z", "Z" + + # Add the external Z-field (hz) to each site + for i in range(L): + builder[i] += -2.0 * hz[i], "Z" + + # Build the local Hamiltonian + exact_model = builder.build_local_ham(L) + + split_opts = { + "max_bond": 10, + "cutoff": 1e-5, + "cutoff_mode": "rel", + "method": "svd", + "renorm": False, + } + + initial_state = MPS_neel_state(L) + + model = setup_dynamic_lse( + [4, 3], + time, + lambda: MPOState(MPO_identity(L)), + partial( + TEBDEvolver, + H=exact_model, + dt=0.05, + order=4, + split_opts=split_opts, + ), + partial( + TEBDEvolver, + H=exact_model, + order=2, + split_opts=split_opts, + ), + initial_state, + ) + np.testing.assert_allclose(model.b, expected_b, rtol=1e-3) + np.testing.assert_allclose(model.A, expected_A, rtol=1e-3) + + prob, coeffs = setup_frobenius_problem(model) + prob.solve() + # NOTE: this particular test converges to fairly different overlaps in the CI on MacOS only. + # Thus, the assertion threshold is so loose. + np.testing.assert_allclose(coeffs.value, expected_coeffs, rtol=0.1) diff --git a/test/backends/quimb_tebd/test_state.py b/test/backends/quimb_tebd/test_state.py new file mode 100644 index 0000000..c4edeb7 --- /dev/null +++ b/test/backends/quimb_tebd/test_state.py @@ -0,0 +1,29 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit_addon_mpf.backends import HAS_QUIMB + +if HAS_QUIMB: + from qiskit_addon_mpf.backends.quimb_tebd import MPOState + from quimb.tensor import MPO_identity + + +@pytest.mark.skipif(not HAS_QUIMB, reason="Quimb is required for these unittests") +class TestMPOState: + def test_unsupported_state(self): + """Test the handling of a non-supported state object.""" + state = MPOState(MPO_identity(5)) + circ = QuantumCircuit(5) + with pytest.raises(TypeError): + state.overlap(circ) diff --git a/test/backends/tenpy_layers/__init__.py b/test/backends/tenpy_layers/__init__.py new file mode 100644 index 0000000..f1a7934 --- /dev/null +++ b/test/backends/tenpy_layers/__init__.py @@ -0,0 +1,11 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/backends/tenpy_layers/test_e2e.py b/test/backends/tenpy_layers/test_e2e.py new file mode 100644 index 0000000..c24ada4 --- /dev/null +++ b/test/backends/tenpy_layers/test_e2e.py @@ -0,0 +1,181 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from functools import partial + +import numpy as np +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import XXPlusYYGate +from qiskit_addon_mpf.backends import HAS_TENPY +from qiskit_addon_mpf.costs import setup_frobenius_problem +from qiskit_addon_mpf.dynamic import setup_dynamic_lse + +if HAS_TENPY: + from qiskit_addon_mpf.backends.tenpy_layers import LayerModel, LayerwiseEvolver + from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state, TEBDEvolver + from tenpy.models import XXZChain2 + + +def gen_ext_field_layer(n, hz): + qc = QuantumCircuit(n) + for q in range(n): + qc.rz(-hz[q], q) + return qc + + +def trotter_step(qc, q, Jxx, Jz): + qc.rzz(Jz, q, q + 1) + qc.append(XXPlusYYGate(2.0 * Jxx), [q, q + 1]) + + +def gen_odd_coupling_layer(n, Jxx, Jz, J): + qc = QuantumCircuit(n) + for q in range(0, n - 1, 2): + trotter_step(qc, q, J[q] * Jxx, J[q] * Jz) + return qc + + +def gen_even_coupling_layer(n, Jxx, Jz, J): + qc = QuantumCircuit(n) + for q in range(1, n - 1, 2): + trotter_step(qc, q, J[q] * Jxx, J[q] * Jz) + return qc + + +@pytest.mark.skipif(not HAS_TENPY, reason="TeNPy is required for these unittests") +class TestEndToEnd: + @pytest.mark.parametrize( + ["time", "expected_A", "expected_b", "expected_coeffs"], + [ + ( + 0.5, + [[1.0, 0.9997562], [0.9997562, 1.0]], + [0.99944125, 0.99857914], + [2.2680558, -1.26805551], + ), + ( + 1.0, + [[1.0, 0.99189288], [0.99189288, 1.0]], + [0.98478594, 0.9676077], + [1.55870387, -0.55870387], + ), + ( + 1.5, + [[1.0, 0.95352741], [0.95352741, 1.0]], + [0.9032996, 0.71967399], + [2.47563369, -1.47563369], + ), + ], + ) + def test_end_to_end(self, time, expected_A, expected_b, expected_coeffs): + np.random.seed(0) + + # constants + L = 10 + W = 0.5 + epsilon = 0.5 + + J = np.random.rand(L - 1) + W * np.ones(L - 1) + # ZZ couplings + Jz = 1.0 + # XX and YY couplings + Jxx = epsilon + + # base coupling + # external field + hz = 0.000000001 * np.array([(-1) ** i for i in range(L)]) + + # This is the full model that we want to simulate. It is used for the "exact" time evolution + # (which is approximated via a fourth-order Suzuki-Trotter formula). + exact_model = XXZChain2( + { + "L": L, + "Jz": 4.0 * Jz * J, + "Jxx": 4.0 * Jxx * J, + "hz": 2.0 * hz, + "bc_MPS": "finite", + "sort_charge": False, + } + ) + + # NOTE: below we are building each layer at a time, but we could also have built a single + # Trotter circuit and sliced it using `qiskit_addon_utils.slicing`. + odd_coupling_layer = LayerModel.from_quantum_circuit( + gen_odd_coupling_layer(L, Jxx, Jz, J), + bc_MPS="finite", + ) + even_coupling_layer = LayerModel.from_quantum_circuit( + gen_even_coupling_layer(L, Jxx, Jz, 2.0 * J), # factor 2 because its the central layer + bc_MPS="finite", + ) + odd_onsite_layer = LayerModel.from_quantum_circuit( + gen_ext_field_layer(L, hz), + keep_only_odd=True, + scaling_factor=0.5, + bc_MPS="finite", + ) + even_onsite_layer = LayerModel.from_quantum_circuit( + gen_ext_field_layer(L, hz), + keep_only_odd=False, + scaling_factor=0.5, + bc_MPS="finite", + ) + # Our layers combine to form a second-order Suzuki-Trotter formula as follows: + layers = [ + odd_coupling_layer, + odd_onsite_layer, + even_onsite_layer, + even_coupling_layer, + even_onsite_layer, + odd_onsite_layer, + odd_coupling_layer, + ] + + options_common = { + "trunc_params": { + "chi_max": 10, + "svd_min": 1e-5, + "trunc_cut": None, + }, + "preserve_norm": False, + } + options_exact = options_common.copy() + options_exact["order"] = 4 + + options_approx = options_common.copy() + + initial_state = MPS_neel_state(exact_model.lat) + + model = setup_dynamic_lse( + [4, 3], + time, + partial(MPOState.initialize_from_lattice, exact_model.lat), + partial( + TEBDEvolver, + model=exact_model, + dt=0.05, + options=options_exact, + ), + partial( + LayerwiseEvolver, + layers=layers, + options=options_approx, + ), + initial_state, + ) + np.testing.assert_allclose(model.b, expected_b, rtol=1e-4) + np.testing.assert_allclose(model.A, expected_A, rtol=1e-4) + + prob, coeffs = setup_frobenius_problem(model) + prob.solve() + np.testing.assert_allclose(coeffs.value, expected_coeffs, rtol=1e-4) diff --git a/test/backends/tenpy_layers/test_model.py b/test/backends/tenpy_layers/test_model.py new file mode 100644 index 0000000..3290b53 --- /dev/null +++ b/test/backends/tenpy_layers/test_model.py @@ -0,0 +1,90 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + + +import numpy as np +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.library import XXPlusYYGate +from qiskit_addon_mpf.backends import HAS_TENPY + +if HAS_TENPY: + from qiskit_addon_mpf.backends.tenpy_layers import LayerModel + + +@pytest.mark.skipif(not HAS_TENPY, reason="Tenpy is required for these unittests") +class TestLayerModel: + def test_from_quantum_circuit(self): + L = 6 + + qc = QuantumCircuit(L) + for i in range(0, L - 1, 2): + qc.rzz(1.0, i, i + 1) + for i in range(L): + qc.rz(1.0, i) + for i in range(1, L - 1, 2): + qc.append(XXPlusYYGate(1.0), [i, i + 1]) + + qc = qc.repeat(2).decompose() + + model = LayerModel.from_quantum_circuit(qc) + expected_H_bonds = [ + np.array( + [ + [4.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [-2.0, 0.0, 0.0, -2.0], + ] + ), + np.array( + [ + [2.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -2.0], + ] + ), + np.array( + [ + [3.0, 0.0, 0.0, -1.0], + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [-1.0, 0.0, 0.0, -1.0], + ] + ), + np.array( + [ + [2.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -2.0], + ] + ), + np.array( + [ + [4.0, 0.0, 0.0, -2.0], + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, -2.0], + ] + ), + ] + assert model.H_bond[0] is None + for expected, actual in zip(expected_H_bonds, model.H_bond[1:]): + np.testing.assert_allclose(expected, actual.to_ndarray().reshape((4, 4))) + + def test_handling_unsupportedown_gate(self): + qc = QuantumCircuit(1) + qc.rx(1.0, 0) + with pytest.raises(NotImplementedError): + LayerModel.from_quantum_circuit(qc) diff --git a/test/backends/tenpy_tebd/__init__.py b/test/backends/tenpy_tebd/__init__.py new file mode 100644 index 0000000..f1a7934 --- /dev/null +++ b/test/backends/tenpy_tebd/__init__.py @@ -0,0 +1,11 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/backends/tenpy_tebd/test_e2e.py b/test/backends/tenpy_tebd/test_e2e.py new file mode 100644 index 0000000..dc5f961 --- /dev/null +++ b/test/backends/tenpy_tebd/test_e2e.py @@ -0,0 +1,120 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +from functools import partial + +import numpy as np +import pytest +from qiskit_addon_mpf.backends import HAS_TENPY +from qiskit_addon_mpf.costs import setup_frobenius_problem +from qiskit_addon_mpf.dynamic import setup_dynamic_lse + +if HAS_TENPY: + from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state, TEBDEvolver + from tenpy.models import XXZChain2 + + +@pytest.mark.skipif(not HAS_TENPY, reason="TeNPy is required for these unittests") +class TestEndToEnd: + @pytest.mark.parametrize( + ["time", "expected_A", "expected_b", "expected_coeffs"], + [ + ( + 0.5, + [[1.0, 0.9997562], [0.9997562, 1.0]], + [0.99944125, 0.99857914], + [2.26805572, -1.26805543], + ), + ( + 1.0, + [[1.0, 0.99189288], [0.99189288, 1.0]], + [0.98478594, 0.9676077], + [1.55870386, -0.55870386], + ), + ( + 1.5, + [[1.0, 0.95352741], [0.95352741, 1.0]], + [0.93918471, 0.71967399], + [2.8617227, -1.8617227], + ), + ], + ) + def test_end_to_end(self, time, expected_A, expected_b, expected_coeffs): + np.random.seed(0) + + # constants + L = 10 + W = 0.5 + epsilon = 0.5 + + J = np.random.rand(L - 1) + W * np.ones(L - 1) + # ZZ couplings + Jz = 1.0 + # XX and YY couplings + Jxx = epsilon + + # base coupling + # external field + hz = 0.000000001 * np.array([(-1) ** i for i in range(L)]) + + # This is the full model that we want to simulate. It is used for the "exact" time evolution + # (which is approximated via a fourth-order Suzuki-Trotter formula). + exact_model = XXZChain2( + { + "L": L, + "Jz": 4.0 * Jz * J, + "Jxx": 4.0 * Jxx * J, + "hz": 2.0 * hz, + "bc_MPS": "finite", + "sort_charge": False, + } + ) + + options_common = { + "trunc_params": { + "chi_max": 10, + "svd_min": 1e-5, + "trunc_cut": None, + }, + "preserve_norm": False, + } + options_exact = options_common.copy() + options_exact["order"] = 4 + + options_approx = options_common.copy() + options_approx["order"] = 2 + + initial_state = MPS_neel_state(exact_model.lat) + + model = setup_dynamic_lse( + [4, 3], + time, + partial(MPOState.initialize_from_lattice, exact_model.lat), + partial( + TEBDEvolver, + model=exact_model, + dt=0.05, + options=options_exact, + ), + partial( + TEBDEvolver, + model=exact_model, + options=options_approx, + ), + initial_state, + ) + np.testing.assert_allclose(model.b, expected_b, rtol=1e-4) + np.testing.assert_allclose(model.A, expected_A, rtol=1e-4) + + prob, coeffs = setup_frobenius_problem(model) + prob.solve() + np.testing.assert_allclose(coeffs.value, expected_coeffs, rtol=1e-4) diff --git a/test/backends/tenpy_tebd/test_state.py b/test/backends/tenpy_tebd/test_state.py new file mode 100644 index 0000000..b95d5d9 --- /dev/null +++ b/test/backends/tenpy_tebd/test_state.py @@ -0,0 +1,32 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +import pytest +from qiskit.circuit import QuantumCircuit +from qiskit_addon_mpf.backends import HAS_TENPY + +if HAS_TENPY: + from qiskit_addon_mpf.backends.tenpy_tebd import MPOState + from tenpy.models.lattice import Chain + from tenpy.networks.site import SpinHalfSite + + +@pytest.mark.skipif(not HAS_TENPY, reason="Tenpy is required for these unittests") +class TestMPOState: + def test_unsupported_state(self): + """Test the handling of a non-supported state object.""" + site = SpinHalfSite() + lattice = Chain(5, site) + state = MPOState.initialize_from_lattice(lattice) + circ = QuantumCircuit(5) + with pytest.raises(TypeError): + state.overlap(circ) diff --git a/test/costs/__init__.py b/test/costs/__init__.py new file mode 100644 index 0000000..f1a7934 --- /dev/null +++ b/test/costs/__init__.py @@ -0,0 +1,11 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. diff --git a/test/costs/test_exact.py b/test/costs/test_exact.py new file mode 100644 index 0000000..ba76d06 --- /dev/null +++ b/test/costs/test_exact.py @@ -0,0 +1,135 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the ``qiskit_addon_mpf.static.exact`` module.""" + +import cvxpy as cp +import numpy as np +import pytest +from qiskit_addon_mpf.costs.exact import setup_exact_problem +from qiskit_addon_mpf.static import setup_static_lse + + +class TestExactCoeffs: + """Tests for the ``qiskit_addon_mpf.static.exact`` module.""" + + def test_setup_exact_problem(self, subtests): + """Tests the :meth:`.setup_exact_problem` method.""" + trotter_steps = [1, 2, 4] + lse = setup_static_lse(trotter_steps, order=2, symmetric=True) + problem, coeffs = setup_exact_problem(lse) + + with subtests.test(msg="final cost"): + final_cost = problem.solve() + pytest.approx(final_cost, 1.888888888888) + + with subtests.test(msg="optimal coefficients"): + expected = np.array([0.02222222, -0.44444444, 1.42222222]) + np.testing.assert_allclose(coeffs.value, expected) + + @pytest.mark.parametrize( + ["trotter_steps", "expected_coeffs", "expected_cost"], + [ + # well-conditioned + ([1, 2], [-1.0, 2.0], 3.0), + ([1, 3], [-0.5, 1.5], 2.0), + ([2, 4], [-1.0, 2.0], 3.0), + ([2, 5], [-2 / 3, 5 / 3], 7 / 3), + ([1, 2, 6], [0.2, -1.0, 1.8], 3.0), + ([1, 2, 7], [1 / 6, -4 / 5, 49 / 30], 2.6), + # ill-conditioned + ([6, 7], [-6, 7], 13.0), + ([3, 4, 5, 6, 7], [27 / 8, -128 / 3, 625 / 4, -216, 2401 / 24], 518.3333332), + ( + [1, 2, 3, 4, 5, 6, 7], + [1 / 720, -8 / 15, 243 / 16, -1024 / 9, 15625 / 48, -3888 / 10, 117649 / 720], + 1007.22220288, + ), + ], + ) + def test_exact_order_1_references( + self, subtests, trotter_steps: list[int], expected_coeffs: list[float], expected_cost: float + ): + """This test ensures the correct behavior of the exact static MPF coefficient model. + + It does so, using the test-cases listed in Appendix E of [1]. + + [1]: A. Carrera Vazquez et al., Quantum 7, 1067 (2023). + https://quantum-journal.org/papers/q-2023-07-25-1067/ + """ + lse = setup_static_lse(trotter_steps, order=1) + problem, coeffs = setup_exact_problem(lse) + + with subtests.test(msg="final cost"): + final_cost = problem.solve() + pytest.approx(final_cost, expected_cost) + + with subtests.test(msg="optimal coefficients"): + np.testing.assert_allclose(coeffs.value, expected_coeffs) + + @pytest.mark.parametrize( + ["order", "trotter_steps", "expected_coeffs", "expected_cost"], + [ + (2, [1, 2], [-1.0 / 3.0, 4.0 / 3.0], 5.0 / 3.0), + (2, [1, 2, 4], [1.0 / 45.0, -4.0 / 9.0, 64.00 / 45.0], 85.0 / 45.0), + (2, [1, 2, 6], [1.0 / 105.0, -1.0 / 6.0, 81.0 / 70.0], 4.0 / 3.0), + (4, [1, 2], [-1.0 / 15.0, 16.0 / 15.0], 17.0 / 15.0), + (4, [1, 2, 3], [1 / 336.0, -32 / 105.0, 729 / 560.0], 1.60952381), + (4, [1, 2, 4], [1 / 945.0, -16 / 189.0, 1024 / 945.0], 1.16931217), + ], + ) + def test_exact_higher_order_references( + self, + subtests, + order: int, + trotter_steps: list[int], + expected_coeffs: list[float], + expected_cost: float, + ): + """This test ensures the correct behavior for higher order formulas. + + It does so, using some of the test-cases listed in Appendix A of [1]. + + [1]: G. H. Low et al, arXiv:1907.11679v2 (2019). + https://arxiv.org/abs/1907.11679v2 + """ + lse = setup_static_lse(trotter_steps, order=order, symmetric=True) + problem, coeffs = setup_exact_problem(lse) + + with subtests.test(msg="final cost"): + final_cost = problem.solve() + pytest.approx(final_cost, expected_cost) + + with subtests.test(msg="optimal coefficients"): + np.testing.assert_allclose(coeffs.value, expected_coeffs) + + def test_setup_exact_problem_params(self, subtests): + """Tests the :meth:`.setup_exact_problem` method with parameters.""" + ks = cp.Parameter(2) + lse = setup_static_lse(ks, order=1) + problem, coeffs = setup_exact_problem(lse) + + for trotter_steps, expected_coeffs, expected_cost in [ + ([1, 2], [-1.0, 2.0], 3.0), + ([1, 3], [-0.5, 1.5], 2.0), + ([2, 4], [-1.0, 2.0], 3.0), + ([2, 5], [-2 / 3, 5 / 3], 7 / 3), + ([6, 7], [-6, 7], 13.0), + ]: + ks.value = trotter_steps + + with subtests.test(msg="final cost"): + final_cost = problem.solve() + pytest.approx(final_cost, expected_cost) + + with subtests.test(msg="optimal coefficients"): + np.testing.assert_allclose(coeffs.value, expected_coeffs) diff --git a/test/costs/test_lse.py b/test/costs/test_lse.py new file mode 100644 index 0000000..2381029 --- /dev/null +++ b/test/costs/test_lse.py @@ -0,0 +1,38 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the ``qiskit_addon_mpf.costs.LSE`` class.""" + +import numpy as np +import pytest +from qiskit_addon_mpf.costs import LSE +from qiskit_addon_mpf.static import setup_static_lse + + +class TestLSE: + """Tests for the ``qiskit_addon_mpf.costs.LSE`` class.""" + + def test_lse_solve(self): + """Tests the :meth:`.LSE.solve` method.""" + trotter_steps = [1, 2, 4] + lse = setup_static_lse(trotter_steps, order=2, symmetric=True) + coeffs = lse.solve() + expected = np.array([0.022222222, -0.444444444, 1.422222222]) + np.testing.assert_allclose(coeffs, expected) + + def test_lse_solve_invalid(self): + """Tests the handling of an invalid model.""" + mat_a = np.asarray([[1.0, 0.9], [0.9, 1.0]]) + vec_b = np.asarray([0.9, 0.9]) + lse = LSE(mat_a, vec_b) + with pytest.raises(ValueError): + lse.solve() diff --git a/test/costs/test_sum_of_squares.py b/test/costs/test_sum_of_squares.py new file mode 100644 index 0000000..8a60c40 --- /dev/null +++ b/test/costs/test_sum_of_squares.py @@ -0,0 +1,75 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2024. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""Tests for the ``qiskit_addon_mpf.static.sum_of_squares`` module.""" + +from typing import Any, ClassVar + +import cvxpy as cp +import numpy as np +import pytest +from qiskit_addon_mpf.costs.sum_of_squares import setup_sum_of_squares_problem +from qiskit_addon_mpf.static import setup_static_lse + + +class TestSumOfSquaresCoeffs: + """Tests for the ``qiskit_addon_mpf.static.sum_of_squares`` module.""" + + CVXPY_SOLVER_SETTINGS: ClassVar[dict[str, Any]] = { + "warm_start": False, + "eps_abs": 1e-7, + "eps_rel": 1e-7, + } + + def test_setup_sum_of_squares_problem(self, subtests): + """Tests the :meth:`.setup_sum_of_squares_problem` method.""" + trotter_steps = [1, 2, 4] + lse = setup_static_lse(trotter_steps, order=2, symmetric=True) + problem, coeffs = setup_sum_of_squares_problem(lse) + + with subtests.test(msg="final cost"): + final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) + pytest.approx(final_cost, 0) + + with subtests.test(msg="optimal coefficients"): + expected = np.array([0.02222225, -0.44444444, 1.42222216]) + np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) + + def test_setup_sum_of_squares_problem_max_l1_norm(self, subtests): + """Tests the :meth:`.setup_sum_of_squares_problem` method with ``max_l1_norm``.""" + trotter_steps = [1, 2, 4] + lse = setup_static_lse(trotter_steps, order=2, symmetric=True) + problem, coeffs = setup_sum_of_squares_problem(lse, max_l1_norm=1.5) + + with subtests.test(msg="final cost"): + final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) + pytest.approx(final_cost, 0.00035765) + + with subtests.test(msg="optimal coefficients"): + expected = np.array([-0.001143293, -0.2488567, 1.25]) + np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) + + def test_setup_sum_of_squares_problem_params(self, subtests): + """Tests the :meth:`.setup_sum_of_squares_problem` method with parameters.""" + ks = cp.Parameter(3) + lse = setup_static_lse(ks, order=2, symmetric=True) + problem, coeffs = setup_sum_of_squares_problem(lse) + + ks.value = [1, 2, 4] + + with subtests.test(msg="final cost"): + final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) + pytest.approx(final_cost, 0) + + with subtests.test(msg="optimal coefficients"): + expected = np.array([0.02222225, -0.44444444, 1.42222216]) + np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) diff --git a/test/static/test_approximate.py b/test/static/test_approximate.py deleted file mode 100644 index 205fd5d..0000000 --- a/test/static/test_approximate.py +++ /dev/null @@ -1,79 +0,0 @@ -# This code is a Qiskit project. -# -# (C) Copyright IBM 2024. -# -# This code is licensed under the Apache License, Version 2.0. You may -# obtain a copy of this license in the LICENSE.txt file in the root directory -# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. -# -# Any modifications or derivative works of this code must retain this -# copyright notice, and modified files need to carry a notice indicating -# that they have been altered from the originals. - -"""Tests for the ``qiskit_addon_mpf.static.approximate`` module.""" - -import unittest -from typing import Any, ClassVar - -import cvxpy as cp -import numpy as np -from qiskit_addon_mpf.static.approximate import setup_approximate_model -from qiskit_addon_mpf.static.lse import setup_lse - - -class TestApproximateCoeffs(unittest.TestCase): - """Tests for the ``qiskit_addon_mpf.static.approximate`` module.""" - - CVXPY_SOLVER_SETTINGS: ClassVar[dict[str, Any]] = { - "warm_start": False, - "eps_abs": 1e-7, - "eps_rel": 1e-7, - } - - def test_setup_approximate_model(self): - """Tests the :meth:`.setup_approximate_model` method.""" - trotter_steps = [1, 2, 4] - lse = setup_lse(trotter_steps, order=2, symmetric=True) - problem, coeffs = setup_approximate_model(lse) - - with self.subTest("final cost"): - final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) - self.assertAlmostEqual(final_cost, 0) - - with self.subTest("optimal coefficients"): - expected = np.array([0.02222225, -0.44444444, 1.42222216]) - np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) - - def test_setup_approximate_model_max_l1_norm(self): - """Tests the :meth:`.setup_approximate_model` method with ``max_l1_norm``.""" - trotter_steps = [1, 2, 4] - lse = setup_lse(trotter_steps, order=2, symmetric=True) - problem, coeffs = setup_approximate_model(lse, max_l1_norm=1.5) - - with self.subTest("final cost"): - final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) - self.assertAlmostEqual(final_cost, 0.00035765) - - with self.subTest("optimal coefficients"): - expected = np.array([-0.001143293, -0.2488567, 1.25]) - np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) - - def test_setup_approximate_model_params(self): - """Tests the :meth:`.setup_approximate_model` method with parameters.""" - ks = cp.Parameter(3) - lse = setup_lse(ks, order=2, symmetric=True) - problem, coeffs = setup_approximate_model(lse) - - ks.value = [1, 2, 4] - - with self.subTest("final cost"): - final_cost = problem.solve(**self.CVXPY_SOLVER_SETTINGS) - self.assertAlmostEqual(final_cost, 0) - - with self.subTest("optimal coefficients"): - expected = np.array([0.02222225, -0.44444444, 1.42222216]) - np.testing.assert_allclose(coeffs.value, expected, rtol=1e-4) - - -if __name__ == "__main__": - unittest.main() diff --git a/test/static/test_exact.py b/test/static/test_exact.py deleted file mode 100644 index 0eeda7b..0000000 --- a/test/static/test_exact.py +++ /dev/null @@ -1,137 +0,0 @@ -# This code is a Qiskit project. -# -# (C) Copyright IBM 2024. -# -# This code is licensed under the Apache License, Version 2.0. You may -# obtain a copy of this license in the LICENSE.txt file in the root directory -# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. -# -# Any modifications or derivative works of this code must retain this -# copyright notice, and modified files need to carry a notice indicating -# that they have been altered from the originals. - -"""Tests for the ``qiskit_addon_mpf.static.exact`` module.""" - -import unittest - -import cvxpy as cp -import numpy as np -from ddt import data, ddt, unpack -from qiskit_addon_mpf.static.exact import setup_exact_model -from qiskit_addon_mpf.static.lse import setup_lse - - -@ddt -class TestExactCoeffs(unittest.TestCase): - """Tests for the ``qiskit_addon_mpf.static.exact`` module.""" - - def test_setup_exact_model(self): - """Tests the :meth:`.setup_exact_model` method.""" - trotter_steps = [1, 2, 4] - lse = setup_lse(trotter_steps, order=2, symmetric=True) - problem, coeffs = setup_exact_model(lse) - - with self.subTest("final cost"): - final_cost = problem.solve() - self.assertAlmostEqual(final_cost, 1.888888888888) - - with self.subTest("optimal coefficients"): - expected = np.array([0.02222222, -0.44444444, 1.42222222]) - np.testing.assert_allclose(coeffs.value, expected) - - @unpack - @data( - # well-conditioned - ([1, 2], [-1.0, 2.0], 3.0), - ([1, 3], [-0.5, 1.5], 2.0), - ([2, 4], [-1.0, 2.0], 3.0), - ([2, 5], [-2 / 3, 5 / 3], 7 / 3), - ([1, 2, 6], [0.2, -1.0, 1.8], 3.0), - ([1, 2, 7], [1 / 6, -4 / 5, 49 / 30], 2.6), - # ill-conditioned - ([6, 7], [-6, 7], 13.0), - ([3, 4, 5, 6, 7], [27 / 8, -128 / 3, 625 / 4, -216, 2401 / 24], 518.3333332), - ( - [1, 2, 3, 4, 5, 6, 7], - [1 / 720, -8 / 15, 243 / 16, -1024 / 9, 15625 / 48, -3888 / 10, 117649 / 720], - 1007.22220288, - ), - ) - def test_exact_order_1_references( - self, trotter_steps: list[int], expected_coeffs: list[float], expected_cost: float - ): - """This test ensures the correct behavior of the exact static MPF coefficient model. - - It does so, using the test-cases listed in Appendix E of [1]. - - [1]: A. Carrera Vazquez et al., Quantum 7, 1067 (2023). - https://quantum-journal.org/papers/q-2023-07-25-1067/ - """ - lse = setup_lse(trotter_steps, order=1) - problem, coeffs = setup_exact_model(lse) - - with self.subTest("final cost"): - final_cost = problem.solve() - self.assertAlmostEqual(final_cost, expected_cost) - - with self.subTest("optimal coefficients"): - np.testing.assert_allclose(coeffs.value, expected_coeffs) - - @unpack - @data( - (2, [1, 2], [-1.0 / 3.0, 4.0 / 3.0], 5.0 / 3.0), - (2, [1, 2, 4], [1.0 / 45.0, -4.0 / 9.0, 64.00 / 45.0], 85.0 / 45.0), - (2, [1, 2, 6], [1.0 / 105.0, -1.0 / 6.0, 81.0 / 70.0], 4.0 / 3.0), - (4, [1, 2], [-1.0 / 15.0, 16.0 / 15.0], 17.0 / 15.0), - (4, [1, 2, 3], [1 / 336.0, -32 / 105.0, 729 / 560.0], 1.60952381), - (4, [1, 2, 4], [1 / 945.0, -16 / 189.0, 1024 / 945.0], 1.16931217), - ) - def test_exact_higher_order_references( - self, - order: int, - trotter_steps: list[int], - expected_coeffs: list[float], - expected_cost: float, - ): - """This test ensures the correct behavior for higher order formulas. - - It does so, using some of the test-cases listed in Appendix A of [1]. - - [1]: G. H. Low et al, arXiv:1907.11679v2 (2019). - https://arxiv.org/abs/1907.11679v2 - """ - lse = setup_lse(trotter_steps, order=order, symmetric=True) - problem, coeffs = setup_exact_model(lse) - - with self.subTest("final cost"): - final_cost = problem.solve() - self.assertAlmostEqual(final_cost, expected_cost) - - with self.subTest("optimal coefficients"): - np.testing.assert_allclose(coeffs.value, expected_coeffs) - - def test_setup_exact_model_params(self): - """Tests the :meth:`.setup_exact_model` method with parameters.""" - ks = cp.Parameter(2) - lse = setup_lse(ks, order=1) - problem, coeffs = setup_exact_model(lse) - - for trotter_steps, expected_coeffs, expected_cost in [ - ([1, 2], [-1.0, 2.0], 3.0), - ([1, 3], [-0.5, 1.5], 2.0), - ([2, 4], [-1.0, 2.0], 3.0), - ([2, 5], [-2 / 3, 5 / 3], 7 / 3), - ([6, 7], [-6, 7], 13.0), - ]: - ks.value = trotter_steps - - with self.subTest("final cost"): - final_cost = problem.solve() - self.assertAlmostEqual(final_cost, expected_cost) - - with self.subTest("optimal coefficients"): - np.testing.assert_allclose(coeffs.value, expected_coeffs) - - -if __name__ == "__main__": - unittest.main() diff --git a/test/static/test_lse.py b/test/test_static.py similarity index 73% rename from test/static/test_lse.py rename to test/test_static.py index bc88abc..65bd536 100644 --- a/test/static/test_lse.py +++ b/test/test_static.py @@ -12,44 +12,43 @@ """Tests for the ``qiskit_addon_mpf.static.lse`` module.""" -import unittest - import cvxpy as cp import numpy as np -from qiskit_addon_mpf.static.lse import setup_lse +import pytest +from qiskit_addon_mpf.static import setup_static_lse -class TestLSE(unittest.TestCase): +class TestLSE: """Tests for the ``qiskit_addon_mpf.static.lse`` module.""" - def test_setup_lse(self): - """Tests the :meth:` setup_lse` method.""" + def test_setup_static_lse(self, subtests): + """Tests the :meth:` setup_static_lse` method.""" trotter_steps = [1, 2, 4] expected_b = np.array([1, 0, 0]) - with self.subTest("order=1"): - lse = setup_lse(trotter_steps, order=1) + with subtests.test(msg="order=1"): + lse = setup_static_lse(trotter_steps, order=1) expected_A = np.array([[1.0, 1.0, 1.0], [1.0, 0.5, 0.25], [1.0, 0.25, 0.0625]]) np.testing.assert_allclose(lse.A, expected_A) np.testing.assert_allclose(lse.b, expected_b) expected_c = np.array([0.33333333, -2.0, 2.66666667]) np.testing.assert_allclose(lse.solve(), expected_c) - with self.subTest("order=2"): - lse = setup_lse(trotter_steps, order=2, symmetric=True) + with subtests.test(msg="order=2"): + lse = setup_static_lse(trotter_steps, order=2, symmetric=True) expected_A = np.array([[1.0, 1.0, 1.0], [1.0, 0.25, 0.0625], [1.0, 0.0625, 0.00390625]]) np.testing.assert_allclose(lse.A, expected_A) np.testing.assert_allclose(lse.b, expected_b) expected_c = np.array([0.022222222, -0.444444444, 1.422222222]) np.testing.assert_allclose(lse.solve(), expected_c) - def test_lse_with_params(self): - """Tests the :meth:` setup_lse` method with parameters.""" + def test_lse_with_params(self, subtests): + """Tests the :meth:` setup_static_lse` method with parameters.""" trotter_steps = cp.Parameter(3) - lse = setup_lse(trotter_steps) + lse = setup_static_lse(trotter_steps) - with self.subTest("assert ValueError"), self.assertRaises(ValueError): + with subtests.test(msg="assert ValueError"), pytest.raises(ValueError): lse.solve() trotter_steps.value = [1, 2, 4] @@ -60,7 +59,3 @@ def test_lse_with_params(self): expected_c = np.array([0.33333333, -2.0, 2.66666667]) np.testing.assert_allclose(lse.solve(), expected_c) - - -if __name__ == "__main__": - unittest.main()