added custom one-site matrix observables for non-qubit local dimensions, including named position observables

CHANGELOG.md

@@ -12,6 +12,10 @@ releases may include breaking changes.
 
 ### Added
 
+- added custom one-site matrix observables for non-qubit local dimensions
+
+  , including named position observables ([#497]) ([**@linusschulte**])
+
 - added direct MPO Process Tensor construction and temporal entanglement
   ([#508]) ([**@aaronleesander**])
 
@@ -229,6 +233,7 @@ changelogs._
 
 <!-- PR links -->
 
+[#497]: https://github.com/munich-quantum-toolkit/yaqs/pull/497
 [#519]: https://github.com/munich-quantum-toolkit/yaqs/pull/519
 [#518]: https://github.com/munich-quantum-toolkit/yaqs/pull/518
 [#516]: https://github.com/munich-quantum-toolkit/yaqs/pull/516

docs/examples/simulation_parameters.md

@@ -33,13 +33,18 @@ do not import gate classes for standard measurements.
 | `"x"`, `"y"`, `"z"`            | Single-qubit Pauli operators                                   | `Observable("z", sites=0)`                  |
 | `"h"`, `"s"`, `"t"`, `"rx"`, … | Other single-qubit gates from the built-in library             | `Observable("h", sites=0)`                  |
 | `"xx"`, `"yy"`, `"zz"`         | Two-qubit Pauli strings                                        | `Observable("zz", sites=[0, 1])`            |
+| `"position"`                   | Position operator for a supplied local position basis          | `Observable("position", 0, positions=grid)` |
 | `"entropy"`                    | Bipartite entanglement entropy across a cut                    | `Observable("entropy", sites=cut)`          |
 | `"schmidt_spectrum"`           | Schmidt spectrum across a cut                                  | `Observable("schmidt_spectrum", sites=cut)` |
 | bitstring / `"pvm"`            | Projection-valued measurement onto a computational basis state | see {doc}`circuit_observables`              |
 
 For custom unitaries and circuit gates, use {doc}`custom_gates` — those
 workflows still use `GateLibrary` or Qiskit circuits directly.
 
+Named observables that require configuration accept keyword-only factory
+arguments. Missing or unknown arguments raise `TypeError`, so misspelled
+parameters are not silently ignored.
+
 ## Start with a preset
 
 You do **not** need to tune every numerical knob before running a simulation.

docs/examples/trapped_ion.md

@@ -30,7 +30,7 @@ reaches the opposite turning point.
 ```{code-cell} ipython3
 import numpy as np
 
-from mqt.yaqs import Hamiltonian, MPO, State
+from mqt.yaqs import Hamiltonian, MPO, Observable, State
 
 omega = 1.0
 initial_displacement = 1.0
@@ -44,6 +44,7 @@ initial_grid_state /= np.linalg.norm(initial_grid_state)
 
 hamiltonian = Hamiltonian.from_mpo(MPO.trapped_ion(positions, masses=[1.0], omega=omega))
 state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim])
+position_observable = Observable("position", 0, positions=positions)
 ```
 
 ## 2. Noiseless evolution to $T/2$
@@ -52,32 +53,34 @@ state = State(length=1, vector=initial_grid_state, physical_dimensions=[grid_dim
 from mqt.yaqs import AnalogSimParams, Simulator
 
 params = AnalogSimParams(
-    observables=[],
+    observables=[position_observable],
     elapsed_time=half_period,
     dt=half_period / 16,
     max_bond_dim=None,
     svd_threshold=1e-12,
     krylov_tol=1e-12,
     preset="exact",
     get_state=True,
-    sample_timesteps=False,
+    sample_timesteps=True,
 )
 
 result = Simulator(show_progress=False).run(state, hamiltonian, params)
 final_state = result.output_state.vector
-final_x = float(np.sum(positions * np.abs(final_state) ** 2))
+position_expectation = np.real(result.expectation_values[0])
+final_x = float(position_expectation[-1])
 ```
 
-The final $\langle x\rangle$ is close to $-x_0$ but not exact because the
-simulation uses a finite grid and a finite-difference kinetic operator.
+The position observable is a custom one-site matrix on the grid basis. The final
+$\langle x\rangle$ is close to $-x_0$ but not exact because the simulation uses
+a finite grid and a finite-difference kinetic operator.
 
 ```{code-cell} ipython3
 print(f"Initial <x>       = {initial_displacement:.6f}")
 print(f"Final <x> at T/2  = {final_x:.6f}")
 print(f"Continuum target  = {-initial_displacement:.6f}")
 ```
 
-## 3. Wavepacket at $t=0$ and $t=T/2$
+## 3. Wavepacket over time
 
 ```{code-cell} ipython3
 ---
@@ -88,24 +91,27 @@ mystnb:
 ---
 import matplotlib.pyplot as plt
 
-prob_initial = np.abs(initial_grid_state) ** 2
-prob_final = np.abs(final_state) ** 2
-
-fig, axes = plt.subplots(1, 2, figsize=(8, 3.2), layout="constrained", sharey=True)
-axes[0].fill_between(positions, prob_initial, alpha=0.35, color="tab:blue")
-axes[0].plot(positions, prob_initial, color="tab:blue", lw=1.5)
-axes[0].set_title(r"$t = 0$")
-axes[0].set_xlabel(r"$x$")
-axes[0].set_ylabel(r"$|\psi(x)|^2$")
-axes[0].grid(alpha=0.3)
-
-axes[1].fill_between(positions, prob_final, alpha=0.35, color="tab:orange")
-axes[1].plot(positions, prob_final, color="tab:orange", lw=1.5)
-axes[1].set_title(rf"$t = T/2$")
-axes[1].set_xlabel(r"$x$")
-axes[1].grid(alpha=0.3)
-
-fig.suptitle("Harmonic wavepacket reflection on a position grid")
+dense_hamiltonian = hamiltonian.to_matrix()
+eigenvalues, eigenvectors = np.linalg.eigh(dense_hamiltonian)
+coefficients = eigenvectors.conj().T @ initial_grid_state
+phases = np.exp(-1j * eigenvalues[:, None] * params.times[None, :])
+states = eigenvectors @ (coefficients[:, None] * phases)
+probability_density = np.abs(states) ** 2
+
+fig, ax = plt.subplots(figsize=(7.2, 3.6), layout="constrained")
+image = ax.imshow(
+    probability_density,
+    aspect="auto",
+    origin="lower",
+    extent=(params.times[0], params.times[-1], positions[0], positions[-1]),
+    cmap="viridis",
+)
+ax.plot(params.times, position_expectation, color="white", lw=1.4, label=r"$\langle x\rangle$")
+ax.set_xlabel(r"$t$")
+ax.set_ylabel(r"$x$")
+ax.set_title("Position-grid wavepacket density")
+ax.legend(loc="upper right")
+fig.colorbar(image, ax=ax, label=r"$|\psi(x,t)|^2$")
 plt.show()
 ```
 

src/mqt/yaqs/core/data_structures/mps.py

@@ -968,32 +968,34 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
         Returns:
             np.complex128: The computed expectation value (typically, its real part is of interest).
 
+        Raises:
+            ValueError: If the observable is not supported or its matrix shape does not match the target site.
+
         Notes:
             A deep copy of the state is used to prevent modifications to the original MPS.
             Requires :meth:`check_covers_sites` to hold for ``sites``; prefer :meth:`expect` for
             gauge-safe evaluation.
         """
         temp_state = copy.deepcopy(self)
-        if operator.gate.matrix.shape[0] == 2:  # Local observable
-            i = None
-            if isinstance(sites, list):
-                i = sites[0]
-            elif isinstance(sites, int):
-                i = sites
-
-            if isinstance(operator.sites, list):
-                assert operator.sites[0] == i, f"Operator sites mismatch {operator.sites[0]}, {i}"
-            elif isinstance(operator.sites, int):
-                assert operator.sites == i, f"Operator sites mismatch {operator.sites}, {i}"
+        sites_list = [sites] if isinstance(sites, int) else list(sites)
+        operator_sites = [operator.sites] if isinstance(operator.sites, int) else list(operator.sites)
 
-            assert i is not None, f"Invalid type for 'sites': expected int or list[int], got {type(sites).__name__}"
+        if operator.gate.interaction == 1:
+            assert len(sites_list) == 1, f"One-site observable requires one site, got {sites_list}."
+            assert operator_sites == sites_list, f"Operator sites mismatch {operator_sites}, {sites_list}"
+            i = sites_list[0]
             a = temp_state.tensors[i]
-            temp_state.tensors[i] = oe.contract("ab, bcd->acd", operator.gate.matrix, a)
+            local_dim = a.shape[0]
+            matrix = np.asarray(operator.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (local_dim, local_dim):
+                msg = f"Local observable matrix shape {matrix.shape} does not match site {i} dimension {local_dim}."
+                raise ValueError(msg)
+            temp_state.tensors[i] = oe.contract("ab, bcd->acd", matrix, a)
 
-        elif operator.gate.matrix.shape[0] == 4:  # Two-site correlator
+        elif operator.gate.interaction == 2:
             assert isinstance(sites, list)
             assert isinstance(operator.sites, list)
-            i, j = sites
+            i, j = sites_list
 
             assert operator.sites[0] == i, "Observable sites mismatch"
             assert operator.sites[1] == j, "Observable sites mismatch"
@@ -1005,14 +1007,18 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
             b = temp_state.tensors[j]
             d_i, left, _ = a.shape
             d_j, _, right = b.shape
+            matrix = np.asarray(operator.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (d_i * d_j, d_i * d_j):
+                msg = f"Two-site observable matrix shape {matrix.shape} does not match site dimensions {d_i} and {d_j}."
+                raise ValueError(msg)
 
             # 1) merge A,B into theta of shape (l, d_i*d_j, r)
             theta = np.tensordot(a, b, axes=(2, 1))  # (d_i, l, d_j, r)
             theta = theta.transpose(1, 0, 2, 3)  # (l, d_i, d_j, r)
             theta = theta.reshape(left, d_i * d_j, right)  # (l, d_i*d_j, r)
 
             # 2) apply operator on the combined phys index
-            theta = oe.contract("ab, cbd->cad", operator.gate.matrix, theta)  # (l, d_i*d_j, r)
+            theta = oe.contract("ab, cbd->cad", matrix, theta)  # (l, d_i*d_j, r)
             theta = theta.reshape(left, d_i, d_j, right)  # back to (l, d_i, d_j, r)
 
             # 3) split via SVD
@@ -1030,6 +1036,9 @@ def local_expect(self, operator: Observable, sites: int | list[int]) -> np.compl
 
             temp_state.tensors[i] = a_new
             temp_state.tensors[j] = b_new
+        else:
+            msg = "Local observable must be one-site or nearest-neighbor two-site."
+            raise ValueError(msg)
 
         return self.scalar_product(temp_state, sites)
 
@@ -1094,26 +1103,41 @@ def bubble_swaps_backward(state: MPS) -> None:
             for i in reversed(range(state.length - 2)):
                 apply_two_site_nn_inplace(state, i, sw)
 
-        sites = [observable.sites] if isinstance(observable.sites, int) else observable.sites
+        sites = [observable.sites] if isinstance(observable.sites, int) else list(observable.sites)
 
-        if observable.gate.matrix.shape[0] == 2:
+        if observable.gate.interaction == 1:
+            if len(sites) != 1:
+                msg = f"One-site local observable requires one site, got {sites}."
+                raise ValueError(msg)
             site = sites[0]
-            self.tensors[site] = oe.contract("ab, bcd->acd", observable.gate.matrix, self.tensors[site])
+            local_dim = self.tensors[site].shape[0]
+            matrix = np.asarray(observable.gate.matrix, dtype=np.complex128)
+            if matrix.shape != (local_dim, local_dim):
+                msg = f"Local observable matrix shape {matrix.shape} does not match site {site} dimension {local_dim}."
+                raise ValueError(msg)
+            self.tensors[site] = oe.contract("ab, bcd->acd", matrix, self.tensors[site])
             return
 
-        if observable.gate.matrix.shape[0] == 4:
+        if observable.gate.interaction == 2:
+            if len(sites) != 2:
+                msg = f"Two-site local observable requires two sites, got {sites}."
+                raise ValueError(msg)
             i, j = int(sites[0]), int(sites[1])
             length = self.length
+            mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
+            d_i = self.tensors[i].shape[0]
+            d_j = self.tensors[j].shape[0]
+            if mat.shape != (d_i * d_j, d_i * d_j):
+                msg = f"Two-site observable matrix shape {mat.shape} does not match site dimensions {d_i} and {d_j}."
+                raise ValueError(msg)
 
             if length == 2:
                 if i == length - 1 and j == 0:
-                    mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
                     g_merged = permuted_periodic_wrap(mat)
                     apply_two_site_nn_inplace(self, 0, g_merged)
                     return
                 i, j = min(i, j), max(i, j)
             elif (i == length - 1 and j == 0) or (i == 0 and j == length - 1):
-                mat = np.asarray(observable.gate.matrix, dtype=np.complex128)
                 bubble_swaps_forward(self)
                 g_merged = permuted_periodic_wrap(mat)
                 apply_two_site_nn_inplace(self, length - 2, g_merged)
@@ -1124,7 +1148,7 @@ def bubble_swaps_backward(state: MPS) -> None:
                 msg = "Only nearest-neighbor two-site observables are currently implemented."
                 raise ValueError(msg)
 
-            apply_two_site_nn_inplace(self, i, np.asarray(observable.gate.matrix, dtype=np.complex128))
+            apply_two_site_nn_inplace(self, i, mat)
             return
 
         msg = "Local observable must be one-site or nearest-neighbor two-site."

src/mqt/yaqs/core/data_structures/simulation_parameters.py

@@ -23,10 +23,10 @@
 
 import numpy as np
 
-from mqt.yaqs.core.libraries.gate_library import GateLibrary
+from mqt.yaqs.core.libraries.gate_library import BaseGate, GateLibrary
 
 if TYPE_CHECKING:
-    from mqt.yaqs.core.libraries.gate_library import BaseGate
+    from numpy.typing import ArrayLike
 
 SimulationPreset = Literal["fast", "balanced", "accurate", "exact"]
 GateMode = Literal["tdvp", "full-tdvp", "swaps", "mpo"]
@@ -233,31 +233,49 @@ class Observable:
         sites: The site or site indices on which this observable is measured.
     """
 
-    def __init__(self, gate: BaseGate | str, sites: int | list[int] | None = None) -> None:
+    def __init__(
+        self,
+        gate: BaseGate | str | ArrayLike,
+        sites: int | list[int] | None = None,
+        **gate_kwargs: object,
+    ) -> None:
         """Initializes an Observable instance.
 
         Args:
-            gate: The gate that will act as the observable.
+            gate: The gate or one-site local matrix that will act as the observable.
             sites: The qubit or site indices on which this observable is measured.
+            **gate_kwargs: Keyword-only arguments for a named gate or observable factory.
+
+        Raises:
+            TypeError: If factory arguments are missing, unexpected, or supplied for a gate instance or matrix.
         """
         if isinstance(gate, str):
-            if gate == "entropy":
-                gate = GateLibrary.entropy()
-            elif gate == "schmidt_spectrum":
-                gate = GateLibrary.schmidt_spectrum()
-            elif gate == "pvm":
-                gate = GateLibrary.pvm(gate)
+            if gate == "pvm":
+                if gate_kwargs:
+                    msg = "'pvm' does not accept observable parameters."
+                    raise TypeError(msg)
+                resolved_gate = GateLibrary.pvm(gate)
             elif hasattr(GateLibrary, gate):
                 attr = getattr(GateLibrary, gate)
-                try:
-                    gate = attr()
-                except TypeError:
-                    gate = GateLibrary.pvm(gate)
+                resolved_gate = attr(**gate_kwargs)
             else:
-                gate = GateLibrary.pvm(gate)
-        assert hasattr(GateLibrary, gate.name), f"Observable {gate.name} not found in GateLibrary."
-        self.gate = copy.deepcopy(gate)
-        if gate.name != "pvm":
+                if gate_kwargs:
+                    msg = f"Unknown observable {gate!r} does not accept observable parameters."
+                    raise TypeError(msg)
+                resolved_gate = GateLibrary.pvm(gate)
+        elif isinstance(gate, BaseGate):
+            if gate_kwargs:
+                msg = "Observable parameters are only supported for named observables."
+                raise TypeError(msg)
+            resolved_gate = gate
+        else:
+            if gate_kwargs:
+                msg = "Observable parameters are only supported for named observables."
+                raise TypeError(msg)
+            resolved_gate = GateLibrary.local(gate)
+        assert hasattr(GateLibrary, resolved_gate.name), f"Observable {resolved_gate.name} not found in GateLibrary."
+        self.gate: BaseGate = copy.deepcopy(resolved_gate)
+        if resolved_gate.name != "pvm":
             assert sites is not None
             self.sites = sites
             self.gate.set_sites(self.sites)