diff --git a/.github/workflows/v2-build-and-deploy.yml b/.github/workflows/v2-build-and-deploy.yml new file mode 100644 index 0000000000..cbfbe1763d --- /dev/null +++ b/.github/workflows/v2-build-and-deploy.yml @@ -0,0 +1,101 @@ +name: V2 Build and Deploy Demos + +permissions: + contents: read + pull-requests: write + actions: read + +on: + workflow_dispatch: + inputs: + environment: + description: SWC environment to deploy to + options: + - swc-staging + - swc-prod + - swc-dev + required: true + type: choice + target: + default: stable + description: PennyLane version to build the demos. Either 'latest' or the most recent 'stable' release. + options: + - latest + - stable + required: true + type: choice + demos: + description: Demos to build and deploy, space-separated list of slugs (e.g. "demo1 demo2 demo3"), or leave empty for all demos. + required: false + type: string + as-previews: + default: false + description: | + Whether to deploy the demos as previews. + + **Please note** that demos built with the latest version cannot be published to swc-staging or swc-prod. + They can only be deployed as previews. + required: false + type: boolean + +jobs: + validate-and-parse-inputs: + runs-on: ubuntu-latest + outputs: + branch: ${{ steps.set-branch.outputs.branch }} + steps: + - name: Set branch + id: set-branch + run: | + if [[ "${{ github.event.inputs.target }}" == "stable" ]]; then + echo "branch=master" >> $GITHUB_OUTPUT + elif [[ "${{ github.event.inputs.target }}" == "latest" ]]; then + echo "branch=dev" >> $GITHUB_OUTPUT + else + echo "branch=" >> $GITHUB_OUTPUT + fi + + - name: Validate preview input + id: validate-preview + run: | + if [[ + ("${{ github.event.inputs.environment }}" == "swc-staging" || + "${{ github.event.inputs.environment }}" == "swc-prod") && + "${{ github.event.inputs.target }}" == "latest" && + "${{ github.event.inputs.as-previews }}" == "false" + ]]; then + echo "=========================" + echo "đźš« Invalid input detected:" + echo "Demos built with the latest version cannot be published to 'swc-staging' or 'swc-prod'." + echo "They can only be deployed as previews." + echo "Please set the 'as-previews' input to 'true' in your workflow configuration." + echo "=========================" + exit 1 + fi + + build: + needs: validate-and-parse-inputs + if: > + (needs.validate-and-parse-inputs.outputs.branch == 'master') || + (needs.validate-and-parse-inputs.outputs.branch == 'dev') + uses: ./.github/workflows/v2-build-demos.yml + with: + ref: ${{ needs.validate-and-parse-inputs.outputs.branch }} + demo-names: ${{ github.event.inputs.demos }} + dev: ${{ github.event.inputs.target == 'latest' }} + save-artifact: true + artifact-name: build-and-deploy-${{ github.event.inputs.target }} + keep-going: false + quiet: false + batch_size: 10 + + deploy: + uses: ./.github/workflows/v2-deploy-demos.yml + needs: + - validate-and-parse-inputs + - build + secrets: inherit + with: + environment: ${{ github.event.inputs.environment }} + artifact-name: build-and-deploy-${{ github.event.inputs.target }} + preview: ${{ github.event.inputs.as-previews }} \ No newline at end of file diff --git a/.github/workflows/v2-build-demos.yml b/.github/workflows/v2-build-demos.yml index 71d857dad3..2f60cc73a8 100644 --- a/.github/workflows/v2-build-demos.yml +++ b/.github/workflows/v2-build-demos.yml @@ -8,13 +8,20 @@ on: type: string demo-names: description: | - Instead of building all demos, only build the demos specified in this list. + Only build the demos specified in a space-separated list. + e.g. demo1 demo2 demo3 + + If not specified, all demos will be built. required: false type: string default: '' + dev: + description: Use development dependencies + required: false + type: boolean + default: false save-artifact: - description: | - Whether to save the built demos as an artifact. + description: Save the built demos as an artifact. required: false type: boolean default: true @@ -24,22 +31,19 @@ on: type: string default: '' keep-going: - description: Whether to keep going if a demo fails to build + description: Keep going if a demo fails to build required: false type: boolean default: false quiet: - description: | - Whether to suppress output from the build process + description: Suppress output from the build process required: false type: boolean default: false batch_size: - description: | - Number of demos to build per job + description: Number of demos to build per job type: number default: 10 - outputs: artifact-name: description: "Name of the artifact containing the built demos" @@ -57,25 +61,28 @@ on: required: false type: string default: '' + dev: + description: Use development dependencies + required: false + type: boolean + default: false save-artifact: - description: | - Whether to save the built demos as an artifact. + description: Save the built demos as an artifact. required: false type: boolean default: true artifact-name: - description: Name of the artifact containing the built demos (defaults to demo-build-{ref}) - required: true + description: Name of the build artifact (defaults to demo-build-{ref}) + required: false type: string default: '' keep-going: - description: Whether to keep going if a demo fails to build + description: Keep going if a demo fails to build required: false type: boolean default: false quiet: - description: | - Whether to suppress output from the build process + description: Suppress output from the build process required: false type: boolean default: false @@ -165,7 +172,8 @@ jobs: --execute \ ${{ inputs.keep-going && '--keep-going' || '--no-keep-going' }} \ ${{ inputs.quiet && '--quiet' || '--no-quiet' }} \ - ${{ steps.demo-name-list.outputs.demo-names }} + ${{ steps.demo-name-list.outputs.demo-names }} \ + ${{ inputs.dev && '--dev' || '--no-dev' }} - name: Upload artifacts id: upload-artifacts @@ -207,5 +215,3 @@ jobs: uses: geekyeggo/delete-artifact@v5 with: name: ${{ needs.generate-build-variables.outputs.artifact-name }}-c* - - diff --git a/.github/workflows/v2-build-pr.yml b/.github/workflows/v2-build-pr.yml new file mode 100644 index 0000000000..6de172691a --- /dev/null +++ b/.github/workflows/v2-build-pr.yml @@ -0,0 +1,82 @@ +name: V2 Build PR + +on: + pull_request: + # remove v2 after testing + branches: [master, dev, v2] + +permissions: + contents: read + +concurrency: + group: build-v2-demos-${{ github.event.pull_request.head.sha }} + cancel-in-progress: true + +jobs: + # Step 1: Identify changed demos + identify-changed-demos: + runs-on: ubuntu-latest + outputs: + updated: ${{ steps.get-changed-demos.outputs.updated }} + deleted: ${{ steps.get-changed-demos.outputs.deleted }} + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Get changed demos + id: get-changed-demos + uses: ./.github/actions/get-changed-demos + + - name: Output changed demos + run: | + echo "Updated Demos: ${{ steps.get-changed-demos.outputs.updated }}" + echo "Deleted Demos: ${{ steps.get-changed-demos.outputs.deleted }}" + + - name: Exit if no changes + if: steps.get-changed-demos.outputs.updated == '' + run: | + echo "No changes found in demos. Exiting workflow." + + # Step 2: Build demos + build: + if: needs.identify-changed-demos.outputs.updated != '' + uses: ./.github/workflows/v2-build-demos.yml + needs: + - identify-changed-demos + with: + ref: ${{ github.event.pull_request.head.sha }} + demo-names: ${{ needs.identify-changed-demos.outputs.updated }} + dev: ${{ github.event.pull_request.base.ref == 'dev' }} + save-artifact: true + artifact-name: demo-build-${{ github.event.pull_request.head.sha }} + keep-going: false + quiet: false + batch_size: 10 + + # Step 3: Save build context + save-build-context: + runs-on: ubuntu-latest + needs: + - build + - identify-changed-demos + steps: + - name: Save Pull Request Event Context + run: | + mkdir -p /tmp/pr + cat >/tmp/pr/pr_info.json <> $GITHUB_OUTPUT + fi + + - name: Unpack Build Information + if: steps.build_context.outputs.result != '' + run: unzip ${{ steps.build_context.outputs.result }} + + - name: Read Build Information + id: read_build_info + if: steps.build_context.outputs.result != '' + uses: actions/github-script@v7 + with: + script: | + const fs = require('fs'); + const buildData = fs.readFileSync('pr_info.json', 'utf8'); + return JSON.parse(buildData); + + - name: Parse Pull Request Event Information + id: pr_info + if: github.event.workflow_run.event == 'pull_request' && steps.build_context.outputs.result != '' + run: | + echo '${{ steps.read_build_info.outputs.result }}' | jq -r '.id' > pr_id.txt + echo '${{ steps.read_build_info.outputs.result }}' | jq -r '.ref' > pr_ref.txt + echo '${{ steps.read_build_info.outputs.result }}' | jq -r '.ref_name' > pr_ref_name.txt + echo '${{ steps.read_build_info.outputs.result }}' | jq -c '.updated_demos' > updated_demos.json + echo '${{ steps.read_build_info.outputs.result }}' | jq -c '.deleted_demos' > deleted_demos.json + + echo "pr_id=$(cat pr_id.txt)" >> $GITHUB_OUTPUT + echo "pr_ref=$(cat pr_ref.txt)" >> $GITHUB_OUTPUT + echo "pr_ref_name=$(cat pr_ref_name.txt)" >> $GITHUB_OUTPUT + echo "updated_demos=$(cat updated_demos.json)" >> $GITHUB_OUTPUT + echo "deleted_demos=$(cat deleted_demos.json)" >> $GITHUB_OUTPUT + + - name: Set job outputs + if: github.event.workflow_run.event == 'pull_request' && steps.build_context.outputs.result != '' + id: set_job_outputs + run: | + echo "pr_id=${{ steps.pr_info.outputs.pr_id }}" >> $GITHUB_OUTPUT + echo "pr_ref=${{ steps.pr_info.outputs.pr_ref }}" >> $GITHUB_OUTPUT + echo "pr_ref_name=${{ steps.pr_info.outputs.pr_ref_name }}" >> $GITHUB_OUTPUT + echo "updated_demos=${{ steps.pr_info.outputs.updated_demos }}" >> $GITHUB_OUTPUT + echo "deleted_demos=${{ steps.pr_info.outputs.deleted_demos }}" >> $GITHUB_OUTPUT + outputs: + pr_id: ${{ steps.set_job_outputs.outputs.pr_id }} + pr_ref: ${{ steps.set_job_outputs.outputs.pr_ref }} + pr_ref_name: ${{ steps.set_job_outputs.outputs.pr_ref_name }} + updated_demos: ${{ steps.set_job_outputs.outputs.updated_demos }} + deleted_demos: ${{ steps.set_job_outputs.outputs.deleted_demos }} + + # Step 2: Deploy the demos to SWC + deploy-preview-demos: + if: github.event.workflow_run.event == 'pull_request' && needs.prepare-build-context.result == 'success' + uses: ./.github/workflows/v2-deploy-demos.yml + needs: prepare-build-context + with: + # TODO: Update SWC environment to "swc-prod" + environment: 'swc-staging' + artifact-name: demo-build-${{ needs.prepare-build-context.outputs.pr_ref }} + preview: true + secrets: inherit + + + ## DISABLED UNTIL RELEASED ## + # Step 3: Create a comment on the PR with the demo links + # generate-comment: + # if: github.event.workflow_run.event == 'pull_request' && needs.prepare-build-context.outputs.pr_id != '' + # runs-on: ubuntu-latest + # needs: prepare-build-context + # steps: + # - name: Create markdown comment from demo names + # id: generate-markdown + # run: | + # demos="${{ needs.prepare-build-context.outputs.updated_demos }}" + + # comment="### Preview(s) are ready! :tada:\n" + # comment+="
\n" + # comment+="Toggle to view preview links\n" + # comment+="\n" + # # TODO: Switch to prod once testing is complete + # for demo in $demos; do + # comment+="- [$demo](https://staging.pennylane.ai/qml/demos/$demo?preview=true)\n" + # done + + # comment+="\n" + # comment+="
" + + # echo "markdown=$comment" >> $GITHUB_OUTPUT + + # - name: Comment on PR + # id: comment-on-pr + # uses: XanaduAI/cloud-actions/create-and-update-pull-request-comment@main + # with: + # github_token: ${{ secrets.github_token }} + # pull_request_number: ${{ needs.prepare-build-context.outputs.pr_id }} + # comment_body: ${{ steps.generate-markdown.outputs.markdown }} diff --git a/README.md b/README.md index 04ad0121c5..c48912f7d0 100644 --- a/README.md +++ b/README.md @@ -33,10 +33,18 @@ and can be downloaded as Jupyter notebooks and Python scripts. You can contribute by submitting a demo via a pull request implementing a recent quantum computing paper/result. + +### Installing the QML tool. + +The `qml` command line tool can be installed by running `pip install .` in the repo root. + ### Adding demos +- Run `qml new` to create a new demo interactively. + - Demos are written in the form of an executable Python script. - - Packages listed in `pyproject.toml` will be available for import during execution. + - Packages listed in `dependencies/requirements-core.in` will be available to all demos by default. + Extra dependencies can be named using a `requirements.in` file in the demo directory. See section below on `Dependency Management` for more details. - Matplotlib plots will be automatically rendered and displayed on the QML website. @@ -47,7 +55,7 @@ quantum computing paper/result. [these steps](https://github.com/PennyLaneAI/qml/tree/master/notebook_converter). - All demos should have a file name beginning with `tutorial_`. - The python files are saved in the `demonstrations` directory. + The python files are saved in the `demonstrations_v2/{demo_name}` directory. - The new demos will avoid using `autograd` or `TensorFlow`, `Jax` and `torch` are recommended instead. Also, if possible, the use of `lightning.qubit` is recommended. - [Restructured Text](http://www.sphinx-doc.org/en/master/usage/restructuredtext/basics.html) @@ -81,10 +89,7 @@ quantum computing paper/result. * To maintain quality and performance, every image should be twice (2X) its visible dimension size on the web page, and at a minimum of `150 ppi/dpi` (preferably `300 ppi/dpi`).
-- Lastly, your demo will need an accompanying _metadata_ file. This file should be named - the same as your python file, but with the `.py` extension replaced with - `.metadata.json`. Check out the `demonstrations_metadata.md` file in this repo for - details on how to format that file and what to include. +- Lastly, your demo will need an accompanying _metadata_ file in the demo directory. - Include the author's information in the `.metadata.json` file. This can be either the author's PennyLane profile `username` or their `name`. If you provide the PennyLane profile username, the author details will be sourced directly from that profile, and the demo will then appear as a contribution on the author's profile. @@ -161,138 +166,23 @@ there are a couple of guidelines to keep in mind. - All submissions must pass code review before being merged into the repository. ## Dependency Management -Due to the large scope of requirements in this repository, the traditional `requirements.txt` file is being phased out -and `pyproject.toml` is being introduced instead, the goal being easier management in regard to adding/updating packages. -To install all the dependencies locally, [poetry](https://python-poetry.org/) needs to be installed. Please follow the -[official installation documentation](https://python-poetry.org/docs/#installation). +Demo dependencies are automatically installed by the `qml` tool during demo +execution. See [dependencies/README.md](./dependencies/README.md) for details +on dependency specifications. ### Installing dependencies -Once poetry has been installed, the dependencies can be installed as follows: -```bash -make environment -``` -Note: This makefile target calls `poetry install` under the hood, you can pass any poetry arguments to this by passing -the `POETRYOPTS` variable. -```bash -make environment POETRYOPTS='--sync --dry-run --verbose' -``` - -The `master` branch of QML uses the latest stable release of PennyLane, whereas the `dev` branch uses the most -up-to-date version from the GitHub repository. If your demo relies on that, install the `dev` dependencies instead -by upgrading all PennyLane and its various plugins to the latest commit from GitHub. -```bash -# Run this instead of running the command above -make environment UPGRADE_PL=true -``` - -#### Installing only the dependencies to build the website without executing demos -It is possible to build the website without executing any of the demo code using `make html-norun` (More details below). - -To install only the base dependencies without the executable dependencies, use: -```bash -make environment BASE_ONLY=true -``` -(This is the equivalent to the previous method of `pip install -r requirements_norun.txt`). - -### Adding / Editing dependencies - -All dependencies need to be added to the pyproject.toml. It is recommended that unless necessary, -all dependencies be pinned to as tight of a version as possible. - -Add the new dependency in the `[tool.poetry.group.executable-dependencies.dependencies]` section of the toml file. - -Once pyproject.toml files have been updated, the poetry.lock file needs to be refreshed: -```bash -poetry lock --no-update -``` -This command will ensure that there are no dependency conflicts with any other package, and everything works. - -The `--no-update` ensures existing package versions are not bumped as part of the locking process. - -If the dependency change is required in prod, open the PR against `master`, or if it's only required in dev, then open -the PR against the `dev` branch, which will be synced to master on the next release of PennyLane. - -#### Adding / Editing PennyLane (or plugin) versions -This process is slightly different from other packages. It is due to the fact that the `master` builds use the stable -releases of PennyLane as stated in the pyproject.toml file. However, for dev builds, we use the latest commit from -GitHub. - -##### Adding a new PennyLane package (plugin) -- Add the package to `pyproject.toml` file with the other pennylane packages and pin it to the latest stable release. -- Add the GitHub installation link to the Makefile, so it is upgraded for dev builds with the other PennyLane packages. - - This should be under the format `$$PYTHON_VENV_PATH/bin/python -m pip install --upgrade git+https://github.com/PennyLaneAI/.git#egg=;\` -- Refresh the poetry lock file by running `poetry lock` +Dependencies are automatically installed when executing demos. A `requirements.txt` file +will be created in the `demo` directory after the build. ## Building -To build the website locally, simply run `make html`. The rendered HTML files -will now be available in `_build/html`. Open `_build/html/index.html` to browse -the built site locally. - -Note that the above command may take some time, as all demos -will be executed and built! Once built, only _modified_ demos will -be re-executed/re-built. - -Alternatively, you may run `make html-norun` to build the website _without_ executing -demos, or build only a single demo using the following command: - -```console -sphinx-build -D sphinx_gallery_conf.filename_pattern=tutorial_QGAN\.py -b html . _build -``` - -where `tutorial_QGAN` should be replaced with the name of the demo to build. - -## Building and running locally on Mac (M1) - -To install dependencies on an M1 Mac and build the QML website, the following instructions may be useful. - -- If python3 is not currently installed, we recommend you install via [Homebrew](https://github.com/conda-forge/miniforge): - - ```bash - brew install python - ``` - -- Follow the steps from `Dependency Management` to setup poetry. - -- Install the base packages by running - - ```bash - make environment BASE_ONLY=true - ``` - - Alternatively, you can do this in a new virtual environment using - - ```bash - python -m venv [venv_name] - cd [venv_name] && source bin/activate - make environment BASE_ONLY=true - ``` - -Once this is complete, you should be able to build the website using `make html-norun`. If this succeeds, the `build` folder should be populated with files. Open `index.html` in your browser to view the built site. - -If you are running into the error message - -``` -command not found: sphinx-build -``` - -you may need to make the following change: - -- In the `Makefile` change `SPHINXBUILD = sphinx-build` to `SPHINXBUILD = python3 -m sphinx.cmd.build`. - -If you are running into the error message - -``` -ModuleNotFoundError: No module named 'the-module-name' -``` +To build and execute demos locally, use `qml build --execute {demo_name} ...`. -you may need to install the module manually: +To build demos using dev dependencies, use `qml build --execute --dev {demo_name}`. -``` -pip3 install the-module-name -``` +Run `qml build --help` for more details. ## Support diff --git a/constraints.txt b/constraints.txt index ad733817c3..d968aba3fe 100644 --- a/constraints.txt +++ b/constraints.txt @@ -1,8 +1,8 @@ -pennylane==0.40.0 -pennylane-cirq==0.40.0 -pennylane-qiskit==0.40.0 -pennylane-qulacs==0.40.0 -pennylane-catalyst==0.10.0 +pennylane==0.41.0 +pennylane-cirq==0.41.0 +pennylane-qiskit==0.41.0 +pennylane-qulacs==0.41.0 +pennylane-catalyst==0.11.0 pennylane-qrack==0.11.1 pyqrack==1.32.12 numpy~=1.24 diff --git a/demonstrations_v2/.gitignore b/demonstrations_v2/.gitignore new file mode 100644 index 0000000000..4414fc1e28 --- /dev/null +++ b/demonstrations_v2/.gitignore @@ -0,0 +1 @@ +requirements.txt diff --git a/demonstrations_v2/quantum_neural_net/demo.py b/demonstrations_v2/quantum_neural_net/demo.py index 020890e5a8..2d8a999233 100644 --- a/demonstrations_v2/quantum_neural_net/demo.py +++ b/demonstrations_v2/quantum_neural_net/demo.py @@ -17,7 +17,7 @@ pytorch_noise PyTorch and noisy devices tutorial_noisy_circuit_optimization Optimizing noisy circuits with Cirq -*Author: Maria Schuld — Posted: 11 October 2019. Last updated: 25 January 2021.* +*Author: Maria Schuld — Posted: 11 October 2019. Last updated: 20 March 2025.* .. warning:: This demo is only compatible with PennyLane version ``0.29`` or below. @@ -84,7 +84,7 @@ def quantum_neural_net(var, x): for v in var: layer(v) - return qml.expval(qml.QuadX(0)) + return qml.expval(qml.X(0)) ############################################################################## @@ -741,4 +741,4 @@ def cost(var, features, labels): ############################################################################## # About the author # ---------------- -# \ No newline at end of file +# diff --git a/demonstrations_v2/quantum_neural_net/metadata.json b/demonstrations_v2/quantum_neural_net/metadata.json index 15ea1533f9..7c1c3b998c 100644 --- a/demonstrations_v2/quantum_neural_net/metadata.json +++ b/demonstrations_v2/quantum_neural_net/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2019-10-11T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2025-03-20T00:00:00+00:00", "categories": [ "Quantum Machine Learning" ], diff --git a/demonstrations_v2/tutorial_annni/demo.py b/demonstrations_v2/tutorial_annni/demo.py new file mode 100644 index 0000000000..95aca197ba --- /dev/null +++ b/demonstrations_v2/tutorial_annni/demo.py @@ -0,0 +1,663 @@ +r"""Supervised and Unsupervised Quantum Machine Learning Models for the Phase Detection of the ANNNI Spin Model +========================================================== + +Quantum Machine Learning (QML) models provide a natural framework for analyzing quantum many-body systems due to the direct mapping between spins and qubits. + +A key property of these systems is their phase diagram, which characterizes different physical phases based on parameters such as magnetization or the strength of interactions. +These diagrams help identify transitions where the system undergoes sudden changes in behaviour, which is essential for understanding the fundamental properties of quantum materials and optimizing quantum technologies. + +Following the same approach as in the paper on `Quantum phase detection generalization from marginal quantum neural network models `_, +we explore the phase diagram of the Axial Next-Nearest-Neighbor Ising (ANNNI) model using both supervised and unsupervised learning methods: + +* **Supervised model**: The Quantum Convolutional Neural Network (QCNN) is trained on a small subset of analytically known phase points, demonstrating its ability to generalize beyond the training data. + +* **Unsupervised model**: Quantum Anomaly Detection (QAD) requires minimal prior knowledge and efficiently identifies potential phase boundaries, making it a valuable tool for systems with unknown structures. + +.. figure:: ../_static/demo_thumbnails/opengraph_demo_thumbnails/OGthumbnail_CERN_ANNNI.png + :align: center + :width: 70% + :target: javascript:void(0) + + Figure 1. Adding an ANNNI Ising to the cake + +The ANNNI model +-------------------------------------------------------------------------------- + +The ANNNI model describes a spin system with three types of competing interactions. Its Hamiltonian is given by + +.. math:: H = -J \sum_{i=1}^{N} \sigma_x^i\sigma_x^{i+1} - \kappa \sigma_x^{i}\sigma_x^{i+2} + h \sigma_z^i, \tag{1} + +where + +* :math:`\sigma_a^i` are the Pauli matrices acting on the :math:`i`-th spin (:math:`a \in \{x, y, z\}`), + +* :math:`J` is the nearest-neighbor coupling constant, which we set to :math:`1` without any loss of generality, + +* :math:`\kappa` controls the strength next-nearest-neighbor interaction, + +* and :math:`h` represents the the transverse magnetic field strength. + +Without loss of generality, we set :math:`J = 1` and consider open boundary conditions for positive :math:`\kappa` and :math:`h`. + +We start by importing the necessary libraries for simulation, optimization, and plotting our results, +as well as setting some important constants: +""" +import pennylane as qml +import numpy as np +from jax import jit, vmap, value_and_grad, random, config +from jax import numpy as jnp +import optax + +import matplotlib.pyplot as plt +from matplotlib.colors import ListedColormap, BoundaryNorm + +config.update("jax_enable_x64", True) + +seed = 123456 + +# Setting our constants +num_qubits = 8 # Number of spins in the Hamiltonian (= number of qubits) +side = 20 # Discretization of the Phase Diagram + +###################################################################### +# Implementing the Hamiltonian +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# In PennyLane, we can easily build the ANNNI's spin Hamiltonian following the same approach as the +# :doc:`demo on spin Hamiltonians `: + +def get_H(num_spins, k, h): + """Construction function the ANNNI Hamiltonian (J=1)""" + + # Interaction between spins (neighbouring): + H = -1 * (qml.PauliX(0) @ qml.PauliX(1)) + for i in range(1, num_spins - 1): + H = H - (qml.PauliX(i) @ qml.PauliX(i + 1)) + + # Interaction between spins (next-neighbouring): + for i in range(0, num_spins - 2): + H = H + k * (qml.PauliX(i) @ qml.PauliX(i + 2)) + + # Interaction of the spins with the magnetic field + for i in range(0, num_spins): + H = H - h * qml.PauliZ(i) + + return H + +###################################################################### +# Defining phase transition lines +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Due to the competition between the three types of interactions, the ANNNI model exhibits a rich and complex phase diagram. +# +# * The **Ising transition line** occurs at +# +# .. math:: h_I(\kappa) \approx \frac{1 - \kappa}{\kappa} \left(1 - \sqrt{\frac{1 - 3 \kappa + 4 \kappa^2 }{1 - \kappa}} \right),\tag{2} +# +# which separates the *ferromagnetic* phase from the *paramagnetic* phase. +# +# * The **Kosterlitz-Thouless (KT) transition line** occurs at +# +# .. math:: h_C(\kappa) \approx 1.05 \sqrt{(\kappa - 0.5) (\kappa - 0.1)}, \tag{3} +# +# which separates the *paramagnetic* phase from the *antiphase*. +# +# +# * Additionally, another phase transition has been numerically addressed but not yet confirmed. The **Berezinskii-Kosterlitz-Thouless (BKT) transition line** occurs at +# +# .. math:: h_{BKT} \approx 1.05 (\kappa - 0.5), \tag{4} +# +# which entirely lies within the *antiphase* region. +# + +def kt_transition(k): + """Kosterlitz-Thouless transition line""" + return 1.05 * np.sqrt((k - 0.5) * (k - 0.1)) + +def ising_transition(k): + """Ising transition line""" + return np.where(k == 0, 1, (1 - k) * (1 - np.sqrt((1 - 3 * k + 4 * k**2) / (1 - k))) / np.maximum(k, 1e-9)) + +def bkt_transition(k): + """Floating Phase transition line""" + return 1.05 * (k - 0.5) + +def get_phase(k, h): + """Get the phase from the DMRG transition lines""" + # If under the Ising Transition Line (Left side) + if k < .5 and h < ising_transition(k): + return 0 # Ferromagnetic + # If under the Kosterlitz-Thouless Transition Line (Right side) + + elif k > .5 and h < kt_transition(k): + return 1 # Antiphase + return 2 # else it is Paramagnetic + +###################################################################### +# .. figure:: ../_static/demonstration_assets/annni/annni_pd_L.png +# :align: center +# :width: 50% +# :target: javascript:void(0) +# +# State preparation +# ----------------- +# +# In this section, we prepare the ground states of the system, which will serve as inputs for both QML models. Several methods can be used, including **Variational Quantum Eigensolver (VQE)**, introduced in [#Peruzzo]_ and demonstrated in the :doc:`demo on VQE `, and **Matrix Product States (MPS)**, illustrated in the :doc:`demo on Constant-depth preparation of MPS with dynamic circuits `. +# For simplicity, in this demo, we compute the ground state directly by finding the *eigenvector* corresponding to the lowest eigenvalue of the Hamiltonian. The resulting states are then loaded into the quantum circuits using PennyLane’s :class:`~pennylane.StatePrep`. +# +# It is important to note that this approach is only feasible within the classically simulable regime, as it becomes quickly intractable for larger system sizes. +# +# Computing the ground states +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# We implement a function to compute the ground state by diagonalizing the Hamiltonian. + +def diagonalize_H(H_matrix): + """Returns the lowest eigenvector of the Hamiltonian matrix.""" + _, psi = jnp.linalg.eigh(H_matrix) # Compute eigenvalues and eigenvectors + return jnp.array(psi[:, 0], dtype=jnp.complex64) # Return the ground state + +# Create meshgrid of the parameter space +ks = np.linspace(0, 1, side) +hs = np.linspace(0, 2, side) +K, H = np.meshgrid(ks, hs) + +# Preallocate arrays for Hamiltonian matrices and phase labels. +H_matrices = np.empty((len(ks), len(hs), 2**num_qubits, 2**num_qubits)) +phases = np.empty((len(ks), len(hs)), dtype=int) + +for x, k in enumerate(ks): + for y, h in enumerate(hs): + H_matrices[y, x] = np.real(qml.matrix(get_H(num_qubits, k, h))) # Get Hamiltonian matrix + phases[y, x] = get_phase(k, h) # Get the respective phase given k and h + +# Vectorized diagonalization +psis = vmap(vmap(diagonalize_H))(H_matrices) + +###################################################################### +# Supervised learning of phases: QCNN +# ----------------------------------- +# +# QCNNs are a class of quantum circuits first introduced in [#Cong]_, inspired by their classical counterparts, Convolutional Neural Networks (CNNs). Like CNNs, QCNNs aim to learn representations from input data by leveraging its local properties. In this implementation, these local properties correspond to the interactions between neighbouring spins. +# +# A QCNN consists of two main components: +# +# * **Convolution layers**: alternating unitaries are applied to pairs of neighbouring spins. +# +# * **Pooling layers**: half of the qubits are measured, and based on the measurement outcome, different rotations are applied to the remaining qubits. +# +# For the output, we consider the model’s probability vector :math:`P(\kappa, h)` over the four computational basis states of the final two-qubit system, obtained using :func:`~pennylane.probs`. Each computational basis state is mapped to a specific phase as follows. +# +# * :math:`\vert 00 \rangle`: Ferromagnetic. +# * :math:`\vert 01 \rangle`: Antiphase. +# * :math:`\vert 10 \rangle`: Paramagnetic. +# * :math:`\vert 11 \rangle`: Trash class. +# +# Circuit definition +# ^^^^^^^^^^^^^^^^^^ +# +# We now define and implement the QCNN circuit. The `qcnn_ansatz` function builds the QCNN architecture by alternating convolution and pooling layers +# until only two qubits remain. The `qcnn_circuit` function embeds an input quantum state through :class:`~pennylane.StatePrep`, applies the QCNN ansatz, and returns +# the probability distribution over the final two-qubit system. Finally, we vectorize the circuit for efficient evaluation. + +def qcnn_ansatz(num_qubits, params): + """Ansatz of the QCNN model + Repetitions of the convolutional and pooling blocks + until only 2 wires are left unmeasured + """ + + # Convolution block + def conv(wires, params, index): + if len(wires) % 2 == 0: + groups = wires.reshape(-1, 2) + else: + groups = wires[:-1].reshape(-1, 2) + qml.RY(params[index], wires=int(wires[-1])) + index += 1 + + for group in groups: + qml.CNOT(wires=[int(group[0]), int(group[1])]) + for wire in group: + qml.RY(params[index], wires=int(wire)) + index += 1 + + return index + + # Pooiling block + def pool(wires, params, index): + # Process wires in pairs: measure one and conditionally rotate the other. + for wire_pool, wire in zip(wires[0::2], wires[1::2]): + m_0 = qml.measure(int(wire_pool)) + qml.cond(m_0 == 0, qml.RX)(params[index], wires=int(wire)) + qml.cond(m_0 == 1, qml.RX)(params[index + 1], wires=int(wire)) + index += 2 + # Remove the measured wire from active wires. + wires = np.delete(wires, np.where(wires == wire_pool)) + + # If an odd wire remains, apply a RX rotation. + if len(wires) % 2 != 0: + qml.RX(params[index], wires=int(wires[-1])) + index += 1 + + return index, wires + + # Initialize active wires and parameter index. + active_wires = np.arange(num_qubits) + index = 0 + + # Initial layer: apply RY to all wires. + for wire in active_wires: + qml.RY(params[index], wires=int(wire)) + index += 1 + + # Repeatedly apply convolution and pooling until there are 2 unmeasured wires + while len(active_wires) > 2: + # Convolution + index = conv(active_wires, params, index) + # Pooling + index, active_wires = pool(active_wires, params, index) + qml.Barrier() + + # Final layer: apply RY to the remaining active wires. + for wire in active_wires: + qml.RY(params[index], wires=int(wire)) + index += 1 + + return index, active_wires + +num_params, output_wires = qcnn_ansatz(num_qubits, [0]*100) + +@qml.qnode(qml.device("default.qubit", wires=num_qubits)) +def qcnn_circuit(params, state): + """QNode with QCNN ansatz and probabilities of unmeasured qubits as output""" + # Input ground state from diagonalization + qml.StatePrep(state, wires=range(num_qubits), normalize = True) + # QCNN + _, output_wires = qcnn_ansatz(num_qubits, params) + + return qml.probs([int(k) for k in output_wires]) + +# Vectorized circuit through vmap +vectorized_qcnn_circuit = vmap(jit(qcnn_circuit), in_axes=(None, 0)) + +# Draw the QCNN Architecture +fig,ax = qml.draw_mpl(qcnn_circuit)(np.arange(num_params), psis[0,0]) + +###################################################################### +# Training of the QCNN +# ^^^^^^^^^^^^^^^^^^^^ +# +# The training is performed by minimizing the **Cross Entropy loss** on the output probabilities +# +# .. math:: \mathcal{L} = -\frac{1}{|S|} \sum_{(\kappa, h) \in S} \sum_{j} y_j^{\frac1T}(\kappa, h) \log(p_j(\kappa, h))^\frac1T, \tag{5} +# +# where +# +# * :math:`S` is the training set, +# +# * :math:`p_j(\kappa, h)` is the model’s predicted probability of the system at :math:`(\kappa, h)` being in the :math:`j`-th phase, +# +# * :math:`y_j(\kappa, h)` represents the one-hot encoded labels for the three phases, +# +# * and :math:`T` is a temperature factor that controls the sharpness of the predicted probability distribution. + +def cross_entropy(pred, Y, T): + """Multi-class cross entropy loss function""" + epsilon = 1e-9 # Small value for numerical stability + pred = jnp.clip(pred, epsilon, 1 - epsilon) # Prevent log(0) + + # Apply sharpening (raise probabilities to the power of 1/T) + pred_sharpened = pred ** (1 / T) + pred_sharpened /= jnp.sum(pred_sharpened, axis=1, keepdims=True) # Re-normalize + + loss = -jnp.sum(Y * jnp.log(pred_sharpened), axis=1) + return jnp.mean(loss) + +###################################################################### +# The analytical points of the ANNNI model correspond to specific regions of the phase diagram where the system simplifies into two well-understood limits: +# +# * **Transverse-field Ising model** at :math:`\kappa = 0` in which we only have the magnetic field and the nearest-neighbor interactions. +# +# * **Quasi classical model** at :math:`h=0` in which we only have the nearest and next-nearest-neighbor interactions. +# +# For these points, we can derive the labels analytically which will then be used for the training of the QCNNs. + +# Mask for the analytical points +analytical_mask = (K == 0) | (H == 0) + +###################################################################### +# .. figure:: ../_static/demonstration_assets/annni/annni_pd_analytical_L.png +# :align: center +# :width: 50% +# :target: javascript:void(0) + +def train_qcnn(num_epochs, lr, T, seed): + """Training function of the QCNN architecture""" + + # Initialize PRNG key + key = random.PRNGKey(seed) + key, subkey = random.split(key) + + # Define the loss function + def loss_fun(params, X, Y): + preds = vectorized_qcnn_circuit(params, X) + return cross_entropy(preds, Y, T) + + # Consider only analytical points for the training + X_train, Y_train = psis[analytical_mask], phases[analytical_mask] + + # Convert labels to one-hot encoding + Y_train_onehot = jnp.eye(4)[Y_train] + + # Randomly initialize the parameters + params = random.normal(subkey, (num_params,)) + + # Initialize Adam optimizer + optimizer = optax.adam(learning_rate=lr) + optimizer_state = optimizer.init(params) + + loss_curve = [] + for epoch in range(num_epochs): + # Compute loss and gradients + loss, grads = value_and_grad(loss_fun)(params, X_train, Y_train_onehot) + + # Update parameters + updates, optimizer_state = optimizer.update(grads, optimizer_state) + params = optax.apply_updates(params, updates) + + loss_curve.append(loss) + + return params, loss_curve + +trained_params, loss_curve = train_qcnn(num_epochs=100, lr=5e-2, T=.5, seed=seed) + +# Plot the loss curve +plt.plot(loss_curve, label="Loss", color="blue", linewidth=2) +plt.xlabel("Epochs"), plt.ylabel("Cross-Entropy Loss") +plt.title("Figure 4. QCNN Training Cross-Entropy Loss Curve") +plt.legend() +plt.grid() +plt.show() + +###################################################################### +# After the training, much alike in [#Monaco]_ and [#Caro]_, we can inspect the model's generalization capabilities, +# by obtaining the predicted phase for every point across the 2D phase diagram and compare these predictions +# with the phase boundaries identified through density matrix renormalization group (DMRG) methods. + +# Take the predicted classes for each point in the phase diagram +predicted_classes = np.argmax( + vectorized_qcnn_circuit(trained_params, psis.reshape(-1, 2**num_qubits)), + axis=1 +) + +colors = ['#80bfff', '#fff2a8', '#80f090', '#da8080',] +phase_labels = ["Ferromagnetic", "Antiphase", "Paramagnetic", "Trash Class",] +cmap = ListedColormap(colors) + +bounds = [-0.5, 0.5, 1.5, 2.5, 3.5] +norm = BoundaryNorm(bounds, cmap.N) + +# Plot the predictions over the phase diagram +plt.figure(figsize=(4,4), constrained_layout=True) +plt.imshow( + predicted_classes.reshape(side, side), + cmap=cmap, + norm=norm, + aspect="auto", + origin="lower", + extent=[0, 1, 0, 2] +) + +# Plot the transition lines (Ising and KT) for reference. +k_vals1 = np.linspace(0.0, 0.5, 50) +k_vals2 = np.linspace(0.5, 1.0, 50) +plt.plot(k_vals1, ising_transition(k_vals1), 'k') +plt.plot(k_vals2, kt_transition(k_vals2), 'k') +plt.plot(k_vals2, bkt_transition(k_vals2), 'k', ls = '--') + +for color, phase in zip(colors, phase_labels[:-1]): + plt.scatter([], [], color=color, label=phase, edgecolors='black') +plt.plot([], [], 'k', label='Transition lines') + +plt.xlabel("k"), plt.ylabel("h") +plt.title("Figure 5. QCNN Classification") +plt.legend() +plt.show() + +###################################################################### +# Despite only being trained on the left and bottom borders of the phase diagram, the QCNN successfully generalizes across the entire diagram, aligning with the phase boundaries predicted by DMRG. +# +# In addition, [#Cea]_ presents an analysis of the QCNN’s performance as the number of qubits (hence the system's size) increases, showing that the overlap between trash class and the expected floating phase becomes more accurate. +# +# Unsupervised learning of phases: Quantum Anomaly Detection +# ---------------------------------------------------------- +# +# Quantum Anomaly Detection (QAD), introduced in [#Kottmann]_, is the quantum version of an autoencoder. However, unlike classical autoencoders, only the encoding (forward) process is trained here. This is because quantum operations are invertible, making a separate decoder unnecessary. +# +# In this method, we start with a single quantum state :math:`|\psi(\kappa, h)\rangle` taken from the ANNNI model. The goal is to optimize the parameters :math:`\theta` of a quantum circuit :math:`V(\theta)` so that it transforms the chosen input state into the following form: +# +# .. math:: V(\theta)|\psi(\kappa, h)\rangle = |\phi\rangle^{N-K} \otimes |0\rangle^{\otimes K}\tag{6} +# +# The equation above means that we are trying to find a unitary transformation that "compresses" the important information of the original state into a smaller number of qubits :math:`(N-K)`, while disentangling and resetting the remaining :math:`K` qubits to the trivial state :math:`|0\rangle`. In other words, the circuit learns to isolate the relevant features of the quantum state into :math:`|\phi\rangle`. +# +# Circuit definition +# ^^^^^^^^^^^^^^^^^^ +# +# We now define and implement the QAD circuit. The `anomaly_ansatz` function builds the QAD architecture. The `anomaly_circuit` function embeds an input quantum state through :class:`~pennylane.StatePrep`, applies the QAD ansatz, and returns the expectation values of the *trash qubits* used to evaluate the compression score. + +def anomaly_ansatz(n_qubit, params): + """Ansatz of the QAD model + Apply multi-qubit gates between trash and non-trash wires + """ + + # Block of gates connecting trash and non-trash wires + def block(nontrash, trash, shift): + # Connect trash wires + for i, wire in enumerate(trash): + target = trash[(i + 1 + shift) % len(trash)] + qml.CZ(wires=[int(wire), int(target)]) + # Connect each nontrash wire to a trash wire + for i, wire in enumerate(nontrash): + trash_idx = (i + shift) % len(trash) + qml.CNOT(wires=[int(wire), int(trash[trash_idx])]) + + depth = 2 # Number of repeated block layers + n_trashwire = n_qubit // 2 + + # Define trash wires as a contiguous block in the middle. + trash = np.arange(n_trashwire // 2, n_trashwire // 2 + n_trashwire) + nontrash = np.setdiff1d(np.arange(n_qubit), trash) + + index = 0 + + # Initial layer: apply RY rotations on all wires. + for wire in np.arange(n_qubit): + qml.RY(params[index], wires=int(wire)) + index += 1 + + # Repeatedly apply blocks of entangling gates and additional rotations. + for shift in range(depth): + block(nontrash, trash, shift) + qml.Barrier() + # In the final layer, only apply rotations on trash wires. + wires_to_rot = np.arange(n_qubit) if shift < depth - 1 else trash + for wire in wires_to_rot: + qml.RY(params[index], wires=int(wire)) + index += 1 + + return index, list(trash) + +num_anomaly_params, trash_wires = qcnn_ansatz(num_qubits, [0]*100) + +@qml.qnode(qml.device("default.qubit", wires=num_qubits)) +def anomaly_circuit(params, state): + """QNode with QAD ansatz and expectation values of the trash wires as output""" + # Input ground state from diagonalization + qml.StatePrep(state, wires=range(num_qubits), normalize = True) + # Quantum Anomaly Circuit + _, trash_wires = anomaly_ansatz(num_qubits, params) + + return [qml.expval(qml.PauliZ(int(k))) for k in trash_wires] + +# Vectorize the circuit using vmap +jitted_anomaly_circuit = jit(anomaly_circuit) +vectorized_anomaly_circuit = vmap(jitted_anomaly_circuit, in_axes=(None, 0)) + +# Draw the QAD Architecture +fig,ax = qml.draw_mpl(anomaly_circuit)(np.arange(num_anomaly_params), psis[0,0]) + +###################################################################### +# Training of the QAD +# ^^^^^^^^^^^^^^^^^^^ +# +# The training process for this architecture follows these steps: +# +# 1. **Selection of Training Event:** a single quantum state is selected as the training event. +# +# 2. **Compression Objective:** the training is performed to achieve the compression of the selected quantum state. This is done by minimizing the following loss function, known as the *compression score*: +# +# .. math:: \mathcal{C} = \frac{1}{2}\sum_{j\in q_T} (1-\left),\tag{7} +# +# where :math:`q_T` refers to the set of trash qubits, which make up :math:`N/2` of the total. +# By doing so, all the information of the input quantum state is compressed in the remaining non-measured qubits. +# +# In this case, the selected quantum state corresponds to the trivial case with :math:`\kappa = 0` and :math:`h = 0`. + +def train_anomaly(num_epochs, lr, seed): + """Training function of the QAD architecture""" + + # Initialize PRNG key + key = random.PRNGKey(seed) + key, subkey = random.split(key) + + # Define the loss function + def loss_fun(params, X): + # Output expectation values of the qubits + score = 1 - jnp.array(jitted_anomaly_circuit(params, X)) + loss_value = jnp.mean(score) + + return loss_value + + # Training set consists only of the k = 0 and h = 0 state + X_train = jnp.array(psis[0, 0]) + + # Randomly initialize parameters + params = random.normal(subkey, (num_anomaly_params,)) + + optimizer = optax.adam(learning_rate=lr) + optimizer_state = optimizer.init(params) + + loss_curve = [] + for epoch in range(num_epochs): + # Get random indices for a batch + key, subkey = random.split(key) + + # Compute loss and gradients + loss, grads = value_and_grad(loss_fun)(params, X_train) + + # Update parameters + updates, optimizer_state = optimizer.update(grads, optimizer_state) + params = optax.apply_updates(params, updates) + + loss_curve.append(loss) + + return params, loss_curve + +trained_anomaly_params, anomaly_loss_curve = train_anomaly(num_epochs=100, lr=1e-1, seed=seed) + +# Plot the loss curve + +plt.plot(anomaly_loss_curve, label="Loss", color="blue", linewidth=2) +plt.xlabel("Epochs"), plt.ylabel("Compression Loss") +plt.title("Figure 6. Anomaly training compression loss curve") +plt.legend() +plt.grid() +plt.show() + + +###################################################################### +# After training the circuit to optimally compress the (0,0) state, we evaluate the compression score for all other input states using the learned parameters. +# +# The model is expected to achieve near-optimal compression for states similar to the training state (namely those belonging to the same phase) resulting in a low compression score. Conversely, states from different phases should exhibit poorer compression, leading to a higher compression score. + +# Evaluate the compression score for each state in the phase diagram +compressions = vectorized_anomaly_circuit(trained_anomaly_params, psis.reshape(-1, 2**num_qubits)) +compressions = jnp.mean(1 - jnp.array(compressions), axis = 0) + +im = plt.imshow(compressions.reshape(side, side), aspect="auto", origin="lower", extent=[0, 1, 0, 2]) + +# Plot transition lines (assuming ising_transition and kt_transition are defined) +plt.plot(np.linspace(0.0, 0.5, 50), ising_transition(np.linspace(0.0, 0.5, 50)), 'k') +plt.plot(np.linspace(0.5, 1.0, 50), kt_transition(np.linspace(0.5, 1.0, 50)), 'k') + +plt.plot([], [], 'k', label='Transition Lines') +plt.scatter([0 +.3/len(ks)], [0 + .5/len(hs)], color='r', marker = 'x', label="Training point", s=50) + +plt.legend(), plt.xlabel("k"), plt.ylabel("h"), plt.title("Figure 7. Phase diagram with QAD") +cbar = plt.colorbar(im) +cbar.set_label(r"Compression Score $\mathcal{C}$") +plt.show() + +###################################################################### +# +# As expected, the compression score is nearly zero within the ferromagnetic phase, while higher scores are observed in the other regions. Surprisingly, the other regions display distinct compression scores that remain consistent within their respective areas. +# +# Using this model, we can clearly identify the three phases and their locations. Furthermore, as observed in [#Cea]_, by increasing the system size (around ~20 spins), a fourth phase emerges in the anticipated floating phase region. These regions are expected to become more sharply defined at larger system sizes, with more pronounced transitions between phases. +# +# Conclusion +# ----------- +# +# Quantum many-body systems present a compelling use case for QML models. In this tutorial, we explored two different approaches: +# +# * **Supervised learning**: the QCNN effectively generalizes phase classification beyond its training data, aligning with phase boundaries identified by classical methods. +# +# * **Unsupervised learning**: QAD successfully distinguishes different quantum phases without requiring labelled training data, making it ideal for systems with unknown structures. +# +# These techniques could be valuable for studying spin systems beyond classical simulability, particularly when non-local interactions are present and tensor network methods become inadequate. +# +# Acknowledgements +# ---------------- +# +# The author would like to acknowledge the contributions and support of Michele Grossi, Sofia Vallecorsa, Oriel Kiss and Enrique Rico Ortega from the CERN Quantum Technology Initiative, and Antonio Mandarino, and María Cea Fernández. The author also gratefully acknowledges the support of DESY (Deutsches Elektronen-Synchrotron). +# +# References +# ----------- +# +# .. [#Monaco] +# +# Saverio Monaco, Oriel Kiss, Antonio Mandarino, Sofia Vallecorsa, Michele Grossi, +# "Quantum phase detection generalization from marginal quantum neural network models", +# `PhysRevB.107.L081105 `__, 2023. +# +# .. [#Cea] +# +# Maria Cea, Michele Grossi, Saverio Monaco, Enrique Rico, Luca Tagliacozzo, Sofia Vallecorsa, +# "Exploring the Phase Diagram of the quantum one-dimensional ANNNI model", +# `arxiv:2402.11022 `__, 2023. +# +# .. [#Peruzzo] +# +# Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, Jeremy L. O'Brien, +# "A variational eigenvalue solver on a quantum processor", +# `Nat. Commun. 5, 4213 `__, 2013. +# +# .. [#Kottmann] +# +# Korbinian Kottmann, Friederike Metz, Joana Fraxanet, Niccolo Baldelli, +# "Variational Quantum Anomaly Detection: Unsupervised mapping of phase diagrams on a physical quantum computer", +# `PhysRevResearch.3.043184 `__, 2022. +# +# .. [#Cong] +# +# Iris Cong, Soonwon Choi, Mikhail D. Lukin, +# "Quantum Convolutional Neural Networks", +# `Nat. Phys. 15, 1273–1278 `__, 2019. +# +# .. [#Caro] +# +# Matthias C. Caro, Hsin-Yuan Huang, Marco Cerezo, Kunal Sharma, Andrew Sornborger, Lukasz Cincio, Patrick J. Coles, +# "Generalization in quantum machine learning from few training data", +# `Nat. Commun. 13, 4919 `__, 2022. +# +# About the author +# ---------------- +# \ No newline at end of file diff --git a/demonstrations_v2/tutorial_annni/metadata.json b/demonstrations_v2/tutorial_annni/metadata.json new file mode 100644 index 0000000000..4857fb7be0 --- /dev/null +++ b/demonstrations_v2/tutorial_annni/metadata.json @@ -0,0 +1,106 @@ +{ + "title": "Supervised and Unsupervised Quantum Machine Learning Models for the Phase Detection of the ANNNI Spin Model", + "authors": [ + { + "username": "SaverioMonaco" + } + ], + "dateOfPublication": "2025-04-28T00:00:00+00:00", + "dateOfLastModification": "2025-05-06T00:00:00+00:00", + "categories": [ + "Quantum Machine Learning" + ], + "tags": [], + "previewImages": [ + { + "type": "thumbnail", + "uri": "/_static/demo_thumbnails/regular_demo_thumbnails/thumbnail_CERN_ANNNI.png" + }, + { + "type": "large_thumbnail", + "uri": "/_static/demo_thumbnails/large_demo_thumbnails/thumbnail_large_CERN_ANNNI.png" + } + ], + "seoDescription": "Based on the paper Quantum phase detection generalization from marginal quantum neural network models, explore the phase diagram of the Axial Next-Nearest-Neighbor Ising (ANNNI) model using both supervised and unsupervised learning methods.", + "doi": "", + "references": [ + { + "id": "Monaco2023", + "type": "article", + "title": "Quantum phase detection generalization from marginal quantum neural network models.", + "authors": "Monaco S., Kiss O., Mandarino A., Vallecorsa S., Michele G.", + "year": "2023", + "journal": "PhysRevB", + "url": "https://journals.aps.org/prb/abstract/10.1103/PhysRevB.107.L081105" + }, + { + "id": "Cea2023", + "type": "preprint", + "title": "Exploring the Phase Diagram of the quantum one-dimensional ANNNI model", + "authors": "Maria Cea, Michele Grossi, Saverio Monaco, Enrique Rico, Luca Tagliacozzo, Sofia Vallecorsa", + "year": "2023", + "journal": "arXiv", + "url": "https://arxiv.org/abs/2402.11022" + }, + { + "id": "Peruzzo2013", + "type": "article", + "title": "A variational eigenvalue solver on a quantum processor", + "authors": "Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J. Love, Alán Aspuru-Guzik, Jeremy L. O'Brien", + "year": "2013", + "journal": "Nat. Commun.", + "url": "https://www.nature.com/articles/ncomms5213" + }, + { + "id": "Kottmann2022", + "type": "article", + "title": "Variational Quantum Anomaly Detection: Unsupervised mapping of phase diagrams on a physical quantum computer", + "authors": "Korbinian Kottmann, Friederike Metz, Joana Fraxanet, Niccolo Baldelli", + "year": "2022", + "journal": "PhysRevResearch", + "url": "https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.3.043184" + }, + { + "id": "Cong2019", + "type": "article", + "title": "Quantum Convolutional Neural Networks", + "authors": "Iris Cong, Soonwon Choi, Mikhail D. Lukin", + "year": "2019", + "journal": "Nat. Phys.", + "url": "https://www.nature.com/articles/s41567-019-0648-8" + }, + { + "id": "Caro2022", + "type": "article", + "title": "Generalization in quantum machine learning from few training data", + "authors": "Matthias C. Caro, Hsin-Yuan Huang, Marco Cerezo, Kunal Sharma, Andrew Sornborger, Lukasz Cincio, Patrick J. Coles", + "year": "2022", + "journal": "Nat. Commun.", + "url": "https://www.nature.com/articles/s41467-022-32550-3" + } + ], + "basedOnPapers": ["10.1103/PhysRevB.107.L081105"], + "referencedByPapers": [], + "relatedContent": [ + { + "type": "demonstration", + "id": "tutorial_how_to_build_spin_hamiltonians", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_constant_depth_mps_prep", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_vqe", + "weight": -1.0 + }, + { + "type": "demonstration", + "id": "tutorial_isingmodel_PyTorch", + "weight": 1.0 + } + ] +} diff --git a/demonstrations_v2/tutorial_classical_expval_estimation/demo.py b/demonstrations_v2/tutorial_classical_expval_estimation/demo.py index 886f40749a..f88b652873 100644 --- a/demonstrations_v2/tutorial_classical_expval_estimation/demo.py +++ b/demonstrations_v2/tutorial_classical_expval_estimation/demo.py @@ -1,6 +1,6 @@ r""" -Classically estimating expectation values from parametrized quantum circuits -============================================================================ +Pauli Propagation: Classically estimating expectation values from parametrized quantum circuits +=============================================================================================== In the race between classical and `quantum computing `__, an important question is whether there exist efficient classical algorithms to simulate quantum diff --git a/demonstrations_v2/tutorial_classical_expval_estimation/metadata.json b/demonstrations_v2/tutorial_classical_expval_estimation/metadata.json index 53e13148c5..cf2b0b6aae 100644 --- a/demonstrations_v2/tutorial_classical_expval_estimation/metadata.json +++ b/demonstrations_v2/tutorial_classical_expval_estimation/metadata.json @@ -1,12 +1,12 @@ { - "title": "Classically estimating expectation values from parametrized quantum circuits", + "title": "Pauli Propagation: Classically estimating expectation values from parametrized quantum circuits", "authors": [ { "username": "dwierichs" } ], "dateOfPublication": "2024-09-10T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2025-04-03T00:00:00+00:00", "categories": [ "Quantum Computing" ], diff --git a/demonstrations_v2/tutorial_classically_boosted_vqe/demo.py b/demonstrations_v2/tutorial_classically_boosted_vqe/demo.py index dd65f9a332..986ec0a7b2 100644 --- a/demonstrations_v2/tutorial_classically_boosted_vqe/demo.py +++ b/demonstrations_v2/tutorial_classically_boosted_vqe/demo.py @@ -449,7 +449,7 @@ def hadamard_test(Uq, Ucl, component="real"): qml.Hadamard(wires=[0]) qml.ControlledQubitUnitary( - Uq.conjugate().T @ Ucl, control_wires=[0], wires=wires[1:] + Uq.conjugate().T @ Ucl, wires=wires ) qml.Hadamard(wires=[0]) diff --git a/demonstrations_v2/tutorial_classically_boosted_vqe/metadata.json b/demonstrations_v2/tutorial_classically_boosted_vqe/metadata.json index 797d86cf63..df8e6990f6 100644 --- a/demonstrations_v2/tutorial_classically_boosted_vqe/metadata.json +++ b/demonstrations_v2/tutorial_classically_boosted_vqe/metadata.json @@ -9,7 +9,7 @@ } ], "dateOfPublication": "2022-10-31T00:00:00+00:00", - "dateOfLastModification": "2024-12-13T00:00:00+00:00", + "dateOfLastModification": "2025-01-23T00:00:00+00:00", "categories": [ "Quantum Chemistry" ], diff --git a/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/demo.py b/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/demo.py index a4f5439a83..cb7ff2379b 100644 --- a/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/demo.py +++ b/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/demo.py @@ -20,7 +20,7 @@ Introduction ------------ -The :doc:`KAK theorem ` is an important result from Lie theory that states that any Lie group element :math:`U` can be decomposed +The :doc:`KAK decomposition ` is an important result from Lie theory that states that any Lie group element :math:`U` can be decomposed as :math:`U = K_1 A K_2,` where :math:`K_{1, 2}` and :math:`A` are elements of two special sub-groups :math:`\mathcal{K}` and :math:`\mathcal{A},` respectively. In special cases, the decomposition simplifies to :math:`U = K A K^\dagger.` @@ -36,14 +36,14 @@ We can use this general result from Lie theory as a powerful circuit decomposition technique. .. note:: We recommend a basic understanding of Lie algebras, see e.g. our :doc:`introduction to (dynamical) Lie algebras for quantum practitioners `. - Otherwise, this demo should be self-contained, though for the mathematically inclined, we further recommend our :doc:`demo on the KAK theorem ` - that dives into the mathematical depths of the theorem and provides more background info. + Otherwise, this demo should be self-contained, though for the mathematically inclined, we further recommend our :doc:`demo on the KAK decomposition ` + that dives into the mathematical depths of the decomposition and provides more background info. Goal: Fast-forwarding time evolutions using the KAK decomposition ----------------------------------------------------------------- Unitary gates in quantum computing are described by the special unitary Lie group :math:`SU(2^n),` so we can use the KAK -theorem to decompose quantum gates into :math:`U = K_1 A K_2.` While the mathematical statement is rather straightforward, +decomposition to factorize quantum gates into :math:`U = K_1 A K_2.` While the mathematical statement is rather straightforward, actually finding this decomposition is not. We are going to follow the recipe prescribed in `Fixed Depth Hamiltonian Simulation via Cartan Decomposition `__ [#Kökcü]_, which tackles this decomposition on the level of the associated Lie algebra via Cartan decomposition. @@ -74,6 +74,7 @@ import numpy as np import pennylane as qml from pennylane import X, Y, Z +from pennylane.liealg import even_odd_involution, cartan_decomp, horizontal_cartan_subalgebra import jax import jax.numpy as jnp @@ -111,10 +112,8 @@ # One common choice of involution is the so-called even-odd involution for Pauli words, # :math:`P = P_1 \otimes P_2 .. \otimes P_n,` where :math:`P_j \in \{I, X, Y, Z\}.` # It essentially counts whether the number of non-identity Pauli operators in the Pauli word is even or odd. - -def even_odd_involution(op): - [pw] = op.pauli_rep - return len(pw) % 2 +# It is readily available in PennyLane as :func:`~.pennylane.liealg.even_odd_involution`, which +# we already imported above. even_odd_involution(X(0)), even_odd_involution(X(0) @ Y(3)) @@ -126,30 +125,9 @@ def even_odd_involution(op): # sort the operators by whether or not they yield a plus or minus sign from the involution function. # This is possible because the operators and involution nicely align with the eigenspace decomposition. -def cartan_decomposition(g, involution): - """Cartan Decomposition g = k + m - - Args: - g (List[PauliSentence]): the (dynamical) Lie algebra to decompose - involution (callable): Involution function :math:`\Theta(\cdot)` to act on PauliSentence ops, should return ``0/1`` or ``True/False``. - - Returns: - k (List[PauliSentence]): the vertical subspace :math:`\Theta(x) = x` - m (List[PauliSentence]): the horizontal subspace :math:`\Theta(x) = -x` """ - m = [] - k = [] - - for op in g: - if involution(op): # vertical space when involution returns True - k.append(op) - else: # horizontal space when involution returns False - m.append(op) - return k, m - -k, m = cartan_decomposition(g, even_odd_involution) +k, m = cartan_decomp(g, even_odd_involution) len(g), len(k), len(m) - ############################################################################## # We have successfully decomposed the 60-dimensional Lie algebra # into a 24-dimensional vertical subspace and a 36-dimensional subspace. @@ -187,51 +165,11 @@ def cartan_decomposition(g, involution): # that commute with it. # # We then obtain a further split of the vector space :math:`\mathfrak{m} = \tilde{\mathfrak{m}} \oplus \mathfrak{h},` -# where :math:`\tilde{\mathfrak{m}}` is just the remainder of :math:`\mathfrak{m}.` +# where :math:`\tilde{\mathfrak{m}}` is just the remainder of :math:`\mathfrak{m}.` The function +# :func:`~.pennylane.liealg.horizontal_cartan_subalgebra` returns some additional information, which we will +# not use here. -def _commutes_with_all(candidate, ops): - r"""Check if ``candidate`` commutes with all ``ops``""" - for op in ops: - com = candidate.commutator(op) - com.simplify() - - if not len(com) == 0: - return False - return True - -def cartan_subalgebra(m, which=0): - """Compute the Cartan subalgebra from the horizontal subspace :math:`\mathfrak{m}` - of the Cartan decomposition - - This implementation is specific for cases of bases of m with pure Pauli words as - detailed in Appendix C in `2104.00728 `__. - - Args: - m (List[PauliSentence]): the horizontal subspace :math:`\Theta(x) = -x - which (int): Choice for the initial element of m from which to construct - the maximal Abelian subalgebra - - Returns: - mtilde (List): remaining elements of :math:`\mathfrak{m}` - s.t. :math:`\mathfrak{m} = \tilde{\mathfrak{m}} \oplus \mathfrak{h}`. - h (List): Cartan subalgebra :math:`\mathfrak{h}`. - - """ - - h = [m[which]] # first candidate - mtilde = m.copy() - - for m_i in m: - if _commutes_with_all(m_i, h): - if m_i not in h: - h.append(m_i) - - for h_i in h: - mtilde.remove(h_i) - - return mtilde, h - -mtilde, h = cartan_subalgebra(m) +g, k, mtilde, h, _ = horizontal_cartan_subalgebra(k, m, tol=1e-8) len(g), len(k), len(mtilde), len(h) ############################################################################## @@ -241,7 +179,7 @@ def cartan_subalgebra(m, which=0): # Variational KhK decomposition # ----------------------------- # -# The KAK theorem is not constructive in the sense that it proves that there exists such a decomposition, but there is no general way of obtaining +# The KAK decomposition is not constructive in the sense that it proves that there exists such a decomposition, but there is no general way of obtaining # it. In particular, there are no linear algebra subroutines implemented in ``numpy`` or ``scipy`` that just compute it for us. # Here, we follow the construction of [#Kökcü]_ for the special case of :math:`H` being in the horizontal space and the decomposition # simplifying to :math:`H = K^\dagger h K`. @@ -282,10 +220,10 @@ def cartan_subalgebra(m, which=0): # def run_opt( - value_and_grad, + loss, theta, - n_epochs=500, - lr=0.1, + n_epochs=1000, + lr=0.05, ): """Boilerplate JAX optimization""" value_and_grad = jax.jit(jax.value_and_grad(loss)) @@ -336,7 +274,7 @@ def loss(theta): A = K_m @ v_m @ K_m.conj().T return jnp.trace(A.conj().T @ H_m).real -theta0 = jnp.ones(len(k), dtype=float) +theta0 = jax.random.normal(jax.random.PRNGKey(0), shape=(len(k),), dtype=float) thetas, energy, _ = run_opt(loss, theta0, n_epochs=1000, lr=0.05) plt.plot(energy - np.min(energy)) @@ -359,13 +297,18 @@ def loss(theta): # .. math:: h_0 = K_c^\dagger H K_c. h_0_m = Kc_m.conj().T @ H_m @ Kc_m -h_0 = qml.pauli_decompose(h_0_m) -print(len(h_0)) +# decompose h_0_m in terms of the basis of h +basis = [qml.matrix(op, wire_order=range(n_wires)) for op in h] +coeffs = qml.pauli.trace_inner_product(h_0_m, basis) + +# ensure that decomposition is correct, i.e. h_0_m is truely an element of just h +h_0_m_recomposed = np.sum([c * op for c, op in zip(coeffs, basis)], axis=0) +print("Decomposition of h_0 is faithful: ", np.allclose(h_0_m_recomposed, h_0_m, atol=1e-10)) + +# sanity check that the horizontal CSA is Abelian, i.e. all its elements commute +print("All elements in h commute with each other: ", qml.liealg.check_abelian(h)) -# assure that h_0 is in \mathfrak{h} -h_vspace = qml.pauli.PauliVSpace(h) -not h_vspace.is_independent(h_0.pauli_rep) ############################################################################## # @@ -393,6 +336,7 @@ def loss(theta): t = 1. U_exact = qml.exp(-1j * t * H) U_exact_m = qml.matrix(U_exact, wire_order=range(n_wires)) +h_0 = qml.dot(coeffs, h) def U_kak(theta_opt, t): qml.adjoint(K)(theta_opt, k) @@ -478,7 +422,7 @@ def compute_res(Us): # Conclusion # ---------- # -# The KAK theorem is a very general mathematical result with far-reaching consequences. +# The KAK decomposition is a very general mathematical result with far-reaching consequences. # While there is no canonical way of obtaining an actual decomposition in practice, we followed # the approach of [#Kökcü]_ which uses a specifically designed loss function and variational # optimization to find the decomposition. diff --git a/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/metadata.json b/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/metadata.json index 1efc93770b..a78b88eea4 100644 --- a/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/metadata.json +++ b/demonstrations_v2/tutorial_fixed_depth_hamiltonian_simulation_via_cartan_decomposition/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2024-12-19T00:00:00+00:00", - "dateOfLastModification": "2025-01-10T00:00:00+00:00", + "dateOfLastModification": "2025-04-11T00:00:00+00:00", "categories": [ "Quantum Computing", "Algorithms" diff --git a/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/demo.py b/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/demo.py new file mode 100644 index 0000000000..8a32264830 --- /dev/null +++ b/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/demo.py @@ -0,0 +1,391 @@ +r"""How to build compressed double-factorized Hamiltonians +========================================================== + +Compressed double factorization offers a powerful approach to overcome key limitations in +quantum chemistry simulations. Specifically, it tackles the runtime's dependency on the +Hamiltonian's one-norm and the shot requirements linked to the number of terms [#cdf]_. +In this tutorial, we will learn how to construct the electronic Hamiltonian in the compressed +double-factorized form using tensor contractions. We will also show how this technique allows +having a linear combination of unitaries (LCU) representation suitable for error-corrected +algorithms, which facilitates efficient simulations via linear-depth circuits with +`Givens rotations `_. + +.. figure:: ../_static/demo_thumbnails/opengraph_demo_thumbnails/OGthumbnail_how_to_build_cdf_hamiltonians.png + :align: center + :width: 70% + :target: javascript:void(0) + +Constructing the electronic Hamiltonian +---------------------------------------- + +The Hamiltonian of a molecular system in the second-quantized form can be expressed as a +sum of the one-body and two-body terms as follows: + +.. math:: H = \mu + \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac{1}{2} \sum_{\sigma \tau, pqrs} g_{pqrs} a^\dagger_{\sigma, p} a^\dagger_{\tau, q} a_{\tau, r} a_{\sigma, s}, + +where the tensors :math:`h_{pq}` and :math:`g_{pqrs}` are the one- and two-body integrals, +:math:`a^\dagger` and :math:`a` are the creation and annihilation operators, :math:`\mu` is the +nuclear repulsion energy constant, :math:`\sigma \in {\uparrow, \downarrow}` represents the spin, +and :math:`p, q, r, s` are the orbital indices. In PennyLane, we can obtain :math:`\mu`, +:math:`h_{pq}` and :math:`g_{pqrs}` using the :func:`~pennylane.qchem.electron_integrals` function: +""" + +import pennylane as qml + +symbols = ["H", "H", "H", "H"] +geometry = qml.math.array([[0., 0., -0.2], [0., 0., -0.1], [0., 0., 0.1], [0., 0., 0.2]]) + +mol = qml.qchem.Molecule(symbols, geometry) +nuc_core, one_body, two_body = qml.qchem.electron_integrals(mol)() + +print(f"One-body and two-body tensor shapes: {one_body.shape}, {two_body.shape}") + +###################################################################### +# In the above expression of :math:`H`, the two-body tensor :math:`g_{pqrs}` +# can be rearranged to define :math:`V_{pqrs}` in the `chemist notation +# `_, +# which leads to a one-body offset term :math:`\sum_{s} V_{pssq}`. This +# allows us to rewrite the Hamiltonian as: +# +# .. math:: H_{\text{C}} = \mu + \sum_{\sigma \in {\uparrow, \downarrow}} \sum_{pq} T_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \sum_{\sigma, \tau \in {\uparrow, \downarrow}} \sum_{pqrs} V_{pqrs} a^\dagger_{\sigma, p} a_{\sigma, q} a^\dagger_{\tau, r} a_{\tau, s}, +# +# with the transformed one-body terms :math:`T_{pq} = h_{pq} - 0.5 \sum_{s} g_{pqss}`. +# We can obtain the :math:`V_{pqrs}` and :math:`T_{pq}` tensors as: +# + +two_chem = 0.5 * qml.math.swapaxes(two_body, 1, 3) # V_pqrs +one_chem = one_body - 0.5 * qml.math.einsum("pqss", two_body) # T_pq + +###################################################################### +# A key feature of this representation is that the modified two-body terms can be factorized +# into a sum of low-rank terms, which can be used to efficiently simulate the Hamiltonian [#cdf2]_. +# We will see how to do this with the double-factorization methods in the next section. +# +# Double factorizing the Hamiltonian +# ----------------------------------- +# +# The double factorization of a Hamiltonian can be described as a Hamiltonian manipulation +# technique based on decomposing :math:`V_{pqrs}` into symmetric tensors :math:`L^{(t)}` +# called *factors*, such that :math:`V_{pqrs} = \sum_t^T L_{pq}^{(t)} L_{rs}^{(t) {\dagger}}` +# and the rank :math:`T \leq N^2`, where :math:`N` is the number of orbitals. We can do this +# by performing an eigenvalue or a pivoted Cholesky decomposition of the modified two-body +# tensor. Moreover, each of the :math:`L^{(t)}` can be further eigendecomposed as +# :math:`L^{(t)}_{pq} = \sum_{i} U_{pi}^{(t)} W_i^{(t)} U_{qi}^{(t)}` to perform a second +# tensor factorization. This enables us to express the two-body tensor :math:`V_{pqrs}` in +# the following double-factorized form in terms of orthonormal core tensors :math:`Z^{(t)}` +# and symmetric leaf tensors :math:`U^{(t)}` [#cdf2]_: +# +# .. math:: V_{pqrs} \approx \sum_t^T \sum_{ij} U_{pi}^{(t)} U_{qi}^{(t)} Z_{ij}^{(t)} U_{rj}^{(t)} U_{sj}^{(t)}, +# +# where :math:`Z_{ij}^{(t)} = W_i^{(t)} W_j^{(t)}`. This decomposition is referred +# to as the *explicit* double factorization (XDF) and decreases the number of terms +# in the qubit basis from :math:`O(N^4)` to :math:`O(N^3)`, assuming the rank of the +# second tensor factorization to be :math:`O(N)`. In PennyLane, this can be done using the +# :func:`~pennylane.qchem.factorize` function, where one can choose the decomposition method for +# the first tensor factorization (``cholesky``), truncate the resulting factors by discarding +# the ones with individual contributions below a specified threshold (``tol_factor``), and +# optionally control the ranks of their second factorization (``tol_eigval``) as shown below: +# + +factors, _, _ = qml.qchem.factorize(two_chem, cholesky=True, tol_factor=1e-5) +print("Shape of the factors: ", factors.shape) + +approx_two_chem = qml.math.tensordot(factors, factors, axes=([0], [0])) +assert qml.math.allclose(approx_two_chem, two_chem, atol=1e-5) + +###################################################################### +# Performing block-invariant symmetry shift +# ------------------------------------------ +# +# We can further improve the double-factorization by employing the block-invariant +# symmetry shift (BLISS) technique, which modifies the Hamiltonian's action on the +# undesired electron-number subspace [#bliss]_. It helps decrease the one-norm and +# the spectral range of the Hamiltonian. In PennyLane, this symmetry shift can be +# done using the :func:`~pennylane.qchem.symmetry_shift` function, which returns +# the symmetry-shifted integral tensors and core constant: +# + +core_shift, one_shift, two_shift = qml.qchem.symmetry_shift( + nuc_core, one_chem, two_chem, n_elec = mol.n_electrons +) # symmetry-shifted terms of the Hamiltonian + +###################################################################### +# Then we can use these shifted terms to obtain a double-factorized representation of +# the Hamiltonian that has a lower one-norm than the original one. For instance, we can +# compare the improvement in the one-norm of the shifted Hamiltonian over the original one +# by accessing the :class:`~.pennylane.resource.DoubleFactorization`'s ``lamb`` attribute: +# + +from pennylane.resource import DoubleFactorization as DF + +DF_chem_norm = DF(one_chem, two_chem, chemist_notation=True).lamb +DF_shift_norm = DF(one_shift, two_shift, chemist_notation=True).lamb +print(f"Decrease in one-norm: {DF_chem_norm - DF_shift_norm}") + +###################################################################### +# Compressing the double factorization +# ------------------------------------- +# +# In many practical scenarios, the double factorization method can be further optimized by +# performing a numerical tensor-fitting of the decomposed two-body terms to obtain :math:`V^\prime` +# such that the approximation error :math:`||V - V^\prime||` remains below a desired threshold +# [#cdf]_. This is referred to as the *compressed* double factorization (CDF) as it reduces the +# number of terms in the factorization of the two-body term from :math:`O(N^3)` to :math:`O(N)` +# and achieves lower approximation errors than the truncated XDF. Compression can be done by +# beginning with :math:`O(N)` random core and leaf tensors and optimizing them based on the +# following cost function :math:`\mathcal{L}` in a greedy manner: +# +# .. math:: \mathcal{L}(U, Z) = \frac{1}{2} \bigg|V_{pqrs} - \sum_t^T \sum_{ij} U_{pi}^{(t)} U_{pj}^{(t)} Z_{ij}^{(t)} U_{qk}^{(t)} U_{ql}^{(t)}\bigg|_{\text{F}} + \rho \sum_t^T \sum_{ij} \bigg|Z_{ij}^{(t)}\bigg|^{\gamma}, +# +# where :math:`|\cdot|_{\text{F}}` denotes the Frobenius norm, :math:`\rho` is a constant +# scaling factor, and :math:`|\cdot|^\gamma` specifies the optional L1 and L2 `regularization +# `_ +# that improves the energy variance of the resulting representation. In PennyLane, this +# compression can be done by using the ``compressed=True`` keyword argument in the +# :func:`~pennylane.qchem.factorize` function. The regularization term will be included +# if the ``regularization`` keyword argument is set to either ``"L1"`` or ``"L2"``: +# + +_, two_body_cores, two_body_leaves = qml.qchem.factorize( + two_shift, tol_factor=1e-2, cholesky=True, compressed=True, regularization="L2" +) # compressed double-factorized shifted two-body terms with "L2" regularization +print(f"Two-body tensors' shape: {two_body_cores.shape, two_body_leaves.shape}") + +approx_two_shift = qml.math.einsum( + "tpk,tqk,tkl,trl,tsl->pqrs", + two_body_leaves, two_body_leaves, two_body_cores, two_body_leaves, two_body_leaves +) # computing V^\prime and comparing it with V below +assert qml.math.allclose(approx_two_shift, two_shift, atol=1e-2) + +###################################################################### +# While the previous shape output for the factors ``(10, 4, 4)`` meant we had :math:`10` two-body +# terms in our factorization, the current shape output ``(6, 4, 4)`` indicates that we have +# :math:`6` terms. This means that the number of terms in the factorization has decreased almost +# by half, which is quite significant! +# +# Constructing the double-factorized Hamiltonian +# ----------------------------------------------- +# +# We can eigendecompose the one-body tensor to obtain similar orthonormal :math:`U^{(0)}` and +# symmetric :math:`Z^{(0)}` tensors for the one-body term and use the compressed factorization +# of the two-body term described in the previous section to express the Hamiltonian in the +# double-factorized form as: +# +# .. math:: H_{\text{CDF}} = \mu + \sum_{\sigma \in {\uparrow, \downarrow}} U^{(0)}_{\sigma} \left( \sum_{p} Z^{(0)}_{p} a^\dagger_{\sigma, p} a_{\sigma, p} \right) U_{\sigma}^{(0)\ \dagger} + \sum_t^T \sum_{\sigma, \tau \in {\uparrow, \downarrow}} U_{\sigma, \tau}^{(t)} \left( \sum_{pq} Z_{pq}^{(t)} a^\dagger_{\sigma, p} a_{\sigma, p} a^\dagger_{\tau, q} a_{\tau, q} \right) U_{\sigma, \tau}^{(t)\ \dagger}. +# +# This Hamiltonian can be easily mapped to the qubit basis via `Jordan-Wigner +# transformation `_ (JWT) using +# :math:`a_p^\dagger a_p = n_p \mapsto 0.5 * (1 - z_p)`, where :math:`n_p` is the number +# operator and :math:`z_p` is the Pauli-Z operation acting on the qubit corresponding to +# orbital :math:`p`. The mapped form naturally gives rise to a measurement grouping, where +# the terms within the basis transformation :math:`U^{(i)}` can be measured simultaneously. +# These can be obtained with the :func:`~pennylane.qchem.basis_rotation` function, which +# performs the double-factorization and JWT automatically. +# +# Another advantage of the double-factorized form is the efficient simulation of the Hamiltonian +# evolution. Before discussing it in the next section, we note that mapping a two-body term to +# the qubit basis will result in two additional one-qubit Pauli-Z terms. We can simplify their +# simulation circuits by accounting for these additional terms directly in the one-body tensor +# using a correction ``one_body_extra``. We can then decompose the corrected one-body terms into +# the orthonormal :math:`U^{\prime(0)}` and symmetric :math:`Z^{\prime(0)}` tensors instead: +# + +two_core_prime = (qml.math.eye(mol.n_orbitals) * two_body_cores.sum(axis=-1)[:, None, :]) +one_body_extra = qml.math.einsum( + 'tpk,tkk,tqk->pq', two_body_leaves, two_core_prime, two_body_leaves +) # one-body correction + +# factorize the corrected one-body tensor to obtain the core and leaf tensors +one_body_eigvals, one_body_eigvecs = qml.math.linalg.eigh(one_shift + one_body_extra) +one_body_cores = qml.math.expand_dims(qml.math.diag(one_body_eigvals), axis=0) +one_body_leaves = qml.math.expand_dims(one_body_eigvecs, axis=0) + +print(f"One-body tensors' shape: {one_body_cores.shape, one_body_leaves.shape}") + +###################################################################### +# We can now specify the Hamiltonian programmatically in the (compressed) +# double-factorized form as a ``dict`` object with the following three keys: +# ``nuc_constant`` (:math:`\mu`), +# ``core_tensors`` (:math:`\left[ Z^{\prime(0)}, Z^{(t_1)}, \ldots, Z^{(t_T)} \right]`), and +# ``leaf_tensors`` (:math:`\left[ U^{\prime(0)}, U^{(t_1)}, \ldots, U^{(t_T)} \right]`): +# + +cdf_hamiltonian = { + "nuc_constant": core_shift[0], + "core_tensors": qml.math.concatenate((one_body_cores, two_body_cores), axis=0), + "leaf_tensors": qml.math.concatenate((one_body_leaves, two_body_leaves), axis=0), +} # CDF Hamiltonian + +###################################################################### +# Simulating the double-factorized Hamiltonian +# --------------------------------------------- +# +# To simulate the time evolution of the CDF Hamiltonian, +# we will first need to learn how to apply the unitary operations +# represented by the exponentiated leaf and core tensors. The former can be done using the +# :class:`~.pennylane.BasisRotation` operation, which implements the unitary transformation +# :math:`\exp \left( \sum_{pq}[\log U]_{pq} (a^\dagger_p a_q - a^\dagger_q a_p) \right)` +# using the `Givens rotation networks +# `_ +# that can be efficiently implemented on quantum hardware. The ``leaf_unitary_rotation`` +# function below does this for a leaf tensor: +# + +def leaf_unitary_rotation(leaf, wires): + """Applies the basis rotation transformation corresponding to the leaf tensor.""" + basis_mat = qml.math.kron(leaf, qml.math.eye(2)) # account for spin + return qml.BasisRotation(unitary_matrix=basis_mat, wires=wires) + +###################################################################### +# Similarly, the unitary transformation for the core tensors can be applied efficiently +# via the ``core_unitary_rotation`` function defined below. The function uses the +# :class:`~.pennylane.RZ` and :class:`~.pennylane.IsingZZ` gates for implementing +# the diagonal and entangling phase rotations for the one- and two-body core tensors, +# respectively, and :class:`~.pennylane.GlobalPhase` for the corresponding global phases: +# + +import itertools as it + +def core_unitary_rotation(core, body_type, wires): + """Applies the unitary transformation corresponding to the core tensor.""" + ops = [] + if body_type == "one_body": # implements one-body term + for wire, cval in enumerate(qml.math.diag(core)): + for sigma in [0, 1]: + ops.append(qml.RZ(-cval, wires=2 * wire + sigma)) + ops.append(qml.GlobalPhase(qml.math.sum(core), wires=wires)) + + if body_type == "two_body": # implements two-body term + for odx1, odx2 in it.product(range(len(wires) // 2), repeat=2): + cval = core[odx1, odx2] / 4.0 + for sigma, tau in it.product(range(2), repeat=2): + if odx1 != odx2 or sigma != tau: + ops.append(qml.IsingZZ(cval, wires=[2*odx1+sigma, 2*odx2+tau])) + gphase = 0.5 * qml.math.sum(core) - 0.25 * qml.math.trace(core) + ops.append(qml.GlobalPhase(-gphase, wires=wires)) + return ops + +###################################################################### +# We can now use these functions to approximate the evolution operator :math:`e^{-iHt}` for +# a time :math:`t` with the Suzuki-Trotter product formula, which uses symmetrized products +# :math:`S_m` defined for an order :math:`m \in [1, 2, 4, \ldots, 2k \in \mathbb{N}]` +# and repeated multiple times [#trotter]_. In general, this can be easily implemented for +# standard non-factorized Hamiltonians using the :class:`~.pennylane.TrotterProduct` operation, +# which defines these products recursively, leading to an exponential scaling in its complexity +# with the number of terms in the Hamiltonian and making it inefficient for larger system sizes. +# +# The exponential scaling can be improved to a great extent by working with the compressed +# double-factorized form of the Hamiltonian as it allows reducing the number of terms in the +# Hamiltonian from :math:`O(N^4)` to :math:`O(N)`. While doing this is not directly supported +# in PennyLane in the form of a template, we can still implement the first-order Trotter +# step using the following :func:`CDFTrotterStep` function that uses the CDF Hamiltonian +# with the ``leaf_unitary_rotation`` and ``core_unitary_rotation`` functions defined earlier. +# We can then use the :func:`~.pennylane.trotterize` function to implement any higher-order +# Suzuki-Trotter products. +# + +def CDFTrotterStep(time, cdf_ham, wires): + """Implements a first-order Trotter step for a CDF Hamiltonian. + + Args: + time (float): time-step for a Trotter step. + cdf_ham (dict): dictionary describing the CDF Hamiltonian. + wires (list): list of integers representing the qubits. + """ + cores, leaves = cdf_ham["core_tensors"], cdf_ham["leaf_tensors"] + for bidx, (core, leaf) in enumerate(zip(cores, leaves)): + # Note: only the first term is one-body, others are two-body + body_type = "two_body" if bidx else "one_body" + qml.prod( + # revert the change-of-basis for leaf tensor + leaf_unitary_rotation(leaf.conjugate().T, wires), + # apply the rotation for core tensor scaled by the time-step + *core_unitary_rotation(time * core, body_type, wires), + # apply the basis rotation for leaf tensor + leaf_unitary_rotation(leaf, wires), + ) # Note: prod applies operations in the reverse order (right-to-left). + +###################################################################### +# We now use this function to simulate the evolution of the :math:`H_4` Hamiltonian +# described in the compressed double-factorized form for a given number of +# steps ``n_steps`` and starting from a Hartree-Fock state ``hf_state`` +# with the following circuit that applies a second-order Trotter product: +# + +time, circ_wires = 1.0, range(2 * mol.n_orbitals) +hf_state = qml.qchem.hf_state(electrons=mol.n_electrons, orbitals=len(circ_wires)) + +@qml.qnode(qml.device("lightning.qubit", wires=circ_wires)) +def cdf_circuit(n_steps, order): + qml.BasisState(hf_state, wires=circ_wires) + qml.trotterize(CDFTrotterStep, n_steps, order)(time, cdf_hamiltonian, circ_wires) + qml.GlobalPhase(cdf_hamiltonian["nuc_constant"], wires=circ_wires) + return qml.state() + +circuit_state = cdf_circuit(n_steps=10, order=2) + +###################################################################### +# We now test the accuracy of the Hamiltonian simulation by evolving the +# Hartree-Fock state analytically ourselves and testing the fidelity +# of the ``evolved_state`` with the ``circuit_state``: +# + +from pennylane.math import fidelity_statevector +from scipy.linalg import expm + +# Evolve the state vector |0...0> to the |HF> state of the system +init_state = qml.math.array([1] + [0] * (2**len(circ_wires) - 1)) +hf_state_vec = qml.matrix(qml.BasisState(hf_state, wires=circ_wires)) @ init_state + +H = qml.qchem.molecular_hamiltonian(mol)[0] # original Hamiltonian +evolved_state = expm(-1j * qml.matrix(H) * time) @ hf_state_vec # e^{-iHt} @ |HF> + +print(f"Fidelity of two states: {fidelity_statevector(circuit_state, evolved_state)}") + +###################################################################### +# As we can see, the fidelity of the evolved state from the circuit is close to +# :math:`1.0`, which indicates that the evolution of the CDF Hamiltonian sufficiently +# matches that of the original one. +# +# Conclusion +# ----------- +# +# Compressed double-factorized representation for the Hamiltonians serves three key purposes. +# First, it provides for a more compact representation of the Hamiltonian that can be stored +# and manipulated easier. Second, it allows more efficient simulations of the Hamiltonian +# time evolution because the number of terms is reduced quadratically. Third, the compact +# representation can be further manipulated to reduce the one-norm of the Hamiltonian, which +# helps reduce the simulation cost when using block encoding or qubitization techniques. +# Overall, employing CDF-based Hamiltonians for quantum chemistry problems provides a +# relatively promising path to reducing the complexity of fault-tolerant quantum algorithms. +# +# References +# ----------- +# +# .. [#cdf] +# +# Oumarou Oumarou, Maximilian Scheurer, Robert M. Parrish, Edward G. Hohenstein, and Christian Gogolin, +# "Accelerating Quantum Computations of Chemistry Through Regularized Compressed Double Factorization", +# `Quantum 8, 1371 `__, 2024. +# +# .. [#cdf2] +# +# Jeffrey Cohn, Mario Motta, and Robert M. Parrish, +# "Quantum Filter Diagonalization with Compressed Double-Factorized Hamiltonians", +# `PRX Quantum 2, 040352 `__, 2021. +# +# .. [#bliss] +# +# Ignacio Loaiza, and Artur F. Izmaylov, +# "Block-Invariant Symmetry Shift: Preprocessing technique for second-quantized Hamiltonians to improve their decompositions to Linear Combination of Unitaries", +# `arXiv:2304.13772 `__, 2023. +# +# .. [#trotter] +# +# Sergiy Zhuk, Niall Robertson, and Sergey Bravyi, +# "Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation", +# `arXiv:2306.12569 `__, 2023. +# +# About the author +# ---------------- +# diff --git a/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/metadata.json b/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/metadata.json new file mode 100644 index 0000000000..a23bde2ef8 --- /dev/null +++ b/demonstrations_v2/tutorial_how_to_build_compressed_double_factorized_hamiltonians/metadata.json @@ -0,0 +1,80 @@ +{ + "title": "How to build compressed double-factorized Hamiltonians", + "authors": [ + { + "username": "whatsis" + } + ], + "dateOfPublication": "2025-03-05T09:00:00+00:00", + "dateOfLastModification": "2025-04-28T09:00:00+00:00", + "categories": [ + "Quantum Chemistry", + "Algorithms", + "How-to" + ], + "tags": [], + "previewImages": [ + { + "type": "thumbnail", + "uri": "/_static/demo_thumbnails/regular_demo_thumbnails/thumbnail_how_to_build_cdf_hamiltonians.png" + }, + { + "type": "large_thumbnail", + "uri": "/_static/demo_thumbnails/large_demo_thumbnails/thumbnail_large_how_to_build_cdf_hamiltonians.png" + } + ], + "seoDescription": "Learn how to build compressed double-factorized Hamiltonians with PennyLane.", + "doi": "", + "references": [ + { + "id": "bliss", + "type": "article", + "title": "Block-Invariant Symmetry Shift: Preprocessing technique for second-quantized Hamiltonians to improve their decompositions to Linear Combination of Unitaries", + "authors": "Ignacio Loaiza, Artur F. Izmaylov", + "year": "2023", + "publisher": "arXiv", + "url": "https://arxiv.org/abs/2304.13772" + }, + { + "id": "cdf", + "type": "article", + "title": "Accelerating Quantum Computations of Chemistry Through Regularized Compressed Double Factorization", + "authors": "Oumarou Oumarou, Maximilian Scheurer, Robert M. Parrish, Edward G. Hohenstein, Christian Gogolin", + "year": "2024", + "publisher": "Quantum", + "url": "https://doi.org/10.22331/q-2024-06-13-1371" + }, + { + "id": "cdf2", + "type": "article", + "title": "Quantum Filter Diagonalization with Compressed Double-Factorized Hamiltonians", + "authors": "OJeffrey Cohn, Mario Motta, Robert M. Parrish", + "year": "2021", + "publisher": "PRX Quantum", + "url": "https://doi.org/10.1103/PRXQuantum.2.040352" + }, + { + "id": "trotter", + "type": "article", + "title": "Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation", + "authors": "Sergiy Zhuk, Niall Robertson, Sergey Bravyi", + "year": "2023", + "publisher": "arXiv", + "url": "https://arxiv.org/pdf/2306.12569" + } + ], + "basedOnPapers": [], + "referencedByPapers": [], + "relatedContent": [ + { + "type": "demonstration", + "id": "tutorial_quantum_chemistry", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_fermionic_operators", + "weight": 1.0 + } + ] +} diff --git a/demonstrations_v2/tutorial_liealgebra/demo.py b/demonstrations_v2/tutorial_liealgebra/demo.py index dcb419dcff..2a482f1fd9 100644 --- a/demonstrations_v2/tutorial_liealgebra/demo.py +++ b/demonstrations_v2/tutorial_liealgebra/demo.py @@ -343,9 +343,9 @@ def IsingGenerators(n, bc="open"): # Now that we know that the Heisenberg model Hamiltonian commutes with any :math:`S_\text{tot}^{\alpha}` for :math:`\alpha \in \{x, y, z\},` we also know that any observable # composed of the total spin components # -# .. math:: \hat{O} = c_x S^x_\text{tot} + c_x S^y_\text{tot} + c_x S^z_\text{tot} +# .. math:: \hat{O} = c_x S^x_\text{tot} + c_y S^y_\text{tot} + c_z S^z_\text{tot} # -# commutes with the Hamiltonian, +# with arbitrary real coefficients :math:`c_x, c_y, c_z \in \mathbb{R}` commutes with the Hamiltonian, # # .. math:: [\hat{O}, H_\text{Heis}] = 0. # diff --git a/demonstrations_v2/tutorial_liealgebra/metadata.json b/demonstrations_v2/tutorial_liealgebra/metadata.json index cbaff49835..85957c653b 100644 --- a/demonstrations_v2/tutorial_liealgebra/metadata.json +++ b/demonstrations_v2/tutorial_liealgebra/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2024-02-27T00:00:00+00:00", - "dateOfLastModification": "2025-02-03T00:00:00+00:00", + "dateOfLastModification": "2025-02-26T00:00:00+00:00", "categories": [ "Quantum Computing", "Getting Started" diff --git a/demonstrations_v2/tutorial_noisy_circuits/demo.py b/demonstrations_v2/tutorial_noisy_circuits/demo.py index 1491215009..af1a7f4d7e 100644 --- a/demonstrations_v2/tutorial_noisy_circuits/demo.py +++ b/demonstrations_v2/tutorial_noisy_circuits/demo.py @@ -144,12 +144,13 @@ def bitflip_circuit(p): qml.CNOT(wires=[0, 1]) qml.BitFlip(p, wires=0) qml.BitFlip(p, wires=1) - return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1)) + return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1)), qml.state() ps = [0.001, 0.01, 0.1, 0.2] for p in ps: - print(f"QNode output for bit flip probability {p} is {bitflip_circuit(p):.4f}") + result = bitflip_circuit(p) + print(f"QNode output for bit flip probability {p} is {result[0]:.4f}") ###################################################################### @@ -158,7 +159,7 @@ def bitflip_circuit(p): # mitigation and error correction are so important. We can use PennyLane to look under the hood and # see the output state of the circuit for the largest noise parameter -print(f"Output state for bit flip probability {p} is \n{np.real(dev.state)}") +print(f"Output state for bit flip probability {p} is \n{result[1]}") ###################################################################### # Besides the bit flip channel, PennyLane supports several other noisy channels that are commonly diff --git a/demonstrations_v2/tutorial_noisy_circuits/metadata.json b/demonstrations_v2/tutorial_noisy_circuits/metadata.json index 92522adae3..0ba26ef0b7 100644 --- a/demonstrations_v2/tutorial_noisy_circuits/metadata.json +++ b/demonstrations_v2/tutorial_noisy_circuits/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2021-02-22T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2024-12-18T00:00:00+00:00", "categories": [ "Getting Started" ], diff --git a/demonstrations_v2/tutorial_odegen/demo.py b/demonstrations_v2/tutorial_odegen/demo.py index 257bb3884c..8894ff68d8 100644 --- a/demonstrations_v2/tutorial_odegen/demo.py +++ b/demonstrations_v2/tutorial_odegen/demo.py @@ -45,10 +45,10 @@ that can be executed on hardware. -SPS & ODEgen ------------- +SPS and ODEgen gradient rules +----------------------------- -Let us start by deriving both the SPS rule and ODEgen. +Let us start by deriving both the SPS and ODEgen rules. We are interested in cost functions of the form @@ -163,6 +163,7 @@ Let us define it in PennyLane and also import some libraries that we are going to need for this demo. """ + import pennylane as qml import numpy as np import jax.numpy as jnp @@ -181,21 +182,21 @@ ############################################################################## # We are going to consider a system of transmon qubits described by the Hamiltonian -# +# # .. math:: H(\theta, t) = - \sum_i \frac{\omega_i}{2} Z_i + \sum_i \Omega_i(t) \sin(\nu_i t + \phi_i(t)) Y_i + \sum_{q, p \in \mathcal{C}} \frac{g_{qp}}{2} (X_i X_p + Y_i Y_p). -# -# The first term describes the single qubits with frequencies :math:`\omega_i.` -# The second term describes the driving with drive amplitudes :math:`\Omega_i,` drive frequencies :math:`\nu_i` and phases :math:`\phi_i.` -# You can check out our :doc:`recent demo on driving qubits on OQC's Lucy ` if +# +# The first term describes the single qubits with frequencies :math:`\omega_i.` +# The second term describes the driving with drive amplitudes :math:`\Omega_i,` drive frequencies :math:`\nu_i` and phases :math:`\phi_i.` +# You can check out our :doc:`recent demo on driving qubits on OQC's Lucy ` if # you want to learn more about the details of controlling transmon qubits. -# The third term describes the coupling between neighboring qubits. We only have two qubits and a simple topology of +# The third term describes the coupling between neighboring qubits. We only have two qubits and a simple topology of # :math:`\mathcal{C} = \{(0, 1)\},` hence only one coupling constant :math:`g_{01}.` -# The coupling is necessary to generate entanglement, which is achieved with cross-resonant driving in fixed-coupling +# The coupling is necessary to generate entanglement, which is achieved with cross-resonant driving in fixed-coupling # transmon systems, as is the case here. -# -# We will use realistic parameters for the transmons, taken from the `coaxmon design paper `_ [#Patterson]_ +# +# We will use realistic parameters for the transmons, taken from the `coaxmon design paper `_ [#Patterson]_ # (this is the blue-print for the transmon qubits in OQC's Lucy that you can :doc:`access on a pulse level in PennyLane `). -# In order to prepare the singlet ground state, we will perform a cross-resonance pulse, i.e. driving one qubit at its coupled neighbor's +# In order to prepare the singlet ground state, we will perform a cross-resonance pulse, i.e. driving one qubit at its coupled neighbor's # frequency for entanglement generation (see [#Patterson]_ or [#Krantz]_) while simultaneously driving the other qubit on resonance. # We choose a gate time of :math:`100 \text{ ns}.` We will use a piecewise constant function :func:`~pennylane.pulse.pwc` to parametrize both # the amplitude :math:`\Omega_i(t)` and the phase :math:`\phi_i(t)` in time, with ``t_bins = 10`` time bins to allow for enough flexibility in the evolution. @@ -221,7 +222,7 @@ def wrapped(p, t): H_pulse += drive_field(T_CR, qubit_freq[0]) * qml.PauliY(wires[1]) # off-resonance driving of qubit 1 ############################################################################## -# We can now define the cost function that computes the expectation value of +# We can now define the cost function that computes the expectation value of # the Heisenberg Hamiltonian after evolving the state with the parametrized pulse Hamiltonian. # We then define the two separate qnodes with ODEgen and SPS as their differentiation methods, respectively. @@ -235,7 +236,14 @@ def qnode0(params): value_and_grad_jax = jax.jit(jax.value_and_grad(qnode_jax)) num_split_times = 8 -qnode_sps = qml.QNode(qnode0, dev, interface="jax", diff_method=qml.gradients.stoch_pulse_grad, use_broadcasting=True, num_split_times=num_split_times) +gradient_kwargs = {"use_broadcasting": True, "num_split_times": num_split_times} +qnode_sps = qml.QNode( + qnode0, + dev, + interface="jax", + diff_method=qml.gradients.stoch_pulse_grad, + gradient_kwargs=gradient_kwargs, +) value_and_grad_sps = jax.value_and_grad(qnode_sps) qnode_odegen = qml.QNode(qnode0, dev, interface="jax", diff_method=qml.gradients.pulse_odegen) @@ -304,12 +312,12 @@ def partial_step(grad_circuit, opt_state, theta): ############################################################################## -# We see that with analytic gradients (ODEgen), we can reach the ground state energy within 100 epochs, whereas with SPS gradients we cannot find the path +# We see that with analytic gradients (ODEgen), we can reach the ground state energy within 100 epochs, whereas with SPS gradients we cannot find the path # towards the minimum due to the stochasticity of the gradient estimates. Note that both optimizations start from the same (random) initial point. # This picture solidifies when repeating this procedure for multiple runs from different random initializations, as was demonstrated in [#Kottmann]_. # # We also want to make sure that this is a fair comparison in terms of quantum resources. In the case of ODEgen, we maximally have :math:`\mathcal{R}_\text{ODEgen} = 2 (4^n - 1) = 30` expectation values. -# For SPS we have :math:`2 N_g N_s = 32` (due to :math:`N_g = 2` and :math:`N_s=8` time samples per gradient that we chose in ``num_split_times`` above). Thus, overall, we require fewer +# For SPS we have :math:`2 N_g N_s = 32` (due to :math:`N_g = 2` and :math:`N_s=8` time samples per gradient that we chose in ``num_split_times`` above). Thus, overall, we require fewer # quantum resources for ODEgen gradients while achieving better performance. # # Conclusion @@ -323,9 +331,8 @@ def partial_step(grad_circuit, opt_state, theta): # Running VQE using ODEgen on hardware has recently been demonstrated in [#Kottmann]_ and you can directly find `the code here `_. - ############################################################################## -# +# # References # ---------- # diff --git a/demonstrations_v2/tutorial_odegen/metadata.json b/demonstrations_v2/tutorial_odegen/metadata.json index ab2eb2a794..df0ad04903 100644 --- a/demonstrations_v2/tutorial_odegen/metadata.json +++ b/demonstrations_v2/tutorial_odegen/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2023-12-12T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2025-01-28T00:00:00+00:00", "categories": [ "Optimization", "Quantum Computing", diff --git a/demonstrations_v2/tutorial_period_finding/demo.py b/demonstrations_v2/tutorial_period_finding/demo.py new file mode 100644 index 0000000000..7c09addfd4 --- /dev/null +++ b/demonstrations_v2/tutorial_period_finding/demo.py @@ -0,0 +1,387 @@ +r""" +Period finding: A problem at the heart of quantum computing +============================================================= + +You might have heard that Shor's algorithm is an instance of "period finding". You might also have heard that, more generally, this is an example of an *Abelian hidden subgroup problem* solved by a quantum computer's stunning ability to implement the Fourier transform efficiently for an intractable number of function values. Hidden subgroup problems and quantum Fourier transforms were all the rage in the quantum computing literature in the 2000s. + +While trends may have moved on, the idea of extracting group structure from the Fourier spectrum is still at the very core of what quantum computers could be useful for. Scott Aaronson, for example, in his 2022 commentary "How Much Structure Is Needed for Huge Quantum Speedup?" presents the following hierarchy: + +.. figure:: ../_static/demonstration_assets/period_finding/aaronson_fig6.png + :align: center + :width: 75% + :target: javascript:void(0); + + Figure 1. Aaronson's 'hierarchy of structure' for quantum algorithms (Fig. 6 in [#Aaronson22]_) + +However, group theory is a huge hurdle for even some of the more seasoned quantum enthusiasts. This demo wants to give a glimpse of what this "Abelian structure" is all about. Luckily, the fruit-fly example of a hidden subgroup problem is just the task of finding the period of an integer valued function—something one can appreciate without much group jargon. A fantastic resource to read up on the basics is the review of hidden subgroup problems by Childs & van Dam (2010) [#Childs2010]_. +""" + +##################################################################### +# The period finding problem +# ------------------- +# +# Consider a function :math:`f(x)` that maps from integers, say the numbers between 0 and 11, to +# some other domain, such as the numbers between 0 and 3. The function is guaranteed to have a +# periodic structure, which means that it repeats after a certain number of integers. We need a +# further technical requirement, which is that the function does not have the same value within a +# period. Here is an example of such a function with a period of 4: +# +# .. figure:: ../_static/demonstration_assets/period_finding/periodic_function.png +# :align: center +# :width: 90% +# :target: javascript:void(0); +# :figwidth: 80% +# +# Figure 2. Example of a discrete periodic function f(x) over the integers x = 0,...,11. +# The function only takes the same value when moving exactly 4 integers on the x-axis. +# +# Importantly, we assume that we have *black-box access* to that function. For a python coder, this +# abstract technical term is actually quite intuitive: imagine some part of your code returns a python function ``f`` +# that you can call by using integers :math:`x \in \{0,..,11\}` as arguments, like ``f(2)``. However, +# you have no knowledge of the definition ``def f(x)``. Your only way to +# learn about the function is to evaluate it at every point. We're ultimately interested in cases +# where there are too many x-values to evaluate each function value and recover the period in a +# brute force manner. +# +# What does period finding have to do with groups? Well, in the language of group theory, the +# integers from 0 to 11 (together with an appropriate operation, like addition modulo 12) form a +# so-called *cyclic group*, which is an example of *Abelian* groups that Aaronson referred to above. +# +# .. admonition:: Abelian group +# :class: note +# +# A group is a set of elements that has: +# 1. an operation that maps two elements a and b of the set into a third element of the set, for example c = a + b, +# 2. an "identity element" e such that e + a = a for any element a, and +# 3. an inverse -a for every element a, such that a + (-a) = e. +# +# A group is called "Abelian" if a + b = b + a for all a and b, otherwise it is called non-Abelian. +# In general, Abelian groups are simpler to work with than non-Abelian groups. +# +# With this definition we can easily see that the group of integers {0,1,...,11} with addition modulo 12 forms an Abelian group: +# there is an addition operation such that a + b = b + a, the element 0 is the identity element since 0 + a = a for +# all a, and the inverse of element a modulo 12 is 12 - a. +# +# The values {0,4,8} form a *subgroup* that is "generated" by the period 4. Finding the period +# means to find the subgroup. The function is said to "hide" the subgroup, since it has the same +# value on all its elements, effectively labeling them. It is also constant on the *cosets* of +# the subgroups, which are shifted sets of the form {0+s, 4+s, 8+s} (where s is an +# element of the group, or a number between 0 and 11). +# +# .. figure:: ../_static/demonstration_assets/period_finding/periodic_function_groups.png +# :align: center +# :width: 70% +# :target: javascript:void(0); +# +# Figure 3. Periodic function from above, from the perspective of group theory. +# +# This jargon hopefully does not scare you, but should illustrate that +# period finding—and the algorithm we will have a look at here—can be generalised to other groups and +# subgroups. It can even be efficiently run for some instances of non-Abelian groups. +# +# The quantum algorithm to find the period of the function above is really simple. Encode the +# function into a quantum state of the form :math:`\sum_x |x \rangle |f(x) \rangle`, apply the quantum +# Fourier transform on the first register, and measure. We then need to do a bit of post-processing: +# The state we measure, written as an integer, is a multiple of the number of periods +# that "fit" into the x-domain. The number of periods is then the biggest number that divides all the measurement outcomes, +# and the period is the number of integers divided by the number of periods. +# In the example above, the quantum algorithm would return random bitstrings that can be +# interpreted as integers {0, 3, 6, 9}. By sampling only two distinct values, say 6 and +# 9, we can determine the largest common divisor as 3, which is the number of periods fitting +# into 12. From there we can recover the number of periods 12/3 = 4. +# +# We'll now implement this recipe, but use a more convenient cyclic group of 16 elements +# {0,..,15}, which can be exactly represented by 4 qubits. + + +##################################################################### +# Implementation of period finding in PennyLane +# ---------------- + +import pennylane as qml +import numpy as np +import matplotlib.pyplot as plt + +##################################################################### +# First we need to define the periodic function :math:`f(x)`. As mentioned, this function is considered +# unknown in principle, but we can call it on a classical computer, and—theoretically, by turning the +# classical logical circuit into a reversible quantum circuit—we assume we can also call it on a quantum computer. The call +# on a quantum computer can be implemented in parallel for a superposition of inputs, which is +# part of the trick of this algorithm. + + +def f(x): + """ + Function whose period we want to find. + + Args: + x (int): integer in {0,..,15} + + Returns: + integer in {0,..,3} + """ + return x % 8 + + +# let's plot this! +x = range(16) +y = [f(x_) for x_ in x] +fig = plt.figure() +ax = fig.add_axes((0.1, 0.2, 0.8, 0.7)) +ax.scatter(x, y) +ax.set_title("Periodic function") +ax.set_ylabel("f(x)") +ax.set_xlabel("x") +fig.text(0.5, 0.05, "Figure 4. Periodic function f(x)", horizontalalignment='center') +plt.show() + +##################################################################### +# +# We will represent the :math:`x` and :math:`f(x)` values of this function as computational basis states +# :math:`| x \rangle | f(x) \rangle`. The amplitudes belonging to that state can be interpreted as a +# "weight" for a specific function value. +# To give an example, the point :math:`f(2) = 2` can be expressed by preparing a state that has a +# nonzero uniform amplitude at :math:`| 0010 \rangle | 11 \rangle`, but zero amplitudes at the states +# :math:`| 0010 \rangle | 00 \rangle`, :math:`| 0010 \rangle | 01 \rangle`, :math:`| 0010 \rangle | 10 \rangle`. +# +# Since we move between integers and computational basis states, we need two utility conversion functions. + + +def to_int(binary): + return int(binary, 2) + + +def to_binary(integer, n): + return format(integer, "b").zfill(n) + + +##################################################################### +# Now we need an oracle that implements the function ``f`` in "quantum parallel". Applied to a +# superposition of inputs, :math:`\sum_x |x \rangle |0 \rangle`, this unitary prepares the state +# :math:`\sum_x |x \rangle |f(x) \rangle`. There are many such unitaries and here we'll somewhat hack it +# together by defining a matrix that does the job. Of course, in a real application, this +# would be a quantum circuit defined by a sequence of gates. + + +def Oracle(f): + """ + Defines the unitary that implements a function f:{0,..,15} -> {0,..,7}. + + Args: + f (func): function to implement + + Returns: + ndarray representing the unitary + """ + + U = np.zeros((2**7, 2**7)) + + for x in range(2**4): + for f_x in range(2**3): + # we know that the initial state has only support on basis states + # of the form |x>|0>, and therefore write all necessary information + # into the entries of the unitary acting on those states + if f_x == 0: + i = to_int(to_binary(x, 4) + "0" * 3) + j = to_int(to_binary(x, 4) + to_binary(f(x), 3)) + U[i, j] = 1 + U[j, i] = 1 + else: + # in all other cases we use trivial entries + i = to_int(to_binary(x, 4) + to_binary(f_x, 3)) + U[i, i] = 1 + + # check that this is a unitary + assert np.allclose(U @ np.linalg.inv(U), np.eye(2**7)) + + return qml.QubitUnitary(U, wires=range(7)) + + +##################################################################### +# Now we're ready to go. Let's implement the famous period finding algorithm. It consists of a +# quantum routine and a classical post-processing step. As mentioned, the quantum part is simple: prepare the +# desired initial state :math:`\sum_x |x \rangle |f(x) \rangle`, apply the quantum Fourier transform +# onto the :math:`|x\rangle`-register, and measure in the computational basis. Effectively, this +# measures in the "Fourier basis", which is where all the magic happens. +# +# We only have to get two unique samples, from which we can compute the period of the function. For this +# reason we define a device with 2 shots. We also add some snapshots to the circuit that we will look at later. + + +dev = qml.device("default.qubit", wires=7, shots=2) + + +@qml.qnode(dev) +def circuit(): + """Circuit to implement the period finding algorithm.""" + + for i in range(4): + qml.Hadamard(wires=i) + + qml.Snapshot("initial_state") + + Oracle(f) + + qml.Snapshot("loaded_function") + + qml.QFT(wires=range(4)) + + qml.Snapshot("fourier_spectrum") + + return qml.sample(wires=range(4)) + + +# take two samples from the circuit +samples = circuit() + +##################################################################### +# The classical post-processing is relatively simple for period finding (whereas +# for general hidden subgroup problems it requires some less trivial algebra). +# If you are curious about the details, you'll find them in +# Childs & van Dam (2013) [#Childs2010]_, Section IVa. + +from fractions import Fraction +from math import lcm + +# convert the bistrings to integers +sample1_int = to_int("".join(str(s) for s in samples[0])) +sample2_int = to_int("".join(str(s) for s in samples[1])) + +# get the denominator of the fraction representation +denominator1 = Fraction(sample1_int / (2**4)).denominator +denominator2 = Fraction(sample2_int / (2**4)).denominator + +# get the least common multiple +result = lcm(denominator1, denominator2) + +print(f"Hidden period: {result}") + +##################################################################### +# Yay, we've found the hidden period! Now, of course this is only impressive if we increase the +# number of qubits. We can find the hidden period of a function defined on :math:`2^n` values in time +# :math:`O(n \mathrm{log}(n))`. + + +##################################################################### +# A peep at the states involved +# ----------------------------- +# +# Somehow the quantum Fourier transform exposed the information we needed. Let's have a closer +# look at the states that were prepared by making use of the snapshots we recorded during the +# circuit simulation. + +dev = qml.device("default.qubit", wires=7, shots=1) +qnode = qml.QNode(circuit, dev) +intermediate_states = qml.snapshots(circuit)() + +##################################################################### +# We can plot them as discrete functions, where the size of a point indicates the absolute value +# of the amplitude. + +import matplotlib.patches as mpatches + +x = [i for i in range(2**4) for _ in range(2**3)] +y = list(range(2**3)) * (2**4) + +states_to_plot = ["initial_state", "loaded_function", "fourier_spectrum"] + +fig, axes = plt.subplots(3, 1, sharex=True, figsize=(7, 6)) + +for state, ax in zip(states_to_plot, axes): + + real = np.real(intermediate_states[state]) + imag = np.imag(intermediate_states[state]) + + # scale the points by the absolute value of the amplitude + sizes_real = [] + for x_, y_ in zip(x, y): + idx = to_int(to_binary(x_, 4) + to_binary(y_, 3)) + sizes_real.append(200 * np.abs(real[idx])) + + sizes_imag = [] + for x_, y_ in zip(x, y): + idx = to_int(to_binary(x_, 4) + to_binary(y_, 3)) + sizes_imag.append(200 * np.abs(imag[idx])) + + ax.scatter(np.array(x), y, s=sizes_real, c="green", alpha=0.5) + ax.scatter(np.array(x), y, s=sizes_imag, c="pink", alpha=0.5) + ax.set_title(state) + ax.set_ylabel("f(x)") + +plt.xlabel("x") +green = mpatches.Patch(color="green", label="abs(real part)") +pink = mpatches.Patch(color="pink", label="abs(imaginary part)") +plt.legend(handles=[green, pink]) +fig.suptitle("Figure 5. Representing the Fourier spectrum", horizontalalignment='center') +fig.tight_layout() +plt.show() + +##################################################################### +# +# First of all, we see that the oracle really prepared a quantum state that we can interpret as +# our periodic function ``f``. +# +# Furthermore, the state representing the Fourier spectrum looks actually quite interesting. But +# the important feature of this state is also clearly visible: Amplitudes concentrate in :math:`x` +# values that are multiples of 2. And 2 was exactly the amount of periods of 8 that fit into 16. + +##################################################################### +# What is the "Abelian structure" exploited by the quantum Fourier transform? +# ---------------------------------------------------------------------------- +# +# Why does the Fourier transform prepare a state with such a convenient structure? Let us have a +# look how it acts on a superposition of inputs :math:`x` for which :math:`f(x)` has the same value, such as +# :math:`\frac{1}{\sqrt{2}} (|2\rangle + |10\rangle)` (where we translated bitstrings to integers for +# better readability). The quantum Fourier transform prepares a new state whose amplitudes are +# proportional to +# +# .. math:: +# +# \sum_k \left(e^{\frac{2 \pi i 2 k}{16}} + e^{\frac{2 \pi i 10 k}{16}} \right) |k \rangle. +# +# In the exponent you find the values 2 and 10, as well as the size of the group, 16. +# Somewhat magically, for some :math:`|k \rangle`, all exponential functions in the sum evaluate to :math:`1`, +# while for all others, the functions cancel each other out and evaluate exactly to zero. + +for k in range(16): + res = 0.5*(np.exp(2 * np.pi * 1j * 2 * k / 16) + np.exp(2 * np.pi * 1j * 10 * k / 16)) + print(f"k={k} --- {np.round(res, 13)}") + +##################################################################### +# This pattern is the same for whichever set of :math:`x` with the same value :math:`f(x)` +# (in other words, which "coset") we picked. It is also +# true for whichever period we encoded, function we chose, and even which (Abelian) group we +# started with. +# +# The "magic" interference, of course, would not surprise a group theorist; it is inherent to the +# structure of the group. The functions :math:`e^{\frac{2 \pi i x k}{16}}`, which we are so used to +# see in quantum theory, are nothing but the *characters* of a group: functions that take group +# elements :math:`x` to complex numbers, and -- a little like eigenvalues in linear algebra -- capture +# the essentials of the group. The destructive interference we see in the quantum Fourier +# transform is nothing other than the orthogonality theorem for these characters. +# +# It is strangely beautiful how quantum interference leverages the structure of Abelian groups in +# the quantum Fourier transform. So far this has only found application for rather abstract +# problems. Even though cryptography, impacted by Shor's algorithm, is what one might consider a +# "real-world" application, it is an outlier, since problems in cryptography are artificially +# *constructed* from abstract mathematical theory. A fascinating question is whether "real" +# real-world applications, like material science or data analysis, could benefit from this +# remarkable confluence of group structure and quantum theory. +# +# +# References +# ------------ +# +# .. [#Aaronson22] +# +# Scott Aaronson, "How Much Structure Is Needed for Huge Quantum Speedups?", +# `arXiv:2209.06930 `__, 2022 +# +# .. [#Childs2010] +# +# Andrew Childs, Vim van Dam, "Quantum algorithms for algebraic problems", +# `Reviews of Modern Physics 82.1: 1-52. `__, 2010 +# +# About the author +# ---------------- +# diff --git a/demonstrations_v2/tutorial_period_finding/metadata.json b/demonstrations_v2/tutorial_period_finding/metadata.json new file mode 100644 index 0000000000..0c216aaa6d --- /dev/null +++ b/demonstrations_v2/tutorial_period_finding/metadata.json @@ -0,0 +1,61 @@ +{ + "title": "Period finding: A problem at the heart of quantum computing", + "authors": [ + { + "username": "mariaschuld" + } + ], + "dateOfPublication": "2025-04-16T10:00:00+00:00", + "dateOfLastModification": "2025-04-29T10:00:00+00:00", + "categories": [ + "Algorithms" + ], + "tags": [], + "previewImages": [ + { + "type": "thumbnail", + "uri": "/_static/demo_thumbnails/regular_demo_thumbnails/thumbnail_period_finding.png" + }, + { + "type": "large_thumbnail", + "uri": "/_static/demo_thumbnails/large_demo_thumbnails/thumbnail_large_period_finding.png" + } + ], + "seoDescription": "Find out how quantum computers exploit group structure and the Fourier transform to uncover hidden patterns in functions.", + "doi": "", + "references": [ + { + "id": "Aaronson2022", + "type": "preprint", + "title": "How Much Structure Is Needed for Huge Quantum Speedups?", + "authors": "Scott Aaronson", + "year": "2022", + "journal": "arXiv", + "url": "https://arxiv.org/pdf/2209.06930" + }, + { + "id": "Childs2010", + "type": "article", + "title": "Quantum algorithms for algebraic problems", + "authors": "Andrew M. Childs, Wim van Dam", + "year": "1997", + "journal": "Reviews of Modern Physics 82, 1", + "doi": "10.1103/RevModPhys.82.1", + "url": "https://doi.org/10.1103/RevModPhys.82.1" + } + ], + "basedOnPapers": [], + "referencedByPapers": [], + "relatedContent": [ + { + "type": "demonstration", + "id": "tutorial_qft", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_qft_arithmetics", + "weight": 1.0 + } + ] +} diff --git a/demonstrations_v2/tutorial_pulse_programming101/demo.py b/demonstrations_v2/tutorial_pulse_programming101/demo.py index b1b4ce33ea..aa050030a1 100644 --- a/demonstrations_v2/tutorial_pulse_programming101/demo.py +++ b/demonstrations_v2/tutorial_pulse_programming101/demo.py @@ -322,15 +322,21 @@ def wrapped(p, t): ############################################################################## # Now we define the ``qnode`` that computes the expectation value of the molecular Hamiltonian. +# We need to wrap the ``qnode`` in a function so that we can convert the expectation value to a real number. +# This will enable use to make use of gradient descent methods that require real-valued loss functions. dev = qml.device("default.qubit", wires=range(n_wires)) -@qml.qnode(dev, interface="jax") def qnode(theta, t=duration): - qml.BasisState(jnp.array(data.tapered_hf_state), wires=H_obj.wires) - qml.evolve(H_pulse)(params=(*theta, *theta), t=t) - return qml.expval(H_obj) + @qml.qnode(dev) + def _qnode_inner(theta, t=duration): + qml.BasisState(jnp.array(data.tapered_hf_state), wires=H_obj.wires) + qml.evolve(H_pulse)(params=(*theta, *theta), t=t) + return qml.expval(H_obj) + + expectation_value = _qnode_inner(theta, t) # Execute the qnode + return jnp.real(expectation_value) # Typecast to real number value_and_grad = jax.jit(jax.value_and_grad(qnode)) diff --git a/demonstrations_v2/tutorial_pulse_programming101/metadata.json b/demonstrations_v2/tutorial_pulse_programming101/metadata.json index 7801031d1f..9c5a8f71c3 100644 --- a/demonstrations_v2/tutorial_pulse_programming101/metadata.json +++ b/demonstrations_v2/tutorial_pulse_programming101/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2023-03-08T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2025-02-14T00:00:00+00:00", "categories": [ "Quantum Hardware", "Quantum Computing" diff --git a/demonstrations_v2/tutorial_qaoa_intro/demo.py b/demonstrations_v2/tutorial_qaoa_intro/demo.py index 789d724ae4..8e73c0d08d 100644 --- a/demonstrations_v2/tutorial_qaoa_intro/demo.py +++ b/demonstrations_v2/tutorial_qaoa_intro/demo.py @@ -38,7 +38,7 @@ When considering quantum circuits, it is often convenient to define them by a series of quantum gates. But there are many instances where it is useful to think of a quantum circuit in terms of a -`Hamiltonian `__. +`Hamiltonian `__. Indeed, gates are physically implemented by performing time evolution under a carefully engineered Hamiltonian. These transformations are described by the time evolution operator, which is a unitary defined as:""" diff --git a/demonstrations_v2/tutorial_qaoa_intro/metadata.json b/demonstrations_v2/tutorial_qaoa_intro/metadata.json index 2fc2f5eb95..a383349c67 100644 --- a/demonstrations_v2/tutorial_qaoa_intro/metadata.json +++ b/demonstrations_v2/tutorial_qaoa_intro/metadata.json @@ -29,5 +29,12 @@ "weight": 1.0 } ], + "hardware": [ + { + "id": "qbraid", + "link": "https://account.qbraid.com?gitHubUrl=https://github.com/PennyLaneAI/pennylane-demo-notebooks.git", + "logo": "/_static/hardware_logos/qbraid.png" + } + ], "discussionForumUrl": "https://discuss.pennylane.ai/t/qaoa-and-optimization/5188" } diff --git a/demonstrations_v2/tutorial_qft/demo.py b/demonstrations_v2/tutorial_qft/demo.py index 487e1a68c6..a85e009dbb 100644 --- a/demonstrations_v2/tutorial_qft/demo.py +++ b/demonstrations_v2/tutorial_qft/demo.py @@ -168,7 +168,7 @@ def circuit(): # This value corresponds to an approximation of :math:`2^nf` where :math:`f` is the frequency and :math:`n` is the # number of qubits. # Once we know the frequency, we invert it to obtain the period :math:`T` of our state. -# In this case, the period is :math:`T = 2^5 / 3 \sim 10.33,` close to the real value of :math:`10.` +# In this case, the period is :math:`T = 2^5 / 3 \sim 10.66,` close to the real value of :math:`10.` # # Conclusion # ---------- diff --git a/demonstrations_v2/tutorial_qft/metadata.json b/demonstrations_v2/tutorial_qft/metadata.json index 99d3b07ae4..99952646e4 100644 --- a/demonstrations_v2/tutorial_qft/metadata.json +++ b/demonstrations_v2/tutorial_qft/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2024-04-16T00:00:00+00:00", - "dateOfLastModification": "2024-11-18T00:00:00+00:00", + "dateOfLastModification": "2025-04-10T00:00:00+00:00", "categories": [ "Algorithms", "Quantum Computing" diff --git a/demonstrations_v2/tutorial_shors_algorithm_catalyst/demo.py b/demonstrations_v2/tutorial_shors_algorithm_catalyst/demo.py new file mode 100644 index 0000000000..5b71df3923 --- /dev/null +++ b/demonstrations_v2/tutorial_shors_algorithm_catalyst/demo.py @@ -0,0 +1,1436 @@ +r""".. role:: html(raw) + :format: html + +Quantum just-in-time compiling Shor's algorithm with Catalyst +============================================================= + +The past few years stimulated a lot of discussion about *hybrid +quantum-classical algorithms*. For a time, this terminology was synonymous +with *variational algorithms*. However, integration with classical +coprocessors is necessary for every quantum algorithm, even ones considered +quintessentially quantum. + +.. figure:: ../_static/demo_thumbnails/opengraph_demo_thumbnails/OGthumbnail_shor_algorithm_catalyst_pennylane.png + :align: center + :width: 70% + :target: javascript:void(0) + +Shor's famous factoring algorithm [#Shor1997]_ is one such example. Consider +an integer :math:`N`, promised to be the product of two primes, :math:`p` and +:math:`q`. Shor's algorithm uses a quantum computer to solve this problem with +exponentially-better scaling than the best-known classical algorithm. But, +have a look at the code below: +""" + +import jax.numpy as jnp + + +def shors_algorithm(N): + p, q = 0, 0 + + while p * q != N: + a = jnp.random.choice(jnp.arange(2, N - 1)) + + if jnp.gcd(N, a) != 1: + p = jnp.gcd(N, a) + return p, N // p + + guess_r = guess_order(N, a) + + if guess_r % 2 == 0: + guess_square_root = (a ** (guess_r // 2)) % N + + if guess_square_root not in [1, N - 1]: + p = jnp.gcd(N, guess_square_root - 1) + q = jnp.gcd(N, guess_square_root + 1) + + return p, q + + +###################################################################### +# If you saw this code out of context, would it even occur to you that it's for +# a quantum algorithm? There are no quantum circuits in sight, and the only "q" +# is a variable name! +# +# As quantum hardware continues to scale up, the way we reason about quantum +# programming is evolving in tandem. Writing circuits gate-by-gate for +# algorithms with hundreds or thousands of qubits is unsustainable. More +# fundamentally, we should consider whether programmers really need deep +# knowledge of quantum computing at all, if the software library can generate +# and compile appropriate quantum code (though, they should probably have at +# least some awareness, since the output of ``guess_order`` is probabilistic!). +# This raises the important questions of what gets compiled, and where, when, +# and how compilation happens. +# +# In PennyLane, classical and quantum code can be compiled *together*, as a +# single workflow, using the `Catalyst +# `_ +# library. This demo leverages their integration to implement Shor's factoring +# algorithm using just-in-time compilation from beginning to end, i.e., classical +# control structure and all. Even better, compilation happens only once per +# distinct *bit-width* of the factored integers, which can lead to huge savings +# in compilation time for realistic problem sizes. +# +# Compilation +# ----------- +# +# Classical compilation +# ^^^^^^^^^^^^^^^^^^^^^ +# +# Compilation is the process of translating operations expressed in a high-level +# language to a low-level language. In compiled languages like C and C++, compilation +# happens offline prior to code execution. A compiler takes a program as input +# and sends it through a sequence of *passes* that perform tasks such as syntax +# analysis, code generation, and optimization. The compiler outputs a new +# program in assembly code, which gets passed to an assembler. The assembler +# translates this code into a machine-executable program that we can feed inputs +# to, then run [#PurpleDragonBook]_. +# +# Compilation is not the only way to execute a program. Python, for example, is +# an *interpreted* language. Both a source program and inputs are fed to the +# interpreter, which processes them line by line and directly yields the program +# output [#PurpleDragonBook]_. +# +# Compilation and interpretation each have strengths and weaknesses. Compilation +# generally leads to faster execution, because optimizations can consider +# the overall structure of a program. However, the executable code is not +# human-readable and thus harder to debug. Interpretation is slower, but +# debugging is often easier because execution can be paused to inspect +# the program state or view diagnostic information [#PurpleDragonBook]_. +# +# In between these extremes lies *just-in-time compilation*. Just-in-time (JIT) +# compilation happens *during* execution. If a programmer specifies a function +# should be JIT compiled, the first time the interpreter sees it, it will spend +# a little more time to construct and cache a compiled version of that +# function. The next time that function is executed, the compiled version can be +# reused, provided the structure of the inputs hasn't changed [#JAXJIT]_. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/compilation_comparison.png +# :scale: 75% +# :align: center +# :alt: The spectrum of compilation techniques: ahead-of-time, just-in-time, and interpretation. +# +# The spectrum of compilation techniques, from ahead-of-time to interpretation. +# +# Quantum compilation +# ^^^^^^^^^^^^^^^^^^^ +# Quantum compilation, like its classical counterpart, lowers an algorithm from +# high-level instructions to low-level instructions. The bulk of +# this process involves converting a circuit expressed as generic, abstract gates +# to a sequence of gates that satisfy the constraints of a particular +# hardware device. Quantum compilation also involves multiple passes through +# one or more intermediate representations, and both +# machine-independent and dependent optimizations, as depicted below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/compilation-stack.svg +# :scale: 75% +# :align: center +# :alt: The quantum compilation stack. +# +# High-level overview of the quantum compilation stack. Each step contains +# numerous subroutines and passes of its own, and many require solving +# computationally hard problems (or very good heuristic techniques). +# +# Developing automated compilation tools is an active area of quantum software +# research. Even if a library contains built-in implementations of quantum +# circuits (e.g., the Fourier transform, or a multi-controlled operation), +# without a proper compiler a user would be left to optimize and map them to +# hardware by hand. This is a laborious (and error-prone!) process, and +# furthermore, is unlikely to be optimal. +# +# Suppose we want to compile and optimize quantum circuits for Shor's algorithm +# to factor an integer :math:`N`. Recalling the pseudocode above, let's break +# the algorithm down into a few distinct steps to identify where quantum +# compilation happens (for a full description of Shor's algorithm, the +# interested reader is referred to the `PennyLane Codebook +# `_). +# +# - Randomly select an integer, :math:`a`, between 2 and +# :math:`N-1` (double check we didn't get lucky and :math:`a` has a common factor with :math:`N`) +# - Using :math:`N` and :math:`a`, generate a circuit for *order finding* on a quantum computer. Execute it, and use the measurement results to obtain a candidate non-trivial square root +# - If the square root is non-trivial, test for valid factors. Otherwise, take more measurement shots, or try a different :math:`a`. +# +# The key thing to note is that for every :math:`N` and :math:`a`, a different +# quantum circuit must be generated, compiled and optimized. Even with a good +# compiler, this will lead to a huge computational overhead! Recall also that in +# a cryptographic context, :math:`N` relates to a public key that is unique for +# every entity. Moreover, for sizes of cryptographic relevance, :math:`N` will +# be a 2048-bit integer (or larger)! It would be ideal if we could reuse and +# share some of this work across different :math:`a` and :math:`N`. To that end, +# JIT compilation is a worthwhile option to explore. +# +# Quantum just-in-time compilation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# In standard PennyLane, quantum circuit execution can be JIT compiled with +# JAX. To learn more, check out the JAX documentation [#JAXJIT]_ and the +# :doc:`PennyLane demo ` on the +# subject. But this only compiles a single circuit. What about all the other +# code around it? What if you also wanted to optimize that quantum circuit, +# based on contextual information? +# +# This is where Catalyst comes in. Catalyst enables quantum JIT (QJIT) +# compilation of the *entire* algorithm from beginning to end. On the surface, it +# looks to be as simple as the following: + +import pennylane as qml + + +@qml.qjit +def shors_algorithm(N): + # Implementation goes here + return p, q + + +###################################################################### +# In practice, it is not so simple, and requires some special-purpose functions +# and data manipulation. But ultimately, we will show how to QJIT the most +# important parts such that the signature can be as minimal as this: + + +@qml.qjit(autograph=True, static_argnums=(1)) +def shors_algorithm(N, n_bits): + # Implementation goes here + return p, q + + +###################################################################### +# Furthermore, along the way we'll leverage the structure of :math:`a` to +# construct more optimal quantum circuits in the QJITted function. +# +# QJIT compiling Shor's algorithm +# ------------------------------- +# +# Quantum subroutines +# ^^^^^^^^^^^^^^^^^^^ +# +# The implementation of the classical portion of Shor's algorithm is +# near-identical to the pseudocode at the beginning, and can be JIT compiled +# essentially as-is. The quantum aspects are where challenges arise. +# +# This section outlines the quantum circuits used in the order-finding +# subroutine. The presented implementation is based on that of Beauregard +# [#Beauregard2003]_. For an integer :math:`N` with an :math:`n = \lceil \log_2 +# N \rceil`-bit representation, we require :math:`2n + 3` qubits, where :math:`n +# + 1` are for computation and :math:`n + 2` are auxiliary. +# +# Order finding is an application of *quantum phase estimation* +# (:doc:`QPE `) for the operator +# +# .. math:: +# +# U_a \vert x \rangle = \vert ax \pmod N \rangle, +# +# where :math:`\vert x \rangle` is the binary representation of integer +# :math:`x`, and :math:`a` is the randomly-generated integer discussed +# above. The full QPE circuit for producing a :math:`t`-bit estimate of a +# phase, :math:`\theta`, is presented below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full.svg +# :width: 600 +# :align: center +# :alt: Quantum phase estimation circuit for order finding. +# +# Initial circuit for order finding with quantum phase estimation (QPE). +# +# This high-level view hides the circuit's complexity; the implementation +# details of :math:`U_a` are not shown, and auxiliary qubits are omitted. In +# what follows, we'll leverage shortcuts afforded by Catalyst and the hybrid +# nature of computation. Specifically, with mid-circuit measurement and reset, we +# can reduce the number of estimation wires to :math:`t=1`. Most arithmetic will +# be performed in the Fourier basis. Finally, with Catalyst, we can vary circuit +# structure based on :math:`a` to save resources, *even though its value isn't +# known until runtime*. +# +# First, we'll use our classical knowledge of :math:`a` to simplify the +# implementation of the controlled :math:`U_a^{2^k}`. Naively, it looks like we +# must apply the controlled :math:`U_a` operations :math:`2^k` times. However, +# note +# +# .. math:: +# +# U_a^{2^k}\vert x \rangle = \vert (a \cdot a \cdots a) x \pmod N \rangle = \vert a^{2^k}x \pmod N \rangle = U_{a^{2^k}} \vert x \rangle. +# +# Since :math:`a` is known, we can classically evaluate :math:`a^{2^k} \pmod +# N` and implement controlled-:math:`U_{a^{2^k}}` instead. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power.svg +# :width: 600 +# :align: center +# :alt: Order finding with controlled operations that take advantage of classical precomputation. +# +# Leveraging knowledge of :math:`a` allows us to precompute powers for the +# controlled :math:`U_a`. +# +# However, there is a tradeoff. The implementation of :math:`U_{a^{2^k}}` will +# be different for each power of :math:`a`, so we +# must compile and optimize more than one circuit. On the other hand, we now run +# only :math:`t` controlled operations instead of :math:`1 + 2 + 4 + \cdots + +# 2^{t-1} = 2^t - 1`. The additional compilation time is likely to be outweighed +# by the reduced number of function calls, especially if we can JIT compile the +# circuit construction. +# +# Next, let's zoom in on an arbitrary controlled :math:`U_a`. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/c-ua.svg +# :width: 700 +# :align: center +# :alt: Quantum phase estimation circuit for order finding. +# +# Implementing a controlled :math:`U_a` with modular arithmetic +# [#Beauregard2003]_. +# +# The control qubit, :math:`\vert c\rangle`, is an estimation qubit. The +# register :math:`\vert x \rangle` and the auxiliary register contain :math:`n` +# and :math:`n + 1` qubits respectively, for reasons described below. +# +# :math:`M_a` multiplies the contents of one register by :math:`a` and adds it to +# another register, in place and modulo :math:`N`, +# +# .. math:: +# +# M_a \vert x \rangle \vert b \rangle \vert 0 \rangle = \vert x \rangle \vert (b + ax) \pmod N \rangle \vert 0 \rangle. +# +# Ignoring the control qubit, we can validate that this circuit implements +# :math:`U_a`: +# +# .. math:: +# +# \begin{eqnarray} +# M_a \vert x \rangle \vert 0 \rangle^{\otimes n + 1} \vert 0 \rangle &=& \vert x \rangle \vert ax \rangle \vert 0 \rangle \\ +# SWAP (\vert x \rangle \vert ax \rangle ) \vert 0 \rangle &=& \vert ax \rangle \vert x \rangle \vert 0 \rangle \\ +# M_{a^{-1}}^\dagger \vert ax \rangle \vert x \rangle \vert 0 \rangle &=& \vert ax\rangle \vert x - a^{-1}(ax) \rangle \vert 0 \rangle \\ +# &=& \vert ax \rangle \vert 0 \rangle^{\otimes n + 1} \vert 0 \rangle, +# \end{eqnarray} +# +# where we've omitted the "mod :math:`N`" for readability, and used the fact +# that the adjoint of addition is subtraction. +# +# A high-level implementation of a controlled :math:`M_a` is shown below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/doubly-controlled-adder-with-control.svg +# :width: 700 +# :align: center +# :alt: Controlled addition of :math:`ax` using a series of double-controlled Fourier adders. +# +# Circuit for controlled multiplication of :math:`ax`. The circuit element labeled :math:`\Phi_+` +# performs addition modulo :math:`N` in the Fourier basis (see main text for a full description) +# [#Beauregard2003]_. +# +# First, note that the controls on the quantum Fourier transforms (QFTs) are not +# needed. If we remove them and :math:`\vert c \rangle = \vert 1 \rangle`, the +# circuit works as expected. If :math:`\vert c \rangle = \vert 0 \rangle`, they +# run, but cancel each other out since none of the operations in between will +# execute (this optimization is broadly applicable to controlled operations, and +# quite useful!). +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/doubly-controlled-adder-with-control-not-on-qft.svg +# :width: 700 +# :align: center +# :alt: Controlled addition of :math:`ax` using a series of double-controlled Fourier adders, controls on QFT removed. +# +# The controls on the QFT can be removed without altering the effect of the circuit +# [#Beauregard2003]_. +# +# At first glance, it's not clear how this produces :math:`a x`. The qubits in +# register :math:`\vert x \rangle` control operations that depend on :math:`a` +# multiplied by various powers of 2. There is a QFT before and after, whose +# purpose is unclear, and we have yet to define the gate labelled :math:`\Phi_+`. +# +# The operation, :math:`\Phi_+`, performs addition modulo :math:`N` in the +# *Fourier basis* [#Draper2000]_. This is another trick that leverages +# prior knowledge of :math:`a`. Rather than performing addition on bits in +# computational basis states, we apply a QFT, adjust the phases based on the +# bits of the number being added, and then an inverse QFT to obtain the result +# in the computational basis. +# +# To understand how Fourier addition works, let's begin with the simpler case of +# non-modular addition. The regular addition circuit, denoted by :math:`\Phi`, +# is shown below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder.svg +# :width: 750 +# :align: center +# :alt: Addition in the Fourier basis. +# +# Circuit for addition in the Fourier basis [#Draper2000]_. Calligraphic +# letters indicate basis states converted to the Fourier basis. +# +# Fourier addition of two :math:`n`-bit numbers uses :math:`n+1` qubits, as it +# accounts for the possibility of *overflow* during addition (this constitutes one +# of our auxiliary qubits). The :math:`\mathbf{R}_k` are rotations that depend +# on the binary representation of :math:`a`, +# +# .. math:: +# +# \mathbf{R}_k = \begin{pmatrix} 1 & 0 \\ 0 & e^{2\pi i\sum_{\ell=0}^{k} \frac{a_\ell}{2^{\ell+1}}} \end{pmatrix}. +# +# A detailed derivation of this circuit is included in the +# :ref:`Appendix `. +# +# Next, we must augment Fourier addition to work modulo :math:`N` (i.e., +# :math:`\Phi_+`). This can be done by adding an auxiliary qubit and a sequence +# of operations to compensate for any overflow that occurs during addition. A +# full explanation of Fourier addition modulo :math:`N` is also provided in the +# :ref:`Appendix `. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_modulo_n.svg +# :width: 800 +# :align: center +# :alt: Addition in the Fourier basis modulo N. +# +# Circuit for doubly-controlled Fourier addition modulo :math:`N` +# [#Beauregard2003]_. Calligraphic letters indicate basis states converted to +# the Fourier basis. +# +# This completes our implementation of the controlled :math:`U_{a^{2^k}}`. The +# full circuit, shown below, uses :math:`t + 2n + 2` qubits. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_combined_multi_est_wires.svg +# :width: 500 +# :align: center +# :alt: Full QPE circuit, all t estimation wires, and decompositions. +# +# Initial circuit for order finding with QPE, and decompositions of its +# constituent subroutines. +# +# +# Taking advantage of classical information +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Above, we incorporated information about :math:`a` in the controlled +# :math:`U_a` and the Fourier adder. In what follows, we identify a few other +# places where classical information can be leveraged. +# +# First, consider the controlled :math:`U_{a^{2^0}} = U_a` at the beginning of +# the algorithm. The only basis state this operation applies to is :math:`\vert +# 1 \rangle`, which gets mapped to :math:`\vert a \rangle`. This is effectively +# just controlled addition of :math:`a - 1` to :math:`1`. Since :math:`a` +# is selected from between :math:`2` and :math:`N-2` inclusive, the addition is +# guaranteed to never overflow. This means we can simply do a controlled Fourier +# addition, and save a significant number of resources! +# +# We can also make some optimizations to the end of the algorithm by keeping +# track of the powers of :math:`a`. If, at iteration :math:`k,` we have +# :math:`a^{2^k} = 1`, no further multiplication is necessary because we would +# be multiplying by :math:`1.` In fact, we can terminate the algorithm early because +# we've found the order of :math:`a` is simply :math:`2^k`. +# +# There are also less-trivial optimizations we can make. Consider the +# sequence of doubly-controlled adders in the controlled :math:`M_a`. Below, we +# show the initial instance where the auxiliary register is in state +# :math:`\vert 0 \rangle`. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/doubly-controlled-adder-simplified.svg +# :width: 600 +# :align: center +# :alt: The controlled multiplier circuit in the context of Shor's algorithm. +# +# The doubly-controlled adders are each controlled on an estimation qubit and a +# qubit in the target register. Consider the state of the system after the +# initial controlled-:math:`U_a` (or, rather, controlled addition of +# :math:`a-1`), +# +# .. math:: +# +# \begin{equation*} +# \vert + \rangle ^{\otimes (t - 1)} \frac{1}{\sqrt{2}} \left( \vert 0 \rangle \vert 1 \rangle + \vert 1 \rangle \vert a \rangle \right). +# \end{equation*} +# +# The next controlled operation is controlled :math:`U_{a^2}`. Since the only +# two basis states present are :math:`\vert 1 \rangle` and :math:`\vert a +# \rangle`, the only doubly-controlled :math:`\Phi_+` that trigger are the first +# one (with the second control on the bottom-most qubit) and those controlled +# on qubits that are :math:`1` in the binary representation of :math:`a`. Thus, +# we only need doubly-controlled operations on qubits where the logical OR of +# the bit representations of :math:`1` and :math:`a` are :math:`1!` We present here an example +# for :math:`a = 5`. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/doubly-controlled-adder-simplified-states.svg +# :width: 800 +# :align: center +# :alt: Leveraging knowledge of :math:`a` to eliminate unnecessary doubly-controlled additions. +# +# Leveraging knowledge of :math:`a` to eliminate unnecessary doubly-controlled additions. +# +# Depending on :math:`a`, this could be major savings, especially at the +# beginning of the algorithm where very few basis states are involved. The same +# trick can be used after the controlled SWAPs, as demonstrated below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/c-ua-basis-states-optimization.svg +# :width: 600 +# :align: center +# :alt: Removing doubly-controlled additions based on expected terms in the register superposition. +# +# Doubly-controlled operations can be removed from controlled multiplications +# by keeping track of basis states involved. +# +# Eventually we expect diminishing returns because each controlled +# :math:`U_{a^{2^k}}` contributes more terms to the superposition. Before the +# :math:`k-` th iteration, the control register contains a superposition of +# :math:`\{ \vert a^j \rangle \}, j = 0, \ldots, 2^{k - 1}` (inclusive), and after the +# controlled SWAPs, the relevant superposition is :math:`\{ \vert a^j \rangle \}, j = +# 2^{k-1}+1, \ldots, 2^{k} - 1`. +# +# The "single-qubit" QPE +# ~~~~~~~~~~~~~~~~~~~~~~ +# +# The optimizations in the previous section focus on reducing the number of +# operations (circuit depth). The number of qubits (circuit width) can be reduced +# using a well-known trick for the QFT. Let's return to the QPE circuit and +# expand the final inverse QFT. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded. +# +# Consider the bottom estimation wire. After the final Hadamard, this qubit only +# applies controlled operations. Rather than preserving its state, we can make a +# measurement and apply rotations classically controlled on the outcome +# (essentially, the reverse of the deferred measurement process). The same can +# be done for the next estimation wire; rotations applied to the remaining +# estimation wires then depend on the previous two measurement outcomes. We can +# repeat this process for all remaining estimation wires. Moreover, we can +# simply reuse the *same* qubit for *all* estimation wires, provided we keep +# track of the measurement outcomes classically, and apply an appropriate +# :math:`RZ` rotation, +# +# .. math:: +# +# \mathbf{M}_{k} = \begin{pmatrix} 1 & 0 \\ 0 & e^{-2\pi i\sum_{\ell=0}^{k} \frac{\theta_{\ell}}{2^{k + 2 - \ell}}} \end{pmatrix}. +# +# This allows us to reformulate the QPE algorithm as: +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-7.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with one estimation qubit. +# +# A step-by-step example is included in the :ref:`Appendix `. +# +# The final Shor circuit +# ~~~~~~~~~~~~~~~~~~~~~~ +# Replacing the controlled-:math:`U_a` with the subroutines derived above, we +# see Shor's algorithm requires :math:`2n + 3` qubits in total, as summarized in +# the graphic below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_combined.svg +# :width: 800 +# :align: center +# :alt: Full implementation of QPE circuit. +# +# From this, we see some additional optimizations are possible. In particular, +# the QFT and inverse QFT are applied before and after each :math:`M_a` block, +# so we can remove the ones that occur between the different controlled +# :math:`U_a` (the only ones remaining are the very first and last, and +# those before and after the controlled SWAPs). + +###################################################################### +# Catalyst implementation +# ^^^^^^^^^^^^^^^^^^^^^^^ +# +# Circuits in hand, we now present the full implementation of Shor's algorithm. +# +# First, we require some utility functions for modular arithmetic: +# exponentiation by repeated squaring (to avoid overflow when exponentiating +# large integers), and computation of inverses modulo :math:`N`. Note that the +# first can be done in regular Python with the built-in ``pow`` method. However, +# it is not JIT-compatible and there is no equivalent in JAX NumPy. + +from jax import numpy as jnp + + +def repeated_squaring(a, exp, N): + """QJIT-compatible function to determine (a ** exp) % N. + + Source: https://en.wikipedia.org/wiki/Modular_exponentiation#Left-to-right_binary_method + """ + bits = jnp.array(jnp.unpackbits(jnp.array([exp]).view("uint8"), bitorder="little")) + total_bits_one = jnp.sum(bits) + + result = jnp.array(1, dtype=jnp.int32) + x = jnp.array(a, dtype=jnp.int32) + + idx, num_bits_added = 0, 0 + + while num_bits_added < total_bits_one: + if bits[idx] == 1: + result = (result * x) % N + num_bits_added += 1 + x = (x**2) % N + idx += 1 + + return result + + +def modular_inverse(a, N): + """QJIT-compatible modular multiplicative inverse routine. + + Source: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Modular_integers + """ + t = jnp.array(0, dtype=jnp.int32) + newt = jnp.array(1, dtype=jnp.int32) + r = jnp.array(N, dtype=jnp.int32) + newr = jnp.array(a, dtype=jnp.int32) + + while newr != 0: + quotient = r // newr + t, newt = newt, t - quotient * newt + r, newr = newr, r - quotient * newr + + if t < 0: + t = t + N + + return t + + +###################################################################### +# We also require a few helper functions for phase estimation. + + +def fractional_binary_to_float(sample): + """Convert an n-bit sample [k1, k2, ..., kn] to a floating point + value using fractional binary representation, + + k = (k1 / 2) + (k2 / 2 ** 2) + ... + (kn / 2 ** n) + """ + powers_of_two = 2 ** (jnp.arange(len(sample)) + 1) + return jnp.sum(sample / powers_of_two) + + +def as_integer_ratio(f): + """QJIT compatible conversion of a floating point number to two 64-bit + integers such that their quotient equals the input to available precision. + """ + mantissa, exponent = jnp.frexp(f) + + i = 0 + while jnp.logical_and(i < 300, mantissa != jnp.floor(mantissa)): + mantissa = mantissa * 2.0 + exponent = exponent - 1 + i += 1 + + numerator = jnp.asarray(mantissa, dtype=jnp.int64) + denominator = jnp.asarray(1, dtype=jnp.int64) + abs_exponent = jnp.abs(exponent) + + if exponent > 0: + num_to_return, denom_to_return = numerator << abs_exponent, denominator + else: + num_to_return, denom_to_return = numerator, denominator << abs_exponent + + return num_to_return, denom_to_return + + +def phase_to_order(phase, max_denominator): + """Given some floating-point phase, estimate integers s, r such that s / r = + phase. Uses a JIT-compatible re-implementation of Fraction.limit_denominator. + """ + numerator, denominator = as_integer_ratio(phase) + + order = 0 + + if denominator <= max_denominator: + order = denominator + + else: + p0, q0, p1, q1 = 0, 1, 1, 0 + + a = numerator // denominator + q2 = q0 + a * q1 + + while q2 < max_denominator: + p0, q0, p1, q1 = p1, q1, p0 + a * p1, q2 + numerator, denominator = denominator, numerator - a * denominator + + a = numerator // denominator + q2 = q0 + a * q1 + + k = (max_denominator - q0) // q1 + bound1 = p0 + k * p1 / q0 + k * q1 + bound2 = p1 / q1 + + loop_res = 0 + + if jnp.abs(bound2 - phase) <= jnp.abs(bound1 - phase): + loop_res = q1 + else: + loop_res = q0 + k * q1 + + order = loop_res + + return order + + +###################################################################### +# Below, we have the implementations of the arithmetic circuits derived in the +# previous section. + +import pennylane as qml +import catalyst +from catalyst import measure + +catalyst.autograph_strict_conversion = True + + +def QFT(wires): + """The standard QFT, redefined because the PennyLane one uses terminal SWAPs.""" + shifts = jnp.array([2 * jnp.pi * 2**-i for i in range(2, len(wires) + 1)]) + + for i in range(len(wires)): + qml.Hadamard(wires[i]) + + for j in range(len(shifts) - i): + qml.ControlledPhaseShift(shifts[j], wires=[wires[(i + 1) + j], wires[i]]) + + +def fourier_adder_phase_shift(a, wires): + """Sends QFT(|b>) -> QFT(|b + a>).""" + n = len(wires) + a_bits = jnp.unpackbits(jnp.array([a]).view("uint8"), bitorder="little")[:n][::-1] + powers_of_two = jnp.array([1 / (2**k) for k in range(1, n + 1)]) + phases = jnp.array([jnp.dot(a_bits[k:], powers_of_two[: n - k]) for k in range(n)]) + + for i in range(len(wires)): + if phases[i] != 0: + qml.PhaseShift(2 * jnp.pi * phases[i], wires=wires[i]) + + +def doubly_controlled_adder(N, a, control_wires, wires, aux_wire): + """Sends |c>|x>QFT(|b>)|0> -> |c>|x>QFT(|b + c x a) mod N>)|0>.""" + qml.ctrl(fourier_adder_phase_shift, control=control_wires)(a, wires) + + qml.adjoint(fourier_adder_phase_shift)(N, wires) + + qml.adjoint(QFT)(wires) + qml.CNOT(wires=[wires[0], aux_wire]) + QFT(wires) + + qml.ctrl(fourier_adder_phase_shift, control=aux_wire)(N, wires) + + qml.adjoint(qml.ctrl(fourier_adder_phase_shift, control=control_wires))(a, wires) + + qml.adjoint(QFT)(wires) + qml.PauliX(wires=wires[0]) + qml.CNOT(wires=[wires[0], aux_wire]) + qml.PauliX(wires=wires[0]) + QFT(wires) + + qml.ctrl(fourier_adder_phase_shift, control=control_wires)(a, wires) + + +def controlled_ua(N, a, control_wire, target_wires, aux_wires, mult_a_mask, mult_a_inv_mask): + """Sends |c>|x>|0> to |c>|ax mod N>|0> if c = 1. + + The mask arguments allow for the removal of unnecessary double-controlled additions. + """ + n = len(target_wires) + + # Apply double-controlled additions where bits of a can be 1. + for i in range(n): + if mult_a_mask[n - i - 1] > 0: + pow_a = (a * (2**i)) % N + doubly_controlled_adder( + N, pow_a, [control_wire, target_wires[n - i - 1]], aux_wires[:-1], aux_wires[-1] + ) + + qml.adjoint(QFT)(wires=aux_wires[:-1]) + + # Controlled SWAP the target and aux wires; note that the top-most aux wire + # is only to catch overflow, so we ignore it here. + for i in range(n): + qml.CNOT(wires=[aux_wires[i + 1], target_wires[i]]) + qml.Toffoli(wires=[control_wire, target_wires[i], aux_wires[i + 1]]) + qml.CNOT(wires=[aux_wires[i + 1], target_wires[i]]) + + # Adjoint of controlled multiplication with the modular inverse of a + a_mod_inv = modular_inverse(a, N) + + QFT(wires=aux_wires[:-1]) + + for i in range(n): + if mult_a_inv_mask[i] > 0: + pow_a_inv = (a_mod_inv * (2 ** (n - i - 1))) % N + qml.adjoint(doubly_controlled_adder)( + N, + pow_a_inv, + [control_wire, target_wires[i]], + aux_wires[:-1], + aux_wires[-1], + ) + + +###################################################################### +# Finally, let's put it all together in the core portion of Shor's algorithm, +# under the ``@qml.qjit`` decorator. + +from jax import random + + +# ``n_bits`` is a static argument because ``jnp.arange`` does not currently +# support dynamically-shaped arrays when jitted. +@qml.qjit(autograph=True, static_argnums=(3)) +def shors_algorithm(N, key, a, n_bits, n_trials): + # If no explicit a is passed (denoted by a = 0), randomly choose a + # non-trivial value of a that does not have a common factor with N. + if a == 0: + while jnp.gcd(a, N) != 1: + key, subkey = random.split(key) + a = random.randint(subkey, (1,), 2, N - 1)[0] + + est_wire = 0 + target_wires = jnp.arange(n_bits) + 1 + aux_wires = jnp.arange(n_bits + 2) + n_bits + 1 + + dev = qml.device("lightning.qubit", wires=2 * n_bits + 3, shots=1) + + @qml.qnode(dev) + def run_qpe(): + meas_results = jnp.zeros((n_bits,), dtype=jnp.int32) + cumulative_phase = jnp.array(0.0) + phase_divisors = 2.0 ** jnp.arange(n_bits + 1, 1, -1) + + a_mask = jnp.zeros(n_bits, dtype=jnp.int64) + a_mask = a_mask.at[0].set(1) + jnp.array( + jnp.unpackbits(jnp.array([a]).view("uint8"), bitorder="little")[:n_bits] + ) + a_inv_mask = a_mask + + # Initialize the target register of QPE in |1> + qml.PauliX(wires=target_wires[-1]) + + # The "first" QFT on the auxiliary register; required here because + # QFT are largely omitted in the modular arithmetic routines due to + # cancellation between adjacent blocks of the algorithm. + QFT(wires=aux_wires[:-1]) + + # First iteration: add a - 1 using the Fourier adder. + qml.Hadamard(wires=est_wire) + + QFT(wires=target_wires) + qml.ctrl(fourier_adder_phase_shift, control=est_wire)(a - 1, target_wires) + qml.adjoint(QFT)(wires=target_wires) + + # Measure the estimation wire and store the phase + qml.Hadamard(wires=est_wire) + meas_results[0] = measure(est_wire, reset=True) + cumulative_phase = -2 * jnp.pi * jnp.sum(meas_results / jnp.roll(phase_divisors, 1)) + + # For subsequent iterations, determine powers of a, and apply controlled + # U_a when the power is not 1. Unnecessary double-controlled operations + # are removed, based on values stored in the two "mask" variables. + powers_cua = jnp.array([repeated_squaring(a, 2**p, N) for p in range(n_bits)]) + + loop_bound = n_bits + if jnp.min(powers_cua) == 1: + loop_bound = jnp.argmin(powers_cua) + + # The core of the QPE routine + for pow_a_idx in range(1, loop_bound): + pow_cua = powers_cua[pow_a_idx] + + if not jnp.all(a_inv_mask): + for power in range(2**pow_a_idx, 2 ** (pow_a_idx + 1)): + next_pow_a = jnp.array([repeated_squaring(a, power, N)]) + a_inv_mask = a_inv_mask + jnp.array( + jnp.unpackbits(next_pow_a.view("uint8"), bitorder="little")[:n_bits] + ) + + qml.Hadamard(wires=est_wire) + + controlled_ua(N, pow_cua, est_wire, target_wires, aux_wires, a_mask, a_inv_mask) + + a_mask = a_mask + a_inv_mask + a_inv_mask = jnp.zeros_like(a_inv_mask) + + # Rotate the estimation wire by the accumulated phase, then measure and reset it + qml.PhaseShift(cumulative_phase, wires=est_wire) + qml.Hadamard(wires=est_wire) + meas_results[pow_a_idx] = measure(est_wire, reset=True) + cumulative_phase = ( + -2 * jnp.pi * jnp.sum(meas_results / jnp.roll(phase_divisors, pow_a_idx + 1)) + ) + + # The adjoint partner of the QFT from the beginning, to reset the auxiliary register + qml.adjoint(QFT)(wires=aux_wires[:-1]) + + return meas_results + + # The "classical" part of Shor's algorithm: run QPE, extract a candidate + # order, then check whether a valid solution is found. We run multiple trials, + # and return a success probability. + p, q = jnp.array(0, dtype=jnp.int32), jnp.array(0, dtype=jnp.int32) + successful_trials = jnp.array(0, dtype=jnp.int32) + + for _ in range(n_trials): + sample = run_qpe() + phase = fractional_binary_to_float(sample) + guess_r = phase_to_order(phase, N) + + # If the guess order is even, we may have a non-trivial square root. + # If so, try to compute p and q. + if guess_r % 2 == 0: + guess_square_root = repeated_squaring(a, guess_r // 2, N) + + if guess_square_root != 1 and guess_square_root != N - 1: + candidate_p = jnp.gcd(guess_square_root - 1, N).astype(jnp.int32) + + if candidate_p != 1: + candidate_q = N // candidate_p + else: + candidate_q = jnp.gcd(guess_square_root + 1, N).astype(jnp.int32) + + if candidate_q != 1: + candidate_p = N // candidate_q + + if candidate_p * candidate_q == N: + p, q = candidate_p, candidate_q + successful_trials += 1 + + return p, q, key, a, successful_trials / n_trials + + +###################################################################### +# Let's ensure the algorithm can successfully factor a small case, :math:`N = +# 15`. We will randomly generate :math:`a` within the function, and do 100 shots +# of the phase estimation procedure to get an idea of the success probability. + + +key = random.PRNGKey(123456789) + +N = 15 +n_bits = int(jnp.ceil(jnp.log2(N))) + +p, q, _, a, success_prob = shors_algorithm(N, key.astype(jnp.uint32), 0, n_bits, 100) + +print(f"Found {N} = {p} x {q} (using random a = {a}) with probability {success_prob:.2f}") + +###################################################################### +# Performance and validation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Next, let's verify QJIT compilation is happening properly. We will run the +# algorithm for different :math:`N` with the same bit width, and different values of +# :math:`a`. We expect the first execution, for the very first :math:`N` and :math:`a`, to +# take longer than the rest. + +import time +import matplotlib.pyplot as plt + +# Some 6-bit numbers +N_values = [33, 39, 51, 55, 57] +n_bits = int(jnp.ceil(jnp.log2(N_values[0]))) + +num_a = 3 + +execution_times = [] + +key = random.PRNGKey(1010101) + +for N in N_values: + unique_a = [] + + while len(unique_a) != num_a: + key, subkey = random.split(key.astype(jnp.uint32)) + a = random.randint(subkey, (1,), 2, N - 1)[0] + if jnp.gcd(a, N) == 1 and a not in unique_a: + unique_a.append(a) + + for a in unique_a: + start = time.time() + p, q, key, _, _ = shors_algorithm(N, key.astype(jnp.uint32), a, n_bits, 1) + end = time.time() + execution_times.append((N, a, end - start)) + + +labels = [f"{ex[0]}, {int(ex[1])}" for ex in execution_times] +times = [ex[2] for ex in execution_times] + +plt.scatter(range(len(times)), times, c=[ex[0] for ex in execution_times]) +plt.xticks(range(len(times)), labels=labels, rotation=80) +plt.xlabel("N, a") +plt.ylabel("Runtime (s)") +plt.tight_layout() +plt.show() + +###################################################################### +# This plot demonstrates exactly what we suspect: changing :math:`N` and +# :math:`a` does not lead to recompilation of the program! This will be +# particularly valuable for large :math:`N`, where traditional circuit +# processing times can grow very large. +# +# To show this more explicitly, let's fix :math:`a = 2`, and generate Shor +# circuits for many different :math:`N` using both the QJIT version, and the +# plain PennyLane version below. Note the standard PennyLane version makes use +# of many of the same subroutines and optimizations, but due to limitations on +# how PennyLane handles mid-circuit measurements, we must use ``qml.cond`` and +# explicit ``qml.PhaseShift`` gates. + + +def shors_algorithm_no_qjit(N, key, a, n_bits, n_trials): + est_wire = 0 + target_wires = list(range(1, n_bits + 1)) + aux_wires = list(range(n_bits + 1, 2 * n_bits + 3)) + + dev = qml.device("lightning.qubit", wires=2 * n_bits + 3, shots=1) + + @qml.qnode(dev) + def run_qpe(): + a_mask = jnp.zeros(n_bits, dtype=jnp.int64) + a_mask = a_mask.at[0].set(1) + jnp.array( + jnp.unpackbits(jnp.array([a]).view("uint8"), bitorder="little")[:n_bits] + ) + a_inv_mask = a_mask + + measurements = [] + + qml.PauliX(wires=target_wires[-1]) + + QFT(wires=aux_wires[:-1]) + + qml.Hadamard(wires=est_wire) + + QFT(wires=target_wires) + qml.ctrl(fourier_adder_phase_shift, control=est_wire)(a - 1, target_wires) + qml.adjoint(QFT)(wires=target_wires) + + qml.Hadamard(wires=est_wire) + measurements.append(qml.measure(est_wire, reset=True)) + + powers_cua = jnp.array([repeated_squaring(a, 2**p, N) for p in range(n_bits)]) + + loop_bound = n_bits + if jnp.min(powers_cua) == 1: + loop_bound = jnp.argmin(powers_cua) + + for pow_a_idx in range(1, loop_bound): + pow_cua = powers_cua[pow_a_idx] + + if not jnp.all(a_inv_mask): + for power in range(2**pow_a_idx, 2 ** (pow_a_idx + 1)): + next_pow_a = jnp.array([repeated_squaring(a, power, N)]) + a_inv_mask = a_inv_mask + jnp.array( + jnp.unpackbits(next_pow_a.view("uint8"), bitorder="little")[:n_bits] + ) + + qml.Hadamard(wires=est_wire) + + controlled_ua(N, pow_cua, est_wire, target_wires, aux_wires, a_mask, a_inv_mask) + + a_mask = a_mask + a_inv_mask + a_inv_mask = jnp.zeros_like(a_inv_mask) + + # The main difference with the QJIT version + for meas_idx, meas in enumerate(measurements): + qml.cond(meas, qml.PhaseShift)( + -2 * jnp.pi / 2 ** (pow_a_idx + 2 - meas_idx), wires=est_wire + ) + + qml.Hadamard(wires=est_wire) + measurements.append(qml.measure(est_wire, reset=True)) + + qml.adjoint(QFT)(wires=aux_wires[:-1]) + + return qml.sample(measurements) + + p, q = jnp.array(0, dtype=jnp.int32), jnp.array(0, dtype=jnp.int32) + successful_trials = jnp.array(0, dtype=jnp.int32) + + for _ in range(n_trials): + sample = jnp.array([run_qpe()]) + phase = fractional_binary_to_float(sample) + guess_r = phase_to_order(phase, N) + + if guess_r % 2 == 0: + guess_square_root = repeated_squaring(a, guess_r // 2, N) + + if guess_square_root != 1 and guess_square_root != N - 1: + candidate_p = jnp.gcd(guess_square_root - 1, N).astype(jnp.int32) + + if candidate_p != 1: + candidate_q = N // candidate_p + else: + candidate_q = jnp.gcd(guess_square_root + 1, N).astype(jnp.int32) + + if candidate_q != 1: + candidate_p = N // candidate_q + + if candidate_p * candidate_q == N: + p, q = candidate_p, candidate_q + successful_trials += 1 + + return p, q, key, a, successful_trials / n_trials + + +###################################################################### +# Let's do the same experiment as before, with three different choices of +# :math:`a` for each :math:`N`. + +execution_times_qjit = [] +execution_times_standard = [] + +key = random.PRNGKey(1010101) + +for N in N_values: + unique_a = [] + + while len(unique_a) != num_a: + key, subkey = random.split(key.astype(jnp.uint32)) + a = random.randint(subkey, (1,), 2, N - 1)[0] + if jnp.gcd(a, N) == 1 and a not in unique_a: + unique_a.append(a) + + for a in unique_a: + # QJIT times + start = time.time() + p, q, _, _, _ = shors_algorithm(N, key.astype(jnp.uint32), a, n_bits, 1) + end = time.time() + execution_times_qjit.append((N, a, end - start)) + + # No QJIT times + start = time.time() + p, q, _, _, _ = shors_algorithm_no_qjit(N, key.astype(jnp.uint32), a, n_bits, 1) + end = time.time() + execution_times_standard.append((N, a, end - start)) + + +labels = [f"{ex[0]}, {int(ex[1])}" for ex in execution_times_qjit] +colours = [ex[0] for ex in execution_times_qjit] + +times_qjit = [ex[2] for ex in execution_times_qjit] +times_standard = [ex[2] for ex in execution_times_standard] + + +plt.scatter(range(len(times)), times_qjit, c=colours, label="QJIT") +plt.scatter(range(len(times)), times_standard, c=colours, marker="v", label="Standard") +plt.xticks(range(0, len(times)), labels=labels, rotation=80) +plt.xlabel("N, a") +plt.ylabel("Runtime (s)") +plt.legend() +plt.tight_layout() +plt.show() + +###################################################################### +# Without QJIT, different values of :math:`a` for the same :math:`N` can have +# wildly different execution times! This is largely due to the +# :math:`a`-specific optimizations. When we use QJIT, we get the benefits of +# that optimization *and* comparable performance across any choice of :math:`a`. +# +# Finally, let's compare different values of N with same choice of :math:`a`. + +N_values = [15, 21, 33, 39, 51, 55, 57, 65] +execution_times_qjit = [] +execution_times_standard = [] + +for N in N_values: + start = time.time() + p, q, key, _, _ = shors_algorithm(N, key.astype(jnp.uint32), 2, n_bits, 1) + end = time.time() + execution_times_qjit.append(end - start) + + start = time.time() + p, q, key, _, _ = shors_algorithm_no_qjit(N, key.astype(jnp.uint32), 2, n_bits, 1) + end = time.time() + execution_times_standard.append(end - start) + +plt.scatter(range(len(N_values)), execution_times_qjit, label="QJIT") +plt.scatter(range(len(N_values)), execution_times_standard, marker="v", label="Standard") +plt.xticks(range(0, len(N_values)), labels=N_values, rotation=80) +plt.xlabel("N") +plt.ylabel("Runtime (s)") +plt.legend() +plt.tight_layout() +plt.show() + + +###################################################################### +# Here we observe that without QJIT, the runtime for different :math:`N`, +# even with the same bit width, may differ greatly. QJIT enables more consistent +# performance, and greatly benefits from reuse of the cached program. +# Preliminary experiments show this to be true even for larger problem sizes! +# +# Conclusions +# ----------- +# +# The ability to leverage a tool like Catalyst means we can quickly generate, +# compile, and optimize very large circuits, even within the context of a larger +# workflow. As a bonus, using JIT compilation means that after the first +# execution, these optimizations come at no extra cost, even though they depend +# on runtime values! +# +# There is still much work to be done, however. For one, the generated +# circuits are not optimized at the individual gate level, so the resource +# counts will be large and impractical. One also observes that, even though +# "higher-level" optimizations become possible, there is still a significant +# amount of manual labour involved in implementing them. Compilation tools that +# can automatically determine how subroutines can be simplified based on +# classical information or on the input quantum superposition would be valuable +# to develop, as they would enable co-optimization of the classical and quantum +# parts of workflows. +# +# *Acknowledgements*: The author thanks the Catalyst team for their support and +# for developing the features needed to carry out this project. Thanks +# especially to David Ittah for developing JIT-compatible implementations of +# utility functions needed for modular exponentiation and phase estimation. +# +# +# References +# ---------- +# +# .. [#Shor1997] +# +# Peter W. Shor (1997) *Polynomial-Time Algorithms for Prime Factorization +# and Discrete Logarithms on a Quantum Computer*. SIAM Journal on Computing, +# Vol. 26, Iss. 5. +# +# .. [#PurpleDragonBook] +# +# Alfred V Aho, Monica S Lam, Ravi Sethi, Jeffrey D Ullman (2007) +# *Compilers: Principles, Techniques, And Tools*. Pearson Education, Inc. +# +# .. [#JAXJIT] +# +# The JAX Authors (2024) *Just-in-time compilation.* `JAX Documentation +# `_. +# +# .. [#Beauregard2003] +# +# Stephane Beauregard (2003) *Circuit for Shor's algorithm using 2n+3 qubits.* +# Quantum Information and Computation, Vol. 3, No. 2 (2003) pp. 175-185. +# +# .. [#Draper2000] +# +# Thomas G. Draper (2000) *Addition on a Quantum Computer.* arXiv preprint, +# arXiv:quant-ph/0008033. +# +# +# Appendix: fun with Fourier transforms +# ------------------------------------- +# +# This appendix provides a detailed derivation of the circuits for addition and +# subtraction in the Fourier basis (both the generic and the modulo :math:`N` +# case), as well as the full explanation of the one-qubit QPE trick. +# +# .. _appendix_fourier_adder: +# +# Standard Fourier addition +# ^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# We recall the circuit from the main text. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder.svg +# :width: 700 +# :align: center +# :alt: Addition in the Fourier basis. +# +# Circuit for addition in the Fourier basis [#Draper2000]_. +# +# In the third circuit, we've absorbed the Fourier transforms into the basis +# states, and denoted this by calligraphic letters. The :math:`\mathbf{R}_k` are +# phase shifts based on :math:`a`, and will be described below. +# +# To understand how Fourier addition works, we begin by revisiting the QFT. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_explanation-1.svg +# :width: 800 +# :align: center +# :alt: The Quantum Fourier Transform. +# +# Circuit for the quantum Fourier transform. Note the big-endian qubit ordering with +# no terminal SWAP gates. +# +# In this circuit we define the phase gate +# +# .. math:: +# +# R_k = \begin{pmatrix} 1 & 0 \\ 0 & e^{\frac{2\pi i}{2^k}} \end{pmatrix}. +# +# The qubits are ordered using big endian notation. For an :math:`n`-bit integer +# :math:`b`, :math:`\vert b\rangle = \vert b_{n-1} \cdots b_0\rangle` and +# :math:`b = \sum_{k=0}^{n-1} 2^k b_k`. +# +# Suppose we wish to add :math:`a` to :math:`b`. We can add a new register +# prepared in :math:`\vert a \rangle`, and use its qubits to control the +# addition of phases to qubits in :math:`\vert b \rangle` (after a QFT is +# applied) in a very particular way: +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_explanation-2.svg +# :width: 800 +# :align: center +# :alt: Adding one integer to another with the Quantum Fourier Transform. +# +# Fourier addition of :math:`a` to :math:`b`. The bit values of :math:`a` +# determine the amount of phase added to the qubits in :math:`b` [#Draper2000]_. +# +# Each qubit in :math:`\vert b \rangle` picks up a phase that depends on the +# bits in :math:`a`. In particular, the :math:`k`'th bit of :math:`b` +# accumulates information about all the bits in :math:`a` with an equal or lower +# index, :math:`a_0, \ldots, a_{k}`. This adds :math:`a_k` to :math:`b_k`, and +# the cumulative effect adds :math:`a` to :math:`b`, up to an inverse QFT. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_explanation-3.svg +# :width: 800 +# :align: center +# :alt: Adding one integer to another with the Quantum Fourier Transform. +# +# However, we must be careful. Fourier basis addition is *not* automatically +# modulo :math:`N`. If the sum :math:`b + a` requires :math:`n+ 1` bits, it will +# overflow. To handle that, one extra qubit is added to the :math:`\vert +# b\rangle` register (initialized to :math:`\vert 0 \rangle`). This is the +# source of one of the auxiliary qubits mentioned earlier. +# +# In our case, a second register of qubits is not required. Since we know +# :math:`a` in advance, we can precompute the amount of phase to apply: qubit +# :math:`\vert b_k \rangle` must be rotated by :math:`\sum_{\ell=0}^{k} +# \frac{a_\ell}{2^{\ell+1}}`. We'll express this as a new gate, +# +# .. math:: +# +# \mathbf{R}_k = \begin{pmatrix} 1 & 0 \\ 0 & e^{2\pi i\sum_{\ell=0}^{k} \frac{a_\ell}{2^{\ell+1}}} \end{pmatrix}. +# +# The final circuit for the Fourier adder is +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_explanation-4.svg +# :width: 500 +# :align: center +# :alt: Full Fourier adder. +# +# +# As one may expect, :math:`\Phi^\dagger` performs subtraction. However, we must +# also consider the possibility of underflow. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_adjoint.svg +# :width: 500 +# :align: center +# :alt: Subtraction in the Fourier basis. +# +# .. _appendix_fourier_adder_modulo_n: +# +# Doubly-controlled Fourier addition modulo :math:`N` +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Let's analyze below circuit, keeping in mind :math:`a < N` and :math:`b < N`, +# and the main register has :math:`n + 1` bits. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/fourier_adder_modulo_n.svg +# :width: 800 +# :align: center +# :alt: Addition in the Fourier basis modulo N. +# +# Circuit for doubly-controlled Fourier addition modulo :math:`N` +# [#Beauregard2003]_. Calligraphic letters indicate states already +# transformed to the Fourier basis. +# +# Suppose both control qubits are :math:`\vert 1 \rangle`. We add :math:`a` to +# :math:`b`, then subtract :math:`N`. If :math:`a + b \geq N`, we get overflow +# on the top qubit (which was included to account for precisely this), but +# subtraction gives us the correct result modulo :math:`N`. After shifting back +# to the computational basis with the :math:`QFT^\dagger`, the topmost qubit is +# in :math:`\vert 0 \rangle` so the CNOT does not trigger. Next, we subtract +# :math:`a` from the register, in state :math:`\vert a + b - N \pmod N \rangle`, +# to obtain :math:`\vert b - N \rangle`. Since :math:`b < N`, there is +# underflow. The top qubit, now in state :math:`\vert 1 \rangle`, does not +# trigger the controlled-on-zero CNOT, and the auxiliary qubit is untouched. +# +# If instead :math:`a + b < N`, we subtracted :math:`N` for no reason, leading +# to underflow. The topmost qubit is :math:`\vert 1 \rangle`, the CNOT will +# trigger, and :math:`N` gets added back. The register then contains +# :math:`\vert b + a \rangle`. Subtracting :math:`a` puts the register in state +# :math:`\vert b \rangle`, which by design will not overflow; the +# controlled-on-zero CNOT triggers, and the auxiliary qubit is returned to +# :math:`\vert 0 \rangle`. +# +# If the control qubits are not both :math:`\vert 1 \rangle`, :math:`N` is +# subtracted and the CNOT triggers (since :math:`b < N`), but :math:`N` is +# always added back. By similar reasoning, the controlled-on-zero CNOT always +# triggers and correctly uncomputes the auxiliary qubit. +# +# Recall that :math:`\Phi_+` is used in :math:`M_a` to add :math:`2^{k}a` to +# :math:`b` (controlled on :math:`x_{k}`) in the Fourier basis. In total, we +# obtain (modulo :math:`N`) +# +# .. math:: +# +# \begin{equation*} +# b + x_{0} \cdot 2^0 a + x_{1} \cdot 2^1 a + \cdots x_{n-1} \cdot 2^{n-1} a = b + a \sum_{k=0}^{n-1} x_{k} 2^k = b + a x. +# \end{equation*} +# +# .. _appendix_single_qubit_qpe: +# +# QPE with one estimation wire +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +# Here we show how the inverse QFT in QPE can be implemented using a single +# estimation wire when mid-circuit measurement and feedforward are available. +# Our starting point is the full algorithm below. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded. +# +# Look carefully at the qubit in :math:`\vert \theta_0\rangle`. After the final +# Hadamard, it is used only for controlled gates. Thus, we can just measure it +# and apply subsequent operations controlled on the classical outcome, +# :math:`\theta_0`. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-2.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded and last estimation qubit measured off. +# +# However, we can do better if we can directly modify the circuit based on the +# measurement outcomes. Instead of applying controlled :math:`R^\dagger_2`, we +# can apply :math:`R^\dagger` where the rotation angle is 0 if :math:`\theta_0 = +# 0`, and :math:`-2\pi i/2^2` if :math:`\theta_0 = 1`, i.e., :math:`R^\dagger_{2 +# \theta_0}`. The same can be done for all other gates controlled on +# :math:`\theta_0`. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-3.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded, last estimation qubit measured, and rotation gates adjusted. +# +# We'll leverage this trick again with the second-last estimation +# qubit. Moreover, we can make a further improvement by noting that once the +# last qubit is measured, we can reset and repurpose it to play the role of the +# second-last qubit. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-4.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded and last estimation qubit reused. +# +# Next, we adjust rotation angles based on measurement values, removing +# the need for classical controls. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-5.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with inverse QFT expanded, last estimation qubit reused, and rotation gates adjusted. +# +# We can do this for all remaining estimation qubits, adding more rotations +# depending on previous measurement outcomes. +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-6.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with one estimation qubit, and unmerged rotation gates. +# +# Finally, since these are all :math:`RZ`, we can merge them. Let +# +# .. math:: +# +# \mathbf{M}_{k} = \begin{pmatrix} 1 & 0 \\ 0 & e^{-2\pi i\sum_{\ell=0}^{k} \frac{\theta_{\ell}}{2^{k + 2 - \ell}}} \end{pmatrix}. +# +# With a bit of index gymnastics, we obtain our final QPE algorithm with a single estimation qubit: +# +# .. figure:: ../_static/demonstration_assets/shor_catalyst/qpe_full_modified_power_with_qft-7.svg +# :width: 800 +# :align: center +# :alt: QPE circuit with one estimation qubit. +# +# About the author +# ---------------- +# diff --git a/demonstrations_v2/tutorial_shors_algorithm_catalyst/metadata.json b/demonstrations_v2/tutorial_shors_algorithm_catalyst/metadata.json new file mode 100644 index 0000000000..92b87a75fe --- /dev/null +++ b/demonstrations_v2/tutorial_shors_algorithm_catalyst/metadata.json @@ -0,0 +1,98 @@ +{ + "title": "Quantum just-in-time compiling Shor's algorithm with Catalyst", + "authors": [ + { + "username": "glassnotes" + } + ], + "dateOfPublication": "2025-04-04T09:00:00+00:00", + "dateOfLastModification": "2025-04-04T09:00:00+00:00", + "categories": [ + "Algorithms", + "Quantum Computing", + "Devices and Performance" + ], + "tags": [], + "previewImages": [ + { + "type": "thumbnail", + "uri": "/_static/demo_thumbnails/regular_demo_thumbnails/thumbnail_shor_algorithm_catalyst_pennylane.png" + }, + { + "type": "large_thumbnail", + "uri": "/_static/demo_thumbnails/large_demo_thumbnails/thumbnail_large_shor_algorithm_catalyst_pennylane.png" + } + ], + "seoDescription": "Just-in-time compile Shor's algorithm from end-to-end with PennyLane and Catalyst.", + "doi": "", + "references": [ + { + "id": "Shor1997", + "type": "article", + "title": "Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer", + "authors": "Peter W. Shor", + "year": "1997", + "journal": "SIAM Journal on Computing Vol. 26, Iss. 5", + "doi": "10.1137/S0097539795293172", + "url": "https://doi.org/10.1137/S0097539795293172" + }, + { + "id": "PurpleDragonBook", + "type": "book", + "title": "Compilers: Principles, Techniques, And Tools", + "authors": "Alfred V Aho, Monica S Lam, Ravi Sethi, Jeffrey D Ullman", + "year": "2007", + "publisher": "Pearson", + "url": "https://www.pearson.com/en-us/subject-catalog/p/compilers-principles-techniques-and-tools/P200000003472/9780133002140" + }, + { + "id": "JAXJIT", + "type": "webpage", + "title": "Just-in-time compilation", + "authors": "The JAX Authors", + "year": "2024", + "publisher": "", + "url": "https://docs.jax.dev/en/latest/jit-compilation.html" + }, + { + "id": "Beauregard2003", + "type": "article", + "title": "Circuit for Shor's algorithm using 2n+3 qubits", + "authors": "Stephane Beauregard", + "year": "2003", + "journal": "Quantum Info. Comput.", + "doi": "10.5555/2011517.2011525", + "url": "https://dl.acm.org/doi/10.5555/2011517.2011525" + }, + { + "id": "Draper2000", + "type": "preprint", + "title": "Addition on a Quantum Computer", + "authors": "Thomas G. Draper", + "year": "2000", + "journal": "", + "doi": "10.48550/arXiv.quant-ph/0008033", + "url": "https://doi.org/10.48550/arXiv.quant-ph/0008033" + } + ], + "basedOnPapers": [ + ], + "referencedByPapers": [], + "relatedContent": [ + { + "type": "demonstration", + "id": "tutorial_jax_transformations", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_qjit_compile_grovers_algorithm_with_catalyst", + "weight": 1.0 + }, + { + "type": "demonstration", + "id": "tutorial_how_to_quantum_just_in_time_compile_vqe_catalyst", + "weight": 1.0 + } + ] +} diff --git a/demonstrations_v2/tutorial_testing_symmetry/demo.py b/demonstrations_v2/tutorial_testing_symmetry/demo.py index 61e0776885..802ac04055 100644 --- a/demonstrations_v2/tutorial_testing_symmetry/demo.py +++ b/demonstrations_v2/tutorial_testing_symmetry/demo.py @@ -308,14 +308,14 @@ def prep_plus(): # Implement controlled symmetry operations on system def CU_sys(): - qml.ControlledQubitUnitary(c_mat @ c_mat, control_wires=[aux[0]], wires=system) - qml.ControlledQubitUnitary(c_mat, control_wires=[aux[1]], wires=system) + qml.ControlledQubitUnitary(c_mat @ c_mat, wires=[aux[0]] + list(system)) + qml.ControlledQubitUnitary(c_mat, wires=[aux[1]] + list(system)) # Implement controlled symmetry operations on copy def CU_cpy(): - qml.ControlledQubitUnitary(c_mat @ c_mat, control_wires=[aux[0]], wires=copy) - qml.ControlledQubitUnitary(c_mat, control_wires=[aux[1]], wires=copy) + qml.ControlledQubitUnitary(c_mat @ c_mat, wires=[aux[0]] + list(copy)) + qml.ControlledQubitUnitary(c_mat, wires=[aux[1]] + list(copy)) ###################################################################### # Let’s combine everything and actually run our circuit! diff --git a/demonstrations_v2/tutorial_testing_symmetry/metadata.json b/demonstrations_v2/tutorial_testing_symmetry/metadata.json index 2a7ca1c3df..b42d1b3179 100644 --- a/demonstrations_v2/tutorial_testing_symmetry/metadata.json +++ b/demonstrations_v2/tutorial_testing_symmetry/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2023-01-24T00:00:00+00:00", - "dateOfLastModification": "2024-11-06T00:00:00+00:00", + "dateOfLastModification": "2025-01-23T00:00:00+00:00", "categories": [ "Algorithms", "Quantum Computing" diff --git a/demonstrations_v2/tutorial_variational_classifier/demo.py b/demonstrations_v2/tutorial_variational_classifier/demo.py index a979f5011f..3255bc29a2 100644 --- a/demonstrations_v2/tutorial_variational_classifier/demo.py +++ b/demonstrations_v2/tutorial_variational_classifier/demo.py @@ -128,7 +128,6 @@ def circuit(weights, x): def variational_classifier(weights, bias, x): return circuit(weights, x) + bias - ############################################################################## # Cost # ~~~~ diff --git a/demonstrations_v2/tutorial_variational_classifier/metadata.json b/demonstrations_v2/tutorial_variational_classifier/metadata.json index 3941e9212a..ccdaee1c42 100644 --- a/demonstrations_v2/tutorial_variational_classifier/metadata.json +++ b/demonstrations_v2/tutorial_variational_classifier/metadata.json @@ -39,5 +39,12 @@ "weight": 1.0 } ], + "hardware": [ + { + "id": "qbraid", + "link": "https://account.qbraid.com?gitHubUrl=https://github.com/PennyLaneAI/pennylane-demo-notebooks.git", + "logo": "/_static/hardware_logos/qbraid.png" + } + ], "discussionForumUrl": "https://discuss.pennylane.ai/t/variational-classifier-demo-question-on-padding/3367" } diff --git a/demonstrations_v2/tutorial_vqe/demo.py b/demonstrations_v2/tutorial_vqe/demo.py index baba6119e6..106fdf4de2 100644 --- a/demonstrations_v2/tutorial_vqe/demo.py +++ b/demonstrations_v2/tutorial_vqe/demo.py @@ -108,7 +108,7 @@ # where :math:`\theta` is the variational parameter to be optimized in order to find # the best approximation to the true ground state. In the Jordan-Wigner [#seeley2012]_ encoding, # the first term :math:`|1100\rangle` represents the `Hartree-Fock (HF) state -# `_ where the two electrons in +# `_ where the two electrons in # the molecule occupy the lowest-energy orbitals. The second term :math:`|0011\rangle` # encodes a double excitation of the HF state where the two particles are excited from # qubits 0, 1 to 2, 3. diff --git a/demonstrations_v2/tutorial_vqe/metadata.json b/demonstrations_v2/tutorial_vqe/metadata.json index 8cc82fc5dd..463adbe8df 100644 --- a/demonstrations_v2/tutorial_vqe/metadata.json +++ b/demonstrations_v2/tutorial_vqe/metadata.json @@ -73,6 +73,11 @@ "id": "aws", "link": "https://github.com/amazon-braket/amazon-braket-examples/blob/main/examples/pennylane/3_Hydrogen_Molecule_geometry_with_VQE/3_Hydrogen_Molecule_geometry_with_VQE.ipynb", "logo": "/_static/hardware_logos/aws.png" + }, + { + "id": "qbraid", + "link": "https://account.qbraid.com?gitHubUrl=https://github.com/PennyLaneAI/pennylane-demo-notebooks.git", + "logo": "/_static/hardware_logos/qbraid.png" } ], "discussionForumUrl": "https://discuss.pennylane.ai/t/a-brief-overview-of-vqe-demo/7333" diff --git a/demonstrations_v2/tutorial_vqe_vqd/requirements.in b/demonstrations_v2/tutorial_vqe_vqd/requirements.in index 7555f688f2..58dcc7fe96 100644 --- a/demonstrations_v2/tutorial_vqe_vqd/requirements.in +++ b/demonstrations_v2/tutorial_vqe_vqd/requirements.in @@ -2,4 +2,3 @@ jax jaxlib optax torch -pennylane==0.39.0 diff --git a/demonstrations_v2/tutorial_zx_calculus/demo.py b/demonstrations_v2/tutorial_zx_calculus/demo.py index a63330a5eb..38aef427ba 100644 --- a/demonstrations_v2/tutorial_zx_calculus/demo.py +++ b/demonstrations_v2/tutorial_zx_calculus/demo.py @@ -547,6 +547,18 @@ def circuit(): manager.canvas.figure = fig fig.set_canvas(manager.canvas) +# Access axes and modify limits +ax = fig.axes[0] +xlim = ax.get_xlim() +ylim = ax.get_ylim() + +# Expand limits by 10% to "zoom out" +zoom_factor = 0.2 +x_margin = (xlim[1] - xlim[0]) * zoom_factor +y_margin = (ylim[1] - ylim[0]) * zoom_factor +ax.set_xlim(xlim[0] - x_margin, xlim[1] + x_margin) +ax.set_ylim(ylim[0] - y_margin, ylim[1] + y_margin) + plt.show() ############################################################################# diff --git a/demonstrations_v2/tutorial_zx_calculus/metadata.json b/demonstrations_v2/tutorial_zx_calculus/metadata.json index f646f23511..6aca7c544f 100644 --- a/demonstrations_v2/tutorial_zx_calculus/metadata.json +++ b/demonstrations_v2/tutorial_zx_calculus/metadata.json @@ -6,7 +6,7 @@ } ], "dateOfPublication": "2023-06-06T00:00:00+00:00", - "dateOfLastModification": "2024-10-07T00:00:00+00:00", + "dateOfLastModification": "2025-04-10T00:00:00+00:00", "categories": [ "Quantum Computing" ], diff --git a/dependencies/README.md b/dependencies/README.md new file mode 100644 index 0000000000..ecef3d624e --- /dev/null +++ b/dependencies/README.md @@ -0,0 +1,8 @@ +# Dependency Specification + +This directory contains files specifying constraints on demo dependencies. + +- [constraints-dev.txt](./constraints-dev.txt) specifies dependency versions for the most recent development builds of Pennylane +- [constraints-stable.txt](./constraints-stable.txt) specifies dependency versions for the most recent stable version of Pennylane +- [requirements-build.txt](./requirements-build.txt) specifies the dependencies for building demo notebooks using sphinx +- [requirements-core.txt](./requirements-core.in) specifies the execution dependencies that are automatically installed for all demos diff --git a/dependencies/constraints-dev.txt b/dependencies/constraints-dev.txt new file mode 100644 index 0000000000..8623e624fb --- /dev/null +++ b/dependencies/constraints-dev.txt @@ -0,0 +1,35 @@ +# Specifies global constraints for the dev build + +--extra-index-url https://test.pypi.org/simple/ +--extra-index-url https://download.pytorch.org/whl/cpu + +pennylane @ git+https://github.com/PennyLaneAI/pennylane.git#egg=pennylane +pennylane-qiskit @ git+https://github.com/PennyLaneAI/pennylane-qiskit.git#egg=pennylane-qiskit +pennylane-qulacs @ git+https://github.com/PennyLaneAI/pennylane-qulacs.git#egg=pennylane-qulacs +iqpot @ git+https://github.com/XanaduAI/iqpopt.git#egg=iqpopt +pennylane-catalyst +pennylane-lightning + +chex<1.0.0 +pyqrack==1.32.12 +numpy~=1.24 +matplotlib==3.7.2 +jax==0.4.28 +jaxlib==0.4.28 +jaxopt==0.8.3 +aiohttp==3.9.5 +fsspec==2024.6.1 +h5py==3.11.0 +openfermionpyscf==0.5 +openqaoa-core==0.2.5 +qiskit>=1.0.0 +qiskit-aer>=0.14.0,<0.16.0 +qiskit_ibm_runtime==0.29.0 +torch==2.1.2+cpu ; sys_platform != 'darwin' +torch==2.1.2 ; sys_platform == 'darwin' +torchvision==0.16.2+cpu ; sys_platform != 'darwin' +torchvision==0.16.2 ; sys_platform == 'darwin' +tensorflow==2.14.1 +optax==0.2.3 +quimb==1.8.2 +kahypar==1.1.7 diff --git a/dependencies/constraints-stable.txt b/dependencies/constraints-stable.txt new file mode 100644 index 0000000000..519f115dc5 --- /dev/null +++ b/dependencies/constraints-stable.txt @@ -0,0 +1,32 @@ +# Specifies global constraints for the stable build +--extra-index-url https://download.pytorch.org/whl/cpu +pennylane==0.41.0 +pennylane-cirq==0.41.0 +pennylane-qiskit==0.41.0 +pennylane-qulacs==0.41.0 +pennylane-catalyst==0.11.0 +pennylane-qrack==0.11.1 + +chex<1.0.0 +pyqrack==1.32.12 +numpy~=1.24 +matplotlib==3.7.2 +jax==0.4.28 +jaxlib==0.4.28 +jaxopt==0.8.3 +aiohttp==3.9.5 +fsspec==2024.6.1 +h5py==3.11.0 +openfermionpyscf==0.5 +openqaoa-core==0.2.5 +qiskit>=1.0.0 +qiskit-aer>=0.14.0,<0.16.0 +qiskit_ibm_runtime==0.29.0 +torch==2.1.2+cpu ; sys_platform != 'darwin' +torch==2.1.2 ; sys_platform == 'darwin' +torchvision==0.16.2+cpu ; sys_platform != 'darwin' +torchvision==0.16.2 ; sys_platform == 'darwin' +tensorflow==2.14.1 +optax==0.2.3 +quimb==1.8.2 +kahypar==1.1.7 diff --git a/dependencies/requirements-build.txt b/dependencies/requirements-build.txt new file mode 100644 index 0000000000..b06e7560fe --- /dev/null +++ b/dependencies/requirements-build.txt @@ -0,0 +1,11 @@ +# Build dependencies +sphinx>=5.0.2 +sphinx_gallery==0.17.1 +Jinja2==3.0.3 +markupsafe==2.1.1 +pyyaml~=6.0.1 +pypandoc==1.5 +pennylane-sphinx-theme @ git+https://github.com/PennyLaneAI/pennylane-sphinx-theme.git@sphinx-update +uv~=0.5 +numpy~=1.24 +pennylane diff --git a/dependencies/requirements-core.in b/dependencies/requirements-core.in new file mode 100644 index 0000000000..54a02bf248 --- /dev/null +++ b/dependencies/requirements-core.in @@ -0,0 +1,13 @@ +# Core dependencies that are always available for every demo +# Note that versions are specified in 'constraints-dev' or 'constraints-stable' +aiohttp +fsspec +h5py +jax +jaxlib +matplotlib +numpy +pennylane +pennylane-lightning +pennylane-catalyst +optax diff --git a/lib/qml/app/app.py b/lib/qml/app/app.py index c40eb4787c..115aa93059 100644 --- a/lib/qml/app/app.py +++ b/lib/qml/app/app.py @@ -1,12 +1,14 @@ import typer from qml.context import Context -from qml.lib import demo, repo +from qml.lib import demo, repo, cli, fs, template import shutil import logging from typing import Annotated, Optional import typing -from pathlib import Path - +import inflection +import re +import json +import rich logging.basicConfig(level=logging.INFO) @@ -36,10 +38,7 @@ def build( keep_going: Annotated[ bool, typer.Option(help="Continue if sphinx-build fails for a demo") ] = False, - overrides_file: Annotated[ - Optional[str], - typer.Option(help="Requirements file containing dependency version overrides"), - ] = None, + dev: Annotated[bool, typer.Option(help="Whether to use dev dependencies")] = False, ) -> None: """Build the named demos.""" ctx = Context() @@ -57,7 +56,73 @@ def build( execute=execute, quiet=quiet, keep_going=keep_going, - overrides_file=Path(overrides_file) if overrides_file else None, + dev=dev, + ) + + +@app.command() +def new(): + """Create a new demo.""" + ctx = Context() + title: str = typer.prompt("Title") + name_default = re.sub(r"\s+", "_", inflection.underscore(title)) + if not name_default.startswith("tutorial_"): + name_default = "tutorial_" + name_default + + while True: + name: str = typer.prompt("Name", name_default) + if demo.get(ctx.demos_dir, name): + print(f"Demo with name '{name}' already exists") + else: + break + + description = typer.prompt("Description", "") + + author_prompt = "Author's pennylane.ai username" + authors: list[str] = [] + authors.append(typer.prompt(author_prompt)) + + while True: + if not typer.confirm("Add another author?"): + break + + authors.append(typer.prompt(author_prompt)) + + small_thumbnail = cli.prompt_path("Thumbnail image") + large_thumbnail = cli.prompt_path("Large thumbnail image") + + demo_dir = ctx.demos_dir / name + demo_dir.mkdir(exist_ok=True) + + if small_thumbnail: + dest = (demo_dir / "thumbnail").with_suffix(small_thumbnail.suffix) + fs.copy_parents(small_thumbnail, dest) + small_thumbnail = str(dest.relative_to(demo_dir)) + + if large_thumbnail: + dest = (demo_dir / "large_thumbnail").with_suffix(large_thumbnail.suffix) + fs.copy_parents(large_thumbnail, dest) + large_thumbnail = str(dest.relative_to(demo_dir)) + + with open(demo_dir / "demo.py", "w") as f: + f.write(template.demo(title)) + + with open(demo_dir / "metadata.json", "w") as f: + json.dump( + template.metadata( + title=title, + description=description, + authors=authors, + thumbnail=small_thumbnail, + large_thumbnail=large_thumbnail, + categories=[], + ), + f, + indent=2, + ) + + rich.print( + f"Created new demo in [bold]{demo_dir.relative_to(ctx.repo_root)}[/bold]" ) diff --git a/lib/qml/context.py b/lib/qml/context.py index 7442f36ce8..0c903100b6 100644 --- a/lib/qml/context.py +++ b/lib/qml/context.py @@ -38,5 +38,17 @@ def repo(self) -> Repo: return Repo.discover() @property - def constraints_file(self) -> Path: - return self.repo_root / "constraints.txt" + def core_requirements_file(self) -> Path: + return self.repo_root / "dependencies" / "requirements-core.in" + + @property + def stable_constraints_file(self) -> Path: + return self.repo_root / "dependencies" / "constraints-stable.txt" + + @property + def dev_constraints_file(self) -> Path: + return self.repo_root / "dependencies" / "constraints-dev.txt" + + @property + def build_requirements_file(self) -> Path: + return self.repo_root / "dependencies" / "requirements-build.txt" diff --git a/lib/qml/lib/cli.py b/lib/qml/lib/cli.py new file mode 100644 index 0000000000..6a1a49a220 --- /dev/null +++ b/lib/qml/lib/cli.py @@ -0,0 +1,27 @@ +import typer +from pathlib import Path +import rich + + +def prompt_path(prompt: str) -> Path | None: + """Prompt user for a file path. + + Args: + prompt: Prompt text + + Returns: + Path: resolved path + None: User provided no inputs + """ + + path = typer.prompt(prompt, "") + if not path: + return None + + path = Path(path).expanduser() + if not path.exists(): + rich.print(f"File [red][bold]{path}[/bold][/red] does not exst.") + + return prompt_path(prompt) + + return path diff --git a/lib/qml/lib/cmds.py b/lib/qml/lib/cmds.py index 73907c5712..5b563bff84 100644 --- a/lib/qml/lib/cmds.py +++ b/lib/qml/lib/cmds.py @@ -51,6 +51,7 @@ def pip_install( constraints: Path | None = None, quiet: bool = True, use_uv: bool = True, + pre: bool = False, ): """Executes `pip install` with the given python interpreter and args. @@ -61,6 +62,8 @@ def pip_install( requirements: Path to a requirements file constraints: Path to a constriants file quiet: Whether to suppress output to stdout + pre: Include pre-release versions of pacakges + Raises: CalledProcessError: The command does not complete successfully @@ -85,6 +88,56 @@ def pip_install( cmd.extend(("--constraint", str(constraints))) if quiet: cmd.append("--quiet") + if pre: + cmd.append("--pre") cmd.extend(str(arg) for arg in args) subprocess.run(cmd).check_returncode() + + +def pip_compile( + python: str | Path, + output_file: Path, + *args: str | Path, + constraints_files: Iterable[Path] | None = None, + quiet: bool = True, + prerelease: bool = False, +) -> None: + """Execute `uv pip compile` with the given python interpreter. + + Args: + python: Path to python interpreter + output_file: Path to output file + args: Input requirements files, or extra arguments to pass + constraints_file: Path to constraints file + quiet: Whether to run in quiet mode + """ + cmd = [ + str(python), + "-m", + "uv", + "pip", + "compile", + "--index-strategy", + "unsafe-best-match", + "--no-header", + "--no-strip-extras", + "--no-strip-markers", + "--no-annotate", + "--emit-index-url", + "--output-file", + str(output_file), + ] + if quiet: + cmd.append("--quiet") + + if constraints_files: + for constraints in constraints_files: + cmd.extend(("--constraints", str(constraints.resolve()))) + + if prerelease: + cmd.append("--prerelease=allow") + + cmd.extend((str(arg) for arg in args)) + + subprocess.run(cmd).check_returncode() diff --git a/lib/qml/lib/demo.py b/lib/qml/lib/demo.py index c0d34c9caf..cde3437108 100644 --- a/lib/qml/lib/demo.py +++ b/lib/qml/lib/demo.py @@ -4,7 +4,6 @@ import shutil from qml.lib import fs, cmds from qml.lib.virtual_env import Virtualenv -from qml.lib.pip_tools import RequirementsGenerator import os import sys from logging import getLogger @@ -99,6 +98,16 @@ def requirements(self) -> frozenset[str]: return frozenset(reqs) +def get(search_dir: Path, name: str) -> Demo | None: + """Get demo with `name`, if it exists.""" + demo = Demo(name=name, path=search_dir / name) + + if not demo.py_file.exists(): + return None + + return demo + + def find(search_dir: Path, *names: str) -> Iterator[Demo]: """Find demos with given names in `search_dir`.""" if not names: @@ -132,7 +141,7 @@ def build( execute: bool, quiet: bool = False, keep_going: bool = False, - overrides_file: Path | None = None, + dev: bool = False, ) -> None: """Build the provided demos using 'sphinx-build', optionally executing them to generate plots and cell outputs. @@ -151,12 +160,8 @@ def build( logger.info("Building %d demos", len(demos)) build_venv = Virtualenv(ctx.build_venv_path) - _install_build_dependencies(build_venv, ctx.build_dir) - - requirements_generator = RequirementsGenerator( - Path(sys.executable), - global_constraints_file=ctx.constraints_file, - overrides_file=overrides_file, + cmds.pip_install( + build_venv.python, requirements=ctx.build_requirements_file, use_uv=False ) for demo in demos: @@ -172,15 +177,14 @@ def build( try: _build_demo( - sphinx_dir=ctx.repo_root, - build_dir=ctx.build_dir, + ctx, build_venv=build_venv, - requirements_generator=requirements_generator, target=target, execute=execute_demo, demo=demo, package=target is BuildTarget.JSON, quiet=quiet, + dev=dev, ) except subprocess.CalledProcessError as exc: if not keep_going: @@ -200,29 +204,52 @@ def build( raise RuntimeError(f"Failed to build {len(failed)} demos", failed) +def generate_requirements( + ctx: Context, demo: Demo, dev: bool, output_file: Path +) -> None: + constraints = [ctx.build_requirements_file] + if dev: + constraints.append(ctx.dev_constraints_file) + else: + constraints.append(ctx.stable_constraints_file) + + requirements_in = [ctx.core_requirements_file] + if demo.requirements_file: + requirements_in.append(demo.requirements_file) + + cmds.pip_compile( + sys.executable, + output_file, + *requirements_in, + constraints_files=constraints, + prerelease=dev, + ) + + def _build_demo( - sphinx_dir: Path, - build_dir: Path, + ctx: Context, build_venv: Virtualenv, demo: Demo, target: BuildTarget, - requirements_generator: "RequirementsGenerator", execute: bool, package: bool, quiet: bool, + dev: bool, ): - out_dir = sphinx_dir / "demos" / demo.name + out_dir = ctx.repo_root / "demos" / demo.name fs.clean_dir(out_dir) - with open(out_dir / "requirements.txt", "w") as f: - f.write(requirements_generator.generate_requirements(demo.requirements)) - + generate_requirements(ctx, demo, dev, out_dir / "requirements.txt") if execute: cmds.pip_install( - build_venv.python, requirements=out_dir / "requirements.txt", quiet=True + build_venv.python, + "--upgrade", + requirements=out_dir / "requirements.txt", + quiet=True, + pre=dev, ) - stage_dir = build_dir / "demonstrations" + stage_dir = ctx.build_dir / "demonstrations" fs.clean_dir(stage_dir) # Need a 'GALLERY_HEADER' file for sphinx-gallery with open(stage_dir / "GALLERY_HEADER.rst", "w"): @@ -240,10 +267,10 @@ def _build_demo( if not execute: cmd.extend(("-D", "plot_gallery=0")) - cmd.extend((str(sphinx_dir), str(build_dir / target.value))) + cmd.extend((str(ctx.repo_root), str(ctx.build_dir / target.value))) sphinx_env = os.environ | { "DEMO_STAGING_DIR": str(stage_dir.resolve()), - "GALLERY_OUTPUT_DIR": str(out_dir.resolve().relative_to(sphinx_dir)), + "GALLERY_OUTPUT_DIR": str(out_dir.resolve().relative_to(ctx.repo_root)), # Make sure demos can find scripts installed in the build venv "PATH": f"{os.environ['PATH']}:{build_venv.path / 'bin'}", } @@ -259,9 +286,9 @@ def _build_demo( if package: _package_demo( demo, - build_dir / "pack", - sphinx_dir / "_static", - build_dir / target.value, + ctx.build_dir / "pack", + ctx.repo_root / "_static", + ctx.build_dir / target.value, out_dir, ) diff --git a/lib/qml/lib/pip_tools.py b/lib/qml/lib/pip_tools.py deleted file mode 100644 index 04e7aaf092..0000000000 --- a/lib/qml/lib/pip_tools.py +++ /dev/null @@ -1,105 +0,0 @@ -import subprocess -from collections import defaultdict -from pathlib import Path -import requirements -import tempfile -import itertools -from collections.abc import Mapping, Sequence - - -class RequirementsGenerator: - """Generates 'requirements.txt' from 'requirements.in' files, with versions constrained - by a global 'constraints.txt' file.""" - - global_constraints: Mapping[str, Sequence[str]] - extra_index_urls: Sequence[str] = ("https://download.pytorch.org/whl/cpu",) - - def __init__( - self, - python_bin: Path, - global_constraints_file: Path, - overrides_file: Path | None = None, - ): - self.python_bin = python_bin - - global_constraints: dict[str, tuple[str, ...]] = defaultdict(tuple) - with open(global_constraints_file, "r") as f: - for req in requirements.parse(f): - global_constraints[req.name] += (req.line,) - - if overrides_file: - overrides: dict[str, tuple[str, ...]] = defaultdict(tuple) - with open(overrides_file, "r") as f: - for req in requirements.parse(f): - overrides[req.name] += (req.line,) - - global_constraints.update(overrides) - - self.global_constraints = global_constraints - self._requirements_in_cache: dict[frozenset[str], str] = {} - - def generate_requirements(self, requirements_in: frozenset[str]): - """Generate a 'requirements.txt' file for the given list of input requirements, with versions from - global constraints. If any of the input requirements are version-qualified, the versions will - override the global constraints. - """ - if cached := self._requirements_in_cache.get(requirements_in): - return cached - - constraints: dict[str, tuple[str, ...]] = defaultdict(tuple) - for req_str in requirements_in: - req = next(requirements.parse(req_str)) - if req.line != req.name: - constraints[req.name] += (req.line,) - - for package in filter( - lambda package: package not in constraints, self.global_constraints - ): - constraints[package] = self.global_constraints[package] - - with tempfile.TemporaryDirectory() as tmpdir: - constraints_file = Path(tmpdir, "constraints.txt") - requirements_file = Path(tmpdir, "requirements.txt") - - with open(constraints_file, "w") as f: - for spec in itertools.chain.from_iterable(constraints.values()): - f.write(spec + "\n") - - with open(requirements_file, "w") as f: - for req_str in requirements_in: - f.write(req_str + "\n") - - cmd = [ - str(self.python_bin), - "-m", - "uv", - "pip", - "compile", - "--index-strategy", - "unsafe-best-match", - "--no-header", - "--no-strip-extras", - "--no-strip-markers", - "--universal", - "--quiet", - "--no-annotate", - "--emit-index-url", - "--constraints", - str(constraints_file), - "--output-file", - str(requirements_file), - ] - for index_url in self.extra_index_urls: - cmd.extend(("--extra-index-url", index_url)) - - cmd.append(str(requirements_file)) - - subprocess.run( - cmd, - ).check_returncode() - - with open(requirements_file, "r") as f: - reqs = f.read() - self._requirements_in_cache[requirements_in] = reqs - - return reqs diff --git a/lib/qml/lib/template.py b/lib/qml/lib/template.py new file mode 100644 index 0000000000..931394ca44 --- /dev/null +++ b/lib/qml/lib/template.py @@ -0,0 +1,61 @@ +from typing import Any +from datetime import timezone, datetime +import inspect + + +def metadata( + title: str, + description: str, + authors: list[str], + thumbnail: str | None, + large_thumbnail: str | None, + categories: list[str], +) -> dict[str, Any]: + today = datetime.now(tz=timezone.utc) + today = today.replace(hour=0, minute=0, second=0, microsecond=0) + + metadata = { + "title": title, + "authors": authors, + "dateOfPublication": today.isoformat(), + "dateOfLastModification": today.isoformat(), + "categores": categories, + "tags": [], + "previewImages": [], + "seoDescription": description, + "doi": "", + "references": [], + "basedOnPapers": [], + "referencedByPapers": [], + "relatedContent": [], + } + + if thumbnail: + metadata["previewImages"].append({"type": "thumbnail", "uri": thumbnail}) + if large_thumbnail: + metadata["previewImages"].append( + {"type": "large_thumbnail", "uri": large_thumbnail} + ) + + return metadata + + +def demo(title: str) -> str: + return inspect.cleandoc(f''' + r""" + {title} + {"="} + + Introduce your demo here! + """ + + print("Hello") + + ############################################################################### + # + # Add comment blocks to separate code blocks + # + + print("World") + + ''') diff --git a/poetry.lock b/poetry.lock index 896d5f8639..ef9867e3c5 100644 --- a/poetry.lock +++ b/poetry.lock @@ -2769,6 +2769,17 @@ enabler = ["pytest-enabler (>=2.2)"] test = ["jaraco.test (>=5.4)", "pytest (>=6,!=8.1.*)", "zipp (>=3.17)"] type = ["pytest-mypy"] +[[package]] +name = "inflection" +version = "0.5.1" +description = "A port of Ruby on Rails inflector to Python" +optional = false +python-versions = ">=3.5" +files = [ + {file = "inflection-0.5.1-py2.py3-none-any.whl", hash = "sha256:f38b2b640938a4f35ade69ac3d053042959b62a0f1076a5bbaa1b9526605a8a2"}, + {file = "inflection-0.5.1.tar.gz", hash = "sha256:1a29730d366e996aaacffb2f1f1cb9593dc38e2ddd30c91250c6dde09ea9b417"}, +] + [[package]] name = "iniconfig" version = "2.0.0" diff --git a/pyproject.toml b/pyproject.toml index e3d9dffd0c..9f68fec8eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,7 @@ dulwich = "<0.22" requirements-parser = "^0.11.0" lxml = "^5.3.0" uv = "^0.5.25" +inflection = "^0.5.1" [tool.poetry.group.base.dependencies] # Base dependencies needed to build the website without any code execution (*-norun)