金融工学の計算に必要な量子リソースの具体的な見積もり — Black–Scholes方程式をLCHS法で解く場合

量子コンピュータの応用先の一つとして、金融工学（Quantitative Finance）が挙げられる。金融工学ではその対象に応じて様々な数理的枠組みが現れるが、そのひとつが偏微分方程式である。偏微分方程式は流体力学や構造解析などCAE（Computer-Aided Engineering）と呼ばれる分野にもしばしば登場し、量子コンピュータによる求解や高速化について多様な手法が知られている。本技術ブログでは、金融工学の典型的な問題である「オプション価格の評価」に用いられる偏微分方程式であるBlack–Scholes方程式に関して、QURI SDKを用いた量子回路実装や必要リソース見積もりを紹介する。具体的には、Black–Scholes方程式を差分法により離散化し、LCHS（Linear Combination of Hamiltonian Simulation）と呼ばれる手法を用いて方程式の解を計算する量子回路を設計・シミュレータ上に実装した。さらに、必要な計算リソースの見積もりも行った。

※ なお、紙幅の都合で詳細な技術的内容は省略している。興味を持たれた方は参考文献等を参照してほしい。

金融工学

金融工学とは、金融商品の価格評価・リスク管理・投資判断などを数学的に扱う分野である。その対象としては以下のような問題が挙げられる。

資産価格の変動をどうモデル化するか — 株価・金利・為替などが時間とともにどう動くかを数式で表す問題
金融派生商品を「公正な価格」でどう評価するか — 原資産（株式など）の価格からそれに依存する契約の価格を決定する問題
リスクをどう管理・最適化するか — 保有資産の組み合わせ（ポートフォリオ）をどう選ぶか、損失リスクをどう測るかという問題

これらの問題に対応して、それぞれ確率微分方程式（資産価格モデル）、偏微分方程式（価格付け方程式）、最適化問題（ポートフォリオ最適化）といった異なる種類の数理的枠組みが用いられる。

本記事で扱うオプションは、将来の決められた時点（満期 $T$ ）に、原資産（株式など）をあらかじめ決めた価格（ストライク $K$ ）で売買する「権利」を持つ金融派生商品の一種である。権利行使のタイミングによって、満期のみで行使できるヨーロピアン型や、満期までの任意の時点で行使できるアメリカン型などに分類される。

以下ではヨーロピアンコールオプション（買う権利）について考える。満期時点の原資産価格を $S_T$ とすると、 $S_T>K$ であれば権利を行使して $S_T-K$ の利益を得て、 $S_T\le K$ であれば権利を放棄すればよい。よって、満期 $T$ での価値は $\max(S_T-K,0)$ と表される。このオプションの現在時刻における価値を評価することが目標となる。

Black–Scholes方程式

最も単純な資産価格モデルとして、時刻 $t$ での原資産価格 $S_t$ が幾何ブラウン運動

dS_t = \mu S_t\, dt + \sigma S_t\, dW_t

（ $W_t$ は標準ブラウン運動、 $\mu$ はドリフト、 $\sigma$ はボラティリティと呼ばれる定数）に従うとする。市場が無摩擦であり取引コストがかからないことや、リスクなし・元手なしで確実に利益を得られる取引機会は存在しないという無裁定条件を仮定すると、時刻 $t$ で原資産価格が $S$ であった時のオプション価格 $V(S,t)$ が満たす偏微分方程式：

\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2S^2\frac{\partial^2V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0

（ $r$ は無リスク金利）を得る。また、終端条件は $V(S,T)=\max(S-K,0)$ と表される。これが1資産のBlack–Scholes方程式であり、複数の原資産を考慮することで多資産へと拡張することも可能である。

1資産の場合には解析解

V(S,t) = S\,\Phi(d_1) - Ke^{-r(T-t)}\,\Phi(d_2) \\ d_1 = \frac{\ln(S/K)+\left(r+\frac{1}{2}\sigma^2\right)(T-t)}{\sigma\sqrt{T-t}}\, \\ d_2 = d_1-\sigma\sqrt{T-t}

（ $\Phi$ は標準正規分布の累積分布関数）が知られている。パラメータを $K=60,\sigma=0.2,r=0.03,T=3$ として、満期までの残り時間ごとの解析解 $V(S,t)$ を図1に示す。満期 $T$ に近づくにつれて、解析解は満期のペイオフ $\max(S-K,0)$ （図中破線）に近づいていく様子が見て取れる。なお、多資産の場合には一般に解析解は存在せず、数値計算を行う必要がある。

図1: 1資産Black–Scholes方程式の解析解

境界条件

解を定めるには $S\to 0,\infty$ の両端で境界条件を与える必要がある。 $S=0$ ではオプションの価値は恒等的に $0$ になる（ $V(0,t)=0$ ）。一方、 $S\to\infty$ ではコールオプションはほぼ確実に権利行使される契約と同等になる。（連続複利である）金利 $r$ の下では、満期 $T$ での価値 $K$ は現在（時刻 $t$ ）の $Ke^{-r(T-t)}$ に相当するため、価格は $V(S,t) = S-Ke^{-r(T-t)}$ となり、特に $S$ 微分は $\partial_S V(S,t)= 1$ に漸近する。これらは解析解の挙動とも整合している。

量子回路実装の方針

Black–Scholes方程式のような線形偏微分方程式は、空間方向（この場合は原資産価格 $S$ ）を $N$ 点に離散化することで、時刻 $t$ での各点での値を並べた $N$ 次元ベクトル $\vec{V}(t)$ に対するアフィン写像

\frac{d\vec{V}(t)}{dt} = A \vec{V}(t) + \vec{b}

（ $A$ は $N$ 次正方行列、 $\vec{b}$ は $N$ 次元ベクトル）として書き直すことができる。満期 $T$ での価格（ペイオフ）は既知であり、知りたいのは現在時刻での価格 $\vec{V}(0)$ であったため、形式的には

\vec{V}(0) = e^{-AT} \vec{V}(T) - \int_0^T d\tau\, e^{-A\tau} \vec{b}

を計算すればよい。この形は流体・構造解析などCAEで現れる線形偏微分方程式と全く同じであり、Schrödingerisationなど既存の量子アルゴリズムがそのまま応用できる[1]。本検討では、非ユニタリな時間発展をユニタリ操作の重み付き和として実装するLCHSと呼ばれる手法を採用した。

LCHSによる量子回路実装

空間の離散化

十分大きな上限 $S_{\max}$ に対し、原資産価格 $S\in[0,S_{\max}]$ を幅 $\Delta S$ で $N$ 点（ $S_0, \dots, S_{N-1}$ ）に離散化し、 $V_j(t):=V(S_j,t)$ とする。Black–Scholes方程式は中心差分近似により

\frac{dV_j}{dt} = rV_j - \frac{1}{2}\sigma^2 S_j^2 \frac{V_{j+1} - 2V_j + V_{j-1}}{(\Delta S)^2} - r S_j \frac{V_{j+1} - V_{j-1}}{2\Delta S}

という常微分方程式系となる（ $j=1,\dots,N-1$ . 引数 $t$ は省略した）。

なお、端点については、 $j=0$ では境界条件 $V(0,t)=0$ より常に $V_0=0$ とすれば良い。 $j=N-1$ では仮想的な格子点 $V_N$ が出現するが、境界条件 $\lim_{S\to\infty}\partial_S V(S,t)=1$ より $V_N=V_{N-2}+2\Delta S$ と定めれば良い。

各格子点での価格を並べたベクトルを

\vec{V}(t) = (V_0(t),\dots,V_{N-1}(t))

とし、この方程式を整理すれば行列 $A$ 、ベクトル $\vec{b}$ が得られる。また、終端（時刻 $T$ ）での値は $V_j(T)=\max(S_j-K,0)$ である。

LCHS

行列 $A$ は一般に（反）エルミートではなく、上記の時間発展 $e^{-At}$ は非ユニタリになる。一方で、量子コンピュータで自然に実装・実現できるのはユニタリ変換であり、非ユニタリ時間発展を計算するには工夫が必要となる。

LCHS（Linear Combination of Hamiltonian Simulation）は、非ユニタリな時間発展 $e^{-At}$ を、ユニタリな時間発展 $e^{-i(kL+H)t}$ の重み付き積分

e^{-At} = \frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}} dk\; f(k)\, e^{-i(kL+H)t} \approx \frac{1}{\sqrt{2\pi}}\sum_{m=-M}^{M} \Delta k\, f(k_m)\, e^{-i(k_mL+H)t} \tag{1}

として表し、量子回路上に実装する手法である。ここで、

$f(k)$ はカーネル関数と呼ばれる。
分割数は $2M+1$ 、幅は $\Delta k$ としており、各分割の代表点は $k_m=m\Delta k$ である。
$L=(A+A^\dagger)/2$ 、 $H=(A-A^\dagger)/(2i)$ はいずれもエルミート行列である。

そして、各 $e^{-i(kL+H)t}$ はQSVT（Quantum Singular Value Transformation; 量子特異値変換）により、 $k$ についての和はLCU（Linear Combination of Unitaries）により実現することができる。（LCUとQSVTについては、Quantum Native Dojo!の書籍や補足ノートも参照）。QSVTでは、ユニタリとは限らない任意の行列を部分行列として埋め込んだ量子回路を用いて、行列の固有値・特異値をさまざまに多項式変換したものを埋め込んだ新たな量子回路を実現することができる。ここでは $e^{-ixt}$ を $x$ の有限次多項式で近似することで、ユニタリな時間発展 $e^{-i(kL+H)t}$ を実現している[2]。カーネル関数 $f(k)$ にはいくつかの数学的な制約があり、ここでは Low-Sommaによって提案された関数

f(k;\gamma)= \sqrt{\frac{2}{\pi}} \frac{e^{1-ik}}{1+k^2} \exp\left(-\frac{k^2+1}{4\gamma^2}\right)

を用いた[3]。なお、 $\gamma$ は許容誤差に応じた実数である。

このLCHSを用いたシミュレーションにおいては、上記(1)式で積分区間を有限にして分割したことによる誤差(LCHSそのものの誤差)や、QSVTで $e^{-ixt}$ を $x$ の有限次多項式として近似したことによる誤差がそれぞれ存在する。実際のシミュレーションにおいては、解全体の許容誤差に応じた適切な分割や次数を選択する必要がある。

また、LCHSではエルミート行列 $L,H$ を量子回路として実現する必要がある。例えば今回の $L$ の対角成分 $(j,j)$ は $r+\sigma^2j^2$ のように $j$ の多項式として与えられており、本検討では、この多項式を古典的に計算して対応する成分を多重制御ゲートで実装する方式（Encoding A）と、量子回路上で多項式を計算して実装する方式（Encoding B）の2通りを実装し、比較した。

スケーリングと量子優位性

格子点数を $N=2^n$ とし、格子点数に依存するコストを考える。詳細は省略するが、各行列 $L, H$ を量子回路に埋め込むために必要なゲート数は、Encoding A では $O(2^n)$ 、Encoding B では $O(n)$ となる（代償として、Encoding Bでは補助量子ビットがより多く必要になる）。以下では、特に指定しない限りゲート数効率の良い Encoding B を仮定する。2階微分項を含むため行列 $A$ のノルムは $\|A\|\sim O(1/(\Delta S)^2)\sim O(2^{2n})$ となり、それに応じてQSVTの呼び出し回数が増えることから、時間発展計算に要するゲート数は $O(2^{2n}n)$ となる。

実際の量子計算においては、ここまで見てきた時間発展部分の実装に加え、初期量子状態の準備と、計算結果の読み出しを行う必要がある。

初期状態準備は、 $t＝T$ における終端条件 $\vec{V}(T)$ に対応する量子状態として、単純な区分線形関数（ $\max(S-K,0)$ ）を状態の振幅に埋め込んだ状態を用意する。これは $O(n)$ のコストで実現される。
LCHSを用いた時間発展シミュレーションによって得られた量子状態（終状態）に測定を行うことで、 $t=0$ におけるオプション価格を読み出すことができる。終状態の各計算基底の振幅の典型的な大きさは $O(1/\sqrt{N})$ であるため、特定の $V_j(0)$ を得るための量子振幅推定（Quantum Native Dojo! 書籍8.3節参照）には $O(2^{n/2})$ 回の演算が必要であると見積もられる。

これらを合わせれば、量子計算において、格子点数に依存するコストは $(O(n)+O(2^{2n}n)) \cdot O(2^{n/2}) = O(2^{5n/2}n)$ となる。

一方、古典計算においては、一般に行列の指数関数 $e^{-AT}$ の計算量は固有値分解などを用いることで $O(2^{3n})$ であることが知られており、量子計算に優位性があるようにも見える。しかし、今回のように疎行列と1つの初期ベクトルに対して $e^{-AT}V$ を求めれば良い場合は、Krylov 部分空間法を用いることが可能で、計算量は $O(2^{2n})$ が支配的であると予想されるから、優位性について確定的なことは言いにくいのが実情である。

このように、漸近的なオーダー解析だけでは、実際の問題に対して量子優位性が得られるかどうかの判断は難しい。量子コンピュータを活用した計算の可能性を信頼性高く見極めるためには、具体的な問題サイズ・誤差に対して必要なゲート数や量子ビット数をきちんと見積もる必要がある。

QURI SDK [4]では、そのような具体的なリソース見積もりを簡単に行うことができる。そこで本検討では、QURI SDKを用いて実際にLCHSの量子回路の実装（設計）を行い、それを誤り訂正量子コンピュータの標準的なゲートセット（Cliffordゲート、 $T$ ゲート、1量子ビットパウリ回転ゲート、トフォリゲート）にコンパイルした場合のリソースについて、定量的な評価を行った。最終的には、これらのゲートをCliffordゲートと1量子ビットパウリ回転ゲートにまで分解し、後者を「パウリ回転ゲート数」としてプロットした。誤り耐性量子計算においては非Cliffordゲートが計算コストの大宗を占めると考えられており、この数値が実際のコストの良い指標となると考えられる。

実装結果

1資産Black–Scholes方程式

パラメータを $S_{\max}=120,K=60,\sigma=0.05,r=0.03,N=2^4,T=1$ とし、LCHSの許容誤差 $\epsilon=0.1$ 、QSVTの許容誤差 $\epsilon_t=0.1$ として計算した結果を図2に示す。左図において、LCHSは量子回路（の古典ベクトルシミュレーション）での計算結果、Classical LCHS は積分の分割などは変えずに各 $e^{-i(k_mL+H)t}$ を古典的に計算した結果、Matrix exponential は $e^{-At}$ を直接古典的に計算した結果、Analyticalは解析解である。また、右図における各工程の誤差はそれらの差を順に取ることで求められる。解析解とよく一致するオプション価格が得られており、誤差もそれぞれ許容値以下となっていることが分かる。

図2: 1資産Black–Scholes方程式の計算結果（左）と誤差（右）

次に、格子点数 $N=2^n$ に対して必要なパウリ回転ゲート数および補助量子ビットを含めた論理qubit数を計測した結果を図3に示す。既に述べたように、理論的なゲート数のオーダーは $O(2^{2n}n)$ である（視覚的に把握しやすいよう図中では $10\cdot 2^{2n}n$ をプロットしている）。Encoding B では Encoding A と比べて qubit数は増えるが、ゲート数は確かに理論的なオーダーに漸近していると言える。

図3: パウリ回転ゲート数（左）とコンパイル後の論理qubit数（右）

2資産Black–Scholes方程式

同様の枠組みは多資産の場合にも拡張できる。2資産Black–Scholes方程式は

\frac{\partial V}{\partial t} + \frac{1}{2} \sigma_1^2 S_1^2 \frac{\partial^2 V}{\partial S_1^2} + \rho \sigma_1 \sigma_2 S_1 S_2 \frac{\partial^2 V}{\partial S_1 \partial S_2} + \frac{1}{2} \sigma_2^2 S_2^2 \frac{\partial^2 V}{\partial S_2^2} + r S_1 \frac{\partial V}{\partial S_1} + r S_2 \frac{\partial V}{\partial S_2} - r V = 0

（ $\sigma_i$ は資産 $i$ のボラティリティ、 $r$ は無リスク金利、 $\rho$ は2資産の相関係数）である。終端条件としては Worst-of Call

V(S_1,S_2,T)=\max(0,\min(S_1,S_2)-K)

を採用する。対応する境界条件は

V(0,S_2,t) = V(S_1,0,t) = 0 \,,\quad \lim_{S_1\to\infty} \partial_{S_1} V(S_1,S_2,t) = 0 \,,\quad \lim_{S_2\to\infty} \partial_{S_2} V(S_1,S_2,t) = 0

となる。

1資産と同様に $S_1\in[0,S_{1,\max}], S_2\in[0,S_{2,\max}]$ をそれぞれ $N_1, N_2$ 点に離散化すると、価格 $V$ は $N_1N_2$ 次元ベクトルとして表される。中心差分近似を取り係数を整理すれば $d\vec{V}/dt=A\vec{V}$ という斉次形になり、 $\vec{V}(0)=e^{-AT}\vec{V}(T)$ を計算すればよい。

パラメータを $S_{1,\max}=S_{2,\max}=120,K=40,\sigma_1=0.05,\sigma_2=0.02,r=0.03,\rho=0.4,N_1=N_2=2^2,T=1$ 、LCHSの許容誤差 $\epsilon=0.1$ 、QSVTの許容誤差 $\epsilon_t=0.1$ として計算した結果を図4に示す。原資産価格が小さい方に価値が依存する Worst-of call の特徴が出ている。また、右図では格子点数 $N_1=2^{n_1}, N_2=2^{n_2}$ に対して必要なパウリ回転ゲート数（Encoding B）及び $\alpha=\sigma_1^2N_1^2+\sigma_2^2N_2^2$ に対するフィッティングを示した。ゲート数の支配的なオーダーは $\alpha^{0.9}$ 程度となっており、これは既知の理論的なオーダー $O(N_1^2)+O(N_2^2)$ と整合する。

図4: 2資産Black–Scholes方程式の計算結果（左）とパウリ回転ゲート数（右）

まとめと今後の展望

本検討では、1資産および2資産のBlack–Scholes方程式の解を求める問題について、LCHSを用いた量子回路実装とシミュレーションを行った。また、求解の際の格子点数を変化させ、必要な計算リソース（非Cliffordゲート数・論理qubit数）を見積もった。

今後の課題としては、金利やボラティリティが時間に依存するようなより現実的な市場を反映したモデルへの拡張が考えられる。

また、本検討で見積もったリソースは、誤り訂正を考慮しない理想的な量子計算機を仮定したものである。実際の誤り耐性量子計算機（FTQC）上での実装を見据えると、surface codeなどの量子誤り訂正符号による符号化を踏まえた物理qubit数・ゲート数・実行時間の見積もりが必要となる。QURI SDKにはこうしたsurface codeベースのアーキテクチャへのコンパイルやリソース見積もりを行うQURI VMという機能があり、これを用いてハードウェアにより近い形でのリソース評価を行うことも本検討の重要なfuture workである。

参考文献

[1] N. Guseynov, N. Liu, C. S. Pun, T. Vaidya, “End-to-End PDE-Based Quantum Algorithms for Multi-Asset Option Pricing under Local and Stochastic Volatility,” arXiv:2605.26610 (2026).

[2] J. M. Martyn, Z. M. Rossi, A. K. Tan, I. L. Chuang, “Grand Unification of Quantum Algorithms,” PRX Quantum 2, 040203 (2021)

[3] G. H. Low, R. D. Somma, “Optimal quantum simulation of linear non-unitary dynamics,” arXiv:2508.19238 (2025).

[4] QunaSys Inc., https://quri-sdk.qunasys.com