QuantumNaturalGradient.cpp 14.9 KB
Newer Older
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
/*******************************************************************************
 * 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:
 *   Thien Nguyen - initial API and implementation
 *******************************************************************************/

#include "QuantumNaturalGradient.hpp"
15
#include "xacc.hpp"
16
#include "xacc_service.hpp"
17
#include <cassert>
18
#include <optional>
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
19
20
using namespace xacc;

21
22
23
24
25
namespace {
std::optional<xacc::quantum::PauliOperator> getGenerator(const InstPtr& in_inst) 
{
    if (in_inst->name() == "Rx")
    {
26
        return xacc::quantum::PauliOperator({{ in_inst->bits()[0], "X" }}, 0.5);
27
28
29
    }
    if (in_inst->name() == "Ry")
    {
30
        return xacc::quantum::PauliOperator({{ in_inst->bits()[0], "Y" }}, 0.5);
31
32
33
    }
    if (in_inst->name() == "Rz")
    {
34
        return xacc::quantum::PauliOperator({{ in_inst->bits()[0], "Z" }}, 0.5);
35
36
37
38
39
    }

    return std::nullopt;
}
}
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
40
41
namespace xacc {
namespace algorithm {
42
bool QuantumNaturalGradient::initialize(const HeterogeneousMap in_parameters)
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
43
{
44
    m_gradientStrategy.reset();
45
46
    m_layers.clear();
    m_nbMetricTensorKernels = 0;
47
    m_metricTermToIdx.clear();
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
    // User can provide a regular gradient strategy.
    // Note: this natural gradient requires a base gradient strategy.
    if (in_parameters.pointerLikeExists<AlgorithmGradientStrategy>("gradient-strategy"))
    {
        m_gradientStrategy = xacc::as_shared_ptr(in_parameters.getPointerLike<AlgorithmGradientStrategy>("gradient-strategy"));
    }

    // No base gradient strategy provided.
    // Must have observable specified.
    if (!m_gradientStrategy)
    {
        if (in_parameters.pointerLikeExists<Observable>("observable"))
        {
            auto observable = xacc::as_shared_ptr(in_parameters.getPointerLike<Observable>("observable"));
            // Create a default gradient strategy based on the observable.
63
            m_gradientStrategy = xacc::getService<AlgorithmGradientStrategy>("central");
64
65
66
67
68
69
70
71
72
            m_gradientStrategy->initialize({ std::make_pair("observable", observable)});
        }
        else
        {
            xacc::error("'observable' must be provided if no 'gradient-strategy' specified.");
            return false;
        }
    }

Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
73
74
75
76
77
    return true;
}

std::vector<std::shared_ptr<CompositeInstruction>> QuantumNaturalGradient::getGradientExecutions(std::shared_ptr<CompositeInstruction> in_circuit, const std::vector<double>& in_x) 
{
78
    auto baseGradientKernels = m_gradientStrategy->getGradientExecutions(in_circuit, in_x);
79
    m_nbMetricTensorKernels = 0;
80
    // Layering the circuit:
81
    m_layers = ParametrizedCircuitLayer::toParametrizedLayers(in_circuit);
82
    m_nbParams = in_x.size();
83
84
85
86
87
88
    std::vector<std::string> paramNames;
    assert(in_circuit->getVariables().size() == in_x.size());
    for (const auto& param : in_circuit->getVariables())
    {
        paramNames.emplace_back(param);
    }
89

90
91
    std::vector<std::shared_ptr<CompositeInstruction>>  metricTensorKernels;
    for (auto& layer : m_layers)
92
    {
93
94
95
96
97
98
99
100
101
102
103
104
105
        auto kernels = constructMetricTensorSubCircuit(layer, paramNames, in_x);
        
        // Insert *non-identity* kernels only
        size_t kernelIdx = 0;
        const auto isIdentityTerm = [](xacc::quantum::PauliOperator& in_pauli){
            return in_pauli.getNonIdentitySubTerms().size() == 0;
        };
        
        for (size_t i = 0; i < layer.kiTerms.size(); ++i)
        {
            if (!isIdentityTerm(layer.kiTerms[i]))
            {
                metricTensorKernels.emplace_back(kernels[kernelIdx]);
106
                m_metricTermToIdx.emplace(layer.kiTerms[i].toString(), metricTensorKernels.size() - 1);
107
108
109
110
111
112
113
114
115
116
117
118
119
            }
            kernelIdx++;
        }        
        
        for (size_t i = 0; i < layer.kikjTerms.size(); ++i)
        {
            if (!isIdentityTerm(layer.kikjTerms[i]))
            {
                metricTensorKernels.emplace_back(kernels[kernelIdx]);
                m_metricTermToIdx.emplace(layer.kikjTerms[i].toString(), metricTensorKernels.size());
            }
            kernelIdx++;
        }  
120
121
    }

122
123
    m_nbMetricTensorKernels = metricTensorKernels.size();
    baseGradientKernels.insert(baseGradientKernels.end(), metricTensorKernels.begin(), metricTensorKernels.end());
124
    return baseGradientKernels;
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
125
126
127
128
}

void QuantumNaturalGradient::compute(std::vector<double>& out_dx, std::vector<std::shared_ptr<AcceleratorBuffer>> in_results)
{
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    assert(m_nbMetricTensorKernels > 0 && in_results.size() > m_nbMetricTensorKernels);
    const auto iterPos = in_results.begin() + (in_results.size() - m_nbMetricTensorKernels);
    // Split the results: regular gradient results + metric tensor results.
    std::vector<std::shared_ptr<AcceleratorBuffer>> baseResults;
    baseResults.assign(in_results.begin(), iterPos);
    std::vector<std::shared_ptr<AcceleratorBuffer>> metricTensorResults;
    metricTensorResults.assign(iterPos, in_results.end());
    std::vector<double> rawDx = out_dx;
    // Calculate the raw gradients (using a regular gradient strategy)
    m_gradientStrategy->compute(rawDx, baseResults);
    // Solve the natural gradient equation:
    const auto gMat = constructMetricTensorMatrix(metricTensorResults);
    arma::dvec gradients(rawDx); 
    arma::dvec newGrads = arma::solve(gMat, gradients);
143
144
    // std::cout << "Regular gradients:\n" << gradients << "\n";
    // std::cout << "Natural gradients:\n" << newGrads << "\n";
145
    out_dx = arma::conv_to<std::vector<double>>::from(newGrads); 
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
146
}
147

148
ObservedKernels QuantumNaturalGradient::constructMetricTensorSubCircuit(ParametrizedCircuitLayer& io_layer, 
149
150
                                                    const std::vector<std::string>& in_varNames, 
                                                    const std::vector<double>& in_varVals) const
151
{
152
    assert(in_varNames.size() == in_varVals.size());
153
154
155
156
    // We need to observe all the generators of parametrized gates plus products of them.
    std::vector<xacc::quantum::PauliOperator> KiTerms;
    std::vector<xacc::quantum::PauliOperator> KiKjTerms;

157
    for (const auto& parOp: io_layer.ops)
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
    {
        auto opGenerator = getGenerator(parOp);
        assert(opGenerator.has_value());
        KiTerms.emplace_back(opGenerator.value());
    }

    for (const auto& genOp1 : KiTerms)
    {
        for (const auto& genOp2 : KiTerms)
        {
            auto kikjOp = genOp1 * genOp2;
            assert(kikjOp.nTerms() == 1);
            KiKjTerms.emplace_back(kikjOp);
        }
    }

    auto gateRegistry = xacc::getService<IRProvider>("quantum");
    auto circuitToObs = gateRegistry->createComposite("__LAYER__COMPOSITE__");
176
177
    std::vector<std::shared_ptr<xacc::CompositeInstruction>> obsComp;
    std::vector<double> resolvedParams;
178
    for (const auto& op : io_layer.preOps)
179
180
181
182
183
184
185
    {
        if (op->isParameterized() && op->getParameter(0).isVariable())
        {
            const auto varName = op->getParameter(0).as<std::string>();
            circuitToObs->addVariable(varName);
            const auto iter = std::find(in_varNames.begin(), in_varNames.end(), varName); 
            assert(iter != in_varNames.end());
186
187
            const size_t idx = std::distance(in_varNames.begin(), iter);
            assert(idx < in_varVals.size());
188
189
190
191
            resolvedParams.emplace_back(in_varVals[idx]);
        }
    }

192
    circuitToObs->addInstructions(io_layer.preOps);
193
194
    // Resolves the pre-ops now that we have passed that layer.
    auto resolvedCirc = resolvedParams.empty() ? circuitToObs :  circuitToObs->operator()(resolvedParams);
195

196
197
198
    io_layer.kiTerms = KiTerms;
    io_layer.kikjTerms = KiKjTerms;
    
199
    for (auto& term : KiTerms)
200
    {   
201
        auto obsKernels = term.observe(resolvedCirc);
202
203
204
205
206
207
        assert(obsKernels.size() == 1);
        obsComp.emplace_back(obsKernels[0]);
    }

    for (auto& term : KiKjTerms)
    {
208
        auto obsKernels = term.observe(resolvedCirc);
209
210
211
212
        assert(obsKernels.size() == 1);
        obsComp.emplace_back(obsKernels[0]);
    }

213
214
    const size_t NUM_KI_TERMS = io_layer.paramInds.size();
    const size_t NUM_KIKJ_TERMS = io_layer.paramInds.size() * io_layer.paramInds.size();
215
216
    // Validate the expected count.
    assert(obsComp.size() == NUM_KI_TERMS + NUM_KIKJ_TERMS);
217
    return obsComp;
218
219
}

220
arma::dmat QuantumNaturalGradient::constructMetricTensorMatrix(const std::vector<std::shared_ptr<xacc::AcceleratorBuffer>>& in_results)
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
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
    arma::dmat gMat(m_nbParams, m_nbParams, arma::fill::zeros);
    size_t blockIdx = 0;
    for (auto& layer : m_layers)
    {
        const auto nbParamsInBlock = layer.paramInds.size();
        // Constructs the block diagonal matrices
        arma::dmat blockMat(nbParamsInBlock, nbParamsInBlock, arma::fill::zeros);

        for (size_t i = 0; i < nbParamsInBlock; ++i)
        {
            for(size_t j = 0; j < nbParamsInBlock; ++j)
            {
                // Entry = <KiKj> - <Ki><Kj>
                // second_order_ev[i, j] - first_order_ev[i] * first_order_ev[j]
                auto firstOrderTerm1 = layer.kiTerms[i];
                const double factor1 = firstOrderTerm1.coefficient().real();
                auto firstOrderTerm2 = layer.kiTerms[j];
                const double factor2 = firstOrderTerm2.coefficient().real();
                auto secondOrderTerm = layer.kiTerms[i] * layer.kiTerms[j];
                
                const auto getExpectationForTerm = [&](xacc::quantum::PauliOperator& in_pauli){
                    if (m_metricTermToIdx.find(in_pauli.toString()) == m_metricTermToIdx.end())
                    {
                        return 1.0;
                    }
                    return in_results[m_metricTermToIdx[in_pauli.toString()]]->getExpectationValueZ();
                };

                const double firstOrderTerm1Exp = getExpectationForTerm(firstOrderTerm1);
                const double firstOrderTerm2Exp = getExpectationForTerm(firstOrderTerm2);
                const double secondOrderTermExp = getExpectationForTerm(secondOrderTerm);
                const double value = factor1*factor2*(secondOrderTermExp - firstOrderTerm1Exp * firstOrderTerm2Exp);
                blockMat(i,j) = value;
            }
        }

        for (size_t i = 0; i < nbParamsInBlock; ++i)
        {
            for(size_t j = 0; j < nbParamsInBlock; ++j)
            {
                gMat(blockIdx + i, blockIdx + j) =  blockMat(i,j);
            }
        }

        blockIdx += nbParamsInBlock;
    }

269
270
271
    return gMat;
}

272
273
274
275
276
277
278
279
280
281
std::vector<ParametrizedCircuitLayer> ParametrizedCircuitLayer::toParametrizedLayers(const std::shared_ptr<xacc::CompositeInstruction>& in_circuit)
{
    const auto variables = in_circuit->getVariables();
    if (variables.empty())
    {
        xacc::error("Input circuit is not parametrized.");
    }

    // Assemble parametrized circuit layer:
    std::vector<ParametrizedCircuitLayer> layers; 
282
    ParametrizedCircuitLayer currentLayer;
283
284
285
286
287
288
    std::set<size_t> qubitsInLayer;
    for (size_t instIdx = 0; instIdx < in_circuit->nInstructions(); ++instIdx)
    {
        const auto& inst = in_circuit->getInstruction(instIdx);
        if (!inst->isParameterized())
        {
289
290
291
292
293
294
295
296
            if (currentLayer.ops.empty())
            {
                currentLayer.preOps.emplace_back(inst);
            }
            else
            {
                currentLayer.postOps.emplace_back(inst);
            }
297
298
299
        }
        else
        {
300
            const auto& instParam = inst->getParameter(0);
301
302
303
304
305
306
307
308
309
            if (instParam.isVariable())
            {
                const auto varName = instParam.as<std::string>();
                const auto iter = std::find(variables.begin(), variables.end(), varName); 
                assert(iter != variables.end());
                assert(inst->bits().size() == 1);
                const auto bitIdx = inst->bits()[0];
                // This qubit line was already acted upon in this layer by a parametrized gate.
                // Hence, start a new layer.
310
                if (!currentLayer.postOps.empty() || xacc::container::contains(qubitsInLayer, bitIdx))
311
312
                {
                    assert(!currentLayer.ops.empty() && !currentLayer.paramInds.empty());
313
                    layers.emplace_back(std::move(currentLayer));
314
315
316
317
318
319
320
321
                    qubitsInLayer.clear();
                }
                currentLayer.ops.emplace_back(inst);
                currentLayer.paramInds.emplace_back(std::distance(iter, variables.begin()));
                qubitsInLayer.emplace(bitIdx);
            } 
            else
            {
322
323
324
325
326
327
328
329
                if (currentLayer.ops.empty())
                {
                    currentLayer.preOps.emplace_back(inst);
                }
                else
                {
                    currentLayer.postOps.emplace_back(inst);
                }       
330
331
332
333
            }
        }
    }

334
335
    assert(!currentLayer.ops.empty() && !currentLayer.paramInds.empty());
    layers.emplace_back(std::move(currentLayer));
336
    
337
338
339
    // Need to make a copy to do two rounds:
    auto layersCopied = layers;

340
341
342
343
344
345
346
    // Add pre-ops and post-ops:
    for (size_t i = 1; i < layers.size(); ++i)
    {
        // Pre-ops:
        auto& thisLayerPreOps = layers[i].preOps;
        const auto& previousLayerOps = layers[i - 1].ops;
        const auto& previousLayerPreOps = layers[i - 1].preOps;
347
348
349
350
351
352
        const auto& previousLayerPostOps = layers[i - 1].postOps;

        // Insert pre-ops, ops, and *raw* post-ops of the previous layer to the begining
        thisLayerPreOps.insert(thisLayerPreOps.begin(), previousLayerPostOps.begin(), previousLayerPostOps.end());
        thisLayerPreOps.insert(thisLayerPreOps.begin(), previousLayerOps.begin(), previousLayerOps.end());
        thisLayerPreOps.insert(thisLayerPreOps.begin(), previousLayerPreOps.begin(), previousLayerPreOps.end());
353
354
355
356
357
358
359
360
    }

    for (int i = layers.size() - 2; i >= 0; --i)
    {
        // Post-ops:
        auto& thisLayerPostOps = layers[i].postOps;
        const auto& nextLayerOps = layers[i + 1].ops;
        const auto& nextLayerPostOps = layers[i + 1].postOps;
361
362
363
364
365
366
        // Get the original pre-ops
        const auto& nextLayerPreOps = layersCopied[i + 1].preOps;

        // Insert next layer pre-ops, ops and its post op
        thisLayerPostOps.insert(thisLayerPostOps.end(), nextLayerPreOps.begin(), nextLayerPreOps.end());
        thisLayerPostOps.insert(thisLayerPostOps.end(), nextLayerOps.begin(), nextLayerOps.end());
367
368
369
370
371
372
373
        thisLayerPostOps.insert(thisLayerPostOps.end(), nextLayerPostOps.begin(), nextLayerPostOps.end());
    }

    // Verify:
    for (const auto& layer: layers)
    {
        assert(layer.preOps.size() + layer.ops.size() + layer.postOps.size() == in_circuit->nInstructions());
374
        assert(!layer.paramInds.empty() && !layer.ops.empty());
375
    }
376

377
    return layers;
378
} 
Nguyen, Thien Minh's avatar
Nguyen, Thien Minh committed
379
380
}
}