Coverage for qml_essentials / operations.py: 84%
845 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-08-18 14:58 +0000
1from typing import Callable, List, Optional, Tuple, Union
2from functools import lru_cache
3import string
4import numpy as np
6import jax
7import jax.numpy as jnp
9from qml_essentials.tape import active_tape, recording # noqa: F401 (re-export)
12def _cdtype():
13 """Return the active JAX complex dtype
14 (complex128 if x64 enabled, else complex64).
15 """
16 return jnp.complex128 if jax.config.x64_enabled else jnp.complex64
19@lru_cache(maxsize=256)
20def _einsum_subscript(
21 n: int,
22 k: int,
23 target_axes: Tuple[int, ...],
24) -> str:
25 """Build an ``einsum`` subscript that fuses contraction + axis restore.
27 Args:
28 n: Total rank of the state tensor (number of qubits for statevectors,
29 ``2 * n_qubits`` for density matrices).
30 k: Number of qubits the gate acts on.
31 target_axes: Tuple of k axis indices in the state tensor that the
32 gate contracts against.
34 Returns:
35 ``einsum`` subscript string, e.g. ``"ab,cBd->cad"`` for a 1-qubit
36 gate on wire 1 of a 3-qubit state.
37 """
38 letters = string.ascii_letters
39 # State indices: one letter per axis
40 state_idx = list(letters[:n])
41 # Contracted indices (the ones being replaced by the gate)
42 contracted = [state_idx[ax] for ax in target_axes]
43 # Gate indices: new output indices + contracted input indices
44 new_out = [letters[n + i] for i in range(k)] # fresh letters for output
45 gate_idx = new_out + contracted # gate shape: (out0, out1, ..., in0, in1, ...)
46 # Result indices: replace target axes with new output letters
47 result_idx = list(state_idx)
48 for i, ax in enumerate(target_axes):
49 result_idx[ax] = new_out[i]
50 return "".join(gate_idx) + "," + "".join(state_idx) + "->" + "".join(result_idx)
53def _contract_and_restore(
54 tensor: jnp.ndarray,
55 gate: jnp.ndarray,
56 k: int,
57 target_axes: List[int],
58) -> jnp.ndarray:
59 """Contract gate against target_axes of tensor and restore axis order.
61 The einsum subscript is cached via :func:`_einsum_subscript` so the
62 string construction only happens once per unique
63 ``(total, k, target_axes)`` combination.
65 Args:
66 tensor: Rank-N tensor (e.g. ``(2,)*n`` for states or ``(2,)*2n``
67 for density matrices).
68 gate: Reshaped gate tensor of shape ``(2,)*2k``.
69 k: Number of qubits the gate acts on (= ``len(target_axes)``).
70 target_axes: The k axes of tensor to contract against.
72 Returns:
73 Updated tensor with the same rank as tensor, with the
74 contracted axes restored to their original positions.
75 """
76 subscript = _einsum_subscript(tensor.ndim, k, tuple(target_axes))
77 return jnp.einsum(subscript, gate, tensor)
80class Operation:
81 """Base class for any quantum operation or observable.
83 Further gates should inherit from this class to realise more specific
84 operations. Generally, operations are created by instantiation inside a
85 circuit function passed to :class:`Script`; the instance is
86 automatically appended to the active tape.
88 An ``Operation`` can also serve as an *observable*: its matrix is used to
89 compute expectation values via ``apply_to_state`` / ``apply_to_density``.
91 Attributes:
92 _matrix: Class-level default gate matrix. Subclasses set this to their
93 fixed unitary. Instances may override it via the *matrix* argument
94 to ``__init__``.
95 _num_wires: Expected number of wires for this gate. Subclasses set
96 this to enforce wire count validation. ``None`` means any number
97 of wires is accepted.
98 _param_names: Tuple of attribute names for the gate parameters.
99 Used by :attr:`parameters` and :meth:`__repr__`.
100 """
102 # Subclasses should set this to the gate's unitary / matrix
103 # Whether this is a controlled operation
104 is_controlled = False
105 # Whether this gate is a Clifford gate (normalises the Pauli group
106 is_clifford = False
108 _matrix: jnp.ndarray = None
109 _num_wires: Optional[int] = None
110 _param_names: Tuple[str, ...] = ()
112 def __init__(
113 self,
114 wires: Union[int, List[int]] = 0,
115 matrix: Optional[jnp.ndarray] = None,
116 record: bool = True,
117 name: Optional[str] = None,
118 ) -> None:
119 """Initialise the operation and optionally register it on the active tape.
121 Args:
122 wires: Qubit index or list of qubit indices this operation acts on.
123 matrix: Optional explicit gate matrix. When provided it overrides
124 the class-level ``_matrix`` attribute.
125 record: If ``True`` (default) and a tape is currently recording,
126 append this operation to the tape. Set to ``False`` for
127 auxiliary objects that should not appear in the circuit
128 (e.g. Hamiltonians used only to build time-dependent
129 evolutions).
130 name: Optional explicit name for this operation. When ``None``
131 (default), the class name is used (e.g. ``"RX"``).
133 Raises:
134 ValueError: If ``_num_wires`` is set and the number of wires
135 doesn't match, or if duplicate wires are provided.
136 """
137 self.name = name or self.__class__.__name__
138 self.wires = list(wires) if isinstance(wires, (list, tuple)) else [wires]
140 if self._num_wires is not None and len(self.wires) != self._num_wires:
141 raise ValueError(
142 f"{self.name} expects {self._num_wires} wire(s), "
143 f"got {len(self.wires)}: {self.wires}"
144 )
145 if len(self.wires) != len(set(self.wires)):
146 raise ValueError(f"{self.name} received duplicate wires: {self.wires}")
148 if matrix is not None:
149 self._matrix = matrix
151 # If a tape is currently recording, append ourselves
152 if record:
153 tape = active_tape()
154 if tape is not None:
155 tape.append(self)
157 @property
158 def parameters(self) -> list:
159 """Return the list of numeric parameters for this operation.
161 Uses the declarative ``_param_names`` tuple to collect parameter
162 values in a canonical order. Non-parametrized gates return an
163 empty list.
165 Returns:
166 List of parameter values (floats or JAX arrays).
167 """
168 return [getattr(self, name) for name in self._param_names]
170 def __repr__(self) -> str:
171 """Return a human-readable representation of this operation.
173 Returns:
174 A string like ``"RX(0.5000, wires=[0])"`` or ``"CX(wires=[0, 1])"``.
175 """
176 params = self.parameters
177 if params:
178 param_str = ", ".join(
179 (
180 f"{float(v):.4f}"
181 if isinstance(v, (float, np.floating, jnp.ndarray))
182 else str(v)
183 )
184 for v in params
185 )
186 return f"{self.name}({param_str}, wires={self.wires})"
187 return f"{self.name}(wires={self.wires})"
189 @property
190 def matrix(self) -> jnp.ndarray:
191 """Return the base matrix of this operation (before lifting).
193 Returns:
194 The gate matrix as a JAX array.
196 Raises:
197 NotImplementedError: If the subclass has not defined ``_matrix``.
198 """
199 if self._matrix is None:
200 raise NotImplementedError(
201 f"{self.__class__.__name__} does not define a matrix."
202 )
203 return self._matrix
205 def decompose(self) -> List["Operation"]:
206 """Decompose this operation into a list of more primitive operations.
208 The returned operations are created with ``record=False`` so the caller
209 controls where they are placed. Reused e.g. by
210 :meth:`~qml_essentials.pauli.PauliCircuit.get_clifford_pauli_gates` to
211 express composite gates in terms of Clifford + Pauli-rotation primitives.
213 Returns:
214 List of :class:`Operation` instances equivalent to this gate.
216 Raises:
217 NotImplementedError: If the gate has no decomposition (it is itself
218 primitive).
219 """
220 raise NotImplementedError(
221 f"{self.__class__.__name__} does not define a decomposition."
222 )
224 @property
225 def wires(self) -> List[int]:
226 """Qubit indices this operation acts on.
228 Returns:
229 List of integer qubit indices.
230 """
231 return self._wires
233 @wires.setter
234 def wires(self, wires: Union[int, List[int]]) -> None:
235 """Set the qubit indices for this operation.
237 Args:
238 wires: A single qubit index or a list of qubit indices.
239 """
240 if isinstance(wires, (list, tuple)):
241 self._wires = list(wires)
242 else:
243 self._wires = [wires]
245 def _update_tape_operation(self, op: "Operation") -> None:
246 """
247 If ``self`` is already on the active tape (the typical case when
248 chaining ``Gate(...).dagger()``), it is replaced by the daggered
249 operation so that only U\\dagger appears on the tape —
250 not both U and ``U\\dagger``.
251 Note that this should only be called immediately after the tape is updated.s
253 Args:
254 op (Operation): New replaced operation on the tape
255 """
256 # If self was recorded on the tape, replace it with the daggered op.
257 tape = active_tape()
258 if tape is not None:
259 if tape and tape[-1] is self:
260 tape[-1] = op
261 else:
262 tape.append(op)
264 def dagger(self) -> "Operation":
265 """Return a new operation, the conjugate transpose (``U\\dagger``)
266 Usage inside a circuit function::
268 RX(0.5, wires=0).dagger()
270 Returns:
271 A new :class:`Operation` with matrix ``U\\dagger`` acting on the same wires.
272 """
273 mat = jnp.conj(self._matrix).T
274 op = Operation(wires=self.wires, matrix=mat, record=False)
276 self._update_tape_operation(op)
278 return op
280 def power(self, power) -> "Operation":
281 """Return a new operation, the power (``U^power``)
282 Usage inside a circuit function::
284 PauliX(wires=0).power(2)
286 Returns:
287 A new :class:`Operation` with matrix ``U\\dagger`` acting on the same wires.
288 """
289 # TODO: support fractional powers
290 mat = jnp.linalg.matrix_power(self._matrix, power)
291 op = Operation(wires=self.wires, matrix=mat, record=False)
293 self._update_tape_operation(op)
295 return op
297 def __mul__(self, other: Union[float, "Operation"]) -> "Operation":
298 """Return a new operation, the product between U and a scalar (``U*x``)
299 or the composition of two operations.
300 Usage inside a circuit function::
302 PauliX(wires=0) * x
303 PauliX(wires=0) * PauliZ(wires=0)
305 Returns:
306 A new :class:`Operation` with matrix ``U*x`` acting on the same wires,
307 or the composed matrix acting on the appropriate wires.
308 """
309 if isinstance(other, Operation):
310 return self.__matmul__(other)
312 mat = other * self._matrix
313 op = Operation(wires=self.wires, matrix=mat, record=False)
315 self._update_tape_operation(op)
317 return op
319 # Also overwrite * for right operands
320 __rmul__ = __mul__
322 def __add__(self, other: "Operation") -> "Operation":
323 """Element-wise addition of two operations on the same wires.
325 Returns:
326 A new :class:`Operation` whose matrix is the sum of both matrices.
328 Raises:
329 ValueError: If the wire sets differ.
330 """
331 if sorted(self.wires) != sorted(other.wires):
332 raise ValueError(
333 f"Can only add operations acting on the same set of wires, "
334 f"got {self.wires} and {other.wires}"
335 )
337 op = Operation(
338 wires=self.wires,
339 matrix=self.matrix + other.matrix,
340 record=False,
341 )
342 return op
344 def prod(self, *ops: "Operation") -> "Operation":
345 """Construct the generalized product (tensor or matrix)
346 of this operation with others.
348 The resulting operation acts on the union of all wire sets.
349 If the wire sets are disjoint, this is a Kronecker product.
350 If the wire sets overlap, the corresponding matrices are multiplied.
352 Usage::
354 res = op1.prod(op2, op3)
355 # or
356 res = Operation.prod(op1, op2, op3)
358 Args:
359 *ops: Variable number of :class:`Operation` instances.
361 Returns:
362 A new :class:`Operation` representing the composed operation.
363 """
364 if not ops:
365 return self
367 all_ops = (self,) + ops
368 all_wires = []
369 for op in all_ops:
370 for w in op.wires:
371 if w not in all_wires:
372 all_wires.append(w)
374 n = len(all_wires)
376 mat = _embed_matrix(all_ops[0].matrix, all_ops[0].wires, all_wires, n)
377 for op in all_ops[1:]:
378 mat_other = _embed_matrix(op.matrix, op.wires, all_wires, n)
379 mat = mat @ mat_other
381 op_names = "*".join(op.name for op in all_ops)
382 return Operation(
383 wires=all_wires, matrix=mat, name=f"Prod({op_names})", record=False
384 )
386 def __matmul__(self, other: "Operation") -> "Operation":
387 """Tensor (Kronecker) product or matrix product of two operations.
389 The resulting operation acts on the union of both wire sets.
390 If the wire sets are disjoint, this is a Kronecker product.
391 If the wire sets overlap, the corresponding matrices are multiplied.
393 Returns:
394 A new :class:`Operation` whose matrix represents the composed
395 operation on the unified wire set.
396 """
397 if not isinstance(other, Operation):
398 return NotImplemented
400 return self.prod(other)
402 def lifted_matrix(self, n_qubits: int) -> jnp.ndarray:
403 """Return the full ``2**n x 2**n`` matrix embedding this gate.
405 Embeds the ``k``-qubit gate matrix into the ``n``-qubit Hilbert space
406 by applying it to the identity matrix via :meth:`apply_to_state`.
407 This is useful for computing ``Tr(O·\\rho )`` directly without vmap.
409 Args:
410 n_qubits: Total number of qubits in the circuit.
412 Returns:
413 The ``(2**n, 2**n)`` matrix of this operation in the full space.
414 """
415 dim = 2**n_qubits
416 # Apply the gate to each basis vector (column of identity)
417 return jax.vmap(lambda col: self.apply_to_state(col, n_qubits))(
418 jnp.eye(dim, dtype=_cdtype())
419 ).T
421 def apply_to_state(self, state: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
422 """Apply this gate to a statevector via tensor contraction.
424 The statevector (shape ``(2**n,)``) is reshaped into a rank-n tensor
425 of shape ``(2,)*n``. The gate (shape ``(2**k, 2**k)``) is reshaped to
426 ``(2,)*2k`` and contracted against the k target wire axes.
428 Memory footprint is O(2**n) and the operation supports arbitrary k.
429 The implementation is fully differentiable through JAX.
431 Args:
432 state: Statevector of shape ``(2**n_qubits,)``.
433 n_qubits: Total number of qubits in the circuit.
435 Returns:
436 Updated statevector of shape ``(2**n_qubits,)``.
437 """
438 k = len(self.wires)
439 gate_tensor = self.matrix.reshape((2,) * 2 * k)
440 psi = state.reshape((2,) * n_qubits)
441 psi_out = _contract_and_restore(psi, gate_tensor, k, self.wires)
442 return psi_out.reshape(2**n_qubits)
444 def apply_to_state_tensor(self, psi: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
445 """Apply this gate to a statevector already in tensor form.
447 Like :meth:`apply_to_state` but expects the state in rank-n tensor
448 form ``(2,)*n`` and returns the result in the same form. This avoids
449 the ``reshape`` calls at the per-gate level when the simulation loop
450 keeps the state in tensor form throughout.
452 Args:
453 psi: Statevector tensor of shape ``(2,)*n_qubits``.
454 n_qubits: Total number of qubits in the circuit.
456 Returns:
457 Updated statevector tensor of shape ``(2,)*n_qubits``.
458 """
459 k = len(self.wires)
460 gate_tensor = self._gate_tensor(k)
461 return _contract_and_restore(psi, gate_tensor, k, self.wires)
463 def _gate_tensor(self, k: int) -> jnp.ndarray:
464 """Return the gate matrix reshaped to ``(2,)*2k`` tensor form.
466 The result is cached on the instance so repeated calls (e.g. from
467 density-matrix simulation which applies U and U*) avoid redundant
468 reshape dispatch.
470 Args:
471 k: Number of qubits the gate acts on.
473 Returns:
474 Gate matrix as a rank-2k tensor of shape ``(2,)*2k``.
475 """
476 cached = getattr(self, "_cached_gate_tensor", None)
477 if cached is not None:
478 return cached
479 gt = self.matrix.reshape((2,) * 2 * k)
480 # Only cache for non-parametrized gates (whose matrix is a class attr)
481 if self._matrix is self.__class__._matrix:
482 object.__setattr__(self, "_cached_gate_tensor", gt)
483 return gt
485 def apply_to_density(self, rho: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
486 """Apply this gate to a density matrix via \\rho -> U\\rho U\\dagger.
488 The density matrix (shape ``(2**n, 2**n)``) is treated as a rank-*2n*
489 tensor with n "ket" axes (0..n-1) and n "bra" axes (n..2n-1).
490 U acts on the ket half; U* acts on the bra half. Both contractions
491 use the shared :func:`_contract_and_restore` helper, keeping the
492 operation allocation-free with respect to building full unitaries.
494 Args:
495 rho: Density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
496 n_qubits: Total number of qubits in the circuit.
498 Returns:
499 Updated density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
500 """
501 k = len(self.wires)
502 U = self._gate_tensor(k)
503 U_conj = jnp.conj(U)
505 rho_t = rho.reshape((2,) * 2 * n_qubits)
507 # Apply U to ket axes, U\\dagger to bra axes
508 rho_t = _contract_and_restore(rho_t, U, k, self.wires)
509 bra_wires = [w + n_qubits for w in self.wires]
510 rho_t = _contract_and_restore(rho_t, U_conj, k, bra_wires)
512 return rho_t.reshape(2**n_qubits, 2**n_qubits)
515class Hermitian(Operation):
516 """A generic Hermitian observable or gate defined by an arbitrary matrix.
518 Example:
519 >>> obs = Hermitian(matrix=my_matrix, wires=0)
520 """
522 def __init__(
523 self,
524 matrix: jnp.ndarray,
525 wires: Union[int, List[int]] = 0,
526 record: bool = True,
527 ) -> None:
528 """Initialise a Hermitian operator.
530 Args:
531 matrix: The Hermitian matrix defining this operator.
532 wires: Qubit index or list of qubit indices this operator acts on.
533 record: If ``True`` (default), record on the active tape. Set to
534 ``False`` when using the Hermitian purely as a Hamiltonian
535 component (e.g. for time-dependent evolution).
536 """
537 super().__init__(
538 wires=wires,
539 matrix=jnp.asarray(matrix, dtype=_cdtype()),
540 record=record,
541 )
543 def __rmul__(self, coeff_fn: Callable) -> "ParametrizedHamiltonian":
544 """Support ``coeff_fn * Hermitian`` -> :class:`ParametrizedHamiltonian`.
546 Args:
547 coeff_fn (Callable): A callable ``(params, t) -> scalar`` giving the
548 time-dependent coefficient.
550 Returns:
551 ParametrizedHamiltonian: A :class:`ParametrizedHamiltonian` pairing
552 *coeff_fn* with this operator's matrix and wires.
554 Raises:
555 TypeError: If *coeff_fn* is not callable.
556 """
557 if not callable(coeff_fn):
558 raise TypeError(
559 f"Left operand of `* Hermitian` must be callable, got {type(coeff_fn)}"
560 )
561 return ParametrizedHamiltonian(terms=[(coeff_fn, self.matrix, self.wires)])
563 def evolve(self, name: Optional[str] = None, **odeint_kwargs) -> Callable:
564 """Return a gate factory for static evolution ``U = exp(-i t H)``.
566 Thin delegator to :meth:`qml_essentials.evolution.Evolution.evolve`.
568 Args:
569 name: Optional name for the produced :class:`Operation`.
570 **odeint_kwargs: Unused for static evolution (accepted for a
571 uniform signature with :meth:`ParametrizedHamiltonian.evolve`).
573 Returns:
574 A callable gate factory ``(t, wires=0) -> Operation``.
575 """
576 from qml_essentials.evolution import Evolution # deferred: circular import
578 return Evolution.evolve(self, name=name, **odeint_kwargs)
581class ParametrizedHamiltonian:
582 """A time-dependent Hamiltonian as a sum of ``coeff * Hermitian`` terms.
584 Mathematically::
586 H(t) = \\sum_i f_i(params_i, t) * H_i
588 Construction is always done from an explicit list of
589 ``(coeff_fn, H_mat, wires)`` triples passed as ``terms``. The
590 common single-term shorthand is the operator form
591 ``coeff_fn * Hermitian(matrix, wires)`` (see
592 :meth:`Hermitian.__rmul__`), which returns a one-term instance.
593 Multi-term Hamiltonians are composed with ``+`` between
594 :class:`ParametrizedHamiltonian` instances::
596 H1 = coeff_x * Hermitian(X, wires=0)
597 H2 = coeff_y * Hermitian(Y, wires=0)
598 H_td = H1 + H2
600 # evolve under the composite Hamiltonian; coeff_args is a list of
601 # parameter sets, one per term, in the order the terms were added:
602 H_td.evolve()([px, py], T=1.0)
604 Attributes:
605 coeff_fns: Tuple of callables ``(params, t) -> scalar``, one per term.
606 H_mats: Tuple of static Hermitian matrices, one per term.
607 wires: Wires this Hamiltonian acts on (union across all terms; for
608 now all terms are required to share the same wire set).
609 """
611 def __init__(
612 self,
613 terms: List[Tuple[Callable, jnp.ndarray, Union[int, List[int]]]],
614 ) -> None:
615 """Build a (possibly multi-term) parametrized Hamiltonian.
617 Args:
618 terms: List of ``(coeff_fn, H_mat, wires)`` triples. Use the
619 ``coeff_fn * Hermitian(...)`` shorthand to build a
620 one-term instance; combine instances with ``+`` to add
621 terms.
623 Raises:
624 ValueError: If the term list is empty, or if terms act on
625 differing wire sets (multi-wire broadcasting is
626 deferred — see :mod:`jaqsi`), or if term matrices have
627 incompatible shapes.
628 """
629 if len(terms) == 0:
630 raise ValueError("ParametrizedHamiltonian needs at least one term.")
632 # Normalise wires (single int -> [int]) and validate consistency.
633 def _wlist(w):
634 return [w] if isinstance(w, int) else list(w)
636 first_wires = _wlist(terms[0][2])
637 for _, _, w in terms[1:]:
638 if _wlist(w) != first_wires:
639 raise ValueError(
640 "All terms of a ParametrizedHamiltonian must currently "
641 "act on the same wires; got "
642 f"{_wlist(w)} vs. {first_wires}. "
643 "Multi-wire broadcasting across terms is not yet supported."
644 )
646 # Validate matrix shape compatibility across terms.
647 first_dim = jnp.asarray(terms[0][1]).shape
648 for _, H, _ in terms[1:]:
649 if jnp.asarray(H).shape != first_dim:
650 raise ValueError(
651 f"All term matrices must have the same shape; got "
652 f"{jnp.asarray(H).shape} vs. {first_dim}."
653 )
655 self._terms: Tuple[Tuple[Callable, jnp.ndarray, List[int]], ...] = tuple(
656 (fn, jnp.asarray(H, dtype=_cdtype()), _wlist(w)) for fn, H, w in terms
657 )
658 self.wires: List[int] = list(first_wires)
660 # --- term accessors -------------------------------------------------
662 @property
663 def coeff_fns(self) -> Tuple[Callable, ...]:
664 """Tuple of coefficient functions, one per term."""
665 return tuple(fn for fn, _, _ in self._terms)
667 @property
668 def H_mats(self) -> Tuple[jnp.ndarray, ...]:
669 """Tuple of Hermitian matrices, one per term."""
670 return tuple(H for _, H, _ in self._terms)
672 @property
673 def n_terms(self) -> int:
674 """Number of terms in the Hamiltonian."""
675 return len(self._terms)
677 # --- composition ---------------------------------------------------
679 def __add__(self, other: "ParametrizedHamiltonian") -> "ParametrizedHamiltonian":
680 """Concatenate term lists: ``H = H1 + H2``."""
681 if not isinstance(other, ParametrizedHamiltonian):
682 return NotImplemented
683 return ParametrizedHamiltonian(terms=list(self._terms) + list(other._terms))
685 def __neg__(self) -> "ParametrizedHamiltonian":
686 """Negate every coefficient: ``-H`` = sum of ``(-f_i) * H_i``."""
687 new_terms = [
688 ((lambda f: lambda p, t: -f(p, t))(fn), H, w) for fn, H, w in self._terms
689 ]
690 return ParametrizedHamiltonian(terms=new_terms)
692 def __sub__(self, other: "ParametrizedHamiltonian") -> "ParametrizedHamiltonian":
693 if not isinstance(other, ParametrizedHamiltonian):
694 return NotImplemented
695 return self + (-other)
697 # --- evolution -----------------------------------------------------
699 def evolve(self, name: Optional[str] = None, **odeint_kwargs) -> Callable:
700 """Return a gate factory for time-dependent evolution.
702 Solves ``dU/dt = -i [sum_i f_i(p_i, t) H_i] U``. Thin delegator to
703 :meth:`qml_essentials.evolution.Evolution.evolve`.
705 Args:
706 name: Optional name for the produced :class:`Operation`.
707 **odeint_kwargs: Solver options forwarded to ``Evolution.evolve``
708 (``atol``, ``rtol``, ``max_steps``, ``throw``, ``solver``,
709 ``magnus_steps``).
711 Returns:
712 A callable gate factory ``(coeff_args, T) -> Operation``.
713 """
714 from qml_essentials.evolution import Evolution # deferred: circular import
716 return Evolution.evolve(self, name=name, **odeint_kwargs)
719class Id(Operation):
720 """Identity gate.
722 Supports an arbitrary number of wires. When more than one wire is
723 given the matrix is the ``2**k x 2**k`` identity (where *k* is the
724 number of wires).
725 """
727 _matrix = jnp.eye(2, dtype=_cdtype())
728 _num_wires = None # accept any number of wires
729 is_clifford = True
731 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
732 """Initialise an identity gate.
734 Args:
735 wires: Qubit index or list of qubit indices this gate acts on.
736 When multiple wires are given the matrix is automatically
737 expanded to the matching ``2**k × 2**k`` identity.
738 """
739 w = list(wires) if isinstance(wires, (list, tuple)) else [wires]
740 k = len(w)
741 if k > 1:
742 kwargs["matrix"] = jnp.eye(2**k, dtype=_cdtype())
743 super().__init__(wires=wires, **kwargs)
746class PauliX(Operation):
747 """Pauli-X gate / observable (bit-flip, \\sigma_x)."""
749 _matrix = jnp.array([[0, 1], [1, 0]], dtype=_cdtype())
750 _num_wires = 1
751 is_clifford = True
753 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
754 """Initialise a Pauli-X gate.
756 Args:
757 wires: Qubit index or list of qubit indices this gate acts on.
758 """
759 super().__init__(wires=wires, **kwargs)
762class PauliY(Operation):
763 """Pauli-Y gate / observable (\\sigma_y)."""
765 _matrix = jnp.array([[0, -1j], [1j, 0]], dtype=_cdtype())
766 _num_wires = 1
767 is_clifford = True
769 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
770 """Initialise a Pauli-Y gate.
772 Args:
773 wires: Qubit index or list of qubit indices this gate acts on.
774 """
775 super().__init__(wires=wires, **kwargs)
778class PauliZ(Operation):
779 """Pauli-Z gate / observable (phase-flip, \\sigma_z)."""
781 _matrix = jnp.array([[1, 0], [0, -1]], dtype=_cdtype())
782 _num_wires = 1
783 is_clifford = True
785 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
786 """Initialise a Pauli-Z gate.
788 Args:
789 wires: Qubit index or list of qubit indices this gate acts on.
790 """
791 super().__init__(wires=wires, **kwargs)
794class H(Operation):
795 """Hadamard gate."""
797 _matrix = jnp.array([[1, 1], [1, -1]], dtype=_cdtype()) / jnp.sqrt(2)
798 _num_wires = 1
799 is_clifford = True
801 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
802 """Initialise a Hadamard gate.
804 Args:
805 wires: Qubit index or list of qubit indices this gate acts on.
806 """
807 super().__init__(wires=wires, **kwargs)
810class S(Operation):
811 """S (phase) gate — a Clifford gate equal to \\sqrt Z.
813 .. math::
814 S = \\begin{pmatrix}1 & 0\\ 0 & i\\end{pmatrix}
815 """
817 _matrix = jnp.array([[1, 0], [0, 1j]], dtype=_cdtype())
818 _num_wires = 1
819 is_clifford = True
821 def __init__(self, wires: Union[int, List[int]] = 0) -> None:
822 """Initialise an S gate.
824 Args:
825 wires: Qubit index or list of qubit indices this gate acts on.
826 """
827 super().__init__(wires=wires)
830class SWAP(Operation):
831 """SWAP gate."""
833 _matrix = jnp.array(
834 [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]], dtype=_cdtype()
835 )
836 _num_wires = 2
837 is_clifford = True
839 def __init__(self, wires: Union[int, List[int]] = 0, **kwargs) -> None:
840 """Initialise a SWAP gate.
842 Args:
843 wires: Qubit index or list of qubit indices this gate acts on.
844 """
845 super().__init__(wires=wires, **kwargs)
848class RandomUnitary(Operation):
849 """Creates a random hermitian matrix and applies it as a gate."""
851 def __init__(
852 self,
853 wires: Union[int, List[int]],
854 key: jax.random.PRNGKey,
855 scale: float = 1.0,
856 record: bool = True,
857 ) -> None:
858 """Initialise a random unitary gate.
860 Args:
861 wires (Union[int, List[int]]): Qubit index or list of qubit indices
862 this gate acts on.
863 key (jax.random.PRNGKey): PRNGKey for randomization.
864 scale (float): Scale of the random unitary (default: 1.0).
865 record (bool): Whether to record this gate on the active tape.
866 """
867 dim = 2 ** len(wires)
868 key_a, key_b = jax.random.split(key)
870 A = (
871 jax.random.normal(key=key_a, shape=(dim, dim))
872 + 1j * jax.random.normal(key=key_b, shape=(dim, dim))
873 ).astype(_cdtype())
874 H = (A + A.conj().T) / 2.0
876 H *= scale / jnp.linalg.norm(H, ord="fro")
878 super().__init__(wires, matrix=H, record=record)
881class DiagonalQubitUnitary(Operation):
882 """A diagonal unitary gate specified by its diagonal entries.
884 Implements ``U = diag(d_0, d_1, ..., d_{2^k-1})`` where each ``d_i`` lies
885 on the unit circle. This is the natural gate for data-encoding
886 Hamiltonians of the form ``S(x) = exp(-i H x)`` where *H* is diagonal in
887 the computational basis (see Peters et al., arXiv:2209.05523).
889 The Golomb encoding strategy uses this gate with diagonal entries
890 ``exp(-i * golomb_marks * x)`` to achieve a maximally non-degenerate
891 Fourier spectrum.
893 Args:
894 diag: 1-D array of ``2**k`` complex values on the unit circle.
895 wires: Qubit indices this gate acts on (s.t. ``2**len(wires) == len(diag)``).
896 **kwargs: Forwarded to :class:`Operation`.
897 """
899 # Do NOT list "diag" in _param_names — the array is not a scalar
900 # parameter and would break drawing helpers that call float(p).
901 _param_names = ()
903 def __init__(
904 self,
905 diag: jnp.ndarray,
906 wires: Union[int, List[int]] = 0,
907 generator: Optional[jnp.ndarray] = None,
908 scale: Optional[float] = None,
909 **kwargs,
910 ) -> None:
911 self.diag = diag
912 # Optional real data-encoding generator: ``diag = exp(-i * generator *
913 # scale)`` with a real diagonal Hamiltonian ``generator`` and real
914 # scalar ``scale``. When present, :meth:`decompose` expands the gate
915 # into commuting Pauli-Z rotations; the complex ``diag`` alone is
916 # insufficient because its phase wraps modulo ``2 pi``.
917 self._generator = generator
918 self._scale = scale
919 wires_list = list(wires) if isinstance(wires, (list, tuple)) else [wires]
920 expected_dim = 2 ** len(wires_list)
921 if diag.shape != (expected_dim,):
922 raise ValueError(
923 f"DiagonalQubitUnitary expects {expected_dim} diagonal entries "
924 f"for {len(wires_list)} wire(s), got shape {diag.shape}"
925 )
926 mat = jnp.diag(diag)
927 # Use a descriptive name for drawing
928 kwargs.setdefault("name", "DiagU")
929 super().__init__(wires=wires, matrix=mat, **kwargs)
931 def decompose(self) -> List["Operation"]:
932 r"""Expand a real-generator diagonal encoding into Pauli-Z rotations.
934 For ``diag = exp(-i H x)`` with a real diagonal Hamiltonian
935 ``H = diag(self._generator)`` and real scalar ``x = self._scale``, the
936 Walsh-Hadamard transform writes ``H = \sum_P \alpha_P P`` over commuting
937 Pauli-Z strings ``P``. Because the strings commute,
939 .. math::
940 e^{-i H x} = \prod_P e^{-i \alpha_P x P}
941 = \prod_P \mathrm{PauliRot}(2 x \alpha_P, P),
943 up to the global phase from the identity term (dropped). Zero-weight
944 strings are omitted.
946 Returns:
947 List of :class:`PauliRot` gates (``record=False``) whose ordered
948 product equals the diagonal gate up to a global phase.
950 Raises:
951 NotImplementedError: If the gate carries no real generator (a
952 generic diagonal unitary has no Pauli-rotation decomposition).
953 """
954 if self._generator is None:
955 return super().decompose()
957 k = len(self.wires)
958 dim = 2**k
959 marks = np.asarray(self._generator, dtype=float).reshape(-1)
960 # Sign vector per qubit position: +1/-1 for basis bit 0/1. Position i
961 # (i = 0 is the most-significant bit of the diagonal index) acts on
962 # ``self.wires[i]`` -- the same order PauliRot uses for its word.
963 signs = np.array(
964 [[1 - 2 * ((j >> (k - 1 - i)) & 1) for j in range(dim)] for i in range(k)]
965 )
967 ops: List["Operation"] = []
968 tol = 1e-12
969 for mask in range(1, dim): # skip mask 0 (identity -> global phase)
970 chi = np.ones(dim)
971 for i in range(k):
972 if (mask >> i) & 1:
973 chi = chi * signs[i]
974 alpha = float(marks @ chi) / dim
975 if abs(alpha) < tol:
976 continue
977 word = "".join("Z" if (mask >> i) & 1 else "I" for i in range(k))
978 theta = 2.0 * alpha * self._scale
979 ops.append(PauliRot(theta, word, wires=self.wires, record=False))
980 return ops
982 def apply_to_state(self, state: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
983 """Apply diagonal gate via element-wise multiplication.
985 For a diagonal unitary, the full ``2^n``-dimensional diagonal is
986 constructed by appropriate Kronecker-product embedding and the gate
987 is applied as an element-wise product, which is significantly cheaper
988 than generic matrix contraction for large qubit counts.
990 Args:
991 state: Statevector of shape ``(2**n_qubits,)``.
992 n_qubits: Total number of qubits in the circuit.
994 Returns:
995 Updated statevector of shape ``(2**n_qubits,)``.
996 """
997 k = len(self.wires)
998 if k == n_qubits and self.wires == list(range(n_qubits)):
999 # Gate acts on all qubits in order — direct element-wise multiply
1000 return state * self.diag
1001 # Fall back to general tensor contraction for arbitrary wire subsets
1002 return super().apply_to_state(state, n_qubits)
1004 def apply_to_density(self, rho: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1005 """Apply diagonal gate to density matrix: rho -> U rho U†.
1007 For diagonal U the transformation is
1008 ``rho_ij -> d_i * conj(d_j) * rho_ij``.
1010 Args:
1011 rho: Density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
1012 n_qubits: Total number of qubits in the circuit.
1014 Returns:
1015 Updated density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
1016 """
1017 k = len(self.wires)
1018 if k == n_qubits and self.wires == list(range(n_qubits)):
1019 d = self.diag
1020 return d[:, None] * jnp.conj(d)[None, :] * rho
1021 return super().apply_to_density(rho, n_qubits)
1024class Barrier(Operation):
1025 """Barrier operation — a no-op used for visual circuit separation.
1027 The barrier does not change the quantum state. It is recorded on the
1028 tape so that drawing backends can insert a visual separator.
1029 """
1031 _matrix = None # not a real gate
1033 def __init__(self, wires: Union[int, List[int]] = 0) -> None:
1034 """Initialise a Barrier.
1036 Args:
1037 wires: Qubit index or list of qubit indices this barrier spans.
1038 """
1039 super().__init__(wires=wires)
1041 def apply_to_state(self, state: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1042 """No-op: return the state unchanged."""
1043 return state
1045 def apply_to_state_tensor(self, psi: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1046 """No-op: return the state tensor unchanged."""
1047 return psi
1049 def apply_to_density(self, rho: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1050 """No-op: return the density matrix unchanged."""
1051 return rho
1054_PAULI_LABELS = ["I", "X", "Y", "Z"]
1055_PAULI_CLASSES = [Id, PauliX, PauliY, PauliZ]
1056_PAULI_MATRICES = {
1057 label: cls._matrix for label, cls in zip(_PAULI_LABELS, _PAULI_CLASSES)
1058}
1059_PAULI_MATS = [_PAULI_MATRICES[label] for label in _PAULI_LABELS]
1062def _make_rotation_gate(pauli_class: type, name: str) -> type:
1063 """Factory for single-qubit rotation gates RX, RY, RZ.
1065 Each gate has the form ``R_P(\\theta) = cos(\\theta/2) I - i sin(\\theta/2) P``.
1067 Args:
1068 pauli_class: One of PauliX, PauliY, PauliZ.
1069 name: Class name for the generated gate (e.g. ``"RX"``).
1071 Returns:
1072 A new :class:`Operation` subclass.
1073 """
1074 pauli_mat = pauli_class._matrix
1076 class _RotationGate(Operation):
1077 # Fancy way of setting docstring to make it generic
1078 __doc__ = (
1079 f"Rotation around the {name[1]} axis: {name}(\\theta) =\n"
1080 f"exp(-i \\theta/2 {name[1]}).\n"
1081 )
1082 _num_wires = 1
1083 _param_names = ("theta",)
1085 def __init__(
1086 self, theta: float, wires: Union[int, List[int]] = 0, **kwargs
1087 ) -> None:
1088 self.theta = theta
1089 c = jnp.cos(theta / 2)
1090 s = jnp.sin(theta / 2)
1091 mat = c * Id._matrix - 1j * s * pauli_mat
1092 super().__init__(wires=wires, matrix=mat, **kwargs)
1094 def generator(self) -> Operation:
1095 """Return the generator as the corresponding Pauli operation."""
1096 return pauli_class(wires=self.wires[0], record=False)
1098 _RotationGate.__name__ = name
1099 _RotationGate.__qualname__ = name
1100 return _RotationGate
1103RX = _make_rotation_gate(PauliX, "RX")
1104RY = _make_rotation_gate(PauliY, "RY")
1105RZ = _make_rotation_gate(PauliZ, "RZ")
1108# Projectors used by controlled-gate factories
1109_P0 = jnp.array([[1, 0], [0, 0]], dtype=_cdtype())
1110_P1 = jnp.array([[0, 0], [0, 1]], dtype=_cdtype())
1113def _make_controlled_gate(target_class: type, name: str) -> type:
1114 """Factory for controlled Pauli gates CX, CY, CZ.
1116 Each gate has the form
1117 ``CP = |0><0| \\otimes I + |1\\langle\\rangle 1| \\otimes P``.
1119 Args:
1120 target_class: The single-qubit gate class (PauliX, PauliY, PauliZ).
1121 name: Class name for the generated gate (e.g. ``"CX"``).
1123 Returns:
1124 A new :class:`Operation` subclass.
1125 """
1126 target_mat = target_class._matrix
1128 class _ControlledGate(Operation):
1129 __doc__ = (
1130 f"Controlled-{target_class.__name__[5:]} gate.\n\n"
1131 f"Applies {target_class.__name__} on the target qubit conditioned "
1132 f"on the control qubit being in state |1\\rangle."
1133 )
1134 _matrix = jnp.kron(_P0, Id._matrix) + jnp.kron(_P1, target_mat)
1135 _num_wires = 2
1136 is_controlled = True
1137 is_clifford = True # CX, CY, CZ are all Clifford gates
1139 def __init__(self, wires: List[int] = [0, 1], **kwargs) -> None:
1140 super().__init__(wires=wires, **kwargs)
1142 def decompose(self) -> List["Operation"]:
1143 # CZ = (H on target) CX (H on target). CX/CY are primitive.
1144 if name != "CZ":
1145 return super().decompose()
1146 c, t = self.wires
1147 return [
1148 H(wires=t, record=False),
1149 CX(wires=[c, t], record=False),
1150 H(wires=t, record=False),
1151 ]
1153 _ControlledGate.__name__ = name
1154 _ControlledGate.__qualname__ = name
1155 return _ControlledGate
1158CX = _make_controlled_gate(PauliX, "CX")
1159CY = _make_controlled_gate(PauliY, "CY")
1160CZ = _make_controlled_gate(PauliZ, "CZ")
1163class CCX(Operation):
1164 """Toffoli (CCX) gate.
1166 The 3-qubit Toffoli gate exercises the arbitrary-k-qubit path in
1167 :meth:`~Operation.apply_to_state` and cannot be expressed as a pair of
1168 2-qubit gates without ancilla, making it a good stress-test for the
1169 simulator.
1170 """
1172 _matrix = jnp.array(
1173 [
1174 [1, 0, 0, 0, 0, 0, 0, 0],
1175 [0, 1, 0, 0, 0, 0, 0, 0],
1176 [0, 0, 1, 0, 0, 0, 0, 0],
1177 [0, 0, 0, 1, 0, 0, 0, 0],
1178 [0, 0, 0, 0, 1, 0, 0, 0],
1179 [0, 0, 0, 0, 0, 1, 0, 0],
1180 [0, 0, 0, 0, 0, 0, 0, 1],
1181 [0, 0, 0, 0, 0, 0, 1, 0],
1182 ],
1183 dtype=_cdtype(),
1184 )
1185 is_controlled = True
1186 _num_wires = 3
1188 def __init__(self, wires: List[int] = [0, 1, 2], **kwargs) -> None:
1189 """Initialise a Toffoli (CCX) gate.
1191 Args:
1192 wires: Three-element list ``[control0, control1, target]``.
1193 """
1194 super().__init__(wires=wires, **kwargs)
1197class CSWAP(Operation):
1198 """Controlled-SWAP (Fredkin) gate.
1200 Swaps the two target qubits conditioned on the control qubit being |1\\rangle.
1202 Args on construction:
1203 wires: ``[control, target0, target1]``.
1204 """
1206 _matrix = jnp.array(
1207 [
1208 [1, 0, 0, 0, 0, 0, 0, 0],
1209 [0, 1, 0, 0, 0, 0, 0, 0],
1210 [0, 0, 1, 0, 0, 0, 0, 0],
1211 [0, 0, 0, 1, 0, 0, 0, 0],
1212 [0, 0, 0, 0, 1, 0, 0, 0],
1213 [0, 0, 0, 0, 0, 0, 1, 0],
1214 [0, 0, 0, 0, 0, 1, 0, 0],
1215 [0, 0, 0, 0, 0, 0, 0, 1],
1216 ],
1217 dtype=_cdtype(),
1218 )
1219 is_controlled = True
1220 _num_wires = 3
1222 def __init__(self, wires: List[int] = [0, 1, 2], **kwargs) -> None:
1223 """Initialise a Controlled-SWAP (Fredkin) gate.
1225 Args:
1226 wires: Three-element list ``[control, target0, target1]``.
1227 """
1228 super().__init__(wires=wires, **kwargs)
1231class ControlledPhaseShift(Operation):
1232 r"""Controlled phase shift gate (CPhase).
1234 Applies a phase shift of ``exp(i * phi)`` to the |11⟩ component of the
1235 two-qubit state, leaving all other computational basis states unchanged.
1236 This is a generalization of the CZ gate: when ``phi = \\pi`` the gate
1237 reduces to CZ.
1239 .. math::
1240 \text{CPhase}(\phi) = \text{diag}(1, 1, 1, e^{i\phi})
1242 which is equivalent to
1243 ``|0⟩⟨0| \\otimes I + |1⟩⟨1| \\otimes P(phi)`` where
1244 ``P(phi) = diag(1, exp(i*phi))``.
1245 """
1247 _num_wires = 2
1248 _param_names = ("phi",)
1249 is_controlled = True
1251 def __init__(self, phi: float, wires: List[int] = [0, 1], **kwargs) -> None:
1252 """Initialise a controlled phase shift gate.
1254 Args:
1255 phi: Phase shift angle in radians.
1256 wires: Two-element list ``[control, target]``.
1257 """
1258 self.phi = phi
1259 phase_gate = jnp.array([[1, 0], [0, jnp.exp(1j * phi)]], dtype=_cdtype())
1260 mat = jnp.kron(_P0, Id._matrix) + jnp.kron(_P1, phase_gate)
1261 super().__init__(wires=wires, matrix=mat, **kwargs)
1264class Rot(Operation):
1265 """General single-qubit rotation:
1266 Rot(\\phi, \\theta, \\omega) = RZ(\\omega) RY(\\theta) RZ(\\phi).
1268 This is the most general SU(2) rotation (up to a global phase). It
1269 decomposes into three successive rotations and has three free parameters.
1270 """
1272 _num_wires = 1
1273 _param_names = ("phi", "theta", "omega")
1275 def __init__(
1276 self,
1277 phi: float,
1278 theta: float,
1279 omega: float,
1280 wires: Union[int, List[int]] = 0,
1281 **kwargs,
1282 ) -> None:
1283 """Initialise a general rotation gate.
1285 Args:
1286 phi: First RZ rotation angle (radians).
1287 theta: RY rotation angle (radians).
1288 omega: Second RZ rotation angle (radians).
1289 wires: Qubit index or list of qubit indices this gate acts on.
1290 """
1291 self.phi = phi
1292 self.theta = theta
1293 self.omega = omega
1294 # Rot(\\phi, \theta, \\omega) = RZ(\\omega) @ RY(\theta) @ RZ(\\phi)
1295 rz_phi = jnp.cos(phi / 2) * Id._matrix - 1j * jnp.sin(phi / 2) * PauliZ._matrix
1296 ry_theta = (
1297 jnp.cos(theta / 2) * Id._matrix - 1j * jnp.sin(theta / 2) * PauliY._matrix
1298 )
1299 rz_omega = (
1300 jnp.cos(omega / 2) * Id._matrix - 1j * jnp.sin(omega / 2) * PauliZ._matrix
1301 )
1302 mat = rz_omega @ ry_theta @ rz_phi
1303 super().__init__(wires=wires, matrix=mat, **kwargs)
1305 def decompose(self) -> List["Operation"]:
1306 """Decompose into ``RZ(phi) RY(theta) RZ(omega)`` (same wire)."""
1307 w = self.wires[0]
1308 return [
1309 RZ(self.phi, wires=w, record=False),
1310 RY(self.theta, wires=w, record=False),
1311 RZ(self.omega, wires=w, record=False),
1312 ]
1315class PauliRot(Operation):
1316 """Multi-qubit Pauli rotation: exp(-i \\theta/2 P) for a Pauli word P.
1318 The Pauli word is given as a string of ``'I'``, ``'X'``, ``'Y'``, ``'Z'``
1319 characters (one per qubit). The rotation matrix is computed as
1320 ``cos(\\theta/2) I - i sin(\\theta/2) P`` where *P* is the tensor product of the
1321 corresponding single-qubit Pauli matrices.
1323 Example::
1325 PauliRot(0.5, "XY", wires=[0, 1])
1326 """
1328 _param_names = ("theta",)
1330 # Map from character to 2x2 matrix (canonical single source of truth)
1331 _PAULI_MAP = _PAULI_MATRICES
1333 def __init__(
1334 self, theta: float, pauli_word: str, wires: Union[int, List[int]] = 0, **kwargs
1335 ) -> None:
1336 """Initialise a PauliRot gate.
1338 Args:
1339 theta: Rotation angle in radians.
1340 pauli_word: A string of ``'I'``, ``'X'``, ``'Y'``, ``'Z'``
1341 characters specifying the Pauli tensor product.
1342 wires: Qubit index or list of qubit indices this gate acts on.
1343 """
1344 from functools import reduce as _reduce
1346 self.theta = theta
1347 self.pauli_word = pauli_word
1349 pauli_matrices = [self._PAULI_MAP[c] for c in pauli_word]
1350 P = _reduce(jnp.kron, pauli_matrices)
1351 dim = P.shape[0]
1352 mat = (
1353 jnp.cos(theta / 2) * jnp.eye(dim, dtype=_cdtype())
1354 - 1j * jnp.sin(theta / 2) * P
1355 )
1356 super().__init__(wires=wires, matrix=mat, **kwargs)
1358 def generator(self) -> Operation:
1359 """Return the generator Pauli tensor product as an :class:`Operation`.
1361 The generator of ``PauliRot(\\theta, word, wires)`` is the tensor product
1362 of single-qubit Pauli matrices specified by *word*. The returned
1363 :class:`Hermitian` wraps that matrix and the gate's wires.
1365 Returns:
1366 :class:`Hermitian` operation representing the Pauli tensor product.
1367 """
1368 from functools import reduce as _reduce
1370 pauli_matrices = [self._PAULI_MAP[c] for c in self.pauli_word]
1371 P = _reduce(jnp.kron, pauli_matrices)
1372 return Hermitian(matrix=P, wires=self.wires, record=False)
1375def _make_pauli_rotation_subclass(name: str, word: str) -> type:
1376 """Build a thin :class:`PauliRot` subclass with the Pauli word fixed.
1378 Used to expose multi-qubit Pauli rotations (``RXX``, ``RYY``, ``RZZ``,
1379 ``RZX``, ...) as standalone classes while sharing PauliRot's matrix
1380 construction and generator logic.
1381 """
1383 sep = " \\otimes "
1384 doc = (
1385 f"{name}(\\theta) = exp(-i \\theta/2\\, {sep.join(word)}).\n\n"
1386 f"Thin :class:`PauliRot` subclass with ``pauli_word={word!r}``."
1387 )
1389 class _PauliRotSubclass(PauliRot):
1390 __doc__ = doc
1391 _num_wires = len(word)
1393 def __init__(
1394 self,
1395 theta: float,
1396 wires: Union[int, List[int]] = None,
1397 **kwargs,
1398 ) -> None:
1399 if wires is None:
1400 wires = list(range(len(word)))
1401 super().__init__(theta, word, wires=wires, **kwargs)
1403 _PauliRotSubclass.__name__ = name
1404 _PauliRotSubclass.__qualname__ = name
1405 return _PauliRotSubclass
1408RXX = _make_pauli_rotation_subclass("RXX", "XX")
1409RYY = _make_pauli_rotation_subclass("RYY", "YY")
1410RZZ = _make_pauli_rotation_subclass("RZZ", "ZZ")
1411RZX = _make_pauli_rotation_subclass("RZX", "ZX")
1414# --- Controlled multi-qubit Pauli rotation ---------------------------------
1417class ControlledPauliRot(Operation):
1418 r"""Multi-controlled multi-qubit Pauli rotation.
1420 Applies ``PauliRot(theta, pauli_word)`` on the *target* wires
1421 conditioned on all *control* wires being in :math:`|1\rangle`.
1423 For a single control wire and a single-character Pauli word this
1424 reduces to the textbook controlled rotations ``CRX``, ``CRY``,
1425 ``CRZ`` — these are exposed below as thin subclasses.
1427 The wire layout is ``[control_0, ..., control_{n_controls-1},
1428 target_0, ..., target_{m-1}]`` where ``m = len(pauli_word)``.
1429 """
1431 _param_names = ("theta",)
1432 is_controlled = True
1434 def __init__(
1435 self,
1436 theta: float,
1437 pauli_word: str,
1438 wires: List[int],
1439 n_controls: int = 1,
1440 **kwargs,
1441 ) -> None:
1442 from functools import reduce as _reduce
1444 self.theta = theta
1445 self.pauli_word = pauli_word
1446 self.n_controls = n_controls
1448 wires_list = [wires] if isinstance(wires, int) else list(wires)
1449 n_targets = len(pauli_word)
1450 if len(wires_list) != n_controls + n_targets:
1451 raise ValueError(
1452 f"ControlledPauliRot expects {n_controls + n_targets} wires "
1453 f"({n_controls} control + {n_targets} target), got "
1454 f"{len(wires_list)}."
1455 )
1457 pauli_matrices = [PauliRot._PAULI_MAP[c] for c in pauli_word]
1458 P = _reduce(jnp.kron, pauli_matrices)
1459 d_t = P.shape[0]
1460 R = (
1461 jnp.cos(theta / 2) * jnp.eye(d_t, dtype=_cdtype())
1462 - 1j * jnp.sin(theta / 2) * P
1463 )
1465 d_c = 2**n_controls
1466 dim = d_c * d_t
1467 # All control patterns except |1...1> act trivially; the active
1468 # block sits in the last d_t x d_t slot.
1469 mat = jnp.eye(dim, dtype=_cdtype())
1470 start = (d_c - 1) * d_t
1471 mat = mat.at[start : start + d_t, start : start + d_t].set(R)
1473 super().__init__(wires=wires_list, matrix=mat, **kwargs)
1475 def generator(self) -> Operation:
1476 """Return the (Hermitian) generator on the full wire set."""
1477 from functools import reduce as _reduce
1479 pauli_matrices = [PauliRot._PAULI_MAP[c] for c in self.pauli_word]
1480 P = _reduce(jnp.kron, pauli_matrices)
1481 d_t = P.shape[0]
1482 d_c = 2**self.n_controls
1483 dim = d_c * d_t
1484 gen = jnp.zeros((dim, dim), dtype=_cdtype())
1485 start = (d_c - 1) * d_t
1486 gen = gen.at[start : start + d_t, start : start + d_t].set(P)
1487 return Hermitian(matrix=gen, wires=self.wires, record=False)
1490def _make_controlled_rotation_subclass(name: str, axis: str) -> type:
1491 """Build a single-control controlled single-qubit rotation subclass.
1493 Reproduces the historical ``CRX``, ``CRY``, ``CRZ`` API as thin
1494 :class:`ControlledPauliRot` subclasses.
1495 """
1497 class _CRotation(ControlledPauliRot):
1498 __doc__ = (
1499 f"Controlled rotation around the {axis} axis.\n\n"
1500 f"Applies R{axis}(\\theta) on the target qubit conditioned on the "
1501 f"control qubit being in state |1\\rangle.\n\n"
1502 f".. math::\n"
1503 f"{name}(\\theta) = |0\\rangle\\langle 0| \\otimes I\n"
1504 f" + |1\\rangle\\langle 1| \\otimes R{axis}(\\theta)"
1505 )
1506 _num_wires = 2
1508 def __init__(self, theta: float, wires: List[int] = [0, 1], **kwargs) -> None:
1509 super().__init__(theta, axis, wires=wires, n_controls=1, **kwargs)
1511 def decompose(self) -> List["Operation"]:
1512 """Decompose into Clifford + single-qubit Pauli rotations."""
1513 c, t = self.wires
1514 theta = self.theta
1515 if axis == "Z":
1516 return [
1517 RZ(theta / 2, wires=t, record=False),
1518 CX(wires=[c, t], record=False),
1519 RZ(-theta / 2, wires=t, record=False),
1520 CX(wires=[c, t], record=False),
1521 ]
1522 if axis == "X":
1523 return [
1524 H(wires=t, record=False),
1525 RZ(theta / 2, wires=t, record=False),
1526 CX(wires=[c, t], record=False),
1527 RZ(-theta / 2, wires=t, record=False),
1528 CX(wires=[c, t], record=False),
1529 H(wires=t, record=False),
1530 ]
1531 # axis == "Y"
1532 return [
1533 RX(-jnp.pi / 2, wires=t, record=False),
1534 RZ(theta / 2, wires=t, record=False),
1535 CX(wires=[c, t], record=False),
1536 RZ(-theta / 2, wires=t, record=False),
1537 RX(jnp.pi / 2, wires=t, record=False),
1538 ]
1540 _CRotation.__name__ = name
1541 _CRotation.__qualname__ = name
1542 return _CRotation
1545CRX = _make_controlled_rotation_subclass("CRX", "X")
1546CRY = _make_controlled_rotation_subclass("CRY", "Y")
1547CRZ = _make_controlled_rotation_subclass("CRZ", "Z")
1550class KrausChannel(Operation):
1551 """Base class for noise channels defined by a set of Kraus operators.
1553 A Kraus channel \\phi(\\rho ) = \\sigma_k K_k \\rho K_k\\dagger
1554 is the most general physical
1555 operation on a quantum state. For a pure unitary gate there is a single
1556 operator K_0 = U satisfying K_0\\daggerK_0 = I; for noisy channels there are
1557 multiple operators.
1559 Subclasses must implement :meth:`kraus_matrices` and return a list of JAX
1560 arrays. :meth:`apply_to_state` is intentionally left unimplemented:
1561 Kraus channels require a density-matrix representation and cannot be
1562 applied to a pure statevector in general.
1563 """
1565 def kraus_matrices(self) -> List[jnp.ndarray]:
1566 """Return the list of Kraus operators for this channel.
1568 Returns:
1569 List of 2-D JAX arrays, each of shape ``(2**k, 2**k)`` where k
1570 is the number of target qubits.
1572 Raises:
1573 NotImplementedError: Subclasses must override this method.
1574 """
1575 raise NotImplementedError
1577 @property
1578 def matrix(self) -> jnp.ndarray:
1579 """Raises TypeError — noise channels have no single unitary matrix.
1581 Raises:
1582 TypeError: Always raised; use :meth:`apply_to_density` instead.
1583 """
1584 raise TypeError(
1585 f"{self.__class__.__name__} is a noise channel and has no single "
1586 "unitary matrix. Use apply_to_density() instead."
1587 )
1589 def apply_to_state(self, state: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1590 """Raises TypeError — noise channels require density-matrix simulation.
1592 Args:
1593 state: Statevector (unused).
1594 n_qubits: Number of qubits (unused).
1596 Raises:
1597 TypeError: Always raised; use ``execute(type='density')`` instead.
1598 """
1599 raise TypeError(
1600 f"{self.__class__.__name__} is a noise channel and cannot be "
1601 "applied to a pure statevector. Use execute(type='density') instead."
1602 )
1604 def apply_to_state_tensor(self, psi: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1605 """Raises TypeError — noise channels require density-matrix simulation."""
1606 raise TypeError(
1607 f"{self.__class__.__name__} is a noise channel and cannot be "
1608 "applied to a pure statevector. Use execute(type='density') instead."
1609 )
1611 def apply_to_density(self, rho: jnp.ndarray, n_qubits: int) -> jnp.ndarray:
1612 """Apply
1613 \\phi(\\rho ) = \\sigma_k K_k \\rho K_k\\dagger using tensor-contraction.
1615 Uses the shared :func:`_contract_and_restore` helper, summing the
1616 result over all Kraus operators.
1618 Args:
1619 rho: Density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
1620 n_qubits: Total number of qubits in the circuit.
1622 Returns:
1623 Updated density matrix of shape ``(2**n_qubits, 2**n_qubits)``.
1624 """
1625 k = len(self.wires)
1626 dim = 2**n_qubits
1627 bra_wires = [w + n_qubits for w in self.wires]
1628 rho_out = jnp.zeros_like(rho)
1630 for K in self.kraus_matrices():
1631 K_t = K.reshape((2,) * 2 * k)
1632 K_conj_t = jnp.conj(K_t)
1633 rho_t = rho.reshape((2,) * 2 * n_qubits)
1634 rho_t = _contract_and_restore(rho_t, K_t, k, self.wires)
1635 rho_t = _contract_and_restore(rho_t, K_conj_t, k, bra_wires)
1636 rho_out = rho_out + rho_t.reshape(dim, dim)
1638 return rho_out
1641class BitFlip(KrausChannel):
1642 r"""Single-qubit bit-flip (Pauli-X) error channel.
1644 .. math::
1645 K_0 = \sqrt{1-p}\,I, \quad K_1 = \sqrt{p}\,X
1647 where *p* \\in [0, 1] is the probability of a bit flip.
1648 """
1650 _num_wires = 1
1651 _param_names = ("p",)
1653 def __init__(self, p: float, wires: Union[int, List[int]] = 0) -> None:
1654 """Initialise a bit-flip channel.
1656 Args:
1657 p: Bit-flip probability, must be in [0, 1].
1658 wires: Qubit index or list of qubit indices this channel acts on.
1660 Raises:
1661 ValueError: If *p* is outside [0, 1].
1662 """
1663 if not 0.0 <= p <= 1.0:
1664 raise ValueError("p must be in [0, 1].")
1665 self.p = p
1666 super().__init__(wires=wires)
1668 def kraus_matrices(self) -> List[jnp.ndarray]:
1669 """Return the two Kraus operators for the bit-flip channel.
1671 Returns:
1672 List ``[K0, K1]`` where K0 = \\sqrt (1-p)·I and K1 = \\sqrt p·X.
1673 """
1674 p = self.p
1675 K0 = jnp.sqrt(1 - p) * Id._matrix
1676 K1 = jnp.sqrt(p) * PauliX._matrix
1677 return [K0, K1]
1680class PhaseFlip(KrausChannel):
1681 r"""Single-qubit phase-flip (Pauli-Z) error channel.
1683 .. math::
1684 K_0 = \sqrt{1-p}\,I, \quad K_1 = \sqrt{p}\,Z
1686 where *p* \\in [0, 1] is the probability of a phase flip.
1687 """
1689 _num_wires = 1
1690 _param_names = ("p",)
1692 def __init__(self, p: float, wires: Union[int, List[int]] = 0) -> None:
1693 """Initialise a phase-flip channel.
1695 Args:
1696 p: Phase-flip probability, must be in [0, 1].
1697 wires: Qubit index or list of qubit indices this channel acts on.
1699 Raises:
1700 ValueError: If *p* is outside [0, 1].
1701 """
1702 if not 0.0 <= p <= 1.0:
1703 raise ValueError("p must be in [0, 1].")
1704 self.p = p
1705 super().__init__(wires=wires)
1707 def kraus_matrices(self) -> List[jnp.ndarray]:
1708 """Return the two Kraus operators for the phase-flip channel.
1710 Returns:
1711 List ``[K0, K1]`` where K0 = \\sqrt (1-p)·I and K1 = \\sqrt p·Z.
1712 """
1713 p = self.p
1714 K0 = jnp.sqrt(1 - p) * Id._matrix
1715 K1 = jnp.sqrt(p) * PauliZ._matrix
1716 return [K0, K1]
1719class DepolarizingChannel(KrausChannel):
1720 r"""Single-qubit depolarizing channel.
1722 .. math::
1723 K_0 = \sqrt{1-p}\,I,\quad K_1 = \sqrt{p/3}\,X,\quad
1724 K_2 = \sqrt{p/3}\,Y,\quad K_3 = \sqrt{p/3}\,Z
1726 where *p* \\in [0, 1]. At p = 3/4 the channel is fully depolarizing.
1727 """
1729 _num_wires = 1
1730 _param_names = ("p",)
1732 def __init__(self, p: float, wires: Union[int, List[int]] = 0) -> None:
1733 """Initialise a depolarizing channel.
1735 Args:
1736 p: Depolarization probability, must be in [0, 1].
1737 wires: Qubit index or list of qubit indices this channel acts on.
1739 Raises:
1740 ValueError: If *p* is outside [0, 1].
1741 """
1742 if not 0.0 <= p <= 1.0:
1743 raise ValueError("p must be in [0, 1].")
1744 self.p = p
1745 super().__init__(wires=wires)
1747 def kraus_matrices(self) -> List[jnp.ndarray]:
1748 """Return the four Kraus operators for the depolarizing channel.
1750 Returns:
1751 List ``[K0, K1, K2, K3]`` corresponding to I, X, Y, Z components.
1752 """
1753 p = self.p
1754 K0 = jnp.sqrt(1 - p) * Id._matrix
1755 K1 = jnp.sqrt(p / 3) * PauliX._matrix
1756 K2 = jnp.sqrt(p / 3) * PauliY._matrix
1757 K3 = jnp.sqrt(p / 3) * PauliZ._matrix
1758 return [K0, K1, K2, K3]
1761class AmplitudeDamping(KrausChannel):
1762 r"""Single-qubit amplitude damping channel.
1764 .. math::
1765 K_0 = \begin{pmatrix}1 & 0\\ 0 & \sqrt{1-\gamma}\end{pmatrix},\quad
1766 K_1 = \begin{pmatrix}0 & \sqrt{\gamma}\\ 0 & 0\end{pmatrix}
1768 where *\\gamma* \\in [0, 1] is the probability of
1769 energy loss (|1\\rangle -> |0\\rangle).
1770 """
1772 _num_wires = 1
1773 _param_names = ("gamma",)
1775 def __init__(self, gamma: float, wires: Union[int, List[int]] = 0) -> None:
1776 """Initialise an amplitude damping channel.
1778 Args:
1779 gamma: Energy-loss probability, must be in [0, 1].
1780 wires: Qubit index or list of qubit indices this channel acts on.
1782 Raises:
1783 ValueError: If *gamma* is outside [0, 1].
1784 """
1785 if not 0.0 <= gamma <= 1.0:
1786 raise ValueError("gamma must be in [0, 1].")
1787 self.gamma = gamma
1788 super().__init__(wires=wires)
1790 def kraus_matrices(self) -> List[jnp.ndarray]:
1791 """Return the two Kraus operators for the amplitude damping channel.
1793 Returns:
1794 List ``[K0, K1]`` as defined in the class docstring.
1795 """
1796 g = self.gamma
1797 K0 = jnp.array([[1.0, 0.0], [0.0, jnp.sqrt(1 - g)]], dtype=_cdtype())
1798 K1 = jnp.array([[0.0, jnp.sqrt(g)], [0.0, 0.0]], dtype=_cdtype())
1799 return [K0, K1]
1802class PhaseDamping(KrausChannel):
1803 r"""Single-qubit phase damping (dephasing) channel.
1805 .. math::
1806 K_0 = \begin{pmatrix}1 & 0\\ 0 & \sqrt{1-\gamma}\end{pmatrix},\quad
1807 K_1 = \begin{pmatrix}0 & 0\\ 0 & \sqrt{\gamma}\end{pmatrix}
1809 where *\\gamma* \\in [0, 1] is the phase damping probability.
1810 """
1812 _num_wires = 1
1813 _param_names = ("gamma",)
1815 def __init__(self, gamma: float, wires: Union[int, List[int]] = 0) -> None:
1816 """Initialise a phase damping channel.
1818 Args:
1819 gamma: Phase-damping probability, must be in [0, 1].
1820 wires: Qubit index or list of qubit indices this channel acts on.
1822 Raises:
1823 ValueError: If *gamma* is outside [0, 1].
1824 """
1825 if not 0.0 <= gamma <= 1.0:
1826 raise ValueError("gamma must be in [0, 1].")
1827 self.gamma = gamma
1828 super().__init__(wires=wires)
1830 def kraus_matrices(self) -> List[jnp.ndarray]:
1831 """Return the two Kraus operators for the phase damping channel.
1833 Returns:
1834 List ``[K0, K1]`` as defined in the class docstring.
1835 """
1836 g = self.gamma
1837 K0 = jnp.array([[1.0, 0.0], [0.0, jnp.sqrt(1 - g)]], dtype=_cdtype())
1838 K1 = jnp.array([[0.0, 0.0], [0.0, jnp.sqrt(g)]], dtype=_cdtype())
1839 return [K0, K1]
1842class ThermalRelaxationError(KrausChannel):
1843 r"""Single-qubit thermal relaxation error channel.
1845 Models simultaneous T_1 energy relaxation and T_2 dephasing. Two regimes
1846 are handled:
1848 T_2 <= T_1 (Markovian dephasing + reset):
1849 Six Kraus operators built from p_z (phase-flip probability), p_r0
1850 (reset-to-|0\\rangle probability) and p_r1 (reset-to-|1\\rangle probability).
1852 T_2 > T_1 (non-Markovian; Choi matrix decomposition):
1853 The Choi matrix is assembled from the relaxation/dephasing rates, then
1854 diagonalised; Kraus operators are K_i = \sqrt \lambda_i · mat(v_i).
1856 Attributes:
1857 pe: Excited-state population (thermal population of |1\\rangle).
1858 t1: T_1 longitudinal relaxation time.
1859 t2: T_2 transverse dephasing time.
1860 tg: Gate duration.
1861 """
1863 _num_wires = 1
1864 _param_names = ("pe", "t1", "t2", "tg")
1866 def __init__(
1867 self,
1868 pe: float,
1869 t1: float,
1870 t2: float,
1871 tg: float,
1872 wires: Union[int, List[int]] = 0,
1873 ) -> None:
1874 """Initialise a thermal relaxation error channel.
1876 Args:
1877 pe: Excited-state population (thermal population of |1\\rangle), in [0, 1].
1878 t1: T_1 longitudinal relaxation time, must be > 0.
1879 t2: T_2 transverse dephasing time, must be > 0 and <= 2·T_1.
1880 tg: Gate duration, must be >= 0.
1881 wires: Qubit index or list of qubit indices this channel acts on.
1883 Raises:
1884 ValueError: If any parameter violates the stated constraints.
1885 """
1886 if not 0.0 <= pe <= 1.0:
1887 raise ValueError("pe must be in [0, 1].")
1888 if t1 <= 0:
1889 raise ValueError("t1 must be > 0.")
1890 if t2 <= 0:
1891 raise ValueError("t2 must be > 0.")
1892 if t2 > 2 * t1:
1893 raise ValueError("t2 must be <= 2·t1.")
1894 if tg < 0:
1895 raise ValueError("tg must be >= 0.")
1896 self.pe = pe
1897 self.t1 = t1
1898 self.t2 = t2
1899 self.tg = tg
1900 super().__init__(wires=wires)
1902 def kraus_matrices(self) -> List[jnp.ndarray]:
1903 """Return the Kraus operators for the thermal relaxation channel.
1905 The number of operators depends on the regime:
1907 * T_2 <= T_1: six operators (identity, phase-flip, two reset-to-|0\\rangle,
1908 two reset-to-|1\\rangle).
1909 * T_2 > T_1: four operators derived from the Choi matrix eigendecomposition.
1911 Returns:
1912 List of 2x2 JAX arrays representing the Kraus operators.
1913 """
1914 pe, t1, t2, tg = self.pe, self.t1, self.t2, self.tg
1916 eT1 = jnp.exp(-tg / t1)
1917 p_reset = 1.0 - eT1
1918 eT2 = jnp.exp(-tg / t2)
1920 if t2 <= t1:
1921 # --- Case T_2 <= T_1: six Kraus operators ---
1922 pz = (1.0 - p_reset) * (1.0 - eT2 / eT1) / 2.0
1923 pr0 = (1.0 - pe) * p_reset
1924 pr1 = pe * p_reset
1925 pid = 1.0 - pz - pr0 - pr1
1927 K0 = jnp.sqrt(pid) * jnp.eye(2, dtype=_cdtype())
1928 K1 = jnp.sqrt(pz) * jnp.array([[1, 0], [0, -1]], dtype=_cdtype())
1929 K2 = jnp.sqrt(pr0) * jnp.array([[1, 0], [0, 0]], dtype=_cdtype())
1930 K3 = jnp.sqrt(pr0) * jnp.array([[0, 1], [0, 0]], dtype=_cdtype())
1931 K4 = jnp.sqrt(pr1) * jnp.array([[0, 0], [1, 0]], dtype=_cdtype())
1932 K5 = jnp.sqrt(pr1) * jnp.array([[0, 0], [0, 1]], dtype=_cdtype())
1933 return [K0, K1, K2, K3, K4, K5]
1935 else:
1936 # --- Case T_2 > T_1: Choi matrix decomposition ---
1937 # Choi matrix (column-major / reshaping convention matching PennyLane)
1938 choi = jnp.array(
1939 [
1940 [1 - pe * p_reset, 0, 0, eT2],
1941 [0, pe * p_reset, 0, 0],
1942 [0, 0, (1 - pe) * p_reset, 0],
1943 [eT2, 0, 0, 1 - (1 - pe) * p_reset],
1944 ],
1945 dtype=_cdtype(),
1946 )
1947 eigenvalues, eigenvectors = jnp.linalg.eigh(choi)
1948 # Each eigenvector (column of eigenvectors) reshaped as 2x2 -> one Kraus op
1949 kraus = []
1950 for i in range(4):
1951 lam = eigenvalues[i]
1952 vec = eigenvectors[:, i]
1953 mat = jnp.sqrt(jnp.abs(lam)) * vec.reshape(2, 2, order="F")
1954 kraus.append(mat.astype(_cdtype()))
1955 return kraus
1958class QubitChannel(KrausChannel):
1959 """Generic Kraus channel from a user-supplied list of Kraus operators.
1961 This replaces PennyLane's ``qml.QubitChannel`` and accepts an arbitrary set
1962 of Kraus matrices satisfying \\sigma_k K_k\\dagger K_k = I.
1964 Example::
1966 kraus_ops = [jnp.sqrt(0.9) * jnp.eye(2), jnp.sqrt(0.1) * PauliX._matrix]
1967 QubitChannel(kraus_ops, wires=0)
1968 """
1970 def __init__(
1971 self, kraus_ops: List[jnp.ndarray], wires: Union[int, List[int]] = 0
1972 ) -> None:
1973 """Initialise a generic Kraus channel.
1975 Args:
1976 kraus_ops: List of Kraus matrices. Each must be a square 2D array
1977 of dimension ``2**k x 2**k`` where k = ``len(wires)``.
1978 wires: Qubit index or list of qubit indices this channel acts on.
1979 """
1980 self._kraus_ops = [jnp.asarray(K, dtype=_cdtype()) for K in kraus_ops]
1981 super().__init__(wires=wires)
1983 def kraus_matrices(self) -> List[jnp.ndarray]:
1984 """Return the stored Kraus operators.
1986 Returns:
1987 List of Kraus operator matrices.
1988 """
1989 return self._kraus_ops
1992def evolve_pauli_with_clifford(
1993 clifford: Operation,
1994 pauli: Operation,
1995 adjoint_left: bool = True,
1996) -> Operation:
1997 """Compute C\\dagger P C (or C P C\\dagger) and
1998 return the result as an Operation.
2000 Both operators are first embedded into the full Hilbert space spanned by
2001 the union of their wire sets. The result is wrapped in a
2002 :class:`Hermitian` so it can be used in further algebra.
2004 Args:
2005 clifford: A Clifford gate.
2006 pauli: A Pauli / Hermitian operator.
2007 adjoint_left: If ``True``, compute C\\dagger P C; otherwise C P C\\dagger.
2009 Returns:
2010 A :class:`Hermitian` wrapping the evolved matrix.
2011 """
2012 all_wires = sorted(set(clifford.wires) | set(pauli.wires))
2013 n = len(all_wires)
2015 C = _embed_matrix(clifford.matrix, clifford.wires, all_wires, n)
2016 P = _embed_matrix(pauli.matrix, pauli.wires, all_wires, n)
2017 Cd = jnp.conj(C).T
2019 if adjoint_left:
2020 result = Cd @ P @ C
2021 else:
2022 result = C @ P @ Cd
2024 return Hermitian(matrix=result, wires=all_wires, record=False)
2027def _embed_matrix(
2028 mat: jnp.ndarray,
2029 op_wires: list,
2030 all_wires: list,
2031 n_total: int,
2032) -> jnp.ndarray:
2033 """Embed a gate matrix into a larger Hilbert space via tensor products.
2035 If the gate already acts on all wires, the matrix is returned as-is.
2036 Otherwise the gate matrix is tensored with identities on the missing
2037 wires, and the resulting matrix rows/columns are permuted so that qubit
2038 ordering matches *all_wires*.
2040 Args:
2041 mat: The gate's unitary matrix of shape ``(2**k, 2**k)`` where
2042 ``k = len(op_wires)``.
2043 op_wires: The wires the gate acts on.
2044 all_wires: The full ordered list of wires.
2045 n_total: ``len(all_wires)``.
2047 Returns:
2048 A ``(2**n_total, 2**n_total)`` matrix.
2049 """
2050 k = len(op_wires)
2051 if k == n_total and list(op_wires) == list(all_wires):
2052 return mat
2054 # Build the full-space matrix by tensoring with identities
2055 # Strategy: tensor I on missing wires, then permute
2056 missing = [w for w in all_wires if w not in op_wires]
2057 # Full matrix = mat \\otimes I_{missing}
2058 full_mat = mat
2059 for _ in missing:
2060 full_mat = jnp.kron(full_mat, jnp.eye(2, dtype=_cdtype()))
2062 # The current ordering is [op_wires..., missing...]
2063 # We need to permute to match all_wires ordering
2064 current_order = list(op_wires) + missing
2065 if current_order != list(all_wires):
2066 perm = [current_order.index(w) for w in all_wires]
2067 full_mat = _permute_matrix(full_mat, perm, n_total)
2069 return full_mat
2072def _permute_matrix(mat: jnp.ndarray, perm: list, n_qubits: int) -> jnp.ndarray:
2073 """Permute the qubit ordering of a matrix.
2075 Given a ``(2**n, 2**n)`` matrix and a permutation of ``[0..n-1]``,
2076 reorder the qubits so that qubit ``i`` moves to position ``perm[i]``.
2078 Args:
2079 mat: Square matrix of dimension ``2**n_qubits``.
2080 perm: Permutation list.
2081 n_qubits: Number of qubits.
2083 Returns:
2084 Permuted matrix of the same shape.
2085 """
2086 dim = 2**n_qubits
2087 # Reshape to tensor, permute axes, reshape back
2088 tensor = mat.reshape([2] * (2 * n_qubits))
2089 # Axes: first n_qubits are row indices, last n_qubits are column indices
2090 row_perm = perm
2091 col_perm = [p + n_qubits for p in perm]
2092 tensor = jnp.transpose(tensor, row_perm + col_perm)
2093 return tensor.reshape(dim, dim)
2096def _dominant_pauli_label(matrix: jnp.ndarray) -> Tuple[complex, str]:
2097 r"""Return the dominant Pauli term ``(coeff, label)`` of a matrix.
2099 Finds the Pauli tensor product :math:`P` (over ``I, X, Y, Z``) with the
2100 largest :math:`|c_P|`, where :math:`c_P = \mathrm{Tr}(P M) / 2^n`. Shared
2101 by :func:`pauli_decompose` and :meth:`PauliWord.from_matrix` so the
2102 brute-force search lives in one place.
2104 Args:
2105 matrix: A ``(2**n, 2**n)`` matrix.
2107 Returns:
2108 ``(coeff, label)`` with *label* a string over ``{I, X, Y, Z}``.
2109 """
2110 from itertools import product as _product
2111 from functools import reduce as _reduce
2113 dim = matrix.shape[0]
2114 n_qubits = int(jnp.round(jnp.log2(dim)))
2116 best_label = "I" * n_qubits
2117 best_coeff = 0.0
2118 for indices in _product(range(4), repeat=n_qubits):
2119 P = _reduce(jnp.kron, [_PAULI_MATS[i] for i in indices])
2120 coeff = jnp.trace(P @ matrix) / dim
2121 if jnp.abs(coeff) > jnp.abs(best_coeff):
2122 best_coeff = coeff
2123 best_label = "".join(_PAULI_LABELS[i] for i in indices)
2124 return best_coeff, best_label
2127def pauli_decompose(matrix: jnp.ndarray, wire_order: Optional[List[int]] = None):
2128 r"""Decompose a Hermitian matrix into a sum of Pauli tensor products.
2130 For an n-qubit matrix (``2**n x 2**n``), returns the dominant Pauli
2131 term (the one with the largest absolute coefficient), wrapped as an
2132 :class:`Operation`. This is sufficient for the Fourier-tree algorithm
2133 which only needs the single non-zero Pauli term produced by Clifford
2134 conjugation of a Pauli operator.
2136 The decomposition uses the trace formula:
2137 ``c_P = Tr(P · M) / 2**n``
2139 Args:
2140 matrix: A ``(2**n, 2**n)`` Hermitian matrix.
2141 wire_order: Optional list of wire indices. If ``None``, defaults
2142 to ``[0, 1, ..., n-1]``.
2144 Returns:
2145 A tuple ``(coeff, op)`` where *coeff* is the complex coefficient and
2146 *op* is the Pauli :class:`Operation` (PauliX, PauliY, PauliZ, I, or
2147 a :class:`Hermitian` for multi-qubit tensor products).
2148 """
2149 from functools import reduce as _reduce
2151 dim = matrix.shape[0]
2152 n_qubits = int(jnp.round(jnp.log2(dim)))
2154 if wire_order is None:
2155 wire_order = list(range(n_qubits))
2157 best_coeff, pauli_label = _dominant_pauli_label(matrix)
2158 label_to_idx = {label: i for i, label in enumerate(_PAULI_LABELS)}
2160 # Build the operation for the dominant term
2161 if sum(1 for ch in pauli_label if ch != "I") <= 1:
2162 # Single-qubit Pauli on one wire (or all-identity)
2163 for q, ch in enumerate(pauli_label):
2164 if ch != "I":
2165 result_op = _PAULI_CLASSES[label_to_idx[ch]](
2166 wires=wire_order[q], record=False
2167 )
2168 result_op._pauli_label = ch
2169 return best_coeff, result_op
2170 result_op = Id(wires=wire_order[0], record=False)
2171 result_op._pauli_label = "I" * n_qubits
2172 return best_coeff, result_op
2173 else:
2174 # Multi-qubit tensor product -> Hermitian with pauli label attached
2175 P = _reduce(jnp.kron, [_PAULI_MATRICES[ch] for ch in pauli_label])
2176 result_op = Hermitian(matrix=P, wires=wire_order, record=False)
2177 result_op._pauli_label = pauli_label
2178 return best_coeff, result_op
2181def pauli_string_from_operation(op: Operation) -> str:
2182 """Extract a Pauli word string from an operation.
2184 Maps ``PauliX`` -> ``"X"``, ``PauliY`` -> ``"Y"``, ``PauliZ`` -> ``"Z"``,
2185 ``I`` -> ``"I"``. For :class:`PauliRot`, returns its stored ``pauli_word``.
2186 For operations produced by :func:`pauli_decompose`, returns the stored
2187 ``_pauli_label`` attribute.
2189 Args:
2190 op: A quantum operation.
2192 Returns:
2193 A string like ``"X"``, ``"ZZ"``, etc.
2194 """
2195 if isinstance(op, PauliRot) and hasattr(op, "pauli_word"):
2196 return op.pauli_word
2197 # Check for label stored by pauli_decompose
2198 if hasattr(op, "_pauli_label"):
2199 return op._pauli_label
2200 name_map = {"PauliX": "X", "PauliY": "Y", "PauliZ": "Z", "I": "I"}
2201 if op.name in name_map:
2202 return name_map[op.name]
2203 # Fall back: decompose the matrix
2204 _, pauli_op = pauli_decompose(op.matrix, wire_order=op.wires)
2205 return pauli_op._pauli_label
2208def prod(*ops: Operation) -> Operation:
2209 """Construct the generalized product (tensor or matrix) of multiple operations.
2211 The resulting operation acts on the union of all wire sets.
2212 If the wire sets are disjoint, this is a Kronecker product.
2213 If the wire sets overlap, the corresponding matrices are multiplied.
2215 Args:
2216 *ops: Variable number of :class:`Operation` instances.
2218 Returns:
2219 A new :class:`Operation` whose matrix represents the composed
2220 operation on the unified wire set.
2221 """
2222 if not ops:
2223 raise ValueError("At least one operation must be provided to prod().")
2224 return ops[0].prod(*ops[1:])
2227# Single-qubit (x, z) bit pattern -> Pauli label, with the convention that a
2228# Pauli word is stored as i^phase * prod_q X_q^{x_q} Z_q^{z_q}.
2229# Under this convention Y = i * X * Z, so the single-qubit Y carries x=z=1.
2230_XZ_TO_LABEL = {(0, 0): "I", (1, 0): "X", (0, 1): "Z", (1, 1): "Y"}
2231_LABEL_TO_XZ = {"I": (0, 0), "X": (1, 0), "Z": (0, 1), "Y": (1, 1)}
2234class PauliWord:
2235 r"""Symbolic n-qubit Pauli operator in the stabilizer-tableau (symplectic)
2236 representation.
2238 A Pauli word is stored as
2240 .. math::
2241 P = i^{\text{phase}} \prod_{q} X_q^{x_q} Z_q^{z_q},
2243 with bit arrays ``x, z \in \{0, 1\}^n`` and an integer ``phase`` taken mod 4
2244 (tracking the scalar ``i^{phase}``). Single-qubit Paulis map as
2245 ``I=(0,0)``, ``X=(1,0)``, ``Z=(0,1)``, ``Y=(1,1)`` (since ``Y = i X Z``).
2247 This replaces the matrix-based Clifford conjugation
2248 (:func:`evolve_pauli_with_clifford` + :func:`pauli_decompose`) with O(n)
2249 symbolic updates, and is shared by both
2250 :class:`~qml_essentials.pauli.PauliCircuit` and the Fourier-tree algorithm.
2252 All operations use NumPy (integer arithmetic), not JAX — this is symbolic
2253 bookkeeping, not numeric computation.
2254 """
2256 __slots__ = ("x", "z", "phase")
2258 def __init__(self, x: np.ndarray, z: np.ndarray, phase: int = 0) -> None:
2259 """Initialise a Pauli word.
2261 Args:
2262 x: Integer/boolean array of X-component bits, length ``n_qubits``.
2263 z: Integer/boolean array of Z-component bits, length ``n_qubits``.
2264 phase: Exponent of the global ``i^{phase}`` scalar (taken mod 4).
2265 """
2266 self.x = np.asarray(x, dtype=np.int8) & 1
2267 self.z = np.asarray(z, dtype=np.int8) & 1
2268 self.phase = int(phase) % 4
2270 # ---- constructors ---------------------------------------------------
2271 @classmethod
2272 def identity(cls, n_qubits: int) -> "PauliWord":
2273 """Return the identity Pauli word on *n_qubits*."""
2274 z = np.zeros(n_qubits, dtype=np.int8)
2275 return cls(z.copy(), z, 0)
2277 @classmethod
2278 def from_pauli_string(
2279 cls, pauli_string: str, wires: List[int], n_qubits: int
2280 ) -> "PauliWord":
2281 """Build a Pauli word from a Pauli string and its wires.
2283 Args:
2284 pauli_string: String over ``{'I', 'X', 'Y', 'Z'}``; one character
2285 per entry of *wires*.
2286 wires: Qubit indices the characters act on.
2287 n_qubits: Total number of qubits in the circuit.
2289 Returns:
2290 The corresponding :class:`PauliWord`.
2291 """
2292 x = np.zeros(n_qubits, dtype=np.int8)
2293 z = np.zeros(n_qubits, dtype=np.int8)
2294 n_y = 0
2295 for ch, w in zip(pauli_string, wires):
2296 xb, zb = _LABEL_TO_XZ[ch]
2297 x[w] = xb
2298 z[w] = zb
2299 if ch == "Y":
2300 n_y += 1
2301 # Each Y contributes a factor i (Y = i X Z), accumulated into phase.
2302 return cls(x, z, n_y % 4)
2304 @classmethod
2305 def from_operation(cls, op: "Operation", n_qubits: int) -> "PauliWord":
2306 """Build a Pauli word from a Pauli-like operation.
2308 Supports :class:`PauliX`/:class:`PauliY`/:class:`PauliZ`/:class:`Id`,
2309 :class:`PauliRot` (via its ``pauli_word``), and any operation carrying a
2310 ``_pauli_label`` (e.g. produced by :func:`pauli_decompose`) or otherwise
2311 decomposable by :func:`pauli_string_from_operation`.
2313 Args:
2314 op: The operation to convert.
2315 n_qubits: Total number of qubits in the circuit.
2317 Returns:
2318 The corresponding :class:`PauliWord`.
2319 """
2320 # Cached symbolic word (e.g. attached to a Clifford-evolved observable).
2321 cached = getattr(op, "_pauli_word", None)
2322 if isinstance(cached, PauliWord) and cached.n_qubits == n_qubits:
2323 return cached
2324 if isinstance(op, PauliRot):
2325 return cls.from_pauli_string(op.pauli_word, op.wires, n_qubits)
2326 # Single-qubit Pauli rotations: generator is the corresponding Pauli.
2327 rot_to_label = {"RX": "X", "RY": "Y", "RZ": "Z"}
2328 if op.name in rot_to_label:
2329 return cls.from_pauli_string(rot_to_label[op.name], op.wires, n_qubits)
2330 name_to_label = {"PauliX": "X", "PauliY": "Y", "PauliZ": "Z", "I": "I"}
2331 if op.name in name_to_label:
2332 return cls.from_pauli_string(name_to_label[op.name], op.wires, n_qubits)
2333 pauli_str = pauli_string_from_operation(op)
2334 return cls.from_pauli_string(pauli_str, op.wires, n_qubits)
2336 @property
2337 def n_qubits(self) -> int:
2338 """Number of qubits this Pauli word spans."""
2339 return self.x.shape[0]
2341 @property
2342 def xy_mask(self) -> np.ndarray:
2343 """Boolean mask of qubits carrying an X or Y (i.e. ``x`` bits set)."""
2344 return self.x.astype(bool)
2346 @property
2347 def is_diagonal(self) -> bool:
2348 """Whether the word is diagonal (only I/Z, i.e. no X component)."""
2349 return not bool(self.x.any())
2351 # ---- algebra --------------------------------------------------------
2352 def commutes_with(self, other: "PauliWord") -> bool:
2353 """Return whether this Pauli word commutes with *other*.
2355 Two Paulis commute iff their symplectic inner product vanishes mod 2.
2356 """
2357 sp = int(np.dot(self.x, other.z) + np.dot(self.z, other.x)) % 2
2358 return sp == 0
2360 def compose(self, other: "PauliWord") -> "PauliWord":
2361 r"""Return the operator product ``self @ other`` as a new Pauli word.
2363 Uses the exact symplectic product rule
2365 .. math::
2366 (X^{x_1} Z^{z_1})(X^{x_2} Z^{z_2})
2367 = (-1)^{z_1 \cdot x_2}\, X^{x_1 \oplus x_2} Z^{z_1 \oplus z_2},
2369 combined with the ``i^{phase}`` scalars (``-1 = i^2``).
2370 """
2371 new_x = self.x ^ other.x
2372 new_z = self.z ^ other.z
2373 cross = int(np.dot(self.z, other.x))
2374 new_phase = (self.phase + other.phase + 2 * cross) % 4
2375 return PauliWord(new_x, new_z, new_phase)
2377 def conjugate_by_clifford(
2378 self, clifford: "Operation", adjoint_left: bool = False
2379 ) -> "PauliWord":
2380 r"""Return the Clifford conjugation of this Pauli word.
2382 Computes ``C P C^\dagger`` (``adjoint_left=False``) or
2383 ``C^\dagger P C`` (``adjoint_left=True``) symbolically, where *C* is one
2384 of the supported Clifford gates ``H, S, CX, CZ`` or a Pauli gate
2385 ``PauliX/PauliY/PauliZ``.
2387 The conjugation is realised by substituting the images of the
2388 single-qubit generators ``X_q`` and ``Z_q`` and re-composing in canonical
2389 order, so all phases are tracked exactly by :meth:`compose`.
2391 Args:
2392 clifford: The Clifford operation to conjugate by.
2393 adjoint_left: If ``True`` compute ``C^\dagger P C``; else
2394 ``C P C^\dagger``.
2396 Returns:
2397 The conjugated :class:`PauliWord`.
2399 Raises:
2400 NotImplementedError: If *clifford* is not a supported gate.
2401 """
2402 n = self.n_qubits
2403 name = clifford.name
2405 # Pauli gates: conjugation is just Q P Q (Q is Hermitian => Q^dagger=Q).
2406 if name in ("PauliX", "PauliY", "PauliZ"):
2407 q = PauliWord.from_operation(clifford, n)
2408 return q.compose(self).compose(q)
2410 try:
2411 images_x, images_z = self._clifford_generator_images(
2412 name, list(clifford.wires), adjoint_left, n
2413 )
2414 except NotImplementedError:
2415 # Any other Clifford (e.g. CY): fall back to the (exact) matrix
2416 # conjugation, which works for arbitrary Cliffords at O(2^n) cost.
2417 return self._conjugate_via_matrix(clifford, adjoint_left)
2419 result = PauliWord.identity(n)
2420 result.phase = self.phase
2421 for q in range(n):
2422 if self.x[q]:
2423 result = result.compose(images_x[q])
2424 if self.z[q]:
2425 result = result.compose(images_z[q])
2426 return result
2428 def _conjugate_via_matrix(
2429 self, clifford: "Operation", adjoint_left: bool
2430 ) -> "PauliWord":
2431 """Matrix-based Clifford conjugation fallback (exact, any Clifford).
2433 Used by :meth:`conjugate_by_clifford` for Cliffords without a symbolic
2434 tableau rule. Reuses :meth:`to_matrix` / :meth:`from_matrix` and the
2435 gate's dense matrix.
2436 """
2437 n = self.n_qubits
2438 C = _embed_matrix(clifford.matrix, clifford.wires, list(range(n)), n)
2439 Cd = jnp.conj(C).T
2440 mat = self.to_matrix()
2441 result = (Cd @ mat @ C) if adjoint_left else (C @ mat @ Cd)
2442 return PauliWord.from_matrix(result)
2444 @staticmethod
2445 def _clifford_generator_images(
2446 name: str, wires: List[int], adjoint_left: bool, n: int
2447 ) -> Tuple[List["PauliWord"], List["PauliWord"]]:
2448 """Images of single-qubit generators ``X_q``/``Z_q`` under a Clifford.
2450 Returns two lists (indexed by qubit) of :class:`PauliWord` giving
2451 ``C X_q C^\\dagger`` and ``C Z_q C^\\dagger`` (or the adjoint direction).
2452 Qubits outside the gate support map to themselves.
2453 """
2455 def single(label: str, q: int) -> "PauliWord":
2456 return PauliWord.from_pauli_string(label, [q], n)
2458 images_x = [single("X", q) for q in range(n)]
2459 images_z = [single("Z", q) for q in range(n)]
2461 if name == "H":
2462 w = wires[0]
2463 images_x[w] = single("Z", w) # H X H = Z
2464 images_z[w] = single("X", w) # H Z H = X
2465 elif name == "S":
2466 w = wires[0]
2467 if adjoint_left: # S^dagger X S = -Y ; S^dagger Z S = Z
2468 images_x[w] = PauliWord.from_pauli_string("Y", [w], n).compose(
2469 PauliWord(np.zeros(n, np.int8), np.zeros(n, np.int8), 2)
2470 )
2471 else: # S X S^dagger = Y ; S Z S^dagger = Z
2472 images_x[w] = single("Y", w)
2473 # images_z[w] unchanged (Z)
2474 elif name == "CX":
2475 c, t = wires
2476 images_x[c] = single("X", c).compose(single("X", t)) # X_c -> X_c X_t
2477 images_z[t] = single("Z", c).compose(single("Z", t)) # Z_t -> Z_c Z_t
2478 # X_t -> X_t and Z_c -> Z_c unchanged ; CX is Hermitian
2479 elif name == "CZ":
2480 c, t = wires
2481 images_x[c] = single("X", c).compose(single("Z", t)) # X_c -> X_c Z_t
2482 images_x[t] = single("Z", c).compose(single("X", t)) # X_t -> Z_c X_t
2483 # Z_c, Z_t unchanged ; CZ is Hermitian
2484 elif name == "SWAP":
2485 a, b = wires
2486 images_x[a], images_x[b] = single("X", b), single("X", a) # swap supports
2487 images_z[a], images_z[b] = single("Z", b), single("Z", a)
2488 else:
2489 raise NotImplementedError(f"No symbolic Clifford rule for gate '{name}'.")
2490 return images_x, images_z
2492 # ---- expectation / conversions -------------------------------------
2493 def zero_expectation(self) -> complex:
2494 r"""Return ``<0|P|0>`` for the all-zero computational basis state.
2496 Non-zero only for diagonal words (I/Z only), in which case it equals the
2497 global phase ``i^{phase}``.
2498 """
2499 if not self.is_diagonal:
2500 return 0.0 + 0.0j
2501 return complex(1j**self.phase)
2503 def to_pauli_string(self) -> str:
2504 """Return the bare Pauli string (ignoring the global phase)."""
2505 return "".join(
2506 _XZ_TO_LABEL[(int(self.x[q]), int(self.z[q]))] for q in range(self.n_qubits)
2507 )
2509 def leading_phase(self) -> complex:
2510 r"""Return the scalar ``c`` such that ``P = c * (bare Pauli string)``.
2512 Because the bare string already contains ``i^{n_Y}`` from its Y factors,
2513 ``c = i^{phase - n_Y}``.
2514 """
2515 n_y = int(((self.x == 1) & (self.z == 1)).sum())
2516 return complex(1j ** ((self.phase - n_y) % 4))
2518 def to_pauli_string_and_phase(self) -> Tuple[str, complex]:
2519 """Return ``(bare Pauli string, leading scalar phase)``."""
2520 return self.to_pauli_string(), self.leading_phase()
2522 def to_matrix(self) -> jnp.ndarray:
2523 r"""Return the dense operator matrix ``i^{phase} \bigotimes_q X^{x_q} Z^{z_q}``.
2525 The per-qubit factor is the symplectic product ``X^{x} Z^{z}`` (so the
2526 ``(1, 1)`` factor is ``XZ = -iY``; the ``Y``-vs-``XZ`` phase is carried by
2527 ``i^{phase}``). Inverse of :meth:`from_matrix`.
2528 """
2529 ident = _PAULI_MATRICES["I"]
2530 xmat = _PAULI_MATRICES["X"]
2531 zmat = _PAULI_MATRICES["Z"]
2532 mat = jnp.array([[1.0 + 0.0j]], dtype=_cdtype())
2533 for q in range(self.n_qubits):
2534 factor = (xmat if self.x[q] else ident) @ (zmat if self.z[q] else ident)
2535 mat = jnp.kron(mat, factor)
2536 return (1j**self.phase) * mat
2538 @classmethod
2539 def from_matrix(cls, matrix: jnp.ndarray) -> "PauliWord":
2540 r"""Build a Pauli word from a matrix that is a single (signed) Pauli.
2542 Recovers the dominant Pauli string and folds its (unit) coefficient
2543 ``c = i^k`` into the word's phase. Intended for matrices that are
2544 exactly a Pauli up to a ``{\pm 1, \pm i}`` scalar (e.g. the result of
2545 Clifford conjugation of a Pauli); the dominant term is returned for
2546 general inputs.
2548 Args:
2549 matrix: A ``(2**n, 2**n)`` matrix proportional to a Pauli string.
2551 Returns:
2552 The corresponding :class:`PauliWord` on ``n`` qubits.
2553 """
2554 coeff, label = _dominant_pauli_label(matrix)
2555 n = len(label)
2556 word = cls.from_pauli_string(label, list(range(n)), n)
2557 # Fold the unit coefficient c = i^k into the phase.
2558 k = int(round(np.angle(complex(coeff)) / (np.pi / 2))) % 4
2559 word.phase = (word.phase + k) % 4
2560 return word
2562 def to_list_repr(self) -> np.ndarray:
2563 """Return the legacy int list representation (I=-1, X=0, Y=1, Z=2)."""
2564 out = np.full(self.n_qubits, -1, dtype=int)
2565 for q in range(self.n_qubits):
2566 label = _XZ_TO_LABEL[(int(self.x[q]), int(self.z[q]))]
2567 out[q] = {"I": -1, "X": 0, "Y": 1, "Z": 2}[label]
2568 return out
2570 def __eq__(self, other: object) -> bool:
2571 if not isinstance(other, PauliWord):
2572 return NotImplemented
2573 return (
2574 self.phase == other.phase
2575 and np.array_equal(self.x, other.x)
2576 and np.array_equal(self.z, other.z)
2577 )
2579 def __repr__(self) -> str:
2580 phase_str = {0: "+", 1: "+i", 2: "-", 3: "-i"}[self.phase]
2581 return f"PauliWord({phase_str}{self.to_pauli_string()})"