注釈

このページは docs/tutorials/qpe_with_sampler.ipynb から生成されました。

Sampler primitive を使用した量子位相推定

量子位相推定（QPE）アルゴリズムは、量子計算における重要なサブルーチンです。ショアの因数分解アルゴリズムや連立方程式の量子アルゴリズム(HHLアルゴリズム) など、多くの量子アルゴリズムの中心的な構成要素として機能しています。このチュートリアルでは、パラメーター化されたQPE回路を構築し、Sampler primitive を使用して実行します。

始める前に

このチュートリアルには、Qiskit Runtimeのサービスインスタンスが必要です。まだの方は、Qiskitの 入門ガイド の指示に従ってセットアップしてください。

背景情報

量子位相推定アルゴリズム

QPE アルゴリズム [1][2] は、ユニタリー演算子 \(U\) が固有ベクトル \(|\psi \rangle\) 、固有値 \(e^{2\pi i \theta}\) を持つとき、thetaの値を推定します。推定のために、状態 \(|\psi \rangle\) を準備し、制御- \(U^{2^j}\) 演算を実行する、ブラックボックス ( オラクル と呼ばれることもあります) を仮定します。

ブラックボックスに依存しているため、QPEアルゴリズムは実際には完全なアルゴリズムではなく、サブルーチンとして考えることができます。このサブルーチンは他のサブルーチンと連携して、ショアの因数分解アルゴリズムやHHLアルゴリズムなど他の計算を実行することができます。

ここではQPEアルゴリズムの詳細を説明はしません。より詳しく知りたい場合は、Qiskitテキストブック [3]のQPEアルゴリズムについての章を読んでみてください。

概要

上で説明したように、QPE アルゴリズムには状態 \(|\psi \rangle\) を準備し、制御- \(U^{2^j}\) 演算を実行するためのブラックボックスがあります。 このチュートリアルでは、さまざまな位相に対応するさまざまなブラックボックスを含む一連の QPE 回路を準備します。 次に、Sampler primitive を使用してこれらの回路を実行し、結果を分析します。 ご覧のとおり、Sampler primitiveを使用すると、パラメーター化された回路を非常に簡単に実行できます。 一連の QPE 回路を作成する代わりに、ブラックボックスが生成する位相を指定するパラメーターとパラメーターの一連の位相値を使用して 1 つ の QPE 回路を作成するだけで済みます。

注記

ショアの素因数分解アルゴリズムのようなQPEアルゴリズムの典型的な使用法では、ブラックボックスはユーザーによって設定されず、他のサブルーチンによって指定されます。しかし、このチュートリアルの目的は、QPEアルゴリズムを使用して、Sampler primitiveを使用したパラメーター化回路の実行の容易さを説明することです。

ステップ1：QPE回路の作成

パラメーター化されたQPE回路を作成する

QPE アルゴリズムの手順は、以下のとおりです。

1. 2 つの量子レジスター (1 つ目は位相を推定するため、2 つ目は固有ベクトル \(|\psi\rangle\) を格納するため) と 1 つの古典レジスターを読み出し用に持つ回路を作成します。

2. 最初のレジスターの量子ビットを \(|0\rangle\) として、2 番目のレジスターを \(|\psi\rangle\) として初期化します。

3. 最初のレジスターに重ね合わせを作成します。

4. 制御-\(U^{2^j}\) ブラックボックスを適用します。

5. 逆量子フーリエ変換 (QFT) を最初のレジスターに適用します。

6. 最初のレジスターを測定します。

次のコードは、キーワード引数 theta と num_qubits を与えて、QPE回路を作成するための関数 create_qpe_circuit を作成します。 theta は \(2\pi\) ファクターを含んでいないことに注意してください。num_qubits 引数は最初のレジスターの量子ビット数を指定します。子ビット数が多いほど、位相推定の精度が良くなります。

import numpy as np

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from qiskit.circuit.library import QFT


def create_qpe_circuit(theta, num_qubits):
    """Creates a QPE circuit given theta and num_qubits."""

    # Step 1: Create a circuit with two quantum registers and one classical register.
    first = QuantumRegister(
        size=num_qubits, name="first"
    )  # the first register for phase estimation
    second = QuantumRegister(
        size=1, name="second"
    )  # the second register for storing eigenvector |psi>
    classical = ClassicalRegister(
        size=num_qubits, name="readout"
    )  # classical register for readout
    qpe_circuit = QuantumCircuit(first, second, classical)

    # Step 2: Initialize the qubits.
    # All qubits are initialized in |0> by default, no extra code is needed to initialize the first register.
    qpe_circuit.x(
        second
    )  # Initialize the second register with state |psi>, which is |1> in this example.

    # Step 3: Create superposition in the first register.
    qpe_circuit.barrier()  # Add barriers to separate each step of the algorithm for better visualization.
    qpe_circuit.h(first)

    # Step 4: Apply a controlled-U^(2^j) black box.
    qpe_circuit.barrier()
    for j in range(num_qubits):
        qpe_circuit.cp(
            theta * 2 * np.pi * (2**j), j, num_qubits
        )  # Theta doesn't contain the 2 pi factor.

    # Step 5: Apply an inverse QFT to the first register.
    qpe_circuit.barrier()
    qpe_circuit.compose(QFT(num_qubits, inverse=True), inplace=True)

    # Step 6: Measure the first register.
    qpe_circuit.barrier()
    qpe_circuit.measure(first, classical)

    return qpe_circuit

create_qpe_circuit 関数の theta 引数に実数を渡すと、位相が固定された回路がエンコードされます。

num_qubits = 4
qpe_circuit_fixed_phase = create_qpe_circuit(
    1 / 2, num_qubits
)  # Create a QPE circuit with fixed theta=1/2.
qpe_circuit_fixed_phase.draw("mpl")

しかし、代わりに Parameter <https://qiskit.org/documentation/stubs/qiskit.circuit.Parameter.html>`__ オブジェクトを渡すと、回路が作成された後にパラメーター値を割り当てることができる、パラメーター化された回路が返されます。このチュートリアルの残りの部分では、パラメーター化されたバージョンのQPE回路を使用することになります。

from qiskit.circuit import Parameter

theta = Parameter(
    "theta"
)  # Create a parameter `theta` whose values can be assigned later.
qpe_circuit_parameterized = create_qpe_circuit(theta, num_qubits)
qpe_circuit_parameterized.draw("mpl")

後で割り当てる位相の値のリストを作成する

パラメーター化されたQPE回路を作成した後、次のステップで回路に割り当てる位相の値のリストを作成します。以下のコードを使用して、 0 から 2 までの21個の位相の値のリストを等間隔で作成します： 0, 0.1, 0.2, …, 1.9, 2.0

number_of_phases = 21
phases = np.linspace(0, 2, number_of_phases)
individual_phases = [
    [ph] for ph in phases
]  # Phases need to be expressed as a list of lists.

ステップ2: 回路をクラウド上の量子コンピューターに送信する

Qiskit Runtime サービスに接続する

まず、 最初のステップ で作成したQiskit Runtimeのサービスインスタンスに接続します。ここでは ibmq_qasm_simulator を使用して、このチュートリアルを実行することにします。

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = "ibmq_qasm_simulator"  # use the simulator

Sampler primitive を使用してパラメーター化された回路を実行する

Qiskit Runtimeサービスが接続され、パラメーター化されたQPE回路、およびパラメーターの値のリストがあれば、Sampler primitiveを使用してQPE回路を実行する準備が整ったことになります。 with ... 構文は、Qiskit Runtime sessionを開きます。Session内では、たった1行のコードで21個の異なるパラメーターの値を持つパラメーター化されたQPE回路を実行することになります。Sampler primitiveは、パラメーター値を受け取り、それを回路に割り当て、異なるパラメーター値を持つ21個の回路として実行し、ビット文字列の (擬似) 確率を含む1つのジョブ結果を返します。これにより、何十行ものコードを書く手間が省けます。

Estimator primitiveは、パラメーター化された回路に対して同様のAPIインターフェイスを持っています。詳しくは APIリファレンス を参照してください。

from qiskit_ibm_runtime import Sampler, Session

with Session(service=service, backend=backend):
    results = (
        Sampler()
        .run(
            [qpe_circuit_parameterized] * len(individual_phases),
            parameter_values=individual_phases,
        )
        .result()
    )

ステップ 3: 結果の分析

1つの回路の結果を分析する

ジョブが完了したら、回路のヒストグラムを見ることで結果の解析を開始することができます。測定結果の（準）確率分布は results.quasi_dists に格納されており、興味のある回路のインデックス（以下の例では \(idx=6\) ）を渡すことで個々の回路の結果にアクセスすることが可能です。

from qiskit.tools.visualization import plot_histogram

idx = 6
plot_histogram(
    results.quasi_dists[idx].binary_probabilities(),
    legend=[f"$\\theta$={phases[idx]:.3f}"],
)

\(\theta\) を推定するために、 \(N_1\) (2進数で 1010 、10進数で 10) を \(2^n\) で割る必要があります。ここで \(n\) は最初のレジスターの量子ビットの数（この例では \(n=4\) ）です。

\[\theta = \frac{N_1}{2^n} = \frac{10}{2^4} = 0.625.\]

これは \(0.6\) の期待値に近いですが、最も可能性の高い結果 (2進数で 1010 、10進数で 10) と2番目に可能性の高い結果 (2進数で 1001 、10進数で 9) を加重平均すれば、もっと良い推定値が得られます。

\[\theta = \frac{P_1 \times \frac{N_1}{2^n} + P_2 \times \frac{N_2}{2^n}}{P_1 + P_2} = \frac{0.554 \times \frac{10}{2^4} + 0.269 \times \frac{9}{2^4}}{0.554 + 0.269} = 0.605,\]

ここで、 \(N_1\) と \(P_1\) は最も可能性の高い結果とその確率、 \(N_2\) と \(P_2\) は2番目に可能性の高い結果とその確率を表します。その結果、 \(0.605\) となり、期待値である \(0.6\) にかなり近くなりました。この方法を用いて、すべての回路の結果を分析することにします。

すべての回路の結果を分析する

以下のヘルパー関数を利用することで、ある位相における \(N_1\), \(P_1\), \(N_2\), \(P_2\) の値を求め、その推定値 \(\theta\) を計算することができます。理想的には \(N_2\) は \(N_1\) の「近傍」であるべきですが（例えば、 1010 の近傍は 1001 と 1011 ）、実際の量子系から得られた結果では、ノイズがあるため 2番目に可能性が高い結果は \(N_1\) の近傍ではない可能性があります。ヘルパー関数は \(N_2\) の値を \(N_1\) の近傍領域から取得します。

def most_likely_bitstring(results_dict):
    """Finds the most likely outcome bit string from a result dictionary."""
    return max(results_dict, key=results_dict.get)


def find_neighbors(bitstring):
    """Finds the neighbors of a bit string.

    Example:
        For bit string '1010', this function returns ('1001', '1011')
    """
    if bitstring == len(bitstring) * "0":
        neighbor_left = len(bitstring) * "1"
    else:
        neighbor_left = format((int(bitstring, 2) - 1), "0%sb" % len(bitstring))

    if bitstring == len(bitstring) * "1":
        neighbor_right = len(bitstring) * "0"
    else:
        neighbor_right = format((int(bitstring, 2) + 1), "0%sb" % len(bitstring))

    return (neighbor_left, neighbor_right)


def estimate_phase(results_dict):
    """Estimates the phase from a result dictionary of a QPE circuit."""

    # Find the most likely outcome bit string N1 and its neighbors.
    num_1_key = most_likely_bitstring(results_dict)
    neighbor_left, neighbor_right = find_neighbors(num_1_key)

    # Get probabilities of N1 and its neighbors.
    num_1_prob = results_dict.get(num_1_key)
    neighbor_left_prob = results_dict.get(neighbor_left)
    neighbor_right_prob = results_dict.get(neighbor_right)

    # Find the second most likely outcome N2 and its probability P2 among the neighbors.
    if neighbor_left_prob is None:
        # neighbor_left doesn't exist
        if neighbor_right_prob is None:
            # both neighbors don't exist, N2 is N1
            num_2_key = num_1_key
            num_2_prob = num_1_prob
        else:
            # If only neighbor_left doesn't exist, N2 is neighbor_right.
            num_2_key = neighbor_right
            num_2_prob = neighbor_right_prob
    elif neighbor_right_prob is None:
        # If only neighbor_right doesn't exist, N2 is neighbor_left.
        num_2_key = neighbor_left
        num_2_prob = neighbor_left_prob
    elif neighbor_left_prob > neighbor_right_prob:
        # Both neighbors exist and neighbor_left has higher probability, so N2 is neighbor_left.
        num_2_key = neighbor_left
        num_2_prob = neighbor_left_prob
    else:
        # Both neighbors exist and neighbor_right has higher probability, so N2 is neighor_right.
        num_2_key = neighbor_right
        num_2_prob = neighbor_right_prob

    # Calculate the estimated phases for N1 and N2.
    num_qubits = len(num_1_key)
    num_1_phase = int(num_1_key, 2) / 2**num_qubits
    num_2_phase = int(num_2_key, 2) / 2**num_qubits

    # Calculate the weighted average phase from N1 and N2.
    phase_estimated = (num_1_phase * num_1_prob + num_2_phase * num_2_prob) / (
        num_1_prob + num_2_prob
    )

    return phase_estimated

estimate_phase ヘルパー関数を使用すると、すべての回路の結果から位相を推定することができます。

qpe_solutions = []
for idx, result_dict in enumerate(results.quasi_dists):
    qpe_solutions.append(estimate_phase(result_dict.binary_probabilities()))

位相 \(\theta\) の理想的な解は、固有値 \(e^{2 \pi i \theta}\) が \(2 \pi\) で周期的なので、周期 1 を持つ周期的な解となります。以下のセルを実行して理想的な解を生成し、QPEアルゴリズムから得られる解と比較することができます。

ideal_solutions = np.append(
    phases[: (number_of_phases - 1) // 2],  # first period
    np.subtract(phases[(number_of_phases - 1) // 2 : -1], 1),  # second period
)
ideal_solutions = np.append(
    ideal_solutions, np.subtract(phases[-1], 2)
)  # starting point of the third period

以下のコードを実行して、解を視覚化できます。

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(10, 6))
plt.plot(phases, ideal_solutions, "--", label="Ideal solutions")
plt.plot(phases, qpe_solutions, "o", label="QPE solutions")

plt.title("Quantum Phase Estimation Algorithm")
plt.xlabel("Input Phase")
plt.ylabel("Output Phase")
plt.legend()

<matplotlib.legend.Legend at 0x130347bb0>

ご覧のように、QPEアルゴリズムで得られた解は、理想的な解に忠実に沿っています。おめでとうございます。QPEアルゴリズムを正常に実行し、良い解を得ることができました!

概要

このチュートリアルでは、パラメーター化されたQPE回路を作成し、Sampler primitiveを使って実行し、結果を解析して、理想的な解に近い解を得ました。Sampler primitiveを使うと、パラメーター化された回路を非常に簡単に実行できることがおわかりいただけると思います。

参考文献

1. Kitaev, A. Y. (1995). Quantum measurements and the Abelian stabilizer problem. arXiv preprint quant-ph/9511026.

2. Nielsen, M., & Chuang, I. (2010). Quantum computation and quantum information (2nd ed., pp. 221-226). Cambridge University Press.

3. Abbas, A., Andersson, S., Asfaw, A., Córcoles, A., Bello, L., Ben-Haim, Y., … & Wootton, J. (2020). Learn quantum computation using qiskit. URL: https://qiskit.org/textbook/preface.html (accessed 07/14/2022).

Qiskit のバージョンと著作権

import qiskit_ibm_runtime

qiskit_ibm_runtime.version.get_version_info()

'0.7.0'

import qiskit.tools.jupyter

%qiskit_version_table
%qiskit_copyright

Version Information

Qiskit Software	Version
qiskit-terra	0.21.0
qiskit-aer	0.10.4
qiskit-ibmq-provider	0.19.2
qiskit	0.37.0
System information
Python version	3.9.12
Python compiler	Clang 12.0.0
Python build	main, Apr 5 2022 01:53:17
OS	Darwin
CPUs	8
Memory (Gb)	16.0
Tue Jul 12 11:40:25 2022 CEST

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.