diff --git a/CHANGELOG.md b/CHANGELOG.md index c75081da..fdbb4bba 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel ### Added +- added MPO zip-up as default long-range gate application method ([#449]) ([**@aaronleesander**]) - parallelized MPO equivalence checker ([#448]) ([**@aaronleesander**]) - added tebd/hybrid/tdvp options for circuit simulation ([#441]) ([**@aaronleesander**]) - added Krylov tolerance as high-level parameter ([#439]) ([**@aaronleesander**]) @@ -126,6 +127,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool +[#449]: https://github.com/munich-quantum-toolkit/yaqs/pull/449 [#448]: https://github.com/munich-quantum-toolkit/yaqs/pull/448 [#447]: https://github.com/munich-quantum-toolkit/yaqs/pull/447 [#446]: https://github.com/munich-quantum-toolkit/yaqs/pull/446 diff --git a/docs/examples/simulation_parameters.md b/docs/examples/simulation_parameters.md index 2ee87561..f3140a66 100644 --- a/docs/examples/simulation_parameters.md +++ b/docs/examples/simulation_parameters.md @@ -175,11 +175,11 @@ Used for noisy strong circuit simulation. Provide observables and optionally ena ### Two-qubit gate mode (`gate_mode`) -Digital circuit simulation on an MPS uses a **TEBD/SVD** two-qubit update by default (`gate_mode="tebd"`). +Digital circuit simulation on an MPS defaults to **`gate_mode="mpo"`** (generic MPO--MPS application): nearest-neighbor gates use the same local TEBD/SVD path as `swaps`, and long-range gates contract an extended gate MPO site-wise (library leg ordering, MPS virtual index before MPO virtual index) followed by compression with `svd_threshold` and `max_bond_dim`. Other modes differ only in how two-qubit gates are applied: -You can also set `gate_mode="hybrid"` to use TEBD/SVD for **nearest-neighbor** gates (adjacent in the internal MPS site order after bit reversal) while keeping the existing **generator MPO + two-site TDVP** path for **long-range** gates. - -Set `gate_mode="tdvp"` to apply TDVP to every two-qubit gate. Set `gate_mode="tebd"` to use TEBD/SVD for all two-qubit gates; long-range gates are implemented by adjacent SWAP insertion before and after the update. +- **`swaps`** — TEBD/SVD for every two-qubit gate; long-range gates are routed with adjacent SWAP insertion before and after the local update. +- **`tdvp`** — TEBD/SVD on nearest-neighbor gates; long-range gates use the generator MPO + two-site TDVP path. +- **`full-tdvp`** — TDVP (generator MPO + two-site TDVP) on every two-qubit gate. ```{code-cell} ipython3 strong = StrongSimParams( diff --git a/src/mqt/yaqs/core/data_structures/mpo.py b/src/mqt/yaqs/core/data_structures/mpo.py index 172621e9..653adc2c 100644 --- a/src/mqt/yaqs/core/data_structures/mpo.py +++ b/src/mqt/yaqs/core/data_structures/mpo.py @@ -10,7 +10,7 @@ from __future__ import annotations import re -from typing import ClassVar +from typing import TYPE_CHECKING, ClassVar, cast, overload import numpy as np import opt_einsum as oe @@ -18,9 +18,19 @@ from numpy.typing import NDArray from .. import linalg -from ..libraries.gate_library import Destroy +from ..libraries.gate_library import BaseGate, Destroy +from .mpo_utils import ( + contract_mpo_site_with_mpo_site, + contract_mpo_site_with_mps_site, + gate_support_mpo_tensors, + identity_mpo_site, +) from .mps import MPS +if TYPE_CHECKING: + from ..methods.decompositions import TruncMode + from .simulation_parameters import StrongSimParams, WeakSimParams + ComplexTensor = NDArray[np.complex128] @@ -49,6 +59,8 @@ class MPO: **Operations** + - ``from_gate(...)``: build an MPO from a two-qubit gate on a chain (optionally identity-padded). + - ``multiply(MPS)`` / ``multiply(MPO)``: apply this MPO to an MPS or left-multiply into another MPO. - ``compress(...)``: SVD-based bond compression sweeps. - ``rotate(...)``: swap physical legs (optionally conjugating). @@ -597,6 +609,58 @@ def identity(cls, length: int, physical_dimension: int = 2) -> MPO: mpo.init_identity(length, physical_dimension=physical_dimension) return mpo + @classmethod + def from_gate(cls, gate: BaseGate, chain_length: int) -> MPO: + """Build an MPO for a two-qubit gate on a chain. + + When ``chain_length`` equals the gate support size, the MPO contains only the + extended gate tensors. When ``chain_length`` is larger, identity MPO sites are + placed outside the support interval ``[min(sites), max(sites)]``. + + Reuses :attr:`~mqt.yaqs.core.libraries.gate_library.BaseGate.mpo_tensors` when + already populated for the gate support. + + Args: + gate: Two-qubit gate with ``sites`` and ``tensor`` (or ``mpo_tensors``) set. + chain_length: Total number of MPO sites (support length or full MPS length). + + Returns: + MPO ready for :meth:`multiply` on an MPS or another MPO. + + Raises: + ValueError: If the gate is not two-qubit or ``chain_length`` is too small. + """ + if gate.interaction != 2: + msg = f"from_gate requires a two-qubit gate, got interaction {gate.interaction}." + raise ValueError(msg) + if len(gate.sites) != 2: + msg = f"from_gate requires exactly two sites, got {len(gate.sites)}." + raise ValueError(msg) + + first_site = min(gate.sites[0], gate.sites[1]) + last_site = max(gate.sites[0], gate.sites[1]) + support_len = last_site - first_site + 1 + if chain_length < support_len: + msg = f"chain_length {chain_length} is smaller than gate support length {support_len}." + raise ValueError(msg) + + support = gate_support_mpo_tensors(gate, first_site=first_site, last_site=last_site) + if chain_length == support_len: + tensors = support + else: + phys_dim = support[0].shape[0] + identity_site = identity_mpo_site(phys_dim) + tensors = [] + for site in range(chain_length): + if site < first_site or site > last_site: + tensors.append(np.array(identity_site, copy=True)) + else: + tensors.append(support[site - first_site]) + + mpo = cls() + mpo.custom(tensors, transpose=False) + return mpo + def init_identity(self, length: int, physical_dimension: int = 2) -> None: """Initialize this MPO in place as the identity operator. @@ -606,14 +670,13 @@ def init_identity(self, length: int, physical_dimension: int = 2) -> None: length: Number of sites. physical_dimension: Local dimension per site (default 2). """ - mat = np.eye(physical_dimension, dtype=np.complex128) - mat = np.expand_dims(mat, (2, 3)) + site = identity_mpo_site(physical_dimension) self.length = length self.physical_dimension = physical_dimension self.tensors = [] for _ in range(length): - self.tensors.append(mat.copy()) + self.tensors.append(np.array(site, copy=True)) def finite_state_machine( self, @@ -920,6 +983,167 @@ def _compress_one_sweep(self, *, direction: str, tol: float, max_bond_dim: int | self.tensors[k] = left self.tensors[k + 1] = right + @overload + def multiply( + self, + other: MPS, + *, + sim_params: StrongSimParams | WeakSimParams | None = None, + compress: bool = True, + ) -> None: ... + + @overload + def multiply( + self, + other: MPO, + *, + start_site: int = 0, + conjugate: bool = False, + compress: bool = True, + tol: float = 1e-12, + max_bond_dim: int | None = None, + n_sweeps: int = 1, + directions: str = "lr_rl", + ) -> None: ... + + def multiply( + self, + other: MPS | MPO, + *, + sim_params: StrongSimParams | WeakSimParams | None = None, + compress: bool = True, + start_site: int = 0, + conjugate: bool = False, + tol: float = 1e-12, + max_bond_dim: int | None = None, + n_sweeps: int = 1, + directions: str = "lr_rl", + ) -> None: + """Left-multiply this MPO into ``other`` (MPS or MPO), updating ``other`` in place. + + For an :class:`~mqt.yaqs.core.data_structures.mps.MPS`, each site is updated by + :func:`~mqt.yaqs.core.data_structures.mpo_utils.contract_mpo_site_with_mps_site`, + optionally followed by a two-site SVD compression sweep driven by ``sim_params``. + + For another :class:`MPO`, each site uses the equivalence-checking + :func:`~mqt.yaqs.core.data_structures.mpo_utils.contract_mpo_site_with_mpo_site` + contraction (``abcd,cefg``), optionally + followed by :meth:`compress` on ``other``. + + Args: + other: Target MPS or MPO to update in place. + sim_params: Truncation settings for MPS compression (required if ``compress`` + is True and ``other`` is an MPS). + compress: Whether to run a compression sweep after contraction. + start_site: When ``len(self) != len(other)``, index on ``other`` where this + MPO is embedded (only for MPO targets). + conjugate: Use the conjugated MPO--MPO contraction (MPO targets only). + tol: MPO compression threshold (MPO targets only). + max_bond_dim: Optional bond-dimension cap for MPO compression. + n_sweeps: Number of MPO compression sweeps. + directions: MPO compression sweep schedule (see :meth:`compress`). + + Raises: + TypeError: If ``other`` is neither an MPS nor an MPO. + """ + if isinstance(other, MPS): + self._multiply_mps( + other, + sim_params=sim_params, + compress=compress, + ) + return + + if not isinstance(other, MPO): + msg = f"multiply expects MPS or MPO, got {type(other).__name__}." + raise TypeError(msg) + + self._multiply_mpo( + other, + start_site=start_site, + conjugate=conjugate, + compress=compress, + tol=tol, + max_bond_dim=max_bond_dim, + n_sweeps=n_sweeps, + directions=directions, + ) + + def _multiply_mps( + self, + state: MPS, + *, + sim_params: StrongSimParams | WeakSimParams | None, + compress: bool, + ) -> None: + """Apply this MPO to ``state`` with optional compression. + + Raises: + ValueError: On length mismatch or missing ``sim_params`` when compressing. + """ + if len(self.tensors) != state.length: + msg = f"MPO length {len(self.tensors)} does not match MPS length {state.length}." + raise ValueError(msg) + + for site, operator in enumerate(self.tensors): + state.tensors[site] = contract_mpo_site_with_mps_site(operator, state.tensors[site]) + + if not compress: + return + if sim_params is None: + msg = "sim_params is required when compress=True for MPO.multiply(MPS)." + raise ValueError(msg) + + state.compress( + sim_params.svd_threshold, + max_bond_dim=sim_params.max_bond_dim, + min_bond_dim=sim_params.min_bond_dim, + trunc_mode=cast("TruncMode", sim_params.trunc_mode), + ) + + def _multiply_mpo( + self, + other: MPO, + *, + start_site: int, + conjugate: bool, + compress: bool, + tol: float, + max_bond_dim: int | None, + n_sweeps: int, + directions: str, + ) -> None: + """Left-multiply this MPO into ``other``. + + Raises: + ValueError: If this MPO cannot be embedded at ``start_site``. + """ + gate_len = len(self.tensors) + target_len = len(other.tensors) + + if gate_len == target_len: + sites = range(target_len) + elif start_site >= 0 and start_site + gate_len <= target_len: + sites = range(start_site, start_site + gate_len) + else: + msg = f"Cannot embed MPO of length {gate_len} at start_site={start_site} into MPO of length {target_len}." + raise ValueError(msg) + + for gate_site, target_site in enumerate(sites): + other.tensors[target_site] = contract_mpo_site_with_mpo_site( + self.tensors[gate_site], + other.tensors[target_site], + conjugate=conjugate, + ) + + if compress: + other.compress( + tol=tol, + max_bond_dim=max_bond_dim, + n_sweeps=n_sweeps, + directions=directions, + ) + def rotate(self, *, conjugate: bool = False) -> None: """Rotates MPO. diff --git a/src/mqt/yaqs/core/data_structures/mpo_utils.py b/src/mqt/yaqs/core/data_structures/mpo_utils.py new file mode 100644 index 00000000..af300bf3 --- /dev/null +++ b/src/mqt/yaqs/core/data_structures/mpo_utils.py @@ -0,0 +1,191 @@ +# Copyright (c) 2025 - 2026 Chair for Design Automation, TUM +# All rights reserved. +# +# SPDX-License-Identifier: MIT +# +# Licensed under the MIT License + +"""Low-level MPO/MPS tensor contractions and gate-MPO construction helpers.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +import opt_einsum as oe +from numpy.typing import NDArray + +from .. import linalg +from ..libraries.gate_library import extend_gate + +if TYPE_CHECKING: + from ..libraries.gate_library import BaseGate + +ComplexTensor = NDArray[np.complex128] + + +def contract_mpo_site_with_mps_site( + mpo_tensor: NDArray[np.complex128], + mps_tensor: NDArray[np.complex128], +) -> NDArray[np.complex128]: + """Contract one MPO site with one MPS site using library leg ordering. + + MPO layout: ``(phys_out, phys_in, chi_left, chi_right)``. + MPS layout: ``(phys, chi_left, chi_right)``. + + Virtual bonds are fused with MPS indices first and MPO indices second on both sides. + + Args: + mpo_tensor: Gate or operator MPO tensor at one site. + mps_tensor: MPS site tensor. + + Returns: + Updated MPS site tensor ``(phys_out, chi_left * mpo_left, chi_right * mpo_right)``. + """ + operator = np.asarray(mpo_tensor, dtype=np.complex128) + site = np.asarray(mps_tensor, dtype=np.complex128) + theta = np.tensordot(operator, site, axes=([1], [0])) + phys_out, mpo_left, mpo_right, mps_left, mps_right = theta.shape + return np.asarray( + theta.transpose(0, 3, 1, 4, 2).reshape( + phys_out, + mps_left * mpo_left, + mps_right * mpo_right, + ), + dtype=np.complex128, + ) + + +def contract_mpo_site_with_mpo_site( + left_mpo_tensor: NDArray[np.complex128], + right_mpo_tensor: NDArray[np.complex128], + *, + conjugate: bool = False, +) -> NDArray[np.complex128]: + """Contract two MPO site tensors (left factor times right factor). + + Uses the transposed layout ``(phys_out, chi_left, phys_in, chi_right)`` and the + ``abcd,cefg`` contraction from equivalence-checking MPO updates. + + Args: + left_mpo_tensor: Left MPO factor at one site, library order + ``(phys_out, phys_in, chi_left, chi_right)``. + right_mpo_tensor: Right MPO factor at the same site. + conjugate: If True, use the conjugated contraction used when updating + the right-hand MPO in equivalence checking. + + Returns: + Product tensor in library order ``(phys_out, phys_in, chi_left, chi_right)``. + """ + tensor1 = np.transpose(np.asarray(left_mpo_tensor, dtype=np.complex128), (0, 2, 1, 3)) + tensor2 = np.transpose(np.asarray(right_mpo_tensor, dtype=np.complex128), (0, 2, 1, 3)) + if conjugate: + theta = oe.contract("abcd,cefg->febagd", tensor1, tensor2) + else: + theta = oe.contract("abcd,cefg->abefdg", tensor1, tensor2) + dims = theta.shape + fused = np.reshape(theta, (dims[0], dims[1] * dims[2], dims[3], dims[4] * dims[5])) + return np.transpose(fused, (0, 2, 1, 3)) + + +def identity_mpo_site(physical_dimension: int) -> NDArray[np.complex128]: + """Single-site identity MPO tensor ``(phys_out, phys_in, 1, 1)``. + + Args: + physical_dimension: Local Hilbert-space dimension per site. + + Returns: + Identity MPO site tensor. + """ + tensor = np.eye(physical_dimension, dtype=np.complex128) + return np.expand_dims(np.expand_dims(tensor, axis=2), axis=3) + + +def gate_tensor_lr_order( + gate: BaseGate, + left_site: int | None = None, + right_site: int | None = None, +) -> NDArray[np.complex128]: + """Return ``gate.tensor`` as ``U[out_l, out_r, in_l, in_r]`` on ascending MPS sites. + + Args: + gate: Two-qubit gate with ``sites`` and ``tensor`` set. + left_site: Lower MPS site index; defaults to ``min(gate.sites)``. + right_site: Higher MPS site index; defaults to ``max(gate.sites)``. + + Returns: + Gate tensor with left-then-right site axis order. + + Raises: + ValueError: If ``gate.sites`` does not match ``(left_site, right_site)``. + """ + if left_site is None or right_site is None: + site0, site1 = gate.sites[0], gate.sites[1] + left_site = min(site0, site1) + right_site = max(site0, site1) + if gate.sites[0] == left_site and gate.sites[1] == right_site: + return np.asarray(gate.tensor, dtype=np.complex128) + if gate.sites[0] == right_site and gate.sites[1] == left_site: + return np.asarray(np.transpose(gate.tensor, (1, 0, 3, 2)), dtype=np.complex128) + msg = f"Gate sites {gate.sites!r} are not consistent with MPS sites ({left_site}, {right_site})." + raise ValueError(msg) + + +def gate_support_mpo_tensors( + gate: BaseGate, + *, + first_site: int, + last_site: int, +) -> list[NDArray[np.complex128]]: + """MPO tensors for the gate support ``[first_site, last_site]`` in library order. + + Args: + gate: Two-qubit gate with optional cached ``mpo_tensors``. + first_site: First site of the support interval (inclusive). + last_site: Last site of the support interval (inclusive). + + Returns: + Support MPO tensors from the gate cache or :func:`~mqt.yaqs.core.libraries.gate_library.extend_gate`. + """ + support_len = last_site - first_site + 1 + try: + cached = gate.mpo_tensors + except AttributeError: + cached = None + if cached is not None and len(cached) == support_len: + return list(cached) + return extend_gate( + gate_tensor_lr_order(gate), + [first_site, last_site], + ) + + +def decompose_theta( + theta: NDArray[np.complex128], threshold: float +) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]: + """SVD-split a fused two-site MPO tensor back into two rank-4 site tensors. + + Args: + theta: High-rank tensor from contracting two MPO sites (and optional gates). + threshold: Singular-value cutoff for truncation. + + Returns: + Left and right MPO site tensors after truncated SVD. + """ + dims = theta.shape + theta = np.transpose(theta, (0, 3, 2, 1, 4, 5)) + theta_matrix = np.reshape(theta, (dims[0] * dims[1] * dims[2], dims[3] * dims[4] * dims[5])) + + u_mat, s_list, v_mat = linalg.svd(theta_matrix, full_matrices=False) + keep = linalg.truncate(s_list, mode="hard_cutoff", threshold=threshold, min_keep=1) + s_list = s_list[:keep] + u_mat = u_mat[:, :keep] + v_mat = v_mat[:keep, :] + + u_tensor = np.reshape(u_mat, (dims[0], dims[1], dims[2], len(s_list))) + + m_mat = np.diag(s_list) @ v_mat + m_tensor = np.reshape(m_mat, (len(s_list), dims[3], dims[4], dims[5])) + m_tensor = np.transpose(m_tensor, (1, 2, 0, 3)) + + return u_tensor, m_tensor diff --git a/src/mqt/yaqs/core/data_structures/mps.py b/src/mqt/yaqs/core/data_structures/mps.py index d72fa609..4d67dc80 100644 --- a/src/mqt/yaqs/core/data_structures/mps.py +++ b/src/mqt/yaqs/core/data_structures/mps.py @@ -24,6 +24,7 @@ if TYPE_CHECKING: from numpy.typing import NDArray + from ..methods.decompositions import TruncMode from .simulation_parameters import AnalogSimParams, Observable, StrongSimParams @@ -597,44 +598,61 @@ def normalize(self, form: str = "B", decomposition: str = "QR") -> None: if form == "B": self.flip_network() - def truncate(self, threshold: float = 1e-12, max_bond_dim: int | None = None) -> None: - """In-place MPS truncation via repeated two-site SVDs.""" - orth_center = self.check_canonical_form()[0] + def compress( + self, + threshold: float, + *, + max_bond_dim: int | None = None, + min_bond_dim: int = 2, + trunc_mode: TruncMode = "discarded_weight", + ) -> None: + """Compress in place via left-to-center and right-to-left two-site SVD sweeps. + + Args: + threshold: SVD truncation threshold (e.g. ``sim_params.svd_threshold``). + max_bond_dim: Optional cap on bond dimension. + min_bond_dim: Minimum bond dimension to retain when truncating. + trunc_mode: ``"discarded_weight"`` or ``"relative"``. + """ if self.length == 1: return - # ——— left­-to-­center sweep ——— - for i in range(orth_center): - a, b = self.tensors[i], self.tensors[i + 1] - merged = merge_two_site(a, b) - a_new, b_new = split_two_site( + canonical = self.check_canonical_form() + orth_center = canonical[0] if canonical and canonical[0] >= 0 else self.length - 1 + + for site in range(orth_center): + left_tensor = self.tensors[site] + right_tensor = self.tensors[site + 1] + merged = merge_two_site(left_tensor, right_tensor) + left_new, right_new = split_two_site( merged, - [a.shape[0], b.shape[0]], + [left_tensor.shape[0], right_tensor.shape[0]], svd_distribution="right", - trunc_mode="discarded_weight", + trunc_mode=trunc_mode, threshold=threshold, max_bond_dim=max_bond_dim, - min_bond_dim=2, + min_bond_dim=min_bond_dim, ) - self.tensors[i], self.tensors[i + 1] = a_new, b_new + self.tensors[site] = left_new + self.tensors[site + 1] = right_new - # flip the network and sweep back self.flip_network() orth_flipped = self.length - 1 - orth_center - for i in range(orth_flipped): - a, b = self.tensors[i], self.tensors[i + 1] - merged = merge_two_site(a, b) - a_new, b_new = split_two_site( + for site in range(orth_flipped): + left_tensor = self.tensors[site] + right_tensor = self.tensors[site + 1] + merged = merge_two_site(left_tensor, right_tensor) + left_new, right_new = split_two_site( merged, - [a.shape[0], b.shape[0]], + [left_tensor.shape[0], right_tensor.shape[0]], svd_distribution="right", - trunc_mode="discarded_weight", + trunc_mode=trunc_mode, threshold=threshold, max_bond_dim=max_bond_dim, - min_bond_dim=2, + min_bond_dim=min_bond_dim, ) - self.tensors[i], self.tensors[i + 1] = a_new, b_new - + self.tensors[site] = left_new + self.tensors[site + 1] = right_new self.flip_network() def scalar_product(self, other: MPS, sites: int | list[int] | None = None) -> np.complex128: diff --git a/src/mqt/yaqs/core/data_structures/simulation_parameters.py b/src/mqt/yaqs/core/data_structures/simulation_parameters.py index 4e275c89..320b1b9f 100644 --- a/src/mqt/yaqs/core/data_structures/simulation_parameters.py +++ b/src/mqt/yaqs/core/data_structures/simulation_parameters.py @@ -29,7 +29,7 @@ from mqt.yaqs.core.libraries.gate_library import BaseGate SimulationPreset = Literal["fast", "balanced", "accurate", "exact"] -GateMode = Literal["hybrid", "tdvp", "tebd"] +GateMode = Literal["tdvp", "full-tdvp", "swaps", "mpo"] class PresetTypes(TypedDict): @@ -102,7 +102,7 @@ def _validate_gate_mode(mode: GateMode) -> GateMode: Raises: ValueError: If ``mode`` is not a supported value. """ - allowed = ("hybrid", "tdvp", "tebd") + allowed = ("tdvp", "full-tdvp", "swaps", "mpo") if mode not in allowed: msg = f"gate_mode must be one of {allowed!r}, got {mode!r}." raise ValueError(msg) @@ -428,7 +428,7 @@ class StrongSimParams(_ObservableOrderingMixin): get_state: If ``True``, request the final state on the returned :class:`~mqt.yaqs.Result`. gate_mode: Two-qubit gate update mode on the MPS digital backend - (``"tebd"``, ``"hybrid"``, or ``"tdvp"``). + (``"swaps"``, ``"tdvp"``, ``"full-tdvp"``, or ``"mpo"``). """ # Properties set as placeholders for code compatibility @@ -449,7 +449,7 @@ def __init__( sample_layers: bool = False, num_mid_measurements: int = 0, random_seed: int | None = None, - gate_mode: GateMode = "tebd", + gate_mode: GateMode = "mpo", ) -> None: """Strong circuit simulation parameters initialization. @@ -475,7 +475,7 @@ def __init__( sample_layers: If ``True``, record observables at sampled circuit layers. num_mid_measurements: Number of mid-circuit measurement barriers when sampling layers. random_seed: If set, makes stochastic trajectories and noise-model sampling reproducible. - gate_mode: Two-qubit gate update mode (default ``"tebd"``). + gate_mode: Two-qubit gate update mode (default ``"mpo"``). """ _validate_random_seed(random_seed) preset_values = SIMULATION_PRESETS[_validate_preset(preset)] @@ -527,7 +527,7 @@ class WeakSimParams: get_state: If ``True``, request the final state on the returned :class:`~mqt.yaqs.Result`. gate_mode: Two-qubit gate update mode on the MPS digital backend - (``"tebd"``, ``"hybrid"``, or ``"tdvp"``). + (``"swaps"``, ``"tdvp"``, ``"full-tdvp"``, or ``"mpo"``). """ # Properties set as placeholders for code compatibility @@ -546,7 +546,7 @@ def __init__( preset: SimulationPreset = "balanced", get_state: bool = False, random_seed: int | None = None, - gate_mode: GateMode = "tebd", + gate_mode: GateMode = "mpo", ) -> None: """Weak circuit simulation initialization. @@ -569,7 +569,7 @@ def __init__( svd_threshold: SVD truncation threshold for bond dimension control. get_state: If ``True``, request the final state on the returned :class:`~mqt.yaqs.Result`. random_seed: If set, makes per-shot jump RNG reproducible. - gate_mode: Two-qubit gate update mode (default ``"tebd"``). + gate_mode: Two-qubit gate update mode (default ``"mpo"``). """ _validate_random_seed(random_seed) preset_values = SIMULATION_PRESETS[_validate_preset(preset)] diff --git a/src/mqt/yaqs/core/methods/bug.py b/src/mqt/yaqs/core/methods/bug.py index 9e95724c..31560c7e 100644 --- a/src/mqt/yaqs/core/methods/bug.py +++ b/src/mqt/yaqs/core/methods/bug.py @@ -220,4 +220,4 @@ def bug(state: MPS, mpo: MPO, sim_params: AnalogSimParams | WeakSimParams | Stro ) state.tensors[0] = updated_tensor # Truncation - state.truncate(sim_params.svd_threshold, sim_params.max_bond_dim) + state.compress(sim_params.svd_threshold, max_bond_dim=sim_params.max_bond_dim) diff --git a/src/mqt/yaqs/digital/digital_tjm.py b/src/mqt/yaqs/digital/digital_tjm.py index b828c073..28ea057c 100644 --- a/src/mqt/yaqs/digital/digital_tjm.py +++ b/src/mqt/yaqs/digital/digital_tjm.py @@ -22,13 +22,14 @@ from qiskit.converters import circuit_to_dag from ..core.data_structures.mpo import MPO +from ..core.data_structures.mpo_utils import gate_tensor_lr_order from ..core.data_structures.mps import MPS from ..core.data_structures.noise_model import NoiseModel from ..core.data_structures.simulation_parameters import ( StrongSimParams, WeakSimParams, ) -from ..core.libraries.gate_library import GateLibrary +from ..core.libraries.gate_library import BaseGate, GateLibrary from ..core.methods.decompositions import merge_two_site, split_two_site from ..core.methods.dissipation import apply_dissipation from ..core.methods.stochastic_process import stochastic_process @@ -223,32 +224,6 @@ def apply_window(state: MPS, mpo: MPO, first_site: int, last_site: int, window_s return short_state, short_mpo, window -def _gate_tensor_left_right_order( - gate: BaseGate, - left_site: int, - right_site: int, -) -> NDArray[np.complex128]: - """Return ``gate.tensor`` with axes ordered as (left site, right site). - - Args: - gate: Two-qubit gate with ``sites`` and ``tensor`` already set. - left_site: Lower MPS site index. - right_site: Higher MPS site index. - - Returns: - Gate tensor ``U[out_l, out_r, in_l, in_r]`` for the left/right ordering. - - Raises: - ValueError: If ``gate.sites`` does not match ``(left_site, right_site)``. - """ - if gate.sites[0] == left_site and gate.sites[1] == right_site: - return gate.tensor - if gate.sites[0] == right_site and gate.sites[1] == left_site: - return np.transpose(gate.tensor, (1, 0, 3, 2)) - msg = f"Gate sites {gate.sites!r} are not consistent with MPS sites ({left_site}, {right_site})." - raise ValueError(msg) - - def apply_two_qubit_gate_tdvp( state: MPS, gate: BaseGate, @@ -318,7 +293,7 @@ def apply_swap(site_left: int) -> None: left_site = min(site0, site1) right_site = max(site0, site1) - u_gate = _gate_tensor_left_right_order(gate, left_site, right_site) + u_gate = gate_tensor_lr_order(gate, left_site, right_site) left_tensor = state.tensors[left_site] right_tensor = state.tensors[right_site] @@ -343,6 +318,28 @@ def apply_swap(site_left: int) -> None: return left_site, right_site +def apply_long_range_gate_mpo( + state: MPS, + gate: BaseGate, + sim_params: StrongSimParams | WeakSimParams, +) -> tuple[int, int]: + """Apply a long-range two-qubit gate via :meth:`~mqt.yaqs.core.data_structures.mpo.MPO.multiply`. + + Args: + state: MPS updated in place. + gate: Two-qubit gate with sites and MPO data populated. + sim_params: Truncation settings for the compression sweep. + + Returns: + ``(first_site, last_site)`` spanning the gate support in MPS order. + """ + site0, site1 = gate.sites[0], gate.sites[1] + first_site = min(site0, site1) + last_site = max(site0, site1) + MPO.from_gate(gate, state.length).multiply(state, sim_params=sim_params, compress=True) + return first_site, last_site + + def apply_two_qubit_gate(state: MPS, node: DAGOpNode, sim_params: StrongSimParams | WeakSimParams) -> tuple[int, int]: """Apply a two-qubit gate using the configured two-qubit gate mode. @@ -360,19 +357,24 @@ def apply_two_qubit_gate(state: MPS, node: DAGOpNode, sim_params: StrongSimParam gate = convert_dag_to_tensor_algorithm(node)[0] site0, site1 = gate.sites[0], gate.sites[1] is_nearest_neighbor = abs(site0 - site1) == 1 - gate_mode: GateMode = getattr(sim_params, "gate_mode", "hybrid") + gate_mode: GateMode = getattr(sim_params, "gate_mode", "mpo") - if gate_mode == "tdvp": + if gate_mode == "full-tdvp": return apply_two_qubit_gate_tdvp(state, gate, sim_params) - if gate_mode == "tebd": + if gate_mode == "swaps": return apply_two_qubit_gate_tebd(state, gate, sim_params) - if gate_mode == "hybrid": + if gate_mode == "tdvp": if is_nearest_neighbor: return apply_two_qubit_gate_tebd(state, gate, sim_params) return apply_two_qubit_gate_tdvp(state, gate, sim_params) + if gate_mode == "mpo": + if is_nearest_neighbor: + return apply_two_qubit_gate_tebd(state, gate, sim_params) + return apply_long_range_gate_mpo(state, gate, sim_params) + msg = f"Unknown gate_mode: {gate_mode!r}" raise ValueError(msg) diff --git a/src/mqt/yaqs/digital/utils/mpo_utils.py b/src/mqt/yaqs/digital/utils/contraction_utils.py similarity index 86% rename from src/mqt/yaqs/digital/utils/mpo_utils.py rename to src/mqt/yaqs/digital/utils/contraction_utils.py index 19561ed9..69586524 100644 --- a/src/mqt/yaqs/digital/utils/mpo_utils.py +++ b/src/mqt/yaqs/digital/utils/contraction_utils.py @@ -5,13 +5,10 @@ # # Licensed under the MIT License -"""Utility functions for MPO-based equivalence checking. +"""DAG-driven MPO update routines for equivalence checking. -This module provides functions to manipulate tensor network representations of quantum operations. -It includes routines for performing SVD-based decompositions (decompose_theta), applying gates to local tensors -(apply_gate), extracting and applying temporal zones from DAGCircuits (apply_temporal_zone), and updating Matrix -Product Operators (MPO) via layer-wise and long-range gate applications. These utilities facilitate the -conversion of quantum circuits into tensor network algorithms for checking equivalence. +Applies temporal zones from paired circuits to an MPO via local contractions, gate updates, +and SVD splits (see :mod:`mqt.yaqs.core.data_structures.mpo_utils` for primitive tensor ops). """ from __future__ import annotations @@ -24,8 +21,8 @@ import opt_einsum as oe from qiskit.converters import dag_to_circuit -from ...core import linalg from ...core.data_structures.mpo import MPO +from ...core.data_structures.mpo_utils import contract_mpo_site_with_mpo_site, decompose_theta from ...parallel_utils import MPContext, available_cpus, limit_worker_threads from .dag_utils import check_longest_gate, convert_dag_to_tensor_algorithm, get_temporal_zone, select_starting_point @@ -40,51 +37,6 @@ from ...core.libraries.gate_library import BaseGate -def decompose_theta( - theta: NDArray[np.complex128], threshold: float -) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]: - """Decompose theta during equivalence checking. - - Perform an SVD-based decomposition of the tensor `theta`, truncating singular values below a given threshold, - and reshape the result into two rank-4 tensors. - - The input tensor is first re-ordered and flattened into a matrix, then decomposed using SVD. - Singular values below the threshold are discarded. The left singular vectors (U) are reshaped into a tensor - with shape (dims[0], dims[1], dims[2], num_sv), and the product of the diagonal singular value matrix with - the right singular vectors (V) is reshaped and transposed into a tensor with shape - (dims[3], dims[4], num_sv, dims[5]). - - Args: - theta (NDArray[np.complex128]): The high-rank tensor to decompose. - threshold (float): The cutoff threshold for singular values. - - Returns: - tuple[NDArray[np.complex128], NDArray[np.complex128]]: - A tuple (U, M) of two reshaped tensors derived from the SVD. - """ - dims = theta.shape - # Reorder indices before flattening - theta = np.transpose(theta, (0, 3, 2, 1, 4, 5)) - theta_matrix = np.reshape(theta, (dims[0] * dims[1] * dims[2], dims[3] * dims[4] * dims[5])) - - u_mat, s_list, v_mat = linalg.svd(theta_matrix, full_matrices=False) - keep = linalg.truncate(s_list, mode="hard_cutoff", threshold=threshold, min_keep=1) - s_list = s_list[:keep] - u_mat = u_mat[:, :keep] - v_mat = v_mat[:keep, :] - - # Reshape U into a tensor of shape (dims[0], dims[1], dims[2], num_sv) - u_tensor = np.reshape(u_mat, (dims[0], dims[1], dims[2], len(s_list))) - - # Compute site tensor M and reshape: first form M = diag(S_list) @ V, - # then reshape to (num_sv, dims[3], dims[4], dims[5]) and transpose to (dims[3], dims[4], num_sv, dims[5]) - m_mat = np.diag(s_list) @ v_mat - m_tensor = np.reshape(m_mat, (len(s_list), dims[3], dims[4], dims[5])) - m_tensor = np.transpose(m_tensor, (1, 2, 0, 3)) - - return u_tensor, m_tensor - - def apply_gate( gate: BaseGate, theta: NDArray[np.complex128], @@ -485,9 +437,7 @@ def apply_long_range_layer(mpo: MPO, dag1: DAGCircuit, dag2: DAGCircuit, thresho and node.qargs[1]._index == gate.qubits[1]._index # noqa: SLF001 ): gate_ = convert_dag_to_tensor_algorithm(node)[0] - mpo_tensors = gate_.mpo_tensors - gate_mpo = MPO() - gate_mpo.custom(mpo_tensors, transpose=False) + gate_mpo = MPO.from_gate(gate_, distance) if conjugate: gate_mpo.rotate(conjugate=True) dag.remove_op_node(node) @@ -533,19 +483,16 @@ def apply_long_range_layer(mpo: MPO, dag1: DAGCircuit, dag2: DAGCircuit, thresho # Process odd-indexed (or hanging) tensor if present. if site_gate_mpo == len(sites) - 1 and any(isinstance(tensor, np.ndarray) for tensor in gate_mpo.tensors): - if not conjugate: - tensor1 = np.transpose(gate_mpo.tensors[site_gate_mpo], (0, 2, 1, 3)) - tensor2 = np.transpose(mpo.tensors[overall_site], (0, 2, 1, 3)) - theta = oe.contract("abcd,cefg->abefdg", tensor1, tensor2) - else: + if conjugate: mpo.rotate() - tensor1 = np.transpose(gate_mpo.tensors[site_gate_mpo], (0, 2, 1, 3)) - tensor2 = np.transpose(mpo.tensors[overall_site], (0, 2, 1, 3)) - theta = oe.contract("abcd,cefg->febagd", tensor1, tensor2) + theta = contract_mpo_site_with_mpo_site( + gate_mpo.tensors[site_gate_mpo], + mpo.tensors[overall_site], + conjugate=conjugate, + ) + if conjugate: mpo.rotate() - - dims = theta.shape - theta = np.reshape(theta, (dims[0], dims[1] * dims[2], dims[3], dims[4] * dims[5])) + theta = np.transpose(theta, (0, 2, 1, 3)) tensor1 = np.transpose(mpo.tensors[overall_site - 1], (0, 2, 1, 3)) theta = oe.contract("abcd, edfg->aebcfg", tensor1, theta) diff --git a/src/mqt/yaqs/equivalence_checker.py b/src/mqt/yaqs/equivalence_checker.py index e3a173a3..b2db3712 100644 --- a/src/mqt/yaqs/equivalence_checker.py +++ b/src/mqt/yaqs/equivalence_checker.py @@ -22,8 +22,8 @@ from qiskit.converters import circuit_to_dag from .core.data_structures.mpo import MPO +from .digital.utils.contraction_utils import iterate from .digital.utils.matrix_utils import check_equivalence_matrix -from .digital.utils.mpo_utils import iterate if TYPE_CHECKING: from qiskit.circuit import QuantumCircuit diff --git a/tests/core/data_structures/test_mpo.py b/tests/core/data_structures/test_mpo.py index c487813d..7c1289ff 100644 --- a/tests/core/data_structures/test_mpo.py +++ b/tests/core/data_structures/test_mpo.py @@ -16,7 +16,8 @@ from mqt.yaqs.core.data_structures.mpo import MPO from mqt.yaqs.core.data_structures.mps import MPS -from mqt.yaqs.core.libraries.gate_library import Destroy, Id +from mqt.yaqs.core.data_structures.simulation_parameters import Observable, StrongSimParams +from mqt.yaqs.core.libraries.gate_library import Destroy, GateLibrary, Id, Z # ---- single-qubit ops ---- _I2 = np.eye(2, dtype=complex) @@ -645,6 +646,101 @@ def test_identity_mpo_tensors_are_independent() -> None: assert mpo.tensors[0] is not mpo.tensors[1] +def test_multiply_mps_identity_preserves_state() -> None: + """Identity MPO.multiply(MPS) leaves the dense state unchanged.""" + length = 4 + state = MPS(length, state="ones") + state.normalize() + expected = np.asarray(state.to_vec(), dtype=np.complex128) + + identity = MPO.identity(length) + identity.multiply(state, compress=False) + + np.testing.assert_allclose(state.to_vec(), expected, atol=1e-10) + + +def test_multiply_mpo_matches_dense_operator_product() -> None: + """Site-wise MPO.multiply(MPO) matches the dense matrix product.""" + length = 4 + left = MPO.ising(length, 1.0, 0.5) + right = MPO.ising(length, 1.0, 0.3) + reference = left.to_matrix() @ right.to_matrix() + + left.multiply(right, compress=False) + np.testing.assert_allclose(right.to_matrix(), reference, atol=1e-10) + + +def test_multiply_mps_with_compression() -> None: + """``multiply(MPS, compress=True)`` runs an SVD sweep using ``sim_params``.""" + length = 3 + state = MPS(length, state="ones") + state.normalize() + gate = GateLibrary.cx() + gate.set_sites(0, 1) + gate_mpo = MPO.from_gate(gate, length) + sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact") + gate_mpo.multiply(state, sim_params=sim_params, compress=True) + state.check_if_valid_mps() + + +def test_multiply_mps_compress_requires_sim_params() -> None: + """Compression without ``sim_params`` raises ``ValueError``.""" + state = MPS(2, state="zeros") + gate = GateLibrary.cx() + gate.set_sites(0, 1) + gate_mpo = MPO.from_gate(gate, 2) + with pytest.raises(ValueError, match="sim_params is required"): + gate_mpo.multiply(state, compress=True) + + +def test_multiply_mps_length_mismatch_raises() -> None: + """MPO and MPS length mismatch raises ``ValueError``.""" + state = MPS(2, state="zeros") + gate = GateLibrary.cx() + gate.set_sites(0, 1) + gate_mpo = MPO.from_gate(gate, 3) + with pytest.raises(ValueError, match="does not match MPS length"): + gate_mpo.multiply(state, compress=False) + + +def test_multiply_invalid_target_type_raises() -> None: + """``multiply`` rejects non-MPS/MPO targets.""" + mpo = MPO.identity(2) + with pytest.raises(TypeError, match="multiply expects MPS or MPO"): + mpo.multiply(object()) # ty: ignore[no-matching-overload] + + +def test_multiply_mpo_embedded_start_site() -> None: + """A shorter gate MPO can be embedded at a non-zero ``start_site``.""" + target = MPO.identity(4) + gate = GateLibrary.cx() + gate.set_sites(0, 1) + support = MPO.from_gate(gate, 2) + support.multiply(target, start_site=1, compress=False) + assert target.length == 4 + target.check_if_valid_mpo() + + +def test_multiply_mpo_with_compression() -> None: + """``multiply(MPO, compress=True)`` runs bond-dimension compression on the target.""" + gate = GateLibrary.cx() + gate.set_sites(0, 1) + support = MPO.from_gate(gate, 2) + target = MPO.identity(2) + support.multiply(target, compress=True, tol=1e-12, n_sweeps=1) + target.check_if_valid_mpo() + + +def test_multiply_mpo_invalid_embed_raises() -> None: + """Embedding a gate MPO outside the target chain raises ``ValueError``.""" + target = MPO.identity(3) + gate = GateLibrary.cx() + gate.set_sites(0, 1) + support = MPO.from_gate(gate, 2) + with pytest.raises(ValueError, match="Cannot embed MPO"): + support.multiply(target, start_site=2, compress=False) + + def test_check_if_identity_non_qubit_physical_dimension() -> None: """Identity check uses the MPO physical dimension, not the qubit default.""" mpo = MPO.identity(2, physical_dimension=3) diff --git a/tests/core/data_structures/test_mpo_utils.py b/tests/core/data_structures/test_mpo_utils.py new file mode 100644 index 00000000..d95b5d86 --- /dev/null +++ b/tests/core/data_structures/test_mpo_utils.py @@ -0,0 +1,146 @@ +# Copyright (c) 2025 - 2026 Chair for Design Automation, TUM +# All rights reserved. +# +# SPDX-License-Identifier: MIT +# +# Licensed under the MIT License + +"""Tests for core MPO tensor utility helpers.""" + +from __future__ import annotations + +from typing import cast + +import numpy as np +import pytest + +from mqt.yaqs.core.data_structures.mpo import MPO +from mqt.yaqs.core.data_structures.mpo_utils import ( + contract_mpo_site_with_mpo_site, + contract_mpo_site_with_mps_site, + decompose_theta, + gate_support_mpo_tensors, + gate_tensor_lr_order, + identity_mpo_site, +) +from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary, extend_gate + + +def test_gate_tensor_lr_order_swaps_sites() -> None: + """Gate tensor axes follow ascending MPS site order.""" + gate = GateLibrary.cx() + gate.set_sites(2, 0) + ordered = gate_tensor_lr_order(gate) + gate.set_sites(0, 2) + direct = gate_tensor_lr_order(gate) + np.testing.assert_allclose(ordered, direct) + + +def test_gate_tensor_lr_order_explicit_sites() -> None: + """Explicit left/right sites select the transpose branch.""" + gate = GateLibrary.cx() + gate.set_sites(3, 1) + tensor = gate_tensor_lr_order(gate, left_site=1, right_site=3) + gate.set_sites(1, 3) + np.testing.assert_allclose(tensor, gate.tensor) + + +def test_gate_tensor_lr_order_inconsistent_sites_raises() -> None: + """Mismatched site indices raise ValueError.""" + gate = GateLibrary.cx() + gate.set_sites(0, 2) + with pytest.raises(ValueError, match="not consistent"): + gate_tensor_lr_order(gate, left_site=0, right_site=3) + + +def test_identity_mpo_site_shape() -> None: + """Identity MPO site has bond dimension one.""" + site = identity_mpo_site(2) + assert site.shape == (2, 2, 1, 1) + np.testing.assert_allclose(site[:, :, 0, 0], np.eye(2)) + + +def test_gate_support_mpo_tensors_uses_cached_mpo_tensors() -> None: + """Cached gate MPO tensors are returned when the support length matches.""" + gate = GateLibrary.cx() + gate.set_sites(1, 3) + first = gate_support_mpo_tensors(gate, first_site=1, last_site=3) + second = gate_support_mpo_tensors(gate, first_site=1, last_site=3) + assert len(first) == len(second) == 3 + np.testing.assert_allclose(first[0], second[0]) + + +def test_gate_support_mpo_tensors_reextends_when_cache_length_mismatches() -> None: + """A cached tensor list with the wrong length triggers ``extend_gate``.""" + gate = GateLibrary.cx() + gate.set_sites(0, 1) + nn_support = gate_support_mpo_tensors(gate, first_site=0, last_site=1) + assert len(nn_support) == 2 + wide = gate_support_mpo_tensors(gate, first_site=0, last_site=2) + assert len(wide) == 3 + + +def test_gate_support_mpo_tensors_calls_extend_gate_without_cache() -> None: + """Gates without ``mpo_tensors`` build support tensors via ``extend_gate``.""" + gate = GateLibrary.cx() + gate.set_sites(0, 1) + tensor = np.asarray(gate.tensor, dtype=np.complex128) + stub = cast("BaseGate", type("GateStub", (), {"interaction": 2, "sites": [0, 1], "tensor": tensor})()) + support = gate_support_mpo_tensors(stub, first_site=0, last_site=1) + expected = extend_gate(gate_tensor_lr_order(gate), [0, 1]) + assert len(support) == 2 + np.testing.assert_allclose(support[0], expected[0]) + np.testing.assert_allclose(support[1], expected[1]) + + +def test_contract_mpo_site_with_mps_site() -> None: + """MPO--MPS site contraction fuses virtual bonds in library order.""" + mpo_site = identity_mpo_site(2) + mps_site = np.ones((2, 1, 1), dtype=np.complex128) + out = contract_mpo_site_with_mps_site(mpo_site, mps_site) + assert out.shape == (2, 1, 1) + + +def test_contract_mpo_site_with_mpo_site_conjugate_flag() -> None: + """MPO--MPO contraction supports the conjugated equivalence-checking path.""" + left = identity_mpo_site(2) + right = identity_mpo_site(2) + plain = contract_mpo_site_with_mpo_site(left, right, conjugate=False) + conj = contract_mpo_site_with_mpo_site(left, right, conjugate=True) + assert plain.shape == conj.shape == (2, 2, 1, 1) + + +def test_decompose_theta_splits_fused_tensor() -> None: + """``decompose_theta`` returns two rank-4 tensors with matching total norm.""" + theta = np.random.default_rng(0).standard_normal((2, 2, 2, 2, 2, 2)) + 1j * np.random.default_rng( + 1 + ).standard_normal((2, 2, 2, 2, 2, 2)) + left, right = decompose_theta(theta, threshold=1e-12) + assert left.ndim == 4 + assert right.ndim == 4 + + +def test_from_gate_rejects_single_qubit_gate() -> None: + """``MPO.from_gate`` requires a two-qubit gate.""" + gate = GateLibrary.x() + gate.set_sites(0) + with pytest.raises(ValueError, match="two-qubit gate"): + MPO.from_gate(gate, 2) + + +def test_from_gate_rejects_chain_shorter_than_support() -> None: + """``chain_length`` must cover the gate support interval.""" + gate = GateLibrary.cx() + gate.set_sites(1, 3) + with pytest.raises(ValueError, match="smaller than gate support"): + MPO.from_gate(gate, 2) + + +def test_from_gate_pads_identity_outside_support() -> None: + """``MPO.from_gate`` inserts identity sites when ``chain_length`` exceeds support width.""" + gate = GateLibrary.cx() + gate.set_sites(1, 3) + mpo = MPO.from_gate(gate, chain_length=6) + assert mpo.length == 6 + assert mpo.tensors[0].shape[-1] == 1 + assert mpo.tensors[5].shape[-1] == 1 diff --git a/tests/core/data_structures/test_mps.py b/tests/core/data_structures/test_mps.py index f503596b..735feaf5 100644 --- a/tests/core/data_structures/test_mps.py +++ b/tests/core/data_structures/test_mps.py @@ -975,7 +975,7 @@ def test_truncate_preserves_orthogonality_center_and_canonicity(center: int) -> assert before_center == center # do a "no-real" truncation (tiny threshold, generous max bond) - mps.truncate(threshold=1e-16, max_bond_dim=100) + mps.compress(threshold=1e-16, max_bond_dim=100) after_center = mps.check_canonical_form()[0] assert after_center == center @@ -1012,7 +1012,7 @@ def test_truncate_reduces_bond_dimensions_and_truncates() -> None: # put it into a known canonical form mps.set_canonical_form(2) # perform a truncation that will cut back to max_bond=3 - mps.truncate(threshold=1e-12, max_bond_dim=3) + mps.compress(threshold=1e-12, max_bond_dim=3) # check validity and that every bond dim <= 3 mps.check_if_valid_mps() @@ -1024,6 +1024,15 @@ def test_truncate_reduces_bond_dimensions_and_truncates() -> None: assert bond_right <= 3 +def test_compress_single_site_returns_immediately() -> None: + """``compress`` is a no-op on a one-site MPS.""" + mps = MPS(1, state="zeros") + before = copy.deepcopy(mps.tensors) + mps.compress(threshold=1e-12) + for before_tensor, after_tensor in zip(before, mps.tensors, strict=True): + np.testing.assert_allclose(before_tensor, after_tensor) + + def _bell_pair_mps() -> MPS: """Auxiliary function to create a Bell-pair MPS. diff --git a/tests/core/data_structures/test_simulation_parameters.py b/tests/core/data_structures/test_simulation_parameters.py index d32ecfe1..fd8e8921 100644 --- a/tests/core/data_structures/test_simulation_parameters.py +++ b/tests/core/data_structures/test_simulation_parameters.py @@ -161,14 +161,15 @@ def test_weak_simparams_default_constructor_uses_balanced() -> None: assert params.svd_threshold == pytest.approx(balanced["svd_threshold"]) assert params.max_bond_dim == balanced["max_bond_dim"] assert params.krylov_tol == pytest.approx(balanced["krylov_tol"]) - assert params.gate_mode == "tebd" + assert params.gate_mode == "mpo" def test_gate_mode_defaults_and_validation() -> None: - """Strong and weak digital params default to TEBD and validate gate_mode names.""" - assert StrongSimParams().gate_mode == "tebd" - assert WeakSimParams(shots=1).gate_mode == "tebd" - assert StrongSimParams(gate_mode="tdvp").gate_mode == "tdvp" + """Strong and weak digital params default to mpo and validate gate_mode names.""" + assert StrongSimParams().gate_mode == "mpo" + assert WeakSimParams(shots=1).gate_mode == "mpo" + assert StrongSimParams(gate_mode="full-tdvp").gate_mode == "full-tdvp" + assert StrongSimParams(gate_mode="mpo").gate_mode == "mpo" with pytest.raises(ValueError, match="gate_mode"): StrongSimParams(gate_mode=cast("GateMode", "invalid")) diff --git a/tests/digital/test_digital_tjm.py b/tests/digital/test_digital_tjm.py index aa90e2a5..a0a57ac1 100644 --- a/tests/digital/test_digital_tjm.py +++ b/tests/digital/test_digital_tjm.py @@ -45,12 +45,14 @@ from mqt.yaqs.core.libraries.circuit_library import create_ising_circuit from mqt.yaqs.core.libraries.gate_library import GateLibrary, X, Y, Z from mqt.yaqs.digital.digital_tjm import ( + apply_long_range_gate_mpo, apply_single_qubit_gate, apply_two_qubit_gate, apply_two_qubit_gate_tebd, apply_window, construct_generator_mpo, create_local_noise_model, + digital_tjm, process_layer, ) from mqt.yaqs.digital.utils.dag_utils import convert_dag_to_tensor_algorithm @@ -656,7 +658,7 @@ def test_ignores_non_mid_barriers_and_handles_measures() -> None: def _run_strong_noiseless( qc: QuantumCircuit, *, - gate_mode: GateMode = "hybrid", + gate_mode: GateMode = "tdvp", max_bond_dim: int | None = None, svd_threshold: float = 1e-12, get_state: bool = False, @@ -685,7 +687,7 @@ def _run_strong_noiseless( return float(np.real(exp[0])) -@pytest.mark.parametrize("gate_mode", ["hybrid", "tebd"]) +@pytest.mark.parametrize("gate_mode", ["tdvp", "swaps"]) def test_nearest_neighbor_gate_modes_agree(gate_mode: str) -> None: """Hybrid and TEBD agree on a purely nearest-neighbor circuit.""" qc = QuantumCircuit(3) @@ -694,7 +696,7 @@ def test_nearest_neighbor_gate_modes_agree(gate_mode: str) -> None: qc.cz(1, 2) qc.rx(0.3, 2) - hybrid_z = _run_strong_noiseless(qc, gate_mode="hybrid") + hybrid_z = _run_strong_noiseless(qc, gate_mode="tdvp") other_z = _run_strong_noiseless(qc, gate_mode=cast("GateMode", gate_mode)) assert hybrid_z == pytest.approx(other_z, abs=1e-12) @@ -705,19 +707,91 @@ def test_long_range_hybrid_matches_tdvp() -> None: qc.h(0) qc.cx(0, 2) - hybrid_z = _run_strong_noiseless(qc, gate_mode="hybrid") - tdvp_z = _run_strong_noiseless(qc, gate_mode="tdvp") + hybrid_z = _run_strong_noiseless(qc, gate_mode="tdvp") + tdvp_z = _run_strong_noiseless(qc, gate_mode="full-tdvp") assert hybrid_z == pytest.approx(tdvp_z, abs=1e-10) +def test_long_range_swap_path_reversed_control_target() -> None: + """TEBD swap routing handles gates with descending site indices (``cx(3, 0)``).""" + qc = QuantumCircuit(4) + qc.h(0) + qc.cx(3, 0) + swaps_z = _run_strong_noiseless(qc, gate_mode="swaps") + mpo_z = _run_strong_noiseless(qc, gate_mode="mpo") + assert swaps_z == pytest.approx(mpo_z, abs=1e-10) + + +def test_weak_noiseless_get_state_returns_mps() -> None: + """Noiseless weak simulation with ``get_state=True`` returns the final MPS.""" + qc = QuantumCircuit(2) + qc.h(0) + qc.measure_all() + sim_params = WeakSimParams(shots=4, gate_mode="tdvp", preset="exact", get_state=True) + state = State(2, initial="zeros") + result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) + assert result.counts is not None + assert sum(result.counts.values()) == 4 + assert result.output_state is not None + assert result.output_state.mps.length == 2 + + +def test_digital_tjm_weak_noisy_get_state_returns_mps() -> None: + """Noisy weak ``digital_tjm`` may return the evolved MPS when ``get_state=True``.""" + mps = MPS(2, state="zeros") + qc = QuantumCircuit(2) + qc.h(0) + noise_model = NoiseModel([{"name": "pauli_x", "sites": [0], "strength": 0.01}]) + sim_params = WeakSimParams(shots=1, gate_mode="tdvp", preset="exact", get_state=True) + counts, _, final = digital_tjm((0, mps, noise_model, sim_params, qc)) + assert isinstance(counts, dict) + assert final is not None + assert final.length == 2 + + +def test_zip_up_nearest_neighbor_matches_tebd() -> None: + """Zip-up uses TEBD on nearest-neighbor gates, matching an all-TEBD run.""" + qc = QuantumCircuit(3) + qc.h(0) + qc.cx(0, 1) + qc.cz(1, 2) + qc.rx(0.3, 2) + + tebd_z = _run_strong_noiseless(qc, gate_mode="swaps") + zip_up_z = _run_strong_noiseless(qc, gate_mode="mpo") + assert zip_up_z == pytest.approx(tebd_z, abs=1e-12) + + +def test_long_range_zip_up_matches_tdvp() -> None: + """Long-range gates use MPO mode and match full-tdvp.""" + qc = QuantumCircuit(4) + qc.h(0) + qc.cx(0, 2) + + zip_up_z = _run_strong_noiseless(qc, gate_mode="mpo") + tdvp_z = _run_strong_noiseless(qc, gate_mode="full-tdvp") + assert zip_up_z == pytest.approx(tdvp_z, abs=1e-10) + + +def test_long_range_zip_up_matches_tebd() -> None: + """Zip-up on long-range gates matches SWAP+TEBD on small circuits.""" + qc = QuantumCircuit(4) + qc.h(0) + qc.cx(0, 2) + + zip_up_z = _run_strong_noiseless(qc, gate_mode="mpo") + tebd_z = _run_strong_noiseless(qc, gate_mode="swaps") + assert zip_up_z == pytest.approx(tebd_z, abs=1e-10) + + def test_long_range_tebd_matches_tdvp() -> None: """TEBD with SWAP insertion matches all-TDVP on a long-range gate.""" qc = QuantumCircuit(4) qc.h(0) qc.cx(0, 2) - tebd_z = _run_strong_noiseless(qc, gate_mode="tebd") - tdvp_z = _run_strong_noiseless(qc, gate_mode="tdvp") + tebd_z = _run_strong_noiseless(qc, gate_mode="swaps") + tdvp_z = _run_strong_noiseless(qc, gate_mode="full-tdvp") assert tebd_z == pytest.approx(tdvp_z, abs=1e-10) @@ -733,7 +807,7 @@ def test_long_range_tebd_matches_qiskit_statevector() -> None: # Use a stricter SVD threshold than the helper default to reduce purely numerical drift # from repeated TEBD split/merge operations (including SWAP insertion). - tebd_vec = _run_strong_noiseless(qc, gate_mode="tebd", get_state=True, svd_threshold=1e-14) + tebd_vec = _run_strong_noiseless(qc, gate_mode="swaps", get_state=True, svd_threshold=1e-14) assert isinstance(tebd_vec, np.ndarray) ref = np.asarray(Statevector(qc).data, dtype=np.complex128) @@ -748,8 +822,8 @@ def test_mixed_circuit_tebd_matches_tdvp() -> None: qc.cx(0, 3) qc.rzz(0.2, 2, 3) - tebd_z = _run_strong_noiseless(qc, gate_mode="tebd") - tdvp_z = _run_strong_noiseless(qc, gate_mode="tdvp") + tebd_z = _run_strong_noiseless(qc, gate_mode="swaps") + tdvp_z = _run_strong_noiseless(qc, gate_mode="full-tdvp") assert tebd_z == pytest.approx(tdvp_z, abs=1e-10) @@ -784,7 +858,7 @@ def test_statevector_matches_qiskit_on_nonsymmetric_probes() -> None: circuits.append(qc) for qc in circuits: - vec = _run_strong_noiseless(qc, gate_mode="hybrid", get_state=True) + vec = _run_strong_noiseless(qc, gate_mode="tdvp", get_state=True) assert isinstance(vec, np.ndarray) ref = np.asarray(Statevector(qc).data, dtype=np.complex128) assert _fidelity(ref, vec) == pytest.approx(1.0, abs=1e-10) @@ -816,7 +890,7 @@ def qiskit_single_pauli_expectation(qc: QuantumCircuit, site: int, op: str) -> f Observable(Z(), 0), Observable(X(), 2), ] - sim_params = StrongSimParams(observables=requested, gate_mode="hybrid", preset="exact") + sim_params = StrongSimParams(observables=requested, gate_mode="tdvp", preset="exact") state = State(3, initial="zeros") result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) @@ -856,7 +930,7 @@ def qiskit_single_pauli_expectation_from_vec(vec: np.ndarray, site: int, op: str Observable(Z(), 0), Observable(X(), 2), ] - sim_params = StrongSimParams(observables=requested, gate_mode="hybrid", preset="exact", get_state=True) + sim_params = StrongSimParams(observables=requested, gate_mode="tdvp", preset="exact", get_state=True) state = State(3, initial="zeros") result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) @@ -882,7 +956,7 @@ def test_expectation_values_align_with_result_observables_order() -> None: # Intentionally not sorted by site: Result order should still match this list. requested = [Observable(Z(), 2), Observable(X(), 0), Observable(Z(), 0)] - sim_params = StrongSimParams(observables=requested, gate_mode="hybrid", preset="exact", get_state=True) + sim_params = StrongSimParams(observables=requested, gate_mode="tdvp", preset="exact", get_state=True) state = State(3, initial="zeros") result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) @@ -913,8 +987,8 @@ def test_mixed_circuit_hybrid_matches_tdvp() -> None: qc.cx(0, 3) qc.rzz(0.2, 2, 3) - hybrid_z = _run_strong_noiseless(qc, gate_mode="hybrid") - tdvp_z = _run_strong_noiseless(qc, gate_mode="tdvp") + hybrid_z = _run_strong_noiseless(qc, gate_mode="tdvp") + tdvp_z = _run_strong_noiseless(qc, gate_mode="full-tdvp") assert hybrid_z == pytest.approx(tdvp_z, abs=1e-10) @@ -924,7 +998,7 @@ def test_bell_state_sanity() -> None: qc.h(0) qc.cx(0, 1) - vec = _run_strong_noiseless(qc, gate_mode="hybrid", get_state=True) + vec = _run_strong_noiseless(qc, gate_mode="tdvp", get_state=True) assert isinstance(vec, np.ndarray) probs = np.abs(vec) ** 2 np.testing.assert_allclose(probs[0], 0.5, atol=1e-10) @@ -943,7 +1017,7 @@ def test_tebd_truncation_respects_max_bond_dim() -> None: state = State(4, initial="zeros") sim_params = StrongSimParams( observables=[Observable(Z(), 0)], - gate_mode="tebd", + gate_mode="swaps", max_bond_dim=2, svd_threshold=1e-6, ) @@ -960,7 +1034,7 @@ def test_weak_simulation_nearest_neighbor_counts() -> None: qc.measure_all() state = State(3, initial="zeros") - sim_params = WeakSimParams(shots=32, gate_mode="hybrid", preset="exact") + sim_params = WeakSimParams(shots=32, gate_mode="tdvp", preset="exact") result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) assert result.counts is not None assert sum(result.counts.values()) == 32 @@ -985,7 +1059,7 @@ def test_weak_counts_match_qiskit_qubit_ordering_deterministic(num_qubits: int, qc.measure_all() shots = 32 - sim_params = WeakSimParams(shots=shots, gate_mode="hybrid", preset="exact") + sim_params = WeakSimParams(shots=shots, gate_mode="tdvp", preset="exact") state = State(num_qubits, initial="zeros") result = Simulator(parallel=False, show_progress=False).run(state, qc, sim_params, None) assert result.counts is not None @@ -1008,7 +1082,7 @@ def test_noisy_nearest_neighbor_smoke() -> None: state = State(2, initial="zeros") sim_params = StrongSimParams( observables=[Observable(Z(), 0)], - gate_mode="hybrid", + gate_mode="tdvp", num_traj=4, random_seed=0, ) @@ -1027,7 +1101,7 @@ def test_apply_two_qubit_gate_tebd_direct_cx_long_range() -> None: dag = circuit_to_dag(qc) node = next(n for n in dag.front_layer() if n.op.name.lower() == "cx") - sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="tebd") + sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="swaps") apply_two_qubit_gate_tebd(mps, convert_dag_to_tensor_algorithm(node)[0], sim_params) mps.normalize(decomposition="SVD") for i, element in enumerate(mps.to_vec()): @@ -1037,6 +1111,27 @@ def test_apply_two_qubit_gate_tebd_direct_cx_long_range() -> None: np.testing.assert_allclose(np.abs(element), 0, atol=1e-10) +def test_apply_long_range_gate_mpo_direct_cx_long_range() -> None: + """Zip-up applies CX(1, 3) on |1111> via extended gate MPO.""" + length = 4 + mps = MPS(length, state="ones") + mps.normalize() + + qc = QuantumCircuit(length) + qc.cx(1, 3) + dag = circuit_to_dag(qc) + node = next(n for n in dag.front_layer() if n.op.name.lower() == "cx") + + sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="mpo") + apply_long_range_gate_mpo(mps, convert_dag_to_tensor_algorithm(node)[0], sim_params) + mps.normalize(decomposition="SVD") + for i, element in enumerate(mps.to_vec()): + if i == 7: + np.testing.assert_allclose(np.abs(element), 1, atol=1e-10) + else: + np.testing.assert_allclose(np.abs(element), 0, atol=1e-10) + + def test_apply_two_qubit_gate_tebd_direct_cnot() -> None: """Direct TEBD helper matches CX on |11> (same check as apply_two_qubit_gate).""" length = 4 @@ -1048,7 +1143,7 @@ def test_apply_two_qubit_gate_tebd_direct_cnot() -> None: dag = circuit_to_dag(qc) node = next(n for n in dag.front_layer() if n.op.name.lower() == "cx") - sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="tebd") + sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="swaps") apply_two_qubit_gate_tebd( mps, convert_dag_to_tensor_algorithm(node)[0], diff --git a/tests/digital/test_mps_utils.py b/tests/digital/test_mps_utils.py new file mode 100644 index 00000000..7823ee2c --- /dev/null +++ b/tests/digital/test_mps_utils.py @@ -0,0 +1,228 @@ +# Copyright (c) 2025 - 2026 Chair for Design Automation, TUM +# All rights reserved. +# +# SPDX-License-Identifier: MIT +# +# Licensed under the MIT License + +"""Tests for generic MPO--MPS gate application via :meth:`~mqt.yaqs.core.data_structures.mpo.MPO.multiply`.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +from qiskit import QuantumCircuit +from qiskit.converters import circuit_to_dag +from qiskit.quantum_info import Statevector + +from mqt.yaqs.core.data_structures.mpo import MPO +from mqt.yaqs.core.data_structures.mps import MPS +from mqt.yaqs.core.data_structures.simulation_parameters import Observable, StrongSimParams +from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary, Z +from mqt.yaqs.digital.digital_tjm import apply_long_range_gate_mpo, apply_two_qubit_gate_tebd +from mqt.yaqs.digital.utils.dag_utils import convert_dag_to_tensor_algorithm + +if TYPE_CHECKING: + from numpy.typing import NDArray + + +def _sim_params() -> StrongSimParams: + return StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="mpo") + + +def _gate_from_circuit(qc: QuantumCircuit, *, op_name: str | None = None) -> BaseGate: + dag = circuit_to_dag(qc) + if op_name is None: + node = next(n for n in dag.op_nodes()) + else: + node = next(n for n in dag.op_nodes() if n.op.name.lower() == op_name.lower()) + return convert_dag_to_tensor_algorithm(node)[0] + + +def _apply_mpo_reference( + length: int, + gate: BaseGate, + *, + compress: bool, +) -> MPS: + """Apply the extended gate MPO on the full chain (reference path). + + Returns: + MPS after explicit :meth:`~mqt.yaqs.core.data_structures.mpo.MPO.multiply`. + """ + state = MPS(length, state="ones") + state.normalize() + MPO.from_gate(gate, length).multiply( + state, + sim_params=_sim_params() if compress else None, + compress=compress, + ) + return state + + +def _qiskit_evolved_vec(qc: QuantumCircuit, label: str) -> NDArray[np.complex128]: + initial = Statevector.from_label(label).data + return np.asarray(Statevector(initial).evolve(qc).data, dtype=np.complex128) + + +def test_identity_mpo_preserves_statevector() -> None: + """Identity MPO on the full chain leaves the dense state unchanged.""" + length = 4 + state = MPS(length, state="ones") + state.normalize() + expected = np.asarray(state.to_vec(), dtype=np.complex128) + + identity_mpo = MPO.identity(length) + identity_mpo.multiply(state, sim_params=_sim_params(), compress=True) + + np.testing.assert_allclose(state.to_vec(), expected, atol=1e-10) + + +def test_nearest_neighbor_cx_mpo_matches_tebd() -> None: + """Adjacent CX via extended MPO matches direct TEBD application.""" + length = 4 + qc = QuantumCircuit(length) + qc.cx(1, 2) + gate = _gate_from_circuit(qc) + sim_params = _sim_params() + + mpo_path = MPS(length, state="ones") + mpo_path.normalize() + MPO.from_gate(gate, length).multiply(mpo_path, sim_params=sim_params, compress=True) + + tebd_path = MPS(length, state="ones") + tebd_path.normalize() + apply_two_qubit_gate_tebd(tebd_path, gate, sim_params) + + np.testing.assert_allclose(mpo_path.to_vec(), tebd_path.to_vec(), atol=1e-10) + + +def test_long_range_cx_matches_statevector_reference() -> None: + """Long-range CX via MPO--MPS matches a dense statevector reference.""" + length = 4 + qc = QuantumCircuit(length) + qc.cx(1, 3) + gate = _gate_from_circuit(qc) + expected = _qiskit_evolved_vec(qc, "1111") + + state = MPS(length, state="ones") + state.normalize() + apply_long_range_gate_mpo(state, gate, _sim_params()) + np.testing.assert_allclose(np.abs(state.to_vec()), np.abs(expected), atol=1e-10) + + +def test_long_range_cx_ones_state() -> None: + """CX(1, 3) on |1111> maps to |1111> (index 7).""" + length = 4 + qc = QuantumCircuit(length) + qc.cx(1, 3) + gate = _gate_from_circuit(qc) + state = MPS(length, state="ones") + state.normalize() + apply_long_range_gate_mpo(state, gate, _sim_params()) + state.normalize(decomposition="SVD") + for index, element in enumerate(state.to_vec()): + if index == 7: + np.testing.assert_allclose(np.abs(element), 1, atol=1e-10) + else: + np.testing.assert_allclose(np.abs(element), 0, atol=1e-10) + + +def test_directional_long_range_vs_nearest_neighbor_cnot() -> None: + """Long-range and nearest-neighbor CNOTs differ and both match TEBD.""" + length = 4 + label = "0110" + sim_params = _sim_params() + + qc_long = QuantumCircuit(length) + qc_long.cx(1, 3) + gate_long = _gate_from_circuit(qc_long) + long_range = MPS(length, state="basis", basis_string=label) + long_range.normalize() + apply_long_range_gate_mpo(long_range, gate_long, sim_params) + + qc_nn = QuantumCircuit(length) + qc_nn.cx(1, 2) + gate_nn = _gate_from_circuit(qc_nn) + nearest = MPS(length, state="basis", basis_string=label) + nearest.normalize() + apply_long_range_gate_mpo(nearest, gate_nn, sim_params) + + tebd_long = MPS(length, state="basis", basis_string=label) + tebd_long.normalize() + apply_two_qubit_gate_tebd(tebd_long, gate_long, sim_params) + + tebd_nn = MPS(length, state="basis", basis_string=label) + tebd_nn.normalize() + apply_two_qubit_gate_tebd(tebd_nn, gate_nn, sim_params) + + np.testing.assert_allclose(np.abs(long_range.to_vec()), np.abs(tebd_long.to_vec()), atol=1e-10) + np.testing.assert_allclose(np.abs(nearest.to_vec()), np.abs(tebd_nn.to_vec()), atol=1e-10) + assert np.max(np.abs(np.abs(long_range.to_vec()) - np.abs(nearest.to_vec()))) > 0.5 + + +def test_apply_long_range_gate_mpo_matches_mpo_reference() -> None: + """Zip-up entry point matches explicit apply-compress MPO reference.""" + length = 4 + qc = QuantumCircuit(length) + qc.cx(1, 3) + gate = _gate_from_circuit(qc) + reference = _apply_mpo_reference(length, gate, compress=True) + + state = MPS(length, state="ones") + state.normalize() + apply_long_range_gate_mpo(state, gate, _sim_params()) + np.testing.assert_allclose(state.to_vec(), reference.to_vec(), atol=1e-10) + + +def test_apply_long_range_gate_mpo_wide_cx_n18() -> None: + """Wide CX(0, n-1) on 18 qubits completes without error.""" + length = 18 + qc = QuantumCircuit(length) + qc.h(0) + qc.cx(0, length - 1) + gate = _gate_from_circuit(qc, op_name="cx") + mps = MPS(length, state="zeros") + apply_long_range_gate_mpo(mps, gate, _sim_params()) + + +def test_apply_long_range_gate_mpo_wide_cx_n32() -> None: + """Wide CX(0, n-1) on 32 qubits completes without label-budget errors.""" + length = 32 + qc = QuantumCircuit(length) + qc.h(0) + qc.cx(0, length - 1) + gate = _gate_from_circuit(qc, op_name="cx") + mps = MPS(length, state="zeros") + apply_long_range_gate_mpo(mps, gate, _sim_params()) + + +def test_swap_via_mpo_matches_tebd() -> None: + """Non-symmetric SWAP on adjacent sites matches TEBD.""" + length = 4 + swap = GateLibrary.swap() + swap.set_sites(1, 2) + sim_params = StrongSimParams(observables=[Observable(Z(), 0)], preset="exact", gate_mode="swaps") + + mpo_path = MPS(length, state="basis", basis_string="1010") + mpo_path.normalize() + MPO.from_gate(swap, length).multiply(mpo_path, sim_params=sim_params, compress=True) + + tebd_path = MPS(length, state="basis", basis_string="1010") + tebd_path.normalize() + apply_two_qubit_gate_tebd(tebd_path, swap, sim_params) + + np.testing.assert_allclose(mpo_path.to_vec(), tebd_path.to_vec(), atol=1e-10) + + +def test_from_gate_reuses_mpo_tensors() -> None: + """from_gate matches extend_gate when gate mpo_tensors are already cached.""" + length = 4 + qc = QuantumCircuit(length) + qc.cx(1, 3) + gate = _gate_from_circuit(qc) + + first = MPO.from_gate(gate, length) + second = MPO.from_gate(gate, length) + np.testing.assert_allclose(first.to_matrix(), second.to_matrix(), atol=1e-12) diff --git a/tests/digital/utils/test_mpo_utils.py b/tests/digital/utils/test_contraction_utils.py similarity index 82% rename from tests/digital/utils/test_mpo_utils.py rename to tests/digital/utils/test_contraction_utils.py index 4b1269ad..f1f81410 100644 --- a/tests/digital/utils/test_mpo_utils.py +++ b/tests/digital/utils/test_contraction_utils.py @@ -5,9 +5,9 @@ # # Licensed under the MIT License -"""Tests for the MPO utility functions used in the equivalence checking framework. +"""Tests for contraction utilities used in the equivalence checking framework. -This module contains unit tests for the MPO utility functions used in the equivalence checking framework. +This module contains unit tests for the digital contraction utilities used in the equivalence checking framework. It verifies the correct functionality of tensor operations including: - SVD-based splitting of MPS tensors (decompose_theta) - Gate application routines (apply_gate, apply_temporal_zone) @@ -25,7 +25,7 @@ import copy from concurrent.futures import ThreadPoolExecutor -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, cast import numpy as np import pytest @@ -33,19 +33,20 @@ from qiskit.converters import circuit_to_dag from mqt.yaqs.core.data_structures.mpo import MPO +from mqt.yaqs.core.data_structures.mpo_utils import decompose_theta from mqt.yaqs.core.libraries.circuit_library import create_ising_circuit -from mqt.yaqs.core.libraries.gate_library import GateLibrary -from mqt.yaqs.digital.utils.dag_utils import convert_dag_to_tensor_algorithm, get_temporal_zone, select_starting_point -from mqt.yaqs.digital.utils.mpo_utils import ( +from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary +from mqt.yaqs.digital.utils.contraction_utils import ( MIN_QUBITS_FOR_MPO_PARALLEL, apply_gate, apply_layer, apply_long_range_layer, apply_temporal_zone, compute_pair_update, - decompose_theta, + iterate, update_mpo, ) +from mqt.yaqs.digital.utils.dag_utils import convert_dag_to_tensor_algorithm, get_temporal_zone, select_starting_point if TYPE_CHECKING: from numpy.typing import NDArray @@ -138,6 +139,21 @@ def test_decompose_theta() -> None: approximate_reconstruction(tensor1, tensor2, theta, atol=1e-5) +def test_apply_gate_identity_is_noop() -> None: + """Identity gates skip contraction and return ``theta`` unchanged.""" + theta = random_theta_6d() + gate = cast( + "BaseGate", + type( + "IdentityGate", + (), + {"name": "I", "interaction": 1, "sites": [0], "matrix": np.eye(2, dtype=np.complex128)}, + )(), + ) + updated = apply_gate(gate, theta, site0=0, site1=1, conjugate=False) + np.testing.assert_allclose(updated, theta) + + @pytest.mark.parametrize("conjugate", [False, True]) def test_apply_single_qubit_gate(*, conjugate: bool) -> None: """Test applying a single-qubit gate (X gate) to a tensor using apply_gate. @@ -290,6 +306,74 @@ def test_apply_long_range_layer() -> None: assert mpo.check_if_identity(1 - 1e-6), "MPO should approximate identity after long-range layer." +def test_apply_long_range_layer_skips_leading_single_qubit_gate() -> None: + """The first layer may begin with a single-qubit gate before a long-range CX.""" + mpo = MPO.identity(3) + circuit = QuantumCircuit(3) + circuit.h(1) + circuit.cx(0, 2) + dag1 = circuit_to_dag(circuit) + dag2 = copy.deepcopy(dag1) + apply_long_range_layer(mpo, dag1, dag2, conjugate=False, threshold=1e-12) + assert mpo.check_if_valid_mpo() + + +def test_apply_long_range_layer_on_wider_mpo() -> None: + """Long-range updates embed gate support inside a longer identity MPO.""" + mpo = MPO.identity(5) + circuit = QuantumCircuit(5) + circuit.cx(0, 2) + dag1 = circuit_to_dag(circuit) + dag2 = copy.deepcopy(dag1) + apply_long_range_layer(mpo, dag1, dag2, conjugate=False, threshold=1e-12) + assert mpo.length == 5 + assert mpo.check_if_valid_mpo() + + +def test_iterate_serial_and_parallel() -> None: + """``iterate`` drives checkerboard sweeps with and without a thread pool.""" + qc = QuantumCircuit(MIN_QUBITS_FOR_MPO_PARALLEL) + qc.cx(0, 2) + dag1 = circuit_to_dag(qc) + + mpo_serial = MPO.identity(MIN_QUBITS_FOR_MPO_PARALLEL) + iterate(mpo_serial, copy.deepcopy(dag1), copy.deepcopy(dag1), threshold=1e-12, parallel=False) + assert mpo_serial.check_if_identity(1 - 1e-6) + + mpo_parallel = MPO.identity(MIN_QUBITS_FOR_MPO_PARALLEL) + iterate( + mpo_parallel, + circuit_to_dag(qc), + circuit_to_dag(qc), + threshold=1e-12, + parallel=True, + max_workers=2, + ) + assert mpo_parallel.check_if_identity(1 - 1e-6) + + +def test_apply_layer_parallel_single_pair_uses_serial_path() -> None: + """Parallel sweeps with one pair fall back to the serial worker path.""" + qc = QuantumCircuit(2) + qc.cx(0, 1) + dag1 = circuit_to_dag(qc) + dag2 = circuit_to_dag(qc) + mpo = MPO.identity(2) + first_iterator, second_iterator = select_starting_point(2, dag1) + with ThreadPoolExecutor(max_workers=2) as pool: + apply_layer( + mpo, + dag1, + dag2, + first_iterator, + second_iterator, + 1e-12, + parallel=True, + thread_pool=pool, + ) + assert mpo.check_if_valid_mpo() + + def _make_n_by_n_circuit(num_qubits: int) -> QuantumCircuit: """Build an ``n`` x ``n`` layered circuit (``n`` qubits, ``n`` repetitions). diff --git a/tests/test_equivalence_checker.py b/tests/test_equivalence_checker.py index d6e71552..9815954f 100644 --- a/tests/test_equivalence_checker.py +++ b/tests/test_equivalence_checker.py @@ -23,7 +23,7 @@ from qiskit.qasm2 import load from mqt.yaqs import DEFAULT_MATRIX_MAX_QUBITS, EquivalenceChecker -from mqt.yaqs.digital.utils.mpo_utils import MIN_QUBITS_FOR_MPO_PARALLEL +from mqt.yaqs.digital.utils.contraction_utils import MIN_QUBITS_FOR_MPO_PARALLEL if TYPE_CHECKING: from mqt.yaqs.equivalence_checker import Representation @@ -229,6 +229,30 @@ def test_checker_rejects_non_int_max_workers() -> None: EquivalenceChecker(max_workers=1.5) # ty: ignore[invalid-argument-type] +def test_checker_rejects_invalid_representation() -> None: + """Unknown ``representation`` strings are rejected at construction.""" + with pytest.raises(ValueError, match="representation must be one of"): + EquivalenceChecker(representation="tensor") # ty: ignore[invalid-argument-type] + + +def test_checker_rejects_bool_matrix_max_qubits() -> None: + """``matrix_max_qubits`` must be a true integer, not a boolean.""" + with pytest.raises(TypeError, match="matrix_max_qubits"): + EquivalenceChecker(matrix_max_qubits=True) + + +def test_checker_rejects_negative_matrix_max_qubits() -> None: + """``matrix_max_qubits`` must be non-negative.""" + with pytest.raises(ValueError, match="non-negative"): + EquivalenceChecker(matrix_max_qubits=-1) + + +def test_check_rejects_mismatched_qubit_counts() -> None: + """``check`` requires both circuits to have the same width.""" + with pytest.raises(ValueError, match="same number of qubits"): + EquivalenceChecker().check(QuantumCircuit(2), QuantumCircuit(3)) + + def test_equivalence_checker_defaults_parallel_true() -> None: """``parallel`` defaults to ``True`` (MPO thread pool still gated by qubit count).""" assert EquivalenceChecker().parallel is True