From ebb2ed9568536c8fbbdb4b181df97b662c30b8f7 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 20 Jul 2025 00:20:00 +0200 Subject: [PATCH 01/19] separated API trotterisation from operations which has the below benefits: - the remaining functions in operations.cpp are precise and do not need to be user-configured for accuracy (i.e. no Trotter hyperparameters) - the remaining functions in operations.cpp merely call the backend and do not include any bespoke logic (i.e. Trotter circuit scheduling) - incoming new Trotter functions for dynamical simulation will be more clearly delineated from the "standard" (and relatively boring) operations - the Trotter logic is isolated in preparation for it becoming more substantial with the introduction of randomisation, commuting groups, and that necessary for Linblad master equation solving --- quest/include/operations.h | 259 -------------------------- quest/include/quest.h | 1 + quest/include/trotterisation.h | 309 +++++++++++++++++++++++++++++++ quest/src/api/CMakeLists.txt | 1 + quest/src/api/operations.cpp | 137 -------------- quest/src/api/trotterisation.cpp | 175 +++++++++++++++++ tests/unit/CMakeLists.txt | 1 + tests/unit/operations.cpp | 17 -- tests/unit/trotterisation.cpp | 36 ++++ 9 files changed, 523 insertions(+), 413 deletions(-) create mode 100644 quest/include/trotterisation.h create mode 100644 quest/src/api/trotterisation.cpp create mode 100644 tests/unit/trotterisation.cpp diff --git a/quest/include/operations.h b/quest/include/operations.h index 0fb06560f..8941fc312 100644 --- a/quest/include/operations.h +++ b/quest/include/operations.h @@ -2461,270 +2461,11 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace); void postMultiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace); -/** @notyettested - * - * Effects (an approximation to) the exponential of @p sum, weighted by @p angle, upon @p qureg, - * via the symmetrized Trotter-Suzuki decomposition (arXiv). - * Increasing @p reps (the number of Trotter repetitions) or @p order (an even, positive integer or one) - * improves the accuracy of the approximation (reducing the "Trotter error" due to non-commuting - * terms of @p sum), though increases the runtime linearly and exponentially respectively. - * - * @formulae - * - * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of - * @f[ - \exp \left(\iu \, \theta \, \hat{H} \right) - * @f] - * via a Trotter-Suzuki decomposition of the specified @p order and number of repetitions (@p reps). - * Simulation is exact, regardless of @p order or @p reps, only when all terms in @p sum commute. - * - * @important - * Note that @f$ \theta @f$ lacks the @f$ -\frac{1}{2} @f$ prefactor present in other functions like - * applyPauliGadget(). - * - * To be precise, let @f$ r = @f$ @p reps and assume @p sum is composed of - * @f$ T @f$-many terms of the form - * @f[ - \hat{H} = \sum\limits_j^T c_j \, \hat{\sigma}_j - * @f] - * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$. - * - * - When @p order=1, this function performs first-order Trotterisation, whereby - * @f[ - \exp(\iu \, \theta \, \hat{H} ) - \approx - \prod\limits^{r} - \prod\limits_{j=1}^{T} - \exp \left( \iu \, \frac{\theta \, c_j}{r} \, \hat\sigma_j \right). - * @f] - * - When @p order=2, this function performs the lowest order "symmetrized" Suzuki decomposition, whereby - * @f[ - \exp(\iu \, \theta \, \hat{H} ) - \approx - \prod\limits^{r} \left[ - \prod\limits_{j=1}^{T} \exp \left( \iu \frac{\theta \, c_j}{2 \, r} \hat\sigma_j \right) - \prod\limits_{j=T}^{1} \exp \left( \iu \frac{\theta \, c_j}{2 \, r} \hat\sigma_j \right) - \right]. - * @f] - * - Greater, even values of @p order (denoted by symbol @f$ n @f$) invoke higher-order symmetrized decompositions - * @f$ S[\theta,n,r] @f$. Letting @f$ p = \left( 4 - 4^{1/(n-1)} \right)^{-1} @f$, these satisfy - * @f{align*} - S[\theta, n, 1] &= - \left( \prod\limits^2 S[p \, \theta, n-2, 1] \right) - S[ (1-4p)\,\theta, n-2, 1] - \left( \prod\limits^2 S[p \, \theta, n-2, 1] \right), - \\ - S[\theta, n, r] &= - \prod\limits^{r} S\left[\frac{\theta}{r}, n, 1\right]. - * @f} - * - * > These formulations are taken from 'Finding Exponential Product Formulas - * > of Higher Orders', Naomichi Hatano and Masuo Suzuki (2005) (arXiv). - * - * @equivalences - * - * - Time evolution of duration @f$ t @f$ under a time-independent Hamiltonian @p sum = @f$ \hat{H} @f$, as - * per the unitary time evolution operator - * @f[ - \hat{U}(t) = \exp(- \iu \, t \,\hat{H} \, / \, \hbar) - * @f] - * is approximated via @f$ \theta = - t / \hbar @f$. - * ``` - qreal time = 3.14; - qreal angle = - time / hbar; - applyTrotterizedPauliStrSumGadget(qureg, sum, angle, order, reps); - * ``` - * - This function is equivalent to applyNonUnitaryTrotterizedPauliStrSumGadget() when passing - * a @p qcomp instance with a zero imaginary component as the @p angle parameter. This latter - * function is useful for generalising dynamical simulation to imaginary-time evolution. - * - * @constraints - * - Unitarity of the prescribed exponential(s) requires that @p sum is Hermitian, ergo containing - * only real coefficients. Validation will check that @p sum is approximately Hermitian, permitting - * coefficients with imaginary components smaller (in magnitude) than epsilon. - * @f[ - \max\limits_{i} \Big|c_i| \le \valeps - * @f] - * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). - * Otherwise, use applyNonUnitaryTrotterizedPauliStrSumGadget() to permit non-Hermitian @p sum - * and ergo effect a non-unitary exponential(s). - * - The @p angle parameter is necessarily real despite the validation epsilon, but can be relaxed - * to an arbitrary complex scalar using applyNonUnitaryTrotterizedPauliStrSumGadget(). - * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly - * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. - * - * @param[in,out] qureg the state to modify. - * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. - * @param[in] angle an effective prefactor of @p sum in the exponent. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions - * - * @throws @validationerror - * - if @p qureg or @p sum are uninitialised. - * - if @p sum is not approximately Hermitian. - * - if @p sum contains non-identities on qubits beyond the size of @p qureg. - * - if @p order is not 1 nor a positive, @b even integer. - * - if @p reps is not a positive integer. - * - * @see - * - applyPauliGadget() - * - applyNonUnitaryTrotterizedPauliStrSumGadget() - * - * @author Tyson Jones - */ -void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps); - - -/// @notyetdoced -/// @notyettested -/// @see -/// - applyTrotterizedPauliStrSumGadget() -/// - applyControlledCompMatr1() -void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps); - - -/// @notyetdoced -/// @notyettested -/// @see -/// - applyTrotterizedPauliStrSumGadget() -/// - applyMultiControlledCompMatr1() -void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps); - - -/// @notyetdoced -/// @notyettested -/// @see -/// - applyTrotterizedPauliStrSumGadget() -/// - applyMultiStateControlledCompMatr1() -void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps); - - -/** @notyettested - * - * A generalisation of applyTrotterizedPauliStrSumGadget() which accepts a complex angle and permits - * @p sum to be non-Hermitian, thereby effecting a potentially non-unitary and non-CPTP operation. - * - * @formulae - * - * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of - * @f[ - \exp \left(\iu \, \theta \, \hat{H} \right) - * @f] - * via a Trotter-Suzuki decomposition of the specified @p order and number of repetitions (@p reps). - * - * See applyTrotterizedPauliStrSumGadget() for more information about the decomposition. - * - * @equivalences - * - * - When @p angle is set to @f$ \theta = \iu \, \tau @f$ and @p sum = @f$ \hat{H} @f$ is Hermitian, - * this function (approximately) evolves @p qureg in imaginary-time. That is, letting - * @f$ \hat{U}(t) = \exp(-\iu \, t \, \hat{H}) @f$ be the normalised unitary evolution operator, this - * function effects the imaginary-time operator - @f[ - \hat{V}(\tau) = \hat{U}(t=-\iu \tau) = \exp(- \tau \hat{H}). - * @f] - * This operation drives the system toward the (unnormalised) groundstate. - * Let @f$ \{ \ket{\phi_i} \} @f$ and @f$ \{ \ket{\lambda_i} \} @f$ be the eigenstates and respective - * eigenvalues of @f$ \hat{H} @f$, which are real due to Hermiticity. - * @f[ - \hat{H} = \sum \limits_i \lambda_i \ket{\phi_i}\bra{\phi_i}, - \;\;\;\;\; \lambda_i \in \mathbb{R}. - * @f] - * - * - When @p qureg is a statevector @f$ \svpsi @f$ and can ergo be expressed in the basis of - * @f$ \{ \ket{\phi_i} \} @f$ as @f$ \svpsi = \sum_i \alpha_i \ket{\phi_i} @f$, - * this function approximates - * @f[ - \svpsi \, \rightarrow \, \hat{V}(\tau) \svpsi = - \sum\limits_i \alpha_i \exp(- \tau \, \lambda_i) \ket{\phi_i}. - * @f] - * - When @p qureg is a density matrix and is ergo expressible as - * @f$ \dmrho = \sum\limits_{ij} \alpha_{ij} \ket{\phi_i}\bra{\phi_j} @f$, this function effects - * @f[ - \dmrho \, \rightarrow \, \hat{V}(\tau) \dmrho \hat{V}(\tau)^\dagger = - \sum\limits_{ij} \alpha_{ij} \exp(-\tau (\lambda_i + \lambda_j)) \ket{\phi_i}\bra{\phi_j}. - * @f] - * - * As @f$ \tau \rightarrow \infty @f$, the resulting unnormalised state approaches statevector - * @f$ \svpsi \rightarrow \alpha_0 \exp(-\tau \lambda_0) \ket{\phi_0} @f$ or density matrix - * @f$ \dmrho \rightarrow \alpha_{0,0} \exp(-2 \tau \lambda_0) \ket{\phi_0}\bra{\phi_0} @f$, - * where @f$ \lambda_0 @f$ is the minimum eigenvalue and @f$ \ket{\phi_0} @f$ is the groundstate. - * Assuming the initial overlap @f$ \alpha_0 @f$ is not zero (or exponentially tiny), - * subsequent renormalisation via setQuregToRenormalized() produces the pure - * ground-state @f$ \ket{\phi_0} @f$. - * - * ``` - // pray for a non-zero initial overlap - initRandomPureState(qureg); // works even for density matrices - - // minimize then renormalise - qreal tau = 10; // impatient infinity - int order = 4; - int reps = 100; - applyNonUnitaryTrotterizedPauliStrSumGadget(qureg, hamil, tau * 1i, order, reps); - setQuregToRenormalized(qureg); - - // ground-state (phi_0) - reportQureg(qureg); - - // lowest lying eigenvalue (lambda_0) - qreal expec = calcExpecPauliStrSum(qureg, hamil); - reportScalar("expec", expec); - * ``` - * - * Note degenerate eigenvalues will yield a pure superposition of the corresponding eigenstates, with - * coefficients informed by the initial, relative populations. - * - * - When @p angle is real and @p sum is Hermitian (has approximately real coefficients), this - * function is equivalent to applyTrotterizedPauliStrSumGadget() - * - * @constraints - * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly - * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. - * - * @param[in,out] qureg the state to modify. - * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. - * @param[in] angle an effective prefactor of @p sum in the exponent. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions - * - * @throws @validationerror - * - if @p qureg or @p sum are uninitialised. - * - if @p sum contains non-identities on qubits beyond the size of @p qureg. - * - if @p order is not 1 nor a positive, @b even integer. - * - if @p reps is not a positive integer. - * - * @author Tyson Jones - */ -void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps); - - // end de-mangler #ifdef __cplusplus } #endif -#ifdef __cplusplus - - -/// @notyettested -/// @notyetvalidated -/// @notyetdoced -/// @cppvectoroverload -/// @see applyMultiControlledTrotterizedPauliStrSumGadget() -void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector controls, PauliStrSum sum, qreal angle, int order, int reps); - - -/// @notyettested -/// @notyetvalidated -/// @notyetdoced -/// @cppvectoroverload -/// @see applyMultiStateControlledTrotterizedPauliStrSumGadget() -void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector controls, std::vector states, PauliStrSum sum, qreal angle, int order, int reps); - - -#endif // __cplusplus - /** @} */ diff --git a/quest/include/quest.h b/quest/include/quest.h index fcc49ab76..58901dcca 100644 --- a/quest/include/quest.h +++ b/quest/include/quest.h @@ -48,6 +48,7 @@ #include "quest/include/debug.h" #include "quest/include/decoherence.h" #include "quest/include/environment.h" +#include "quest/include/trotterisation.h" #include "quest/include/initialisations.h" #include "quest/include/channels.h" #include "quest/include/operations.h" diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h new file mode 100644 index 000000000..6c16e3d9a --- /dev/null +++ b/quest/include/trotterisation.h @@ -0,0 +1,309 @@ +/** @file + * API signatures for effecting Trotterised operators which + * approximate the action of exponentials of PauliStrSum + * + * @author Tyson Jones + * + * @defgroup trotterisation Trotterisation + * @ingroup api + * @brief Functions for Trottersing operations upon Quregs. + * @{ + */ + +#ifndef TROTTERISATION_H +#define TROTTERISATION_H + +#include "quest/include/qureg.h" +#include "quest/include/paulis.h" +#include "quest/include/matrices.h" + +#ifdef __cplusplus + #include +#endif + + + +/** + * @defgroup trotter_paulistrsum PauliStrSum gadgets + * @brief Functions for using Trotterisation to approximate the action of + * exponentials of weighted sums of Pauli tensors upon Quregs. + * @{ + */ + + +#ifdef __cplusplus +extern "C" { +#endif + + +/** @notyettested + * + * Effects (an approximation to) the exponential of @p sum, weighted by @p angle, upon @p qureg, + * via the symmetrized Trotter-Suzuki decomposition (arXiv). + * Increasing @p reps (the number of Trotter repetitions) or @p order (an even, positive integer or one) + * improves the accuracy of the approximation (reducing the "Trotter error" due to non-commuting + * terms of @p sum), though increases the runtime linearly and exponentially respectively. + * + * @formulae + * + * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of + * @f[ + \exp \left(\iu \, \theta \, \hat{H} \right) + * @f] + * via a Trotter-Suzuki decomposition of the specified @p order and number of repetitions (@p reps). + * Simulation is exact, regardless of @p order or @p reps, only when all terms in @p sum commute. + * + * @important + * Note that @f$ \theta @f$ lacks the @f$ -\frac{1}{2} @f$ prefactor present in other functions like + * applyPauliGadget(). + * + * To be precise, let @f$ r = @f$ @p reps and assume @p sum is composed of + * @f$ T @f$-many terms of the form + * @f[ + \hat{H} = \sum\limits_j^T c_j \, \hat{\sigma}_j + * @f] + * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$. + * + * - When @p order=1, this function performs first-order Trotterisation, whereby + * @f[ + \exp(\iu \, \theta \, \hat{H} ) + \approx + \prod\limits^{r} + \prod\limits_{j=1}^{T} + \exp \left( \iu \, \frac{\theta \, c_j}{r} \, \hat\sigma_j \right). + * @f] + * - When @p order=2, this function performs the lowest order "symmetrized" Suzuki decomposition, whereby + * @f[ + \exp(\iu \, \theta \, \hat{H} ) + \approx + \prod\limits^{r} \left[ + \prod\limits_{j=1}^{T} \exp \left( \iu \frac{\theta \, c_j}{2 \, r} \hat\sigma_j \right) + \prod\limits_{j=T}^{1} \exp \left( \iu \frac{\theta \, c_j}{2 \, r} \hat\sigma_j \right) + \right]. + * @f] + * - Greater, even values of @p order (denoted by symbol @f$ n @f$) invoke higher-order symmetrized decompositions + * @f$ S[\theta,n,r] @f$. Letting @f$ p = \left( 4 - 4^{1/(n-1)} \right)^{-1} @f$, these satisfy + * @f{align*} + S[\theta, n, 1] &= + \left( \prod\limits^2 S[p \, \theta, n-2, 1] \right) + S[ (1-4p)\,\theta, n-2, 1] + \left( \prod\limits^2 S[p \, \theta, n-2, 1] \right), + \\ + S[\theta, n, r] &= + \prod\limits^{r} S\left[\frac{\theta}{r}, n, 1\right]. + * @f} + * + * > These formulations are taken from 'Finding Exponential Product Formulas + * > of Higher Orders', Naomichi Hatano and Masuo Suzuki (2005) (arXiv). + * + * @equivalences + * + * - Time evolution of duration @f$ t @f$ under a time-independent Hamiltonian @p sum = @f$ \hat{H} @f$, as + * per the unitary time evolution operator + * @f[ + \hat{U}(t) = \exp(- \iu \, t \,\hat{H} \, / \, \hbar) + * @f] + * is approximated via @f$ \theta = - t / \hbar @f$. + * ``` + qreal time = 3.14; + qreal angle = - time / hbar; + applyTrotterizedPauliStrSumGadget(qureg, sum, angle, order, reps); + * ``` + * - This function is equivalent to applyNonUnitaryTrotterizedPauliStrSumGadget() when passing + * a @p qcomp instance with a zero imaginary component as the @p angle parameter. This latter + * function is useful for generalising dynamical simulation to imaginary-time evolution. + * + * @constraints + * - Unitarity of the prescribed exponential(s) requires that @p sum is Hermitian, ergo containing + * only real coefficients. Validation will check that @p sum is approximately Hermitian, permitting + * coefficients with imaginary components smaller (in magnitude) than epsilon. + * @f[ + \max\limits_{i} \Big|c_i| \le \valeps + * @f] + * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). + * Otherwise, use applyNonUnitaryTrotterizedPauliStrSumGadget() to permit non-Hermitian @p sum + * and ergo effect a non-unitary exponential(s). + * - The @p angle parameter is necessarily real despite the validation epsilon, but can be relaxed + * to an arbitrary complex scalar using applyNonUnitaryTrotterizedPauliStrSumGadget(). + * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly + * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. + * + * @param[in,out] qureg the state to modify. + * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. + * @param[in] angle an effective prefactor of @p sum in the exponent. + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) + * @param[in] reps the number of Trotter repetitions + * + * @throws @validationerror + * - if @p qureg or @p sum are uninitialised. + * - if @p sum is not approximately Hermitian. + * - if @p sum contains non-identities on qubits beyond the size of @p qureg. + * - if @p order is not 1 nor a positive, @b even integer. + * - if @p reps is not a positive integer. + * + * @see + * - applyPauliGadget() + * - applyNonUnitaryTrotterizedPauliStrSumGadget() + * + * @author Tyson Jones + */ +void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps); + + +/// @notyetdoced +/// @notyettested +/// @see +/// - applyTrotterizedPauliStrSumGadget() +/// - applyControlledCompMatr1() +void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps); + + +/// @notyetdoced +/// @notyettested +/// @see +/// - applyTrotterizedPauliStrSumGadget() +/// - applyMultiControlledCompMatr1() +void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps); + + +/// @notyetdoced +/// @notyettested +/// @see +/// - applyTrotterizedPauliStrSumGadget() +/// - applyMultiStateControlledCompMatr1() +void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps); + + +/** @notyettested + * + * A generalisation of applyTrotterizedPauliStrSumGadget() which accepts a complex angle and permits + * @p sum to be non-Hermitian, thereby effecting a potentially non-unitary and non-CPTP operation. + * + * @formulae + * + * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of + * @f[ + \exp \left(\iu \, \theta \, \hat{H} \right) + * @f] + * via a Trotter-Suzuki decomposition of the specified @p order and number of repetitions (@p reps). + * + * See applyTrotterizedPauliStrSumGadget() for more information about the decomposition. + * + * @equivalences + * + * - When @p angle is set to @f$ \theta = \iu \, \tau @f$ and @p sum = @f$ \hat{H} @f$ is Hermitian, + * this function (approximately) evolves @p qureg in imaginary-time. That is, letting + * @f$ \hat{U}(t) = \exp(-\iu \, t \, \hat{H}) @f$ be the normalised unitary evolution operator, this + * function effects the imaginary-time operator + @f[ + \hat{V}(\tau) = \hat{U}(t=-\iu \tau) = \exp(- \tau \hat{H}). + * @f] + * This operation drives the system toward the (unnormalised) groundstate. + * Let @f$ \{ \ket{\phi_i} \} @f$ and @f$ \{ \ket{\lambda_i} \} @f$ be the eigenstates and respective + * eigenvalues of @f$ \hat{H} @f$, which are real due to Hermiticity. + * @f[ + \hat{H} = \sum \limits_i \lambda_i \ket{\phi_i}\bra{\phi_i}, + \;\;\;\;\; \lambda_i \in \mathbb{R}. + * @f] + * + * - When @p qureg is a statevector @f$ \svpsi @f$ and can ergo be expressed in the basis of + * @f$ \{ \ket{\phi_i} \} @f$ as @f$ \svpsi = \sum_i \alpha_i \ket{\phi_i} @f$, + * this function approximates + * @f[ + \svpsi \, \rightarrow \, \hat{V}(\tau) \svpsi = + \sum\limits_i \alpha_i \exp(- \tau \, \lambda_i) \ket{\phi_i}. + * @f] + * - When @p qureg is a density matrix and is ergo expressible as + * @f$ \dmrho = \sum\limits_{ij} \alpha_{ij} \ket{\phi_i}\bra{\phi_j} @f$, this function effects + * @f[ + \dmrho \, \rightarrow \, \hat{V}(\tau) \dmrho \hat{V}(\tau)^\dagger = + \sum\limits_{ij} \alpha_{ij} \exp(-\tau (\lambda_i + \lambda_j)) \ket{\phi_i}\bra{\phi_j}. + * @f] + * + * As @f$ \tau \rightarrow \infty @f$, the resulting unnormalised state approaches statevector + * @f$ \svpsi \rightarrow \alpha_0 \exp(-\tau \lambda_0) \ket{\phi_0} @f$ or density matrix + * @f$ \dmrho \rightarrow \alpha_{0,0} \exp(-2 \tau \lambda_0) \ket{\phi_0}\bra{\phi_0} @f$, + * where @f$ \lambda_0 @f$ is the minimum eigenvalue and @f$ \ket{\phi_0} @f$ is the groundstate. + * Assuming the initial overlap @f$ \alpha_0 @f$ is not zero (or exponentially tiny), + * subsequent renormalisation via setQuregToRenormalized() produces the pure + * ground-state @f$ \ket{\phi_0} @f$. + * + * ``` + // pray for a non-zero initial overlap + initRandomPureState(qureg); // works even for density matrices + + // minimize then renormalise + qreal tau = 10; // impatient infinity + int order = 4; + int reps = 100; + applyNonUnitaryTrotterizedPauliStrSumGadget(qureg, hamil, tau * 1i, order, reps); + setQuregToRenormalized(qureg); + + // ground-state (phi_0) + reportQureg(qureg); + + // lowest lying eigenvalue (lambda_0) + qreal expec = calcExpecPauliStrSum(qureg, hamil); + reportScalar("expec", expec); + * ``` + * + * Note degenerate eigenvalues will yield a pure superposition of the corresponding eigenstates, with + * coefficients informed by the initial, relative populations. + * + * - When @p angle is real and @p sum is Hermitian (has approximately real coefficients), this + * function is equivalent to applyTrotterizedPauliStrSumGadget() + * + * @constraints + * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly + * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. + * + * @param[in,out] qureg the state to modify. + * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. + * @param[in] angle an effective prefactor of @p sum in the exponent. + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) + * @param[in] reps the number of Trotter repetitions + * + * @throws @validationerror + * - if @p qureg or @p sum are uninitialised. + * - if @p sum contains non-identities on qubits beyond the size of @p qureg. + * - if @p order is not 1 nor a positive, @b even integer. + * - if @p reps is not a positive integer. + * + * @author Tyson Jones + */ +void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps); + + +// end de-mangler +#ifdef __cplusplus +} +#endif + +#ifdef __cplusplus + + +/// @notyettested +/// @notyetvalidated +/// @notyetdoced +/// @cppvectoroverload +/// @see applyMultiControlledTrotterizedPauliStrSumGadget() +void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector controls, PauliStrSum sum, qreal angle, int order, int reps); + + +/// @notyettested +/// @notyetvalidated +/// @notyetdoced +/// @cppvectoroverload +/// @see applyMultiStateControlledTrotterizedPauliStrSumGadget() +void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vector controls, std::vector states, PauliStrSum sum, qreal angle, int order, int reps); + + +#endif // __cplusplus + +/** @} */ + + + +#endif // TROTTERISATION_H + +/** @} */ // (end file-wide doxygen defgroup) diff --git a/quest/src/api/CMakeLists.txt b/quest/src/api/CMakeLists.txt index a6d3137b6..afb298df5 100644 --- a/quest/src/api/CMakeLists.txt +++ b/quest/src/api/CMakeLists.txt @@ -11,5 +11,6 @@ target_sources(QuEST operations.cpp paulis.cpp qureg.cpp + trotterisation.cpp types.cpp ) \ No newline at end of file diff --git a/quest/src/api/operations.cpp b/quest/src/api/operations.cpp index 1b57f7d2a..0ec30c72d 100644 --- a/quest/src/api/operations.cpp +++ b/quest/src/api/operations.cpp @@ -1376,145 +1376,8 @@ void postMultiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace) { // workspace -> qureg, and qureg -> sum * qureg } -void internal_applyFirstOrderTrotterRepetition( - Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, bool reverse -) { - // apply each sum term as a gadget, in forward or reverse order - for (qindex i=0; i -> exp(i angle * sum)|psi> - qcomp arg = angle * coeff; - localiser_statevec_anyCtrlPauliGadget(qureg, ketCtrls, states, str, arg); - - if (!qureg.isDensityMatrix) - continue; - - // effect rho -> rho dagger(i angle * sum) - arg *= paulis_hasOddNumY(str) ? 1 : -1; - str = paulis_getShiftedPauliStr(str, qureg.numQubits); - localiser_statevec_anyCtrlPauliGadget(qureg, braCtrls, states, str, arg); - } -} - -void internal_applyHigherOrderTrotterRepetition( - Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, int order -) { - if (order == 1) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, false); - - } else if (order == 2) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, false); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, true); - - } else { - qreal p = 1. / (4 - std::pow(4, 1./(order-1))); - qcomp a = p * angle; - qcomp b = (1-4*p) * angle; - - int lower = order - 2; - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - } -} - -void internal_applyAllTrotterRepetitions( - Qureg qureg, int* controls, int* states, int numControls, - PauliStrSum sum, qcomp angle, int order, int reps -) { - // exp(i angle sum) = identity when angle=0 - if (angle == qcomp(0,0)) - return; - - // prepare control-qubit lists once for all invoked gadgets below - auto ketCtrlsVec = util_getVector(controls, numControls); - auto braCtrlsVec = (qureg.isDensityMatrix)? util_getBraQubits(ketCtrlsVec, qureg) : vector{}; - auto statesVec = util_getVector(states, numControls); - - qcomp arg = angle / reps; - - // perform carefully-ordered sequence of gadgets - for (int r=0; r controls, PauliStrSum sum, qreal angle, int order, int reps) { - - applyMultiControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), controls.size(), sum, angle, order, reps); -} - -void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, vector controls, vector states, PauliStrSum sum, qreal angle, int order, int reps) { - validate_controlsMatchStates(controls.size(), states.size(), __func__); - - applyMultiStateControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), states.data(), controls.size(), sum, angle, order, reps); -} - /* diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp new file mode 100644 index 000000000..43d371059 --- /dev/null +++ b/quest/src/api/trotterisation.cpp @@ -0,0 +1,175 @@ +/** @file + * API definitions for functions which involve Trotterising + * exponential operators, such as PauliStrSum gadgets, and + * so are inherently approximate. + * + * @author Tyson Jones + */ + +#include "quest/include/qureg.h" +#include "quest/include/paulis.h" +#include "quest/include/matrices.h" + +#include "quest/src/core/validation.hpp" +#include "quest/src/core/utilities.hpp" +#include "quest/src/core/localiser.hpp" + +#include + +using std::vector; + + + +/* + * INTERNAL UTILS + */ + +extern bool paulis_hasOddNumY(PauliStr str); +extern PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift); + +void internal_applyFirstOrderTrotterRepetition( + Qureg qureg, vector& ketCtrls, vector& braCtrls, + vector& states, PauliStrSum sum, qcomp angle, bool reverse +) { + // apply each sum term as a gadget, in forward or reverse order + for (qindex i=0; i -> exp(i angle * sum)|psi> + qcomp arg = angle * coeff; + localiser_statevec_anyCtrlPauliGadget(qureg, ketCtrls, states, str, arg); + + if (!qureg.isDensityMatrix) + continue; + + // effect rho -> rho dagger(i angle * sum) + arg *= paulis_hasOddNumY(str) ? 1 : -1; + str = paulis_getShiftedPauliStr(str, qureg.numQubits); + localiser_statevec_anyCtrlPauliGadget(qureg, braCtrls, states, str, arg); + } +} + +void internal_applyHigherOrderTrotterRepetition( + Qureg qureg, vector& ketCtrls, vector& braCtrls, + vector& states, PauliStrSum sum, qcomp angle, int order +) { + if (order == 1) { + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, false); + + } else if (order == 2) { + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, true); + + } else { + qreal p = 1. / (4 - std::pow(4, 1./(order-1))); + qcomp a = p * angle; + qcomp b = (1-4*p) * angle; + + int lower = order - 2; + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + } +} + +void internal_applyAllTrotterRepetitions( + Qureg qureg, int* controls, int* states, int numControls, + PauliStrSum sum, qcomp angle, int order, int reps +) { + // exp(i angle sum) = identity when angle=0 + if (angle == qcomp(0,0)) + return; + + // prepare control-qubit lists once for all invoked gadgets below + auto ketCtrlsVec = util_getVector(controls, numControls); + auto braCtrlsVec = (qureg.isDensityMatrix)? util_getBraQubits(ketCtrlsVec, qureg) : vector{}; + auto statesVec = util_getVector(states, numControls); + + qcomp arg = angle / reps; + + // perform carefully-ordered sequence of gadgets + for (int r=0; r controls, PauliStrSum sum, qreal angle, int order, int reps) { + + applyMultiControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), controls.size(), sum, angle, order, reps); +} + +void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, vector controls, vector states, PauliStrSum sum, qreal angle, int order, int reps) { + validate_controlsMatchStates(controls.size(), states.size(), __func__); + + applyMultiStateControlledTrotterizedPauliStrSumGadget(qureg, controls.data(), states.data(), controls.size(), sum, angle, order, reps); +} diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 033751832..76642be01 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -12,5 +12,6 @@ target_sources(tests operations.cpp paulis.cpp qureg.cpp + trotterisation.cpp types.cpp ) \ No newline at end of file diff --git a/tests/unit/operations.cpp b/tests/unit/operations.cpp index b8522321a..7c6c39096 100644 --- a/tests/unit/operations.cpp +++ b/tests/unit/operations.cpp @@ -2117,20 +2117,3 @@ TEST_CASE( "applyNonUnitaryPauliGadget", TEST_CATEGORY ) { /** @} (end defgroup) */ - - - -/** - * @todo - * UNTESTED FUNCTIONS - */ - -void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps); - -void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps); - -void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps); - -void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps); - -void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps); diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp new file mode 100644 index 000000000..174b7c66e --- /dev/null +++ b/tests/unit/trotterisation.cpp @@ -0,0 +1,36 @@ +/** @file + * Unit tests of the trotterisation module. + * + * @author Tyson Jones + * + * @defgroup unittrotter Trotterisation + * @ingroup unittests + */ + +#include "quest.h" + + + +/* + * UTILITIES + */ + +#define TEST_CATEGORY \ + LABEL_UNIT_TAG "[trotterisation]" + + + +/** + * @todo + * UNTESTED FUNCTIONS + */ + +void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps); + +void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps); + +void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps); + +void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps); + +void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps); From 5cf77b5acca9087963da227d0337fe52ad9b17bc Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 20 Jul 2025 00:59:14 +0200 Subject: [PATCH 02/19] added closed evolution merely aliasing the existing gadget Trotterisation --- quest/include/operations.h | 2 +- quest/include/trotterisation.h | 34 ++++++++++++++++++++++++++++++++ quest/src/api/trotterisation.cpp | 29 +++++++++++++++++++++++++++ 3 files changed, 64 insertions(+), 1 deletion(-) diff --git a/quest/include/operations.h b/quest/include/operations.h index 8941fc312..c6ef231d9 100644 --- a/quest/include/operations.h +++ b/quest/include/operations.h @@ -116,7 +116,7 @@ void multiplyCompMatr1(Qureg qureg, int target, CompMatr1 matrix); * * In general, this function will break the normalisation of @p qureg and result in a * non-physical state, and is useful for preparing sub-expressions of formulae like - * the Linbladian. + * the Lindbladian. * * @myexample * ``` diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 6c16e3d9a..a4aae0756 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -304,6 +304,40 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vec + +/** + * @defgroup trotter_timeevol Time evolution + * @brief Functions for approximate dynamical simulation. + * @{ + */ + + +#ifdef __cplusplus +extern "C" { +#endif + + +/// @notyettested +/// @notyetvalidated +/// @notyetdoced +void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamiltonian, qreal time, int order, int reps); + + +/// @notyettested +/// @notyetvalidated +/// @notyetdoced +void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamiltonian, qreal tau, int order, int reps); + + +// end de-mangler +#ifdef __cplusplus +} +#endif + +/** @} */ + + + #endif // TROTTERISATION_H /** @} */ // (end file-wide doxygen defgroup) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 43d371059..cc1b6b586 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -173,3 +173,32 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, vector Date: Mon, 21 Jul 2025 20:04:14 +0200 Subject: [PATCH 03/19] added demo applyTrotterizedPauliNoisyTimeEvolution though it is unsure this should be the final interface; we can also accept jump operators given as projectors and diagonal matrices without much greater complexity, and SHOULD also deal with exponentiating general dense matrices. This commit also includes a Trotter patch we'll commit separately to dev for publicity --- quest/include/trotterisation.h | 10 +- quest/src/api/trotterisation.cpp | 155 +++++++++++++++++++++++++------ tests/unit/trotterisation.cpp | 8 ++ 3 files changed, 144 insertions(+), 29 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index a4aae0756..cda730395 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -320,13 +320,19 @@ extern "C" { /// @notyettested /// @notyetvalidated /// @notyetdoced -void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamiltonian, qreal time, int order, int reps); +void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps); /// @notyettested /// @notyetvalidated /// @notyetdoced -void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamiltonian, qreal tau, int order, int reps); +void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); + + +/// @notyettested +/// @notyetvalidated +/// @notyetdoced +void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); // end de-mangler diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index cc1b6b586..99ec32e8a 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -29,7 +29,7 @@ extern PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift); void internal_applyFirstOrderTrotterRepetition( Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, bool reverse + vector& states, PauliStrSum sum, qcomp angle, bool postmultiply, bool reverse ) { // apply each sum term as a gadget, in forward or reverse order for (qindex i=0; i rho dagger(i angle * sum) arg *= paulis_hasOddNumY(str) ? 1 : -1; str = paulis_getShiftedPauliStr(str, qureg.numQubits); @@ -53,14 +58,14 @@ void internal_applyFirstOrderTrotterRepetition( void internal_applyHigherOrderTrotterRepetition( Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, int order + vector& states, PauliStrSum sum, qcomp angle, int order, bool postmultiply ) { if (order == 1) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, postmultiply, false); } else if (order == 2) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, false); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, true); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, postmultiply, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, postmultiply, true); } else { qreal p = 1. / (4 - std::pow(4, 1./(order-1))); @@ -68,17 +73,17 @@ void internal_applyHigherOrderTrotterRepetition( qcomp b = (1-4*p) * angle; int lower = order - 2; - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); // angle -> a + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower, postmultiply); // angle -> b + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); } } void internal_applyAllTrotterRepetitions( Qureg qureg, int* controls, int* states, int numControls, - PauliStrSum sum, qcomp angle, int order, int reps + PauliStrSum sum, qcomp angle, int order, int reps, bool postmultiply ) { // exp(i angle sum) = identity when angle=0 if (angle == qcomp(0,0)) @@ -94,7 +99,7 @@ void internal_applyAllTrotterRepetitions( // perform carefully-ordered sequence of gadgets for (int r=0; r newStrings; + vector newCoeffs; + + // premature optimisation + qindex numNewTerms = 2 * hamil.numTerms + numJumps; + newStrings.reserve(numNewTerms); + newCoeffs.reserve(numNewTerms); + + // collect -i[H,rho] terms + for (qindex n=0; n Date: Wed, 23 Jul 2025 23:40:11 +0200 Subject: [PATCH 04/19] added doc for closed evolution --- quest/include/trotterisation.h | 350 ++++++++++++++++++++++++--------- 1 file changed, 260 insertions(+), 90 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index cda730395..74e420f2c 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -38,15 +38,16 @@ extern "C" { /** @notyettested * - * Effects (an approximation to) the exponential of @p sum, weighted by @p angle, upon @p qureg, + * Effects an approximation to the exponential of @p sum, weighted by @p angle times @f$ i @f$, upon @p qureg, * via the symmetrized Trotter-Suzuki decomposition (arXiv). * Increasing @p reps (the number of Trotter repetitions) or @p order (an even, positive integer or one) - * improves the accuracy of the approximation (reducing the "Trotter error" due to non-commuting - * terms of @p sum), though increases the runtime linearly and exponentially respectively. + * improves the accuracy of the approximation by reducing the "Trotter error" due to non-commuting + * terms of @p sum, though increases the runtime linearly and exponentially respectively. * * @formulae * - * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of + * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle @f$ \in \mathbb{R} @f$. This function approximates + * the action of * @f[ \exp \left(\iu \, \theta \, \hat{H} \right) * @f] @@ -54,7 +55,7 @@ extern "C" { * Simulation is exact, regardless of @p order or @p reps, only when all terms in @p sum commute. * * @important - * Note that @f$ \theta @f$ lacks the @f$ -\frac{1}{2} @f$ prefactor present in other functions like + * Observe that @f$ \theta @f$ lacks the @f$ -\frac{1}{2} @f$ prefactor present in other functions like * applyPauliGadget(). * * To be precise, let @f$ r = @f$ @p reps and assume @p sum is composed of @@ -64,7 +65,8 @@ extern "C" { * @f] * where @f$ c_j @f$ is the coefficient of the @f$ j @f$-th PauliStr @f$ \hat{\sigma}_j @f$. * - * - When @p order=1, this function performs first-order Trotterisation, whereby + * - When @p order=1, this function performs first-order Trotterisation, where the terms of @p sum + * are effected in a repeated, arbitrary but fixed order. * @f[ \exp(\iu \, \theta \, \hat{H} ) \approx @@ -72,7 +74,9 @@ extern "C" { \prod\limits_{j=1}^{T} \exp \left( \iu \, \frac{\theta \, c_j}{r} \, \hat\sigma_j \right). * @f] - * - When @p order=2, this function performs the lowest order "symmetrized" Suzuki decomposition, whereby + * + * - When @p order=2, this function performs the lowest order "symmetrized" Suzuki decomposition, whereby + * each repetition effects the terms of @p sum forward then in reverse. * @f[ \exp(\iu \, \theta \, \hat{H} ) \approx @@ -81,8 +85,11 @@ extern "C" { \prod\limits_{j=T}^{1} \exp \left( \iu \frac{\theta \, c_j}{2 \, r} \hat\sigma_j \right) \right]. * @f] + * * - Greater, even values of @p order (denoted by symbol @f$ n @f$) invoke higher-order symmetrized decompositions - * @f$ S[\theta,n,r] @f$. Letting @f$ p = \left( 4 - 4^{1/(n-1)} \right)^{-1} @f$, these satisfy + * @f$ S[\theta,n,r] @f$. These see the lower order Trotter circuits repeated twice forward, then reversed, then + * twice forward again, recursively. To be precise, letting @f$ p = \left( 4 - 4^{1/(n-1)} \right)^{-1} @f$, these + * satisfy * @f{align*} S[\theta, n, 1] &= \left( \prod\limits^2 S[p \, \theta, n-2, 1] \right) @@ -98,22 +105,20 @@ extern "C" { * * @equivalences * - * - Time evolution of duration @f$ t @f$ under a time-independent Hamiltonian @p sum = @f$ \hat{H} @f$, as - * per the unitary time evolution operator + * - By passing @f$ \theta = - \Delta t / \hbar @f$, this function approximates unitary time evolution of a closed + * system under the time-independent Hamiltonian @p sum = @f$ \hat{H} @f$ over a duration of @f$ \Delta t @f$, as + * described by propagator * @f[ - \hat{U}(t) = \exp(- \iu \, t \,\hat{H} \, / \, \hbar) + \hat{U}(\Delta t) = \exp(- \iu \, \Delta t \,\hat{H} \, / \, \hbar), * @f] - * is approximated via @f$ \theta = - t / \hbar @f$. - * ``` - qreal time = 3.14; - qreal angle = - time / hbar; - applyTrotterizedPauliStrSumGadget(qureg, sum, angle, order, reps); - * ``` + * as utilised by the function applyTrotterizedUnitaryTimeEvolution(). + * * - This function is equivalent to applyNonUnitaryTrotterizedPauliStrSumGadget() when passing * a @p qcomp instance with a zero imaginary component as the @p angle parameter. This latter * function is useful for generalising dynamical simulation to imaginary-time evolution. * * @constraints + * * - Unitarity of the prescribed exponential(s) requires that @p sum is Hermitian, ergo containing * only real coefficients. Validation will check that @p sum is approximately Hermitian, permitting * coefficients with imaginary components smaller (in magnitude) than epsilon. @@ -123,14 +128,18 @@ extern "C" { * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). * Otherwise, use applyNonUnitaryTrotterizedPauliStrSumGadget() to permit non-Hermitian @p sum * and ergo effect a non-unitary exponential(s). - * - The @p angle parameter is necessarily real despite the validation epsilon, but can be relaxed - * to an arbitrary complex scalar using applyNonUnitaryTrotterizedPauliStrSumGadget(). + * + * - The @p angle parameter is necessarily real to retain unitarity, but can be relaxed to an arbitrary + * complex scalar (i.e. a @p qcomp) using applyNonUnitaryTrotterizedPauliStrSumGadget(). This permits + * cancelling the complex unit @f$ i @f$ to effect non-unitary @f$ \exp(\theta \, \hat{H}) @f$ as + * is useful for imaginary-time evolution. + * * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly - * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. + * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute, or @p reps @f$ \rightarrow \infty @f$. * * @param[in,out] qureg the state to modify. * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. - * @param[in] angle an effective prefactor of @p sum in the exponent. + * @param[in] angle the prefactor of @p sum times @f$ i @f$ in the exponent. * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) * @param[in] reps the number of Trotter repetitions * @@ -144,6 +153,7 @@ extern "C" { * @see * - applyPauliGadget() * - applyNonUnitaryTrotterizedPauliStrSumGadget() + * - applyTrotterizedUnitaryTimeEvolution() * * @author Tyson Jones */ @@ -176,12 +186,13 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con /** @notyettested * - * A generalisation of applyTrotterizedPauliStrSumGadget() which accepts a complex angle and permits + * A generalisation of applyTrotterizedPauliStrSumGadget() which accepts a complex @p angle and permits * @p sum to be non-Hermitian, thereby effecting a potentially non-unitary and non-CPTP operation. * * @formulae * - * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle. This function approximates the action of + * Let @f$ \hat{H} = @f$ @p sum and @f$ \theta = @f$ @p angle @f$ \in \mathbb{C} @f$. This function + * approximates the action of * @f[ \exp \left(\iu \, \theta \, \hat{H} \right) * @f] @@ -191,69 +202,19 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con * * @equivalences * - * - When @p angle is set to @f$ \theta = \iu \, \tau @f$ and @p sum = @f$ \hat{H} @f$ is Hermitian, - * this function (approximately) evolves @p qureg in imaginary-time. That is, letting - * @f$ \hat{U}(t) = \exp(-\iu \, t \, \hat{H}) @f$ be the normalised unitary evolution operator, this - * function effects the imaginary-time operator - @f[ - \hat{V}(\tau) = \hat{U}(t=-\iu \tau) = \exp(- \tau \hat{H}). - * @f] - * This operation drives the system toward the (unnormalised) groundstate. - * Let @f$ \{ \ket{\phi_i} \} @f$ and @f$ \{ \ket{\lambda_i} \} @f$ be the eigenstates and respective - * eigenvalues of @f$ \hat{H} @f$, which are real due to Hermiticity. - * @f[ - \hat{H} = \sum \limits_i \lambda_i \ket{\phi_i}\bra{\phi_i}, - \;\;\;\;\; \lambda_i \in \mathbb{R}. - * @f] - * - * - When @p qureg is a statevector @f$ \svpsi @f$ and can ergo be expressed in the basis of - * @f$ \{ \ket{\phi_i} \} @f$ as @f$ \svpsi = \sum_i \alpha_i \ket{\phi_i} @f$, - * this function approximates - * @f[ - \svpsi \, \rightarrow \, \hat{V}(\tau) \svpsi = - \sum\limits_i \alpha_i \exp(- \tau \, \lambda_i) \ket{\phi_i}. - * @f] - * - When @p qureg is a density matrix and is ergo expressible as - * @f$ \dmrho = \sum\limits_{ij} \alpha_{ij} \ket{\phi_i}\bra{\phi_j} @f$, this function effects - * @f[ - \dmrho \, \rightarrow \, \hat{V}(\tau) \dmrho \hat{V}(\tau)^\dagger = - \sum\limits_{ij} \alpha_{ij} \exp(-\tau (\lambda_i + \lambda_j)) \ket{\phi_i}\bra{\phi_j}. - * @f] - * - * As @f$ \tau \rightarrow \infty @f$, the resulting unnormalised state approaches statevector - * @f$ \svpsi \rightarrow \alpha_0 \exp(-\tau \lambda_0) \ket{\phi_0} @f$ or density matrix - * @f$ \dmrho \rightarrow \alpha_{0,0} \exp(-2 \tau \lambda_0) \ket{\phi_0}\bra{\phi_0} @f$, - * where @f$ \lambda_0 @f$ is the minimum eigenvalue and @f$ \ket{\phi_0} @f$ is the groundstate. - * Assuming the initial overlap @f$ \alpha_0 @f$ is not zero (or exponentially tiny), - * subsequent renormalisation via setQuregToRenormalized() produces the pure - * ground-state @f$ \ket{\phi_0} @f$. - * - * ``` - // pray for a non-zero initial overlap - initRandomPureState(qureg); // works even for density matrices - - // minimize then renormalise - qreal tau = 10; // impatient infinity - int order = 4; - int reps = 100; - applyNonUnitaryTrotterizedPauliStrSumGadget(qureg, hamil, tau * 1i, order, reps); - setQuregToRenormalized(qureg); - - // ground-state (phi_0) - reportQureg(qureg); - - // lowest lying eigenvalue (lambda_0) - qreal expec = calcExpecPauliStrSum(qureg, hamil); - reportScalar("expec", expec); - * ``` - * - * Note degenerate eigenvalues will yield a pure superposition of the corresponding eigenstates, with - * coefficients informed by the initial, relative populations. + * - When @p angle is set to @f$ \theta = \iu \, \Delta \tau @f$ and @p sum = @f$ \hat{H} @f$ is Hermitian, + * this function (approximately) evolves @p qureg in imaginary-time for duration @f$ \Delta \tau @f$, + * effecting non-unitary propagator + @f[ + \exp(- \Delta \tau \hat{H}) + * @f] + * as utilised by applyTrotterizedImaginaryTimeEvolution(). * - * - When @p angle is real and @p sum is Hermitian (has approximately real coefficients), this - * function is equivalent to applyTrotterizedPauliStrSumGadget() + * - When @p angle is real and @p sum is Hermitian (i.e. has approximately real coefficients), the effected + * operation is unitary and this function becomes equivalent to applyTrotterizedPauliStrSumGadget(). * * @constraints + * * - This function only ever effects @f$ \exp \left(\iu \, \theta \, \hat{H} \right) @f$ exactly * when all PauliStr in @p sum = @f$ \hat{H} @f$ commute. * @@ -304,7 +265,6 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vec - /** * @defgroup trotter_timeevol Time evolution * @brief Functions for approximate dynamical simulation. @@ -317,15 +277,225 @@ extern "C" { #endif -/// @notyettested -/// @notyetvalidated -/// @notyetdoced +/** @notyettested + * + * Unitarily time evolves @p qureg for the duration @p time under the time-independent Hamiltonian @p hamil, + * as approximated by symmetrized Trotterisation of the specified @p order and number of cycles @p reps. + * See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * + * @formulae + * + * Let @f$ \hat{H} = @f$ @p hamil and @f$ t = @f$ @p time @f$ \in \mathbb{R} @f$. This function approximates + * the action of the unitary-time evolution operator/propagator + * @f[ + \hat{U}(t) = \exp \left(- \iu \, t \, \hat{H} \right), + * @f] + * as solves the time-independent Schrödinger equation. When @p qureg is a statevector @f$ \svpsi @f$, the + * resulting state approximates + * @f[ + \approx U(t) \svpsi + * @f] + * while when @p qureg is a density matrix @f$ \dmrho @f$, the result approximates + * @f[ + \approx U(t) \, \dmrho \, U(t)^\dagger. + * @f] + * + * @equivalences + * + * - This function merely wraps applyTrotterizedPauliStrSumGadget() which effects @f$ \exp(\iu \theta \hat{H}) @f$, + * passing @f$ \theta = - t @f$. + * + * @constraints + * + * - Unitarity requires that @p hamil is Hermitian and ergo contains only real coefficients. Validation will check that + * @p hamil is approximately Hermitian, permitting coefficients with imaginary components smaller (in magnitude) than + * epsilon. + * @f[ + \max\limits_{i} \Big|c_i| \le \valeps + * @f] + * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). The imaginary components + * of the Hamiltonian _are_ considered during simulation. + * + * - The @p time parameter is necessarily real to retain unitarity. It can be substituted for a strictly imaginary + * scalar to perform imaginary-time evolution (as per Wick rotation @f$ t \rightarrow - \iu \tau @f$) via + * applyTrotterizedImaginaryTimeEvolution(), or generalised to an arbitrary complex number through direct use of + * applyNonUnitaryTrotterizedPauliStrSumGadget(). + * + * - The simulated system is _closed_ with dynamics described fully by the Hamiltonian @p hamil. Open or otherwise noisy + * system dynamics can be simulated with applyTrotterizedPauliNoisyTimeEvolution(). + * + * - Simulation is exact such that the effected operation is precisely @f$ \exp(-\iu t \hat{H}) @f$ only when + * @p reps @f$ \rightarrow \infty @f$ or all terms in @p hamil commute with one another. + * + * @myexample + * + * ``` + Qureg qureg = createDensityQureg(10); + PauliStrSum hamil = createInlinePauliStrSum(R"( + 1 ZZI + 2 IZZ + 3 ZIZ + 1.5 XII + 2.5 IXI + 3.5 IIX + )"); + + qreal time = 0.8 * hbar; + int order = 4; + int reps = 20; + applyTrotterizedUnitaryTimeEvolution(qureg, hamil, time, order, reps); + * ``` + * + * @see + * - applyTrotterizedImaginaryTimeEvolution() + * - applyTrotterizedPauliNoisyTimeEvolution() + * - applyNonUnitaryTrotterizedPauliStrSumGadget() + * + * @param[in,out] qureg the state to modify. + * @param[in] hamil the Hamiltonian as a a weighted sum of Pauli strings. + * @param[in] time the duration over which to simulate evolution. + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) + * @param[in] reps the number of Trotter repetitions + * + * @throws @validationerror + * - if @p qureg or @p hamil are uninitialised. + * - if @p hamil contains non-identities on qubits beyond the size of @p qureg. + * - if @p hamil is not approximately Hermitian. + * - if @p order is not 1 nor a positive, @b even integer. + * - if @p reps is not a positive integer. + * + * @author Tyson Jones + */ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps); -/// @notyettested -/// @notyetvalidated -/// @notyetdoced +/** @notyetdoced + * + * Simulates imaginary-time evolution of @p qureg for the duration @p tau under the time-independent + * Hamiltonian @p hamil, as approximated by symmetrized Trotterisation of the specified @p order and + * number of cycles @p reps. See applyTrotterizedPauliStrSumGadget() for information about the Trotter + * method. + * + * > [!IMPORTANT] + * > This is a non-physical operation and breaks the normalisation of state which can be restored + * > via setQuregToRenormalized(). + * + * @formulae + * + * Let @f$ \hat{H} = @f$ @p hamil and @f$ \tau = @f$ @p tau @f$ \in \mathbb{R} @f$. This function + * approximates the action of the non-unitary imaginary-time propagator + * @f[ + \hat{V}(\tau) = \exp \left(- \tau \, \hat{H} \right), + * @f] + * as prescribed by Wick rotating (substituting time @f$ t @f$ for @f$ t \rightarrow -\iu \tau @f$) + * the time-independent Schrödinger equation. When @p qureg is a statevector @f$ \svpsi @f$, the + * resulting state approximates + * @f[ + \approx V(\tau) \svpsi + * @f] + * while when @p qureg is a density matrix @f$ \dmrho @f$, the result approximates + * @f[ + \approx V(\tau) \, \dmrho \, V(\tau)^\dagger. + * @f] + * + * @par Utility + * + * Imaginary-time evolution drives the system toward the (unnormalised) groundstate of the Hamiltonian. + * Let @f$ \{ \ket{\phi_i} \} @f$ and @f$ \{ \ket{\lambda_i} \} @f$ be the eigenstates and respective + * eigenvalues of @f$ \hat{H} @f$, which are real due to Hermiticity. + * @f[ + \hat{H} = \sum \limits_i \lambda_i \ket{\phi_i}\bra{\phi_i}, + \;\;\;\;\; \lambda_i \in \mathbb{R}. + * @f] + * + * - When @p qureg is a statevector @f$ \svpsi @f$ and can ergo be expressed in the basis of + * @f$ \{ \ket{\phi_i} \} @f$ as @f$ \svpsi = \sum_i \alpha_i \ket{\phi_i} @f$, + * this function approximates + * @f[ + \svpsi \, \rightarrow \, \hat{V}(\tau) \svpsi = + \sum\limits_i \alpha_i \exp(- \tau \, \lambda_i) \ket{\phi_i}. + * @f] + * - When @p qureg is a density matrix and is ergo expressible as + * @f$ \dmrho = \sum\limits_{ij} \alpha_{ij} \ket{\phi_i}\bra{\phi_j} @f$, this function effects + * @f[ + \dmrho \, \rightarrow \, \hat{V}(\tau) \dmrho \hat{V}(\tau)^\dagger = + \sum\limits_{ij} \alpha_{ij} \exp(-\tau (\lambda_i + \lambda_j)) \ket{\phi_i}\bra{\phi_j}. + * @f] + * + * As @f$ \tau \rightarrow \infty @f$, the resulting unnormalised state approaches statevector + * @f$ \svpsi \rightarrow \alpha_0 \exp(-\tau \lambda_0) \ket{\phi_0} @f$ or density matrix + * @f$ \dmrho \rightarrow \alpha_{0,0} \exp(-2 \tau \lambda_0) \ket{\phi_0}\bra{\phi_0} @f$, + * where @f$ \lambda_0 @f$ is the minimum eigenvalue and @f$ \ket{\phi_0} @f$ is the groundstate. + * Assuming the initial overlap @f$ \alpha_0 @f$ is not zero (or exponentially tiny), + * subsequent renormalisation via setQuregToRenormalized() produces the pure + * ground-state @f$ \ket{\phi_0} @f$ or @f$ \ket{\phi_0}\bra{\phi_0} @f$. + * + * Note degenerate minimum eigenvalues will yield a pure superposition of the corresponding + * eigenstates, with coefficients informed by the initial, relative populations. + * + * @equivalences + * + * - This function merely wraps applyNonUnitaryTrotterizedPauliStrSumGadget() which effects @f$ \exp(\iu \theta \hat{H}) @f$, + * passing @f$ \theta = \tau \iu @f$. + * + * @constraints + * + * - While the process of imaginary-time evolution is non-unitary (and non-physical), Hermiticity of @p hamil is still + * assumed, requiring it contains only real coefficients. Validation will check that @p hamil is _approximately_ Hermitian, + * permitting coefficients with imaginary components smaller (in magnitude) than epsilon. + * @f[ + \max\limits_{i} \Big|c_i| \le \valeps + * @f] + * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). Beware however that + * imaginary-time evolution under a non-Hermitian Hamiltonian will _not_ necessarily approach the lowest lying eigenstate + * (the eigenvalues may be non-real) so is likely of limited utility. + * + * - The @p tau parameter is necessarily real such that evolution approaches the groundstate (modulo renormalisation). + * It can generalised to an arbitrary complex number through direct use of applyNonUnitaryTrotterizedPauliStrSumGadget(). + * + * - Simulation is exact such that the effected operation is precisely @f$ \exp(-\tau \hat{H}) @f$ only when + * @p reps @f$ \rightarrow \infty @f$ or all terms in @p hamil commute with one another. + * + * @myexample + * + * ``` + // pray for a non-zero initial overlap + initRandomPureState(qureg); // works even for density matrices + + // minimize then renormalise + qreal tau = 10; // impatient infinity + int order = 4; + int reps = 100; + applyTrotterizedImaginaryTimeEvolution(qureg, hamil, tau, order, reps); + setQuregToRenormalized(qureg); + + // ground-state (phi_0) + reportQureg(qureg); + + // lowest lying eigenvalue (lambda_0) + qreal expec = calcExpecPauliStrSum(qureg, hamil); + reportScalar("expec", expec); + * ``` + * + * @see + * - applyTrotterizedUnitaryTimeEvolution() + * - applyNonUnitaryTrotterizedPauliStrSumGadget() + * + * @param[in,out] qureg the state to modify. + * @param[in] hamil the Hamiltonian as a a weighted sum of Pauli strings. + * @param[in] tau the duration over which to simulate imaginary-time evolution. + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) + * @param[in] reps the number of Trotter repetitions + * + * @throws @validationerror + * - if @p qureg or @p hamil are uninitialised. + * - if @p hamil contains non-identities on qubits beyond the size of @p qureg. + * - if @p hamil is not approximately Hermitian. + * - if @p order is not 1 nor a positive, @b even integer. + * - if @p reps is not a positive integer. + * + * @author Tyson Jones + */ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); From 19e18b2f9cae8d8f775746f038d693a22bc4b6d3 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Wed, 23 Jul 2025 23:40:39 +0200 Subject: [PATCH 05/19] updated dynamics examples to use time-evol in lieu of direct gadgets, for clarity --- examples/extended/dynamics.c | 7 +++---- examples/extended/dynamics.cpp | 7 +++---- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/examples/extended/dynamics.c b/examples/extended/dynamics.c index 8b94e0a7d..0d6763c1d 100644 --- a/examples/extended/dynamics.c +++ b/examples/extended/dynamics.c @@ -1,6 +1,5 @@ /** @file - * An example of using QuEST (primarily function - * applyTrotterizedPauliStrSumGadget()) to perform + * An example of using QuEST to perform closed * dynamical simulation via Trotterisation of the * unitary-time evolution operator. * @@ -158,7 +157,7 @@ int main() { for (int i=0; i qreal time = dt * (i+1); @@ -188,7 +187,7 @@ int main() { // verify results by uninterrupted higher-order simulation to target time initPlusState(qureg); - applyTrotterizedPauliStrSumGadget(qureg, hamil, -dt*steps, order+2, reps*steps); + applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps); reportScalar("final ", calcExpecPauliStrSum(qureg, observ)); // clean up diff --git a/examples/extended/dynamics.cpp b/examples/extended/dynamics.cpp index ae232e348..95245fb87 100644 --- a/examples/extended/dynamics.cpp +++ b/examples/extended/dynamics.cpp @@ -1,6 +1,5 @@ /** @file - * An example of using QuEST (primarily function - * applyTrotterizedPauliStrSumGadget()) to perform + * An example of using QuEST to perform closed * dynamical simulation via Trotterisation of the * unitary-time evolution operator. * @@ -155,7 +154,7 @@ int main() { for (int i=0; i qreal time = dt * (i+1); @@ -182,7 +181,7 @@ int main() { // verify results by uninterrupted higher-order simulation to target time initPlusState(qureg); - applyTrotterizedPauliStrSumGadget(qureg, hamil, -dt*steps, order+2, reps*steps); + applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt*steps, order+2, reps*steps); reportScalar("final ", calcExpecPauliStrSum(qureg, observ)); // clean up From 315ea414a54ab09ab572553cde1e33361caefc3b Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Wed, 23 Jul 2025 23:41:25 +0200 Subject: [PATCH 06/19] made imaginary-time evolution assert Hermiticity since it really only counts as imaginary-time evolution if the Hamiltonian is Hermitian - otherwise it's just sparkling non-unitary propagator Trotterisation --- quest/src/api/trotterisation.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 99ec32e8a..e90696728 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -207,6 +207,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea validate_quregFields(qureg, __func__); validate_pauliStrSumFields(hamil, __func__); validate_pauliStrSumTargets(hamil, qureg, __func__); + validate_pauliStrSumIsHermitian(hamil, __func__); validate_trotterParams(qureg, order, reps, __func__); // exp(-tau H) = exp(x i H) | x=tau*i From 01c51e137af0d79096545145eb55b0744b0bde9c Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Wed, 23 Jul 2025 23:57:35 +0200 Subject: [PATCH 07/19] patched non-unitary Trotter on density matrix Trotterisation of a PauliStrSum gadget sees a density matrix operated upon by a sequence of PauliStr gadgets, each like rho -> exp(i angle coeff term) rho exp(i angle coeff term)^dagger where angle is a user-parameter, term is a PauliStr, and coeff is its coefficient in the PauliStrSum. In the linearised setting, this becomes: ||rho> -> conj[exp(i angle coeff term)] (x) exp(i angle coeff term) ||rho> When the gadget is unitary (i.e. the PauliStrSum is Hermitian and the angle parameter is real), the dagger of the LHS operator above simplifies to exp(+- i angle coeff term) where the +- sign is admitted by conj(term), depending on the parity of the Y operators therein. This simplification was erroneously used in the general case, when the PauliStrSum is not necessarily Hermitian and angle can be complex. In that general case, we have instead that the LHS operator's conjugation simplifies to: exp(+- i conj(angle*coeff) term) This patch corrects the argument of the prescribed gadget. --- quest/src/api/trotterisation.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index e90696728..10a36ddac 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -37,7 +37,7 @@ void internal_applyFirstOrderTrotterRepetition( qcomp coeff = sum.coeffs[j]; PauliStr str = sum.strings[j]; - // effect |psi> -> exp(i angle * sum)|psi> + // effect |psi> -> exp(i angle * coeff * term)|psi> qcomp arg = angle * coeff; localiser_statevec_anyCtrlPauliGadget(qureg, ketCtrls, states, str, arg); @@ -49,8 +49,10 @@ void internal_applyFirstOrderTrotterRepetition( if (!postmultiply) continue; - // effect rho -> rho dagger(i angle * sum) - arg *= paulis_hasOddNumY(str) ? 1 : -1; + // effect rho -> rho exp(i angle * coeff * term)^dagger via linearised + // ||rho>> -> conj(exp(i angle * coeff * term)) (x) I ||rho>> + // = exp(+- i conj(angle) conj(coeff) term) (x) I ||rho>> + arg = std::conj(arg) * (paulis_hasOddNumY(str) ? 1 : -1); str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliGadget(qureg, braCtrls, states, str, arg); } From d661ce766225fb7b1d80b7533a9e23655e53cb8e Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Wed, 23 Jul 2025 23:59:58 +0200 Subject: [PATCH 08/19] added memo to generalise Linblad to PauliStrSum --- quest/src/api/trotterisation.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 10a36ddac..b340b711e 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -229,6 +229,11 @@ extern PauliStr paulis_getKetAndBraPauliStr(PauliStr str, Qureg qureg); extern "C" { +/// @todo +/// we will not expose the below function; instead, we will trivially generalise it +/// to where each jump operator is a PauliStrSum, hackily empowering (sub-optimal) +/// projectors including off-diagonal, which is much more useful + void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps) { validate_quregFields(qureg, __func__); validate_quregIsDensityMatrix(qureg, __func__); From 55222ad3b9c6513322d7bc69f25f1092b5b6f8ab Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sat, 16 Aug 2025 22:45:50 +0200 Subject: [PATCH 09/19] tidying use of paulis_hasOddNumY by replacing it with a function which directly answers the caller's question; what sign when multiplied upon the PauliStr produces its conjugate? --- quest/src/api/operations.cpp | 21 +++++++++++---------- quest/src/api/paulis.cpp | 7 +++++++ quest/src/api/trotterisation.cpp | 10 +++++----- 3 files changed, 23 insertions(+), 15 deletions(-) diff --git a/quest/src/api/operations.cpp b/quest/src/api/operations.cpp index 0ec30c72d..0a5b75235 100644 --- a/quest/src/api/operations.cpp +++ b/quest/src/api/operations.cpp @@ -29,8 +29,9 @@ using std::vector; * PRVIATE UTILITIES */ +extern int paulis_getSignOfPauliStrConj(PauliStr str); + extern bool paulis_isIdentity(PauliStr str); -extern bool paulis_hasOddNumY(PauliStr str); extern PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift); extern PauliStr paulis_getKetAndBraPauliStr(PauliStr str, Qureg qureg); @@ -1250,7 +1251,7 @@ void postMultiplyPauliStr(Qureg qureg, PauliStr str) { validate_quregIsDensityMatrix(qureg, __func__); validate_pauliStrTargets(qureg, str, __func__); - qcomp factor = paulis_hasOddNumY(str)? -1 : 1; // undo transpose + qcomp factor = paulis_getSignOfPauliStrConj(str); // undo transpose str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliTensor(qureg, {}, {}, str, factor); } @@ -1292,7 +1293,7 @@ void applyMultiStateControlledPauliStr(Qureg qureg, int* controls, int* states, // operation sinto a single tensor, i.e. +- (shift(str) (x) str), to // avoid superfluous re-enumeration of the state if (qureg.isDensityMatrix && numControls == 0) { - factor = paulis_hasOddNumY(str)? -1 : 1; + factor = paulis_getSignOfPauliStrConj(str); ctrlVec = util_getConcatenated(ctrlVec, util_getBraQubits(ctrlVec, qureg)); stateVec = util_getConcatenated(stateVec, stateVec); str = paulis_getKetAndBraPauliStr(str, qureg); @@ -1302,7 +1303,7 @@ void applyMultiStateControlledPauliStr(Qureg qureg, int* controls, int* states, // but density-matrix control qubits require two distinct operations if (qureg.isDensityMatrix && numControls > 0) { - factor = paulis_hasOddNumY(str)? -1 : 1; + factor = paulis_getSignOfPauliStrConj(str); ctrlVec = util_getBraQubits(ctrlVec, qureg); str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliTensor(qureg, ctrlVec, stateVec, str, factor); @@ -1366,7 +1367,7 @@ void postMultiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace) { // post-multiply each term in-turn, mixing into output qureg, then undo using idempotency for (qindex i=0; i prefixZ) { int sign = 1; diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index b340b711e..82b1e4d19 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -24,7 +24,7 @@ using std::vector; * INTERNAL UTILS */ -extern bool paulis_hasOddNumY(PauliStr str); +extern int paulis_getSignOfPauliStrConj(PauliStr str); extern PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift); void internal_applyFirstOrderTrotterRepetition( @@ -51,8 +51,8 @@ void internal_applyFirstOrderTrotterRepetition( // effect rho -> rho exp(i angle * coeff * term)^dagger via linearised // ||rho>> -> conj(exp(i angle * coeff * term)) (x) I ||rho>> - // = exp(+- i conj(angle) conj(coeff) term) (x) I ||rho>> - arg = std::conj(arg) * (paulis_hasOddNumY(str) ? 1 : -1); + // = exp(- i conj(angle coeff) sign term) (x) I ||rho>> + arg = - std::conj(arg) * paulis_getSignOfPauliStrConj(str); str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliGadget(qureg, braCtrls, states, str, arg); } @@ -279,7 +279,7 @@ void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qre // term of i conj(H) (x) I newStrings.push_back(paulis_getShiftedPauliStr(oldStr, qureg.numQubits)); - newCoeffs.push_back(1_i * (paulis_hasOddNumY(oldStr) ? -1 : 1) * std::conj(oldCoeff)); + newCoeffs.push_back(1_i * paulis_getSignOfPauliStrConj(oldStr) * std::conj(oldCoeff)); } // collect jump terms @@ -287,7 +287,7 @@ void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qre // gamma_k conj(L_k) (x) L_k newStrings.push_back(paulis_getKetAndBraPauliStr(jumps[n], qureg)); - newCoeffs.push_back(damps[n] * (paulis_hasOddNumY(jumps[n]) ? -1 : 1)); + newCoeffs.push_back(damps[n] * paulis_getSignOfPauliStrConj(jumps[n])); } // spoof a PauliStrSum to avoid superfluous alloc From b2240e17c58b795f9c69b5ecdb354ea43e424da9 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 18 Aug 2025 20:51:52 +0200 Subject: [PATCH 10/19] separated PauliStr and PauliStrSum sections --- quest/src/api/paulis.cpp | 59 ++++++++++++++++++++++-------------- quest/src/core/localiser.cpp | 10 +++--- 2 files changed, 41 insertions(+), 28 deletions(-) diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp index 397f63259..32b6bb408 100644 --- a/quest/src/api/paulis.cpp +++ b/quest/src/api/paulis.cpp @@ -87,7 +87,7 @@ void freeAllMemoryIfAnyAllocsFailed(PauliStrSum sum) { /* - * INTERNAL UTILITIES + * INTERNAL PauliStr UTILITIES * * callable by other internal files but which are not exposed in the header * because we do not wish to make them visible to users. Ergo other internal @@ -139,12 +139,6 @@ int paulis_getIndOfLefmostNonIdentityPauli(PauliStr* strings, qindex numStrings) } -int paulis_getIndOfLefmostNonIdentityPauli(PauliStrSum sum) { - - return paulis_getIndOfLefmostNonIdentityPauli(sum.strings, sum.numTerms); -} - - bool paulis_containsXOrY(PauliStr str) { int maxInd = paulis_getIndOfLefmostNonIdentityPauli(str); @@ -160,16 +154,6 @@ bool paulis_containsXOrY(PauliStr str) { } -bool paulis_containsXOrY(PauliStrSum sum) { - - for (qindex i=0; i,3> paulis_getSeparateInds(PauliStr str, Qureg qureg) { +array,3> paulis_getSeparateInds(PauliStr str) { vector iXYZ = paulis_getTargetInds(str); vector iX, iY, iZ; @@ -286,18 +270,25 @@ PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift) { } -PauliStr paulis_getKetAndBraPauliStr(PauliStr str, Qureg qureg) { +PauliStr paulis_getTensorProdOfPauliStr(PauliStr left, PauliStr right, int numQubits) { + + // computes left (tensor) right, assuming right is smaller than numQubits + PauliStr shifted = paulis_getShiftedPauliStr(left, numQubits); - PauliStr shifted = paulis_getShiftedPauliStr(str, qureg.numQubits); - // return a new stack PauliStr instance (avoiding C++20 initialiser) PauliStr out; - out.lowPaulis = str.lowPaulis | shifted.lowPaulis; - out.highPaulis = str.highPaulis | shifted.highPaulis; + out.lowPaulis = right.lowPaulis | shifted.lowPaulis; + out.highPaulis = right.highPaulis | shifted.highPaulis; return out; } +PauliStr paulis_getKetAndBraPauliStr(PauliStr str, Qureg qureg) { + + return paulis_getTensorProdOfPauliStr(str, str, qureg.numQubits); +} + + PAULI_MASK_TYPE paulis_getKeyOfSameMixedAmpsGroup(PauliStr str) { PAULI_MASK_TYPE key = 0; @@ -319,6 +310,28 @@ PAULI_MASK_TYPE paulis_getKeyOfSameMixedAmpsGroup(PauliStr str) { } + +/* + * INTERNAL PauliStrSum UTILITIES + */ + + +int paulis_getIndOfLefmostNonIdentityPauli(PauliStrSum sum) { + + return paulis_getIndOfLefmostNonIdentityPauli(sum.strings, sum.numTerms); +} + + +bool paulis_containsXOrY(PauliStrSum sum) { + + for (qindex i=0; i, vecto extern bool paulis_containsXOrY(PauliStr str); extern vector paulis_getTargetInds(PauliStr str); -extern std::array,3> paulis_getSeparateInds(PauliStr str, Qureg qureg); +extern std::array,3> paulis_getSeparateInds(PauliStr str); extern int paulis_getPrefixZSign(Qureg qureg, vector prefixZ) ; extern qcomp paulis_getPrefixPaulisElem(Qureg qureg, vector prefixY, vector prefixZ); @@ -1298,7 +1298,7 @@ void anyCtrlPauliTensorOrGadget(Qureg qureg, vector ctrls, vector ctrl // - prefix X,Y determine communication, because they apply bit-not to rank // - prefix Y,Z determine node-wide coefficient, because they contain rank-determined !=1 elements // - suffix X,Y,Z determine local amp coefficients - auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str, qureg); + auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str); auto [prefixX, suffixX] = util_getPrefixAndSuffixQubits(targsX, qureg); auto [prefixY, suffixY] = util_getPrefixAndSuffixQubits(targsY, qureg); auto [prefixZ, suffixZ] = util_getPrefixAndSuffixQubits(targsZ, qureg); @@ -1994,7 +1994,7 @@ qcomp getDensMatrExpecPauliStrTermOfOnlyThisNode(Qureg qureg, PauliStr str) { // caller must reduce the returned value between nodes if necessary // all ket-paulis are in the suffix state - auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str, qureg); + auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str); // optimised scenario when str = I if (targsX.empty() && targsY.empty() && targsZ.empty()) @@ -2019,7 +2019,7 @@ qcomp localiser_statevec_calcExpecPauliStr(Qureg qureg, PauliStr str) { // - prefix Y,Z determine node-wide coefficient, because they contain rank-determined !=1 elements // - suffix X,Y,Z determine local amp coefficients // noting that when !qureg.isDistributed, all paulis will be in suffix - auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str, qureg); + auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str); auto [prefixX, suffixX] = util_getPrefixAndSuffixQubits(targsX, qureg); auto [prefixY, suffixY] = util_getPrefixAndSuffixQubits(targsY, qureg); auto [prefixZ, suffixZ] = util_getPrefixAndSuffixQubits(targsZ, qureg); @@ -2107,7 +2107,7 @@ qcomp localiser_statevec_calcExpecPauliStrSum(Qureg qureg, PauliStrSum sum) { // for each term within the current group... for (auto& [str, coeff] : terms) { - auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str, qureg); + auto [targsX, targsY, targsZ] = paulis_getSeparateInds(str); auto [prefixX, suffixX] = util_getPrefixAndSuffixQubits(targsX, qureg); auto [prefixY, suffixY] = util_getPrefixAndSuffixQubits(targsY, qureg); auto [prefixZ, suffixZ] = util_getPrefixAndSuffixQubits(targsZ, qureg); From ca6e9f824ec8578f916e28ab9e4b5a16b76a98e9 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 18 Aug 2025 20:53:25 +0200 Subject: [PATCH 11/19] made createPauliStrSum validate enough RAM since the same facilities are anyway needed (with much better motivation) by Lindblad evolution --- quest/src/core/memory.cpp | 9 +++++++++ quest/src/core/memory.hpp | 2 ++ quest/src/core/validation.cpp | 28 +++++++++++++++++++++------- tests/unit/paulis.cpp | 6 ++++++ 4 files changed, 38 insertions(+), 7 deletions(-) diff --git a/quest/src/core/memory.cpp b/quest/src/core/memory.cpp index c8b81fc88..7f11494a1 100644 --- a/quest/src/core/memory.cpp +++ b/quest/src/core/memory.cpp @@ -388,6 +388,15 @@ bool mem_canSuperOpFitInMemory(int numQubits, qindex numBytesPerNode) { } +bool mem_canPauliStrSumFitInMemory(qindex numTerms, qindex numBytesPerNode) { + + // awkwardly arranged to avoid overflow when numTerms is too large + size_t numBytesPerTerm = sizeof(PauliStr) + sizeof(qcomp); + qindex maxNumTerms = numBytesPerNode / numBytesPerTerm; // floors + return numTerms <= maxNumTerms; +} + + /* * MEMORY ALLOCATION SUCCESS diff --git a/quest/src/core/memory.hpp b/quest/src/core/memory.hpp index 7b112b027..b624d6c03 100644 --- a/quest/src/core/memory.hpp +++ b/quest/src/core/memory.hpp @@ -101,6 +101,8 @@ bool mem_canMatrixFitInMemory(int numQubits, bool isDense, int numNodes, qindex bool mem_canSuperOpFitInMemory(int numQubits, qindex numBytesPerNode); +bool mem_canPauliStrSumFitInMemory(qindex numTerms, qindex numBytesPerNode); + /* diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index a8223cb3f..78d305bac 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -693,6 +693,9 @@ namespace report { string NEW_PAULI_STR_SUM_DIFFERENT_NUM_STRINGS_AND_COEFFS = "Given a different number of Pauli strings (${NUM_STRS}) and coefficients ${NUM_COEFFS}."; + string NEW_PAULI_STR_SUM_CANNOT_FIT_INTO_CPU_MEM = + "A PauliStrSum containing ${NUM_TERMS} terms cannot fit in the available RAM of ${NUM_BYTES} bytes."; + string NEW_PAULI_STR_SUM_STRINGS_ALLOC_FAILED = "Attempted allocation of the PauliStrSum's ${NUM_TERMS}-term array of Pauli strings (${NUM_BYTES} bytes total) unexpectedly failed."; @@ -3190,14 +3193,25 @@ void validate_controlAndPauliStrTargets(Qureg qureg, int ctrl, PauliStr str, con void validate_newPauliStrSumParams(qindex numTerms, const char* caller) { - // note we do not bother checking whether RAM has enough memory to contain - // the new Pauli sum, because the caller to this function has already - // been passed data of the same size (and it's unlikely the user is about - // to max RAM), and the memory requirements scale only linearly with the - // parameters (e.g. numTerms), unlike the exponential scaling of the memory - // of Qureg and CompMatr, for example - assertThat(numTerms > 0, report::NEW_PAULI_STR_SUM_NON_POSITIVE_NUM_STRINGS, {{"${NUM_TERMS}", numTerms}}, caller); + + // attempt to fetch RAM, and simply return if we fail; if we unknowingly + // didn't have enough RAM, then alloc validation will trigger later + size_t memPerNode = 0; + try { + memPerNode = mem_tryGetLocalRamCapacityInBytes(); + } catch(mem::COULD_NOT_QUERY_RAM e) { + return; + } + + // pedantically check whether the PauliStrSum fits in memory. This seems + // ridiculous/pointless because the user is expected to have already prepared + // data of an equivalent size (which is passed), but checking means we catch + // when the user has passed an erroneous 'numTerms' which is way too large, + // avoiding a seg fault + + bool fits = mem_canPauliStrSumFitInMemory(numTerms, memPerNode); + assertThat(fits, report::NEW_PAULI_STR_SUM_CANNOT_FIT_INTO_CPU_MEM, {{"${NUM_TERMS}", numTerms}, {"${NUM_BYTES}", memPerNode}}, caller); } void validate_newPauliStrSumMatchingListLens(qindex numStrs, qindex numCoeffs, const char* caller) { diff --git a/tests/unit/paulis.cpp b/tests/unit/paulis.cpp index 255d30c4b..8b6fc4f90 100644 --- a/tests/unit/paulis.cpp +++ b/tests/unit/paulis.cpp @@ -328,6 +328,12 @@ TEST_CASE( "createPauliStrSum", TEST_CATEGORY ) { REQUIRE_THROWS_WITH( createPauliStrSum(nullptr, nullptr, numTerms), ContainsSubstring("number of terms must be a positive integer") ); } + SECTION( "exceeds memory" ) { + + // can choose even a number of terms so large that its size (in bytes) overflows + REQUIRE_THROWS_WITH( createPauliStrSum(nullptr, nullptr, 1LL << 60), ContainsSubstring("cannot fit in the available RAM") ); + } + SECTION( "mismatching lengths" ) { // specific to the C++ interface From f0a31e0a38c59e67cb04e5a76f3548ff203a28a5 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 18 Aug 2025 21:03:40 +0200 Subject: [PATCH 12/19] added PauliStrSum jump operators allowing the (potentially sub-optimal) expression of arbitrary operators like projectors, etc --- quest/include/trotterisation.h | 2 +- quest/src/api/paulis.cpp | 108 +++++++++++++++++++++ quest/src/api/trotterisation.cpp | 156 ++++++++++++++++++++----------- quest/src/core/errors.cpp | 36 +++++++ quest/src/core/errors.hpp | 18 ++++ quest/src/core/utilities.cpp | 35 +++++++ quest/src/core/utilities.hpp | 5 + quest/src/core/validation.cpp | 82 ++++++++++++++++ quest/src/core/validation.hpp | 18 ++++ 9 files changed, 403 insertions(+), 57 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 74e420f2c..594db6da9 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -502,7 +502,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea /// @notyettested /// @notyetvalidated /// @notyetdoced -void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); +void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps); // end de-mangler diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp index 32b6bb408..b07c065fd 100644 --- a/quest/src/api/paulis.cpp +++ b/quest/src/api/paulis.cpp @@ -20,6 +20,7 @@ #include "quest/src/comm/comm_routines.hpp" #include +#include #include #include #include @@ -310,6 +311,32 @@ PAULI_MASK_TYPE paulis_getKeyOfSameMixedAmpsGroup(PauliStr str) { } +std::pair paulis_getPauliStrProd(PauliStr strA, PauliStr strB) { + + // a . b = coeff * (a ^ b) + PauliStr strOut; + strOut.lowPaulis = strA.lowPaulis ^ strB.lowPaulis; + strOut.highPaulis = strA.highPaulis ^ strB.highPaulis; + + // coeff = product of single-site product coeffs + qcomp coeff = 1; + for (int i=0; i= numQubits) + error_pauliStrSumHasMoreQubitsThanSpecifiedInTensorProd(); + if (out.numTerms != in.numTerms * in.numTerms) + error_pauliStrSumTensorProdHasIncorrectNumTerms(); + + // conj(in) (x) in = sum_jk conj(c_j) c_k conj(P_j) (x) P_k... + qindex i = 0; + for (qindex j=0; j= numQubits) + error_pauliStrSumHasMoreQubitsThanSpecifiedInConjShift(); + if (out.numTerms != in.numTerms) + error_pauliStrSumConjHasIncorrectNumTerms(); + + // where conj(c P) = conj(c) sign P + for (qindex i=0; i @@ -26,6 +27,10 @@ using std::vector; extern int paulis_getSignOfPauliStrConj(PauliStr str); extern PauliStr paulis_getShiftedPauliStr(PauliStr str, int pauliShift); +extern void paulis_setPauliStrSumToScaledTensorProdOfConjWithSelf(PauliStrSum out, qreal factor, PauliStrSum in, int numQubits); +extern void paulis_setPauliStrSumToScaledProdOfAdjointWithSelf(PauliStrSum out, qreal factor, PauliStrSum in); +extern void paulis_setPauliStrSumToShiftedConj(PauliStrSum out, PauliStrSum in, int numQubits); +extern qindex paulis_getNumTermsInPauliStrSumProdOfAdjointWithSelf(PauliStrSum in); void internal_applyFirstOrderTrotterRepetition( Qureg qureg, vector& ketCtrls, vector& braCtrls, @@ -109,6 +114,37 @@ void internal_applyAllTrotterRepetitions( /// implement these above or into another function? } +qindex internal_getNumTotalSuperPropagatorTerms(PauliStrSum hamil, PauliStrSum* jumps, int numJumps) { + + // this function returns 0 to indicate an overflow, which will never + // be confused for the correct non-overflowed output because hamil.numTerms>0 + int OVERFLOW_FLAG = 0; + + if (util_willProdOverflow({2,hamil.numTerms})) + return OVERFLOW_FLAG; + + // I (x) H + conj(H) (x) I + qindex numTerms = 2 * hamil.numTerms; + + for (int i=0; i newStrings; - vector newCoeffs; + // validate memory allocations for all super-propagator terms + vector superStrings; + vector superCoeffs; + auto callbackString = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(PauliStr), __func__); }; + auto callbackCoeff = [&]() { validate_tempAllocSucceeded(false, numSuperTerms, sizeof(qcomp), __func__); }; + util_tryAllocVector(superStrings, numSuperTerms, callbackString); + util_tryAllocVector(superCoeffs, numSuperTerms, callbackCoeff); - // premature optimisation - qindex numNewTerms = 2 * hamil.numTerms + numJumps; - newStrings.reserve(numNewTerms); - newCoeffs.reserve(numNewTerms); + qindex superTermInd = 0; // collect -i[H,rho] terms for (qindex n=0; n #include #include +#include #include #include #include @@ -339,6 +340,39 @@ qcomp util_getPowerOfI(size_t exponent) { return values[exponent % 4]; } +bool util_willProdOverflow(vector terms) { + + qindex max = std::numeric_limits::max(); + qindex prod = 1; + + for (auto x : terms) { + + // division floors, so prod strictly exceed + if (prod > max / x) + return true; + + prod *= x; + } + + return false; +} + +bool util_willSumOverflow(vector terms) { + + qindex max = std::numeric_limits::max(); + qindex sum = 0; + + for (auto x : terms) { + + if (sum >= max - x) + return true; + + sum += x; + } + + return false; +} + /* @@ -1107,6 +1141,7 @@ void util_tryAllocVector(vector &vec, qindex size, std::function &vec, qindex size, std::function errFunc) { tryAllocVector(vec, size, errFunc); } void util_tryAllocVector(vector &vec, qindex size, std::function errFunc) { tryAllocVector(vec, size, errFunc); } void util_tryAllocVector(vector &vec, qindex size, std::function errFunc) { tryAllocVector(vec, size, errFunc); } +void util_tryAllocVector(vector &vec, qindex size, std::function errFunc) { tryAllocVector(vec, size, errFunc); } // cuQuantum needs a vector overload, which we additionally define when qreal!=double. Gross! #if FLOAT_PRECISION != 2 diff --git a/quest/src/core/utilities.hpp b/quest/src/core/utilities.hpp index c514821bf..a55e94ab5 100644 --- a/quest/src/core/utilities.hpp +++ b/quest/src/core/utilities.hpp @@ -105,6 +105,10 @@ bool util_isApproxReal(qcomp num, qreal epsilon); qcomp util_getPowerOfI(size_t exponent); +bool util_willProdOverflow(vector terms); + +bool util_willSumOverflow(vector terms); + /* @@ -401,6 +405,7 @@ void util_tryAllocVector(vector &vec, qindex size, std::function &vec, qindex size, std::function errFunc); void util_tryAllocVector(vector &vec, qindex size, std::function errFunc); void util_tryAllocVector(vector &vec, qindex size, std::function errFunc); +void util_tryAllocVector(vector &vec, qindex size, std::function errFunc); // cuQuantum needs a vector overload, which we additionally define when qreal!=double. Gross! #if FLOAT_PRECISION != 2 diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index 78d305bac..51b61b7fc 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -877,6 +877,10 @@ namespace report { "Cannot perform this ${NUM_TARGS}-target operation upon a ${NUM_QUREG_QUBITS}-qubit density-matrix distributed between ${NUM_NODES} nodes, since each node's communication buffer (with capacity for ${NUM_QUREG_AMPS_PER_NODE} amps) cannot simultaneously store the ${NUM_TARG_AMPS} mixed remote amplitudes."; + /* + * TROTTERISATION PARAMETERS + */ + string INVALID_TROTTER_ORDER = "Invalid Trotter order (${ORDER}). The order parameter must be positive and even, or unity."; @@ -884,6 +888,23 @@ namespace report { "Invalid number of Trotter repetitions (${REPS}). The number of repetitions must be positive."; + /* + * TIME EVOLUTION PARAMETERS + */ + + string NEGATIVE_NUM_LINDBLAD_JUMP_OPS = + "The number of jump operators must be zero or positive."; + + string NEGATIVE_LINDBLAD_DAMPING_RATE = + "One or more damping rates were negative, beyond the tolerance set by the validation epsilon."; + + string NUM_LINDBLAD_SUPER_PROPAGATOR_TERMS_OVERFLOWED = + "The given Hamiltonian and jump operators suggest a Lindbladian superpropagator with more weighted Pauli strings than can be stored in a qindex type."; + + string NEW_LINDBLAD_SUPER_PROPAGATOR_CANNOT_FIT_INTO_CPU_MEM = + "The Lindbladian superpropagator would contain ${NUM_TERMS} weighted Pauli strings and exceed the available RAM of ${NUM_BYTES} bytes."; + + /* * CHANNEL PARAMETERS */ @@ -3765,6 +3786,12 @@ void validate_mixedAmpsFitInNode(Qureg qureg, int numTargets, const char* caller assertThat(qureg.numAmpsPerNode >= numTargAmps, msg, vars, caller); } + + +/* + * TROTTERISATION PARAMETERS + */ + void validate_trotterParams(Qureg qureg, int order, int reps, const char* caller) { bool isEven = (order % 2) == 0; @@ -3774,6 +3801,61 @@ void validate_trotterParams(Qureg qureg, int order, int reps, const char* caller +/* + * TIME EVOLUTION PARAMETERS + */ + +void validate_lindbladJumpOps(PauliStrSum* jumps, int numJumps, Qureg qureg, const char* caller) { + + assertThat(numJumps >= 0, report::NEGATIVE_NUM_LINDBLAD_JUMP_OPS, caller); + + // @todo + // these error messages report as if each jump operator was "the" PauliStrSum + // to a function expecting one, and should be tailored to them being "a" jump op + for (int n=0; n= 0, report::NEGATIVE_NUM_LINDBLAD_JUMP_OPS, caller); + + if (isNumericalValidationDisabled()) + return; + + // in lieu of asserting positivity, we somewhat unusually permit small negative + // damping rates just for consistency with other numerical validation tolerances + for (int n=0; n= - global_validationEpsilon, report::NEGATIVE_LINDBLAD_DAMPING_RATE, caller); +} + +void validate_numLindbladSuperPropagatorTerms(qindex numSuperTerms, const char* caller) { + + assertThat(numSuperTerms != 0, report::NUM_LINDBLAD_SUPER_PROPAGATOR_TERMS_OVERFLOWED, caller); + + // attempt to fetch RAM, and simply return if we fail; if we unknowingly + // didn't have enough RAM, then alloc validation will trigger later + size_t memPerNode = 0; + try { + memPerNode = mem_tryGetLocalRamCapacityInBytes(); + } catch(mem::COULD_NOT_QUERY_RAM e) { + return; + } + + // check whether the superpropagator fits in memory + bool fits = mem_canPauliStrSumFitInMemory(numSuperTerms, memPerNode); + assertThat(fits, report::NEW_LINDBLAD_SUPER_PROPAGATOR_CANNOT_FIT_INTO_CPU_MEM, {{"${NUM_TERMS}", numSuperTerms}, {"${NUM_BYTES}", memPerNode}}, caller); + +} + + + /* * CHANNEL PARAMETERS */ diff --git a/quest/src/core/validation.hpp b/quest/src/core/validation.hpp index 92baac843..e15bca4e9 100644 --- a/quest/src/core/validation.hpp +++ b/quest/src/core/validation.hpp @@ -414,10 +414,28 @@ void validate_rotationAxisNotZeroVector(qreal x, qreal y, qreal z, const char* c void validate_mixedAmpsFitInNode(Qureg qureg, int numTargets, const char* caller); + + +/* + * TROTTERISATION PARAMETERS + */ + void validate_trotterParams(Qureg qureg, int order, int reps, const char* caller); +/* + * TIME EVOLUTION PARAMETERS + */ + +void validate_lindbladJumpOps(PauliStrSum* jumps, int numJumps, Qureg qureg, const char* caller); + +void validate_lindbladDampingRates(qreal* damps, int numJumps, const char* caller); + +void validate_numLindbladSuperPropagatorTerms(qindex numSuperTerms, const char* caller); + + + /* * DECOHERENCE */ From f4edce350e2ebb95b97576e08a42777b1b70878b Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 18 Aug 2025 22:01:28 +0200 Subject: [PATCH 13/19] eep --- quest/src/api/trotterisation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 4e4cb3ff2..f9adb6633 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -142,7 +142,7 @@ qindex internal_getNumTotalSuperPropagatorTerms(PauliStrSum hamil, PauliStrSum* } // indicate no overflow - return OVERFLOW_FLAG; + return numTerms; } From af839067ec2eadf2172c5f51e7c2ace9479af5bc Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 18 Aug 2025 23:36:40 +0200 Subject: [PATCH 14/19] clarified 'postmultiply' to inverse 'onlyApplyLeft' --- quest/src/api/trotterisation.cpp | 59 ++++++++++++++++---------------- tests/unit/trotterisation.cpp | 2 -- 2 files changed, 30 insertions(+), 31 deletions(-) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index f9adb6633..64c0dc19a 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -34,7 +34,7 @@ extern qindex paulis_getNumTermsInPauliStrSumProdOfAdjointWithSelf(PauliStrSum i void internal_applyFirstOrderTrotterRepetition( Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, bool postmultiply, bool reverse + vector& states, PauliStrSum sum, qcomp angle, bool onlyLeftApply, bool reverse ) { // apply each sum term as a gadget, in forward or reverse order for (qindex i=0; i rho exp(i angle * coeff * term)^dagger via linearised @@ -65,14 +65,14 @@ void internal_applyFirstOrderTrotterRepetition( void internal_applyHigherOrderTrotterRepetition( Qureg qureg, vector& ketCtrls, vector& braCtrls, - vector& states, PauliStrSum sum, qcomp angle, int order, bool postmultiply + vector& states, PauliStrSum sum, qcomp angle, int order, bool onlyLeftApply ) { if (order == 1) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, postmultiply, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, onlyLeftApply, false); } else if (order == 2) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, postmultiply, false); - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, postmultiply, true); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, onlyLeftApply, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, onlyLeftApply, true); } else { qreal p = 1. / (4 - std::pow(4, 1./(order-1))); @@ -80,17 +80,17 @@ void internal_applyHigherOrderTrotterRepetition( qcomp b = (1-4*p) * angle; int lower = order - 2; - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); // angle -> a - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower, postmultiply); // angle -> b - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, postmultiply); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, onlyLeftApply); // angle -> a + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, onlyLeftApply); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower, onlyLeftApply); // angle -> b + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, onlyLeftApply); + internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower, onlyLeftApply); } } void internal_applyAllTrotterRepetitions( Qureg qureg, int* controls, int* states, int numControls, - PauliStrSum sum, qcomp angle, int order, int reps, bool postmultiply + PauliStrSum sum, qcomp angle, int order, int reps, bool onlyLeftApply ) { // exp(i angle sum) = identity when angle=0 if (angle == qcomp(0,0)) @@ -106,7 +106,7 @@ void internal_applyAllTrotterRepetitions( // perform carefully-ordered sequence of gadgets for (int r=0; r -> U |psi>, rho -> U rho U^dagger + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply); } void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) { @@ -171,8 +172,8 @@ void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle validate_pauliStrSumIsHermitian(sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps, onlyLeftApply); } void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, PauliStrSum sum, qreal angle, int order, int reps) { @@ -182,8 +183,8 @@ void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, Pauli validate_controlAndPauliStrSumTargets(qureg, control, sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps, onlyLeftApply); } void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int numControls, PauliStrSum sum, qreal angle, int order, int reps) { @@ -193,8 +194,8 @@ void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls validate_controlsAndPauliStrSumTargets(qureg, controls, numControls, sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps, onlyLeftApply); } void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls, int* states, int numControls, PauliStrSum sum, qreal angle, int order, int reps) { @@ -205,8 +206,8 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con validate_controlStates(states, numControls, __func__); // permits states==nullptr validate_trotterParams(qureg, order, reps, __func__); - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps, onlyLeftApply); } } // end de-mangler @@ -239,8 +240,8 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal // exp(-i t H) = exp(x i H) | x=-t qcomp angle = - time; - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply); } void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps) { @@ -252,8 +253,8 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea // exp(-tau H) = exp(x i H) | x=tau*i qcomp angle = qcomp(0, tau); - bool postmultiply = true; - internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, postmultiply); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, hamil, angle, order, reps, onlyLeftApply); } } // end de-mangler @@ -350,8 +351,8 @@ void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qre // effect exp(t S) = exp(x i S) | x=-i*time, left-multiplying only qcomp angle = qcomp(0, -time); - bool postmultiply = false; - internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, superSum, angle, order, reps, postmultiply); + bool onlyLeftApply = true; + internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, superSum, angle, order, reps, onlyLeftApply); } } // end de-mangler diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp index 85e9b47f0..3ac362fc3 100644 --- a/tests/unit/trotterisation.cpp +++ b/tests/unit/trotterisation.cpp @@ -39,6 +39,4 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); - -// TODO: just a draft function void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); From eb6cf1213bfe3e188cf179af1e0596f648970ba3 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Tue, 19 Aug 2025 17:03:13 +0200 Subject: [PATCH 15/19] removed Pauli from applyTrotterizedPauliNoisyTimeEvolution --- quest/include/trotterisation.h | 2 +- quest/src/api/trotterisation.cpp | 2 +- tests/unit/trotterisation.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 594db6da9..9a6a0bd4c 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -502,7 +502,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea /// @notyettested /// @notyetvalidated /// @notyetdoced -void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps); +void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps); // end de-mangler diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 375794c59..760432ee7 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -267,7 +267,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea extern "C" { -void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps) { +void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps) { validate_quregFields(qureg, __func__); validate_quregIsDensityMatrix(qureg, __func__); validate_pauliStrSumFields(hamil, __func__); diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp index 3ac362fc3..bc2c7aab6 100644 --- a/tests/unit/trotterisation.cpp +++ b/tests/unit/trotterisation.cpp @@ -39,4 +39,4 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); -void applyTrotterizedPauliNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); +void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); From 3551128d6334248ac667726d6181626a51f8d327 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Tue, 19 Aug 2025 17:03:37 +0200 Subject: [PATCH 16/19] added doc for applyTrotterizedNoisyTimeEvolution --- quest/include/trotterisation.h | 145 ++++++++++++++++++++++++++++++++- 1 file changed, 142 insertions(+), 3 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 9a6a0bd4c..7c20c24b6 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -499,9 +499,148 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); -/// @notyettested -/// @notyetvalidated -/// @notyetdoced +/** @notyettested + * + * Simulates open dynamics of @p qureg as per the Lindblad master equation, under the time-independent + * Hamiltonian @p hamil and jump operators @p jumps with corresponding damping rates @p damps, with + * evolution approximated by symmetrized Trotterisation of the specified @p order and number of cycles + * @p reps. See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * + * @formulae + * + * Let @f$ \rho = @f$ @p qureg, @f$ \hat{H} = @f$ @p hamil, @f$ t = @f$ @p time, and denote the @f$ i @f$-th + * element of @p damps and @p jumps as @f$ \gamma_i @f$ and @f$ \hat{J}_i @f$ respectively. The Lindblad + * master equation prescribes that @f$ \rho @f$ time-evolves according to + * @f[ + \frac{\mathrm{d}}{\mathrm{d}t} \rho = -\iu [\hat{H}, \rho] + \sum\limits_i \gamma_i \left( + \hat{J}_i \rho \hat{J}_i^\dagger - \frac{1}{2} \left\{ \hat{J}_i^\dagger \hat{J}_i, \rho \right\} + \right). + * @f] + * This function works by building a superoperator of the right-hand-side which acts upon the space of + * linearised @f$\rho@f$, + * @f[ + \boldsymbol{L} = -\iu \left( \hat{\id} \otimes \hat{H} - \hat{H}^* \otimes \hat{\id} \right) + + \sum\limits_i \gamma_i \left( + \hat{J}_i^* \otimes \hat{J}_i - \frac{1}{2} \hat{\id} \otimes (\hat{J}^\dagger J_i) + - \frac{1}{2} (\hat{J}^\dagger J_i)^* \otimes \hat{\id} + \right), + * @f] + * as a non-Hermitian weighted sum of Pauli strings (a PauliStrSum). The superoperator @f$ \boldsymbol{L} @f$ + * informs a superpropagator which exactly solves evolution as: + * @f[ + \ket{\rho(t)} = \exp\left( t \boldsymbol{L} \right) \ket{\rho(0)}. + * @f] + * This function approximates the superpropagator @f$ \exp\left( t \boldsymbol{L} \right) @f$ using a higher-order + * symmetrized Suzuki-Trotter decomposition, as informed by parameters @p order and @p reps. + * + * @par Utility + * + * This function simulates time evolution of an open system, where the jump operators model interactions with + * the environment. This can capture sophisticated decoherence processes of the quantum state which are untenable + * to model as discrete operations with functions like mixKrausMap(). This function also proves useful for + * preparing realistic, physical input states to quantum metrological circuits, or the general high-performance + * simulation of digital time evolution of condensed matter systems. + * + * @equivalences + * + * - When `numJumps = 0`, evolution is unitary and the Lindblad master equation simplifes to the Liouville–von Neumann + * equation, which is equivalently (and more efficiently) simulated via applyTrotterizedUnitaryTimeEvolution(). + * + * @constraints + * + * - Each damping rate in @p damps is expected to be a zero or positive number, in order for evolution to be trace + * preserving. Validation will assert that each damping rate @f$ \gamma_i @f$ satisfies + * @f[ + \min\limits_{i} \gamma_i \ge - \valeps + * @f] + * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). Non-trace-preserving, + * negative damping rates can be simulated by disabling numerical validation via `setValidationEpsilon(0)`. + * + * - The @p time parameter is necessarily real, and cannot be generalised to imaginary or complex like in other + * functions. Generalisation is trivially numerically possible, but has no established physical meaning and so + * is not exposed in the API. Please open an issue on Github for advice on complex-time simulation. + * + * - Simulation is exact only when @p reps @f$ \rightarrow \infty @f$ or all terms in the superoperator + * @f$ \boldsymbol{L} @f$ incidentally commute with one another, and otherwise incorporates Trotter error. + * Unlike for unitary evolution, Trotter error _does_ break normalisation of the state and so this function + * is generally non-trace-preserving. In theory, normalisation can be restored with setQuregToRenormalized() + * though noticable norm-breaking indicates evolution was inaccurate, and should instead be repeated with + * increased @p order or @p reps parameters. + * + * - The function instantiates superoperator @f$ \boldsymbol{L} @f$ above as a temporary PauliStrSum, incurring a + * memory and time overhead which grows quadratically with the number of terms in @p hamil, plus quadratically + * with the number in each jump operator. These overheads may prove prohibitively costly for PauliStrSum + * containing very many terms. + * + * @myexample + * + * ``` + // |+><+| + Qureg qureg = createDensityQureg(3); + initPlusState(qureg); + + PauliStrSum hamil = createInlinePauliStrSum(R"( + 1 IIX + 2 IYI + 3 ZZZ + )"); + + // |0><0| + PauliStrSum jump1 = createInlinePauliStrSum(R"( + 0.5 I + 0.5 Z + )"); + + // |1><0| + PauliStrSum jump2 = createInlinePauliStrSum(R"( + 0.5 X + -0.5i Y + )"); + + // "noisiness" + qreal damps[] = {.3, .4}; + PauliStrSum jumps[] = {jump1, jump2}; + int numJumps = 2; + + reportScalar("initial energy", calcExpecPauliStrSum(qureg, hamil)); + + // time and accuracy + qreal time = 0.5; + int order = 4; + int reps = 100; + applyTrotterizedNoisyTimeEvolution(qureg, hamil, damps, jumps, numJumps, time, order, reps); + + reportScalar("final energy", calcExpecPauliStrSum(qureg, hamil)); + * ``` + * + * @see + * - applyTrotterizedUnitaryTimeEvolution() + * - applyTrotterizedImaginaryTimeEvolution() + * + * @param[in,out] qureg the density-matrix state to evolve and modify. + * @param[in] hamil the Hamiltonian of the qubit system (excludes any environment). + * @param[in] damps the damping rates of each jump operator in @p jumps. + * @param[in] jumps the jump operators specified as PauliStrSum. + * @param[in] numJumps the length of list @p jumps (and @p damps). + * @param[in] time the duration through which to evolve the state. + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...). + * @param[in] reps the number of Trotter repetitions. + * + * @throws @validationerror + * - if @p qureg, @p hamil or any element of @p jumps are uninitialised. + * - if @p qureg is not a density matrix. + * - if @p hamil or any element of @p jumps contains non-identities on qubits beyond the size of @p qureg. + * - if @p hamil is not approximately Hermitian. + * - if @p numJumps is negative. + * - if any element of @p damps is not approximately positive. + * - if the total number of Lindbladian superoperator terms overflows the `qindex` type. + * - if all Lindbladian superoperator terms cannot simultaneously fit into CPU memory. + * - if memory allocation of the Lindbladian superoperator terms unexpectedly fails. + * - if @p order is not 1 nor a positive, @b even integer. + * - if @p reps is not a positive integer. + * + * @author Tyson Jones + */ void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStrSum* jumps, int numJumps, qreal time, int order, int reps); From 0dc6454445d9f4edc8f9729b32517cc27a6068fa Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Tue, 19 Aug 2025 17:03:49 +0200 Subject: [PATCH 17/19] patched time-evol doc --- quest/include/trotterisation.h | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 7c20c24b6..84b34b57c 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -123,7 +123,7 @@ extern "C" { * only real coefficients. Validation will check that @p sum is approximately Hermitian, permitting * coefficients with imaginary components smaller (in magnitude) than epsilon. * @f[ - \max\limits_{i} \Big|c_i| \le \valeps + \max\limits_{i} |c_i| \le \valeps * @f] * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). * Otherwise, use applyNonUnitaryTrotterizedPauliStrSumGadget() to permit non-Hermitian @p sum @@ -140,8 +140,8 @@ extern "C" { * @param[in,out] qureg the state to modify. * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. * @param[in] angle the prefactor of @p sum times @f$ i @f$ in the exponent. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...). + * @param[in] reps the number of Trotter repetitions. * * @throws @validationerror * - if @p qureg or @p sum are uninitialised. @@ -221,8 +221,8 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con * @param[in,out] qureg the state to modify. * @param[in] sum a weighted sum of Pauli strings to approximately exponentiate. * @param[in] angle an effective prefactor of @p sum in the exponent. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...). + * @param[in] reps the number of Trotter repetitions. * * @throws @validationerror * - if @p qureg or @p sum are uninitialised. @@ -311,7 +311,7 @@ extern "C" { * @p hamil is approximately Hermitian, permitting coefficients with imaginary components smaller (in magnitude) than * epsilon. * @f[ - \max\limits_{i} \Big|c_i| \le \valeps + \max\limits_{i} |c_i| \le \valeps * @f] * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). The imaginary components * of the Hamiltonian _are_ considered during simulation. @@ -322,10 +322,11 @@ extern "C" { * applyNonUnitaryTrotterizedPauliStrSumGadget(). * * - The simulated system is _closed_ with dynamics described fully by the Hamiltonian @p hamil. Open or otherwise noisy - * system dynamics can be simulated with applyTrotterizedPauliNoisyTimeEvolution(). + * system dynamics can be simulated with applyTrotterizedNoisyTimeEvolution(). * * - Simulation is exact such that the effected operation is precisely @f$ \exp(-\iu t \hat{H}) @f$ only when - * @p reps @f$ \rightarrow \infty @f$ or all terms in @p hamil commute with one another. + * @p reps @f$ \rightarrow \infty @f$ or all terms in @p hamil commute with one another. Conveniently, Trotter error + * does _not_ break normalisation of the state since the approximating circuit remains unitary. * * @myexample * @@ -348,14 +349,14 @@ extern "C" { * * @see * - applyTrotterizedImaginaryTimeEvolution() - * - applyTrotterizedPauliNoisyTimeEvolution() + * - applyTrotterizedNoisyTimeEvolution() * - applyNonUnitaryTrotterizedPauliStrSumGadget() * * @param[in,out] qureg the state to modify. * @param[in] hamil the Hamiltonian as a a weighted sum of Pauli strings. * @param[in] time the duration over which to simulate evolution. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...). + * @param[in] reps the number of Trotter repetitions. * * @throws @validationerror * - if @p qureg or @p hamil are uninitialised. @@ -369,7 +370,7 @@ extern "C" { void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps); -/** @notyetdoced +/** @notyettested * * Simulates imaginary-time evolution of @p qureg for the duration @p tau under the time-independent * Hamiltonian @p hamil, as approximated by symmetrized Trotterisation of the specified @p order and @@ -444,7 +445,7 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal * assumed, requiring it contains only real coefficients. Validation will check that @p hamil is _approximately_ Hermitian, * permitting coefficients with imaginary components smaller (in magnitude) than epsilon. * @f[ - \max\limits_{i} \Big|c_i| \le \valeps + \max\limits_{i} |c_i| \le \valeps * @f] * where the validation epsilon @f$ \valeps @f$ can be adjusted with setValidationEpsilon(). Beware however that * imaginary-time evolution under a non-Hermitian Hamiltonian will _not_ necessarily approach the lowest lying eigenstate @@ -484,8 +485,8 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal * @param[in,out] qureg the state to modify. * @param[in] hamil the Hamiltonian as a a weighted sum of Pauli strings. * @param[in] tau the duration over which to simulate imaginary-time evolution. - * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...) - * @param[in] reps the number of Trotter repetitions + * @param[in] order the order of the Trotter-Suzuki decomposition (e.g. @p 1, @p 2, @p 4, ...). + * @param[in] reps the number of Trotter repetitions. * * @throws @validationerror * - if @p qureg or @p hamil are uninitialised. From dc874ddf564e979e869bb562916d29f3a8a5f118 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Tue, 19 Aug 2025 17:04:12 +0200 Subject: [PATCH 18/19] documented failed mathcal-rendering attempts for posterity! --- utils/docs/Doxyfile | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/utils/docs/Doxyfile b/utils/docs/Doxyfile index daa1f0143..1d1086fe3 100644 --- a/utils/docs/Doxyfile +++ b/utils/docs/Doxyfile @@ -1927,7 +1927,10 @@ MATHJAX_RELPATH = # MATHJAX_EXTENSIONS = ams # This tag requires that the tag USE_MATHJAX is set to YES. -MATHJAX_EXTENSIONS = +### Tyson: +### these were added to attemptedly fix \mathcal rendering +### but they did nothing, alas! I retain for posterity +MATHJAX_EXTENSIONS = TeX/AMSmath TeX/AMSsymbols # The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces # of code that will be used on startup of the MathJax code. See the MathJax site @@ -2101,7 +2104,10 @@ PAPER_TYPE = a4 # If left blank no extra packages will be included. # This tag requires that the tag GENERATE_LATEX is set to YES. -EXTRA_PACKAGES = times, amsmath, xr, amsfonts, tikz +### Tyson: +### \mathcal{L} does not render despite amssymb appearing below +### (and MathJax given the AMSsymbols extension elsewhere), grr! +EXTRA_PACKAGES = xr times amsmath amssymb amsfonts tikz # The LATEX_HEADER tag can be used to specify a user-defined LaTeX header for # the generated LaTeX document. The header should contain everything until the From ad21e011d3aeeb224e779242b752305734cb06a2 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Tue, 19 Aug 2025 17:44:52 +0200 Subject: [PATCH 19/19] reformatted Trotter ref in doc --- quest/include/trotterisation.h | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 84b34b57c..ca234c5d5 100644 --- a/quest/include/trotterisation.h +++ b/quest/include/trotterisation.h @@ -198,7 +198,7 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con * @f] * via a Trotter-Suzuki decomposition of the specified @p order and number of repetitions (@p reps). * - * See applyTrotterizedPauliStrSumGadget() for more information about the decomposition. + * > See applyTrotterizedPauliStrSumGadget() for more information about the decomposition. * * @equivalences * @@ -281,7 +281,6 @@ extern "C" { * * Unitarily time evolves @p qureg for the duration @p time under the time-independent Hamiltonian @p hamil, * as approximated by symmetrized Trotterisation of the specified @p order and number of cycles @p reps. - * See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. * * @formulae * @@ -299,6 +298,8 @@ extern "C" { * @f[ \approx U(t) \, \dmrho \, U(t)^\dagger. * @f] + * + * > See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. * * @equivalences * @@ -374,8 +375,7 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal * * Simulates imaginary-time evolution of @p qureg for the duration @p tau under the time-independent * Hamiltonian @p hamil, as approximated by symmetrized Trotterisation of the specified @p order and - * number of cycles @p reps. See applyTrotterizedPauliStrSumGadget() for information about the Trotter - * method. + * number of cycles @p reps. * * > [!IMPORTANT] * > This is a non-physical operation and breaks the normalisation of state which can be restored @@ -398,6 +398,8 @@ void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal * @f[ \approx V(\tau) \, \dmrho \, V(\tau)^\dagger. * @f] + * + * > See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. * * @par Utility * @@ -505,7 +507,7 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea * Simulates open dynamics of @p qureg as per the Lindblad master equation, under the time-independent * Hamiltonian @p hamil and jump operators @p jumps with corresponding damping rates @p damps, with * evolution approximated by symmetrized Trotterisation of the specified @p order and number of cycles - * @p reps. See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * @p reps. * * @formulae * @@ -534,6 +536,8 @@ void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qrea * This function approximates the superpropagator @f$ \exp\left( t \boldsymbol{L} \right) @f$ using a higher-order * symmetrized Suzuki-Trotter decomposition, as informed by parameters @p order and @p reps. * + * > See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * * @par Utility * * This function simulates time evolution of an open system, where the jump operators model interactions with