quantumlib · Acciaccatura · Jul 11, 2026 · Jul 11, 2026 · Jul 11, 2026 · Jul 11, 2026
@@ -319,7 +319,7 @@ def sparse(self, qubits: Sequence[cirq.Qid] | None = None) -> cirq.PauliString:
 
         return pauli_string.PauliString(
             coefficient=self.coefficient,
-            qubit_pauli_map={q: PAULI_GATES[p] for q, p in zip(qubits, self.pauli_mask) if p},
+            qubit_pauli_map={q: PAULI_GATES[p] for q, p in zip(qubits, self.pauli_mask)},
         )
 
     def __str__(self) -> str:

@@ -35,7 +35,7 @@
 if TYPE_CHECKING:
     import cirq
 
-UnitPauliStringT = frozenset[tuple[raw_types.Qid, pauli_gates.Pauli]]
+UnitPauliStringT = frozenset[tuple[raw_types.Qid, pauli_gates.Pauli | identity.IdentityGate]]
 PauliSumLike = Union[
     complex, PauliString, 'PauliSum', pauli_string.SingleQubitPauliStringGateOperation
 ]
@@ -336,7 +336,7 @@ def _is_linear_dict_of_unit_pauli_string(linear_dict: value.LinearDict[UnitPauli
         for qid, pauli in k:
             if not isinstance(qid, raw_types.Qid):
                 return False
-            if not isinstance(pauli, pauli_gates.Pauli):
+            if not isinstance(pauli, pauli_gates.Pauli) and pauli != identity.I:
                 return False
 
     return True
@@ -468,7 +468,10 @@ def from_pauli_strings(cls, terms: PauliString | list[PauliString]) -> PauliSum:
             terms = [terms]
         termdict: defaultdict[UnitPauliStringT, value.Scalar] = defaultdict(lambda: 0)
         for pstring in terms:
-            key = frozenset(pstring._qubit_pauli_map.items())
+            key = frozenset(
+                list(pstring._qubit_pauli_map.items())
+                + [(qubit, identity.I) for qubit in pstring.identity_qubits]
+            )
             termdict[key] += pstring.coefficient
         return cls(linear_dict=value.LinearDict(termdict))
 

@@ -154,7 +154,7 @@ class PauliString(raw_types.Operation, Generic[TKey]):
     def __init__(
         self,
         *contents: cirq.PAULI_STRING_LIKE,
-        qubit_pauli_map: dict[TKey, cirq.Pauli] | None = None,
+        qubit_pauli_map: Mapping[TKey, cirq.Pauli | cirq.IdentityGate] | None = None,
         coefficient: cirq.TParamValComplex = 1,
     ):
         """Initializes a new `PauliString` operation.
@@ -166,9 +166,7 @@ def __init__(
                 values is given, they are each individually converted and then
                 multiplied from left to right in order.
             qubit_pauli_map: Initial dictionary mapping qubits to pauli
-                operations. Defaults to the empty dictionary. Note that, unlike
-                dictionaries passed to contents, this dictionary must not
-                contain any identity gate values. Further note that this
+                operations. Defaults to the empty dictionary. Note that this
                 argument specifies values that are logically *before* factors
                 specified in `contents`; `contents` are *right* multiplied onto
                 the values in this dictionary.
@@ -177,12 +175,22 @@ def __init__(
         Raises:
             TypeError: If the `qubit_pauli_map` has values that are not Paulis.
         """
+        self._qubit_pauli_map: dict[TKey, cirq.Pauli] = {}
+        # Maintain a list of identity qubits.
+        # Historically this class retains only qubits with Pauli operations applied
+        # to them. However, downstream calculations (e.g. applying exponentiation, etc.)
+        # relies on knowledge of where identity qubits are placed as well. To maintain
+        # existing behavior, track identity qubits under a separate field.
+        self._identity_qubits: Sequence[cirq.Qid] = []
         if _compat.__cirq_debug__.get() and qubit_pauli_map is not None:
-            for v in qubit_pauli_map.values():
-                if not isinstance(v, pauli_gates.Pauli):
-                    raise TypeError(f'{v} is not a Pauli')
+            for k, v in qubit_pauli_map.items():
+                if isinstance(v, pauli_gates.Pauli):
+                    self._qubit_pauli_map[k] = v
+                elif v == identity.I:
+                    self._identity_qubits.append(k)
+                else:
+                    raise TypeError(f'{v} is not a Pauli or identity')
 
-        self._qubit_pauli_map: dict[TKey, cirq.Pauli] = qubit_pauli_map or {}
         self._coefficient: cirq.TParamValComplex | sympy.Expr = (
             coefficient if isinstance(coefficient, sympy.Expr) else complex(coefficient)
         )
@@ -329,6 +337,11 @@ def qubits(self) -> tuple[TKey, ...]:
         """Returns a tuple of qubits on which this pauli string acts."""
         return tuple(self.keys())
 
+    @property
+    def identity_qubits(self) -> tuple[cirq.Qid, ...]:
+        """Returns a tuple of qubits that this pauli string is applying I on."""
+        return tuple(self._identity_qubits)
+
     def _circuit_diagram_info_(self, args: cirq.CircuitDiagramInfoArgs) -> list[str]:
         if not len(self._qubit_pauli_map):
             return NotImplemented

@@ -108,9 +108,35 @@ def matrix(self) -> np.ndarray:
         """
         if protocols.is_parameterized(self._exponent):
             raise ValueError("Exponent should not parameterized.")
-        ret = np.ones(1)
-        for pauli_string_exp in self:
-            ret = np.kron(ret, protocols.unitary(pauli_string_exp))
+
+        # Retrieve all qubits that are involved in producing the PauliSum used by
+        # this class, including identity qubits. These qubits will be used to
+        # compute the appropriate total size of the final unitary.
+        qubits = set()
+        for pauli_string in self._pauli_sum:
+            qubits.update(pauli_string.qubits)
+            qubits.update(pauli_string.identity_qubits)
+
+        # Calculate the identity e^{i \theta \sum P} = \prod e^{i \theta P_n}
+        ret = np.eye(2 ** len(qubits))
+        for pauli_string in self._pauli_sum:
+            # For each e^{i \theta P_n}, compute its value using the identity
+            # e^{i \theta P_n} = I \cos{(\theta)} + i P_n \sin{(\theta)}
+            # in the context of the qubits involved in the PauliSum.
+            pauli_string_qubits = set(pauli_string.qubits)
+            pauli_string_unitary = np.ones(1, dtype=float)
+            for qubit in sorted(qubits):
+                if qubit in pauli_string_qubits:
+                    pauli_string_unitary = np.kron(
+                        pauli_string_unitary, protocols.unitary(pauli_string[qubit])
+                    )
+                else:
+                    pauli_string_unitary = np.kron(pauli_string_unitary, np.eye(2))
+            theta = pauli_string.coefficient * self._multiplier
+            theta *= self._exponent
+            cos_term = np.cos(theta) * np.eye(2 ** len(qubits))
+            sin_term = np.sin(theta) * 1j * pauli_string_unitary
+            ret = ret @ (cos_term + sin_term)
         return ret
 
     @_compat.cached_method

@@ -102,6 +102,14 @@ def test_pauli_sum_exponential_parameterized_matrix_raises() -> None:
             cirq.PauliSumExponential(2j * cirq.X(q0) + 3j * cirq.Z(q1), np.pi / 2),
             np.array([[1j, 0, 0, 0], [0, -1j, 0, 0], [0, 0, 1j, 0], [0, 0, 0, -1j]]),
         ),
+        (
+            cirq.PauliSumExponential(cirq.DensePauliString('XI')(q0, q1), exponent=np.pi / 4),
+            np.array([[1, 0, 1j, 0], [0, 1, 0, 1j], [1j, 0, 1, 0], [0, 1j, 0, 1]]) / (2**0.5),
+        ),
+        (
+            cirq.PauliSumExponential(cirq.DensePauliString('IY')(q0, q1), exponent=3 * np.pi / 4),
+            np.array([[-1, 1, 0, 0], [-1, -1, 0, 0], [0, 0, -1, 1], [0, 0, -1, -1]]) / (2**0.5),
+        ),
     ),
 )
 def test_pauli_sum_exponential_has_correct_unitary(psum_exp, expected_unitary) -> None: