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 diff --git a/quest/include/trotterisation.h b/quest/include/trotterisation.h index 6c16e3d9a..ca234c5d5 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,41 +105,43 @@ 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. * @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 * 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] 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] 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. * * @throws @validationerror * - if @p qureg or @p sum are uninitialised. @@ -144,6 +153,7 @@ extern "C" { * @see * - applyPauliGadget() * - applyNonUnitaryTrotterizedPauliStrSumGadget() + * - applyTrotterizedUnitaryTimeEvolution() * * @author Tyson Jones */ @@ -176,92 +186,43 @@ 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] * 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 * - * - 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. * * @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. @@ -304,6 +265,399 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, std::vec +/** + * @defgroup trotter_timeevol Time evolution + * @brief Functions for approximate dynamical simulation. + * @{ + */ + + +#ifdef __cplusplus +extern "C" { +#endif + + +/** @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. + * + * @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] + * + * > See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * + * @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} |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 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. Conveniently, Trotter error + * does _not_ break normalisation of the state since the approximating circuit remains unitary. + * + * @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() + * - 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. + * + * @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 + * + * 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. + * + * > [!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] + * + * > See applyTrotterizedPauliStrSumGadget() for information about the Trotter method. + * + * @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} |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); + + +/** @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. + * + * @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. + * + * > 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 + * 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); + + +// end de-mangler +#ifdef __cplusplus +} +#endif + +/** @} */ + + + #endif // TROTTERISATION_H /** @} */ // (end file-wide doxygen defgroup) diff --git a/quest/src/api/operations.cpp b/quest/src/api/operations.cpp index ce5edb579..5c415fbc5 100644 --- a/quest/src/api/operations.cpp +++ b/quest/src/api/operations.cpp @@ -30,8 +30,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); @@ -966,7 +967,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); @@ -976,7 +977,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); @@ -1230,8 +1231,8 @@ void applyNonUnitaryPauliGadget(Qureg qureg, PauliStr str, qcomp angle) { if (!qureg.isDensityMatrix) return; - // conj(e^i(a)XZ) = e^(-i conj(a)XZ) but conj(Y)=-Y, so odd-Y undoes phase negation - phase = std::conj(phase) * (paulis_hasOddNumY(str) ? 1 : -1); + // conj(e^i(a)P) = e^(-i s conj(a) P) + phase = - std::conj(phase) * paulis_getSignOfPauliStrConj(str); str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliGadget(qureg, {}, {}, str, phase); } @@ -1273,8 +1274,8 @@ void applyMultiStateControlledPauliGadget(Qureg qureg, int* controls, int* state if (!qureg.isDensityMatrix) return; - // conj(e^iXZ) = e^(-iXZ), but conj(Y)=-Y, so odd-Y undoes phase negation - phase *= paulis_hasOddNumY(str) ? 1 : -1; + // conj(e^(i a P)) = e^(-i s a P) + phase *= - paulis_getSignOfPauliStrConj(str); ctrlVec = util_getBraQubits(ctrlVec, qureg); str = paulis_getShiftedPauliStr(str, qureg.numQubits); localiser_statevec_anyCtrlPauliGadget(qureg, ctrlVec, stateVec, str, phase); diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp index c770e5fb9..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 @@ -87,7 +88,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 +140,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 +155,6 @@ bool paulis_containsXOrY(PauliStr str) { } -bool paulis_containsXOrY(PauliStrSum sum) { - - for (qindex i=0; i prefixZ) { int sign = 1; @@ -239,7 +231,7 @@ qindex paulis_getTargetBitMask(PauliStr str) { } -array,3> paulis_getSeparateInds(PauliStr str, Qureg qureg) { +array,3> paulis_getSeparateInds(PauliStr str) { vector iXYZ = paulis_getTargetInds(str); vector iX, iY, iZ; @@ -279,18 +271,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; @@ -312,6 +311,54 @@ 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 @@ -24,12 +25,16 @@ 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); +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, - vector& states, PauliStrSum sum, qcomp angle, 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 -> 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); + // term finished upon statevector if (!qureg.isDensityMatrix) continue; - // effect rho -> rho dagger(i angle * sum) - arg *= paulis_hasOddNumY(str) ? 1 : -1; + // Linbladian propagator is only ever pre-multiplied + if (onlyLeftApply) + continue; + + // 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 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); } @@ -53,14 +65,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 onlyLeftApply ) { if (order == 1) { - internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle, onlyLeftApply, 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, onlyLeftApply, false); + internal_applyFirstOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, angle/2, onlyLeftApply, true); } else { qreal p = 1. / (4 - std::pow(4, 1./(order-1))); @@ -68,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); // angle -> a - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, b, lower); // angle -> b - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); - internal_applyHigherOrderTrotterRepetition(qureg, ketCtrls, braCtrls, states, sum, a, lower); + 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 + PauliStrSum sum, qcomp angle, int order, int reps, bool onlyLeftApply ) { // exp(i angle sum) = identity when angle=0 if (angle == qcomp(0,0)) @@ -94,7 +106,7 @@ void internal_applyAllTrotterRepetitions( // perform carefully-ordered sequence of gadgets for (int r=0; r0 + 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 -> 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) { @@ -127,7 +172,8 @@ void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle validate_pauliStrSumIsHermitian(sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - internal_applyAllTrotterRepetitions(qureg, nullptr, nullptr, 0, sum, angle, order, reps); + 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) { @@ -137,7 +183,8 @@ void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, Pauli validate_controlAndPauliStrSumTargets(qureg, control, sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - internal_applyAllTrotterRepetitions(qureg, &control, nullptr, 1, sum, angle, order, reps); + 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) { @@ -147,7 +194,8 @@ void applyMultiControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* controls validate_controlsAndPauliStrSumTargets(qureg, controls, numControls, sum, __func__); validate_trotterParams(qureg, order, reps, __func__); - internal_applyAllTrotterRepetitions(qureg, controls, nullptr, numControls, sum, angle, order, reps); + 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) { @@ -158,7 +206,8 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, int* con validate_controlStates(states, numControls, __func__); // permits states==nullptr validate_trotterParams(qureg, order, reps, __func__); - internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps); + bool onlyLeftApply = false; + internal_applyAllTrotterRepetitions(qureg, controls, states, numControls, sum, angle, order, reps, onlyLeftApply); } } // end de-mangler @@ -173,3 +222,138 @@ void applyMultiStateControlledTrotterizedPauliStrSumGadget(Qureg qureg, 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); + + qindex superTermInd = 0; + + // collect -i[H,rho] terms + for (qindex n=0; n, 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); 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/utilities.cpp b/quest/src/core/utilities.cpp index e30aa86d6..207b76d7f 100644 --- a/quest/src/core/utilities.cpp +++ b/quest/src/core/utilities.cpp @@ -28,6 +28,7 @@ #include #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 a8223cb3f..51b61b7fc 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."; @@ -874,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."; @@ -881,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 */ @@ -3190,14 +3214,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) { @@ -3751,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; @@ -3760,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 */ 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 diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp index 174b7c66e..bc2c7aab6 100644 --- a/tests/unit/trotterisation.cpp +++ b/tests/unit/trotterisation.cpp @@ -34,3 +34,9 @@ void applyControlledTrotterizedPauliStrSumGadget(Qureg qureg, int control, Pauli 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); + +void applyTrotterizedUnitaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal time, int order, int reps); + +void applyTrotterizedImaginaryTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal tau, int order, int reps); + +void applyTrotterizedNoisyTimeEvolution(Qureg qureg, PauliStrSum hamil, qreal* damps, PauliStr* jumps, int numJumps, qreal time, int order, int reps); 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