量子コンピュータを用いた高速数値積分
宇野 隼平iFast Numerical Integration Algorithms using Quantum Computer
Shumpei UNO
近年, 高速計算が可能な次世代のコンピュータの候補として, 量子コンピュータに大きな注目が集まっている. 本稿では, 量子コンピュータの今後の重要な応用先のひとつとして期待される, 数値積分を実行する量子アルゴリズムについて紹介する. 量子コンピュータを使用することで, これまでの代表的な数値積分手法 (モンテカルロ積分) に較べて, ほ
ぼ 2 乗の高速化が果たされる. 本稿では, 近年の著者らの研究を含めて, 3 種類の量子積分アルゴリズムを紹介する.

(キーワード): 量子コンピュータ, 量子アルゴリズム, 数値積分, 量子振幅推定, モンテカルロ積分

1 はじめに行っている.
量子コンピュータの重要な応用先の一つとして, 多

近年, これまでのコンピュータとは異なる原理に

次元数値積分が期待されている. 数値積分は, 例えば,より動作する量子コンピュータに大きな注目が集

金融工学 [5, 6], 物理シミュレーション [7, 8, 9], 安全

まっている. これまで, IBM, Intel, Google 等の大手

工学 [10, 11], 量子化学計算 [12, 13] を始めとして, 社

IT メーカーの他, 多くのベンチャー企業も量子コン
ピュータ開発への参入を表明しており, 開発競争が本格化し始めて来ている [1]. 既に, 小規模な量子コン

会的・学術的に幅広い分野で必要不可欠な計算である. 数値積分は, これまでのコンピュータ (以下, 「従来型のコンピュータ」という) では, 多くの場合, モン

ピュータをクラウドを通じて利用する環境が提供さ

テカルロ法を用いて計算されるが, 量子コンピュータ

れており [2, 3], 今後の急速な発展に備えて, 試験的を用いることで, モンテカルロ法と較べて 2 乗の速度に量子コンピュータを利用することも可能である.

向上が見込まれるアルゴリズムが知られている [14].

現在利用可能な量子コンピュータは, まだ小規模

近年, これらの数値積分の量子アルゴリズムを, 金融

でノイズの多いものであるが, 今後, もし大規模な工学を中心として, 具体的な問題へ適用するための量子コンピュータが実現すれば, これまで計算時間提案が行われている [15, 16, 17, 18]. このアルゴリの関係により, 解くことが難しかった問題を解決す

ズムの一部は, IBM がオープンソースとしている公ることが期待されており, 現在, 科学技術の様々な分開している量子コンピュータ向けの SDK(Software

野で, 量子コンピュータの活用に向けた検討が盛ん

Development Kit) の Qiskit [19] に, ライブラリとして組み込まれており, データを与えるだけで, 手軽に,に行われている. 例えば, 慶應義塾大学には, 量子コ
ンピュータを用いた実用的なアプリケーションの研

量子コンピュータの金融工学への応用を試験するこ究開発を, 産学連携により行うことを目的に掲げたとが出来る状況である.

IBM Q ネットワークハブが設立され [4], 金融業界からみずほフィナンシャルグループ, 三菱 UFJ フィ
ナンシャル・グループ, 化学業界から JSR, 三菱ケミ
カルの 4 社が参画している. 著者もみずほフィナン
シャルグループの一員として本活動に参画し, 様々なバックグラウンドを持つ研究者と議論することで, 量これまで提案されている多くの量子数値積分アルゴリズムでは Brassard らの量子振幅推定法 [20] が用いられている. Brassard らの量子振幅推定法は,比較的多くの演算操作が必要なアルゴリズムであり,非常に高い量子コンピュータの性能が要求される.
一方で, 現在の量子コンピュータでは, ビットの量子子コンピュータの将来的な応用方法について検討を状態を保つことの困難さ等に起因して, 可能な演算操
作の回数に制限がある. このため, なるべく早い時期

i サイエンスソリューション部 チーフコンサルタント 博士に量子コンピュータを活用するためには, Brassard

(理学)

1

らの量子振幅推定法 [20] に較べて演算操作の回数のまた, p(x) は ρ(x) を, I˜d を M d 個に等分割した各超少ないアルゴリズムを開発することが重要な課題で立方体上で積分した関数であり,ある. この課題に対する筆者らの最近の取り組み [21]満たす. f (x) は, h(x) の格子点上での値である:

本稿の構成は以下のとおりである. まず, 2 節では,
従来型のコンピュータにおいて代表的な数値積分手

3 節において, 量子コンピュータを用いた 3 種類の数
値積分手法を紹介する.

2 従来型のコンピュータによる積分計算
本節では, 数値積分問題の定式化及び従来型のコンピュータによる計算手法の概要について述べる. な
お, ここでの数値積分問題の定式化は主に [22, 23] に従う.

d 次元の数値積分問題は, 有界な実関数 g : Rd → R
の積分

∫
(1)

Id

が解析的に評価出来ない場合に, 数値的に近似値を求める問題である. ここで Id は d 次元空間内の有界な部分空間 Id ⊂ Rd である.
式 (1) は一般性を失わずに, 上界, 下界から決まる
既知の定数によりスケールすることで, 以下の形に書き直すことが出来る.

∫
I(h) =

ρ(x)h(x)dx

(2)

I˜d

ここで, I˜d = [0, 1]d は d 次元超立方体, h は 0 から 1の間に値域ii をとる実関数 h : I˜d → [0, 1] であり, また, ρ は以下を満たす実関数である:

1 ∑h(xi )
I(h) ∼
N i=1

ρ(x)dx = 1

(3)

I˜d

リズムは, 積分の各次元を M 個に分割し, 矩形近似
することである. この時, 積分 I(h) の近似は

p(x)f (x)(6)

Chebyshev の不等式により, サンプルの数 N (被積分関数 h の評価回数) と誤差 ε の関係は, N ∼ O(ϵ−2 )となることが知られている [6]. モンテカルロ法は式(4) のような直接的な和の評価に較べて, 遥かに効率的な手法であるが, 一方で, 実用上は, N ∼ O(ϵ−2 )の値は大きくなることがあり (例えば, 誤差を ε ∼O(10−5 ) にしようと思えば, N ∼ O(1010 ) 程度のサンプル数が必要になる), モンテカルロ法による数値

式 (2) の近似値を求めるための最も単純なアルゴ

∑

(5)

N

∫

I(h) ∼ S(f ) =1

1 次元の積分区間を M 等分に分割した場合に, 誤()
差が ε ∼ O M −1 で表せることから, 式 (4) を用いて誤差 ε を達成するためには, 被積分関数 f をO(ε−d ) 回程度評価する必要がある. 典型的な金融の問題では, 積分の次元 d(金融資産の種類) は数百〜数千次元であるため, 例えば, d = 100 として, 誤差をε ∼ O(10−3 ) 程度まで評価するとすると, 積分の近似
値 S(f ) を計算するために, 被積分関数 f を O(10300 )回程度評価する必要がある. これは, f が 1 回の浮動小数点演算で計算できたとしても, 1PFlops(1 秒
間に 1015 回の浮動小数点演算が可能) のスパコンを使って, O(10285 ) 秒 ∼ O(10278 ) 年かかる計算である. 現在の宇宙年齢が O(1010 ) 年であることを鑑みると, 矩形近似 (4) を使って, 積分計算を実行することは実質上不可能である. この困難は主に, 演算回数
が次元 d の増加に対して, 指数関数的に増加する事により生じている.
以上の事情により, 高次元の積分を評価する場合には, 式 (4) を直接評価することは無く, モンテカルロ
法により評価を行うのが一般的である. モンテカル
ロ法では, 式 (2) の p(x) を確率密度関数とみなすことで, 式 (2) を統計的期待値と見なして計算を行う.
例えば, p(x) に従う確率変数 x の N 個のサンプル
x1 , x2 , · · · , xN が得られたとき, (2) の近似値は以下の式により与えられる.

法であるモンテカルロ法について紹介する. その後,

g(x)dxx∈Jd p(x) = 1 をx2
xd ), ,··· ,f (x1 , x2 , · · · , xd ) = h
M M
M
(x

についても本稿で簡単に紹介する.

I(g) =

∑

(4)

x∈Jd

と和の形で計算することが出来る. ここで, Jd は M dd個の整数格子点 Jd = {0, 1, 2, · · · , M − 1} である.

積分評価は様々な分野で多くの時間を占める計算の一つである.

ii この段階では h の値域を [0, 1] の範囲に制限する必要は無

以下で説明するように, 量子コンピュータを用い

いが, 後の量子アルゴリズムの際の便利のため, 値域を設定した.

たアルゴリズムでは, 被積分関数 f の評価回数 N と2

誤差 ε の関係は N ∼ O(ϵ−1 ) となり, 従来型のコン

量子コンピュータの演算は, 量子コンピュータの状

ピュータのモンテカルロ計算と較べて, 効率的に数値

態 |ψ⟩n を別の状態 |ψ ′ ⟩n に遷移させるユニタリ演算

積分を行うことが可能である.

子 U として表される.

3 量子アルゴリズム|ψ ′ ⟩n = U |ψ⟩n(9)

本節では量子コンピュータを用いて数値積分を行う 3 種類のアルゴリズムについて紹介する.

本稿では, 量子コンピュータの計算は, 初期状態 |0⟩n

まず, 3.1 節において量子計算計算モデルの基本事に対して, 様々なユニタリ演算子 U1 , U2 , · · · を作用

項をごく簡単に述べる. 次に, 3.2 節において, 現在広させることで進んでいくこととする. なお, 本稿で用く用いられている Brassard らの手法 [20] を用いたいるユニタリ演算子は付録 A にまとめる.

数値積分手法について述べる. その後, 3.3 節におい最後に, 量子コンピュータの計算結果の読み出して, Brassard らの手法に較べて演算数が少ないと考(測定) について述べる. 本稿では, 計算終了時の量子∑状態が |ψ⟩n = x∈{0,1}n αx |x⟩n で表されるとした
とき, 計算結果の読み出しを行うことで, n ビット列x ∈ {0, 1}n が確率 |αx |2 で得られることとする.
以下では, この計算モデルを用いた数値積分計算の概要を述べる.

えられる Abrams の手法 [22] について述べる. 最後に, 3.4 節において, 本稿の著者を含めた共同研究である, 鈴木らの手法 [21] について述べる.

3.1

量子計算モデルの基本事項

ここでは, 以下で用いる (ゲート型) 量子コンピ
ュータの計算モデルについて, ごく簡単に述べる. よ3.2 Brassard らの量子振幅推定法を用いた数値積分
本小節では, 数値積分が量子振幅推定の問題に書き直せることを利用して, 数値積分を行う量子アルゴ
リズムについて述べる. 本小節の記載は主に [15] に従う.
数値積分の量子アルゴリズムでは, 式 (4) の数値
積分の矩形近似値を算出することを目的とする. なお, 以下では, 表記の簡潔さのため, 数値積分の次元d を 1 次元として取り扱い, 式 (4) の積分区間の分割
数 M = 2n 個とした近似値

り詳細な内容については, 教科書 [24, 25] 等を参照されたい.
従来型のコンピュータでは, ビットは 0 又は 1 の確定した状態を取るが, 量子コンピュータのビット(以下, 「量子ビット」という) は重ね合わせ状態を取ることが出来る. より具体的には, 量子ビットの状態|ϕ⟩ は, 2 次元ヒルベルト空間 H2 = C2 内の単位ベク
トル|ϕ⟩ = α |0⟩ + β |1⟩(7)で表される. ここで, α, β は |α|2 + |β|2 = 1 を満たす

S(f ) =

任意の複素数であり, H2 の正規直交基底 |0⟩ , |1⟩ は

n
2∑
−1

p(x)f (x)

(10)

x=0

計算基底と呼ばれる.
量子コンピュータの状態は, 複数の量子ビットを合を求めることとするが, 任意の有限次元 d の積分へのわせた状態として表される. 具体的には, n 個の量子拡張は容易である.

ビットを組み合わせた状態 |ψ⟩n として, n 個の H2

付録 A の式 (58) 及び (59) から構成される演算子のテンソル積の状態空間 H2⊗n = C2 内の単位ベクnトルとして表される.

∑|ψ⟩n =A = R(P ⊗ I1 )αx |x⟩n(8)を計算の初期状態 |0⟩n+1 に作用させると

x∈{0,1}n∑x∈{0,1}n |αx |= 1 を満たす 2n 個の複素数の組であり, 確率振幅と呼ばれる. また, |x⟩は, 計算基底 |0⟩ , |1⟩ の n 個の積状態である. 以下,
x ∈ {0, 1}n を, 非負の整数 {0, 1, 2, · · · , 2n−1 } の 2
進数表記と見なすこととする.
ここで, αx は(11)2|Ψ⟩n+1 = A |0⟩n+1
=

n
2∑
−1 √

p(x) |x⟩n(√)√
f (x) |0⟩ + 1 − f (x) |1⟩

x=0

(12)
3

が得られる. 正規直交基底 |Ψ̃0 ⟩n+1 及び |Ψ̃1 ⟩n+1 を|Ψ̃0 ⟩n+1 =

n
2∑
−1 √

x=0

|Ψ̃1 ⟩n+1 =

√

p(x)

n
2∑
−1 √

p(x)

x=0

√

f (x)
|x⟩n |0⟩
S(f )
(13)

Q = UΨ UΨ̃0
UΨ̃0 = In ⊗ Z()
UΨ = A In+1 − 2 |0⟩n+1 ⟨0|n+1 A†

1 − f (x)|x⟩n |1⟩
1 − S(f )と導入すると, 式 (12) の |Ψ⟩n+1 は|Ψ⟩n+1 =√S(f ) |Ψ̃0 ⟩n+1 +√√S(f ) = sin θ, (0 ≤ θ ≤ π2 ) と置く. 式 (14)の |Ψ⟩n+1 に作用する, 量子振幅増幅法の演算子 Q
は, 以下のように定義される.

ため,

(15)

振幅増幅演算子 Q を式 (14) の |Ψ⟩n+1 に j 回作用させると,

1 − S(f ) |Ψ̃1 ⟩n+1
(14)

Qj |Ψ⟩n+1と書き直すことが出来る. 式 (14) では, 数値積分の= sin((2j + 1)θ) |Ψ̃0 ⟩n+1 + cos((2j + 1)θ) |Ψ̃1 ⟩n+1

近似値 S(f ) の平方根が |Ψ̃0 ⟩n+1 の確率振幅として(16)現れており, もし確率振幅を十分な精度で求めることが出来れば, 数値積分の近似値を求められることにとなり, Q は式 (12) で張られる部分空間 HΨ 内の回なる.

転演算子として作用することがわかる. 部分空間 HΨ
内における Q の固有状態 |Ψ± ⟩ は,今, 仮に, 式 (14) の |Ψ⟩n+1 を繰り返し測定することで確率振幅を推定することを考える (以下, 2 種類)
1 (|Ψ± ⟩n+1 = √ |Ψ̃1 ⟩n+1 ∓ i |Ψ̃0 ⟩n+12の測定結果を取る独立な試行を「ベルヌーイ試行」という). この時, ベルヌーイ分布の標準偏差の式によ

(

(17)

)

と構成することが可能であり, それぞれ固有値 e±2iθであり, 従来型のコンピュータによるモンテカルロ法に対応する. もし, 振幅増幅演算子 Q の固有値 e±2iθ

と同程度のサンプル数が必要となる. このため, 量子のどちらかを推定することが出来れば, 積分の近似値

コンピュータを用いることによる加速を得るためにS(f ) = sin2 θ を計算出来ることとなる.

り, サンプル数 N と誤差 ε の関係は, N ∼ O ε−2

量子コンピュータでユニタリ演算子の固有値をは, 確率振幅を推定するための何らかの工夫が必要と求める方法として, 量子位相推定法が知られている

なる.

[24, 27]. 量子位相推定法は, ユニタリ演算子とその固有ベクトルが与えられた時に, 効率的に固有値を算出する手法である. 量子位相推定法を適用するため,
式 (14) の |Ψ⟩n+1 を, Q の固有ベクトル |Ψ± ⟩n+1 で展開すると,

単純なサンプリングに比べて効率的に確率振幅を推定する方法としては, Brassard らによって提案された手法が知られている [20]. 以下, モンテカルロ法
における被積分関数 f の評価回数 N に対応して, 量子計算における関数の呼び出し回数 N を式 (11) で定義される演算子 A の演算回数と定義する. この時,)
1 (|Ψ⟩n+1 = √ eiθ |Ψ+ ⟩n+1 + e−iθ |Ψ− ⟩n+1 (18)2

Brassard らの手法を用いた確率振幅の推定では, 誤
差 ε を達成するための演算子 A の呼び出し回数 N()
は, N ∼ O ε−1 iii であり, モンテカルロ法のほぼ平方根程度の関数の呼び出し回数で, 同程度の誤差を達成することが出来る. これは, 誤差を ε ∼ O(10−5 )にしようと思えば, N ∼ O(105 ) 程度の関数の呼び出
し回数で良いこととなり, モンテカルロ法と較べて大幅な速度向上が見込まれる.
以下に, Brassard らの振幅推定法を用いた S(f ) の推定方法の概要を述べる.
Brassard らの振幅推定法は, 量子振幅増幅法 [20,26] と位相推定法 [24] から構成される. 以下, 便宜のと書き直せる.
次 に, |Ψn+1 ⟩ に 使 用 す る n + 1 ビ ッ ト に 加 えて, 新たに m 個の補助ビットを追加した量子状態|0⟩m |Ψ⟩n+1 を考える. 量子位相推定法では, この追加した m 個の補助ビットに位相の推定結果が出力される. まず, 状態 |0⟩m |Ψ⟩n+1 の m 個の補助ビットに対して, Hadamard 演算子を作用させると, 量子状
態はH ⊗m |0⟩m |Ψ⟩n+1
2 −1
(
)
1 ∑
|y⟩m eiθ |Ψ+ ⟩n+1 + e−iθ |Ψ− ⟩n+1=√
2m y=0m

iii 本稿の計算量オーダーは, 対数の因子を除いたものである

4

の誤差で得られたとした時, S(f ) の推定値 Ŝ(f ) の(19)

推定誤差は, 誤差の伝播により,となる. 式 (19) の量子状態に, m 量子ビット |y⟩m を√制御ビットとする制御ユニタリ演算子 c Q を作用さ|S(f )−Ŝ(f )| ≤ 2πk

せると,
2 −1
(
)
1 ∑
Q√
|y⟩m eiθ |Ψ+ ⟩n+1 + e−iθ |Ψ− ⟩n+1m
2 y=0

( π )2
S(f ) (1 − S(f ))
2
+k
2m
2m
(22)

m

c

m
−1
iθ 2∑

e
=√

2 y=0

と求められる.
以上の手法において, 演算子 A の呼び出し回数 N
は, 演算子 Q の中に A が 2 回含まれることを考慮す

e2iyθ |y⟩m |Ψ+ ⟩n+1

ると, N = 2m+1 −1 回となる. よって, この手法における誤差 ε = |S(f ) − Ŝ(f )| と関数 f の呼出回数 N

(

2m −1

e−iθ ∑ −2iyθe|y⟩m |Ψ− ⟩n+1
+ √
2 y=0

以上の説明では簡単のため, 数値積分の次元 d = 1としたが, 以上の計算過程から明らかなように, 量子(20)振幅増幅法においても, モンテカルロ法と同様に, 関が得られる. ここで, 周波数解析等を思い起こして,数の呼出回数 N は次元 d に依存しないことに留意

離散フーリエ変換を行うことで振動数を取り出せるする.

ことの類推から, 式 (20) の m 量子ビットに付録 A

以上で概要を与えた Brassard ら [20] の振幅推定−1の式 (60) の量子逆フーリエ変換演算子 Fmを作用法では, 多くの制御ユニタリ演算子 c Q を用いていさせた量子状態る. 現在, 利用可能な量子コンピュータ [2, 3] では Q
等の単なるユニタリ演算子に較べて, c Q 等の制御ユ

2 −1eiθ ∑ 2iyθ −1√eFm |y⟩m |Ψ+ ⟩n+1
2 y=0
m

m
−1
−iθ 2∑

e
+ √

2

)

の関係は N ∼ O ε−1 であることがわかる. なお,

ニタリ演算子は演算数が多く, ノイズが生じやすい傾
向にあり, この傾向は今後もしばらくは続くと考えられる. このため, 近い将来の量子コンピュータの活用−1e−2iyθ Fm|y⟩m |Ψ− ⟩n+1

y=0に向けて, 制御ユニタリ演算子 c Q をなるべく使用し

(21)

2m −1

eiθ ∑
=√
α+ (x) |x⟩m |Ψ+ ⟩n+1
2 x=0ないアルゴリズムの開発は重要な課題である. 以下,
制御ユニタリ演算子 c Q を使わずに量子振幅を行う

2 種類のアルゴリズムを紹介する.

2m −1e−iθ ∑α− (x) |x⟩m |Ψ− ⟩n+1
+ √
2 x=0
(

の係数 α+ (x) 及び α− (x) は, 2m πθ 及び 2m 1 − πθ

3.3 Abrams らの量子振幅推定法を用いた数値積分
本小節では, Abrams らの方法 [22] に基づき, 量子
位相推定法を使用せずに, 量子振幅増幅法のみで振幅推定を行い, 式 (10) の S(f ) を求める手法について述べる. この手法では, 量子位相推定を使わないため,
3.2 節の方法で使われた制御ユニタリ演算子 c Q や量−1を使う必要が無く, 演算数の子逆フーリエ変換 Fm削減が期待される. ただし, 後述するように, Abrams
らの方法では, 1 種類の振幅増幅演算子 Q のみでな
く, 測定結果に依存して, 様々な演算子を用いて振幅増幅を行う必要があることに留意する.
Abrams らの振幅推定法の大まかな概念を図 1 に示す. この手法は, Iteration によって, 徐々に正確な振幅を推定していくアルゴリズムである. まず, S(f )を埋め込んだ初期状態に対して, ベルヌーイ試行を繰り返し, 振幅 S(f ) の大きさを適当な誤差範囲で定める. 次に, この誤差範囲を振幅増幅法で拡大した上)付近にピークを持つ関数であると予想できる. もし,
この付近に十分高いピークを持つとすれば, m 個の補助ビットの量子状態 |x⟩m を観測すれば, θ または π − θ の近似値を得ることが出来るため, S(f ) =sin2 (θ) = sin2 (π − θ) から, 数値積分の近似値を推定することが出来る.
実際, 式 (21) の状態に対して, m ビット量子状態|x⟩ を観測したときに, θ または π − θ が ±π/2m 以()
内の誤差で得られる確率の下限値 π82 ∼ 81% が,Brassard ら [20] により求められている. 同様に, kを適当な整数とした時, θ または π − θ が ±πk/2m 以∑k
1内の誤差で得られる確率の下限は, π82
i=1 (2i−1)2と求めることが出来る (例えば, k = 5 とすれば成功確率を 95% 以上). θ または π − θ が ±πk/2m 以内
5

図 1 Abrams らの振幅推定法 [22] の概念図. ベルヌーイ試行による誤差範囲の設定と, 振幅増幅法による誤差範囲の拡大を繰り返す Iteration を行うことで, 正確な振幅の推定値を求めていく.

で, ベルヌーイ試行を行い, より正確な誤差範囲を定

ここで, 式 (25) の量子状態 |Ψ1 ⟩ に対して, N1 回のめる. Abrams らの手法は, この誤差範囲の設定と振ベルヌーイ試行をした場合に, S(f ) の推定値 Ŝ1 (f )幅増幅法を繰り返す Iteration を行うことで, 徐々にがどの程度の誤差範囲で求まるかを算出する. 試行

正確な振幅の推定値を求めていく手法である. なお,

回数 N , 確率 p のベルヌーイ分布に従う確率変数の後述するように, この手法において, 関数の評価回数

誤差範囲 εp が, N が十分大きい場合に, 適当な定数√()
N と誤差 ε の関係は, N ∼ O ε−1 と Brassard らの振幅推定法と同程度である.
本小節では, まず, 3.3.1 節において, Iteration の 1
回目のアルゴリズムを示し, 続いて, 3.3.2 節において k 回目の Iteration のアルゴリズムを示す. 最後に, 3.3.3 節において, 関数の評価回数 N と誤差 ε の関係式を導く.

c1 と標準偏差√
εp = c1

Rf |x⟩n |0⟩ = |x⟩n f (x) |0⟩ +

√

√
c1

Ŝ1 (f ) の推定誤差 ε1 は√√1
c1 1 − S(f )2
c1
ε1 =
≤
2
N1
2 N1

(28)

3.3.2 節で述べるように, Amrams らの手法では,
こうして求められた誤差範囲 Ŝ(f ) − ϵ21 ≤ S(f ) ≤Ŝ1 (f ) + ϵ21 を振幅増幅法を用いて拡大することで, 更に精度のよい推定を行う手法である.
式 (28) により, ε1 の値と c1 の値を定めれば, Iteration 内で必要なベルヌーイ試行の回数 N1 を求める
ことが出来る. c1 の値は, 3.3.3 節で議論するように(24)を計算の初期状態 |0⟩n+1 に作用させると|Ψ1 ⟩n+1 = Bf |0⟩n+1
p(x)f (x) |0⟩n+1 + · · ·(27)と求められる.

式 (23) と付録 A の式 (58) から構成される演算子=

S(f )2 (1 − S(f )2 )N1で与えられる. 誤差の伝播公式を考慮すれば, 推定値

(23)

n
2∑
−1

(26)

回測定した場合の S(f )2 の推定誤差は,)
1 − f (x)2 |1⟩

Bf = (P † ⊗ I1 )Rf (P ⊗ I1 )p(1 − p)Nと表されることを考慮すると, 量子状態 |Ψ1 ⟩ を N1

3.3.1 Iteration の 1 回目
本手法では, 関数 f (x) を計算するため, 付録 A の式 (59) と類似した以下の演算子を仮定する.
(

p(1−p)
を用いて,
N

(25)

x=0

本手法の成功確率が有限になるように定める必要がが得られ, |0⟩n+1 の確率振幅として, S(f ) が得らある. また, 設定する誤差 ε1 は, 最終的に求めたい誤れる.

差 ε に較べて遥かに大きい値をとることに留意する

6

(例えば, 最終的に求めたい誤差 ε ∼ 10−5 に対して
ε1 = 12 など).
なお, 以下では簡単のため, 全ての Iteration でベ
ルヌーイ試行の推定誤差を一定の値 ε0 に設定するとして議論を進めていく.

が得られ, |0⟩n+1 の確率振幅として, Dk (f ) が得られる. 以下, 被積分関数 f (x) の計算が引き算に較べて遥かに複雑であるとし, Bf と Bfk の総演算回数を被積分関数 f の評価回数 N として数えることとする.
ここで, 0 ≤ D(f ) ≤ εk−1であることを考慮して,0
式 (33) の |0⟩n+1 の振幅増幅を行う. 以下, 便宜のた√
Dk (f ) = sin θk , (0 ≤ θk ≤ π2 ) と置く. 振幅増幅演算子 Qk を

Iteration の k 回目
ここでは, Abrams らの振幅推定手法 [22] の k 回
目の Iteration の手続きを示す.
Iteration の k − 1 回 目 ま で で, S(f ) の 推 定 値Ŝk−1 (f ) が誤差 εk−1の範囲で求まっていると仮定す
0
る. この時, 以下の手続きを行うことで, Iteration の
k 回目において, ほぼ定数回のベルヌーイ試行を行う
ことで, 振幅 S(f ) の推定値 Ŝk (f ) を誤差 εk0 の範囲で求めることが出来ることを示す.
Iteration の k 回 目 の 手 続 き で は, Ŝk−1 (f ) −
3.3.2

εk−1
0
2

≤ S(f ) ≤ Ŝk−1 (f ) +

εk−1
0
2

め,

Qk = UΨk UΨ̃k0
UΨ̃k0 = In+1 − 2 |0⟩n+1 ⟨0|n+1()
UΨk = Bk In+1 − 2 |0⟩n+1 ⟨0|n+1 Bk†と 定 義 す る と, 振 幅 増 幅 演 算 子 Qk を 式 (33) の|Ψk ⟩n+1 に jk 回作用させることで,
Qjkk |Ψ⟩n+1 = sin((2jk + 1)θk ) |0⟩n+1 + · · · (35)の範囲の確率振幅を振幅増幅法で拡大することで推定値 Ŝk (f ) を求と, |0⟩n+1 の振幅を増幅させることが出来る. jk のめる. この目的のため, まず, 以下の量 fk (x) 及び値は, 振幅が最も増幅するように,

Dk (f ) を定義する:(
fk (x) = f (x) −Dk (f ) =

n
2∑
−1

εk−1
Ŝk−1 (f ) − 0
2

(34)

)

(2jk + 1) sin−1 εk0 ∼
(29)

π
2

(36)

と選ぶこととする. 式 (35) の量子状態に対して Nk
回の測定を行った時の sin2 ((2jk + 1)θk ) の推定値がp(x)fk (x)ベルヌーイ分布に従うことから, 誤差の伝播公式を考

x=0

もし, こうして定義した Dk (f ) の推定値 D̂k (f ) を,慮すると, Dk (f ) = sin2 θk の推定誤差は, 適当な定誤差 εk0 の範囲で求めることが出来れば, S(f ) の推定

数 ck を用いて,(

値 Ŝk (f ) を誤差 εk0 で以下のようにして求めること

ck cos

が出来る:

εk =

Ŝk (f ) = Ŝk−1 (f ) −εk−1
0+ D̂k (f )
2

(30)

(38)

なので, 誤差範囲を εk = εk0 にするために必要なベル
ヌーイ試行の回数 Nk は, 式 (37) から, k が大きい場と付録 A の式 (58) から構成される演算子

合には, ck の因子を除いて, ほぼ定数iv になることが

(32)

わかる.
以上により, ほぼ定数回のベルヌーイ試行により,を定義する. 式 (32) を計算の初期状態 |0⟩n+1 に作用

k 回目の Iteration において, S(f ) の推定値 Ŝk (f ) を
誤差 εk0 で求める手続きを示した.

させると|Ψk ⟩n+1 = Bfk |0⟩n+1
p(x)fk (x) |0⟩n+1 + · · ·

ck
√
(37)
2(2jk + 1) Nk()
jk ∼ O ε−k
0

(31)

=

2(2jk + 1) Nk≤

式 (36) から, εk0 が十分小さければ,()√
Rfk |x⟩n |0⟩ = |x⟩n fk (x) |0⟩ + 1 − fk (x)2 |1⟩

n
2∑
−1

√

)

と表される.

式 (29) の Dk (f ) を推定するため,

Bfk = (P † ⊗ I1 )Rfk (P ⊗ I1 )sin−1 (Dk (f ))2jk +1(33)
iv ここでは, k や ε に依存しないという意味での定数

x=07

3.3.3 関数の呼出回数と誤差の関係
ここでは, 以上により概要を示した Abrams の振幅推定法 [22] における関数 f の呼出回数 N と誤差()
ε の関係が , (対数因子を除いて)N ∼ O ε−1 により与えられることを示す.
m 回目の Iteration により, 最終的に要求される精度 εm0 = ε に達したとする. この時, Abrams の振幅推定法における関数 f の総呼び出し回数は, k 回目のIteration における振幅増幅の回数 jk とベルヌーイ
試行の回数 Nk により, 以下のように与えられる.
N=m∑(2jk + 1)Nkを用いると式 (42) から, ベルヌーイ試行の回数 Nkと誤差 ε の関係が Nk ∼ O(ln(ε−1 )) 程度であることがわかる. 以上により, f の全評価回数 N と最終的な誤差 ε の関係は, 対数の因子を除いて N ∼ O(ε−1 )程度であることがわかる.
以上に示した Abrams らの振幅推定法 [22] は, 量子位相推定法を使わない手法であるため, Brassardらの手法 [20] で必要であった制御ユニタリ演算子
c

Q や逆フーリエ変換演算子 Fm が不要となっている. このため, Abrams らの手法は, 近い将来の量子コンピュータの活用に向けて, Brassard らの手法と較べて, 有効な手法である. 一方で, (32) のように,測定結果に応じて, 振幅増幅演算子 Qk を変更する必要があるため, 実用上, 演算子の構成に伴う困難が生じる可能性がある. 以下, 3.4 節では, 1 種類の振幅増幅演算子のみを用いて, 振幅推定を行う鈴木らの最近の研究 [21] について述べる.

(39)

k=1

こ こ で, jk の ε 依 存 性 は, 式 (38) に よ り, jk ∼()
O εm−kε−1 で与えられる. 以下, 残りの部分, ベル
0
ヌーイ試行の回数 Nk と誤差 ε の関係を示す.
この目的のため, m 回目の Iteration の後の S(f )の推測値 Sm (f ) が, 誤差 ε 以内に入る確率が一定の
値 Cp v 以上であることを要請するとする. 簡単のため, 各 Iteration 内の失敗確率を一定の値 Pfail となるとする. この時, m 回の Iteration 後の成功確率が,Cp 以上であるとすると,
m(1 − Pfail )≥ Cp

鈴木らの振幅推定法を用いた数値積分

3.4

この小節では, 鈴木らによる, 振幅増幅法と最尤推定を組み合わせた振幅推定法 [21] について述べる.
なお, 本研究は, 慶應義塾大学の IBM Q ネットワー
クハブにおいて, 本稿の著者を含めた共同研究として行った研究である.

(40)

鈴木らの振幅推定法の概念図を図 2 に示す. 大まが成り立つ. ここで, 平均値 µ, 標準偏差 σ に従う確

かに言えば, 鈴木らの手法は, 推定したい量 S(f ) と率変数 X に対する Chebyshev の不等式

相関するような様々な量子状態について測定を行い,

1
Pfail (|X − µ| ≥ ck σ) ≤ 2
ck

得られる情報を組み合わせることで, S(f ) を推定する手法である. より具体的には, 振幅増幅法を用いて(41)

S(f ) に相関する様々な量子状態を作り, 得られた測定結果から最尤推定により, S(f ) の推定を行う手法
である.
以下, 3.4.1 節において本手法の概要を, 3.4.2 節において本手法の関数の呼出回数と誤差の関係を示す.

を考慮すると, 式 (40) を満たすためには次の関係式が成り立てば良いことがわかる.

1
ln
m

(

1
Cp

)
∼

1
c2k

(42)

3.4.1

ここで, 各 Iteration の誤差 ε0 と最終的な誤差 ε の手法の概要

以下, 鈴木らの手法を用いて, 式 (14) の |Ψ̃0 ⟩n+1

関係式√
S(f ) = sin θ を推定する方法を述べる.
本手法においても, 3.3 節において概要を説明した
Abrams らの手法 [22] と同様に, Iteration により,√S(f ) を求める. k 回目の Iterationvi では, 式 (14)

の確率振幅

εm
0 =ε

(43)

及び k 回目の Iteration におけるベルヌーイ試行の誤差√
ε0 = cksin2 ((2jk + 1)θk ) cos2 ((2jk + 1)θk )(44)
Nkvi k 回目の振幅増幅回路の実行に際しては k − 1 回目より前の結果を使用しないため, 独立に実行することが可能であり, 正確には Iteration では無いが, ここでは Abrams らの手法と対比させるため, Iteration と呼ぶこととした.

v 例えば 1 等の適当な定数38図 2 鈴木ら [21] の振幅増幅法の概念図. 各 Iteration では, jk 回の振幅増幅を行った量子状態 Qjk |Ψ⟩ に対する測定結果をもとに, 尤度関数 Lk を作成する. 全ての Iteration の尤度関数の積から全体の尤度関数
L を求めることで, L を最大にする引数として, 推定値 Ŝ(f ) が得られる.

定値 θ̂ を用いて, S(f ) の推定値を Ŝ(f ) = sin2 θ̂ とに対して, jk 回の振幅増幅を行った状態

Qjk |Ψ⟩n+1 = sin((2jk + 1)θ) |Ψ̃0 ⟩n+1
+ cos((2jk + 1)θ) |Ψ̃1 ⟩n+1

求めることが出来る.
ここで, 0 ≤ θ ≤ π/2 の範囲では, θ と S(f ) =sin2 θ が 1 対 1 対応することから, 尤度関数 L をS(f ) の関数としても見なすことが出来ることを考慮
して, L(⃗h; S(f )) = L(⃗h; θ) と書くとすると, 式 (48)(45)

に対して, Nk 回のベルヌーイ試行を行うとする. この時 |Ψ̃0 ⟩n+1 が観測される回数が hk 回であったとによる推定は

すると, θ に対する尤度関数 Lk (hk ; θ) は (定数項を除いて) 以下の式で表される.

Ŝ(f ) = arg max ln L(⃗h; S(f ))

(49)

S(f )

Lk (hk ; θ)[
]hk [ 2
]Nk −hk= sin2 ((2jk + 1)θ)cos ((2jk + 1)θ)と見なすことも出来ることに留意する.
以上により, 最尤推定法による S(f ) の予測値 Ŝ(f )(46)の推定方法について述べた. 以下では, 本手法における関数 f の呼出回数 N と誤差 ε の関係について議論

以上の手続きを k = 1 から m まで, 適当な振幅増幅の回数 ⃗j = {j1 , j2 , · · · , jm } 及びベルヌーイ試行の回する.

⃗ = {N1 , N2 , · · · , Nm } に対して行い, |Ψ̃0 ⟩数N
n+1

3.4.2 関数の呼出回数と誤差の関係
ここでは, 鈴木らの手法 [21] における, 関数 f の呼出回数 N と誤差 ε の関係について, フィッシャー情報量を用いて議論する. なお, フィッシャー情報量の詳細については, 統計学の教科書 [28] 等を参照されたい.
一般に, 確率変数 x が, 母数 a により決まる確率分
布 L(x; a) に従うとき, フィッシャー情報量 I(a) は
以下により与えられる.
[(
)2 ]
∂
I(a) = E(50)
ln L(x; a)∂aが観測された回数が ⃗h = {h1 , h2 , · · · , hm } と得られたとする. この時, 全観測結果に対する尤度関数は式(46) を用いて, 以下のように構成することが出来る.
L(⃗h; θ) =m∏

Lk (hk ; θ)

(47)

k=1

このように構成された尤度関数を用いて, θ の最尤推定を行う. 最尤推定値 θ̂ は, 尤度関数を最大化する θとして与えられる:θ̂ = arg max L(⃗h; θ)θ= arg max ln L(⃗h; θ)(48)

ここで, 期待値 E は確率変数 x について取るものθとする. a の任意の不偏推定量 â の推定誤差 ε =√E[(â − a)2 ] について, 推定誤差の下限は, フィッシャー情報量を用いて, 以下のクラメール–ラオ不等ここで, 2 番目の等式では対数関数 ln が単純増加関数であることを用いた. 以上により得られた最尤推9紹介した. 数値積分で誤差 ε を達成するために, モン

式 [28] により与えられることが知られている:

ε2 = E[(â − a)2 ] ≥

1
I(a)

(

)

テカルロ法では O ε−2 回の関数の演算が必要だったのに対して, 量子コンピュータを使うことで, 演算

(51)

(

)

回数が O ε−1 回程度に減少することを示した.

なお, 最尤推定法を用いた推定では, サンプル数が大

現在, 金融を始めとして, 各分野に量子コンピュー

きくなるに連れて, 推定誤差は漸近的に, 式 (51) の下タを使う取り組みが始まっており, 今後, 将来の量子限に近づくことが知られている.

コンピュータの活用に向けて, 本稿で紹介した手法を鈴木らの手法 [21] において, 式 (46), (47) 及び (50)

基本として, 各分野に特化したような更に洗練されたからフィッシャー情報量は手法の開発が期待される.

∑1
Nk (2jk + 1)2 (52)I(a) =
S(f )(1 − S(f ))
m

付録 A

k=1

本付録では, 本文中で使用するユニタリ演算子をまとめる.
Hadamard 演算子 H の 1 量子ビット計算基底 |0⟩ への作用をと与えられる. また, 関数 f の呼出回数 N は, 以下の式により与えられる.

N=

m
∑

Nk (2jk + 1)

本文中で使用するユニタリ演算子(53)

H |0⟩ =

k=1

以上により, もし, 振幅増幅の回数 ⃗j = {j1 , · · · , jm }|0⟩ + |1⟩
√
2

(54)

と定義する.

Hadamard 演算子 H ⊗n の n 量子ビット計算基底 |0⟩n

⃗ = {N1 , · · · , Nm } を具及びベルヌーイ試行の回数 Nへの作用を体的に決めれば, 最尤推定法による推定誤差 ε の下n限と関数の呼出回数 N の関係を, 式 (51), (52) 及びH(53) を用いて導くことが出来る.
たとえば, 全ての k = 1, · · · , m に対して, ベル
ヌーイ試行の回数を一定回数 Nk = N0 として, また, jk を指数関数的に増加する場合 m1 = 0, m2 =
20 , m3 = 21 , · · · , mm = 2(m−2) を考えると, 推定
誤差 ε の下限と関数の呼出回数 N の関係は N ≥()
O ε−1 となることが示される. 最尤推定法では, サ

⊗n

2 −1
1 ∑
|x⟩n
|0⟩n = √2n x=0(55)と定義する.
1 量子ビット計算基底 |0⟩ , |1⟩ への Z の作用を,Z |0⟩ = |0⟩(56)

Z |1⟩ = − |1⟩

n ビット量子状態 |Ψ⟩n に作用する恒等演算子 In を

ンプル数が十分多いときにこの不等式の下限に漸近

In |Ψ⟩n = |Ψ⟩n

することが知られているが, 鈴木ら [21] は, 数値シ(57)と定義する.
式 (10) の確率関数 p(x) を計算するための n 量子ビットに作用する演算子をミュレーションにより, N0 = 100 回程度と選ぶことで, 下限をほぼ達成することを確かめている.
本手法は, 制御ユニタリ演算子 c Q や逆フーリエ変換演算子 Fm を使用しないため, 3.2 節において概説

P |0⟩n =した Brassard らの手法に較べて, 演算数が減少する

n
2∑
−1

√
p(x) |x⟩n

(58)

x=0

ことが期待される. 実際に, 鈴木ら [21] は, IBM の公と定義し, f (x) を計算するための n + 1 ビットに作用する演算子を開している量子コンピュータ向けの SDK(Software

Development Kit) の Qiskit[19] を用いて, 簡単な三角関数の積分について, 本手法と Brassard らの手法の量子回路を実装し, C-NOT 演算子の個数が 5% 程度に削減することを確認している.

R |x⟩n |0⟩ = |x⟩n(√f (x) |0⟩ +√1 − f (x) |1⟩)(59)と定義する.
量子逆フーリエ変換を行う n ビットに作用する演算子を

4 まとめ

1
Fn−1 |x⟩n = √2n

本稿では, 量子コンピュータの応用先の例として,
数値積分を行う 3 種類の量子アルゴリズムについてと定義する.

10

n
2∑
−1

y=0

n

e−2πixy/2 |y⟩n(60)

n 量子ビットに作用するユニタリ演算子 U が与えられた
とき, m ビットを制御ビットとして, m + n ビットに作用する制御ユニタリ演算子 c U をc

U |x⟩m |y⟩n = |x⟩m U x |y⟩nJournal of chemical physics, 149(17):174112, 2018.
[14] A. Montanaro. Quantum speedup of Monte Carlo
methods. Proceedings of the Royal Society A:
Mathematical, Physical and Engineering Sciences,471(2181):20150301, 2015.
[15] P. Rebentrost, B. Gupt and T. R. Bromley. Quantum computational ﬁnance: Monte Carlo pricing of ﬁnancial derivatives. Physical Review A,
98:022321, 2018.
[16] S. Woerner and D. J. Egger. Quantum risk analysis. npj Quantum Information, 5(1):15, 2019.
[17] A. Martin et al. Towards pricing ﬁnancial derivatives with an ibm quantum computer. arXiv
preprint arXiv:1904.05803, 2019.
[18] N. Stamatopoulos et al.
Option pricing
using quantum computers.
arXiv preprint
arXiv:1905.02666, 2019.
[19] G. Aleksandrowicz et al. Qiskit: An open-source
framework for quantum computing, 2019.
[20] G. Brassard et al. Quantum amplitude ampliﬁcation and estimation. Contemporary Mathematics,
305:53–74, 2002.
[21] Y. Suzuki et al.
Amplitude estimation
without phase estimation.
arXiv preprint
arXiv:1904.10246, 2019.
[22] D. S. Abrams and C. P. Williams. Fast quantum
algorithms for numerical integrals and stochastic
processes. arXiv:quant-ph/9908083, 1999.
[23] A. Papageorgiou and J. F. Traub. Quantum Algorithms for Continuous Problems and Their Applications, pages 151–178. John Wiley & Sons, Ltd,
2014.
[24] M. A. Nielsen and I. Chuang. Quantum computation and quantum information, 2002.
[25] 石坂 智, 小川 朋宏, 河内 亮周, 木村 元, 林 正人. 量子
情報科学入門. 共立出版, 2012.
[26] L. K. Grover. A fast quantum mechanical algorithm for database search. Proceedings of 28th
Annual ACM Symposium on Theory of Computing, pages 212–219, 1996.
[27] A. Y. Kitaev. Quantum measurements and the
Abelian stabilizer problem. Electronic Colloquium on Computational Complexity, 3, 1996.
[28] C. R. Rao. Linear statistical inference and its applications, volume 2. Wiley New York, 1973.

(61)と定義する.

引用文献[1] QUANTUM COMPUTING REPORT.
https://quantumcomputingreport.com/
scorecards/qubit-count. Accessed:2019-05-29.
[2] IBM Q Experience.
https://quantumexperience.ng.bluemix.net/
qx/editor. Accessed:2019-05-29.
[3] Quantum Cloud Services - Rigetti.
https://www.rigetti.com/qcs. Accessed:201905-29.
[4] 慶應義塾大学.
https://www.keio.ac.jp/ja/news/2018/5/22/
27-44149/. Accessed:2019-05-29.
[5] P. Glasserman. Monte Carlo Methods in Financial
Engineering. Springer, 2004.
[6] 湯前 祥二, 鈴木 輝好. モンテカルロ法の金融工学への応用. 朝倉書店, 2000.
[7] W. Krauth. Statistical Mechanics: Algorithms
and Computations. Oxford University Press,
2006.
[8] M. Creutz, L. Jacobs and C. Rebbi. Monte Carlocomputations in lattice gauge theories. Physics
Reports, 95(4):201–282, 1983.
[9] 青木 慎也. 格子上の場の理論. シュプリンガー・ジャ
パン, 2012.
[10] M. EricksonKirk et al. Technical basis for revisionof the pressurized thermal shock (PTS) screening
limit in the PTS rule (10 CFR 50.61). NUREG1806, US Nuclear Regulatory Commission, 2007.
[11] J. Katsuyama et al. Guideline on a structural
integrity assessment for reactor pressure vessel
based on probabilistic fracture mechanics. JAEAResearch, 2017.
[12] C. M. Johnson et al.
Monte Carlo explicitly correlated second-order many-body perturbation theory. The Journal of chemical physics,
145(15):154115, 2016.
[13] C. M. Johnson et al. Monte Carlo explicitly correlated many-body Green’s function theory. The

11