Merge pull request #204 from tnguyen-ornl/tnguyen/qite

Implements Quantum Imaginary Time Evolution algorithm

docs/source/extensions.rst

@@ -1220,6 +1220,122 @@ or in Python
    # i.e. phase value of 1/2^3 = 1/8.
    print(buffer)
 
+QITE
+++++
+The Quantum Imaginary Time Evolution (QITE) Algorithm requires the following input information:
+(`Motta et al. (2020) <https://arxiv.org/pdf/1901.07653.pdf>`_)
+
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|  Algorithm Parameter   |                  Parameter Description                          |             type                     |
++========================+=================================================================+======================================+
+|    observable          | The hermitian operator represents the cost Hamiltonian.         | std::shared_ptr<Observable>          |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|    accelerator         | The Accelerator backend to target                               | std::shared_ptr<Accelerator>         |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|    steps               | The number of Trotter steps.                                    | int                                  |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|    step-size           | The Trotter step size.                                          | double                               |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+
+Optionally, users can provide these parameters:
+
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|  Algorithm Parameter   |                  Parameter Description                          |             type                     |
++========================+=================================================================+======================================+
+|    ansatz              | State preparation circuit.                                      | std::shared_ptr<CompositeInstruction>|
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|    analytical          | If true, perform an analytical run rather than                  | boolean                              |
+|                        | executing quantum circuits on the Accelerator backend.          |                                      |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+|    initial-state       | For `analytical` mode only, select the initial state.           | int                                  |
++------------------------+-----------------------------------------------------------------+--------------------------------------+
+
+This Algorithm will add ``opt-val`` (``double``) which is the energy value at the final Trotter step to the provided ``AcceleratorBuffer``.
+The results of the algorithm are therefore retrieved via these keys (see snippet below). 
+Also, energy values at each Trotter step are stored in the ``exp-vals`` field (``vector<double>``).
+
+Note: during execution, the following line may be logged to the output console:
+
+.. code:: cpp
+
+   warning: solve(): system seems singular; attempting approx solution
+
+This is completely normal and can be safely ignored.
+
+.. code:: cpp
+
+   #include "xacc.hpp"
+   #include "xacc_observable.hpp"
+   #include "xacc_service.hpp"
+
+   int main(int argc, char **argv) {
+      xacc::Initialize(argc, argv);
+      // Use the Qpp simulator as the accelerator
+      auto acc = xacc::getAccelerator("qpp");
+      
+      auto buffer = xacc::qalloc(1);
+      auto observable = xacc::quantum::getObservable(
+            "pauli",
+            std::string("0.7071067811865475 X0 + 0.7071067811865475 Z0"));
+      
+      auto qite = xacc::getService<xacc::Algorithm>("qite");
+      const int nbSteps = 25;
+      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)
+      });
+
+      qite->execute(buffer);
+      std::cout << "Min Energy: " << (*buffer)["opt-val"].as<double>() << "\n";
+   }
+
+In Python:
+
+.. code:: python
+
+   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 Hamiltonian as an XACC PauliOperator
+   ham = xacc.getObservable('pauli', '0.70710678118 X0 + 0.70710678118 Z0')
+
+   # We just need 1 qubit
+   buffer = xacc.qalloc(1)
+
+   # Horizontal axis: 0 -> 2.5
+   # The number of Trotter steps 
+   nbSteps = 25
+
+   # The Trotter step size
+   stepSize = 0.1
+
+   # Create the QITE algorithm
+   qite = xacc.getAlgorithm('qite', {
+                           'accelerator': qpu,
+                           'observable': ham,
+                           'step-size': stepSize,
+                           'steps': nbSteps
+                           })
+
+   result = qite.execute(buffer)
+
+   # Expected result: ~ -1
+   print('Min energy value = ', buffer.getInformation('opt-val'))
+   E = buffer.getInformation('exp-vals')
+   # Plot energy vs. beta
+   plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro', label = 'XACC QITE')
+   plt.grid()
+   plt.show()
+
+
 Accelerator Decorators
 ----------------------
 ROErrorDecorator

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

python/examples/qite_example.py

+# Simple 1-qubit demonstration of the Quatum Imaginary Time Evolution algorithm
+# Reference: https://www.nature.com/articles/s41567-019-0704-4
+# Target H = 1/sqrt(2)(X + Z)
+# Expected minimum value: -1.0
+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 Hamiltonian as an XACC PauliOperator
+ham = xacc.getObservable('pauli', '0.70710678118 X0 + 0.70710678118 Z0')
+
+# We just need 1 qubit
+buffer = xacc.qalloc(1)
+
+# See Fig. 2 (e) of https://www.nature.com/articles/s41567-019-0704-4
+# Horizontal axis: 0 -> 2.5
+# The number of Trotter steps 
+nbSteps = 25
+
+# The Trotter step size
+stepSize = 0.1
+
+# Create the QITE algorithm
+qite = xacc.getAlgorithm('qite', {
+                        'accelerator': qpu,
+                        'observable': ham,
+                        'step-size': stepSize,
+                        'steps': nbSteps
+                        })
+
+result = qite.execute(buffer)
+
+# Expected result: ~ -1
+print('Min energy value = ', buffer.getInformation('opt-val'))
+E = buffer.getInformation('exp-vals')
+# Reproduce Fig. 2(e) of https://www.nature.com/articles/s41567-019-0704-4
+plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro', label = 'XACC QITE')
+plt.grid()
+plt.show()
\ No newline at end of file

quantum/plugins/algorithms/CMakeLists.txt

@@ -19,6 +19,7 @@ add_subdirectory(qpt)
 add_subdirectory(qaoa)
 add_subdirectory(qpe)
 add_subdirectory(gradient_strategies)
+add_subdirectory(qite)
 
 file(GLOB PYDECORATORS ${CMAKE_CURRENT_SOURCE_DIR}/vqe/python/*.py
                        ${CMAKE_CURRENT_SOURCE_DIR}/ml/ddcl/python/*.py)

quantum/plugins/algorithms/qite/CMakeLists.txt

+# *******************************************************************************
+# Copyright (c) 2019 UT-Battelle, LLC.
+# All rights reserved. This program and the accompanying materials
+# are made available under the terms of the Eclipse Public License v1.0
+# and Eclipse Distribution License v.10 which accompany this distribution.
+# The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
+# and the Eclipse Distribution License is available at
+# https://eclipse.org/org/documents/edl-v10.php
+#
+# Contributors:
+#   Thien Nguyen - initial API and implementation
+# *******************************************************************************/
+set(LIBRARY_NAME xacc-algorithm-qite)
+
+file(GLOB SRC *.cpp)
+
+usfunctiongetresourcesource(TARGET ${LIBRARY_NAME} OUT SRC)
+usfunctiongeneratebundleinit(TARGET ${LIBRARY_NAME} OUT SRC)
+
+add_library(${LIBRARY_NAME} SHARED ${SRC})
+
+target_include_directories(
+  ${LIBRARY_NAME}
+  PUBLIC .
+  ${CMAKE_SOURCE_DIR}/tpls/armadillo)
+
+target_link_libraries(${LIBRARY_NAME} PUBLIC xacc CppMicroServices xacc-quantum-gate)
+
+set(_bundle_name xacc_algorithm_qite)
+set_target_properties(${LIBRARY_NAME}
+                      PROPERTIES COMPILE_DEFINITIONS
+                                 US_BUNDLE_NAME=${_bundle_name}
+                                 US_BUNDLE_NAME
+                                 ${_bundle_name})
+
+usfunctionembedresources(TARGET
+                         ${LIBRARY_NAME}
+                         WORKING_DIRECTORY
+                         ${CMAKE_CURRENT_SOURCE_DIR}
+                         FILES
+                         manifest.json)
+
+if(APPLE)
+  set_target_properties(${LIBRARY_NAME}
+                        PROPERTIES INSTALL_RPATH "@loader_path/../lib")
+  set_target_properties(${LIBRARY_NAME}
+                        PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
+else()
+  set_target_properties(${LIBRARY_NAME}
+                        PROPERTIES INSTALL_RPATH "$ORIGIN/../lib")
+  set_target_properties(${LIBRARY_NAME} PROPERTIES LINK_FLAGS "-shared")
+  # Armadillo solver needs LAPACK
+  find_package(LAPACK)
+  if(LAPACK_FOUND)
+   target_link_libraries(${LIBRARY_NAME} PRIVATE ${LAPACK_LIBRARIES})
+  else()
+    message(STATUS "LAPACK NOT FOUND. QITE plugin may not work.")
+  endif()
+endif()
+
+if(XACC_BUILD_TESTS)
+  add_subdirectory(tests)
+endif()
+
+install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)

quantum/plugins/algorithms/qite/QITEActivator.cpp

+/*******************************************************************************
+ * Copyright (c) 2019 UT-Battelle, LLC.
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v1.0 which accompanies this
+ * distribution. The Eclipse Public License is available at
+ * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
+ *License is available at https://eclipse.org/org/documents/edl-v10.php
+ *
+ * Contributors:
+ *   Thien Nguyen - initial API and implementation
+ *******************************************************************************/
+#include "qite.hpp"
+
+#include "cppmicroservices/BundleActivator.h"
+#include "cppmicroservices/BundleContext.h"
+#include "cppmicroservices/ServiceProperties.h"
+
+#include <memory>
+#include <set>
+
+using namespace cppmicroservices;
+
+namespace {
+
+/**
+ */
+class US_ABI_LOCAL QITEActivator : public BundleActivator {
+
+public:
+  QITEActivator() {}
+
+  /**
+   */
+  void Start(BundleContext context) {
+    auto c = std::make_shared<xacc::algorithm::QITE>();
+    context.RegisterService<xacc::Algorithm>(c);
+  }
+
+  /**
+   */
+  void Stop(BundleContext /*context*/) {}
+};
+
+} // namespace
+
+CPPMICROSERVICES_EXPORT_BUNDLE_ACTIVATOR(QITEActivator)

quantum/plugins/algorithms/qite/manifest.json

+{
+  "bundle.symbolic_name" : "xacc_algorithm_qite",
+  "bundle.activator" : true,
+  "bundle.name" : "XACC QITE Algorithm",
+  "bundle.description" : ""
+}

quantum/plugins/algorithms/qite/qite.cpp

+/*******************************************************************************
+ * Copyright (c) 2019 UT-Battelle, LLC.
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License v1.0
+ * and Eclipse Distribution License v1.0 which accompanies this
+ * distribution. The Eclipse Public License is available at
+ * http://www.eclipse.org/legal/epl-v10.html and the Eclipse Distribution
+ *License is available at https://eclipse.org/org/documents/edl-v10.php
+ *
+ * Contributors:
+ *   Thien Nguyen - initial API and implementation
+ *******************************************************************************/
+#include "qite.hpp"
+
+#include "Observable.hpp"
+#include "xacc.hpp"
+#include "xacc_service.hpp"
+#include "PauliOperator.hpp"
+#include "Circuit.hpp"
+#include <memory>
+#include <armadillo>
+#include <cassert>
+
+namespace {
+const std::complex<double> I{ 0.0, 1.0};
+int findMatchingPauliIndex(const std::vector<std::string>& in_OpList, const std::string& in_obsStr)
+{
+  // Returns true if the *operators* of the two terms are identical
+  // e.g. a Z0X1 and b Z0X1 -> true 
+  const auto comparePauliString = [](const std::string& in_a, const std::string& in_b) -> bool {
+    // Strip the coefficient part
+    auto opA = in_a.substr(in_a.find_last_of(")") + 1);
+    auto opB = in_b.substr(in_b.find_last_of(")") + 1);
+    opA.erase(std::remove(opA.begin(), opA.end(), ' '), opA.end()); 
+    opB.erase(std::remove(opB.begin(), opB.end(), ' '), opB.end()); 
+    return opA == opB;
+  };
+
+  for (int i = 0; i < in_OpList.size(); ++i)
+  {
+    std::shared_ptr<xacc::Observable> obs = std::make_shared<xacc::quantum::PauliOperator>();
+    const std::string pauliObsStr = "1.0 " + in_OpList[i];
+    obs->fromString(pauliObsStr);
+    
+    if (comparePauliString(obs->toString(), in_obsStr))
+    {
+      return i;
+    }
+  }
+  // Failed!
+  return -1;
+}
+
+// Project/flatten the target observable into the full list of all
+// possible Pauli operator combinations.
+// e.g. H = a X + b Z (1 qubit)
+// -> { 0.0, a, 0.0, b } (the ordering is I, X, Y, Z)
+std::vector<double> observableToVec(std::shared_ptr<xacc::Observable> in_observable, const std::vector<std::string>& in_pauliObsList) 
+{
+  std::vector<double> obsProjCoeffs(in_pauliObsList.size(), 0.0);
+  for (const auto& term: in_observable->getNonIdentitySubTerms())
+  {
+    const auto index = findMatchingPauliIndex(in_pauliObsList, term->toString());
+    assert(index >= 0);
+    obsProjCoeffs[index] = term->coefficient().real();
+  }
+  return obsProjCoeffs;
+};
+
+// Helper to generate all permutation of Pauli observables:
+// e.g.
+// 1-qubit: I, X, Y, Z
+// 2-qubit: II, IX, IY, IZ, XI, XX, XY, XZ, YI, YX, YY, YZ, ZI, ZX, ZY, ZZ
+std::vector<std::string> generatePauliPermutation(int in_nbQubits)
+{
+  assert(in_nbQubits > 0);
+  const int nbPermutations = std::pow(4, in_nbQubits);
+  std::vector<std::string> opsList;
+  opsList.reserve(nbPermutations);
+  
+  const std::vector<std::string> pauliOps { "X", "Y", "Z" };
+  const auto addQubitPauli = [&opsList, &pauliOps](int in_qubitIdx){
+    const auto currentOpListSize = opsList.size();
+    for (int i = 0; i < currentOpListSize; ++i)
+    {
+      auto& currentOp = opsList[i];
+      for (const auto& pauliOp : pauliOps)
+      {
+        const auto newOp = currentOp + pauliOp + std::to_string(in_qubitIdx);
+        opsList.emplace_back(newOp);
+      }
+    }
+  };
+  
+  opsList = { "", "X0", "Y0", "Z0" };
+  for (int i = 1; i < in_nbQubits; ++i) 
+  {
+    addQubitPauli(i);
+  }
+
+  assert(opsList.size() == nbPermutations);
+  std::sort(opsList.begin(), opsList.end());
+
+  return opsList;
+};
+
+arma::cx_mat createSMatrix(const std::vector<std::string>& in_pauliOps, const std::vector<double>& in_tomoExp) 
+{
+  const auto sMatDim = in_pauliOps.size();
+  arma::cx_mat S_Mat(sMatDim, sMatDim, arma::fill::zeros);
+  arma::cx_vec b_Vec(sMatDim, arma::fill::zeros); 
+  
+  const auto calcSmatEntry = [&](int in_row, int in_col) -> std::complex<double> {
+    // Map the tomography expectation to the S matrix
+    // S(i, j) = <psi|sigma_dagger(i)sigma(j)|psi>
+    // sigma_dagger(i)sigma(j) will produce another Pauli operator with an additional coefficient.
+    // e.g. sigma_x * sigma_y = i*sigma_z
+    const auto leftOp = "1.0 " + in_pauliOps[in_row];
+    const auto rightOp = "1.0 " + in_pauliOps[in_col];
+    xacc::quantum::PauliOperator left(leftOp);
+    xacc::quantum::PauliOperator right(rightOp);
+    auto product = left * right;
+    const auto index = findMatchingPauliIndex(in_pauliOps, product.toString());
+    return in_tomoExp[index]*product.coefficient();
+  };
+
+  // S matrix:
+  for (int i = 0; i < sMatDim; ++i)
+  {
+    for (int j = 0; j < sMatDim; ++j)
+    {
+      S_Mat(i, j) = calcSmatEntry(i, j);
+    }
+  }
+  return S_Mat;
+}
+}
+
+using namespace xacc;
+namespace xacc {
+namespace algorithm {
+bool QITE::initialize(const HeterogeneousMap &parameters) 
+{
+  bool initializeOk = true;
+  if (!parameters.pointerLikeExists<Accelerator>("accelerator")) 
+  {
+    std::cout << "'accelerator' is required.\n";
+    initializeOk = false;
+  }
+
+  if (!parameters.keyExists<int>("steps")) 
+  {
+    std::cout << "'steps' is required.\n";
+    initializeOk = false;
+  }
+
+  if (!parameters.keyExists<double>("step-size")) 
+  {
+    std::cout << "'step-size' is required.\n";
+    initializeOk = false;
+  }
+
+  if (!parameters.pointerLikeExists<Observable>("observable")) 
+  {
+    std::cout << "'observable' is required.\n";
+    initializeOk = false;
+  }
+  
+  if (initializeOk)
+  {
+    m_accelerator = xacc::as_shared_ptr(parameters.getPointerLike<Accelerator>("accelerator"));
+    m_nbSteps = parameters.get<int>("steps");
+    m_dBeta = parameters.get<double>("step-size");
+    m_observable = xacc::as_shared_ptr(parameters.getPointerLike<Observable>("observable"));
+  }
+
+  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");