量子計算機用SDK,pytketを利用したVQE

今回はCambridge Quantum Computingから出ている量子計算機用アルゴリズム開発キット、pytketを紹介します。そしてそれを用いて変分量子アルゴリズムの一つであるVariational Quantum Eigensolver(VQE)の実行プログラムを解説します。このpytketは、blueqatと同様に様々なバックエンドを利用でき、珍しくHoneywellのマシンにもジョブ投入できます。しかしながら、このpytketはバージョンアップごとに仕様が二転三転するため、すぐに前の書き方が通用しなくなります。この記事においてはこれらの環境で計算した結果を載せます。

python環境:python3.9.10(64-bit)

pytketバージョン:pytket0.17.0, pytket-qiskit0.20.0

実行環境:Visual Studio 2017

pytketにおいては、量子回路は次のように記述します。

from pytket import Circuit
from pytket.extensions.qiskit.backends import AerBackend, IBMQBackend, AerStateBackend

このバージョンにおいては、バックエンドは次のように記述し、回路の実行結果をこのように取り出します。

circ = Circuit(2, 2)
circ.Rx(0.3, 0).Ry(0.5, 1).Rz(-0.6, 1).measure_all()

b = AerBackend()
bs = AerStateBackend()
b.compile_circuit(circ)        # performs the minimal compilation to satisfy the device/simulator constraints
handle = b.process_circuit(circ, 1024)   # run the circuit 1024 times
bs.compile_circuit(circ)
from pytket.utils import probs_from_counts
print(b.get_result(handle))
print(bs)
state = bs.process_circuit(circ).get_state()
nvector = state.round(16)
print(probs_from_counts(b.get_result(handle).get_counts())) # retrieve and return the readouts
print(nvector)

バックエンドを変えると様々な量子計算機やそのシミュレーターで計算が可能です。しかし、残念ながらblueqatのようにそれを介さない計算は出来ません。ここからはVQEの実行法を解説します。もし、VQEの方法がわからないなら、こちら を参照してください。下記のソースはGitHubにおいて公式から配布されているソースコード ですが、これにおいて計算するのは原子間距離を0.75(Å)とした水素分子の基底エネルギーです。ハミルトニアンはopenfermionのものを使い、基底はSTO-3Gです。クラスターはUCCSDのものを使い、変数は一電子励起2個と二電子励起1個の3つです。それらを回路に登録しているのが120行目です。122行目がエネルギーを回路から計算する部分とクーロン斥力ポテンシャルの和を返す部分です。

 
from pytket import Circuit
from pytket.extensions.qiskit.backends import AerBackend, IBMQBackend, AerStateBackend
import sys
from openfermion import *
from scipy import *
from scipy.optimize import minimize as scipy_minimizer
from scipy import interpolate
import numpy as np
import numpy as np
import copy
from collections import Counter, defaultdict
from functools import reduce
import itertools
import random
import warnings
import time
import csv


# # VQE for Unitary Coupled Cluster using tket


# In this tutorial, we will focus on:
# - building parameterised ansätze for variational algorithms;
# - compilation tools for UCC-style ansätze.


# This example assumes the reader is familiar with the Variational Quantum Eigensolver and its application to electronic structure problems through the Unitary Coupled Cluster approach.
#
# To run this example, you will need [object Object] and [object Object], as well as [object Object], [object Object], and [object Object].
#
# We will start with a basic implementation and then gradually modify it to make it faster, more general, and less noisy. The final solution is given in full at the bottom of the notebook.
#
# Suppose we have some electronic configuration problem, expressed via a physical Hamiltonian. (The Hamiltonian and excitations in this example were obtained using [object Object] version 0.5.2 and [object Object] for H2, bond length 0.75A, sto3g basis, Jordan-Wigner encoding, with no qubit reduction or orbital freezing.)


from openfermion import QubitOperator


from openfermion import QubitOperator
from scipy.optimize import minimize
from sympy import symbols


from pytket.extensions.qiskit import AerBackend
from pytket.circuit import Circuit, Qubit
from pytket.partition import PauliPartitionStrat
from pytket.passes import GuidedPauliSimp, FullPeepholeOptimise
from pytket.pauli import Pauli, QubitPauliString
from pytket.utils import get_operator_expectation_value, gen_term_sequence_circuit
from pytket.utils.operators import QubitPauliOperator


# Obtain electronic Hamiltonian:


hamiltonian = (
    -0.8153001706270075 * QubitOperator("")
    + 0.16988452027940318 * QubitOperator("Z0")
    + -0.21886306781219608 * QubitOperator("Z1")
    + 0.16988452027940323 * QubitOperator("Z2")
    + -0.2188630678121961 * QubitOperator("Z3")
    + 0.12005143072546047 * QubitOperator("Z0 Z1")
    + 0.16821198673715723 * QubitOperator("Z0 Z2")
    + 0.16549431486978672 * QubitOperator("Z0 Z3")
    + 0.16549431486978672 * QubitOperator("Z1 Z2")
    + 0.1739537877649417 * QubitOperator("Z1 Z3")
    + 0.12005143072546047 * QubitOperator("Z2 Z3")
    + 0.04544288414432624 * QubitOperator("X0 X1 X2 X3")
    + 0.04544288414432624 * QubitOperator("X0 X1 Y2 Y3")
    + 0.04544288414432624 * QubitOperator("Y0 Y1 X2 X3")
    + 0.04544288414432624 * QubitOperator("Y0 Y1 Y2 Y3")
)
nuclear_repulsion_energy = 0.70556961456


hamiltonian_op = QubitPauliOperator.from_OpenFermion(hamiltonian)


# Obtain terms for single and double excitations:


q = [Qubit(i) for i in range(4)]
xyii = QubitPauliString([q[0], q[1]], [Pauli.X, Pauli.Y])
yxii = QubitPauliString([q[0], q[1]], [Pauli.Y, Pauli.X])
iixy = QubitPauliString([q[2], q[3]], [Pauli.X, Pauli.Y])
iiyx = QubitPauliString([q[2], q[3]], [Pauli.Y, Pauli.X])
xxxy = QubitPauliString(q, [Pauli.X, Pauli.X, Pauli.X, Pauli.Y])
xxyx = QubitPauliString(q, [Pauli.X, Pauli.X, Pauli.Y, Pauli.X])
xyxx = QubitPauliString(q, [Pauli.X, Pauli.Y, Pauli.X, Pauli.X])
yxxx = QubitPauliString(q, [Pauli.Y, Pauli.X, Pauli.X, Pauli.X])
yyyx = QubitPauliString(q, [Pauli.Y, Pauli.Y, Pauli.Y, Pauli.X])
yyxy = QubitPauliString(q, [Pauli.Y, Pauli.Y, Pauli.X, Pauli.Y])
yxyy = QubitPauliString(q, [Pauli.Y, Pauli.X, Pauli.Y, Pauli.Y])
xyyy = QubitPauliString(q, [Pauli.X, Pauli.Y, Pauli.Y, Pauli.Y])


# Symbolic UCC ansatz generation:


syms = symbols("p0 p1 p2")
singles_syms = {xyii: syms[0], yxii: -syms[0], iixy: syms[1], iiyx: -syms[1]}
doubles_syms = {
    xxxy: 0.25 * syms[2],
    xxyx: -0.25 * syms[2],
    xyxx: 0.25 * syms[2],
    yxxx: -0.25 * syms[2],
    yyyx: -0.25 * syms[2],
    yyxy: 0.25 * syms[2],
    yxyy: -0.25 * syms[2],
    xyyy: 0.25 * syms[2],
}
print(singles_syms, doubles_syms)
excitation_op = QubitPauliOperator({**singles_syms, **doubles_syms})
ucc_ref = Circuit(4).X(0).X(2)
ucc = gen_term_sequence_circuit(excitation_op, ucc_ref)


# Circuit simplification:


GuidedPauliSimp().apply(ucc)
FullPeepholeOptimise().apply(ucc)


# Connect to a simulator/device:


backend = AerBackend()


# Objective function:




def objective(params):
    circ = ucc.copy()
    sym_map = dict(zip(syms, params)) #assign parameters.
    circ.symbol_substitution(sym_map)
    return (
        get_operator_expectation_value(
            circ,
            hamiltonian_op,
            backend,
            n_shots=4000,
            partition_strat=PauliPartitionStrat.CommutingSets,
        )
        + nuclear_repulsion_energy
    ).real




# Optimise against the objective function:


initial_params = [1e-4, 1e-4, 4e-1]
result = minimize(objective, initial_params, method="Nelder-Mead")
print("Final parameter values", result.x)
print("Final energy value", result.fun)


# Exercises:
# - Replace the [object Object] call with its implementation and use this to pull the analysis for measurement reduction outside of the objective function, so our circuits can be fully determined and compiled once. This means that the [object Object] method will need to be applied to each measurement circuit instead of just the state preparation circuit.
# - Use the [object Object] class to add some mitigation of the measurement errors. Start by running the characterisation circuits first, before your main VQE loop, then apply the mitigation to each of the circuits run within the objective function.
# - Change the [object Object] by passing in a [object Object] [object Object] to simulate a noisy device. Compare the accuracy of the objective function both with and without the circuit simplification. Try running a classical optimiser over the objective function and compare the convergence rates with different noise models. If you have access to a QPU, try changing the [object Object] to connect to that and compare the results to the simulator.

これを実行した結果は大体-1.13(Hartree)になります。おまけにこのコードの最後には宿題もついてます。気が向いたらやってみるといいでしょう。今 後はこれを計算した結果を載せていきたいと思います。