adapt_vqe.cpp 13.7 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

    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
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
      int nTermsCommutator = std::dynamic_pointer_cast<PauliOperator>(commutators[operatorIdx])->getTerms().size();
      if(nTermsCommutator != 0){

        // Print number of instructions for computing <observable>
        xacc::info("Number of instructions for commutator calculation: " 
                    + std::to_string(nTermsCommutator));
        // 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){
          ss << std::setprecision(12) << "[H," << operatorIdx << "] = " << commutatorValue << "\n";
          xacc::info(ss.str());
          ss.str(std::string());
        }

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

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

    gradientNorm = std::sqrt(gradientNorm);
257
258
259
    ss << std::setprecision(12) << "Norm of gradient vector: " << gradientNorm << " a.u.\n";
    xacc::info(ss.str());
    ss.str(std::string());
260
261

    if (gradientNorm < _gradThreshold) { // ADAPT-VQE converged
262
263
264
265
266
267
      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";
268
      for (auto param : x){
269
        ss << std::setprecision(12) << param << " ";
270
      }
271
272
273
      xacc::info(ss.str() + "\n");
      ss.str(std::string());

274
275
      xacc::info("Final ADAPT-VQE circuit\n" + ansatzInstructions->toString());

276
277
278
279
      return; 

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

280
      xacc::info("\nVQE optimization of current ansatz.\n\n");
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
315
316
317
318
319
320
321

      // 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(),
322
                      {std::make_pair(optimizer->name() + "-optimizer", optimizer->get_algorithm()),
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
                      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;

342
343
344
345
346
      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";
347
      for (auto param : x){
348
        ss << param << " ";
349
      }
350
351
      xacc::info(ss.str() + "\n");
      ss.str(std::string());
352

353
      ss << "Ansatz at ADAPT iteration " << std::to_string(iter + 1) << ": \n";
354
      for (auto op : ansatzOps){
355
        ss << op << " ";
356
      }
357
358
359
      xacc::info(ss.str() + "\n");
      ss.str(std::string());

360
      if(_printOps){
361
        xacc::info("Printing operators at ADAPT iteration " + std::to_string(iter + 1) + "\n\n");
362
        for(auto op : ansatzOps){
363
364
          xacc::info("Operator index: " + std::to_string(op) + "\n");
          xacc::info(operators[op]->toString() + "\n\n");
365
366
367
368
        }
      }

    } else {
369
      xacc::info("ADAPT-VQE did not converge in " + std::to_string(_maxIter) + " iterations.\n");
370
371
372
373
374
375
376
377
      return;
    }

  }
}

} // namespace adapt_vqe
} // namespace xacc