Fix PauliEvolutionGate (using product formulas) for all-identity Pauli terms (Qiskit#13634)

* fix PauliEvo for all identities

* fix rustiq

Co-authored-by: Alexander Ivrii <[email protected]>

* fix docs

* regression test 13644

---------

Co-authored-by: Alexander Ivrii <[email protected]>

Author: Cryoris

crates/accelerate/src/circuit_library/pauli_evolution.rs

@@ -192,22 +192,22 @@ fn multi_qubit_evolution(
 /// followed by a CX-chain and then a single Pauli-Z rotation on the last qubit. Then the CX-chain
 /// is uncomputed and the inverse basis transformation applied. E.g. for the evolution under the
 /// Pauli string XIYZ we have the circuit
-///                 ┌───┐┌───────┐┌───┐
-/// 0: ─────────────┤ X ├┤ Rz(2) ├┤ X ├───────────
-///    ┌──────┐┌───┐└─┬─┘└───────┘└─┬─┘┌───┐┌────┐
-/// 1: ┤ √Xdg ├┤ X ├──■─────────────■──┤ X ├┤ √X ├
-///    └──────┘└─┬─┘                   └─┬─┘└────┘
-/// 2: ──────────┼───────────────────────┼────────
-///     ┌───┐    │                       │  ┌───┐
-/// 3: ─┤ H ├────■───────────────────────■──┤ H ├─
-///     └───┘                               └───┘
+///
+///        ┌───┐      ┌───┐┌───────┐┌───┐┌───┐
+///     0: ┤ H ├──────┤ X ├┤ Rz(2) ├┤ X ├┤ H ├────────
+///        └───┘      └─┬─┘└───────┘└─┬─┘└───┘
+///     1: ─────────────┼─────────────┼───────────────
+///        ┌────┐┌───┐  │             │  ┌───┐┌──────┐
+///     2: ┤ √X ├┤ X ├──■─────────────■──┤ X ├┤ √Xdg ├
+///        └────┘└─┬─┘                   └─┬─┘└──────┘
+///     3: ────────■───────────────────────■──────────
 ///
 /// Args:
 ///     num_qubits: The number of qubits in the Hamiltonian.
 ///     sparse_paulis: The Paulis to implement. Given in a sparse-list format with elements
-///         ``(pauli_string, qubit_indices, coefficient)``. An element of the form
-///         ``("IXYZ", [0,1,2,3], 0.2)``, for example, is interpreted in terms of qubit indices as
-///          I_q0 X_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 0.4.
+///         ``(pauli_string, qubit_indices, rz_rotation_angle)``. An element of the form
+///         ``("XIYZ", [0,1,2,3], 2)``, for example, is interpreted in terms of qubit indices as
+///         X_q0 I_q1 Y_q2 Z_q3 and will use a RZ rotation angle of 2.
 ///     insert_barriers: If ``true``, insert a barrier in between the evolution of individual
 ///         Pauli terms.
 ///     do_fountain: If ``true``, implement the CX propagation as "fountain" shape, where each
@@ -244,7 +244,7 @@ pub fn py_pauli_evolution(
         }
 
         paulis.push(pauli);
-        times.push(time); // note we do not multiply by 2 here, this is done Python side!
+        times.push(time); // note we do not multiply by 2 here, this is already done Python side!
         indices.push(tuple.get_item(1)?.extract::<Vec<u32>>()?)
     }
 
@@ -266,12 +266,12 @@ pub fn py_pauli_evolution(
         },
     );
 
-    // When handling all-identity Paulis above, we added the time as global phase.
-    // However, the all-identity Paulis should add a negative phase, as they implement
-    // exp(-i t I). We apply the negative sign here, to only do a single (-1) multiplication,
-    // instead of doing it every time we find an all-identity Pauli.
+    // When handling all-identity Paulis above, we added the RZ rotation angle as global phase,
+    // meaning that we have implemented of exp(i 2t I). However, what we want it to implement
+    // exp(-i t I). To only use a single multiplication, we apply a factor of -0.5 here.
+    // This is faster, in particular as long as the parameter expressions are in Python.
     if modified_phase {
-        global_phase = multiply_param(&global_phase, -1.0, py);
+        global_phase = multiply_param(&global_phase, -0.5, py);
     }
 
     CircuitData::from_packed_operations(py, num_qubits as u32, 0, evos, global_phase)

crates/accelerate/src/synthesis/evolution/pauli_network.rs

@@ -211,7 +211,7 @@ fn inject_rotations(
         if pauli_support_size == 0 {
             // in case of an all-identity rotation, update global phase by subtracting
             // the angle
-            global_phase = radd_param(global_phase, multiply_param(&angles[i], -1.0, py), py);
+            global_phase = radd_param(global_phase, multiply_param(&angles[i], -0.5, py), py);
             hit_paulis[i] = true;
             dag.remove_node(i);
         } else if pauli_support_size == 1 && dag.is_front_node(i) {

crates/circuit/src/operations.rs

@@ -2345,8 +2345,14 @@ pub fn add_param(param: &Param, summand: f64, py: Python) -> Param {
 }
 
 pub fn radd_param(param1: Param, param2: Param, py: Python) -> Param {
-    match [param1, param2] {
+    match [&param1, &param2] {
         [Param::Float(theta), Param::Float(lambda)] => Param::Float(theta + lambda),
+        [Param::Float(theta), Param::ParameterExpression(_lambda)] => {
+            add_param(&param2, *theta, py)
+        }
+        [Param::ParameterExpression(_theta), Param::Float(lambda)] => {
+            add_param(&param1, *lambda, py)
+        }
         [Param::ParameterExpression(theta), Param::ParameterExpression(lambda)] => {
             Param::ParameterExpression(
                 theta

qiskit/synthesis/evolution/suzuki_trotter.py

@@ -140,7 +140,7 @@ def expand(
 
         .. code-block:: text
 
-            ("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], 2)
+            ("X", [0], t), ("ZZ", [0, 1], 2t), ("X", [0], t)
 
         Note that the rotation angle contains a factor of 2, such that that evolution
         of a Pauli :math:`P` over time :math:`t`, which is :math:`e^{itP}`, is represented

releasenotes/notes/fix-pauli-evo-all-identity-b129acd854d8c391.yaml

@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    The :class:`.PauliEvolutionGate`, if used with a product formula synthesis (this is the default),
+    did not correctly handle all-identity terms in the operator. The all-identity term
+    should introduce a global phase equal to ``-evolution_time``, but was off by a factor of 2
+    and could break for parameterized times. This behavior is now fixed.
+    Fixed `#13625 <https://github.com/Qiskit/qiskit/issues/13625>`__.

test/python/circuit/library/test_evolution_gate.py

@@ -479,6 +479,36 @@ def atomic_evolution(pauli, time):
         decomposed = evo_gate.definition.decompose()
         self.assertEqual(decomposed.count_ops()["cx"], reps * 3 * 4)
 
+    def test_all_identity(self):
+        """Test circuit with all identity Paulis works correctly."""
+        evo = PauliEvolutionGate(I ^ I, time=1).definition
+        expected = QuantumCircuit(2, global_phase=-1)
+        self.assertEqual(expected, evo)
+
+    def test_global_phase(self):
+        """Test a circuit with parameterized global phase terms.
+
+        Regression test of #13625.
+        """
+        pauli = (X ^ X) + (I ^ I) + (I ^ X)
+        time = Parameter("t")
+        evo = PauliEvolutionGate(pauli, time=time)
+
+        expected = QuantumCircuit(2, global_phase=-time)
+        expected.rxx(2 * time, 0, 1)
+        expected.rx(2 * time, 0)
+
+        with self.subTest(msg="check circuit"):
+            self.assertEqual(expected, evo.definition)
+
+        # since all terms in the Pauli operator commute, we can compare to an
+        # exact matrix exponential
+        time_value = 1.76123
+        bound = evo.definition.assign_parameters([time_value])
+        exact = scipy.linalg.expm(-1j * time_value * pauli.to_matrix())
+        with self.subTest(msg="check correctness"):
+            self.assertEqual(Operator(exact), Operator(bound))
+
     def test_sympify_is_real(self):
         """Test converting the parameters to sympy is real.
 

test/python/circuit/library/test_evolved_op_ansatz.py

@@ -199,6 +199,21 @@ def test_detect_commutation(self):
         # this Hamiltonian should be split into 2 commuting groups, hence we get 2 parameters
         self.assertEqual(2, circuit.num_parameters)
 
+    def test_evolution_with_identity(self):
+        """Test a Hamiltonian containing an identity term.
+
+        Regression test of #13644.
+        """
+        hamiltonian = SparsePauliOp(["III", "IZZ", "IXI"])
+        ansatz = hamiltonian_variational_ansatz(hamiltonian, reps=1)
+        bound = ansatz.assign_parameters([1, 1])  # we have two non-commuting groups, hence 2 params
+
+        expected = QuantumCircuit(3, global_phase=-1)
+        expected.rzz(2, 0, 1)
+        expected.rx(2, 1)
+
+        self.assertEqual(expected, bound)
+
 
 def evolve(pauli_string, time):
     """Get the reference evolution circuit for a single Pauli string."""