Implementation of ADAPT-VQE and Parameter Shift

Signed-off-by: Daniel Claudino <6d3@ornl.gov>

python/examples/adapt_vqe_h2.py

+import xacc
+
+qpu = xacc.getAccelerator('tnqvm')
+optimizer = xacc.getOptimizer('nlopt',{'nlopt-optimizer':'l-bfgs'})
+buffer = xacc.qalloc(4)
+
+geom = '''
+H  0.000000   0.0      0.0
+H   0.0        0.0  .7474
+'''
+H = xacc.getObservable('pyscf', {'basis': 'sto-3g', 'geometry': geom})
+print(H.toString())
+
+adapt = xacc.getAlgorithm('adapt-vqe', {'accelerator': qpu,
+                                  'optimizer': optimizer,
+                                  'observable': H,
+                                  'nElectrons': 2,
+                                  'gradient-strategy': 'parameter-shift-gradient',
+                                  'print-operators': True,
+                                  'threshold': 1.0e-2,
+                                  'print-threshold': 1.0e-10,
+                                  'maxiter': 2,
+                                  'pool': "uccsd"})
+# execute
+adapt.execute(buffer)

quantum/observable/fermion/FermionOperator.cpp

@@ -24,23 +24,26 @@ namespace quantum {
 FermionTerm &FermionTerm::operator*=(const FermionTerm &v) noexcept {
   coeff() *= std::get<0>(v);
 
-//   std::cout << "FermionTerm: " << id() << ", " << FermionTerm::id(std::get<1>(v)) << "\n";
+   //std::cout << "FermionTerm: " << id() << ", " << FermionTerm::id(std::get<1>(v)) << "\n";
   auto otherOps = std::get<1>(v);
   for (auto &kv : otherOps) {
     auto site = kv.first;
     auto c_or_a = kv.second;
-    // std::cout << "HELLO: " << site << ", " << std::boolalpha << c_or_a << "\n";
+    //std::cout << "\n\nHELLO: " << site << ", " << std::boolalpha << c_or_a << "\n";
     Operators o = ops();    
     if (!o.empty()) {
+
       auto it = std::find_if(o.begin(), o.end(),
                            [&](const std::pair<int, bool> &element) {
                              return element.first == site;
                            });
-    // std::cout << it->first << ", " << std::boolalpha << it->second << "\n";
+      //std::cout << it->first << ", " << std::boolalpha << it->second << "\n";
       if (it->first == site) {
-      if (it->second && c_or_a) {
+        if (it->second == c_or_a) {
         // zero out this FermionTerm
           ops().clear();
+        } else {//this adds the adjoint of operators whose sites are already in the product
+          ops().push_back({site, c_or_a});
         }
 
       } else {
@@ -104,6 +107,8 @@ const std::string FermionOperator::toString() {
     s << c << " " << std::get<2>(kv.second) << " ";
     Operators ops = std::get<1>(kv.second);
 
+    // I changed this to get the operators printed in the order they're built
+    /*
     std::vector<int> creations, annhilations;
     for (auto &t : ops) {
       if (t.second) {
@@ -112,18 +117,25 @@ const std::string FermionOperator::toString() {
         annhilations.push_back(t.first);
       }
     }
-
     for (auto &t : creations) {
       s << t << "^" << std::string(" ");
     }
     for (auto &t : annhilations) {
       s << t << std::string(" ");
+    }*/
+
+    for (auto &t : ops) {
+      if (t.second) {
+        s << t.first << "^" << std::string(" ");
+      }
+      else{
+        s << t.first << std::string(" ");
+      }
     }
 
     s << "+ ";
   }
 
-//   std::cout << "tostring " << s.str() << "\n";
   auto r = s.str().substr(0, s.str().size() - 2);
   xacc::trim(r);
   return r;
@@ -205,20 +217,31 @@ FermionOperator::operator*=(const FermionOperator &v) noexcept {
   std::unordered_map<std::string, FermionTerm> newTerms;
   for (auto &kv : terms) {
     for (auto &vkv : v.terms) {
-      auto multTerm = kv.second * vkv.second;
+
+      FermionTerm multTerm;
+      if (kv.first == "I"){
+        multTerm = vkv.second;
+      } else {
+        multTerm = kv.second * vkv.second;
+      }
+
       if (!multTerm.ops().empty()) {
         auto id = multTerm.id();
-
-        if (!newTerms.insert({id, multTerm}).second) {
+        if (!newTerms.insert({id, multTerm}).second && kv.first != "I") {
           newTerms.at(id).coeff() += multTerm.coeff();
         }
 
+        if (!newTerms.insert({id, multTerm}).second && kv.first == "I") {
+          newTerms.at(id).coeff() = std::get<0>(kv.second) * std::get<0>(vkv.second);
+        }
+
         if (std::abs(newTerms.at(id).coeff()) < 1e-12) {
           newTerms.erase(id);
         }
       }
     }
   }
+
   terms = newTerms;
   return *this;
 }
@@ -258,6 +281,13 @@ FermionOperator::operator*=(const std::complex<double> v) noexcept {
   return *this;
 }
 
+std::shared_ptr<Observable> FermionOperator::commutator(std::shared_ptr<Observable> op) {
+
+  FermionOperator& A = *std::dynamic_pointer_cast<FermionOperator>(op);
+  std::shared_ptr<FermionOperator> commutatorHA = std::make_shared<FermionOperator>((*this) * A - A * (*this));
+  return std::dynamic_pointer_cast<Observable>(commutatorHA);
+}
+
 } // namespace quantum
 } // namespace xacc
 

quantum/observable/fermion/FermionOperator.hpp

@@ -204,6 +204,9 @@ void fromOptions(const HeterogeneousMap& options) override {
       return;
   }
 
+  std::shared_ptr<Observable>
+  commutator(std::shared_ptr<Observable> obs) override;
+
 };
 
 } // namespace quantum

quantum/observable/fermion/tests/FermionOperatorTester.cpp

@@ -69,6 +69,30 @@ TEST(FermionOperatorTester, checkFromStr2) {
     FermionOperator op(src);
     std::cout << op.toString() << "\n";
 }
+
+TEST(FermionOperatorTester, checkMult) {
+    //FermionOperator op1("(-0.5,0)  3^ 2^ 0 1 + (0.5,-0)  1^ 0^ 2 3 + (0.5,-0)  0^ 1^ 3 2 + (0.5,0)  2^ 3^ 0 1 + (0.5,0)  3^ 3^ 1 1 + (-0.5,0)  1^ 1^ 3 3 + (-0.5,0)  1^ 0^ 3 2 + (-0.5,0)  2^ 3^ 1 0"), op2("(0.708024, 0)");
+    FermionOperator op1("3^ 2^ 3 2"), op2("0^ 3 1^ 2");
+
+    for (auto b : op1.getTerms()){std::cout << "First " << b.second.id() << "\n";}
+    for (auto b : op2.getTerms()){std::cout << "Second " << b.second.id() << "\n";}
+
+    FermionTerm mult;
+    for (auto a : op1.getTerms()){
+      for (auto b : op2.getTerms()){
+        mult = a.second * b.second;
+      }
+    }
+
+    for (auto k : mult.ops()){
+      std::cout << k.first << " " << k.second << "\n";
+    }
+
+    //auto mult = op1.getTerms().second * op2.getTerms().second;
+
+    std::cout << "Print op =" << (op1*op2).toString() << "\n\n";
+}
+
 int main(int argc, char** argv) {
     xacc::Initialize(argc,argv);
    ::testing::InitGoogleTest(&argc, argv);

quantum/observable/pauli/PauliOperator.cpp

@@ -760,6 +760,14 @@ void PauliOperator::fromXACCIR(std::shared_ptr<IR> ir) {
   }
 }
 
+std::shared_ptr<Observable> PauliOperator::commutator(std::shared_ptr<Observable> op) {
+
+  PauliOperator& A = *std::dynamic_pointer_cast<PauliOperator>(op);
+  std::shared_ptr<PauliOperator> commutatorHA =  std::make_shared<PauliOperator>((*this) * A - A * (*this));
+  return std::dynamic_pointer_cast<Observable>(commutatorHA);
+
+}
+
 } // namespace quantum
 } // namespace xacc
 

quantum/observable/pauli/PauliOperator.hpp

@@ -339,6 +339,9 @@ public:
   const std::string name() const override { return "pauli"; }
   const std::string description() const override { return ""; }
   void fromOptions(const HeterogeneousMap &options) override { return; }
+  
+  std::shared_ptr<Observable>
+  commutator(std::shared_ptr<Observable> obs) override;
 };
 } // namespace quantum
 

quantum/plugins/algorithms/CMakeLists.txt

@@ -10,6 +10,7 @@
 # Contributors:
 #   Alexander J. McCaskey - initial API and implementation
 # *******************************************************************************/
+add_subdirectory(adapt_vqe)
 add_subdirectory(vqe)
 add_subdirectory(rdm)
 add_subdirectory(ml)
@@ -17,6 +18,7 @@ add_subdirectory(rotoselect)
 add_subdirectory(qpt)
 add_subdirectory(qaoa)
 add_subdirectory(qpe)
+add_subdirectory(gradient_strategies)
 
 file(GLOB PYDECORATORS ${CMAKE_CURRENT_SOURCE_DIR}/vqe/python/*.py
                        ${CMAKE_CURRENT_SOURCE_DIR}/ml/ddcl/python/*.py)

quantum/plugins/algorithms/adapt_vqe/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:
+#   Alexander J. McCaskey - initial API and implementation
+# *******************************************************************************/
+set(LIBRARY_NAME xacc-algorithm-adapt_vqe)
+
+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 .)
+
+target_link_libraries(${LIBRARY_NAME} PUBLIC xacc CppMicroServices PRIVATE xacc-quantum-gate)
+
+set(_bundle_name xacc_algorithm_adapt_vqe)
+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")
+endif()
+
+if(XACC_BUILD_TESTS)
+  add_subdirectory(tests)
+endif()
+
+install(TARGETS ${LIBRARY_NAME} DESTINATION ${CMAKE_INSTALL_PREFIX}/plugins)
\ No newline at end of file

quantum/plugins/algorithms/adapt_vqe/adapt_vqe.cpp

+/*******************************************************************************
+ * Copyright (c) 2020 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:
+ *   Daniel Claudino - initial API and implementation
+ *******************************************************************************/
+
+#include "adapt_vqe.hpp"
+
+#include "FermionOperator.hpp"
+#include "ObservableTransform.hpp"
+#include "xacc.hpp"
+#include "xacc_service.hpp"
+#include "xacc_observable.hpp"
+#include "Circuit.hpp"
+#include "AlgorithmGradientStrategy.hpp"
+
+#include <memory>
+#include <iomanip>
+
+using namespace xacc;
+using namespace xacc::quantum;
+
+namespace xacc {
+namespace algorithm {
+
+bool ADAPT_VQE::initialize(const HeterogeneousMap &parameters) {
+
+  if (!parameters.pointerLikeExists<Observable>("observable")) {
+    std::cout << "Obs was false\n";
+    return false;
+  } else if (!parameters.pointerLikeExists<Accelerator>(
+                 "accelerator")) {
+    std::cout << "Acc was false\n";
+    return false;
+  } else if(!parameters.stringExists("pool")){
+    return false;
+  }
+
+  optimizer = parameters.get<std::shared_ptr<Optimizer>>("optimizer");
+  observable = parameters.get<std::shared_ptr<Observable>>("observable");
+  accelerator = parameters.get<std::shared_ptr<Accelerator>>("accelerator");
+  nElectrons = parameters.get<int>("nElectrons");
+
+  if (parameters.keyExists<int>("maxiter")) {
+    _maxIter = parameters.get<int>("maxiter");
+  } 
+
+  if (parameters.keyExists<double>("threshold")) {
+    _gradThreshold = parameters.get<double>("threshold");
+  }
+
+  if (parameters.keyExists<double>("print-threshold")) {
+    _printThreshold = parameters.get<double>("print-threshold");
+  }
+
+  if (parameters.keyExists<bool>("print-operators")) {
+    _printOps = parameters.get<bool>("print-operators");
+  }
+
+  pool = parameters.getString("pool");
+
+  // Check if Observable is Fermion or Pauli and manipulate accordingly
+  //
+  // if string has ^, it's FermionOperator
+  if (observable->toString().find("^") != std::string::npos){
+
+    auto jw = xacc::getService<ObservableTransform>("jw");
+    if (std::dynamic_pointer_cast<FermionOperator>(observable)) {
+      observable = jw->transform(observable);
+    } else {
+      auto fermionObservable = xacc::quantum::getObservable("fermion", observable->toString());
+      observable = jw->transform(std::dynamic_pointer_cast<Observable>(fermionObservable));      
+    }
+
+  // observable is PauliOperator, but does not cast down to it
+  // Not sure about the likelihood of this happening, but want to cover all bases
+  } else if (observable->toString().find("X") != std::string::npos
+            || observable->toString().find("Y") != std::string::npos
+            || observable->toString().find("Z") != std::string::npos
+            && !std::dynamic_pointer_cast<PauliOperator>(observable)){
+
+    auto pauliObservable = xacc::quantum::getObservable("pauli", observable->toString());
+    observable = std::dynamic_pointer_cast<Observable>(pauliObservable);
+
+  }
+
+  if (parameters.keyExists<std::vector<double>>("initial-parameters")) {
+    initialParameters= parameters.get<std::vector<double>>("initial-parameters");
+  }
+
+  if (parameters.keyExists<std::vector<int>>("initial-ansatz")) {
+    initialAnsatz = parameters.get<std::vector<int>>("initial-ansatz");
+  }
+
+  if (parameters.stringExists("gradient-strategy")) {
+    gradStrategyName = parameters.getString("gradient-strategy");
+  }
+
+  return true;
+}
+
+const std::vector<std::string> ADAPT_VQE::requiredParameters() const {
+  return {"observable", "optimizer", "accelerator", "nElectrons", "pool"};
+}
+
+void ADAPT_VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
+
+  auto ansatzRegistry = xacc::getIRProvider("quantum");
+  auto ansatzInstructions = ansatzRegistry->createComposite("ansatzCircuit");
+  auto operatorPool = xacc::getService<OperatorPool>(pool);
+  auto operators = operatorPool->generate(buffer->size(), nElectrons);
+  std::vector<std::shared_ptr<Observable>> pauliOps;
+  std::vector<int> ansatzOps;
+  auto jw = xacc::getService<ObservableTransform>("jw");
+
+  // Mean-field state
+  std::size_t j;
+  for (int i = 0; i < nElectrons/2; i++) {
+    j = (std::size_t) i;
+    auto alphaXGate = ansatzRegistry->createInstruction("X", std::vector<std::size_t>{j});
+    ansatzInstructions->addInstruction(alphaXGate);
+    j = (std::size_t) (i + buffer->size()/2);
+    auto betaXGate = ansatzRegistry->createInstruction("X", std::vector<std::size_t>{j});
+    ansatzInstructions->addInstruction(betaXGate);
+  }
+
+  // Vector of commutators, need to compute them only once
+  std::vector<std::shared_ptr<Observable>> commutators;
+  for (auto op : operators){
+    std::shared_ptr<Observable> pauliOp;
+    if(std::dynamic_pointer_cast<PauliOperator>(op)){
+      pauliOp = op;
+    } else {
+      pauliOp = jw->transform(op);
+    }
+    pauliOps.push_back(pauliOp);
+    commutators.push_back(observable->commutator(pauliOp));
+  }
+
+  int initialIter = 0;
+  double oldEnergy = 0.0;
+  std::vector<double> x; // these are the variational parameters
+
+  // Resume from a previously optimized ansatz
+  if (!initialAnsatz.empty()){
+
+    if (!initialParameters.empty()){
+      x = initialParameters;
+    } else {
+      x.resize(initialAnsatz.size());
+    }
+
+    initialIter = initialAnsatz.size();
+    for (int i = 0; i < initialAnsatz.size(); i++){
+
+      auto exp_i_theta = std::dynamic_pointer_cast<quantum::Circuit>(
+          xacc::getService<Instruction>("exp_i_theta"));
+
+      exp_i_theta->expand(
+          {std::make_pair("pauli", pauliOps[initialAnsatz[i]]->toString()),
+          std::make_pair("param_id", std::string("x") + std::to_string(i)),
+          std::make_pair("no-i", true)});
+
+      ansatzInstructions->addVariable(std::string("x") + std::to_string(i));
+
+      for (auto& inst : exp_i_theta->getInstructions()){  
+        ansatzInstructions->addInstruction(inst);
+      }
+    }
+
+    auto newOptimizer = xacc::getOptimizer(optimizer->name(),
+                  {std::make_pair(optimizer->name() + "-optimizer", optimizer->optimizer_name()),
+                  std::make_pair("initial-parameters", x)});
+
+    auto init_vqe = xacc::getAlgorithm(
+        "vqe", {std::make_pair("observable", observable),
+                std::make_pair("optimizer", optimizer),
+                std::make_pair("accelerator", accelerator),
+                std::make_pair("ansatz", ansatzInstructions)});
+    auto tmp_buffer = xacc::qalloc(buffer->size());
+    oldEnergy = init_vqe->execute(tmp_buffer, x)[0];
+    std::cout << "Initial energy = " << oldEnergy << "\n";
+
+  }
+
+  std::cout << "Operator pool: " << operatorPool->name() << "\n";
+  std::cout << "Number of operators in the pool: " << pauliOps.size() << "\n\n";
+
+  // start ADAPT loop
+  for (int iter = initialIter; iter < _maxIter; iter++){
+
+    std::cout << "Iteration: " << iter + 1 << "\n";
+    std::cout << "Computing [H, A]\n" << "\n";
+    std::cout << "Printing commutators with absolute value above " << _printThreshold << "\n";
+
+    int maxCommutatorIdx = 0;
+    double maxCommutator = 0.0;
+    double gradientNorm = 0.0;
+
+    // Loop over non-vanishing commutators and select the one with largest magnitude
+    for (int operatorIdx = 0; operatorIdx < commutators.size(); operatorIdx++){
+
+      // only compute commutators if they aren't zero
+      if(std::dynamic_pointer_cast<PauliOperator>(commutators[operatorIdx])->getTerms().size() != 0){
+      // observe the commutators with the updated circuit ansatz
+      auto grad_vqe = xacc::getAlgorithm(
+          "vqe", {std::make_pair("observable", commutators[operatorIdx]),
+                  std::make_pair("optimizer", optimizer),
+                  std::make_pair("accelerator", accelerator),
+                  std::make_pair("ansatz", ansatzInstructions)});
+      auto tmp_buffer = xacc::qalloc(buffer->size());
+      auto commutatorValue = std::real(grad_vqe->execute(tmp_buffer, x)[0]);
+
+      if(abs(commutatorValue) > _printThreshold){
+        std::cout << std::setprecision(12) << "[H," << operatorIdx << "] = " << commutatorValue << "\n";
+      }
+
+      // update maxCommutator
+      if(abs(commutatorValue) > abs(maxCommutator)){
+        maxCommutatorIdx = operatorIdx;
+        maxCommutator = commutatorValue;
+      }
+
+      gradientNorm += commutatorValue * commutatorValue;
+      }
+    }
+
+    gradientNorm = std::sqrt(gradientNorm);
+    std::cout << "Max gradient component: [H, " << maxCommutatorIdx << "] = " << maxCommutator << " a.u.\n";
+    std::cout << "Norm of gradient vector: " << gradientNorm << " a.u.\n";
+
+    if (gradientNorm < _gradThreshold) { // ADAPT-VQE converged
+      std::cout << "\nADAPT-VQE converged in " << iter << " iterations.\n";
+      std::cout << "ADAPT-VQE energy: " << oldEnergy << " a.u.\n";
+      std::cout << "Optimal parameters: \n";
+      for (auto param : x){
+        std::cout << param << " ";
+      }
+      std::cout << "\n";
+      return; 
+
+    } else if (iter < _maxIter) { // Add operator and reoptimize 
+
+      std::cout << "\nVQE optimization of current ansatz.\n\n";
+
+      // keep track of growing ansatz
+      ansatzOps.push_back(maxCommutatorIdx);
+
+      // Instruction service for the operator to be added to the ansatz
+      auto maxCommutatorGate = std::dynamic_pointer_cast<quantum::Circuit>(
+          xacc::getService<Instruction>("exp_i_theta"));
+
+      // Create instruction for new operator
+      maxCommutatorGate->expand(
+          {std::make_pair("pauli", pauliOps[maxCommutatorIdx]->toString()),
+          std::make_pair("param_id", std::string("x") + std::to_string(iter)),
+          std::make_pair("no-i", true)});
+
+      // Add instruction for new operator to the current ansatz
+      ansatzInstructions->addVariable(std::string("x") + std::to_string(iter));
+
+      // Append new instructions to current circuit
+      for (auto& inst : maxCommutatorGate->getInstructions()){  
+        ansatzInstructions->addInstruction(inst);
+      }
+
+      // Instantiate gradient class
+      std::shared_ptr<AlgorithmGradientStrategy> gradientStrategy;
+      if(!gradStrategyName.empty()){
+        gradientStrategy = xacc::getService<AlgorithmGradientStrategy>(gradStrategyName);
+        gradientStrategy->optionalParameters({std::make_pair("observable", observable),
+                                            std::make_pair("commutator", commutators[maxCommutatorIdx]),
+                                            std::make_pair("jw", true)
+                                            });
+      } else {
+        gradientStrategy = nullptr;
+      }
+
+      // If using gradient-based optimizer 
+      // Convergence is improved if passing initial parameters to optimizer
+      // so we create a new instance of Optimizer with them
+      std::shared_ptr<Optimizer> newOptimizer;
+      if(gradientStrategy){
+        x.insert(x.begin(), 0.0);
+        newOptimizer = xacc::getOptimizer(optimizer->name(),
+                      {std::make_pair(optimizer->name() + "-optimizer", optimizer->optimizer_name()),
+                      std::make_pair("initial-parameters", x)});
+      } else {
+        newOptimizer = optimizer;
+      }
+
+      // Start VQE optimization 
+      auto sub_vqe = xacc::getAlgorithm(
+          "vqe", {std::make_pair("observable", observable),