From b56b07a56d34bb684e386b32412db7d40fe63128 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 28 Jul 2026 15:53:22 -0600 Subject: [PATCH 1/3] Starting to fix some Trotter issues --- analysis/config_types.py | 7 +++++++ analysis/unitary.py | 9 ++++++++- common/trotter_flattened.py | 2 ++ 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/analysis/config_types.py b/analysis/config_types.py index 9947b460..66efcdc7 100644 --- a/analysis/config_types.py +++ b/analysis/config_types.py @@ -131,7 +131,13 @@ def encode_ramped_trotter(self, **kwargs): self.trotter_combine_terms = kwargs.get("trotter_combine_terms", True) self.ordering_method = kwargs.get("ordering_method", None) self.trotter_order = kwargs.get("trotter_order", None) + self.trotter_steps = kwargs.get("trotter_steps", None) self.phase_scale_factor = kwargs.get("phase_scale_factor", 1.01) + # Validate trotter_steps + if self.trotter_steps is not None: + if not isinstance(self.trotter_steps, int) or self.trotter_steps < 1: + raise ValueError( + f"trotter_steps must be a positive integer, got {self.trotter_steps}") # Validate phase_scale_factor if self.phase_scale_factor <= 0: raise ValueError( @@ -159,6 +165,7 @@ def _generate_TOML_table(self): self.save_if_present(table, "trotter_combine_terms") self.save_if_present(table, "ordering_method") self.save_if_present(table, "trotter_order") + self.save_if_present(table, "trotter_steps") self.save_if_present(table, "phase_scale_factor") return table diff --git a/analysis/unitary.py b/analysis/unitary.py index 928afca7..515c9c3f 100644 --- a/analysis/unitary.py +++ b/analysis/unitary.py @@ -292,9 +292,16 @@ def encode_ramped_trotter( logger.verbose("-- auto selecting order based on unoptimized ramp count") Nsteps = max(1, math.ceil(Nsteps0)) - logger.info(f"-- using {method} Trotter formula with {Nsteps} steps ({Nsteps0})") + logger.info(f"-- Trotter recommendation: use {method} formula with {Nsteps} steps ({Nsteps0})") logger.verbose(f"-- unoptimized ramp count = {Nsteps * ramps_per_step}") + # Check for manual override of Trotter steps + if hasattr(config_unitary, 'trotter_steps') and config_unitary.trotter_steps is not None: + Nsteps = config_unitary.trotter_steps + logger.warning(f"-- User override for number of Trotter steps: using {Nsteps} steps.") + else: + logger.info(f"-- Using recommended number of Trotter steps: using {Nsteps} steps.") + pauli_strings = hamiltonian.get_all_pauli_strings(return_as='strings') pauli_strings = reorder_paulis(pauli_strings, config_unitary.ordering_method) diff --git a/common/trotter_flattened.py b/common/trotter_flattened.py index c17e0877..74ec3bc6 100644 --- a/common/trotter_flattened.py +++ b/common/trotter_flattened.py @@ -69,6 +69,8 @@ def get_trotterization_coefficients(method): s2 = 0.5 / (4.0 - cbrt(4.0)) k = 0.5 - 4.0 * s2 return (s2, s2, s2, s2, k, k, s2, s2, s2, s2) + elif m in (6, "sixth order"): + raise NotImplementedError("Sixth order not yet available.") elif m in (8, "eighth order", "morales 2022", "morales 2025"): # Eighth-order method from Morales et al. (2022), recommended by Ostmeyer (2023); # paper modified and updated on arXiv as Morales et al. (2025) From 57baaac01ca92c5f2c9edda2710414a405348bf2 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 28 Jul 2026 17:01:43 -0600 Subject: [PATCH 2/3] Allow the use of any available Trotter method if step count is provided (no need to automatically derive step count). --- analysis/unitary.py | 199 ++++++++++++++++++++++++-------------------- 1 file changed, 108 insertions(+), 91 deletions(-) diff --git a/analysis/unitary.py b/analysis/unitary.py index 515c9c3f..c957b8eb 100644 --- a/analysis/unitary.py +++ b/analysis/unitary.py @@ -205,102 +205,119 @@ def encode_ramped_trotter( assert tevol_hbar >= 0.0 logger.info(f"-- using timestep = {tevol_hbar}") - # Compute the number of Trotter steps based on the Trotter error from arXiv:1912.08854v3 - # Using fast implementation for better performance (100-150x more samples per second) - from qhat.analysis.trotter_coefficients_fast import trotter_error_estimator_fast - - # Read error coefficient computation mode from config (with backward-compatible defaults) - error_coeff_mode = getattr(config_unitary, 'error_coeff_mode', 'monte_carlo') - error_coeff_auto_exact = getattr(config_unitary, 'error_coeff_auto_exact', False) - error_coeff_time_limit = getattr(config_unitary, 'error_coeff_time_limit', 60) - - # Validate user configuration - if error_coeff_mode not in ['monte_carlo', 'exact']: - logger.warning(f"Invalid error_coeff_mode='{error_coeff_mode}', using 'monte_carlo'") - error_coeff_mode = 'monte_carlo' - - logger.info(f"-- Error coefficient computation mode: {error_coeff_mode}") - if error_coeff_mode == 'monte_carlo' and error_coeff_auto_exact: - logger.info(f"-- Auto-switch to exact coefficient computation: enabled") - logger.info(f"-- Error coefficient time limit: {error_coeff_time_limit}s") - - # Compute error coefficients with user-specified mode - c1, c2 = trotter_error_estimator_fast( - hamiltonian.get_grouped_terms(), - error_coeff_time_limit, - mode=error_coeff_mode, - auto_exact=error_coeff_auto_exact - ) - logger.info(f"-- Trotter error coefficients: C1 = {c1}, C2 = {c2}") - assert config_unitary.energy_error is not None - logger.info(f"-- allowable energy error = {config_unitary.energy_error} Hartrees") - logger.info(f"-- effective energy range = {2 * math.pi / tevol_hbar} Hartrees") - eps_trotter = config_unitary.energy_error * tevol_hbar / (2 * math.pi) - logger.info(f"-- allowable fractional Trotter error = {eps_trotter}") - # use config_unitary.energy_error here instead of eps_trotter, because these equations are - # derived in terms of the absolute energy error rather than the fractional energy error - s1 = tevol_hbar * config_unitary.error_scale * c1 / config_unitary.energy_error - s2 = tevol_hbar * math.sqrt(config_unitary.error_scale * c2 / config_unitary.energy_error) - logger.info(f"-- Trotter step count: s1 = {s1}, s2 = {s2}") - - # Cost factors: r1=1 (single ramp), r2=2 (forward+backward) - r1 = 1 - r2 = 2 - - # Check if user explicitly requested a specific Trotter order + # Check if user provided manual step count + user_provided_steps = (hasattr(config_unitary, 'trotter_steps') and + config_unitary.trotter_steps is not None) + + # Get the requested Trotter order/method requested_order = config_unitary.trotter_order - if requested_order is not None: - # User explicitly requested a specific order - if requested_order.lower() == "first order": - method = "first order" - Nsteps0 = s1 - ramps_per_step = 1 - elif requested_order.lower() == "second order": + if user_provided_steps: + # User provided step count manually - allow any method from the implementation + Nsteps = config_unitary.trotter_steps + logger.info(f"-- User provided number of Trotter steps: {Nsteps}") + + if requested_order is None: + raise ValueError( + "When providing trotter_steps manually, you must also specify trotter_order." + ) + + # When steps are provided, we use the requested method directly + method = requested_order.lower() + logger.info(f"-- Using user-specified Trotter method: {method}") + + else: + # Auto-calculate optimal step count - only certain methods support this + # Compute the number of Trotter steps based on the Trotter error from arXiv:1912.08854v3 + # Using fast implementation for better performance (100-150x more samples per second) + from qhat.analysis.trotter_coefficients_fast import trotter_error_estimator_fast + + # Read error coefficient computation mode from config (with backward-compatible defaults) + error_coeff_mode = getattr(config_unitary, 'error_coeff_mode', 'monte_carlo') + error_coeff_auto_exact = getattr(config_unitary, 'error_coeff_auto_exact', False) + error_coeff_time_limit = getattr(config_unitary, 'error_coeff_time_limit', 60) + + # Validate user configuration + if error_coeff_mode not in ['monte_carlo', 'exact']: + logger.warning(f"Invalid error_coeff_mode='{error_coeff_mode}', using 'monte_carlo'") + error_coeff_mode = 'monte_carlo' + + logger.info(f"-- Error coefficient computation mode: {error_coeff_mode}") + if error_coeff_mode == 'monte_carlo' and error_coeff_auto_exact: + logger.info(f"-- Auto-switch to exact coefficient computation: enabled") + logger.info(f"-- Error coefficient time limit: {error_coeff_time_limit}s") + + # Compute error coefficients with user-specified mode + c1, c2 = trotter_error_estimator_fast( + hamiltonian.get_grouped_terms(), + error_coeff_time_limit, + mode=error_coeff_mode, + auto_exact=error_coeff_auto_exact + ) + logger.info(f"-- Trotter error coefficients: C1 = {c1}, C2 = {c2}") + assert config_unitary.energy_error is not None + logger.info(f"-- allowable energy error = {config_unitary.energy_error} Hartrees") + logger.info(f"-- effective energy range = {2 * math.pi / tevol_hbar} Hartrees") + eps_trotter = config_unitary.energy_error * tevol_hbar / (2 * math.pi) + logger.info(f"-- allowable fractional Trotter error = {eps_trotter}") + # use config_unitary.energy_error here instead of eps_trotter, because these equations are + # derived in terms of the absolute energy error rather than the fractional energy error + s1 = tevol_hbar * config_unitary.error_scale * c1 / config_unitary.energy_error + s2 = tevol_hbar * math.sqrt(config_unitary.error_scale * c2 / config_unitary.energy_error) + logger.info(f"-- Trotter step count: s1 = {s1}, s2 = {s2}") + + # Cost factors: r1=1 (single ramp), r2=2 (forward+backward) + r1 = 1 + r2 = 2 + + if requested_order is not None: + # User explicitly requested a specific order for auto-calculation + if requested_order.lower() == "first order": + method = "first order" + Nsteps0 = s1 + ramps_per_step = 1 + elif requested_order.lower() == "second order": + method = "second order" + Nsteps0 = s2 + ramps_per_step = 2 + elif requested_order.lower() == "fourth order": + # Fourth-order step count using C2-based heuristic + # NOTE: This is an approximation provided by Claude (not human verified) since we don't + # compute the exact C4 coefficient + x = config_unitary.error_scale * c2 / config_unitary.energy_error + s4 = tevol_hbar * x ** (1.0/3.0) + logger.warning( + "Fourth-order step count (s4) uses an APPROXIMATION based on C2, not an exact C4 " + "coefficient. Unlike first-order (C1) and second-order (C2), this is a heuristic " + "estimate and may be less accurate." + ) + logger.info(f"-- Fourth-order step count: s4 = {s4}") + method = "fourth order" + Nsteps0 = s4 + ramps_per_step = 10 + else: + raise ValueError( + f"Invalid trotter_order '{requested_order}' for automatic step calculation. " + f"Only 'first order', 'second order', and 'fourth order' have step auto-calculation. " + f"To use other methods (third, eighth order, etc.), provide trotter_steps manually." + ) + logger.info(f"-- User requested {method}") + else: + # Auto-select between first and second order based on cost + # -- Fourth order is excluded because the step estimate is not as reliable. It should be + # added here once we have a better calculation of C4 and therefore s4. method = "second order" Nsteps0 = s2 ramps_per_step = 2 - elif requested_order.lower() == "fourth order": - # Fourth-order step count using C2-based heuristic - # NOTE: This is an approximation provided by Claude (not human verified) since we don't - # compute the exact C4 coefficient - x = config_unitary.error_scale * c2 / config_unitary.energy_error - s4 = tevol_hbar * x ** (1.0/3.0) - logger.warning( - "Fourth-order step count (s4) uses an APPROXIMATION based on C2, not an exact C4 " - "coefficient. Unlike first-order (C1) and second-order (C2), this is a heuristic " - "estimate and may be less accurate." - ) - logger.info(f"-- Fourth-order step count: s4 = {s4}") - method = "fourth order" - Nsteps0 = s4 - ramps_per_step = 10 - else: - raise ValueError(f"Invalid trotter_order '{requested_order}'. Must be 'first order', 'second order', or 'fourth order'.") - logger.info(f"-- User requested {method}") - else: - # Auto-select between first and second order based on cost - # -- Fourth order is excluded because the step estimate is not as reliable. It should be - # added here once we have a better calculation of C4 and therefore s4. - method = "second order" - Nsteps0 = s2 - ramps_per_step = 2 - if r1 * s1 < r2 * s2: - method = "first order" - Nsteps0 = s1 - ramps_per_step = 1 - logger.verbose("-- auto selecting order based on unoptimized ramp count") - - Nsteps = max(1, math.ceil(Nsteps0)) - logger.info(f"-- Trotter recommendation: use {method} formula with {Nsteps} steps ({Nsteps0})") - logger.verbose(f"-- unoptimized ramp count = {Nsteps * ramps_per_step}") - - # Check for manual override of Trotter steps - if hasattr(config_unitary, 'trotter_steps') and config_unitary.trotter_steps is not None: - Nsteps = config_unitary.trotter_steps - logger.warning(f"-- User override for number of Trotter steps: using {Nsteps} steps.") - else: - logger.info(f"-- Using recommended number of Trotter steps: using {Nsteps} steps.") + if r1 * s1 < r2 * s2: + method = "first order" + Nsteps0 = s1 + ramps_per_step = 1 + logger.verbose("-- auto selecting order based on unoptimized ramp count") + + Nsteps = max(1, math.ceil(Nsteps0)) + logger.info(f"-- Trotter recommendation: use {method} formula with {Nsteps} steps ({Nsteps0})") + logger.verbose(f"-- unoptimized ramp count = {Nsteps * ramps_per_step}") pauli_strings = hamiltonian.get_all_pauli_strings(return_as='strings') pauli_strings = reorder_paulis(pauli_strings, config_unitary.ordering_method) From f9dd771552018bac0d5991a8eb3b0992297f212c Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 29 Jul 2026 11:00:46 -0600 Subject: [PATCH 3/3] Additional Trotter methods, additional names for Trotter methods, additional tests. --- common/tests/test_trotter_flattened.py | 264 +++++++++++++++++++++++++ common/trotter_flattened.py | 145 ++++++++++++-- 2 files changed, 390 insertions(+), 19 deletions(-) diff --git a/common/tests/test_trotter_flattened.py b/common/tests/test_trotter_flattened.py index d4d2e021..ad816d77 100644 --- a/common/tests/test_trotter_flattened.py +++ b/common/tests/test_trotter_flattened.py @@ -657,6 +657,70 @@ def test_get_coefficients_case_insensitive(self): coeffs3 = get_trotterization_coefficients("second order") assert coeffs1 == coeffs2 == coeffs3 + def test_get_coefficients_suzuki_aliases(self): + """Test Suzuki fourth-order method aliases.""" + coeffs1 = get_trotterization_coefficients(4) + coeffs2 = get_trotterization_coefficients("fourth order") + coeffs3 = get_trotterization_coefficients("suzuki5") + coeffs4 = get_trotterization_coefficients("suzuki 5") + coeffs5 = get_trotterization_coefficients("suzuki 1990") + coeffs6 = get_trotterization_coefficients("suzuki 1990 5-term") + assert coeffs1 == coeffs2 == coeffs3 == coeffs4 == coeffs5 == coeffs6 + assert len(coeffs1) == 10 # 5 cycles + + def test_get_coefficients_blanes_moan_fourth(self): + """Test Blanes & Moan fourth-order method (equation 46).""" + coeffs1 = get_trotterization_coefficients("bm4") + coeffs2 = get_trotterization_coefficients("blanes and moan 2002 fourth order") + coeffs3 = get_trotterization_coefficients("blanes moan 4") + assert coeffs1 == coeffs2 == coeffs3 + # Should have 12 coefficients (6 cycles, ramped format) + assert len(coeffs1) == 12 + # Verify all are real numbers + assert all(isinstance(c, (int, float)) for c in coeffs1) + + def test_get_coefficients_ostmeyer_fourth(self): + """Test Ostmeyer 2023 fourth-order method (equation 40).""" + coeffs1 = get_trotterization_coefficients("opt4") + coeffs2 = get_trotterization_coefficients("ostmeyer 2023 fourth order") + coeffs3 = get_trotterization_coefficients("optimised fourth order") + assert coeffs1 == coeffs2 == coeffs3 + # Should have 10 coefficients (5 cycles, ramped format) + assert len(coeffs1) == 10 + # Verify all are real numbers + assert all(isinstance(c, (int, float)) for c in coeffs1) + + def test_get_coefficients_blanes_moan_sixth(self): + """Test Blanes & Moan sixth-order method (equation 53).""" + coeffs1 = get_trotterization_coefficients(6) + coeffs2 = get_trotterization_coefficients("sixth order") + coeffs3 = get_trotterization_coefficients("bm6") + coeffs4 = get_trotterization_coefficients("blanes and moan 2002 sixth order") + coeffs5 = get_trotterization_coefficients("blanes moan 6") + assert coeffs1 == coeffs2 == coeffs3 == coeffs4 == coeffs5 + # Should have 20 coefficients (10 cycles, ramped format) + assert len(coeffs1) == 20 + # Verify all are real numbers + assert all(isinstance(c, (int, float)) for c in coeffs1) + + def test_get_coefficients_morales_eighth(self): + """Test Morales eighth-order method aliases.""" + coeffs1 = get_trotterization_coefficients(8) + coeffs2 = get_trotterization_coefficients("eighth order") + coeffs3 = get_trotterization_coefficients("morales") + coeffs4 = get_trotterization_coefficients("morales 2022") + coeffs5 = get_trotterization_coefficients("morales 2025") + assert coeffs1 == coeffs2 == coeffs3 == coeffs4 == coeffs5 + assert len(coeffs1) == 34 # 17 cycles + + def test_get_coefficients_verlet_aliases(self): + """Test second-order Verlet/Leapfrog aliases.""" + coeffs1 = get_trotterization_coefficients(2) + coeffs2 = get_trotterization_coefficients("second order") + coeffs3 = get_trotterization_coefficients("verlet") + coeffs4 = get_trotterization_coefficients("leapfrog") + assert coeffs1 == coeffs2 == coeffs3 == coeffs4 + # ================================================================================== # Test: Convenience constructors @@ -854,3 +918,203 @@ def test_equivalence_old_and_new_interface(self): U_old = bloq_old.tensor_contract() U_new = bloq_new.tensor_contract() assert np.allclose(U_old, U_new) + + +# ================================================================================== +# Test: New Ostmeyer 2023 methods +# ================================================================================== + +class TestOstmeyerMethods: + """Test the newly implemented methods from Ostmeyer 2023 paper.""" + + def test_blanes_moan_fourth_instantiation(self): + """Test that Blanes & Moan fourth-order method can be instantiated.""" + trotter = Trotterization.from_method( + pauli_terms=[("X", 1.0), ("Y", 1.0)], + method="bm4", + time=1.0, + num_steps=10 + ) + assert len(trotter.coefficients) == 12 + assert trotter.num_steps == 10 + + def test_blanes_moan_fourth_unitarity(self): + """Test that Blanes & Moan fourth produces unitary matrices.""" + trotter = Trotterization.from_method( + pauli_terms=[("XY", 0.5), ("YZ", 0.3)], + method="bm4", + time=1.0, + num_steps=5 + ) + U = trotter.tensor_contract() + + # Check U† U = I + U_dag_U = U.conj().T @ U + identity = np.eye(U.shape[0]) + assert np.allclose(U_dag_U, identity, atol=1e-10), \ + "Blanes & Moan fourth does not produce unitary matrix" + + def test_ostmeyer_fourth_instantiation(self): + """Test that Ostmeyer fourth-order method can be instantiated.""" + trotter = Trotterization.from_method( + pauli_terms=[("X", 1.0), ("Y", 1.0)], + method="opt4", + time=1.0, + num_steps=10 + ) + assert len(trotter.coefficients) == 10 + assert trotter.num_steps == 10 + + def test_ostmeyer_fourth_unitarity(self): + """Test that Ostmeyer fourth produces unitary matrices.""" + trotter = Trotterization.from_method( + pauli_terms=[("XY", 0.5), ("YZ", 0.3)], + method="opt4", + time=1.0, + num_steps=5 + ) + U = trotter.tensor_contract() + + # Check U† U = I + U_dag_U = U.conj().T @ U + identity = np.eye(U.shape[0]) + assert np.allclose(U_dag_U, identity, atol=1e-10), \ + "Ostmeyer fourth does not produce unitary matrix" + + def test_blanes_moan_sixth_instantiation(self): + """Test that Blanes & Moan sixth-order method can be instantiated.""" + trotter = Trotterization.from_method( + pauli_terms=[("X", 1.0), ("Y", 1.0)], + method="bm6", + time=1.0, + num_steps=10 + ) + assert len(trotter.coefficients) == 20 + assert trotter.num_steps == 10 + + def test_blanes_moan_sixth_unitarity(self): + """Test that Blanes & Moan sixth produces unitary matrices.""" + trotter = Trotterization.from_method( + pauli_terms=[("XY", 0.5), ("YZ", 0.3)], + method="bm6", + time=1.0, + num_steps=5 + ) + U = trotter.tensor_contract() + + # Check U† U = I + U_dag_U = U.conj().T @ U + identity = np.eye(U.shape[0]) + assert np.allclose(U_dag_U, identity, atol=1e-10), \ + "Blanes & Moan sixth does not produce unitary matrix" + + def test_fourth_order_methods_comparison(self): + """Compare all fourth-order methods - verify they produce valid unitaries.""" + # Test case: two non-commuting Hamiltonians + # Note: As discussed in Ostmeyer (2023), different fourth-order methods have + # different error accumulation properties. Ostmeyer's "optimised" method has + # high theoretical efficiency but unfavorable error accumulation in practice, + # often being outperformed by Blanes & Moan's method. + h_x = 1.0 + h_y = 1.0 + time = 1.0 + num_steps = 10 + + # Exact evolution + H = h_x * pauli_string_to_matrix("X") + h_y * pauli_string_to_matrix("Y") + U_exact = expm(-1j * H * time) + + # Test that all methods can be instantiated and produce valid results + methods = [ + "suzuki5", + "bm4", + "opt4" + ] + + for method_name in methods: + trotter = Trotterization.from_method( + pauli_terms=[("X", h_x), ("Y", h_y)], + method=method_name, + time=time, + num_steps=num_steps + ) + U = trotter.tensor_contract() + + # Verify unitary + U_dag_U = U.conj().T @ U + identity = np.eye(U.shape[0]) + assert np.allclose(U_dag_U, identity, atol=1e-10), \ + f"{method_name} does not produce unitary matrix" + + # Verify it's at least approximating the correct evolution + error = np.linalg.norm(U - U_exact) + assert error < 1.0, f"{method_name} error is unreasonably large" + + def test_sixth_order_better_than_fourth(self): + """Test that sixth-order is more accurate than fourth-order.""" + h_x = 1.0 + h_y = 1.0 + time = 1.0 + num_steps = 5 + + # Exact evolution + H = h_x * pauli_string_to_matrix("X") + h_y * pauli_string_to_matrix("Y") + U_exact = expm(-1j * H * time) + + # Fourth-order (Suzuki) + trotter4 = Trotterization.from_method( + pauli_terms=[("X", h_x), ("Y", h_y)], + method=4, + time=time, + num_steps=num_steps + ) + U4 = trotter4.tensor_contract() + error4 = np.linalg.norm(U4 - U_exact) + + # Sixth-order (Blanes & Moan) + trotter6 = Trotterization.from_method( + pauli_terms=[("X", h_x), ("Y", h_y)], + method=6, + time=time, + num_steps=num_steps + ) + U6 = trotter6.tensor_contract() + error6 = np.linalg.norm(U6 - U_exact) + + # Sixth-order should be more accurate + assert error6 < error4 + + def test_new_methods_resource_estimation(self): + """Test that resource estimation works for all new methods.""" + terms = [("XY", 0.5), ("YZ", 0.3)] + time = 1.0 + num_steps = 10 + + for method_name in ["bm4", "opt4", "bm6"]: + trotter = Trotterization.from_method( + pauli_terms=terms, + method=method_name, + time=time, + num_steps=num_steps + ) + tc = t_complexity(trotter) + + # Should have non-negative resource counts + assert tc.t >= 0 + assert tc.rotations >= 0 + assert tc.clifford >= 0 + + def test_coefficient_symmetry_properties(self): + """Test that new methods produce symmetric coefficient sequences.""" + # For symmetric Trotter methods, the coefficient sequence should be palindromic + # when considering the full (c,d) pairs + + for method_name in ["bm4", "opt4", "bm6"]: + coeffs = get_trotterization_coefficients(method_name) + + # Check that we have an even number (pairs of c,d) + assert len(coeffs) % 2 == 0 + + # For these methods, verify basic properties + # All coefficients should be real + assert all(isinstance(c, (int, float)) for c in coeffs) diff --git a/common/trotter_flattened.py b/common/trotter_flattened.py index 74ec3bc6..f825995e 100644 --- a/common/trotter_flattened.py +++ b/common/trotter_flattened.py @@ -26,12 +26,20 @@ def get_trotterization_coefficients(method): Args: method: Either a string/int naming a method, or a sequence of coefficients. Supported names: - - 1, "first order": Minimal first-order method (1.0,) - - 2, "second order": Minimal second-order method (0.5, 0.5) - - 3, "third order", "ruth 1983": Third-order Ruth (1983) - - "symmetrized ruth 1983": Fourth-order symmetrized Ruth - - 4, "fourth order", "suzuki 1990": Fourth-order Suzuki (1990) - - 8, "eighth order", "morales 2022", "morales 2025": Eighth-order Morales et al. + - 1, "first order": First-order (1 cycle) + - 2, "second order", "verlet", "leapfrog": Second-order Verlet/Leapfrog (1 cycle) + - 3, "third order", "ruth 1983": Third-order Ruth (1983) (5 cycles) + - "symmetrized ruth 1983": Fourth-order symmetrized Ruth (10 cycles) + - 4, "fourth order", "suzuki5", "suzuki 5", "suzuki 1990", "suzuki 1990 5-term": + Fourth-order Suzuki (1990) 5-term - DEFAULT (5 cycles) + - "bm4", "blanes and moan 2002 fourth order", "blanes moan 4": + Fourth-order Blanes & Moan (2002) - recommended by Ostmeyer (6 cycles) + - "opt4", "ostmeyer 2023 fourth order", "optimised fourth order": + Fourth-order Ostmeyer (2023) - high theoretical efficiency (5 cycles) + - 6, "sixth order", "bm6", "blanes and moan 2002 sixth order", "blanes moan 6": + Sixth-order Blanes & Moan (2002) (10 cycles) + - 8, "eighth order", "morales", "morales 2022", "morales 2025": + Eighth-order Morales et al. (2022) (17 cycles) Returns: Tuple of coefficients for the ramped Trotterization @@ -50,30 +58,129 @@ def get_trotterization_coefficients(method): m = method.lower() if m in (1, "first order"): - # Minimal first-order method + # First-order method (1 cycle) return (1.0,) - elif m in (2, "second order"): - # Minimal second-order method + elif m in (2, "second order", "verlet", "leapfrog"): + # Second-order Verlet/Leapfrog method (1 cycle) return (0.5, 0.5) elif m in (3, "third order", "ruth 1983"): - # Third-order Ruth (1983) + # Third-order Ruth (1983) (5 cycles) return (7./24., 3./8., 3./8., -25./24., 1.0) elif m == "symmetrized ruth 1983": - # Symmetrized version of Ruth (1983) that raises it to fourth order but uses twice as - # many terms; expected to be less optimal than "suzuki 1990" + # Symmetrized version of Ruth (1983) that raises it to fourth order (10 cycles) + # Less efficient than Suzuki 1990 5-term return (7./48., 3./16., 3./16., -25./48., 0.5, 0.5, -25./48., 3./16., 3./16., 7./48.) - elif m in (4, "fourth order", "suzuki 1990"): + elif m in (4, "fourth order", "suzuki5", "suzuki 5", "suzuki 1990", "suzuki 1990 5-term"): + # Fourth-order Suzuki (1990) 5-term recursion (5 cycles) - DEFAULT # Apply the recursion relation from Suzuki (1990) to the standard second-order method # to get a fourth-order method; also discussed in Hatano and Suzuki (2005) and - # Ostmeyer (2023) + # Ostmeyer (2023). Ostmeyer 2023 recommends Suzuki 5-term recursion for + # low-precision Hamiltonians with more than two terms in the Hamiltonian. s2 = 0.5 / (4.0 - cbrt(4.0)) k = 0.5 - 4.0 * s2 return (s2, s2, s2, s2, k, k, s2, s2, s2, s2) - elif m in (6, "sixth order"): - raise NotImplementedError("Sixth order not yet available.") - elif m in (8, "eighth order", "morales 2022", "morales 2025"): - # Eighth-order method from Morales et al. (2022), recommended by Ostmeyer (2023); - # paper modified and updated on arXiv as Morales et al. (2025) + elif m in ("bm4", "blanes and moan 2002 fourth order", "blanes moan 4"): + # Fourth-order Blanes & Moan (2002) from equation 46 of Ostmeyer (2023) (6 cycles) + # Ostmeyer 2023 recommends this as a reasonable default and notes it is optimal when + # there are exactly two terms in the Hamiltonian. Has favorable error accumulation. + # From paper equation 46: + a1 = 0.07920369643119569 + b1 = 0.209515106613362 + a2 = 0.353172906049774 + b2 = -0.143851773179818 + a3 = -0.04206508035771955 + b3 = 0.5 - b1 - b2 + a4 = 1.0 - 2.0 * (a1 + a2 + a3) + # By symmetry: a5 = a3, a6 = a2, a7 = a1 and b4 = b3, b5 = b2, b6 = b1 + + # Convert (a,b) to (c,d) using equations 4-7: + c1 = a1 + d1 = b1 - c1 + c2 = a2 - d1 + d2 = b2 - c2 + c3 = a3 - d2 + d3 = b3 - c3 + c4 = a4 - d3 + d4 = b3 - c4 # b4 = b3 by symmetry + c5 = a3 - d4 # a5 = a3 by symmetry + d5 = b2 - c5 # b5 = b2 by symmetry + c6 = a2 - d5 # a6 = a2 by symmetry + d6 = b1 - c6 # b6 = b1 by symmetry + + return (c1, d1, c2, d2, c3, d3, c4, d4, c5, d5, c6, d6) + elif m in ("opt4", "ostmeyer 2023 fourth order", "optimised fourth order"): + # Fourth-order Ostmeyer (2023) from equation 40 (5 cycles) + # This is the theoretically most efficient q = 5 cycles fourth-order scheme. + # However, empirically it has unfavorable error accumulation over time and is + # usually outperformed by Blanes & Moan's scheme (46). + # From paper equation 40: + a1 = 0.09257547473195787 + b1 = 0.2540996315529392 + a2 = 0.4627160310210738 + b2 = -0.1676517240119692 + a3 = 0.5 - (a1 + a2) + b3 = 1.0 - 2.0 * (b1 + b2) + # By symmetry: a4 = a2, a5 = a1 and b4 = b2, b5 = b1 + + # Convert (a,b) to (c,d) using equations 4-7: + c1 = a1 + d1 = b1 - c1 + c2 = a2 - d1 + d2 = b2 - c2 + c3 = a3 - d2 + d3 = b3 - c3 + c4 = a2 - d3 # a4 = a2 by symmetry + d4 = b2 - c4 # b4 = b2 by symmetry + c5 = a1 - d4 # a5 = a1 by symmetry + d5 = b1 - c5 # b5 = b1 by symmetry + + return (c1, d1, c2, d2, c3, d3, c4, d4, c5, d5) + elif m in (6, "sixth order", "bm6", "blanes and moan 2002 sixth order", "blanes moan 6"): + # Sixth-order Blanes & Moan (2002) from equation 53 of Ostmeyer (2023) (10 cycles) + # Ostmeyer 2023 recommends this for high precision (relative error < 10^(-4)). + # Most efficient known order n = 6 scheme. + # From paper equation 53: + a1 = 0.0502627644003922 + a2 = 0.413514300428344 + a3 = 0.04507988979439977 + a4 = -0.188054853819569 + a5 = 0.54196067845078 + a6 = 1.0 - 2.0 * (a1 + a2 + a3 + a4 + a5) + + b1 = 0.148816447901042 + b2 = -0.132385865767784 + b3 = 0.067307604692185 + b4 = 0.432666402578175 + b5 = 0.5 - (b1 + b2 + b3 + b4) + # By symmetry: b6 = b5, b7 = b4, b8 = b3, b9 = b2, b10 = b1 + + # Convert (a,b) to (c,d) using equations 4-7: + c1 = a1 + d1 = b1 - c1 + c2 = a2 - d1 + d2 = b2 - c2 + c3 = a3 - d2 + d3 = b3 - c3 + c4 = a4 - d3 + d4 = b4 - c4 + c5 = a5 - d4 + d5 = b5 - c5 + c6 = a6 - d5 + d6 = b5 - c6 # b6 = b5 by symmetry + c7 = a5 - d6 # a7 = a5 by symmetry + d7 = b4 - c7 # b7 = b4 by symmetry + c8 = a4 - d7 # a8 = a4 by symmetry + d8 = b3 - c8 # b8 = b3 by symmetry + c9 = a3 - d8 # a9 = a3 by symmetry + d9 = b2 - c9 # b9 = b2 by symmetry + c10 = a2 - d9 # a10 = a2 by symmetry + d10 = b1 - c10 # b10 = b1 by symmetry + + return (c1, d1, c2, d2, c3, d3, c4, d4, c5, d5, + c6, d6, c7, d7, c8, d8, c9, d9, c10, d10) + elif m in (8, "eighth order", "morales", "morales 2022", "morales 2025"): + # Eighth-order Morales et al. (2022) from Ostmeyer (2023) (17 cycles) + # Recommended by Ostmeyer (2023); paper modified and updated on arXiv as Morales et al. (2025) b1 = 0.12783360986284110837857554950443 b2 = 0.56148845266356446893590729572808 b3 = -0.38400573301491401473462588779099