Cirqを利用した古典コンピュータ上での量子回路シミュレーション

kenya-sk
17 min readJan 31, 2024

古典コンピュータ上でのSimulator

Cirqには、古典コンピュータ上で利用可能な2種類のシミュレータが用意されています。これらのシミュレータでは、比較的小規模の量子回路の動作をテストすることが可能です。より大規模な量子回路をシミュレーションしたい場合は、外部の高性能なシミュレータの利用を検討してください。これらのシミュレータはCirqのインターフェイスを提供しているため、容易に切り替えることが可能となっています。

  • cirq.Simulator
    - 純粋状態 (pure state)のシミュレーションを実行できる
    - 純粋状態とは、量子システムが単一の量子状態にあることを意味し、一つの純粋状態はベクトルで表現され、ユニタリ行列を用いて他の純粋状態へ変換できる
  • cirq.DensityMatrixSimulator
    - 混合状態 (mixed state)のシミュレーションを実行できる
    - 混合状態とは、異なる純粋状態が統計的な確率で混合した状態 (量子力学の重ね合わせとは異なる)を表し、一般的に密度行列を用いて表現される
    - 実際の量子デバイスで実行した時のように、外部環境の影響等のノイズを考慮したシミュレーション実行時に適している

今回は、cirq.Simulatorを中心にまとめております。cirq.DensityMatrixSimulatorについて後日公開させていただきます。

純粋状態のシミュレーション

純粋状態のシンプルな量子回路でシミュレーションを実行してみます。Simulatorの利用方法については、前の記事で紹介したものと基本的には同様となっております。今回のコードでは、基本となる回路を他の章でも活用するため、量子回路の構築部分を関数として切り出しています。状態ベクトルを出力するsimulateメソッドと観測結果のみを出力するrunメソッドに対応できるように、引数measで観測を行うかどうかを制御できるようになっています。

# 2量子ビットを定義
q0 = cirq.GridQubit(0, 0)
q1 = cirq.GridQubit(1, 0)

# 1. 各qubitにXゲートの平方根を適用
# 2. q0をcontrol、q1をtargetとして制御Zゲートを適用
# 3. 各qubitにXゲートの平方根を適用
# 4. meas=Trueの場合、各qubitを観測
def basic_circuit(meas=True):
sqrt_x = cirq.X**0.5
yield sqrt_x(q0), sqrt_x(q1)
yield cirq.CZ(q0, q1)
yield sqrt_x(q0), sqrt_x(q1)
if meas:
yield cirq.measure(q0, key="q0"), cirq.measure(q1, key="q1")

# 上記の回路を構築
circuit = cirq.Circuit()
circuit.append(basic_circuit())
print(circuit)

# 出力
# (0, 0): ───X^0.5───@───X^0.5───M('q0')───
# │
# (1, 0): ───X^0.5───@───X^0.5───M('q1')───

上記の回路をrunメソッドを利用して、シミュレーションを行い観測結果を確認してみます。runメソッドは結果をResultクラスとして返します。また、シミュレーションは量子コンピュータの挙動を模倣しているため、得られる観測結果は確率的となります。確率的な挙動は、numpyのrandomモジュールで制御されています。デバック等の目的で一時的に乱数のseedを固定したい場合は、Simulatorインスタンス作成時にnp.random.RandomStateで指定を行うことができます。

simulator = cirq.Simulator()
result = simulator.run(circuit)
print(result)

# 出力
# q0=1
# q1=0


# シミュレーションの乱数を固定した場合
simulator = cirq.Simulator(seed=np.random.RandomState(42))
result = simulator.run(circuit, repetitions=3)
print(result)

# 出力
# q0=011
# q1=110

量子回路の状態ベクトルを確認したい場合には、simulateメソッドを利用します。simulateメソッドは結果をSimulationTrialResultクラスとして返し、final_state_vectorプロパティに状態ベクトルを格納しています。量子デバイス上では結果として得られるのは観測結果のみとなるため、状態ベクトルはSimulator上でのみ確認できます。

# 量子回路の状態ベクトルを確認
circuit = cirq.Circuit()
circuit.append(basic_circuit(meas=False))
result = simulator.simulate(circuit, qubit_order=[q0, q1])
print(result.final_state_vector)

# 出力
# [0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ]

シミュレータはデフォルトで実数をnp.float32、複素数をnp.complex64のデータ型で確保します。より高精度な計算が必要な場合は、np.float64やnp.complex128にキャストしてシミュレータを実行することも可能です。

simulate_expectation_valuesメソッドを利用して、Operationの期待値を得ることも可能です。そのため、一度に複数のOperationの結果を扱いたい場合や量子アルゴリズムの評価時に有用です。

XX_obs = cirq.X(q0) * cirq.X(q1)
ZZ_obs = cirq.Z(q0) * cirq.Z(q1)
ev_list = simulator.simulate_expectation_values(
cirq.Circuit(basic_circuit(meas=False)), observables=[XX_obs, ZZ_obs]
)
# eval_list[0]->XX_obsの期待値、eval_list[1]->ZZ_obsの期待値
print(ev_list)

# 出力
# [(1+0j), 0j]

段階的なシミュレーション

デバックの際には、回路の最終的な結果だけではなく、各ステップでの挙動を確認する必要が出てきます。そこで、cirqではMoment毎の結果を返すsimulate_moment_stepsイテレータが用意されています。こちらを利用することでMomentを1ステップとみなし、段階的な中間結果を確認しながらシミュレーションを行えます。イテレータからは、各ステップの状態ベクトルを持つStepResultクラスが返されます。

# Momentを1stepとしてシミュレーション可能
circuit = cirq.Circuit()
circuit.append(basic_circuit(True))
print(circuit)
print("\n---------------- simulate step by step ----------------")
for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
print(f"state at step {i}: {np.round(step.state_vector(copy=True), 3)}")

# 出力
# (0, 0): ───X^0.5───@───X^0.5───M('q0')───
# │
# (1, 0): ───X^0.5───@───X^0.5───M('q1')───
#
# ---------------- simulate step by step ----------------
# state at step 0: [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j]
# state at step 1: [0. +0.5j 0.5+0.j 0.5+0.j 0. +0.5j]
# state at step 2: [0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ]
# state at step 3: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]

もし利用しているシミュレータがsimulate_moment_stepsをサポートしていない場合は、回路をチャンクに分割して実行することで同様のことを実現できます。しかし、コードの記述量が増え、可読性も下がるため、simulate_moment_stepsが使えない場合の代替手段とする。

# simulate_moment_stepsを利用した場合と同様のstepにチャンクを作成
chunks = [cirq.Circuit(moment) for moment in basic_circuit()]

# 前のチャンクで得られた状態ベクトルを保持し、次のチャンク実行時の初期状態として利用
prev_state = 0
for i, chunk in enumerate(chunks):
result = simulator.simulate(chunk, initial_state=prev_state)
prev_state = result.final_state_vector
print(f"state at step {i}: {np.around(prev_state, 3)}")

# 出力
# state at step 0: [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j]
# state at step 1: [0. +0.5j 0.5+0.j 0.5+0.j 0. +0.5j]
# state at step 2: [0.5+0.j 0. +0.5j 0. +0.5j 0.5+0.j ]
# state at step 3: [0.+0.j 0.+1.j 0.+0.j 0.+0.j]

パラメータを含むシミュレーション

ここまでは固定値を持つゲートを利用して回路を構築してきました。シミュレータの場合は、cirq.ParamResolverを利用することでパラメータを含む回路を定義することができます。cirq.ParamResolverはsympyライブラリを利用して定義されたSymbolを値にマップする役割を担います。このことにより、実験の管理が容易になります。

q0 = cirq.GridQubit(0, 0)
q1 = cirq.GridQubit(1, 0)

# パウリX演算子をx乗し、xをパラメータとして定義
root_w_gate = cirq.X ** sympy.Symbol('x')
circuit = cirq.Circuit()
circuit.append([root_w_gate(q0), root_w_gate(q1)])
print(circuit)

# パラメータxに0.0, 0.25, 0.5, 0.75, 1.0を実行時に代入してシミュレーション
for y in range(5):
resolver = cirq.ParamResolver({'x': y / 4.0})
result = simulator.simulate(circuit, resolver)
print(f"params: {result.params}, state vector: {np.around(result.final_state_vector, 2)}")

# 出力
# 0: ───X^x───
#
# 1: ───X^x───
# params: cirq.ParamResolver({'x': 0.0}), state vector: [1.+0.j 0.+0.j 0.+0.j 0.+0.j]
# params: cirq.ParamResolver({'x': 0.25}), state vector: [ 0.6 +0.6j 0.25-0.25j 0.25-0.25j -0.1 -0.1j ]
# params: cirq.ParamResolver({'x': 0.5}), state vector: [0. +0.5j 0.5+0.j 0.5+0.j 0. -0.5j]
# params: cirq.ParamResolver({'x': 0.75}), state vector: [-0.1 +0.1j 0.25+0.25j 0.25+0.25j 0.6 -0.6j ]
# params: cirq.ParamResolver({'x': 1.0}), state vector: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]

特定のパラメータ値の集合に対する一連の実験は、sweepと呼ばれる機能で定義することでより便利になります。sweepを利用してシミュレータを実行するときには、run_sweepメソッドを利用します。run_sweepはパラメータの値と観測結果を格納したResultクラスを全実験について保持したList形式のResultsを返します。

resolvers = [cirq.ParamResolver({'x': y / 2.0}) for y in range(3)]
circuit = cirq.Circuit()
circuit.append([rot_w_gate(q0), rot_w_gate(q1)])
circuit.append([cirq.measure(q0, key='q0'), cirq.measure(q1, key='q1')])
print(circuit)

# パラメータの集合はparams引数で渡し、
# 各パラメータについて、実験を2回繰り返す (repetitions引数で指定)
results = simulator.run_sweep(program=circuit, params=resolvers, repetitions=2)
for result in results:
print(f"params:{result.params}")
print(f"measurements:\n{result}\n")

# 出力
# 0: ───X^x───M('q0')───
#
# 1: ───X^x───M('q1')───
# params:cirq.ParamResolver({'x': 0.0})
# measurements:
# q0=00
# q1=00
#
# params:cirq.ParamResolver({'x': 0.5})
# measurements:
# q0=01
# q1=00
#
# params:cirq.ParamResolver({'x': 1.0})
# measurements:
# q0=11
# q1=11

上記ではパラメータのListを作成しsweepを構築しましたが、Cirqでは以下の方法でもsweepを構築することができます。より詳細な利用方法については、こちらの公式ドキュメントを参考にしてください。

  • cirq.Linespace
    - 特定の範囲を等間隔で分割した値の集合を作成
  • cirq.Points
    - 特定の値の集合を作成
  • cirq.Product
    - 複数のsweepに含まれる要素の可能な全ての組み合わせで構成される集合を作成
    - 「*」演算子で定義することも可能
  • cirq.Zip
    - 複数のsweepに含まれる要素を要素毎 (index毎)のペアとして集合を作成
    - 「+」演算子で定義することも可能
# 0.0 - 1.5までを4等分した値の集合を作成
params = cirq.Linspace(key="theta", start=0, stop=1.5, length=4)
for param in params:
print(param)

# 出力
# cirq.ParamResolver({'theta': 0.0})
# cirq.ParamResolver({'theta': 0.5})
# cirq.ParamResolver({'theta': 1.0})
# cirq.ParamResolver({'theta': 1.5})

# ----------------------------------
# 特定の値の集合を作成
params = cirq.Points(key="theta", points=[0, 1, 3])
for param in params:
print(param)

# 出力
# cirq.ParamResolver({'theta': 0})
# cirq.ParamResolver({'theta': 1})
# cirq.ParamResolver({'theta': 3})

# ----------------------------------
# パラメータの組み合わせから集合を作成
sweep1 = cirq.Linspace("theta", 0, 1.0, 3)
sweep2 = cirq.Points("gamma", [0, 3])
# for param in sweep1 * sweep2:と定義することも可能
for param in cirq.Product(sweep1, sweep2):
print(param)

# 出力
# cirq.ParamResolver({'theta': 0.0, 'gamma': 0})
# cirq.ParamResolver({'theta': 0.0, 'gamma': 3})
# cirq.ParamResolver({'theta': 0.5, 'gamma': 0})
# cirq.ParamResolver({'theta': 0.5, 'gamma': 3})
# cirq.ParamResolver({'theta': 1.0, 'gamma': 0})
# cirq.ParamResolver({'theta': 1.0, 'gamma': 3})

# ----------------------------------
# パラメータの要素積から集合を作成
sweep1 = cirq.Points("theta", [1, 2, 3])
sweep2 = cirq.Points("gamma", [0, 3, 4])
# for param in sweep1 + sweep2:と定義することも可能
for param in cirq.Zip(sweep1, sweep2):
print(param)

# 出力
# cirq.ParamResolver({'theta': 1, 'gamma': 0})
# cirq.ParamResolver({'theta': 2, 'gamma': 3})
# cirq.ParamResolver({'theta': 3, 'gamma': 4})

References

--

--