Fixed QITE algorithm impl

This mostly completed the impl.
Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

Also, added python examples

python/examples/qite_deuteron_h2_example.py

+# Simple 2-qubit demonstration of the Quatum Imaginary Time Evolution algorithm
+# for Deuteron H2 Hamiltonian
+# Reference: https://arxiv.org/pdf/1912.06226.pdf
+# Expected minimum value: -1.74886
+import xacc,sys, numpy as np
+import matplotlib.pyplot as plt
+
+# Get access to the desired QPU and
+# allocate some qubits to run on
+qpu = xacc.getAccelerator('qpp')
+
+# Construct the H2 Hamiltonian
+ham = xacc.getObservable('pauli', '5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1')
+
+xacc.qasm('''.compiler xasm
+.circuit prep
+.qbit q
+X(q[0]);
+''')
+prep_circuit = xacc.getCompiled('prep')
+
+# We need 2 qubits for this case (H2)
+buffer = xacc.qalloc(2)
+
+# Horizontal axis: 0 -> 0.3 (step = 0.05)
+# The number of Trotter steps 
+nbSteps = 6
+
+# The Trotter step size
+stepSize = 0.05
+
+# Create the QITE algorithm
+qite = xacc.getAlgorithm('qite', {
+                        'accelerator': qpu,
+                        'observable': ham,
+                        'step-size': stepSize,
+                        'steps': nbSteps,
+                        'ansatz': prep_circuit
+                        })
+
+result = qite.execute(buffer)
+
+# Expected result: ~ -1.74886
+print('Min energy value = ', buffer.getInformation('opt-val'))
+E = buffer.getInformation('exp-vals')
+# Reproduce Fig. 3(a) of https://arxiv.org/pdf/1912.06226.pdf
+plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro-', label = 'XACC QITE')
+plt.ylim((-2.0, -0.25))
+plt.yticks(np.arange(-2.0, -0.25, step=0.25)) 
+plt.grid()
+plt.show()
\ No newline at end of file

python/examples/qite_deuteron_h3_example.py

+# Simple 3-qubit demonstration of the Quatum Imaginary Time Evolution algorithm
+# for Deuteron H3 Hamiltonian
+# Reference: https://arxiv.org/pdf/1912.06226.pdf
+# Expected minimum value: -2.04482
+import xacc,sys, numpy as np
+import matplotlib.pyplot as plt
+
+# Get access to the desired QPU and
+# allocate some qubits to run on
+qpu = xacc.getAccelerator('qpp')
+
+# Construct the H3 Hamiltonian
+ham = xacc.getObservable('pauli', '5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1 + 9.625 - 9.625 Z2 - 3.91 X1 X2 - 3.91 Y1 Y2')
+
+xacc.qasm('''.compiler xasm
+.circuit prep
+.qbit q
+X(q[0]);
+''')
+prep_circuit = xacc.getCompiled('prep')
+
+# We need 3 qubits for this case (H3)
+buffer = xacc.qalloc(3)
+
+# Horizontal axis: 0 -> 0.3 (step = 0.05)
+# The number of Trotter steps 
+nbSteps = 6
+
+# The Trotter step size
+stepSize = 0.05
+
+# Create the QITE algorithm
+qite = xacc.getAlgorithm('qite', {
+                        'accelerator': qpu,
+                        'observable': ham,
+                        'step-size': stepSize,
+                        'steps': nbSteps,
+                        'ansatz': prep_circuit
+                        })
+
+result = qite.execute(buffer)
+
+# Expected result: ~ -2.04482
+print('Min energy value = ', buffer.getInformation('opt-val'))
+E = buffer.getInformation('exp-vals')
+# Reproduce Fig. 3(b) of https://arxiv.org/pdf/1912.06226.pdf
+plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro-', label = 'XACC QITE')
+plt.ylim((-2.5, -0.0))
+plt.yticks(np.arange(-2.5, 0.0, step=0.5)) 
+plt.grid()
+plt.show()
\ No newline at end of file

quantum/plugins/algorithms/qite/qite.cpp

@@ -120,7 +120,6 @@ arma::cx_mat createSMatrix(const std::vector<std::string>& in_pauliOps, const st
     xacc::quantum::PauliOperator left(leftOp);
     xacc::quantum::PauliOperator right(rightOp);
     auto product = left * right;
-    // std::cout << left.toString() << " * " << right.toString() << " = " << product.toString() << "\n";
     const auto index = findMatchingPauliIndex(in_pauliOps, product.toString());
     return in_tomoExp[index]*product.coefficient();
   };
@@ -175,10 +174,27 @@ bool QITE::initialize(const HeterogeneousMap &parameters)
     m_observable = xacc::as_shared_ptr(parameters.getPointerLike<Observable>("observable"));
   }
 
-  m_analytical = true;
+  m_analytical = false;
   if (parameters.keyExists<bool>("analytical")) 
   {
     m_analytical = parameters.get<bool>("analytical");
+    if (m_analytical)
+    {
+      // Default initial state is 0
+      m_initialState = 0;
+      if (parameters.keyExists<int>("initial-state"))
+      {
+        m_initialState = parameters.get<int>("initial-state");
+      } 
+    }
+  }
+
+  m_ansatz = nullptr;
+  // Ansatz here is just a state preparation circuit:
+  // e.g. if we want to start in state |01>, not |00>
+  if (parameters.pointerLikeExists<CompositeInstruction>("ansatz")) 
+  {
+    m_ansatz = parameters.getPointerLike<CompositeInstruction>("ansatz");
   }
 
   m_approxOps.clear();
@@ -294,8 +310,6 @@ double QITE::calcCurrentEnergy(int in_nbQubits) const
 std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuffer>& in_buffer, std::shared_ptr<CompositeInstruction> in_kernel, std::shared_ptr<Observable> in_hmTerm) const
 {
   const auto pauliOps = generatePauliPermutation(in_buffer->size());
-  assert(in_hmTerm->getSubTerms().size() == 1);
-  assert(in_hmTerm->getNonIdentitySubTerms().size() == 1);
 
   // Step 1: Observe the kernels using the various Pauli
   // operators to calculate S and b.
@@ -340,18 +354,13 @@ std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuff
   {
     for (int j = 0; j < sMatDim; ++j)
     {
-      S_Mat(i, j) = calcSmatEntry(sigmaExpectation, i, j);
+      const auto entry = calcSmatEntry(sigmaExpectation, i, j);
+      S_Mat(i, j) = std::abs(entry) > 1e-12 ? entry : 0.0;
     }
   }
 
   // b vector:
   const auto obsProjCoeffs = observableToVec(in_hmTerm, pauliOps);
-  // std::cout << "Observable Pauli Vec: [";
-  // for (const auto& elem: obsProjCoeffs)
-  // {
-  //   std::cout << elem << ", ";
-  // }
-  // std::cout << "]\n";
 
   // Calculate c: Eq. 3 in https://arxiv.org/pdf/1901.07653.pdf
   double c = 1.0;
@@ -360,7 +369,6 @@ std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuff
     c -= 2.0 * m_dBeta * obsProjCoeffs[i] * sigmaExpectation[i];
   }
   
-  // std::cout << "c = " << c << "\n";
   for (int i = 0; i < sMatDim; ++i)
   {
     std::complex<double> b = (sigmaExpectation[i]/ std::sqrt(c) - sigmaExpectation[i])/ m_dBeta;
@@ -373,22 +381,12 @@ std::shared_ptr<Observable> QITE::calcAOps(const std::shared_ptr<AcceleratorBuff
     }
     b = I*b - I*std::conj(b);
     // Set b_Vec
-    b_Vec(i) = b;
+    b_Vec(i) = std::abs(b) > 1e-12 ? b : 0.0;
   }
-
-  // std::cout << "S Matrix: \n" << S_Mat << "\n"; 
-  // std::cout << "B Vector: \n" << b_Vec << "\n"; 
-  // Add regularizer
-  arma::cx_mat dalpha(sMatDim, sMatDim, arma::fill::eye); 
-  dalpha = 0.1 * dalpha;
-
-  auto lhs = S_Mat + S_Mat.st() + dalpha;
+  
+  auto lhs = S_Mat + S_Mat.st();
   auto rhs = -b_Vec;
-
-  // std::cout << "LHS Matrix: \n" << lhs << "\n"; 
-  // std::cout << "RHS Vector: \n" << rhs << "\n"; 
   arma::cx_vec a_Vec = arma::solve(lhs, rhs);
-  // std::cout << "Result A Vector: \n" << a_Vec << "\n"; 
 
   // Now, we have the decomposition of A observable in the basis of
   // all possible Pauli combinations.
@@ -421,19 +419,20 @@ void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
     // Run on hardware/simulator using quantum gates/measure
     // Initial energy
     m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
-    
+    auto hamOp = std::make_shared<xacc::quantum::PauliOperator>();
+    for (const auto& hamTerm : m_observable->getNonIdentitySubTerms())
+    {
+      *hamOp = *hamOp + *(std::dynamic_pointer_cast<xacc::quantum::PauliOperator>(hamTerm));
+    }
+
     // Time stepping:
     for (int i = 0; i < m_nbSteps; ++i)
     {
-      for (const auto& hamTerm : m_observable->getNonIdentitySubTerms())
-      {
-        // Propagates the state via Trotter steps:
-        auto kernel = constructPropagateCircuit();
-        // Optimizes/calculates next A ops
-        auto nextAOps = calcAOps(buffer, kernel, hamTerm);
-        m_approxOps.emplace_back(nextAOps);
-      }
-
+      // Propagates the state via Trotter steps:
+      auto kernel = constructPropagateCircuit();
+      // Optimizes/calculates next A ops
+      auto nextAOps = calcAOps(buffer, kernel, hamOp);
+      m_approxOps.emplace_back(nextAOps); 
       m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
     }
     assert(m_energyAtStep.size() == m_nbSteps + 1);
@@ -496,38 +495,36 @@ void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
 
     // Initial state
     arma::cx_vec psiVec(1 << buffer->size(), arma::fill::zeros);
-    psiVec(1) = 1.0;
-    // Time stepping:
-    for (int i = 0; i < m_nbSteps; ++i)
+    psiVec(m_initialState) = 1.0;
+    arma::cx_mat hMat(1 << buffer->size(), 1 << buffer->size(), arma::fill::zeros);
+    auto identityTerm = m_observable->getIdentitySubTerm();
+    if (identityTerm)
     {
-      arma::cx_mat hMat(1 << buffer->size(), 1 << buffer->size(), arma::fill::zeros);
-      
-      auto identityTerm = m_observable->getIdentitySubTerm();
-      if (identityTerm)
-      {
-        arma::cx_mat idTerm(1 << buffer->size(), 1 << buffer->size(), arma::fill::eye);
-        hMat += identityTerm->coefficient() * idTerm;
-      }
+      arma::cx_mat idTerm(1 << buffer->size(), 1 << buffer->size(), arma::fill::eye);
+      hMat += identityTerm->coefficient() * idTerm;
+    }
 
-      double stateVecNorm = 1.0;
-      for (const auto& hamTerm : m_observable->getNonIdentitySubTerms())
+    for (const auto& hamTerm : m_observable->getNonIdentitySubTerms())
+    {
+      auto pauliCast = std::dynamic_pointer_cast<xacc::quantum::PauliOperator>(hamTerm);
+      const auto hamMat = pauliCast->toDenseMatrix(buffer->size());
+    
+      for (int i = 0; i < hMat.n_rows; ++i)
       {
-        auto pauliCast = std::dynamic_pointer_cast<xacc::quantum::PauliOperator>(hamTerm);
-        const auto hamMat = pauliCast->toDenseMatrix(buffer->size());
-      
-        for (int i = 0; i < hMat.n_rows; ++i)
+        for (int j = 0; j < hMat.n_cols; ++j)
         {
-          for (int j = 0; j < hMat.n_cols; ++j)
-          {
-            const int index = i*hMat.n_rows + j;
-            hMat(i, j) += hamMat[index];
-          }
+          const int index = i*hMat.n_rows + j;
+          hMat(i, j) += hamMat[index];
         }
       }
+    }    
+
+    // Time stepping:
+    for (int i = 0; i < m_nbSteps; ++i)
+    {
       double normAfter = 0.0;
       arma::cx_vec dPsiVec(1 << buffer->size(), arma::fill::zeros);
       std::tie(dPsiVec, normAfter) = expMinusHamTerm(hMat, psiVec, m_dBeta);
-      stateVecNorm *= normAfter;
       // Eq. 8, SI of https://arxiv.org/pdf/1901.07653.pdf
       dPsiVec = dPsiVec - psiVec;
       std::vector<std::complex<double>> pauliExp;

quantum/plugins/algorithms/qite/qite.hpp

@@ -61,6 +61,10 @@ private:
   mutable std::vector<double> m_energyAtStep;
   // If a pure analytical run is requested.
   bool m_analytical;
+  // For analytical solver only: the initial state
+  // For accelerator-based simulation, the Ansatz is used to
+  // prepare the initial state.
+  int m_initialState;
 };
 } // namespace algorithm
 } // namespace xacc

quantum/plugins/algorithms/qite/tests/QITETester.cpp

@@ -54,14 +54,18 @@ TEST(QITETester, checkDeuteuronH2)
   std::shared_ptr<xacc::Observable> observable = std::make_shared<xacc::quantum::PauliOperator>();
   observable->fromString("5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1");
   auto acc = xacc::getAccelerator("qpp");
-  const int nbSteps = 10;
+  const int nbSteps = 5;
   const double stepSize = 0.1;
+  auto compiler = xacc::getCompiler("xasm");
+  auto ir = compiler->compile(R"(__qpu__ void f(qbit q) { X(q[0]); })", nullptr);
+  auto x = ir->getComposite("f");
 
   const bool initOk =  qite->initialize({
     std::make_pair("accelerator", acc),
     std::make_pair("steps", nbSteps),
     std::make_pair("observable", observable),
-    std::make_pair("step-size", stepSize)
+    std::make_pair("step-size", stepSize),
+    std::make_pair("ansatz", x),
   });
 
   EXPECT_TRUE(initOk);
@@ -71,9 +75,56 @@ TEST(QITETester, checkDeuteuronH2)
   const double finalEnergy = (*buffer)["opt-val"].as<double>();
   std::cout << "Final Energy: " << finalEnergy << "\n";
   EXPECT_NEAR(finalEnergy, -1.74886, 1e-3);
+}
+
+TEST(QITETester, checkDeuteuronH2Analytical) 
+{
+  auto qite = xacc::getService<xacc::Algorithm>("qite");
+  std::shared_ptr<xacc::Observable> observable = std::make_shared<xacc::quantum::PauliOperator>();
+  observable->fromString("5.907 - 2.1433 X0X1 - 2.1433 Y0Y1 + .21829 Z0 - 6.125 Z1");
+  auto acc = xacc::getAccelerator("qpp");
+  const int nbSteps = 10;
+  const double stepSize = 0.1;
+
+  const bool initOk =  qite->initialize({
+    std::make_pair("accelerator", acc),
+    std::make_pair("steps", nbSteps),
+    std::make_pair("observable", observable),
+    std::make_pair("step-size", stepSize),
+    std::make_pair("analytical", true),
+    std::make_pair("initial-state", 1)
+  });
+
+  EXPECT_TRUE(initOk);
+  auto buffer = xacc::qalloc(2);
+  qite->execute(buffer);
+  const double finalEnergy = (*buffer)["opt-val"].as<double>();
+  std::cout << "Final Energy: " << finalEnergy << "\n";
+  EXPECT_NEAR(finalEnergy, -1.74886, 1e-3);
   buffer->print();
+  const std::string Aop_Analytical = (*buffer)["A-op"].as<std::string>();
+  std::shared_ptr<xacc::quantum::PauliOperator> aOpObservable = std::make_shared<xacc::quantum::PauliOperator>();
+  aOpObservable->fromString(Aop_Analytical);
+  // Expected: recover UCC ansatz: "X0Y1 − X1Y0"
+  EXPECT_EQ(aOpObservable->getSubTerms().size(), 2);
+  auto term1 = aOpObservable->getSubTerms()[0];
+  auto term2 = aOpObservable->getSubTerms()[1];
+  const auto term1Coeff = term1->coefficient();
+  const auto term2Coeff = term2->coefficient();
+  // Two terms have opposite signs
+  const auto diff = term1Coeff + term2Coeff;
+  EXPECT_NEAR(std::norm(diff), 0.0, 1e-6);
+  const auto endsWith = [](const std::string& in_string, const std::string& in_ending) -> bool {
+    return std::equal(in_ending.rbegin(), in_ending.rend(), in_string.rbegin());
+  };
+
+  // two UCC's terms: "Y0 X1" "X0 Y1"
+  const bool check1 = endsWith(term1->toString(), "Y0 X1") || endsWith(term1->toString(), "X0 Y1");
+  const bool check2 = endsWith(term2->toString(), "Y0 X1") || endsWith(term2->toString(), "X0 Y1");
+  EXPECT_TRUE(check1 && check2);
 }
 
+
 TEST(QITETester, checkDeuteuronH3) 
 {
   auto qite = xacc::getService<xacc::Algorithm>("qite");
@@ -82,12 +133,16 @@ TEST(QITETester, checkDeuteuronH3)
   auto acc = xacc::getAccelerator("qpp");
   const int nbSteps = 10;
   const double stepSize = 0.1;
+  auto compiler = xacc::getCompiler("xasm");
+  auto ir = compiler->compile(R"(__qpu__ void f1(qbit q) { X(q[0]); })", nullptr);
+  auto x = ir->getComposite("f1");
 
   const bool initOk =  qite->initialize({
     std::make_pair("accelerator", acc),
     std::make_pair("steps", nbSteps),
     std::make_pair("observable", observable),
-    std::make_pair("step-size", stepSize)
+    std::make_pair("step-size", stepSize),
+    std::make_pair("ansatz", x)
   });
 
   EXPECT_TRUE(initOk);
@@ -96,7 +151,7 @@ TEST(QITETester, checkDeuteuronH3)
   qite->execute(buffer);
   const double finalEnergy = (*buffer)["opt-val"].as<double>();
   std::cout << "Final Energy: " << finalEnergy << "\n";
-  EXPECT_NEAR(finalEnergy, -2.04482, 1e-3);
+  EXPECT_NEAR(finalEnergy, -2.04482, 0.01);
 }
 
 int main(int argc, char **argv) {