🎉

量子コンペやるぜ(IBM Quantum Challenge Fall 2021 4c回答編)

2021/11/19に公開

IBM Quantum Challenge Fall 2021

IBM Quantum Challenge Fall 2021略して、量コンをやったぜ。
青バッチをゲトしたけど、最終問題の4cは答えが合わなかったぜ。

Lukas Botschさんの回答

量コンのGithubページ見ているとPRが何件かあったので、その中の一つの回答を勝手に解説するぜ。
今日のパクリ先は、これだ。

QAOA

4cは、QAOAをライブラリを使わずに量子回路を自力で作って、そのコストを競う勝負だ。
まずは、QAOAが何かというのがわからないと話にならない。
4cに辿りつくまでにライブラリを使ったQAOAをやっているはずだが、中身がよくわからない。
そんな野郎はコレを見るとよい。
https://www.youtube.com/watch?v=nLY8pvcASnk

QRAM

QRAMを理解しておいた方がよい。
過去問(2020年)とかIMPLEMENTING QRAM IN QISKIT WITH CODEがわかりやすい。

全体のコード

まずはコードを全部のせてしまおう。
Lukas Botschさんのコードをちょっとだけ修正している。
引数を元々の引数の名前に直したのと、回路図を出力するようにしたのと、ローカルで実行できるようにしてある程度だ。

from typing import List, Union
import math
import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, assemble
from qiskit.circuit import Gate
from qiskit.circuit.library.standard_gates import *
#from qiskit.circuit.library import QFT
from qiskit.visualization import plot_histogram

from qiskit import execute, Aer

instance_examples = [
    {
        'L1': [3, 7, 3, 4, 2, 6, 2, 2, 4, 6, 6],
        'L2': [7, 8, 7, 6, 6, 9, 6, 7, 6, 7, 7],
        'C1': [2, 2, 2, 3, 2, 4, 2, 2, 2, 2, 2],
        'C2': [4, 3, 3, 4, 4, 5, 3, 4, 4, 3, 4],
        'C_max': 33
    },
    # {
    #     'L1': [4, 2, 2, 3, 5, 3, 6, 3, 8, 3, 2],
    #     'L2': [6, 5, 8, 5, 6, 6, 9, 7, 9, 5, 8],
    #     'C1': [3, 3, 2, 3, 4, 2, 2, 3, 4, 2, 2],
    #     'C2': [4, 4, 3, 5, 5, 3, 4, 5, 5, 3, 5],
    #     'C_max': 38
    # },
    # {
    #     'L1': [5, 4, 3, 3, 3, 7, 6, 4, 3, 5, 3],
    #     'L2': [9, 7, 5, 5, 7, 8, 8, 7, 5, 7, 9],
    #     'C1': [2, 2, 4, 2, 3, 4, 2, 2, 2, 2, 2],
    #     'C2': [3, 4, 5, 4, 4, 5, 3, 3, 5, 3, 5],
    #     'C_max': 35
    # }
]

def solver_function(L1: list, L2: list, C1: list, C2: list, C_max: int) -> QuantumCircuit:

    # the number of qubits representing answers
    index_qubits = len(L1)

    # the maximum possible total cost
    max_c = sum([max(l0, l1) for l0, l1 in zip(C1, C2)])

    # the number of qubits representing data values can be defined using the maximum possible total cost as follows:
    data_qubits = math.ceil(math.log(max_c, 2)) + 1 if not max_c & (max_c - 1) == 0 else math.ceil(math.log(max_c, 2)) + 2

    ### Phase Operator ###
    def phase_return(index_qubits: int, gamma: float, L1: list, L2: list, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qc = QuantumCircuit(qr_index)
        ### U_1(gamma * (lambda2 - lambda1)) for each qubit ###
        for i in range(len(L1)):
            qc.p(-gamma * (L2[i] - L1[i]) / 2, qr_index[i])

        qc.draw(output='mpl', filename="phase_return.png")
        return qc.to_gate(label=" Penalty Dephasing ") if to_gate else qc

    def qft(data_qubits: int, inverse=False, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qc = QuantumCircuit(data_qubits)
        N = data_qubits-1
        qc.h(N)
        for n in range(1, N+1):
            qc.cp(2*np.pi / 2**(n+1), N-n, N)

        for i in range(1, N):
            qc.h(N-i)
            for n in range(1, N-i+1):
                qc.cp(2*np.pi / 2**(n+1), N-(n+i), N-i)

        qc.h(0)

        # Calculate IQFT
        qc = qc.inverse() if inverse else qc

        qc.draw(output='mpl', filename="qft.png")
        return qc.to_gate(label=" QFT ") if to_gate and not inverse else qc.to_gate(label=" IQFT ") if to_gate else qc

    def qft_add(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qc = QuantumCircuit(data_qubits)
        N = data_qubits-1
        bin_const = [int(x) for x in bin(const)[2:]] # Big endian
        bin_const = [0]*(N-len(bin_const)) + bin_const
        for n in range(1, N+1):
            if bin_const[n-1]:
                qc.p(2*np.pi / 2**(n+1), N)

        for i in range(1, N+1):
            for n in range(N-i+1):
                if bin_const[n+i-1]:
                    qc.p(2*np.pi / 2**(n+1), N-i)

        qc.draw(output='mpl', filename="qft_add.png")
        return qc.to_gate(label=f" [ +{const} ] ") if to_gate else qc

    def cost_calculation(index_qubits: int, cost_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qr_data = QuantumRegister(data_qubits, "data")
        qc = QuantumCircuit(qr_index, qr_data)

        qc.append(qft(cost_qubits), qr_data)
        for i, (c1, c2) in enumerate(zip(list1, list2)):
            qc.append(qft_add(cost_qubits, c2).control(1), [qr_index[i]] + qr_data[:])
            qc.x(qr_index[i])
            qc.append(qft_add(cost_qubits, c1).control(1), [qr_index[i]] + qr_data[:])
            qc.x(qr_index[i])
        qc.append(qft(cost_qubits, inverse=True), qr_data)

        qc.draw(output='mpl', filename="cost_calculation.png")
        return qc.to_gate(label=" Cost Calculation ") if to_gate else qc

    def constraint_testing(data_qubits: int, C_max: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)

        c = data_qubits - 1
        w = 2**c - C_max - 1

        if w > 0:
            # Add w to cost(z)
            qc.append(qft(data_qubits), qr_data)
            qc.append(qft_add(data_qubits, w), qr_data)
            qc.append(qft(data_qubits, inverse=True), qr_data)
            # Set F if the condition cost(z) > C_max holds
            qc.cx(qr_data[c], qr_f)

        qc.draw(output='mpl', filename="constraint_testing.png")
        return qc.to_gate(label=" constraint_testing ") if to_gate else qc

    def penalty_dephasing(data_qubits: int, alpha: float, gamma: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)
        for i in range(data_qubits):
            qc.cp((2**i) * alpha * gamma, qr_f, qr_data[i])
        
        # Add phase -alpha * (Cmax + w) = -alpha * 2^(cost_qubits-1) to all
        # infeasible states
        C_max_plus_w = 2**(data_qubits-1)-1
        qc.p(-C_max_plus_w * alpha * gamma, qr_f)

        qc.draw(output='mpl', filename="penalty_dephasing.png")
        return qc.to_gate(label=" Penalty Dephasing ") if to_gate else qc

    def reinitialization(index_qubits: int, data_qubits: int, C1: list, C2: list, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_index, qr_data, qr_f)

        qc.append(constraint_testing(data_qubits, C_max, False).inverse().to_gate(label="  Constraint Testing +  "), qr_data[:] + qr_f[:])
        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2, False).inverse().to_gate(label="  Cost Calculation +  "), qr_index[:] + qr_data[:])

        qc.draw(output='mpl', filename="reinitialization.png")
        return qc.to_gate(label=" Reinitialization ") if to_gate else qc

    def mixing_operator(index_qubits: int, beta: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qc = QuantumCircuit(qr_index)

        qc.rx(2*beta, qr_index)

        qc.draw(output='mpl', filename="mixing_operator.png")
        return qc.to_gate(label=" Mixing Operator ") if to_gate else qc

    qr_index = QuantumRegister(index_qubits, "index") # index register
    qr_data = QuantumRegister(data_qubits, "data") # data register
    qr_f = QuantumRegister(1, "flag") # flag register
    cr_index = ClassicalRegister(index_qubits, "c_index") # classical register storing the measurement result of index register
    qc = QuantumCircuit(qr_index, qr_data, qr_f, cr_index)

    ### initialize the index register with uniform superposition state ###
    qc.h(qr_index)

    ### DO NOT CHANGE THE CODE BELOW
    p = 5
    alpha = 1
    for i in range(p):

        ### set fixed parameters for each round ###
        beta = 1 - (i + 1) / p
        gamma = (i + 1) / p

        ### return part ###
        qc.append(phase_return(index_qubits, gamma, L1, L2), qr_index)

        ### step 1: cost calculation ###
        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2), qr_index[:] + qr_data[:])

        ### step 2: Constraint testing ###
        qc.append(constraint_testing(data_qubits, C_max), qr_data[:] + qr_f[:])

        ### step 3: penalty dephasing ###
        qc.append(penalty_dephasing(data_qubits, alpha, gamma), qr_data[:] + qr_f[:])

        ### step 4: reinitialization ###
        qc.append(reinitialization(index_qubits, data_qubits, C1, C2, C_max), qr_index[:] + qr_data[:] + qr_f[:])

        ### mixing operator ###
        qc.append(mixing_operator(index_qubits, beta), qr_index)

    ### measure the index ###
    ### since the default measurement outcome is shown in big endian, it is necessary to reverse the classical bits in order to unify the endian ###
    qc.measure(qr_index, cr_index[::-1])

    return qc

def grade_ex4c():
    backend = Aer.get_backend('qasm_simulator')
    for ex in instance_examples:
        test_qc = solver_function(ex['L1'], ex['L2'], ex['C1'], ex['C2'], ex['C_max'])
        test_qc.draw(output='mpl', filename="4c.png")
        job = execute(test_qc, backend, shots=1000)
        result = job.result()
        print(result)
        count =result.get_counts()
        print(count)
        hist = plot_histogram(count)
        hist.savefig('4c-histogram.png')

grade_ex4c()

実行すると、以下のような感じで回答が表示される。

{'01101111110': 1, '11100011000': 1, '10100010000': 1, '10001111001': 1, '10101010100': 1, '10100101110': 1, '11110011100': 2, '00110110010': 1, '10010110010': 1, '00111011101': 1, '10111110100': 1, '01011110011': 1, '11101110110': 1, '11101110100': 1, '11110011000': 2, '10100110010': 1, '11111110011': 1, '00101110110': 2, '10010100000': 1, '10011111100': 1, '01110101100': 1, '00100101000': 1, '00011111010': 2, '10110011011': 2, '10110010101': 1, '01011111100': 1, '10101011110': 1, '00101110001': 1, '00101101000': 1, '11101101000': 1, '00110111011': 1, '00101011010': 2, '01111101100': 1, '01101101000': 1, '10011111110': 1, '10001111010': 1, '00110011010': 1, '11001111010': 1, '10110011100': 1, '11011111010': 1, '11111110010': 2, '11101100010': 1, '11111110110': 2, '10110100000': 1, '10110001001': 1, '01101111101': 1, '11010111010': 1, '11110110000': 2, '10101011001': 2, '10011110010': 1, '11001011010': 1, '11111010010': 1, '10100011000': 2, '11010111000': 1, '10101001001': 1, '11101111011': 1, '01110111110': 1, '00101101001': 1, '00110010000': 1, '10111100110': 1, '10101010001': 1, '10111011001': 1, '00100010100': 1, '11110110100': 1, '00111011001': 1, '10000111000': 1, '01110110110': 1, '01101011010': 1, '00101101100': 1, '01110110010': 1, '00101100010': 1, '00101110011': 3, '11101011100': 1, '11111110100': 1, '11101101001': 1, '11111110101': 3, '01001111100': 1, '01001111110': 1, '00111101000': 1, '10100011110': 1, '10111101010': 5, '10111110000': 12, '01111111110': 4, '00000111000': 1, '11111011010': 2, '01110011110': 1, '00110100000': 1, '11100111110': 1, '01111111010': 12, '11011011000': 1, '00011111100': 1, '11100110001': 1, '00110111000': 2, '11110111111': 1, '10111011100': 6, '11111111100': 7, '10111110111': 2, '01111111100': 6, '10100111000': 6, '10111011010': 6, '00000111010': 1, '10110101000': 1, '10110011010': 4, '00111101001': 2, '10101110110': 1, '11101111001': 1, '10100110000': 1, '10101110010': 2, '00010111000': 1, '00100111110': 2, '10111110101': 1, '10110111010': 12, '01101001010': 1, '00100111011': 1, '10001111000': 4, '11101101010': 1, '10011111000': 3, '11110111011': 2, '01110011000': 1, '00110111001': 4, '11111110000': 5, '01101111000': 12, '10111101000': 6, '10101101010': 1, '01111111000': 27, '01111111101': 1, '10100111110': 1, '10101011010': 6, '10111011000': 22, '11100111011': 1, '10101011100': 2, '10110101100': 1, '10110101001': 1, '11111111010': 17, '11100010010': 1, '00111110100': 2, '00001111010': 1, '10111111010': 28, '10101110000': 4, '00100010101': 1, '10110001000': 2, '10101101110': 1, '11001011100': 1, '00101111000': 11, '01111110000': 2, '11110101010': 1, '00110111010': 3, '11101100110': 1, '10001111110': 2, '01111001010': 2, '10110111111': 2, '10001011000': 2, '10100111001': 3, '11001111000': 4, '11111011000': 12, '10111111000': 70, '10111111100': 7, '11110110010': 3, '10100011011': 1, '10100111100': 7, '00111011000': 6, '01110111100': 4, '01000111010': 1, '10110111100': 9, '01111110010': 4, '10101111100': 10, '00110110000': 1, '11101111100': 6, '00111001100': 1, '10110111000': 30, '01111000010': 1, '11100111000': 9, '01101110010': 1, '00111110001': 1, '01010010010': 1, '11111010000': 2, '10110110100': 3, '00011111000': 1, '11111111000': 18, '10110110010': 1, '10100111010': 7, '10100011100': 1, '00101111010': 9, '00110101100': 1, '11101111000': 25, '10101111000': 42, '01101111100': 1, '00100110000': 1, '01011111000': 2, '10101101000': 3, '00101111100': 6, '11101111010': 6, '11101011000': 8, '10110110110': 1, '01101111010': 4, '01011111011': 1, '01111011010': 3, '11111101000': 5, '10111111110': 7, '11100011010': 1, '00010111011': 1, '01101110100': 1, '11111011100': 2, '00100111000': 2, '00101110100': 1, '01100111110': 1, '10000111100': 1, '00110111110': 2, '01000111110': 1, '10110111110': 6, '10010111110': 1, '10100011010': 1, '11100110101': 1, '00100111010': 1, '10010011000': 1, '00111111001': 2, '10111111011': 4, '10101110100': 4, '00110111100': 6, '00101111001': 2, '01110111000': 5, '00111101110': 1, '10101111010': 17, '00001011101': 1, '00111101100': 2, '00111111111': 1, '10110011000': 3, '00111110011': 1, '01111011000': 3, '01010111000': 1, '10101111110': 5, '10010111000': 1, '11110111001': 2, '01110111001': 1, '10111110001': 2, '01011110100': 1, '00100011000': 1, '11101011010': 2, '10111010010': 1, '11110011010': 6, '11110111010': 11, '11100111010': 5, '00111100110': 1, '11111111011': 1, '00011101001': 1, '01110110100': 1, '10101011101': 1, '00111001001': 1, '10101001010': 3, '00011110101': 1, '11011010110': 1, '10101011000': 10, '00111111100': 4, '10101111101': 1, '10011110100': 1, '00111111110': 6, '11010001010': 1, '00110101010': 2, '11001111111': 1, '11011111000': 1, '11010000011': 1, '10101101001': 1, '00111111000': 34, '01110111010': 4, '10111111001': 2, '10110010110': 1, '10100111011': 1, '00101011001': 2, '00111110000': 2, '11111100001': 1, '10001110100': 1, '00101011110': 1, '11111110111': 3, '01100111100': 2, '10001001010': 1, '01111100011': 1, '11110111100': 2, '00111101010': 1, '11110111000': 14, '11110111110': 6, '00010011000': 1, '00100111100': 1, '11101110010': 2, '00111011010': 3, '10111110010': 2, '11101110000': 3, '11111111110': 2, '00111011100': 1, '00111111010': 8, '01111011001': 2, '00101111110': 2}

ヒストグラムにすると、こんな感じになる。

"00111111000"と"10110111000"が最適解のはずで、このケースの値が多くなれば正解だ。

phase_return

では、関数毎に見ていこう。

    def phase_return(index_qubits: int, gamma: float, L1: list, L2: list, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qc = QuantumCircuit(qr_index)
        ### U_1(gamma * (lambda2 - lambda1)) for each qubit ###
        for i in range(len(L1)):
            qc.p(-gamma * (L2[i] - L1[i]) / 2, qr_index[i])

        qc.draw(output='mpl', filename="phase_return.png")
        return qc.to_gate(label=" Penalty Dephasing ") if to_gate else qc

収益計算パートで、L2からL1を引けばよい。
注意点としては、u1ゲートかpゲートを使うのだが、u1ゲートは非推奨だからpゲートを使うのがよいだろう。
初期コメントでは、U1に「gamma * (lambda2 - lambda1)」を設定すると教えてくれているが微妙に違う。

公式ページのjupyterの方に書いてある式に従って設定する必要がある。
回転ゲート U_1\left(\gamma_i \left(\lambda_{2}^{t}-\lambda_{1}^{t}\right)\right)= e^{-i \frac{\gamma_i \left(\lambda_{2}^{t}-\lambda_{1}^{t}\right)} 2}

回路図はこうなる。

qiskitの場合、量子回路を生成して実行する形となる。
そのため、回路を生成したタイミングでpゲートに指定する値は確定済みになっている。

qft

量子フーリエ変換の関数。
組み込みのqiskit.circuit.library.QFTがあるが自分で定義している。
組み込みのQFTだとエンディアンの都合が悪いからだ。
詳しくは、(qiskitのQFTモジュールを使う際のエンディアンに気を付ける)[https://qiita.com/ohken322/items/b1006db91861d4a54d22]を見るとよい。

    def qft(data_qubits: int, inverse=False, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qc = QuantumCircuit(data_qubits)
        N = data_qubits-1
        qc.h(N)
        for n in range(1, N+1):
            qc.cp(2*np.pi / 2**(n+1), N-n, N)

        for i in range(1, N):
            qc.h(N-i)
            for n in range(1, N-i+1):
                qc.cp(2*np.pi / 2**(n+1), N-(n+i), N-i)

        qc.h(0)

        # Calculate IQFT
        qc = qc.inverse() if inverse else qc

        qc.draw(output='mpl', filename="qft.png")
        return qc.to_gate(label=" QFT ") if to_gate and not inverse else qc.to_gate(label=" IQFT ") if to_gate else qc

回路図は、こうなります。

qft_add

QFT加算器。QFT Const Adderを見ると詳しく説明してくれている。

    def qft_add(data_qubits: int, const: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qc = QuantumCircuit(data_qubits)
        N = data_qubits-1
        bin_const = [int(x) for x in bin(const)[2:]] # Big endian
        bin_const = [0]*(N-len(bin_const)) + bin_const
        for n in range(1, N+1):
            if bin_const[n-1]:
                qc.p(2*np.pi / 2**(n+1), N)

        for i in range(1, N+1):
            for n in range(N-i+1):
                if bin_const[n+i-1]:
                    qc.p(2*np.pi / 2**(n+1), N-i)

        qc.draw(output='mpl', filename="qft_add.png")
        return qc.to_gate(label=f" [ +{const} ] ") if to_gate else qc

cost_calculation

コストを足していく。そのときにindexビットを見て、どちらが選ばれているか判断する。
qft_add(cost_qubits, c2).control(1)のcontrol(1)がそれ。

    def cost_calculation(index_qubits: int, cost_qubits: int, list1: list, list2: list, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qr_data = QuantumRegister(data_qubits, "data")
        qc = QuantumCircuit(qr_index, qr_data)

        qc.append(qft(cost_qubits), qr_data)
        for i, (c1, c2) in enumerate(zip(list1, list2)):
            qc.append(qft_add(cost_qubits, c2).control(1), [qr_index[i]] + qr_data[:])
            qc.x(qr_index[i])
            qc.append(qft_add(cost_qubits, c1).control(1), [qr_index[i]] + qr_data[:])
            qc.x(qr_index[i])
        qc.append(qft(cost_qubits, inverse=True), qr_data)

        qc.draw(output='mpl', filename="cost_calculation.png")
        return qc.to_gate(label=" Cost Calculation ") if to_gate else qc

回路図を見ると、制御がわかりやすい。

Xゲートを使って、打ち消しを入れている。
例えば、0番目のindexビットが0の場合は、最初のCNOTの先にある+4は適用されず、次にXゲートが適用されて、1になり、その次のCNOTの先にある+2が適用される。そして、またXゲートで元に戻る。
このゲートの挟み込みは、よく使われる回路パターンだ。

constraint_testing

    def constraint_testing(data_qubits: int, C_max: int, to_gate=True) -> Union[Gate, QuantumCircuit]:
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)

        c = data_qubits - 1
        w = 2**c - C_max - 1

        if w > 0:
            # Add w to cost(z)
            qc.append(qft(data_qubits), qr_data)
            qc.append(qft_add(data_qubits, w), qr_data)
            qc.append(qft(data_qubits, inverse=True), qr_data)
            # Set F if the condition cost(z) > C_max holds
            qc.cx(qr_data[c], qr_f)

        qc.draw(output='mpl', filename="constraint_testing.png")
        return qc.to_gate(label=" constraint_testing ") if to_gate else qc

劣化コストの値がC_{max}より大きいindexにペナルティを追加する。
コストを足していった結果、QRAMの上位ビットがキャリーオーバしたかどうかで判定する。
まず、そのために以下でQRAMのビット数を計算する。

c = data_qubits - 1
w = 2**c - C_max - 1

QRAMに足した結果、qc.cx(qr_data[c], qr_f)で上位ビットを見てコストオーバだとflagをたてる。

penalty_dephasing

ビット毎に見ていってペナルティを適用する。

\alpha\left(cost(z)-C_{m a x}\right)=\sum_{j=0}^{k-1} 2^{j} \alpha A_{1}[j]-2^{c} \alpha
    def penalty_dephasing(data_qubits: int, alpha: float, gamma: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_data, qr_f)
        for i in range(data_qubits):
            qc.cp((2**i) * alpha * gamma, qr_f, qr_data[i])
        
        # Add phase -alpha * (Cmax + w) = -alpha * 2^(cost_qubits-1) to all
        # infeasible states
        C_max_plus_w = 2**(data_qubits-1)-1
        qc.p(-C_max_plus_w * alpha * gamma, qr_f)

        qc.draw(output='mpl', filename="penalty_dephasing.png")
        return qc.to_gate(label=" Penalty Dephasing ") if to_gate else qc

reinitialization

再初期化は、constraint_testing()cost_calculation()だけ行う。
実行する順番に注意。

    def reinitialization(index_qubits: int, data_qubits: int, C1: list, C2: list, C_max: int, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qr_data = QuantumRegister(data_qubits, "data")
        qr_f = QuantumRegister(1, "flag")
        qc = QuantumCircuit(qr_index, qr_data, qr_f)

        qc.append(constraint_testing(data_qubits, C_max, False).inverse().to_gate(label="  Constraint Testing +  "), qr_data[:] + qr_f[:])
        qc.append(cost_calculation(index_qubits, data_qubits, C1, C2, False).inverse().to_gate(label="  Cost Calculation +  "), qr_index[:] + qr_data[:])

        qc.draw(output='mpl', filename="reinitialization.png")
        return qc.to_gate(label=" Reinitialization ") if to_gate else qc

mixing_operator

この演算子は、index registerの各qubitへのR_x(2\beta_i)ゲートによって実現できます。

の説明どおりになやればよい。

    def mixing_operator(index_qubits: int, beta: float, to_gate = True) -> Union[Gate, QuantumCircuit]:
        qr_index = QuantumRegister(index_qubits, "index")
        qc = QuantumCircuit(qr_index)

        qc.rx(2*beta, qr_index)

        qc.draw(output='mpl', filename="mixing_operator.png")
        return qc.to_gate(label=" Mixing Operator ") if to_gate else qc

全体の回路

以上!

Discussion