diff --git a/CHANGELOG.md b/CHANGELOG.md index 655415780..360462c9f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel ### Changed +- Changed Pauli summation to use finite state machine construction ([`#308`]) ([**@aaronleesander**]) - 🔧 Replace `mypy` with `ty` ([#304]) ([**@denialhaag**]) ## [0.3.3] - 2026-01-12 @@ -60,6 +61,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool +[#308]: https://github.com/munich-quantum-toolkit/yaqs/pull/308 [#304]: https://github.com/munich-quantum-toolkit/yaqs/pull/304 [#298]: https://github.com/munich-quantum-toolkit/yaqs/pull/298 [#216]: https://github.com/munich-quantum-toolkit/yaqs/pull/216 diff --git a/src/mqt/yaqs/core/data_structures/networks.py b/src/mqt/yaqs/core/data_structures/networks.py index 33e5c0926..d3f5a5b0d 100644 --- a/src/mqt/yaqs/core/data_structures/networks.py +++ b/src/mqt/yaqs/core/data_structures/networks.py @@ -1376,13 +1376,15 @@ def from_pauli_sum( Each term is given as ``(coeff, spec)`` where ``spec`` is a string like ``"Z0 Z1"``, ``"X7"``, or ``""`` for the identity. Terms are assembled by - direct-summing their MPO bond spaces and then compressed in-place. + constructing a finite state machine (FSM) that represents the sum of terms + directly, resulting in an optimal or near-optimal bond dimension without + intermediate compression steps. Args: terms: List of ``(coefficient, spec)`` Pauli terms. length: Number of sites (L). physical_dimension: Local dimension (only ``2`` is supported). - tol: SVD truncation threshold used during compression. + tol: SVD truncation threshold used during final compression. max_bond_dim: Optional hard cap on the kept MPO bond dimension. n_sweeps: Number of compression sweeps (>= 0). @@ -1392,7 +1394,8 @@ def from_pauli_sum( Notes: The resulting MPO represents the sum of all provided terms (including - coefficients). Compression is performed via ``self.compress(..., directions="lr_rl")``. + coefficients). The construction uses an FSM approach which is significantly + faster than summing individual MPOs for large numbers of terms. """ if physical_dimension != 2: msg = "Only physical_dimension=2 is supported by this Pauli MPO builder." @@ -1404,27 +1407,107 @@ def from_pauli_sum( self.length = length self.physical_dimension = physical_dimension - mpo_terms: list[list[np.ndarray]] = [] + if not terms: + self.tensors = [np.zeros((2, 2, 1, 1), dtype=complex) for _ in range(length)] + return + + # 1. Parse terms into dense lists of operator names. + # Structure: terms list of (coeff, [op_at_site_0, op_at_site_1, ...]) + parsed_terms: list[tuple[complex | float, list[str]]] = [] for coeff, spec in terms: - ops = self._parse_pauli_string(spec) # {site: "X|Y|Z|I"} - # validate sites + ops - for site, lab in ops.items(): + ops_map = self._parse_pauli_string(spec) + # Validate sites + for site, lab in ops_map.items(): if not (0 <= site < length): msg = f"Site index {site} outside [0, {length - 1}]." raise ValueError(msg) if lab not in self._VALID: msg = f"Invalid local op {lab!r}; expected one of {sorted(self._VALID)}." raise ValueError(msg) - mpo_terms.append(self._create_term(length, ops, coeff)) - if not mpo_terms: - self.tensors = [np.zeros((2, 2, 1, 1), dtype=complex) for _ in range(length)] - else: - acc = mpo_terms[0] - for nxt in mpo_terms[1:]: - acc = self._sum_terms(acc, nxt) - self.tensors = acc + # Fill missing sites with Identity "I" + op_list = [ops_map.get(i, "I") for i in range(length)] + parsed_terms.append((coeff, op_list)) + + # 2. Assign State IDs (Right-to-Left) + # We identify unique "suffix states" needed at each bond. + # A state at bond i is uniquely defined by the pair (Operator at site i, State at bond i+1). + + # `term_trajectories[term_idx][i]` stores the State ID at bond `i` for `term_idx`. + # Bond indices range from 0 (left of site 0) to L (right of site L-1). + term_trajectories = [[0] * (length + 1) for _ in range(len(parsed_terms))] + + # Initialize right boundary (Bond L): All terms end at the "sink" state (ID 0). + for t_idx in range(len(parsed_terms)): + term_trajectories[t_idx][length] = 0 + + # bond_state_maps[i] stores the mapping: (Op_str, Next_State_ID) -> Current_State_ID + bond_state_maps: list[dict[tuple[str, int], int]] = [{} for _ in range(length + 1)] + + # Sweep Right-to-Left (sites L-1 down to 1) to build the FSM transitions. + # We stop at bond 1. Bond 0 is always the single "Start" state. + for i in range(length - 1, 0, -1): + next_bond = i + 1 + current_bond = i + + unique_states_map = bond_state_maps[current_bond] + next_id = 0 + + for t_idx, (_, ops) in enumerate(parsed_terms): + op = ops[i] + next_state = term_trajectories[t_idx][next_bond] + signature = (op, next_state) + + if signature not in unique_states_map: + unique_states_map[signature] = next_id + next_id += 1 + + term_trajectories[t_idx][current_bond] = unique_states_map[signature] + + # 3. Build Tensors (Left-to-Right) + self.tensors = [] + paulis = self._PAULI_2 + + for i in range(length): + # Determine bond dimensions based on number of unique states at boundaries + if i == 0: + d_left = 1 + d_right = 1 if length == 1 else len(bond_state_maps[1]) + # Handle edge case where d_right is 0 (should not happen if terms exist) + if length > 1 and d_right == 0: + d_right = 1 + else: + d_left = len(bond_state_maps[i]) + d_right = 1 if i == length - 1 else len(bond_state_maps[i + 1]) + + # Allocate tensor: (phys_out, phys_in, left, right) + tensor = np.zeros((2, 2, d_left, d_right), dtype=complex) + + if i == 0: + # First site: Accumulate coefficients and split into initial branches. + for t_idx, (coeff, ops) in enumerate(parsed_terms): + op_name = ops[i] + op_mat = paulis[op_name] + target_state = term_trajectories[t_idx][1] + + # Accumulate contribution. Multiple terms may map to the same target state. + tensor[:, :, 0, target_state] += coeff * op_mat + else: + # Internal sites: deterministic transitions. + # Each row (current_id) in the tensor corresponds to a unique state from Step 2. + # This state maps to exactly one (op, next_id) pair. + map_i = bond_state_maps[i] + for (op_name, next_id), current_id in map_i.items(): + op_mat = paulis[op_name] + tensor[:, :, current_id, next_id] = op_mat + + self.tensors.append(tensor) + + # 4. Final Compression + # The FSM construction is optimal for one-sided (suffix) uniqueness. + # A standard two-sweep compression ("lr_rl") puts the MPO in canonical form + # and removes any remaining redundancies (e.g., common prefixes). self.compress(tol=tol, max_bond_dim=max_bond_dim, n_sweeps=n_sweeps, directions="lr_rl") assert self.check_if_valid_mpo(), "MPO initialized wrong" @@ -1632,39 +1715,6 @@ def check_if_identity(self, fidelity: float) -> bool: # Checks if trace is not a singular values for partial trace return not np.round(np.abs(trace), 1) / 2**self.length < fidelity - @classmethod - def _create_term(cls, length: int, ops: dict[int, str], coeff: complex) -> list[np.ndarray]: - """Construct an MPO tensor list for a single Pauli-string term. - - Builds a length-`length` MPO corresponding to a single Pauli operator - product specified by a mapping from site indices to Pauli labels. - - Each site tensor has bond dimension 1, and the overall MPO represents: - coeff ⨂_i P_i - where P_i ∈ {I, X, Y, Z}. - - Args: - length: Total number of sites in the MPO. - ops: Mapping from site index to Pauli label - ('I', 'X', 'Y', or 'Z'). Sites not present in the dictionary - default to identity. - coeff: Scalar coefficient multiplying the operator. - - Returns: - list[np.ndarray]: List of MPO tensors, one per site, each of shape - (2, 2, 1, 1). - """ - tensors: list[np.ndarray] = [] - p = cls._PAULI_2 - for i in range(length): - lab = ops.get(i, "I") - term = np.zeros((2, 2, 1, 1), dtype=complex) - term[:, :, 0, 0] = p[lab] - tensors.append(term) - tensors[0] = tensors[0].copy() - tensors[0][:, :, 0, 0] *= coeff - return tensors - @classmethod def _parse_pauli_string(cls, spec: str) -> dict[int, str]: """Parse a Pauli-string specification into a site-to-operator mapping. @@ -1706,48 +1756,3 @@ def _parse_pauli_string(cls, spec: str) -> dict[int, str]: msg = f"Invalid token(s) in spec '{spec}'. Use forms like 'X0 Y2 Z5'." raise ValueError(msg) return out - - @staticmethod - def _sum_terms(tensors1: list[np.ndarray], tensors2: list[np.ndarray]) -> list[np.ndarray]: - """Form the direct-sum MPO representing the sum of two MPO terms. - - Given two MPOs `tensors1` and `tensors2` of equal length, constructs a new MPO - representing their sum by block-diagonal concatenation of virtual - bond spaces. - - Args: - tensors1: First MPO term as a list of tensors. - tensors2: Second MPO term as a list of tensors. - - Returns: - list[np.ndarray]: MPO tensors representing tensors1 + tensors2. - - Notes: - - This operation increases the MPO bond dimension additively. - - No compression is performed here; callers are expected to - compress the result afterward if desired. - """ - length = len(tensors1) - assert length == len(tensors2) - if length == 1: - return [tensors1[0] + tensors2[0]] - out: list[np.ndarray] = [] - for k in range(length): - a, b = tensors1[k], tensors2[k] - d = a.shape[0] - la, ra = a.shape[2], a.shape[3] - lb, rb = b.shape[2], b.shape[3] - if k == 0: - c = np.zeros((d, d, 1, ra + rb), dtype=complex) - c[:, :, 0, :ra] = a[:, :, 0, :] - c[:, :, 0, ra:] = b[:, :, 0, :] - elif k == length - 1: - c = np.zeros((d, d, la + lb, 1), dtype=complex) - c[:, :, :la, 0] = a[:, :, :, 0] - c[:, :, la:, 0] = b[:, :, :, 0] - else: - c = np.zeros((d, d, la + lb, ra + rb), dtype=complex) - c[:, :, :la, :ra] = a - c[:, :, la:, ra:] = b - out.append(c) - return out diff --git a/tests/core/data_structures/test_networks.py b/tests/core/data_structures/test_networks.py index cfd12be56..349272a94 100644 --- a/tests/core/data_structures/test_networks.py +++ b/tests/core/data_structures/test_networks.py @@ -1279,34 +1279,6 @@ def test_from_pauli_sum_empty_terms_builds_zero_mpo() -> None: assert np.allclose(t, 0.0) -def test_from_pauli_sum_calls_create_term_for_each_term(monkeypatch: pytest.MonkeyPatch) -> None: - """Pauli-sum MPO construction: each term is converted via _create_term before summation.""" - mpo = MPO() - - # Make parsing deterministic and valid. - monkeypatch.setattr(mpo, "_parse_pauli_string", lambda _spec: {0: "Z"}) - - calls: list[tuple[int, dict[int, str], complex | float]] = [] - - def fake_create_term(length: int, ops: dict[int, str], coeff: complex) -> list[np.ndarray]: - calls.append((length, ops, coeff)) - # Minimal valid term tensor list: identity-ish MPO with bond dim 1 - return [np.zeros((2, 2, 1, 1), dtype=complex) for _ in range(length)] - - monkeypatch.setattr(mpo, "_create_term", fake_create_term) - - # Avoid depending on compress/check behavior in this unit test. - monkeypatch.setattr(mpo, "compress", lambda **_kwargs: None) - monkeypatch.setattr(mpo, "check_if_valid_mpo", lambda: True) - - mpo.from_pauli_sum(terms=[(2.0, "Z0"), (3.0, "Z0")], length=2, n_sweeps=0) - - assert calls == [ - (2, {0: "Z"}, 2.0), - (2, {0: "Z"}, 3.0), - ] - - def test_compress_raises_on_negative_n_sweeps() -> None: """MPO compress input validation: negative n_sweeps must raise.""" mpo = MPO()