patched Trotter sign and added non-unitary Trotter (#647)

- applyTrotterizedPauliStrSumGadget() was documented to effect exp(i t H) but actually effected exp(-i t H), eep! Thankfully the doc warned the function was untested. It now correctly effects exp(i t H).
- defensively removed hardcoding of the scalar responsible for the above bug, which undoes the coefficient convention of the applyPauliGadget() functions
- added applyNonUnitaryTrotterizedPauliStrSumGadget() which permits a non-Hermitian Hamiltonian (i.e. with non-negligible imaginary components of the coefficients) and a complex angle parameter. Among other things, this permits simple imaginary-time evolution
- added the full doc for both functions
- fixed defunct doxygen command (@cppoverload)

Note that both these Trotter functions remain without unit tests since impending new functions are anticipated which will generalise the tests.

Author: TysonRayJones

quest/include/calculations.h

@@ -336,7 +336,7 @@ qcomp calcExpecNonHermitianFullStateDiagMatrPower(Qureg qureg, FullStateDiagMatr
 /// @notyettested
 /// @notyetdoced
 /// @notyetvalidated
-/// @cppoverload
+/// @cppvectoroverload
 /// @see calcProbOfMultiQubitOutcome()
 qreal calcProbOfMultiQubitOutcome(Qureg qureg, std::vector<int> qubits, std::vector<int> outcomes);
 
@@ -354,7 +354,7 @@ std::vector<qreal> calcProbsOfAllMultiQubitOutcomes(Qureg qureg, std::vector<int
 /// @notyettested
 /// @notyetdoced
 /// @notyetvalidated
-/// @cppoverload
+/// @cppvectoroverload
 /// @see calcPartialTrace()
 Qureg calcPartialTrace(Qureg qureg, std::vector<int> traceOutQubits);
 
@@ -363,7 +363,7 @@ Qureg calcPartialTrace(Qureg qureg, std::vector<int> traceOutQubits);
 /// @notyettested
 /// @notyetdoced
 /// @notyetvalidated
-/// @cppoverload
+/// @cppvectoroverload
 /// @see calcReducedDensityMatrix()
 Qureg calcReducedDensityMatrix(Qureg qureg, std::vector<int> retainQubits);
 

quest/include/operations.h

@@ -1834,6 +1834,7 @@ void multiplyPauliGadget(Qureg qureg, PauliStr str, qreal angle);
      qcomp factor = cexp(- theta / 2 * 1.i);
      setQuregToSuperposition(factor, qureg, 0,qureg,0,qureg);
  *   ```
+ * - Passing @p angle=0 is equivalent to effecting the identity, leaving the state unchanged.
  *
  * @myexample
  * ```
@@ -1850,13 +1851,20 @@ void multiplyPauliGadget(Qureg qureg, PauliStr str, qreal angle);
     // concisely
     applyPauliGadget(qureg, getInlinePauliStr("XYZ",{0,1,7}), theta);
  * ```
- * - Passing @p angle=0 is equivalent to effecting the identity, leaving the state unchanged.
+ *
+ * @see
+ *  - applyNonUnitaryPauliGadget()
  */
 void applyPauliGadget(Qureg qureg, PauliStr str, qreal angle);
 
-/// @notyetdoced
+
+/** @notyetdoced
+ * 
+ * This function generalises applyPauliGadget() to accept a complex angle.
+ */
 void applyNonUnitaryPauliGadget(Qureg qureg, PauliStr str, qcomp angle);
 
+
 /// @notyetdoced
 void applyControlledPauliGadget(Qureg qureg, int control, PauliStr str, qreal angle);
 
@@ -2265,8 +2273,13 @@ extern "C" {
 void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace);
 
 
-/** @notyetdoced
- * @notyettested
+/** @notyettested
+ * 
+ * Effects (an approximation to) the exponential of @p sum, weighted by @p angle, upon @p qureg,
+ * via the symmetrized Trotter-Suzuki decomposition (<a href="https://arxiv.org/abs/math-ph/0506007">arXiv</a>).
+ * 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 
  * 
@@ -2275,14 +2288,18 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace);
       \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 (necessarily real) coefficient of the @f$ j @f$-th PauliStr @f$ \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[
@@ -2315,10 +2332,161 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace);
  * 
  * > These formulations are taken from 'Finding Exponential Product Formulas
  * > of Higher Orders', Naomichi Hatano and Masuo Suzuki (2005) (<a href="https://arxiv.org/abs/math-ph/0506007">arXiv</a>).
+ * 
+ * @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);
 
 
+/** @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
 }

quest/src/api/operations.cpp

@@ -1130,18 +1130,20 @@ void multiplyPauliStrSum(Qureg qureg, PauliStrSum sum, Qureg workspace) {
     // workspace -> qureg, and qureg -> sum * qureg
 }
 
-void applyFirstOrderTrotter(Qureg qureg, PauliStrSum sum, qreal angle, bool reverse) {
+void applyFirstOrderTrotter(Qureg qureg, PauliStrSum sum, qcomp angle, bool reverse) {
 
     // (internal, invoked by applyTrotterizedPauliStrSumGadget)
 
     for (qindex i=0; i<sum.numTerms; i++) {
         int j = reverse? sum.numTerms - i - 1 : i;
-        qreal arg = 2 * angle * std::real(sum.coeffs[j]); // 2 undoes Gadget convention
-        applyPauliGadget(qureg, sum.strings[j], arg); // caller disabled valiation therein
+
+        // effect exp(i angle * sum) by undoing gadget pre-factor
+        qcomp arg = angle * sum.coeffs[j] / util_getPhaseFromGateAngle(1);
+        applyNonUnitaryPauliGadget(qureg, sum.strings[j], arg); // caller disabled valiation therein
     }
 }
 
-void applyHigherOrderTrotter(Qureg qureg, PauliStrSum sum, qreal angle, int order) {
+void applyHigherOrderTrotter(Qureg qureg, PauliStrSum sum, qcomp angle, int order) {
 
     // (internal, invoked by applyTrotterizedPauliStrSumGadget)
 
@@ -1154,8 +1156,8 @@ void applyHigherOrderTrotter(Qureg qureg, PauliStrSum sum, qreal angle, int orde
 
     } else {
         qreal p = 1. / (4 - std::pow(4, 1./(order-1)));
-        qreal a = p * angle;
-        qreal b = (1-4*p) * angle;
+        qcomp a = p * angle;
+        qcomp b = (1-4*p) * angle;
 
         int lower = order - 2;
         applyFirstOrderTrotter(qureg, sum, a, lower);
@@ -1166,15 +1168,15 @@ void applyHigherOrderTrotter(Qureg qureg, PauliStrSum sum, qreal angle, int orde
     }
 }
 
-void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) {
+void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps) {
     validate_quregFields(qureg, __func__);
     validate_pauliStrSumFields(sum, __func__);
     validate_pauliStrSumTargets(sum, qureg, __func__);
-    validate_pauliStrSumIsHermitian(sum, __func__);
     validate_trotterParams(qureg, order, reps, __func__);
+    // sum is permitted to be non-Hermitian
 
     // exp(i angle sum) = identity when angle=0
-    if (angle == 0)
+    if (angle == qcomp(0,0))
         return;
 
     // record validation state then disable to avoid repeated
@@ -1196,6 +1198,20 @@ void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle
     /// implement these above or into another function?
 }
 
+void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps) {
+
+    // validate inputs here despite re-validation below so that func name is correct in error message
+    validate_quregFields(qureg, __func__);
+    validate_pauliStrSumFields(sum, __func__);
+    validate_pauliStrSumTargets(sum, qureg, __func__);
+    validate_trotterParams(qureg, order, reps, __func__);
+
+    // crucially, require that sum coefficients are real
+    validate_pauliStrSumIsHermitian(sum, __func__);
+
+    applyNonUnitaryTrotterizedPauliStrSumGadget(qureg, sum, angle, order, reps);
+}
+
 } // end de-mangler
 
 

quest/src/core/utilities.cpp

@@ -912,12 +912,13 @@ qreal util_getPhaseFromGateAngle(qreal angle) {
 
     return - angle / 2;
 }
-
 qcomp util_getPhaseFromGateAngle(qcomp angle) {
 
-    return -angle / qcomp(2.0, 0.0);
+    return - angle / 2;
 }
 
+
+
 /*
  * DECOHERENCE FACTORS
  */

quest/src/core/utilities.hpp

@@ -349,7 +349,6 @@ util_VectorIndexRange util_getLocalIndRangeOfVectorElemsWithinNode(int rank, qin
  */
 
 qreal util_getPhaseFromGateAngle(qreal angle);
-
 qcomp util_getPhaseFromGateAngle(qcomp angle);
 
 

tests/unit/operations.cpp

@@ -200,7 +200,7 @@ void TEST_ON_CACHED_QUREG_AND_MATRIX(quregCache quregs, matrixCache matrices, au
  * - applyDiagMatrPower:    diagpower
  * - applyCompMatr:         compmatr
  * - applyPauliStr:         paulistr
- * - applyPauliGadgt:       pauligad
+ * - applyPauliGadget:      pauligad
 */
 
 enum ArgsFlag { none, scalar, axisrots, diagmatr, diagpower, compmatr, paulistr, pauligad };
@@ -1959,3 +1959,5 @@ TEST_CASE( "applyNonUnitaryPauliGadget", TEST_CATEGORY ) {
  */
 
 void applyTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qreal angle, int order, int reps);
+
+void applyNonUnitaryTrotterizedPauliStrSumGadget(Qureg qureg, PauliStrSum sum, qcomp angle, int order, int reps);