adapt_vqe.cpp 13.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/*******************************************************************************
 * 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>
26
27
#include <sstream>
#include <string>
28
29
30
31
32
33
34
35
36
37

using namespace xacc;
using namespace xacc::quantum;

namespace xacc {
namespace algorithm {

bool ADAPT_VQE::initialize(const HeterogeneousMap &parameters) {

  if (!parameters.pointerLikeExists<Observable>("observable")) {
38
    xacc::info("Obs was false\n");
39
    return false;
40
41
42
43
  } 
  
  if (!parameters.pointerLikeExists<Accelerator>("accelerator")) {
    xacc::info("Acc was false\n");
44
    return false;
45
46
47
  } 
  
  if(!parameters.stringExists("pool")){
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    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);

  }

99
100
  if (parameters.keyExists<std::vector<double>>("checkpoint-parameters")) {
    checkpointParams= parameters.get<std::vector<double>>("checkpoint-parameters");
101
102
  }

103
104
  if (parameters.keyExists<std::vector<int>>("checkpoint-ops")) {
    checkpointOps = parameters.get<std::vector<int>>("checkpoint-ops");
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
  }

  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");
127
  std::stringstream ss;
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

  // 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
158
  if (!checkpointOps.empty()){
159

160
161
    if (!checkpointParams.empty()){
      x = checkpointParams;
162
    } else {
163
      x.resize(checkpointOps.size());
164
165
    }

166
167
    initialIter = checkpointOps.size();
    for (int i = 0; i < checkpointOps.size(); i++){
168
169
170
171
172

      auto exp_i_theta = std::dynamic_pointer_cast<quantum::Circuit>(
          xacc::getService<Instruction>("exp_i_theta"));

      exp_i_theta->expand(
173
          {std::make_pair("pauli", pauliOps[checkpointOps[i]]->toString()),
174
175
176
177
178
179
180
181
182
183
184
          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(),
185
                  {std::make_pair(optimizer->name() + "-optimizer", optimizer->get_algorithm()),
186
187
188
189
190
191
192
193
194
                  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];
195
196
197
198

    ss << std::setprecision(12) << oldEnergy << "\n";
    xacc::info(ss.str());
    ss.str(std::string());
199
200
201

  }

202
203
  xacc::info("Operator pool: " + operatorPool->name() + "\n");
  xacc::info("Number of operators in the pool: " + std::to_string(pauliOps.size()) + "\n\n");
204
205
206
207

  // start ADAPT loop
  for (int iter = initialIter; iter < _maxIter; iter++){

208
209
210
    xacc::info("Iteration: " + std::to_string(iter + 1) + "\n");
    xacc::info("Computing [H, A]\n\n");
    xacc::info("Printing commutators with absolute value above " + std::to_string(_printThreshold) + "\n");
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

    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){
231
232
233
        ss << std::setprecision(12) << "[H," << operatorIdx << "] = " << commutatorValue << "\n";
        xacc::info(ss.str());
        ss.str(std::string());
234
235
236
237
238
239
240
241
242
243
244
      }

      // update maxCommutator
      if(abs(commutatorValue) > abs(maxCommutator)){
        maxCommutatorIdx = operatorIdx;
        maxCommutator = commutatorValue;
      }

      gradientNorm += commutatorValue * commutatorValue;
      }
    }
245
246
247
248
249
    
    ss << std::setprecision(12) << "Max gradient component: [H, " 
        << maxCommutatorIdx << "] = " << maxCommutator << " a.u.\n";
    xacc::info(ss.str());
    ss.str(std::string());
250
251

    gradientNorm = std::sqrt(gradientNorm);
252
253
254
    ss << std::setprecision(12) << "Norm of gradient vector: " << gradientNorm << " a.u.\n";
    xacc::info(ss.str());
    ss.str(std::string());
255
256

    if (gradientNorm < _gradThreshold) { // ADAPT-VQE converged
257
258
259
260
261
262
      xacc::info("\nADAPT-VQE converged in " + std::to_string(iter) + " iterations.\n");
      ss << std::setprecision(12) << "ADAPT-VQE energy: " << oldEnergy << " a.u.\n";
      xacc::info(ss.str());
      ss.str(std::string());

      ss << "Optimal parameters: \n";
263
      for (auto param : x){
264
        ss << std::setprecision(12) << param << " ";
265
      }
266
267
268
      xacc::info(ss.str() + "\n");
      ss.str(std::string());

269
270
271
272
      return; 

    } else if (iter < _maxIter) { // Add operator and reoptimize 

273
      xacc::info("\nVQE optimization of current ansatz.\n\n");
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314

      // 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(),
315
                      {std::make_pair(optimizer->name() + "-optimizer", optimizer->get_algorithm()),
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
                      std::make_pair("initial-parameters", x)});
      } else {
        newOptimizer = optimizer;
      }

      // Start VQE optimization 
      auto sub_vqe = xacc::getAlgorithm(
          "vqe", {std::make_pair("observable", observable),
                  std::make_pair("optimizer", newOptimizer),
                  std::make_pair("accelerator", accelerator),
                  std::make_pair("gradient_strategy", gradientStrategy),
                  std::make_pair("ansatz", ansatzInstructions)
                  });
      sub_vqe->execute(buffer);

      auto newEnergy = (*buffer)["opt-val"].as<double>();
      x = (*buffer)["opt-params"].as<std::vector<double>>();
      oldEnergy = newEnergy;

335
336
337
338
339
      ss << std::setprecision(12) << "\nEnergy at ADAPT iteration " << iter + 1 << ": " << newEnergy << "\n";
      xacc::info(ss.str());
      ss.str(std::string());

      ss << std::setprecision(12) << "Parameters at ADAPT iteration " << iter + 1 << ": \n";
340
      for (auto param : x){
341
        ss << param << " ";
342
      }
343
344
      xacc::info(ss.str() + "\n");
      ss.str(std::string());
345

346
      ss << "Ansatz at ADAPT iteration " << std::to_string(iter + 1) << ": \n";
347
      for (auto op : ansatzOps){
348
        ss << op << " ";
349
      }
350
351
352
      xacc::info(ss.str() + "\n");
      ss.str(std::string());

353
      if(_printOps){
354
        xacc::info("Printing operators at ADAPT iteration " + std::to_string(iter + 1) + "\n\n");
355
        for(auto op : ansatzOps){
356
357
          xacc::info("Operator index: " + std::to_string(op) + "\n");
          xacc::info(operators[op]->toString() + "\n\n");
358
359
360
361
        }
      }

    } else {
362
      xacc::info("ADAPT-VQE did not converge in " + std::to_string(_maxIter) + " iterations.\n");
363
364
365
366
367
368
369
370
      return;
    }

  }
}

} // namespace adapt_vqe
} // namespace xacc