# Copyright 2025 Qilimanjaro Quantum Tech
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from __future__ import annotations
import copy
import re
from abc import ABC
from collections import defaultdict
from itertools import combinations, product
from typing import TYPE_CHECKING, Callable, ClassVar, TypeAlias
import numpy as np
from loguru import logger
from scipy.sparse import csr_matrix, kron, spmatrix
from qilisdk.core.expression import Expression
from qilisdk.core.parameterizable import Parameterizable
from qilisdk.core.qtensor import QTensor
from qilisdk.core.types import Number
from qilisdk.settings import get_settings
from qilisdk.utils.hashing import hash as qili_hash
from qilisdk.utils.visualization.style import HamiltonianStyle
from qilisdk.yaml import yaml
from .exceptions import InvalidHamiltonianOperation
if TYPE_CHECKING:
from collections.abc import Callable, Iterator
from qilisdk.core.variables import Parameter
_DIVISION_BY_OPERATORS_MESSAGE = "Division by operators is not supported"
_GENERIC_VARIABLE_IN_TERM_MESSAGE = "Expression provided contains generic variables that are not Parameter."
_GENERIC_VARIABLE_IN_HAMILTONIAN_MESSAGE = (
"Only Parameters are allowed to be used in hamiltonians. Generic Variables are not supported"
)
def _complex_dtype() -> np.dtype:
return np.dtype(get_settings().complex_precision.dtype)
###############################################################################
# Public Factory Functions
###############################################################################
[documentos]
def Z(qubit: int) -> Hamiltonian:
return PauliZ(qubit).to_hamiltonian()
[documentos]
def X(qubit: int) -> Hamiltonian:
return PauliX(qubit).to_hamiltonian()
[documentos]
def Y(qubit: int) -> Hamiltonian:
return PauliY(qubit).to_hamiltonian()
[documentos]
def I(qubit: int = 0) -> Hamiltonian: # ruff: ignore[ambiguous-function-name]
return PauliI(qubit).to_hamiltonian()
###############################################################################
# Abstract Base PauliOperator
###############################################################################
[documentos]
class PauliOperator(ABC):
"""
A generic abstract Pauli operator that acts on one qubit.
Example:
.. code-block:: python
from qilisdk.analog import PauliX
op = PauliX(0)
Note: You can also use the factory functions X(q), Y(q), Z(q), I(q) to get a Hamiltonian object.
"""
_NAME: ClassVar[str]
_MATRIX: ClassVar[np.ndarray]
_MATRIX_CACHE: ClassVar[dict[np.dtype, np.ndarray]]
def __init_subclass__(cls) -> None:
super().__init_subclass__()
cls._MATRIX_CACHE = {}
def __init__(self, qubit: int) -> None:
# QSDK-05: reject negative qubit indices at construction (the upper bound
# depends on the Hamiltonian / circuit and is enforced at execution).
if qubit < 0:
raise ValueError(f"Qubit index must be non-negative, got {qubit}.")
self._qubit = qubit
@property
@property
@classmethod
[documentos]
def matrix_for_dtype(cls, dtype: np.dtype) -> np.ndarray:
cached = cls._MATRIX_CACHE.get(dtype)
if cached is None:
cached = cls._MATRIX.astype(dtype, copy=False)
cls._MATRIX_CACHE[dtype] = cached
return cached
@property
[documentos]
def matrix(self) -> np.ndarray:
return self.matrix_for_dtype(_complex_dtype())
[documentos]
def to_hamiltonian(self) -> Hamiltonian:
"""Convert this single operator to a Hamiltonian with one term.
Returns:
Hamiltonian: The converted Hamiltonian.
"""
return Hamiltonian({(self,): 1})
def __hash__(self) -> int:
return qili_hash(self._NAME, self._qubit)
def __eq__(self, other: object) -> bool:
if isinstance(other, Hamiltonian):
return other == self
if not isinstance(other, PauliOperator):
return False
return (self._NAME == other._NAME) and (self._qubit == other._qubit)
def __repr__(self) -> str:
return f"{self.name}({self.qubit})"
def __str__(self) -> str:
return f"{self.name}({self.qubit})"
# ----------- Arithmetic Operators ------------
def __add__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return self.to_hamiltonian() + other
__radd__ = __add__
__iadd__ = __add__
def __sub__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return self.to_hamiltonian() - other
def __rsub__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return other - self.to_hamiltonian()
__isub__ = __sub__
def __mul__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return self.to_hamiltonian() * other
def __rmul__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return other * self.to_hamiltonian()
__imul__ = __mul__
def __truediv__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
return self.to_hamiltonian() / other
def __rtruediv__(self, _: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
raise InvalidHamiltonianOperation(_DIVISION_BY_OPERATORS_MESSAGE)
__itruediv__ = __truediv__
###############################################################################
# Concrete Operator Classes
###############################################################################
@yaml.register_class
[documentos]
class PauliZ(PauliOperator):
_NAME: ClassVar[str] = "Z"
_MATRIX: ClassVar[np.ndarray] = np.array([[1, 0], [0, -1]], dtype=_complex_dtype())
@yaml.register_class
[documentos]
class PauliX(PauliOperator):
_NAME: ClassVar[str] = "X"
_MATRIX: ClassVar[np.ndarray] = np.array([[0, 1], [1, 0]], dtype=_complex_dtype())
@yaml.register_class
[documentos]
class PauliY(PauliOperator):
_NAME: ClassVar[str] = "Y"
_MATRIX: ClassVar[np.ndarray] = np.array([[0, -1j], [1j, 0]], dtype=_complex_dtype())
@yaml.register_class
[documentos]
class PauliI(PauliOperator):
_NAME: ClassVar[str] = "I"
_MATRIX: ClassVar[np.ndarray] = np.array([[1, 0], [0, 1]], dtype=_complex_dtype())
_PAULI_CLASS_BY_NAME: dict[str, type[PauliOperator]] = {cls._NAME: cls for cls in (PauliI, PauliX, PauliY, PauliZ)}
# Wrapping a lattice direction into a ring only adds a new bond once it spans at least 3 sites;
# below that the wrap-around bond would duplicate one the open lattice already has.
_MIN_PERIODIC_EXTENT = 3
# Coefficients of the model constructors, either one value shared by every term or one value per term
def _validate_nqubits(nqubits: int, minimum: int = 1, label: str = "") -> None:
"""
Validate the number of qubits given for a certain Hamiltonian constructor.
Args:
nqubits (int): the qubit count to validate.
minimum (int, optional): the smallest count the model can be built on, e.g. 2 for a model
with two-qubit couplings. Defaults to 1.
label (str, optional): the model name, used in the error message. Defaults to "".
Raises:
ValueError: if ``nqubits`` is not greater than zero.
ValueError: if ``nqubits`` is below ``minimum``.
"""
if nqubits <= 0:
raise ValueError(f"nqubits must be greater than zero, got {nqubits}.")
if nqubits < minimum:
raise ValueError(f"{label} Hamiltonians need at least {minimum} qubits, got {nqubits}.")
def _parse_coefficients(coefficient: Coefficient, count: int, name: str) -> list[complex]:
"""
Expand the coefficient of a group of Hamiltonian terms into one value per term.
A plain number gives the same value to every term, while a list gives each term the value at its
own position.
Args:
coefficient (Coefficient): the value shared by every term, or one value per term.
count (int): the number of terms the coefficient weights.
name (str): the argument name, used in the error message.
Returns:
list[complex]: the coefficient of each term, in the order the terms are generated.
Raises:
ValueError: if a list is given whose length does not match the number of terms.
"""
if isinstance(coefficient, (list, tuple)):
if len(coefficient) != count:
raise ValueError(f"{name} must hold one coefficient per term, got {len(coefficient)} for {count} terms.")
return [complex(value) for value in coefficient]
return [complex(coefficient)] * count
# Cache sparse single-qubit matrices once per dtype to avoid rebuilding them for every term.
_SINGLE_QUBIT_SPARSE: dict[tuple[str, np.dtype], csr_matrix] = {}
def _get_single_qubit_sparse_matrix(name: str) -> csr_matrix:
dtype = _complex_dtype()
key = (name, dtype)
cached = _SINGLE_QUBIT_SPARSE.get(key)
if cached is None:
pauli_cls = _PAULI_CLASS_BY_NAME[name]
cached = csr_matrix(pauli_cls.matrix_for_dtype(dtype), dtype=dtype)
_SINGLE_QUBIT_SPARSE[key] = cached
return cached
@yaml.register_class
[documentos]
class Hamiltonian(Parameterizable):
"""
Represent a Hamiltonian expressed as a linear combination of Pauli operators.
Example:
.. code-block:: python
from qilisdk.analog.hamiltonian import Hamiltonian, X, Z
H = X(0) * X(1) + Z(1)
"""
_EPS: float = 1e-14
_PAULI_PRODUCT_TABLE: ClassVar[dict[tuple[str, str], tuple[complex, Callable[..., PauliOperator]]]] = {
("X", "X"): (1, PauliI),
("X", "Y"): (1j, PauliZ),
("X", "Z"): (-1j, PauliY),
("Y", "X"): (-1j, PauliZ),
("Y", "Y"): (1, PauliI),
("Y", "Z"): (1j, PauliX),
("Z", "X"): (1j, PauliY),
("Z", "Y"): (-1j, PauliX),
("Z", "Z"): (1, PauliI),
}
def __init__(
self, elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] | None = None
) -> None:
"""
Build a Hamiltonian from a mapping of Pauli operator products to coefficients.
Args:
elements (dict[tuple[PauliOperator, ...], complex | Expression | Parameter], optional):
Mapping from operator tuples to numerical coefficients or symbolic parameters. For example:
.. code-block:: python
{
(Z(0), Y(1)): 1.0,
(X(1),): 1j,
}
Defaults to None, which creates an empty Hamiltonian.
Raises:
ValueError: If the provided coefficients include generic variables instead of parameters.
"""
super().__init__()
self._elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = defaultdict(complex)
if elements:
for key, val in elements.items():
if isinstance(val, Expression):
if not val.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_HAMILTONIAN_MESSAGE)
for parameter in val.free_parameters():
self._add_parameter(parameter.label, parameter)
self._elements[key] += val
self.simplify()
@property
[documentos]
def nqubits(self) -> int:
"""Number of qubits on which the Hamiltonian acts."""
qubits = {op.qubit for key in self._elements for op in key}
return max(qubits) + 1 if qubits else 0
@property
[documentos]
def elements(self) -> dict[tuple[PauliOperator, ...], complex]:
"""Return the stored operator-coefficient mapping with symbolic terms evaluated."""
return {k: (v if isinstance(v, (int, float, complex)) else v.evaluate({})) for k, v in self._elements.items()}
[documentos]
def simplify(self) -> Hamiltonian:
"""Simplify the Hamiltonian expression by removing near-zero terms and accumulating constant terms.
Returns:
Hamiltonian: Simplified Hamiltonian
"""
# 1) Remove near-zero
keys_to_remove = [
key for key, value in self._elements.items() if isinstance(value, complex) and abs(value) < Hamiltonian._EPS
]
for key in keys_to_remove:
del self._elements[key]
# 2) Accumulate identities that do NOT act on qubit=0 => I(0)
to_accumulate = [
(key, value)
for key, value in self._elements.items()
if len(key) == 1 and key[0].name == "I" and key[0].qubit != 0
]
for key, value in to_accumulate:
del self._elements[key]
self._elements[PauliI(0),] += value
return self
def _apply_operator_on_qubit(self, terms: list[PauliOperator], padding: int = 0) -> spmatrix:
"""Get the matrix representation of a single term by taking the tensor product
of operators acting on each qubit. For qubits with no operator in `terms`,
the identity is used.
Args:
terms (list[PauliOperator]): A list of Pauli operators in the term.
Returns:
spmatrix: The full matrix representation of the term.
Raises:
ValueError: If multiple operators act on the same qubit.
"""
total_qubits = self.nqubits + padding
if total_qubits == 0:
return csr_matrix((1, 1), dtype=_complex_dtype())
ordered_terms = sorted(terms, key=lambda op: op.qubit)
# Check we don't have multiple operators on the same qubit
qubit_indices = [op.qubit for op in ordered_terms]
if len(qubit_indices) != len(set(qubit_indices)):
raise ValueError("The list should not contain multiple operators acting on the same qubit.")
identity_single = _get_single_qubit_sparse_matrix("I")
idx = 0
next_op = ordered_terms[0] if ordered_terms else None
result: spmatrix | None = None
for qubit in range(total_qubits):
if next_op is not None and next_op.qubit == qubit:
single = _get_single_qubit_sparse_matrix(next_op.name)
idx += 1
next_op = ordered_terms[idx] if idx < len(ordered_terms) else None
else:
single = identity_single
result = single if result is None else kron(result, single, format="csr")
# Added for type safety, but should never occur
return result if result is not None else csr_matrix((2**total_qubits, 2**total_qubits), dtype=_complex_dtype())
[documentos]
def to_matrix(self) -> spmatrix:
"""Return the full matrix representation of the Hamiltonian by summing over all terms.
Returns:
spmatrix: The sparse matrix representation of the Hamiltonian.
"""
logger.debug("[Hamiltonian] Building sparse matrix representation over {} qubits", self.nqubits)
dim = 2**self.nqubits
# Initialize a zero matrix of the appropriate dimension.
result = csr_matrix((dim, dim), dtype=_complex_dtype())
for coeff, term in self:
result += coeff * self._apply_operator_on_qubit(term)
return result
[documentos]
def to_qtensor(self, total_nqubits: int | None = None) -> QTensor:
"""Return the Hamiltonian as a ``QTensor`` built from the sparse matrix representation.
Args:
total_nqubits (int, optional): Specify the total number of qubits that this hamiltonian acts on. Defaults to None.
Returns:
QTensor: The QTensor object representation of the Hamiltonian.
Raises:
ValueError: If the total_nqubits provided is lower than the number of qubits effected by the hamiltonian.
"""
nqubits = total_nqubits or self.nqubits
padding = nqubits - self.nqubits
if nqubits < self.nqubits:
raise ValueError(
f"The total number of qubits can't be less than the number of the qubits effected by this hamiltonian ({self.nqubits})"
)
logger.debug("[Hamiltonian] Building QTensor representation over {} qubits", nqubits)
dim = 2 ** (nqubits)
# Initialize a zero matrix of the appropriate dimension.
result = csr_matrix((dim, dim), dtype=_complex_dtype())
for coeff, term in self:
result += coeff * self._apply_operator_on_qubit(term, padding=padding)
return QTensor(result)
[documentos]
def get_static_hamiltonian(self) -> Hamiltonian:
"""Return a Hamiltonian containing only constant coefficients."""
out = Hamiltonian()
for pauli, value in self.elements.items():
aux: Hamiltonian | PauliOperator = pauli[0]
for p in list(pauli)[1:]:
aux *= p
out += aux * value
return out
[documentos]
def draw(self, style: HamiltonianStyle | None = None, filepath: str | None = None) -> None:
"""Render this Hamiltonian as an interaction graph and optionally save it to a file.
Every qubit is drawn as a node whose disc is split into one slice per local field acting
on it (labelled with the Pauli type), and every two-qubit term is drawn as an edge between
the qubits it couples, with a line style per coupling type. Slice and edge colours encode
the coefficient of the corresponding term, as described by the accompanying colour bar.
Expressions acting on three or more qubits are drawn as star-shaped hyperedges joined at their
centroid, and a constant (identity) term is annotated as an energy offset.
If ``filepath`` is given, the resulting figure is saved to disk (the output format is
inferred from the file extension, e.g. ``.png``, ``.pdf``, ``.svg``); otherwise the figure
is shown.
Args:
style (HamiltonianStyle | None, optional): Customization options for the plot appearance.
Defaults to :class:`HamiltonianStyle`.
filepath (str | None, optional): If provided, saves the plot to the specified file path.
Example:
.. code-block:: python
from qilisdk.analog import X, Z
H = X(0) + 2 * Z(1) + 0.5 * Z(0) * Z(1)
H.draw()
"""
logger.debug("[Hamiltonian] Drawing Hamiltonian with style: {} and filepath: {}", style, filepath)
from qilisdk.utils.visualization.hamiltonian_renderers import ( # ruff: ignore[import-outside-top-level]
MatplotlibHamiltonianRenderer,
)
renderer = MatplotlibHamiltonianRenderer(self, style=style or HamiltonianStyle())
renderer.plot()
if filepath:
renderer.save(filepath)
else:
renderer.show()
def __iter__(self) -> Iterator[tuple[complex, list[PauliOperator]]]:
for key, value in self.elements.items():
yield value, list(key)
# ------- Equality & hashing --------
def __eq__(self, other: object) -> bool:
if other == Hamiltonian.ZERO:
return bool(
len(self._elements) == 0
or (len(self._elements) == 1 and (PauliI(0),) in self._elements and self._elements[PauliI(0),] == 0)
)
if isinstance(other, Number):
return bool(
len(self._elements) == 1 and (PauliI(0),) in self._elements and self._elements[PauliI(0),] == other
)
if isinstance(other, PauliOperator):
other = other.to_hamiltonian()
if isinstance(other, QTensor):
other = Hamiltonian.from_qtensor(other)
if not isinstance(other, Hamiltonian):
return False
return dict(self._elements) == dict(other._elements)
def __ne__(self, other: object) -> bool:
return not self.__eq__(other)
def __hash__(self) -> int:
return qili_hash(self._elements)
def __copy__(self) -> Hamiltonian:
return Hamiltonian(elements=self._elements.copy())
# ------- String representation --------
def __repr__(self) -> str:
return str(self)
def __str__(self) -> str:
# Return "0" if there are no terms
if not self._elements:
return "0"
def _format_coeff(c: complex) -> str:
re, im = c.real, c.imag
# 1) Purely real?
if abs(im) < Hamiltonian._EPS:
re_int = np.round(re)
if abs(re - re_int) < Hamiltonian._EPS:
return str(int(re_int)) # e.g. '2' instead of '2.0'
return str(re) # e.g. '2.5'
# 2) Purely imaginary?
if abs(re) < Hamiltonian._EPS:
im_int = np.round(im)
if abs(im - im_int) < Hamiltonian._EPS:
return f"{int(im_int)}j" # e.g. 2 => '2j', -3 => '-3j'
return f"{im}j" # e.g. '2.5j'
# 3) General complex with nonzero real & imag
s = str(c) # e.g. '(3+2j)'
return s
# We want to place the single identity term (I(0),) at the front if it exists
items = list(self.elements.items())
try:
i = next(idx for idx, (key, _) in enumerate(items) if len(key) == 1 and key[0] == (PauliI(0)))
item = items.pop(i)
items.insert(0, item)
except StopIteration:
pass
parts = []
for idx, (operator, coeff) in enumerate(items):
base_str = _format_coeff(coeff)
if idx == 0:
# first term
if len(operator) == 1 and operator[0].name == "I":
coeff_str = base_str
elif base_str == "1":
coeff_str = ""
elif base_str == "-1":
coeff_str = "-"
else:
coeff_str = base_str
elif base_str == "1":
coeff_str = "+"
elif base_str == "-1":
coeff_str = "-"
elif base_str.startswith("-"):
coeff_str = f"- {base_str[1:]}"
else:
coeff_str = f"+ {base_str}"
# Operators string
ops_str = " ".join(str(op) for op in operator if op.name != "I")
if coeff_str and ops_str:
parts.append(f"{coeff_str} {ops_str}")
else:
parts.append(coeff_str + ops_str)
return " ".join(parts)
[documentos]
def get_commuting_partitions(self) -> list[dict[tuple[PauliOperator, ...], complex | Expression | Parameter]]:
"""
Split the Hamiltonian into a list of partitions, each containing commuting terms.
For now this is a greedy algorithm, but a smarter graph-coloring approach could be used later.
Returns:
list[dict[tuple[PauliOperator, ...], complex | Expression | Parameter]]:
A list of dictionaries, each representing a partition of the Hamiltonian containing commuting terms.
"""
logger.debug("[Hamiltonian] Partitioning Hamiltonian into commuting groups")
partitions: list[dict[tuple[PauliOperator, ...], complex | Expression | Parameter]] = []
# Check each term with each partition
for term, coeff in self.elements.items():
placed = False
for partition in partitions:
# Check if the term commutes with all terms in the current partition.
# Both are single Pauli strings, so the parity check is exact and avoids
# building intermediate Hamiltonians.
if all(self._pauli_strings_commute(term, other_term) for other_term in partition):
# If so, add it to this partition
partition[term] = coeff
placed = True
break
# Otherwise create a new partition for this term
if not placed:
partitions.append({term: coeff})
return partitions
@classmethod
[documentos]
def transverse_field(cls, nqubits: int, x_coefficient: Coefficient = 1.0) -> Hamiltonian:
"""
Build a transverse field, :math:`\\sum_i h_i\\, X_i`.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
H = Hamiltonian.transverse_field(nqubits=2, x_coefficient=1.3)
# 1.3 X(0) + 1.3 X(1)
H = Hamiltonian.transverse_field(nqubits=2, x_coefficient=[1.3, 0.7])
# 1.3 X(0) + 0.7 X(1)
Args:
nqubits (int): the number of qubits the Hamiltonian acts on.
x_coefficient (Coefficient, optional): the field strength on every qubit, or a list
holding the strength of each qubit in turn. Defaults to 1.0.
Returns:
Hamiltonian: the transverse-field Hamiltonian.
Raises:
ValueError: if ``nqubits`` is not greater than zero.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits)
logger.debug("[Hamiltonian] Building transverse field over {} qubits", nqubits)
x = _parse_coefficients(x_coefficient, nqubits, "x_coefficient")
return cls({(PauliX(qubit),): x[qubit] for qubit in range(nqubits)})
@classmethod
[documentos]
def longitudinal_field(cls, nqubits: int, z_coefficient: Coefficient = 1.0) -> Hamiltonian:
"""
Build a longitudinal field, :math:`\\sum_i h_i\\, Z_i`.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
H = Hamiltonian.longitudinal_field(nqubits=2, z_coefficient=1.3)
# 1.3 Z(0) + 1.3 Z(1)
H = Hamiltonian.longitudinal_field(nqubits=2, z_coefficient=[1.3, 0.7])
# 1.3 Z(0) + 0.7 Z(1)
Args:
nqubits (int): the number of qubits the Hamiltonian acts on.
z_coefficient (Coefficient, optional): the field strength on every qubit, or a list
holding the strength of each qubit in turn. Defaults to 1.0.
Returns:
Hamiltonian: the longitudinal-field Hamiltonian.
Raises:
ValueError: if ``nqubits`` is not greater than zero.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits)
logger.debug("[Hamiltonian] Building longitudinal field over {} qubits", nqubits)
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
return cls({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
@classmethod
[documentos]
def ising(
cls,
nqubits: int,
zz_coefficient: Coefficient = 1.0,
z_coefficient: Coefficient = 0.0,
) -> Hamiltonian:
"""
Build an all-to-all Ising Hamiltonian, :math:`\\sum_{i<j} J_{ij}\\, Z_i Z_j + \\sum_i h_i\\, Z_i`.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
H = Hamiltonian.ising(nqubits=2, zz_coefficient=2.0)
# 2 Z(0) Z(1)
Args:
nqubits (int): the number of qubits the Hamiltonian acts on. Must be at least 2.
zz_coefficient (Coefficient, optional): the coupling on every pair of qubits, or a list
holding the coupling of each pair in turn, ordered as ``(0, 1), (0, 2), ..., (1, 2),
...``. Defaults to 1.0.
z_coefficient (Coefficient, optional): the longitudinal field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 0.0, which leaves the field out
entirely.
Returns:
Hamiltonian: the Ising Hamiltonian.
Raises:
ValueError: if ``nqubits`` is less than 2.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits, minimum=2, label="Ising")
logger.debug("[Hamiltonian] Building Ising model over {} qubits", nqubits)
pairs = list(combinations(range(nqubits), 2))
zz = _parse_coefficients(zz_coefficient, len(pairs), "zz_coefficient")
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = {
(PauliZ(first), PauliZ(second)): coefficient for (first, second), coefficient in zip(pairs, zz, strict=True)
}
elements.update({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
return cls(elements)
@classmethod
[documentos]
def ising_chain(
cls,
nqubits: int,
zz_coefficient: Coefficient = 1.0,
z_coefficient: Coefficient = 0.0,
periodic: bool = False,
) -> Hamiltonian:
"""
Build a 1D nearest-neighbour Ising chain, :math:`\\sum_i J_i\\, Z_i Z_{i+1} + \\sum_i h_i\\, Z_i`.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
H = Hamiltonian.ising_chain(nqubits=3, zz_coefficient=2.0)
# 2 Z(0) Z(1) + 2 Z(1) Z(2)
# Closing the chain into a ring adds the bond between the two ends.
H = Hamiltonian.ising_chain(nqubits=3, zz_coefficient=2.0, periodic=True)
# 2 Z(0) Z(1) + 2 Z(1) Z(2) + 2 Z(0) Z(2)
Args:
nqubits (int): the number of qubits in the chain. Must be at least 2.
zz_coefficient (Coefficient, optional): the coupling on every bond, or a list holding the
coupling of each bond in turn, ordered along the chain and ending with the
wrap-around bond when ``periodic`` is set. Defaults to 1.0.
z_coefficient (Coefficient, optional): the longitudinal field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 0.0, which leaves the field out
entirely.
periodic (bool, optional): whether to close the chain into a ring by coupling the last
qubit back to the first. Ignored for fewer than 3 qubits, where the ring bond would
duplicate the one bond the chain already has. Defaults to False.
Returns:
Hamiltonian: the Ising chain Hamiltonian.
Raises:
ValueError: if ``nqubits`` is less than 2.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits, minimum=2, label="Ising chain")
logger.debug("[Hamiltonian] Building Ising chain over {} qubits (periodic={})", nqubits, periodic)
bonds = [(qubit, qubit + 1) for qubit in range(nqubits - 1)]
if periodic and nqubits >= _MIN_PERIODIC_EXTENT:
bonds.append((0, nqubits - 1))
zz = _parse_coefficients(zz_coefficient, len(bonds), "zz_coefficient")
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = {
(PauliZ(first), PauliZ(second)): coefficient for (first, second), coefficient in zip(bonds, zz, strict=True)
}
elements.update({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
return cls(elements)
@classmethod
[documentos]
def ising_grid(
cls,
rows: int,
columns: int,
zz_coefficient: Coefficient = 1.0,
z_coefficient: Coefficient = 0.0,
periodic: bool = False,
) -> Hamiltonian:
"""
Build a 2D nearest-neighbour Ising model on a ``rows`` x ``columns`` square lattice.
Qubits are numbered in row-major order, so the qubit at ``(row, column)`` has index
``row * columns + column``, and the Hamiltonian acts on ``rows * columns`` qubits. Each qubit
is coupled to its neighbour to the right and its neighbour below, giving the usual square
lattice.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
# A 2x2 lattice: qubits 0 1 on the top row, 2 3 on the bottom.
H = Hamiltonian.ising_grid(rows=2, columns=2)
# Z(0) Z(1) + Z(0) Z(2) + Z(1) Z(3) + Z(2) Z(3)
Args:
rows (int): the number of lattice rows.
columns (int): the number of lattice columns.
zz_coefficient (Coefficient, optional): the coupling on every bond, or a list holding the
coupling of each bond in turn, ordered by the site each bond starts from in row-major
order and, within a site, giving its horizontal bond before its vertical one.
Defaults to 1.0.
z_coefficient (Coefficient, optional): the longitudinal field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 0.0, which leaves the field out
entirely.
periodic (bool, optional): whether to wrap the lattice into a torus by coupling each edge
back to the opposite one. Wrapping is skipped along any direction shorter than 3
sites, where it would duplicate an existing bond. Defaults to False.
Returns:
Hamiltonian: the Ising grid Hamiltonian.
Raises:
ValueError: if ``rows`` or ``columns`` is not greater than zero.
ValueError: if the lattice holds fewer than 2 qubits.
ValueError: if a list of coefficients does not hold one value per term.
"""
if rows <= 0:
raise ValueError(f"rows must be greater than zero, got {rows}.")
if columns <= 0:
raise ValueError(f"columns must be greater than zero, got {columns}.")
nqubits = rows * columns
_validate_nqubits(nqubits, minimum=2, label="Ising grid")
logger.debug("[Hamiltonian] Building {}x{} Ising grid (periodic={})", rows, columns, periodic)
def index(row: int, column: int) -> int:
return row * columns + column
bonds: list[tuple[int, int]] = []
for row in range(rows):
for column in range(columns):
if column + 1 < columns:
bonds.append((index(row, column), index(row, column + 1)))
elif periodic and columns >= _MIN_PERIODIC_EXTENT:
bonds.append((index(row, 0), index(row, column)))
if row + 1 < rows:
bonds.append((index(row, column), index(row + 1, column)))
elif periodic and rows >= _MIN_PERIODIC_EXTENT:
bonds.append((index(0, column), index(row, column)))
zz = _parse_coefficients(zz_coefficient, len(bonds), "zz_coefficient")
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = {
(PauliZ(first), PauliZ(second)): coefficient for (first, second), coefficient in zip(bonds, zz, strict=True)
}
elements.update({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
return cls(elements)
@classmethod
[documentos]
def transverse_field_ising(
cls,
nqubits: int,
x_coefficient: Coefficient = 1.0,
zz_coefficient: Coefficient = 1.0,
z_coefficient: Coefficient = 0.0,
) -> Hamiltonian:
"""
Build an all-to-all transverse-field Ising Hamiltonian.
.. math::
H = \\sum_{i<j} J_{ij}\\, Z_i Z_j + \\sum_i h^x_i\\, X_i + \\sum_i h^z_i\\, Z_i
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
H = Hamiltonian.transverse_field_ising(nqubits=2, x_coefficient=1.3, zz_coefficient=-2)
# 1.3 X(0) + 1.3 X(1) - 2 Z(0) Z(1)
Args:
nqubits (int): the number of qubits the Hamiltonian acts on. Must be at least 2.
x_coefficient (Coefficient, optional): the transverse field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 1.0.
zz_coefficient (Coefficient, optional): the coupling on every pair of qubits, or a list
holding the coupling of each pair in turn, ordered as ``(0, 1), (0, 2), ..., (1, 2),
...``. Defaults to 1.0.
z_coefficient (Coefficient, optional): the longitudinal field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 0.0, which leaves the field out
entirely.
Returns:
Hamiltonian: the transverse-field Ising Hamiltonian.
Raises:
ValueError: if ``nqubits`` is less than 2.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits, minimum=2, label="Transverse-field Ising")
logger.debug("[Hamiltonian] Building transverse-field Ising model over {} qubits", nqubits)
pairs = list(combinations(range(nqubits), 2))
x = _parse_coefficients(x_coefficient, nqubits, "x_coefficient")
zz = _parse_coefficients(zz_coefficient, len(pairs), "zz_coefficient")
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = {
(PauliX(qubit),): x[qubit] for qubit in range(nqubits)
}
elements.update(
{
(PauliZ(first), PauliZ(second)): coefficient
for (first, second), coefficient in zip(pairs, zz, strict=True)
}
)
elements.update({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
return cls(elements)
@classmethod
[documentos]
def heisenberg(
cls,
nqubits: int,
xx_coefficient: Coefficient = 1.0,
yy_coefficient: Coefficient | None = None,
zz_coefficient: Coefficient | None = None,
z_coefficient: Coefficient = 0.0,
) -> Hamiltonian:
"""
Build an all-to-all Heisenberg Hamiltonian.
.. math::
H = \\sum_{i<j} \\left( J^x_{ij} X_i X_j + J^y_{ij} Y_i Y_j + J^z_{ij} Z_i Z_j \\right)
+ \\sum_i h_i\\, Z_i
Leaving ``yy_coefficient`` and ``zz_coefficient`` at their defaults gives the isotropic XXX
model; setting ``zz_coefficient`` alone gives the XXZ model; setting all three independently
gives the fully anisotropic XYZ model.
Example:
.. code-block:: python
from qilisdk.analog import Hamiltonian
# Isotropic XXX model.
H = Hamiltonian.heisenberg(nqubits=2, xx_coefficient=0.5)
# 0.5 X(0) X(1) + 0.5 Y(0) Y(1) + 0.5 Z(0) Z(1)
# XXZ model with an anisotropic ZZ coupling.
H = Hamiltonian.heisenberg(nqubits=2, xx_coefficient=1.0, zz_coefficient=0.3)
# A disordered field, set qubit by qubit.
H = Hamiltonian.heisenberg(nqubits=3, xx_coefficient=1.0, z_coefficient=[-1.2, 0.4, 2.7])
Args:
nqubits (int): the number of qubits the Hamiltonian acts on. Must be at least 2.
xx_coefficient (Coefficient, optional): the ``XX`` coupling on every pair of qubits, or a
list holding the coupling of each pair in turn, ordered as ``(0, 1), (0, 2), ...,
(1, 2), ...``. Defaults to 1.0.
yy_coefficient (Coefficient, optional): the ``YY`` coupling, given the same way as
``xx_coefficient``. Defaults to None, which reuses ``xx_coefficient``.
zz_coefficient (Coefficient, optional): the ``ZZ`` coupling, given the same way as
``xx_coefficient``. Defaults to None, which reuses ``xx_coefficient``.
z_coefficient (Coefficient, optional): the longitudinal field on every qubit, or a list
holding the field of each qubit in turn. Defaults to 0.0, which leaves the field out
entirely.
Returns:
Hamiltonian: the Heisenberg Hamiltonian.
Raises:
ValueError: if ``nqubits`` is less than 2.
ValueError: if a list of coefficients does not hold one value per term.
"""
_validate_nqubits(nqubits, minimum=2, label="Heisenberg")
if yy_coefficient is None:
yy_coefficient = xx_coefficient
if zz_coefficient is None:
zz_coefficient = xx_coefficient
logger.debug("[Hamiltonian] Building Heisenberg model over {} qubits", nqubits)
pairs = list(combinations(range(nqubits), 2))
xx = _parse_coefficients(xx_coefficient, len(pairs), "xx_coefficient")
yy = _parse_coefficients(yy_coefficient, len(pairs), "yy_coefficient")
zz = _parse_coefficients(zz_coefficient, len(pairs), "zz_coefficient")
z = _parse_coefficients(z_coefficient, nqubits, "z_coefficient")
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = {}
for axis, couplings in (("X", xx), ("Y", yy), ("Z", zz)):
elements.update(
{
(_PAULI_CLASS_BY_NAME[axis](first), _PAULI_CLASS_BY_NAME[axis](second)): coefficient
for (first, second), coefficient in zip(pairs, couplings, strict=True)
}
)
elements.update({(PauliZ(qubit),): z[qubit] for qubit in range(nqubits)})
return cls(elements)
@classmethod
[documentos]
def from_qtensor(cls, tensor: QTensor, tol: float | None = None, prune: float | None = None) -> Hamiltonian:
"""
Expand a qtensor (dense operator) on n qubits into a sum of Pauli strings,
returning a qilisdk.analog.Hamiltonian.
Args:
tol (float): Hermiticity check tolerance. Defaults to global zero tolerance setting.
prune (float): Drop coefficients whose absolute value satisfies ``abs(c) < prune`` to reduce numerical noise. Defaults to global zero tolerance setting.
Returns:
Hamiltonian: Sum_{P in {I,X,Y,Z}^{⊗ n}} c_P * P with c_P = Tr(qt * P) / 2^n
Raises:
ValueError: If the input is not square, not a power-of-two dimension, or not Hermitian w.r.t. `tol`.
"""
if tol is None:
tol = get_settings().atol
if prune is None:
prune = get_settings().atol
A = np.asarray(tensor.dense())
dim = tensor.shape[0]
n = round(np.log2(dim))
logger.debug("[Hamiltonian] Expanding {}-qubit dense operator into Pauli basis", n)
if not tensor.is_hermitian():
raise ValueError("Matrix is not Hermitian within tolerance; cannot form a Hamiltonian.")
# QiliSDK Pauli constructors indexed by qubit id
pauli_for: dict[int, Callable[[int], Hamiltonian]] = {0: I, 1: X, 2: Y, 3: Z}
# Normalization from orthonormality: Tr(P_a P_b) = 2^n δ_ab
norm = 1.0 / (2**n)
# Prebuild per-qubit operator “letters” so we can construct full strings quickly
H = Hamiltonian()
# Full Pauli basis (includes identity on any subset automatically)
for word in product((0, 1, 2, 3), repeat=n):
# Compose the word into a QiliSDK operator acting on the proper qubits
# Example word=(3,1,0) for n=3 -> Z(0)*X(1)*I(2)
op = Hamiltonian()
for q, letter in enumerate(word):
# multiply by the operator on qubit q
op = op * pauli_for[letter](q) if op != 0 else pauli_for[letter](q)
# Convert to dense once; no padding needed because it spans all n qubits
dense_P = op.to_qtensor(n).dense()
# Coefficient c_P = Tr(A P) / 2^n (P is Hermitian)
c = norm * np.trace(A @ dense_P)
# Numerical safety: coefficients should be real for Hermitian A and P
if abs(c.imag) < tol:
c = c.real
if abs(c) > prune:
H += c * op
# Optional: verify round-trip (use a slightly looser atol to tolerate pruning)
if not np.allclose(H.to_qtensor(n).dense(), A, atol=max(10 * prune, 1e-9)):
# If this triggers, consider lowering `prune` or raising `tol`.
raise ValueError("Pauli expansion failed round-trip check; try adjusting tolerances.")
return H
@classmethod
[documentos]
def parse(cls, hamiltonian_str: str) -> Hamiltonian:
logger.debug("[Hamiltonian] Parsing Hamiltonian from string")
hamiltonian_str = hamiltonian_str.strip()
# 1) remove *all* spaces inside any ( … ) group (coefficients or indices)
hamiltonian_str = re.sub(
r"\(\s*([0-9A-Za-z.+\-j\s]+?)\s*\)",
lambda m: "(" + re.sub(r"\s+", "", str(m.group(1))) + ")",
hamiltonian_str,
)
# 2) collapse multiple spaces down to one (outside the parens now)
hamiltonian_str = re.sub(r"\s+", " ", hamiltonian_str)
# 3) ensure a single space between a closing “)” and the next operator token like X(0)/Y(1)/etc.
hamiltonian_str = re.sub(r"\)\s*(?=[XYZI]\()", ") ", hamiltonian_str)
# Special case: "0" => empty Hamiltonian
if hamiltonian_str == "0":
return cls({})
elements: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = defaultdict(complex)
# If there's no initial +/- sign, prepend '+ ' for easier splitting
if not hamiltonian_str.startswith("+") and not hamiltonian_str.startswith("-"):
hamiltonian_str = "+ " + hamiltonian_str
# Replace " - " with " + - " so each term is split on " + "
hamiltonian_str = hamiltonian_str.replace(" - ", " + - ")
# Split on " + "
tokens = hamiltonian_str.split(" + ")
# Remove any empty tokens (can happen if the string started "+ ")
tokens = [t.strip() for t in tokens if t.strip()]
# Regex to match operator tokens like "Z(0)", "X(1)", "I(0)"
operator_pattern = re.compile(r"([XYZI])\((\d+)\)")
def parse_token(token: str) -> tuple[complex, list[PauliOperator]]:
def looks_like_number(text: str) -> bool:
# Make sure it's not empty
if text:
# If the first char is digit, '(', '.', '+', '-', or '0',
# or if 'j' is present, assume it's numeric
first = text[0]
if first.isdigit() or first in {"(", ".", "+", "-"}:
return True
return "j" in text
sign = 1
# Check leading sign
if token.startswith("-"):
sign = -1
token = token[1:].strip()
elif token.startswith("+"):
# optional leading '+'
token = token[1:].strip()
words = token.split()
if not words:
# e.g. just "-" or "+"
# means coefficient = ±1, no operators
return complex(sign), []
# Attempt to parse the first word as a numeric coefficient
maybe_coeff = words[0]
# Decide if 'maybe_coeff' is numeric or an operator
if looks_like_number(maybe_coeff):
# parse as a complex number
coeff_str = maybe_coeff
# If it's e.g. '(2.5+3j)', remove parentheses
if coeff_str.startswith("(") and coeff_str.endswith(")"):
coeff_str = coeff_str[1:-1]
coeff_val = complex(coeff_str) * sign
# consume this word
words = words[1:]
else:
# No explicit coefficient => ±1
coeff_val = complex(sign)
# Now parse the remaining words as operators
ops = []
for w in words:
match = operator_pattern.fullmatch(w)
if not match:
raise ValueError(f"Unrecognized operator format: '{w}'")
name, qubit_str = match.groups()
qubit = int(qubit_str)
op = _PAULI_CLASS_BY_NAME[name](qubit)
ops.append(op)
return coeff_val, ops
for token in tokens:
coeff, op_list = parse_token(token)
if not op_list:
# purely scalar => store as (I(0),)
elements[PauliI(0),] += coeff
else:
# Sort operators by qubit for canonical ordering
op_list.sort(key=lambda op: op.qubit)
elements[tuple(op_list)] += coeff
hamiltonian = cls(elements)
hamiltonian.simplify()
return hamiltonian
[documentos]
def commutator(self, h: Hamiltonian) -> Hamiltonian:
"""compute the commutator of the current hamiltonian with another hamiltonian (h)
Args:
h (Hamiltonian): the second hamiltonian.
Returns:
Hamiltonian: the commutator.
"""
return self * h - h * self
[documentos]
def anticommutator(self, h: Hamiltonian) -> Hamiltonian:
"""compute the anticommutator of the current hamiltonian with another hamiltonian (h)
Args:
h (Hamiltonian): the second hamiltonian.
Returns:
Hamiltonian: the anticommutator.
"""
return self * h + h * self
[documentos]
def commutes_with(self, h: Hamiltonian) -> bool:
"""Check whether this Hamiltonian commutes with another one (``[self, h] == 0``).
This is faster than materialising the full commutator ``self * h - h * self``:
every pair of Pauli strings either commutes or anticommutes, and which one
holds can be decided with a cheap qubit-overlap parity check instead of a
matrix/operator product. Commuting pairs contribute nothing to the commutator
and are skipped entirely, so only the anticommuting pairs (for which
``[P, Q] = 2 P Q``) are ever multiplied out and accumulated.
Args:
h (Hamiltonian): the Hamiltonian to test commutation against.
Returns:
bool: ``True`` if the two Hamiltonians commute, ``False`` otherwise.
"""
residual: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = defaultdict(complex)
for ops1, c1 in self._elements.items():
for ops2, c2 in h._elements.items():
if Hamiltonian._pauli_strings_commute(ops1, ops2):
continue
# Anticommuting strings: [P, Q] = P Q - Q P = 2 P Q.
phase, new_ops = self._multiply_sets(ops1, ops2)
residual[new_ops] += 2 * phase * c1 * c2
return Hamiltonian(residual) == Hamiltonian.ZERO
[documentos]
def vector_norm(self) -> float:
"""
Returns:
float: the vector norm of the hamiltonian.
"""
s = 0
for coeff, _ in self:
s += np.conj(coeff) * coeff
return np.real(np.sqrt(s))
[documentos]
def frobenius_norm(self) -> float:
"""
Returns:
float: the forbenius norm of the hamiltonian.
"""
n = self.nqubits
s = 0
for coeff, _ in self:
s += np.conj(coeff) * coeff
return np.real(np.sqrt(s) * np.sqrt(2**n))
[documentos]
def trace(self) -> Number:
"""
Returns:
float: the trace of the hamiltonian.
"""
n = self.nqubits
d = 2**n
t = self._elements.get((PauliI(0),), 0)
if isinstance(t, Expression):
return t.evaluate({}) * d
return t * d
# ------- Internal multiplication helpers --------
@staticmethod
def _multiply_sets(
set1: tuple[PauliOperator, ...], set2: tuple[PauliOperator, ...]
) -> tuple[complex, tuple[PauliOperator, ...]]:
# Combine all operators into a single list
combined = list(set1) + list(set2)
# Group by qubit
combined.sort(key=lambda op: op.qubit)
sum_dict: dict[int, list[PauliOperator]] = defaultdict(list)
for op in combined:
sum_dict[op.qubit].append(op)
accumulated_phase = complex(1)
final_ops: list[PauliOperator] = []
for qubit_ops in sum_dict.values():
op1 = qubit_ops[0]
phase = complex(1)
# Multiply together all operators on the same qubit
for op2 in qubit_ops[1:]:
aux_phase, op1 = Hamiltonian._multiply_pauli(op1, op2)
phase *= aux_phase
if op1.name != "I":
final_ops.append(op1)
accumulated_phase *= phase
# If everything simplified to identity, we store I(0)
if not final_ops:
final_ops = [PauliI(0)]
# Sort again by qubit (to keep canonical form)
final_ops.sort(key=lambda op: op.qubit)
return accumulated_phase, tuple(final_ops)
@staticmethod
def _pauli_strings_commute(ops1: tuple[PauliOperator, ...], ops2: tuple[PauliOperator, ...]) -> bool:
"""Return whether two Pauli strings commute, using only a qubit-overlap parity check.
Two Pauli strings always either commute or anticommute. They anticommute iff they
disagree (both non-identity and a different Pauli) on an odd number of qubits, since
single-qubit Paulis anticommute exactly when they are distinct and non-identity.
"""
names2 = {op.qubit: op.name for op in ops2 if op.name != "I"}
if not names2:
return True
anticommuting = 0
for op in ops1:
if op.name == "I":
continue
other = names2.get(op.qubit)
if other is not None and other != op.name:
anticommuting += 1
return anticommuting % 2 == 0
@staticmethod
def _multiply_pauli(op1: PauliOperator, op2: PauliOperator) -> tuple[complex, PauliOperator]:
if op1.qubit != op2.qubit:
raise ValueError("Operators must act on the same qubit for multiplication.")
# If either is identity, no phase
if op1.name == "I":
return (1, op2)
if op2.name == "I":
return (1, op1)
# Look up the product in the table
key = (op1.name, op2.name)
result = Hamiltonian._PAULI_PRODUCT_TABLE.get(key)
if result is None:
raise InvalidHamiltonianOperation(f"Multiplying {op1} and {op2} not supported.")
phase, op_cls = result
# By convention, an I operator is always I(0) in this code
if op_cls is PauliI:
return phase, PauliI(0)
# Otherwise, keep the same qubit
return phase, op_cls(op1.qubit)
# ------- Public arithmetic operators --------
def __add__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
out = copy.copy(self)
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
out._add_inplace(other)
return out.simplify()
def __radd__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
return self.__add__(other)
def __sub__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
out = copy.copy(self)
out._sub_inplace(other)
return out.simplify()
def __rsub__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
# (other - self)
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
out = copy.copy(other if isinstance(other, Hamiltonian) else Hamiltonian() + other)
out._sub_inplace(self)
return out.simplify()
def __neg__(self) -> Hamiltonian:
return -1 * self
def __mul__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
out = copy.copy(self)
out._mul_inplace(other)
return out.simplify()
def __rmul__(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> Hamiltonian:
if isinstance(other, Expression) and not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_TERM_MESSAGE)
if isinstance(other, Hamiltonian):
out = copy.copy(other)
out._mul_inplace(self)
return out.simplify()
return self.__mul__(other)
def __truediv__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
out = copy.copy(self)
out._div_inplace(other)
return out.simplify()
def __rtruediv__(self, other: Number | PauliOperator | Hamiltonian) -> Hamiltonian:
# (other / self)
raise InvalidHamiltonianOperation(_DIVISION_BY_OPERATORS_MESSAGE)
__iadd__ = __add__
__isub__ = __sub__
__imul__ = __mul__
__itruediv__ = __truediv__
def _add_inplace(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> None:
if isinstance(other, Hamiltonian):
# If it's empty, do nothing
if not other.elements:
return
# Otherwise, add each term
for key, val in other._elements.items(): # ruff: ignore[private-member-access]
self._elements[key] += val
self._update_parameters(other._parameters) # ruff: ignore[private-member-access]
elif isinstance(other, PauliOperator):
# Just add 1 to that single operator key
self._elements[other,] += 1
elif isinstance(other, (int, float, complex)):
if abs(other) < get_settings().atol:
return
# Add the scalar to (I(0),)
self._elements[PauliI(0),] += other
elif isinstance(other, Expression):
if not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_HAMILTONIAN_MESSAGE)
for parameter in other.free_parameters():
self._add_parameter(parameter.label, parameter)
self._elements[PauliI(0),] += other
else:
raise InvalidHamiltonianOperation(f"Invalid addition between Hamiltonian and {other.__class__.__name__}.")
def _sub_inplace(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> None:
if isinstance(other, Hamiltonian):
for key, val in other._elements.items(): # ruff: ignore[private-member-access]
self._elements[key] -= val
self._update_parameters(other._parameters) # ruff: ignore[private-member-access]
elif isinstance(other, PauliOperator):
self._elements[other,] -= 1
elif isinstance(other, (int, float, complex)):
if abs(other) < get_settings().atol:
return
self._elements[PauliI(0),] -= other
elif isinstance(other, Expression):
if not other.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_HAMILTONIAN_MESSAGE)
for parameter in other.free_parameters():
self._add_parameter(parameter.label, parameter)
self._elements[PauliI(0),] -= other
else:
raise InvalidHamiltonianOperation(
f"Invalid subtraction between Hamiltonian and {other.__class__.__name__}."
)
def _scale_inplace(self, factor: Number | Expression | Parameter) -> None:
"""Multiply every coefficient by a scalar factor.
Raises:
ValueError: if ``factor`` is a symbolic expression that is not fully parameterized.
"""
if isinstance(factor, (int, float, complex)):
# 0 and 1 short-circuits
if abs(factor) < get_settings().atol:
self._elements.clear()
return
if factor == 1:
return
else:
if not factor.is_parameterized():
raise ValueError(_GENERIC_VARIABLE_IN_HAMILTONIAN_MESSAGE)
for parameter in factor.free_parameters():
self._add_parameter(parameter.label, parameter)
for k in self._elements:
self._elements[k] *= factor
@staticmethod
def _identity_coefficient(other: Hamiltonian) -> complex | Expression | Parameter | None:
"""Return the coefficient of ``other`` if it is a scalar times the identity, ``None`` otherwise."""
if len(other._elements) != 1:
return None
((ops, coefficient),) = other._elements.items()
if len(ops) == 1 and ops[0].name == "I" and ops[0].qubit == 0:
return coefficient
return None
def _mul_inplace(self, other: Number | PauliOperator | Hamiltonian | Expression | Parameter) -> None:
if isinstance(other, (int, float, complex, Expression)):
self._scale_inplace(other)
return
if isinstance(other, PauliOperator):
# Convert single PauliOperator -> Hamiltonian with 1 key
other = other.to_hamiltonian()
if not isinstance(other, Hamiltonian):
raise InvalidHamiltonianOperation(
f"Invalid multiplication between Hamiltonian and {other.__class__.__name__}."
)
if not other._elements: # ruff: ignore[private-member-access]
# Multiply by "0" Hamiltonian => 0
self._elements.clear()
return
# A scalar multiple of the identity is just a rescaling
coefficient = self._identity_coefficient(other)
if coefficient is not None:
self._scale_inplace(coefficient)
return
# Otherwise, we do the general multiply
new_dict: dict[tuple[PauliOperator, ...], complex | Expression | Parameter] = defaultdict(complex)
for ops1, c1 in self._elements.items():
for ops2, c2 in other._elements.items(): # ruff: ignore[private-member-access]
phase, new_ops = self._multiply_sets(ops1, ops2)
new_dict[new_ops] += phase * c1 * c2
self._elements = new_dict
self._update_parameters(other._parameters) # ruff: ignore[private-member-access]
def _div_inplace(self, other: Number | PauliOperator | Hamiltonian) -> None:
# Only valid for scalars
if not isinstance(other, (int, float, complex)):
raise InvalidHamiltonianOperation(_DIVISION_BY_OPERATORS_MESSAGE)
if abs(other) < get_settings().atol:
raise ZeroDivisionError("Cannot divide by zero.")
self._mul_inplace(1 / other)
