Qiskit: Qiskit Aquaを使わないQAOAの実装

はじめに

私の最初の記事で，TSPをQAOAで解いてみたというものを投稿しました．しかし，これはqiskit aquaという自分で回路を設計しなくても可能ないわゆるモジュール化されたものを使用していました．
そこで今回はqiskit aquaを使わずにQAOAを使ってみたいと思います．まず，最初の段階としてmax cutを用います．正直，他の問題もハミルトニアンを変更するだけでなんとかなると思いますので，qiskitで量子プログラミングをしていく人たちの参考になればなと思います．
日本語で量子プログラミングを解説いているサイトは複数ありますが，どれも別の量子プログラミング言語を用いていることから，qiskitで記述している私の記事にも意味はあると考えています．

なお，今回の記事でもQAOA自体の解説は他の方にお任せするとしてqiskitのコードの解説をしていきたいと思います．

参考: Quantum Native Dojo

code

まずは全体のcodeを

# coding: utf-8

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute
from qiskit.aqua.utils import tensorproduct
from qiskit import BasicAer
from qiskit.quantum_info.analysis import average_data

from scipy.optimize import minimize
import numpy as np
import time


def classical_minimize(cost_func, initial_params, options, method='powell'):
    print('classical minimize is starting now... ')
    start_time = time.time()
    result = minimize(cost_func, initial_params, options=options, method=method)
    print('running time: {}'.format(time.time() - start_time))
    print('opt_cost: {}'.format(result.fun))
    print('opt_params: {}'.format(result.x))
    return result.x


class MAXCUT:

    # 初期値設定
    def __init__(self, n_qubits, weight_matrix, p=1, num_steps=1):
        self.n_qubits = n_qubits
        self.weight_matrix = weight_matrix
        self.P = p
        self.num_steps = num_steps
        self.edge_list = []
        self._make_edge_list()

    # edge list の作成 from weight matrix
    def _make_edge_list(self):
        for i in range(self.weight_matrix.shape[0]):
            for j in range(i+1, self.weight_matrix.shape[1]):
                if self.weight_matrix[i][j] != 0:
                    self.edge_list.append([i, j])

    '------------------------------------------------------------------------------------------------------------------'

    def Zi_Zj(self, q1, q2):
        I_mat = np.array([[1, 0], [0, 1]])
        Z_mat = np.array([[1, 0], [0, -1]])
        if q1 == 0 or q2 == 0:
            tensor = Z_mat
        else:
            tensor = I_mat
        for i in range(1, self.n_qubits):
            if i == q1 or i == q2:
                tensor = tensorproduct(tensor, Z_mat)
            else:
                tensor = tensorproduct(tensor, I_mat)
        return tensor

    def observable(self):
        obs = np.zeros((2**self.n_qubits, 2**self.n_qubits))
        for a_edge in self.edge_list:
            q1, q2 = a_edge
            obs = obs + 0.5 * self.Zi_Zj(q1, q2)
        return obs

    '------------------------------------------------------------------------------------------------------------------'

    # U_C(gamma)
    def add_U_C(self, qc, gamma):
        for a_edge in self.edge_list:
            q1, q2 = a_edge
            qc.cx(q1, q2)
            qc.rz(-2**gamma, q2)
            qc.cx(q1, q2)
        return qc

    # U_X(beta)
    def add_U_X(self, qc, beta):
        for i in range(self.n_qubits):
            qc.rx(-2*beta, i)
        return qc

    def QAOA_output_onelayer(self, params, run=False):
        beta, gamma = params
        qr = QuantumRegister(self.n_qubits)
        cr = ClassicalRegister(self.n_qubits)
        qc = QuantumCircuit(qr, cr)
        qc.h(range(self.n_qubits))
        qc = self.add_U_C(qc, gamma)
        qc = self.add_U_X(qc, beta)
        qc.measure(range(self.n_qubits), range(self.n_qubits))
        NUM_SHOTS = 10000
        seed = 1234
        backend = BasicAer.get_backend('qasm_simulator')
        results = execute(qc, backend, shots=NUM_SHOTS, seed_simulator=seed).result()
        counts = results.get_counts(qc)
        expectation = average_data(counts, self.observable())
        return expectation

    '------------------------------------------------------------------------------------------------------------------'

    def minimize(self):
        initial_params = np.array([0.1, 0.1])
        opt_params = classical_minimize(self.QAOA_output_onelayer, initial_params,
                                        options={'maxiter':500}, method='powell')
        return opt_params

    def run(self):
        opt_params = self.minimize()
        beta_opt, gamma_opt = opt_params
        qr = QuantumRegister(self.n_qubits)
        cr = ClassicalRegister(self.n_qubits)
        qc = QuantumCircuit(qr, cr)
        qc.h(range(self.n_qubits))
        qc = self.add_U_C(qc, gamma_opt)
        qc = self.add_U_X(qc, beta_opt)
        qc.measure(range(self.n_qubits), range(self.n_qubits))
        NUM_SHOTS = 10000
        seed = 1234
        backend = BasicAer.get_backend('qasm_simulator')
        results = execute(qc, backend, shots=NUM_SHOTS, seed_simulator=seed).result()
        print(results.get_counts())


if __name__ == '__main__':
    weight_matrix = np.array([[0, 1, 0, 1],
                             [1, 0, 1, 0],
                             [0, 1, 0, 1],
                             [1, 0, 1, 0]])
    cut = MAXCUT(4, weight_matrix, p=1)
    cut.run()

気になる関数だけをピックアップして解説していきます．

C(Z)の作成

    def Zi_Zj(self, q1, q2):
        I_mat = np.array([[1, 0], [0, 1]])
        Z_mat = np.array([[1, 0], [0, -1]])
        if q1 == 0 or q2 == 0:
            tensor = Z_mat
        else:
            tensor = I_mat
        for i in range(1, self.n_qubits):
            if i == q1 or i == q2:
                tensor = tensorproduct(tensor, Z_mat)
            else:
                tensor = tensorproduct(tensor, I_mat)
        return tensor

    def observable(self):
        obs = np.zeros((2**self.n_qubits, 2**self.n_qubits))
        for a_edge in self.edge_list:
            q1, q2 = a_edge
            obs = obs + 0.5 * self.Zi_Zj(q1, q2)
        return obs

<\beta,\gamma|C(Z)|\beta,\gamma>

を計算する際に必要なC(Z)をここで作成しています．
この際にテンソル積を計算する必要がありますが，qiskit.aqua.utilsのtensorproductで計算可能です．

期待値の計算

qiskitで期待値を計算するには

        results = execute(qc, backend, shots=NUM_SHOTS, seed_simulator=seed).result()
        counts = results.get_counts(qc)
        expectation = average_data(counts, self.observable())

の方法でできます．

結果

実行結果は以下の通りになりました．

{'0101': 2578, '1110': 171, '1101': 146, '1001': 793, '1111': 141, '0011': 815, '0111': 122, '0010': 159, '0000': 161, '0100': 170, '0001': 151, '0110': 802, '1000': 160, '1100': 811, '1010': 2682, '1011': 138}

値が大きい，0101と1010が解として取り出すことができます．