サンプリング推定

量子コンピュータの演算子の期待値を推定するためには、サンプリング測定が必要です。サンプリング測定では、量子回路の実行とそれに続く量子ビットの測定を複数回繰り返します。この繰り返し測定の統計量を用いて演算子の期待値を推定します。

実際の量子コンピュータでサンプリング推定を行う前に、量子回路シミュレータでどのようにシミュレーションできるかを見てみましょう。このチュートリアルでは、以下の方法を紹介します。

• Operatorのグループピング
• グルーピング結果から測定回路を構築
• 測定回路にショットを割り当て、サンプリング
• サンプリング結果から期待値を再構築

前提条件

このチュートリアルで使用するQURI Partsモジュール: quri-parts-circuit,quri-parts-core,quri-parts-qulacs

以下のようにインストールすることができます:

!pip install "quri-parts[qulacs]"

概要

このチュートリアルの主な目的は、sampling estimatorを用いてサンプリング推定を行う方法を紹介することです。sampling estimatorとは、サンプリング実験によって期待値を計算するQuantumEstimatorのことです。準備として、期待値を推定したい演算子と状態を作成します。

from quri_parts.core.operator import Operator, pauli_label, PAULI_IDENTITY
from quri_parts.core.state import quantum_state

op = Operator({
    pauli_label("Z0"): 0.25,
    pauli_label("Z1 Z2"): 2.0,
    pauli_label("X1 X2"): 0.5 + 0.25j,
    pauli_label("Z1 Y3"): 1.0j,
    pauli_label("Z2 Y3"): 1.5 + 0.5j,
    pauli_label("X1 Y3"): 2.0j,
    PAULI_IDENTITY: 3.0,
})

print("Operator:")
print(op)


state = quantum_state(4, bits=0b0101)
print("")
print("State:")
print(state)

#output
    Operator:
    0.25*Z0 + 2.0*Z1 Z2 + (0.5+0.25j)*X1 X2 + 1j*Z1 Y3 + (1.5+0.5j)*Z2 Y3 + 2j*X1 Y3 + 3.0*I
    
    State:
    ComputationalBasisState(qubit_count=4, bits=0b101, phase=0π/2)

ここでは，sampling estimatorを作成し，使用するためのサンプルコードを提供します。以下の仕様でsampling estimatorを作成します:

• 総ショット数: 5000
• 測定ファクトリ: ビットごとに可換なパウリへのグルーピング
• ショットアロケータ: 比例ショット配分

「測定ファクトリ」と「ショットアロケータ」については後のセクションで説明します。とりあえず、一時的にこれらを当然と考えて、それらを使ってsampling estimatorを作ってみましょう:

from quri_parts.core.estimator.sampling import create_sampling_estimator
from quri_parts.core.measurement import bitwise_commuting_pauli_measurement
from quri_parts.core.sampling.shots_allocator import create_proportional_shots_allocator
from quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_sampler

total_shots = 5000
concurrent_sampler = create_qulacs_vector_concurrent_sampler()
measurement_factory = bitwise_commuting_pauli_measurement
shots_allocator = create_proportional_shots_allocator()

sampling_estimator = create_sampling_estimator(
    total_shots,
    concurrent_sampler,
    measurement_factory,
    shots_allocator
)

サンプリングの結果を見てみましょう:

estimate = sampling_estimator(op, state)
print("Estimated expectation value:", estimate.value)
print("Standard error of this sampling estimation:", estimate.error)

#output
    Estimated expectation value: (0.6527269558463302+0.02582389806961259j)
    Standard error of this sampling estimation: 0.070823557319059

サンプリング推定の結果を、vector estimatorによる正確な結果と比較することができます。

from quri_parts.qulacs.estimator import create_qulacs_vector_estimator

vector_estimator = create_qulacs_vector_estimator()
vector_estimator(op, state)

#output
    _Estimate(value=(0.75+0j), error=0.0)

演算子を可換なパウリ文字列にグループ化

サンプリング実験では、演算子の期待値を見積もる方法の一つとして、各パウリ文字列の期待値を見積もり、それらを合計する方法があります。しかし、複数のパウリ文字列が可換であれば、一度に測定した方が効率的です。従って、最初のステップは、パウリ文字列をいくつかの可換パウリ文字列のセットにグループ化することです。このパウリ文字列のグループ化は、演算子推定の文脈で重要な研究課題です。

次に、グループ化の結果から、測定とサンプリングのための各グループの回路を構成する必要があり、これを測定回路と呼びます。測定回路に含まれる明示的なゲートは、グループ内のパウリ文字列に依存します。

ここでは、演算子と上記で構築した状態を使用します:

print("Operator:")
for pl, coeff in op.items():
    print(f"   + {coeff} * {pl}")

print("")
print("State:", state)

#output
    Operator:
    + 0.25 * Z0
    + 2.0 * Z1 Z2
    + (0.5+0.25j) * X1 X2
    + 1j * Z1 Y3
    + (1.5+0.5j) * Z2 Y3
    + 2j * X1 Y3
    + 3.0 * I
    
    State: ComputationalBasisState(qubit_count=4, bits=0b101, phase=0π/2)

QURI Partsでどのように

• 演算子のグループ化
• 測定回路の構築
• 測定結果によるパウリ文字列の再構築

を行うかを説明します。最後に、便利なオブジェクトを紹介して、このセクションを締めくくります:CommutablePauliSetMeasurementという便利なオブジェクトを紹介します。

パウリのグループ化

QURI Partsでは、グループ化関数はPauliGroupingで表されます。これらは、OperatorまたはPauliLabelのセットをCommutablePauliSetのセットに変換する関数です。関数のシグネチャは次式で与えられます:

from typing import Callable, Set, Union, Iterable
from typing_extensions import TypeAlias
from quri_parts.core.operator import PauliLabel

CommutablePauliSet: TypeAlias = Set[PauliLabel]

#: Represents a grouping function
PauliGrouping = Callable[
    [Union[Operator, Iterable[PauliLabel]]], Iterable[CommutablePauliSet]
]

例として、最も単純なパウリのグループ化の1つであるビット単位のグループ化を使用する。この場合、グループはパウリ文字列のビット単位の可換性に基づいて決定されます。このグループ分けは次のようにテストできます:

from quri_parts.core.operator.grouping import bitwise_pauli_grouping
pauli_sets = bitwise_pauli_grouping(op)
print(pauli_sets)

#output
    frozenset({frozenset({PauliLabel({(1, <SinglePauli.X: 1>), (3, <SinglePauli.Y: 2>)})}), frozenset({PauliLabel()}), frozenset({PauliLabel({(1, <SinglePauli.X: 1>), (2, <SinglePauli.X: 1>)})}), frozenset({PauliLabel({(0, <SinglePauli.Z: 3>)}), PauliLabel({(2, <SinglePauli.Z: 3>), (1, <SinglePauli.Z: 3>)})}), frozenset({PauliLabel({(2, <SinglePauli.Z: 3>), (3, <SinglePauli.Y: 2>)}), PauliLabel({(3, <SinglePauli.Y: 2>), (1, <SinglePauli.Z: 3>)})})})

グループ化メソッドはパウリラベルのセットのfrozensetを返すので、上記は少し複雑に見えます。

print(f"Number of groups: {len(pauli_sets)}")
for i, pauli_set in enumerate(pauli_sets):
    labels = ", ".join([str(pauli) for pauli in pauli_set])
    print(f"Group {i} contains: {labels}")

#output
    Number of groups: 5
    Group 0 contains: X1 Y3
    Group 1 contains: I
    Group 2 contains: X1 X2
    Group 3 contains: Z0, Z1 Z2
    Group 4 contains: Z2 Y3, Z1 Y3

ここで、QURI Partsで使用可能なグループ化関数の一覧を示します:

グループ化関数	説明
individual_pauli_grouping	グループ化なし
bitwise_pauli_grouping	ビットごとの可換性に基づくパウリ文字列の可換性のチェック
sorted_injection_grouping	arXiv:1908.06942で説明

測定回路

可換パウリセットの測定を行うには、測定前に適用する回路を構築する必要があります。ビット単位のパウリグループ化の場合、対応する回路はbitwise_commuting_pauli_measurement_circuit関数で構築することができます:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurement_circuit
from quri_parts.circuit import QuantumCircuit
from quri_parts.circuit.utils.circuit_drawer import draw_circuit

pauli_set = {pauli_label("X0 Z2 Y3"), pauli_label("X0 Z1 Y3")}
measurement_circuit = bitwise_commuting_pauli_measurement_circuit(pauli_set)
draw_circuit(QuantumCircuit(4, gates=measurement_circuit))

#output
         ___          
        | H |         
    0 --|0  |---------
        |___|         
                      
                      
    1 ----------------
                      
                      
                      
    2 ----------------
                      
         ___     ___  
        |Sdg|   | H | 
    3 --|1  |---|2  |-
        |___|   |___| 

状態のサンプリング

量子状態の可換パウリセットを測定するためには、以下のステップが必要です:

• 状態を準備する回路を構成
• 状態を準備した回路の後にパウリセットの測定回路を連結
• 結合した回路でサンプリング

ここでは簡単のために初期状態としてComputationalBasisStateを使用していますが、CircuitQuantumStateであれば何でも使用できます。

from quri_parts.core.state import quantum_state
from quri_parts.qulacs.sampler import create_qulacs_vector_sampler

sampler = create_qulacs_vector_sampler()

initial_state = quantum_state(4, bits=0b0101)

# Circuit for state preparation
state_prep_circuit = initial_state.circuit

# Measurement circuit
pauli_set = {pauli_label("Z2 Y3"), pauli_label("Z1 Y3")}
measurement_circuit = bitwise_commuting_pauli_measurement_circuit(pauli_set)

# Concatenate measurement circuit
sampled_circuit = state_prep_circuit + measurement_circuit
print("State preparation circuit concatenated with measurement circuit:")
draw_circuit(sampled_circuit)

# Sampling
sampling_result = sampler(sampled_circuit, shots=1000)
print("\n")
print("Sampling result:")
print({bin(bits): count for bits, count in sampling_result.items()})

#output
    State preparation circuit concatenated with measurement circuit:
         ___          
        | X |         
    0 --|0  |---------
        |___|         
                      
                      
    1 ----------------
                      
         ___          
        | X |         
    2 --|1  |---------
        |___|         
         ___     ___  
        |Sdg|   | H | 
    3 --|2  |---|3  |-
        |___|   |___| 


    Sampling result:
    {'0b101': 501, '0b1101': 499}

サンプリング結果からパウリ文字列の期待値を再構成

次に、測定されたパウリ文字列の固有値を測定結果から読み出す必要があります。QURI Partsでは、これをPauliReconstructorFactoryとPauliReconstructorで行います。PauliReconstructorは測定結果に対応するPauli文字列の固有値を返す関数で、PauliReconstructorFactoryはPauliLabelからPauliReconstructorを生成する関数です。QURI Partsで定義されているこれらの関数のシグネチャを以下に示します。

#: PauliReconstructor represents a function that reconstructs a value of a Pauli
#: operator from a measurement result of its measurement circuit.
PauliReconstructor: TypeAlias = Callable[[int], int]

#: PauliReconstructorFactory represents a factory function that returns
#: a :class:`~PauliReconstructor` for a given Pauli operator.
PauliReconstructorFactory: TypeAlias = Callable[[PauliLabel], PauliReconstructor]

上のセクションの例では、$Z_2 Y_3$と$Z_1 Y_3$でサンプリング測定を行い、0b1101と0b0101の$2$種類のビットパターンを得ました。ビットごとのパウリグループ化の場合、パウリ演算子の値はサンプリングされたビットパターンから以下のように再構成できます:

from quri_parts.core.measurement import bitwise_pauli_reconstructor_factory

# Obtain a reconstructor for Z2 Y3
reconstructor = bitwise_pauli_reconstructor_factory(pauli_label("Z2 Y3"))

# Reconstruct a value of Z2 Y3 from a sampled bit pattern 0b1101
pauli_value = reconstructor(0b1101)
print("Measurement result: 0b1101.", f"Corresponding eigenvalue: {pauli_value}")

# Reconstruct from 0b0101
pauli_value = reconstructor(0b0101)
print("Measurement result: 0b0101.", f"Corresponding eigenvalue: {pauli_value}")

#output
    Measurement result: 0b1101. Corresponding eigenvalue: 1
    Measurement result: 0b0101. Corresponding eigenvalue: -1

$Z_2 Y_3$の期待値を以下のように計算することができます:

pauli_exp = (
    reconstructor(0b1101) * sampling_result[0b1101] +
    reconstructor(0b0101) * sampling_result[0b0101]
) / 1000
print(pauli_exp)

# Equivalent to above
pauli_exp = sum(
    reconstructor(bits) * count for bits, count in sampling_result.items()
) / sum(sampling_result.values())
print(pauli_exp)

#output
    -0.002
    -0.002

QURI Partsは、サンプリング結果からPauliLabelの期待値を求める便利な方法として、trivial_pauli_expectation_estimator関数も提供しています。

from quri_parts.core.estimator.sampling import trivial_pauli_expectation_estimator

pauli_exp = trivial_pauli_expectation_estimator(sampling_result, pauli_label("Z2 Y3"))
print(pauli_exp)

#output
    -0.002

ここではビットごとのパウリグループ化を使っているので、trivial_pauli_expectation_estimatorが使えます。より一般的なケースでは、general_pauli_expectation_estimatorをPauliReconstructorFactoryと一緒に使う必要があります:

from quri_parts.core.estimator.sampling import general_pauli_expectation_estimator
pauli_exp = general_pauli_expectation_estimator(
    sampling_result, pauli_label("Z2 Y3"), bitwise_pauli_reconstructor_factory
)
print(pauli_exp)

#output
    -0.002

最後に、各パウリ項の寄与を合計することにより、元の演算子の期待値を見積もることができます。以下のように$Z_2 Y_3$の寄与を計算します:

# Coefficient of Z2 Y3 in op
coef = op[pauli_label("Z2 Y3")]
pauli_contrib = coef * pauli_exp
print(pauli_contrib)

#output
    (-0.003-0.001j)

すべてのパウリ文字列についてこの手順を繰り返すと、元の演算子の期待値を推定することができます。

CommutablePauliSetMeasurementオブジェクト

上記の手順は少し複雑なので、ショートカットの方法が用意されています。それを使うために、まずCommutablePauliSetMeasurementについて紹介します:これは以下の要素を含むデータ構造です:

• pauli_set: 一緒に測定される可換パウリ演算子のセット
• measurement_circuit: 与えられたpauli_setを測定するために必要な回路
• pauli_reconstructor_factory: サンプリング結果からパウリ演算子の値を再構成する関数のファクトリ

CommutablePauliSetMeasurementのセットは、特定の測定スキームに応じて構築する必要があります。それらは通常、CommutablePauliSetMeasurementFactoryによって取得されます。CommutablePauliSetMeasurementFactoryは、OperatorまたはPauliLabelsのシーケンスを取り込んで、CommutablePauliSetMeasurementのシーケンスを返します。

例として、ビットごとのパウリグループ化は、以下のようにできます:

from quri_parts.core.measurement import bitwise_commuting_pauli_measurement

measurements = bitwise_commuting_pauli_measurement(op)
print(f"Number of CommutablePauliSetMeasurement: {len(measurements)}")
measurement = measurements[-1]

#output
    Number of CommutablePauliSetMeasurement: 5

• 可換パウリセット

print(", ".join(map(lambda s: str(s), measurement.pauli_set)))

#output
    Z2 Y3, Z1 Y3

• 状態を準備した回路＋測定回路

  $\langle Z_2 Y_3 \rangle$と$\langle Z_1 Y_3 \rangle$を推定する回路を示し、$1000$ショットサンプリングを行う。

full_circuit = state.circuit + measurement.measurement_circuit
draw_circuit(full_circuit)

samling_result = sampler(full_circuit, 1000)
print("\n")
print(f"Measurement Result: {samling_result}")

#output
         ___          
        | X |         
    0 --|0  |---------
        |___|         
                      
                      
    1 ----------------
                      
         ___          
        | X |         
    2 --|1  |---------
        |___|         
         ___     ___  
        |Sdg|   | H | 
    3 --|2  |---|3  |-
        |___|   |___| 
    
    
    Measurement Result: Counter({13: 501, 5: 499})

• 再構築ファクトリ

  measurementからPauliReconstructorFactoryを取り出し、$\langle Z_2 Y_3 \rangle$と$\langle Z_1 Y_3 \rangle$を計算します。

print(f"Reconstructor factory: {measurement.pauli_reconstructor_factory}")

for pl in measurement.pauli_set:
    estimate = general_pauli_expectation_estimator(
        sampling_result, pl, measurement.pauli_reconstructor_factory
    )
    print(f" <{pl}> = {estimate}")
    

#output
    Reconstructor factory: <function bitwise_pauli_reconstructor_factory at 0x127af98b0>
     <Z2 Y3> = -0.002
     <Z1 Y3> = 0.002

ここで、QURI Partsで提供されているPauliReconstructorFactoryを列挙します:

• bitwise_commuting_pauli_measurement: ビットごとの可換性にもとづいてグループ化
• individual_pauli_measurement: グループ化なし(全てのPauliLabelが別のグループとなる)

ショットアロケータ

ここまでで、OperatorをいくつかのCommutablePauliSetにグループ化する方法を学びました。また、測定回路からサンプリング結果を読み出して、個々のPauliLabelの期待値を評価する方法も学びました。最後に推定に必要なのがPauliSamplingShotsAllocatorです。これは、各パウリセットの測定にどのようにサンプリングショットを割り当てるかを指定します。

インターフェース

QURI Partsでは、サンプリング推定に使用されるショットアロケータは、PauliSamplingShotsAllocator ^[1]で表されます。これは、対象となるOperatorまたはCommutablePauliSetのセットに従って、総ショット数をCommutablePauliSetに分配します。関数のシグネチャは以下で与えられます。

from typing import Collection
from quri_parts.core.sampling import PauliSamplingSetting

CommutablePauliSetMeasurement: TypeAlias = Callable[
    [Operator, Collection[CommutablePauliSet], int], Collection[PauliSamplingSetting]
]

ここで、返されるPauliSamplingSettingは、CommutablePauliSetとそれに割り当てられたショット数を保持する名前付きタプルです。

[1] ^ QURI Partsが提供する別のタイプのショットアロケータWeightedSamplingShotsAllocatorがあります。しかし、これらはサンプリング推定には使われません。そのため、このチュートリアルでは紹介しません。

使用可能なショットアロケータ

最後に、QURI Partsで使用可能なショットアロケータを紹介します:

from quri_parts.core.sampling.shots_allocator import (
    create_equipartition_shots_allocator,
    create_proportional_shots_allocator,
    create_weighted_random_shots_allocator,
)
# Allocates shots equally among the Pauli sets
allocator = create_equipartition_shots_allocator()
# Allocates shots proportional to Pauli coefficients in the operator
allocator = create_proportional_shots_allocator()
# Allocates shots using random weights
allocator = create_weighted_random_shots_allocator(seed=777)

以下にアロケータをまとめます:

アロケータ	生成関数	説明	リファレンス
equipartition shots allocator	create_equipartition_shots_allocator	パウリセットに均等にショットを割り当てます	-
proportional shots allocator	create_proportional_shots_allocator	演算子のパウリ係数に比例したショットを割り当てます	- arXiv:2112.07416 (2021). 付録 D. - Phys. Rev. A 92, 042303
weighted random shots allocator	create_weighted_random_shots_allocator	ランダムな重みを使用してショットを割り当てます	arXiv:2004.06252 (2020)

サンプリング推定

サンプリング推定を行うために必要なすべてのコンポーネントを紹介するのは長い道のりでした。QURI Partsでは、サンプリング推定計算をカスタマイズできるsampling_estimate関数を提供しています。CircuitQuantumStateの演算子の期待値を推定するために必要な要素は以下の通りです:

• 総ショット: 期待値の推定に使用するショットの総数
• Concurrent Sampler: 全ての (状態準備の回路＋測定回路) のサンプリングを行うConcurrentSampler
• 測定ファクトリ: 入力演算子に対してグループ化を実行し、関係するパウリ文字列のすべての期待値を評価するCommutablePauliSetMeasurementFactory
• ショットアロケータ: 測定ファクトリによって生成された CommutablePauliSetの間で合計ショットを分配する PauliSamplingShotsAllocator

どのように使用するかを示します。

from quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_sampler
from quri_parts.core.estimator.sampling import sampling_estimate

concurrent_sampler = create_qulacs_vector_concurrent_sampler()

estimate = sampling_estimate(
    op,            # Operator to estimate
    state,         # A circuit state
    5000,          # Total sampling shots
    concurrent_sampler, # ConcurrentSampler
    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement
    allocator,     # PauliSamplingShotsAllocator
)
print(f"Estimated expectation value: {estimate.value}")
print(f"Standard error of estimation: {estimate.error}")

#output
    Estimated expectation value: (0.7812513815547097+0.002275914633594121j)
    Standard error of estimation: 0.0707060360741753

ほとんどの場合、アルゴリズムにEstimatorやConcurrentEstimatorを組み込みたいと思うかもしれません。QURI Partsではcreate_sampling_estimatorとcreate_sampling_concurrent_estimatorを提供し、推定チュートリアルで紹介したvector estimatorと同じように使用できる（Concurrent）QuantumEstimatorを作成することができます。

from quri_parts.core.estimator.sampling import create_sampling_estimator
estimator = create_sampling_estimator(
    5000,          # Total sampling shots
    concurrent_sampler, # ConcurrentSampler
    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement
    allocator,     # PauliSamplingShotsAllocator
)
estimate = estimator(op, state)
print(f"Estimated expectation value: {estimate.value}")
print(f"Standard error of estimation: {estimate.error}")

#output
    Estimated expectation value: (0.7533110100454796-0.0441731501615594j)
    Standard error of estimation: 0.07035284315335183

concurrent sampling estimatorsでは、演算子と状態のセットをもう1つ作成します:

from quri_parts.core.state import GeneralCircuitQuantumState
from quri_parts.circuit import H, QuantumCircuit

op_1 = op
op_2 = Operator({
    PAULI_IDENTITY: 8,
    pauli_label("X0 Y1 Z2 X3"): 9 + 8j
})

state_1 = state
state_2 = GeneralCircuitQuantumState(4, QuantumCircuit(4, gates=[H(i) for i in range(4)]))

from quri_parts.core.estimator.sampling import create_sampling_concurrent_estimator

concurrent_estimator = create_sampling_concurrent_estimator(
    5000,          # Total sampling shots
    concurrent_sampler, # ConcurrentSampler
    bitwise_commuting_pauli_measurement, # Factory function for CommutablePauliSetMeasurement
    allocator,     # PauliSamplingShotsAllocator
)
estimates = concurrent_estimator([op_1, op_2], [state_1, state_2])

print(f"Estimated expectation value 1: {estimates[0].value}")
print(f"Standard error of estimation 1: {estimates[0].error}")
print("")
print(f"Estimated expectation value 2: {estimates[1].value}")
print(f"Standard error of estimation 2: {estimates[1].error}")

#output
    Estimated expectation value 1: (0.7328563810878067+0.02033531023801298j)
    Standard error of estimation 1: 0.0711939414351656
    
    Estimated expectation value 2: (7.8092-0.1696j)
    Standard error of estimation 2: 0.1702555909214144