diff --git a/libs/solvers/include/cudaq/solvers/hamiltonian_simulation/trotter.h b/libs/solvers/include/cudaq/solvers/hamiltonian_simulation/trotter.h new file mode 100644 index 00000000..86ce42e4 --- /dev/null +++ b/libs/solvers/include/cudaq/solvers/hamiltonian_simulation/trotter.h @@ -0,0 +1,64 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2026 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ +#pragma once + +#include "cudaq.h" +#include "cudaq/qis/pauli_word.h" +#include "cudaq/spin_op.h" + +#include +#include + +namespace cudaq::solvers { + +inline constexpr int first_order_trotter = 1; +inline constexpr int second_order_trotter = 2; +inline constexpr int fourth_order_trotter = 4; + +inline constexpr double forest_ruth_w1 = 1.3512071919596578; +inline constexpr double forest_ruth_w0 = -1.7024143839193153; + +/// @brief Flattened Hamiltonian representation used by Trotter routines. +/// @details Identity terms are stored separately. For H = c I + H', +/// apply_trotter() implements the product-formula circuit for H' and omits the +/// phase exp(-i c t). That phase is globally unobservable for a single +/// unconditioned state evolution, but it can matter for controlled evolution, +/// overlaps, phase estimation, Krylov/QEL moments, and other interference-based +/// algorithms. Callers in those contexts must account for identity_coefficient. +struct trotter_terms { + std::vector coefficients; + std::vector words; + double identity_coefficient = 0.0; + std::size_t num_qubits = 0; +}; + +/// @brief Extract real non-identity Pauli terms from a spin operator. +/// @details The returned identity_coefficient is not consumed by +/// apply_trotter(). It is returned so higher-level algorithms can preserve the +/// omitted exp(-i * identity_coefficient * time) phase when that phase becomes +/// observable as a relative phase. +/// @throws std::invalid_argument for coefficients with non-negligible imaginary +/// parts. +trotter_terms make_trotter_terms(const cudaq::spin_op &hamiltonian, + double coefficient_tolerance = 1e-12); + +/// @brief Apply Suzuki-Trotter evolution to a live quantum register. +/// @details This QPU-facing primitive consumes flattened Pauli terms rather +/// than cudaq::spin_op so that kernels only receive device-lowerable data. +/// Identity terms are intentionally omitted. Thus, if the original Hamiltonian +/// was H = c I + H', this applies a product-formula approximation to +/// exp(-i H' t), not the full exp(-i H t). The omitted exp(-i c t) phase is +/// harmless for ordinary expectation values of a single unconditioned evolved +/// state, but must be tracked or reintroduced by controlled/interference-based +/// algorithms. +__qpu__ void apply_trotter(const std::vector &coefficients, + const std::vector &words, + double time, std::size_t steps, int order, + cudaq::qview<> qubits); + +} // namespace cudaq::solvers diff --git a/libs/solvers/lib/CMakeLists.txt b/libs/solvers/lib/CMakeLists.txt index 4b13c0be..ab78dfbf 100644 --- a/libs/solvers/lib/CMakeLists.txt +++ b/libs/solvers/lib/CMakeLists.txt @@ -31,6 +31,11 @@ add_library(cudaq-solvers SHARED version.cpp ) +cudaqx_add_device_code(cudaq-solvers + SOURCES + hamiltonian_simulation/trotter.cpp +) + add_subdirectory(adapt) add_subdirectory(optimizers) add_subdirectory(stateprep) diff --git a/libs/solvers/lib/hamiltonian_simulation/trotter.cpp b/libs/solvers/lib/hamiltonian_simulation/trotter.cpp new file mode 100644 index 00000000..5b44d855 --- /dev/null +++ b/libs/solvers/lib/hamiltonian_simulation/trotter.cpp @@ -0,0 +1,91 @@ +/******************************************************************************* + * Copyright (c) 2026 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "cudaq/solvers/hamiltonian_simulation/trotter.h" + +#include +#include + +namespace cudaq::solvers { + +trotter_terms make_trotter_terms(const cudaq::spin_op &hamiltonian, + double coefficient_tolerance) { + if (coefficient_tolerance < 0.0) + throw std::invalid_argument( + "trotter error - coefficient tolerance must be non-negative."); + + trotter_terms terms; + terms.num_qubits = hamiltonian.num_qubits(); + + for (const auto &term : hamiltonian) { + const auto coefficient = term.evaluate_coefficient(); + if (std::abs(coefficient.imag()) > coefficient_tolerance) + throw std::invalid_argument( + "trotter error - only real Hamiltonian coefficients are supported."); + + const auto real_coefficient = coefficient.real(); + if (term.is_identity()) { + // Preserve identity contributions for callers that need the omitted + // exp(-i c t) phase in controlled or interference-based algorithms. + terms.identity_coefficient += real_coefficient; + continue; + } + + terms.coefficients.push_back(real_coefficient); + terms.words.push_back(term.get_pauli_word(terms.num_qubits)); + } + + return terms; +} + +__qpu__ void apply_trotter(const std::vector &coefficients, + const std::vector &words, + double time, std::size_t steps, int order, + cudaq::qview<> qubits) { + if (steps == 0 || coefficients.size() != words.size()) + return; + + if (order != 1 && order != 2 && order != 4) + return; + + const double dt = time / static_cast(steps); + for (std::size_t step = 0; step < steps; ++step) { + if (order == 1) { + for (std::size_t i = 0; i < words.size(); ++i) + exp_pauli(-dt * coefficients[i], qubits, words[i]); + } else if (order == 4) { + for (std::size_t i = 0; i < words.size(); ++i) + exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i], qubits, + words[i]); + for (std::size_t i = words.size(); i > 0; --i) + exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i - 1], qubits, + words[i - 1]); + + for (std::size_t i = 0; i < words.size(); ++i) + exp_pauli(-0.5 * -1.7024143839193153 * dt * coefficients[i], qubits, + words[i]); + for (std::size_t i = words.size(); i > 0; --i) + exp_pauli(-0.5 * -1.7024143839193153 * dt * coefficients[i - 1], qubits, + words[i - 1]); + + for (std::size_t i = 0; i < words.size(); ++i) + exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i], qubits, + words[i]); + for (std::size_t i = words.size(); i > 0; --i) + exp_pauli(-0.5 * 1.3512071919596578 * dt * coefficients[i - 1], qubits, + words[i - 1]); + } else { + for (std::size_t i = 0; i < words.size(); ++i) + exp_pauli(-0.5 * dt * coefficients[i], qubits, words[i]); + for (std::size_t i = words.size(); i > 0; --i) + exp_pauli(-0.5 * dt * coefficients[i - 1], qubits, words[i - 1]); + } + } +} + +} // namespace cudaq::solvers diff --git a/libs/solvers/python/bindings/solvers/py_solvers.cpp b/libs/solvers/python/bindings/solvers/py_solvers.cpp index ca2e8ae6..d0122ee3 100644 --- a/libs/solvers/python/bindings/solvers/py_solvers.cpp +++ b/libs/solvers/python/bindings/solvers/py_solvers.cpp @@ -5,6 +5,8 @@ * This source code and the accompanying materials are made available under * * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ +#include +#include #include #include #include @@ -20,6 +22,7 @@ #include "cudaq/python/PythonCppInterop.h" #include "cudaq/solvers/adapt.h" +#include "cudaq/solvers/hamiltonian_simulation/trotter.h" #include "cudaq/solvers/qaoa.h" #include "cudaq/solvers/stateprep/ceo.h" #include "cudaq/solvers/stateprep/uccgsd.h" @@ -818,9 +821,67 @@ operator pool. Only integer and list configuration values are currently supporte )#"); } +void bindTrotter(nb::module_ &mod) { + auto terms_to_tuple = [](const cudaq::spin_op &hamiltonian, + double coefficient_tolerance) { + auto terms = + cudaq::solvers::make_trotter_terms(hamiltonian, coefficient_tolerance); + return nb::make_tuple(terms.coefficients, terms.words, + terms.identity_coefficient, terms.num_qubits); + }; + + mod.def( + "make_trotter_terms", terms_to_tuple, nb::arg("hamiltonian"), + nb::arg("coefficient_tolerance") = 1e-12, + R"(Return flattened non-identity coefficients, Pauli words, identity coefficient, and qubit count. + +The identity coefficient is returned separately because apply_trotter() omits +that global phase. For H = c I + H', apply_trotter() approximates exp(-i H' t). +Callers using controlled evolution, overlaps, phase estimation, Krylov/QEL +moments, or other interference-based algorithms must account for the omitted +exp(-i c t) phase when it becomes a relative phase.)"); + + mod.def( + "make_trotter_terms", + [terms_to_tuple](const cudaq::spin_op_term &hamiltonian, + double coefficient_tolerance) { + return terms_to_tuple(cudaq::spin_op(hamiltonian), + coefficient_tolerance); + }, + nb::arg("hamiltonian"), nb::arg("coefficient_tolerance") = 1e-12); + + mod.def( + "apply_trotter", + [](const std::vector &, const std::vector &, + double, std::size_t, int, cudaq::qview<>) {}, + R"(Apply Suzuki-Trotter evolution inside a CUDA-Q kernel. + +The QPU-facing form consumes flattened terms: + + apply_trotter(coefficients, pauli_words, time, steps, order, qubits) + +Use make_trotter_terms(H) on the host to extract coefficients and words from a +SpinOperator before passing them to a kernel. Identity terms are not applied by +this primitive; use the returned identity coefficient to track or reintroduce +the omitted phase in controlled or interference-based algorithms.)"); + + cudaq::python::registerDeviceKernel( + nb::cast(mod.attr("__name__")), "apply_trotter", + cudaq::python::getMangledArgsString< + const std::vector &, const std::vector &, + double, std::size_t, int, cudaq::qview<>>()); + + cudaq::python::addDeviceKernelInterop< + const std::vector &, const std::vector &, + double, std::size_t, int, cudaq::qview<>>( + mod, "hamiltonian_simulation", "apply_trotter", + "Apply Suzuki-Trotter evolution inside a CUDA-Q kernel."); +} + void bindSolvers(nb::module_ &mod) { addStatePrepKernels(mod); + bindTrotter(mod); auto solvers = mod; //.def_submodule("solvers"); bindOperators(solvers); diff --git a/libs/solvers/python/cudaq_solvers/__init__.py b/libs/solvers/python/cudaq_solvers/__init__.py index e9c2397f..c0e7b930 100644 --- a/libs/solvers/python/cudaq_solvers/__init__.py +++ b/libs/solvers/python/cudaq_solvers/__init__.py @@ -9,10 +9,82 @@ try: from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import * from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import __version__ + from ._pycudaqx_solvers_the_suffix_matters_cudaq_solvers import make_trotter_terms as _cpp_make_trotter_terms except ImportError as _e: from ._system_dep_check import raise_if_missing_system_dep raise_if_missing_system_dep(_e, "cudaq-solvers") raise + + +def _maybe_call(value): + return value() if callable(value) else value + + +def _is_python_spin_operator(value): + return hasattr(value, "term_count") and hasattr(value, "__iter__") + + +def _is_python_spin_term(value): + return hasattr(value, "evaluate_coefficient") and hasattr( + value, "get_pauli_word") + + +def _term_coefficient(term, coefficient_tolerance): + coefficient = term.evaluate_coefficient() + if abs(coefficient.imag) > coefficient_tolerance: + raise ValueError( + "trotter error - only real Hamiltonian coefficients are supported.") + return float(coefficient.real) + + +def _term_qubit_extent(term): + max_degree = _maybe_call(getattr(term, "max_degree", -1)) + return max_degree + 1 if max_degree >= 0 else 0 + + +def make_trotter_terms(hamiltonian, coefficient_tolerance=1e-12): + """Return flattened terms for Suzuki-Trotter circuit primitives. + + Returns ``(coefficients, words, identity_coefficient, num_qubits)`` where + ``words`` are padded Pauli strings suitable for CUDA-Q kernel arguments. + + ``apply_trotter`` omits identity terms. For ``H = c I + H'``, it applies a + product-formula approximation to ``exp(-i H' t)`` and leaves the phase + ``exp(-i c t)`` to the caller. This phase cancels in ordinary expectation + values of one unconditioned evolved state, but it can matter for controlled + evolution, overlaps, phase estimation, Krylov/QEL moments, and other + interference-based algorithms. + """ + if coefficient_tolerance < 0.0: + raise ValueError( + "trotter error - coefficient tolerance must be non-negative.") + + if _is_python_spin_term(hamiltonian): + num_qubits = _term_qubit_extent(hamiltonian) + terms = [hamiltonian] + elif _is_python_spin_operator(hamiltonian): + num_qubits = int(_maybe_call(getattr(hamiltonian, "qubit_count", 0))) + terms = list(hamiltonian) + else: + return _cpp_make_trotter_terms(hamiltonian, coefficient_tolerance) + + coefficients = [] + words = [] + identity_coefficient = 0.0 + + for term in terms: + coefficient = _term_coefficient(term, coefficient_tolerance) + term_extent = _term_qubit_extent(term) + num_qubits = max(num_qubits, term_extent) + if term.is_identity(): + identity_coefficient += coefficient + continue + words.append(term.get_pauli_word(num_qubits)) + coefficients.append(coefficient) + + return coefficients, words, identity_coefficient, num_qubits + + try: from .gqe_algorithm.gqe import gqe except ImportError: diff --git a/libs/solvers/python/tests/test_trotter.py b/libs/solvers/python/tests/test_trotter.py new file mode 100644 index 00000000..dd4ed6ec --- /dev/null +++ b/libs/solvers/python/tests/test_trotter.py @@ -0,0 +1,305 @@ +# ============================================================================ # +# Copyright (c) 2026 NVIDIA Corporation & Affiliates. # +# All rights reserved. # +# # +# This source code and the accompanying materials are made available under # +# the terms of the Apache License 2.0 which accompanies this distribution. # +# ============================================================================ # + +import numpy as np +import pytest + +import cudaq +from cudaq import spin +import cudaq_solvers as solvers + +FOREST_RUTH_W1 = 1.3512071919596578 +FOREST_RUTH_W0 = -1.7024143839193153 + + +def _apply_pauli_to_vector(word, vector): + word = str(word) + result = np.zeros_like(vector, dtype=np.complex128) + n_qubits = len(word) + for basis, amplitude in enumerate(vector): + target = basis + phase = 1.0 + 0.0j + for qubit, op in enumerate(word): + bit = (basis >> qubit) & 1 + if op == "I": + continue + if op == "X": + target ^= 1 << qubit + elif op == "Y": + target ^= 1 << qubit + phase *= -1.0j if bit else 1.0j + elif op == "Z": + phase *= -1.0 if bit else 1.0 + else: + raise ValueError(op) + result[target] += phase * amplitude + return result + + +def _apply_pauli_rotation(vector, word, angle): + return (np.cos(angle) * vector - + 1.0j * np.sin(angle) * _apply_pauli_to_vector(word, vector)) + + +def _second_order_step(vector, coefficients, words, tau): + state = vector + for coefficient, word in zip(coefficients, words): + state = _apply_pauli_rotation(state, word, 0.5 * tau * coefficient) + for coefficient, word in reversed(list(zip(coefficients, words))): + state = _apply_pauli_rotation(state, word, 0.5 * tau * coefficient) + return state + + +def _simulate_trotter(coefficients, + words, + identity, + num_qubits, + time, + steps, + order, + ket, + include_identity=True): + if steps == 0: + raise ValueError("steps must be greater than zero") + if order not in (1, 2, 4): + raise ValueError("order must be one of {1, 2, 4}") + if len(coefficients) != len(words): + raise ValueError("coefficients and words must have equal length") + if ket.size != 2**num_qubits: + raise ValueError("ket length does not match num_qubits") + + state = np.array(ket, dtype=np.complex128, copy=True) + dt = time / steps + for _ in range(steps): + if order == 1: + for coefficient, word in zip(coefficients, words): + state = _apply_pauli_rotation(state, word, dt * coefficient) + elif order == 2: + state = _second_order_step(state, coefficients, words, dt) + else: + state = _second_order_step(state, coefficients, words, + FOREST_RUTH_W1 * dt) + state = _second_order_step(state, coefficients, words, + FOREST_RUTH_W0 * dt) + state = _second_order_step(state, coefficients, words, + FOREST_RUTH_W1 * dt) + + if include_identity and identity != 0.0: + state *= np.exp(-1.0j * identity * time) + return state + + +def _pauli_matrix(word): + dim = 2**len(word) + matrix = np.zeros((dim, dim), dtype=np.complex128) + for basis in range(dim): + vector = np.zeros(dim, dtype=np.complex128) + vector[basis] = 1.0 + matrix[:, basis] = _apply_pauli_to_vector(word, vector) + return matrix + + +def _exact_evolve(coefficients, words, identity, time, ket): + matrix = identity * np.eye(ket.size, dtype=np.complex128) + for coefficient, word in zip(coefficients, words): + matrix += coefficient * _pauli_matrix(str(word)) + eigenvalues, eigenvectors = np.linalg.eigh(matrix) + return eigenvectors @ (np.exp(-1.0j * time * eigenvalues) * + (eigenvectors.conj().T @ ket)) + + +def _two_qubit_product_state(rx_angle, ry_angle): + q0 = np.array([np.cos(0.5 * rx_angle), -1.0j * np.sin(0.5 * rx_angle)], + dtype=np.complex128) + q1 = np.array( + [np.cos(0.5 * ry_angle), np.sin(0.5 * ry_angle)], dtype=np.complex128) + return np.array( + [q0[basis & 1] * q1[(basis >> 1) & 1] for basis in range(4)], + dtype=np.complex128) + + +def _phase_align_error(actual, expected): + overlap = np.vdot(expected, actual) + if abs(overlap) > 0.0: + actual = actual * np.exp(-1.0j * np.angle(overlap)) + return np.linalg.norm(actual - expected) + + +def test_no_public_host_statevector_trotter_api(): + assert not hasattr(solvers, "trotter") + + +def test_make_trotter_terms_extracts_python_spin_operator(): + hamiltonian = (0.7 * spin.x(0) + 0.4 * spin.z(1) - + 0.2 * cudaq.SpinOperator.from_word("II")) + + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + + by_word = { + str(word): coefficient + for coefficient, word in zip(coefficients, words) + } + assert num_qubits == 2 + assert identity == pytest.approx(-0.2) + assert by_word["XI"] == pytest.approx(0.7) + assert by_word["IZ"] == pytest.approx(0.4) + + +def test_make_trotter_terms_accepts_single_python_spin_term(): + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + 0.5 * spin.y(2)) + + assert coefficients == pytest.approx([0.5]) + assert [str(word) for word in words] == ["IIY"] + assert identity == pytest.approx(0.0) + assert num_qubits == 3 + + +def test_product_formula_reference_improves_with_order(): + hamiltonian = (0.7 * spin.x(0) + 0.4 * spin.z(1) + + 0.31 * spin.x(0) * spin.z(1) + 0.23 * spin.y(0) * spin.y(1)) + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + + rng = np.random.default_rng(7) + ket = rng.normal(size=4) + 1.0j * rng.normal(size=4) + ket = ket / np.linalg.norm(ket) + + time = 0.8 + steps = 2 + exact = _exact_evolve(coefficients, words, identity, time, ket) + + first = _simulate_trotter(coefficients, words, identity, num_qubits, time, + steps, 1, ket) + second = _simulate_trotter(coefficients, words, identity, num_qubits, time, + steps, 2, ket) + fourth = _simulate_trotter(coefficients, words, identity, num_qubits, time, + steps, 4, ket) + + first_error = _phase_align_error(first, exact) + second_error = _phase_align_error(second, exact) + fourth_error = _phase_align_error(fourth, exact) + + assert second_error < first_error + assert fourth_error < second_error + + +def test_apply_trotter_kernel_interop_with_flattened_terms(): + hamiltonian = spin.x(0) + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + assert identity == 0.0 + assert num_qubits == 1 + + @cudaq.kernel + def evolve(coeffs: list[float], paulis: list[cudaq.pauli_word], t: float): + q = cudaq.qvector(1) + solvers.hamiltonian_simulation.apply_trotter(coeffs, paulis, t, 1, 2, q) + + state = np.asarray(cudaq.get_state(evolve, coefficients, words, 0.25), + dtype=np.complex128) + expected = _simulate_trotter(coefficients, words, identity, num_qubits, + 0.25, 1, 2, + np.array([1.0, 0.0], dtype=np.complex128)) + + np.testing.assert_allclose(state, expected, atol=1e-12) + + +def test_apply_trotter_kernel_matches_reference_for_orders(): + hamiltonian = (0.7 * spin.x(0) + 0.4 * spin.z(1) + + 0.31 * spin.x(0) * spin.z(1) + 0.23 * spin.y(0) * spin.y(1)) + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + ket = np.array([1.0, 0.0, 0.0, 0.0], dtype=np.complex128) + + @cudaq.kernel + def evolve(coeffs: list[float], paulis: list[cudaq.pauli_word], t: float, + steps: int, order: int): + q = cudaq.qvector(2) + solvers.hamiltonian_simulation.apply_trotter(coeffs, paulis, t, steps, + order, q) + + for order in (1, 2, 4): + state = np.asarray(cudaq.get_state(evolve, coefficients, words, 0.8, 3, + order), + dtype=np.complex128) + expected = _simulate_trotter(coefficients, words, identity, num_qubits, + 0.8, 3, order, ket) + np.testing.assert_allclose(state, expected, atol=1e-6) + + +def test_apply_trotter_kernel_orders_track_exact_evolution(): + hamiltonian = (0.37 * spin.x(0) - 0.22 * spin.z(1) + + 0.19 * spin.x(0) * spin.x(1) + 0.41 * spin.y(0) * spin.y(1) + + 0.13 * spin.z(0) * spin.x(1)) + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + assert identity == pytest.approx(0.0) + assert num_qubits == 2 + + time = 0.7 + steps = 4 + rx_angle = 0.37 + ry_angle = -0.52 + ket = _two_qubit_product_state(rx_angle, ry_angle) + exact = _exact_evolve(coefficients, words, identity, time, ket) + + @cudaq.kernel + def evolve(coeffs: list[float], paulis: list[cudaq.pauli_word], t: float, + n_steps: int, order: int, theta0: float, theta1: float): + q = cudaq.qvector(2) + rx(theta0, q[0]) + ry(theta1, q[1]) + solvers.hamiltonian_simulation.apply_trotter(coeffs, paulis, t, n_steps, + order, q) + + errors = {} + for order in (1, 2, 4): + state = np.asarray(cudaq.get_state(evolve, coefficients, words, time, + steps, order, rx_angle, ry_angle), + dtype=np.complex128) + errors[order] = _phase_align_error(state, exact) + + assert errors[2] < errors[1] + assert errors[4] < errors[2] + assert errors[1] < 2.0e-2 + assert errors[2] < 6.0e-4 + assert errors[4] < 5.0e-6 + + +def test_apply_trotter_kernel_handles_four_qubit_hamiltonian_with_many_terms(): + hamiltonian = (0.11 * spin.x(0) - 0.17 * spin.y(1) + 0.23 * spin.z(2) - + 0.29 * spin.x(3) + 0.31 * spin.x(0) * spin.x(1) + + 0.37 * spin.y(1) * spin.z(2) - 0.41 * spin.z(0) * spin.x(3) + + 0.43 * spin.x(0) * spin.y(2) * spin.z(3) - + 0.47 * spin.y(0) * spin.y(1) * spin.x(2) + + 0.53 * spin.z(0) * spin.x(1) * spin.y(2) * spin.z(3)) + coefficients, words, identity, num_qubits = solvers.make_trotter_terms( + hamiltonian) + + assert num_qubits == 4 + assert len(coefficients) > 8 + assert len(coefficients) == len(words) + + ket = np.zeros(16, dtype=np.complex128) + ket[0] = 1.0 + + @cudaq.kernel + def evolve(coeffs: list[float], paulis: list[cudaq.pauli_word], t: float, + steps: int, order: int): + q = cudaq.qvector(4) + solvers.hamiltonian_simulation.apply_trotter(coeffs, paulis, t, steps, + order, q) + + state = np.asarray(cudaq.get_state(evolve, coefficients, words, 0.37, 2, 2), + dtype=np.complex128) + expected = _simulate_trotter(coefficients, words, identity, num_qubits, + 0.37, 2, 2, ket) + + np.testing.assert_allclose(state, expected, atol=1e-6) diff --git a/libs/solvers/unittests/CMakeLists.txt b/libs/solvers/unittests/CMakeLists.txt index bb80834f..77796189 100644 --- a/libs/solvers/unittests/CMakeLists.txt +++ b/libs/solvers/unittests/CMakeLists.txt @@ -88,6 +88,11 @@ add_executable(test_qaoa test_qaoa.cpp) target_link_libraries(test_qaoa PRIVATE GTest::gtest_main cudaq-solvers test-kernels cudaq::cudaq) gtest_discover_tests(test_qaoa) +add_executable(test_trotter test_trotter.cpp) +target_link_libraries(test_trotter PRIVATE GTest::gtest_main cudaq-solvers cudaq::cudaq) +add_dependencies(CUDAQXSolversUnitTests test_trotter) +gtest_discover_tests(test_trotter) + add_executable(test_adapt_mpi test_adapt_mpi.cpp) target_link_libraries(test_adapt_mpi PRIVATE GTest::gtest cudaq-solvers test-kernels cudaq::cudaq) gtest_discover_tests(test_adapt_mpi) diff --git a/libs/solvers/unittests/test_trotter.cpp b/libs/solvers/unittests/test_trotter.cpp new file mode 100644 index 00000000..697576dd --- /dev/null +++ b/libs/solvers/unittests/test_trotter.cpp @@ -0,0 +1,329 @@ +/******************************************************************************* + * Copyright (c) 2026 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "cudaq/solvers/hamiltonian_simulation/trotter.h" + +#include "cudaq.h" +#include +#include +#include +#include +#include + +using namespace cudaq::spin; + +namespace { +void expect_close(std::complex actual, std::complex expected, + double tol = 1e-6) { + EXPECT_NEAR(actual.real(), expected.real(), tol); + EXPECT_NEAR(actual.imag(), expected.imag(), tol); +} + +std::vector basis_bits(std::size_t basis, std::size_t num_qubits) { + std::vector bits(num_qubits, 0); + for (std::size_t qubit = 0; qubit < num_qubits; ++qubit) + bits[num_qubits - 1 - qubit] = (basis >> qubit) & 1U; + return bits; +} + +void normalize(std::vector> &state) { + double norm = 0.0; + for (const auto &litude : state) + norm += std::norm(amplitude); + + const auto scale = 1.0 / std::sqrt(norm); + for (auto &litude : state) + amplitude *= scale; +} + +std::size_t infer_num_qubits_from_state_size(std::size_t state_size) { + if (state_size == 0) + throw std::invalid_argument("input ket cannot be empty."); + + if ((state_size & (state_size - 1)) != 0) + throw std::invalid_argument("input ket length must be a power of two."); + + std::size_t num_qubits = 0; + while ((std::size_t{1} << num_qubits) < state_size) + ++num_qubits; + return num_qubits; +} + +void validate_trotter_inputs(const cudaq::solvers::trotter_terms &terms, + std::size_t steps, int order, + std::size_t state_size) { + if (steps == 0) + throw std::invalid_argument("steps must be greater than zero."); + + if (order != cudaq::solvers::first_order_trotter && + order != cudaq::solvers::second_order_trotter && + order != cudaq::solvers::fourth_order_trotter) + throw std::invalid_argument("order must be one of {1, 2, 4}."); + + if (terms.coefficients.size() != terms.words.size()) + throw std::invalid_argument( + "coefficients and Pauli words must have equal length."); + + const auto state_num_qubits = infer_num_qubits_from_state_size(state_size); + const auto num_qubits = + terms.num_qubits == 0 ? state_num_qubits : terms.num_qubits; + if (state_size != (std::size_t{1} << num_qubits)) + throw std::invalid_argument( + "input ket length does not match Hamiltonian qubits."); +} + +void apply_pauli_rotation(std::vector> &state, + const cudaq::pauli_word &word, double angle, + std::size_t num_qubits) { + const auto word_string = word.str(); + if (word_string.size() != num_qubits) + throw std::invalid_argument( + "Pauli word length does not match qubit count."); + + const auto c = std::cos(angle); + const auto s = std::sin(angle); + const std::complex minus_i_s{0.0, -s}; + std::vector> rotated(state.size(), {0.0, 0.0}); + + for (std::size_t basis = 0; basis < state.size(); ++basis) { + std::size_t target = basis; + std::complex phase{1.0, 0.0}; + + for (std::size_t qubit = 0; qubit < num_qubits; ++qubit) { + const bool bit = ((basis >> qubit) & 1U) != 0; + switch (word_string[num_qubits - 1 - qubit]) { + case 'I': + break; + case 'X': + target ^= (std::size_t{1} << qubit); + break; + case 'Y': + target ^= (std::size_t{1} << qubit); + phase *= bit ? std::complex{0.0, -1.0} + : std::complex{0.0, 1.0}; + break; + case 'Z': + phase *= bit ? -1.0 : 1.0; + break; + default: + throw std::invalid_argument("unsupported Pauli word."); + } + } + + rotated[basis] += c * state[basis]; + rotated[target] += minus_i_s * phase * state[basis]; + } + + state = std::move(rotated); +} + +void apply_first_order_step(std::vector> &state, + const cudaq::solvers::trotter_terms &terms, + double tau) { + for (std::size_t i = 0; i < terms.words.size(); ++i) + apply_pauli_rotation(state, terms.words[i], tau * terms.coefficients[i], + terms.num_qubits); +} + +void apply_second_order_step(std::vector> &state, + const cudaq::solvers::trotter_terms &terms, + double tau) { + for (std::size_t i = 0; i < terms.words.size(); ++i) + apply_pauli_rotation(state, terms.words[i], + 0.5 * tau * terms.coefficients[i], terms.num_qubits); + + for (std::size_t i = terms.words.size(); i > 0; --i) + apply_pauli_rotation(state, terms.words[i - 1], + 0.5 * tau * terms.coefficients[i - 1], + terms.num_qubits); +} + +void apply_ordered_step(std::vector> &state, + const cudaq::solvers::trotter_terms &terms, double tau, + int order) { + if (order == cudaq::solvers::first_order_trotter) { + apply_first_order_step(state, terms, tau); + return; + } + + if (order == cudaq::solvers::second_order_trotter) { + apply_second_order_step(state, terms, tau); + return; + } + + apply_second_order_step(state, terms, cudaq::solvers::forest_ruth_w1 * tau); + apply_second_order_step(state, terms, cudaq::solvers::forest_ruth_w0 * tau); + apply_second_order_step(state, terms, cudaq::solvers::forest_ruth_w1 * tau); +} + +std::vector> +simulate_trotter_statevector(cudaq::solvers::trotter_terms terms, double time, + std::size_t steps, int order, + const std::vector> &ket) { + validate_trotter_inputs(terms, steps, order, ket.size()); + + const auto state_num_qubits = infer_num_qubits_from_state_size(ket.size()); + if (terms.num_qubits == 0) + terms.num_qubits = state_num_qubits; + + auto state = ket; + const auto dt = time / static_cast(steps); + for (std::size_t step = 0; step < steps; ++step) + apply_ordered_step(state, terms, dt, order); + + if (terms.identity_coefficient != 0.0) { + const std::complex phase = + std::exp(std::complex{0.0, -terms.identity_coefficient * time}); + for (auto &litude : state) + amplitude *= phase; + } + + return state; +} + +struct apply_trotter_test_kernel { + void operator()(std::size_t num_qubits, + const std::vector &coefficients, + const std::vector &words, double time, + std::size_t steps, int order) __qpu__ { + cudaq::qvector q(num_qubits); + cudaq::solvers::apply_trotter(coefficients, words, time, steps, order, q); + } +}; +} // namespace + +TEST(TrotterTest, MakeTrotterTermsDropsIdentityTerms) { + cudaq::spin_op hamiltonian = + 1.2 * x(0) + 0.4 * z(1) - 0.7 * cudaq::spin_op::from_word("II"); + + auto terms = cudaq::solvers::make_trotter_terms(hamiltonian); + EXPECT_EQ(terms.coefficients.size(), 2); + EXPECT_EQ(terms.words.size(), 2); + EXPECT_NEAR(terms.identity_coefficient, -0.7, 1e-12); + EXPECT_EQ(terms.num_qubits, 2); +} + +TEST(TrotterTest, MakeTrotterTermsRejectsInvalidInputs) { + EXPECT_THROW(cudaq::solvers::make_trotter_terms(z(0), -1.0), + std::invalid_argument); + + cudaq::spin_op complex_hamiltonian = std::complex{0.0, 0.25} * x(0); + EXPECT_THROW(cudaq::solvers::make_trotter_terms(complex_hamiltonian), + std::invalid_argument); +} + +TEST(TrotterTest, TestReferenceRejectsInvalidInputs) { + auto terms = cudaq::solvers::make_trotter_terms(z(0)); + std::vector> ket{{1.0, 0.0}, {0.0, 0.0}}; + + EXPECT_THROW(simulate_trotter_statevector(terms, 0.25, 0, 2, ket), + std::invalid_argument); + EXPECT_THROW(simulate_trotter_statevector(terms, 0.25, 1, 3, ket), + std::invalid_argument); + EXPECT_THROW( + simulate_trotter_statevector( + terms, 0.25, 1, 2, std::vector>{{1.0, 0.0}}), + std::invalid_argument); + + auto bad_terms = terms; + bad_terms.coefficients.push_back(0.5); + EXPECT_THROW(simulate_trotter_statevector(bad_terms, 0.25, 1, 2, ket), + std::invalid_argument); +} + +TEST(TrotterTest, TestReferenceIncludesIdentityGlobalPhase) { + auto terms = + cudaq::solvers::make_trotter_terms(2.0 * cudaq::spin_op::from_word("II")); + std::vector> ket{ + {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}; + const double time = 0.3; + + auto evolved = simulate_trotter_statevector( + terms, time, 1, cudaq::solvers::second_order_trotter, ket); + const auto phase = std::exp(std::complex{0.0, -2.0 * time}); + + expect_close(evolved[0], phase); + expect_close(evolved[1], {0.0, 0.0}); + expect_close(evolved[2], {0.0, 0.0}); + expect_close(evolved[3], {0.0, 0.0}); +} + +TEST(TrotterTest, ApplyTrotterKernelMatchesStatevectorReferenceForOrders) { + cudaq::spin_op hamiltonian = 0.7 * cudaq::spin_op::from_word("XI") + + 0.4 * cudaq::spin_op::from_word("IZ") + + 0.31 * cudaq::spin_op::from_word("XZ") + + 0.23 * cudaq::spin_op::from_word("YY"); + auto terms = cudaq::solvers::make_trotter_terms(hamiltonian); + + std::vector> ket{ + {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}; + + const double time = 0.8; + const std::size_t steps = 3; + for (int order : {cudaq::solvers::first_order_trotter, + cudaq::solvers::second_order_trotter, + cudaq::solvers::fourth_order_trotter}) { + auto expected = + simulate_trotter_statevector(terms, time, steps, order, ket); + auto actual = + cudaq::get_state(apply_trotter_test_kernel{}, terms.num_qubits, + terms.coefficients, terms.words, time, steps, order); + + for (std::size_t i = 0; i < expected.size(); ++i) + expect_close(actual.amplitude(basis_bits(i, terms.num_qubits)), + expected[i]); + } +} + +TEST(TrotterTest, ApplyTrotterHandlesFourQubitHamiltonianWithManyTerms) { + cudaq::spin_op hamiltonian = 0.11 * cudaq::spin_op::from_word("XIII") - + 0.17 * cudaq::spin_op::from_word("IYII") + + 0.23 * cudaq::spin_op::from_word("IIZI") - + 0.29 * cudaq::spin_op::from_word("IIIX") + + 0.31 * cudaq::spin_op::from_word("XXII") + + 0.37 * cudaq::spin_op::from_word("IYZI") - + 0.41 * cudaq::spin_op::from_word("ZIIX") + + 0.43 * cudaq::spin_op::from_word("XIYZ") - + 0.47 * cudaq::spin_op::from_word("YYXI") + + 0.53 * cudaq::spin_op::from_word("ZXYZ"); + auto terms = cudaq::solvers::make_trotter_terms(hamiltonian); + + ASSERT_EQ(terms.num_qubits, 4); + ASSERT_GT(terms.coefficients.size(), 8); + ASSERT_EQ(terms.coefficients.size(), terms.words.size()); + + std::vector> ket(16, {0.0, 0.0}); + ket[0] = {1.0, 0.0}; + + const double time = 0.37; + const std::size_t steps = 2; + const int order = cudaq::solvers::second_order_trotter; + + auto expected = simulate_trotter_statevector(terms, time, steps, order, ket); + auto actual = + cudaq::get_state(apply_trotter_test_kernel{}, terms.num_qubits, + terms.coefficients, terms.words, time, steps, order); + + for (std::size_t i = 0; i < expected.size(); ++i) + expect_close(actual.amplitude(basis_bits(i, terms.num_qubits)), + expected[i]); +} + +TEST(TrotterTest, ApplyTrotterTreatsIdentityOnlyHamiltonianAsNoOp) { + auto terms = + cudaq::solvers::make_trotter_terms(1.5 * cudaq::spin_op::from_word("II")); + std::vector> ket{ + {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}}; + + auto actual = cudaq::get_state(apply_trotter_test_kernel{}, terms.num_qubits, + terms.coefficients, terms.words, 0.25, 4, + cudaq::solvers::second_order_trotter); + + for (std::size_t i = 0; i < ket.size(); ++i) + expect_close(actual.amplitude(basis_bits(i, terms.num_qubits)), ket[i]); +}