量子アルゴリズム Grover Adaptive Search でポートフォリオ最適化問題を計算する方法の紹介

はじめに

本記事では論文^[1]で提案されている Grover Adaptive Search の QUBO (Quadratic Unconstrained Binary Optimization) 問題への応用を使ってポートフォリオ最適化問題を解く方法を紹介する. QUBO 問題は 量子アニーリング方式を使った計算が有名であるが, 本記事の方法は量子ゲート方式のアルゴリズムである. Grover Adaptive Search は Qiskit Optimization の Grover Optimizer^[2] での計算実装を紹介する. 量子ゲート方式の実機を使っての計算は試していない.

対象読者

• 量子ゲート方式の量子アルゴリズムで QUBO 問題を解く方法に興味のある方.
• ポートフォリオ最適化問題を量子アルゴリズムで解く方法に興味のある方.
• Grover Search アルゴリズムの応用に興味のある方.

Grover Adaptive Search

計算基底の表記

本記事では $n$ 量子ビットを使った計算を考え, 対応する計算基底 ($N = 2^n$ 個のベクトルからなる) を以下のように表す.

$$\{ \ket{0}_n, \ket{1}_n, \dots, \ket{2^n -1}_n \}.$$

以下, 量子ゲートを構成するユニタリ作用素はこの計算基底によって張られる複素ベクトル空間 $C^N$ 上で考える.

Grover Search に使われる作用素

Grover Adaptive Search の基礎となる Grover Search アルゴリズムは, 以下の $3$ つのユニタリ作用素を繰り返し適用して目的の番号列 $I \subset \{ 0, \dots, 2^n -1 \}$ を高確率で抽出できる量子状態を得るアルゴリズムであった.

1. $A$: 量子状態を用意する作用素.
   Grover Search では $A$ として $n$ 個の Hadamard ゲートのテンソル積 $H^{\otimes n}$ を考え, 未知のどんな探索対象番号列にもバランスよく対応できるようにするため, アルゴリズムの最初にすべての基底ベクトルの状態の等しい重ね合わせ状態を用意する.

$$A\ket{0}_n = H^{\otimes n}\ket{0}_n = \frac{1}{\sqrt{2^n}}\sum_{i=0}^{2^n -1}\ket{i}_n.$$

2. $O$: オラクル作用素.
   探索対象番号列 $I \subset \{ 0, \dots, 2^n -1 \}$ に対してオラクル作用素 $O$ は $I$ に属する番号の基底ベクトルは $-1$ 倍し, 属さない基底ベクトルは変化させないユニタリ作用素である. オラクル作用素をすべての基底ベクトルの状態の等しい重ね合わせ状態に適用すると以下のようになる.

$$O (A\ket{0}_n) = \frac{1}{\sqrt{2^n}}\sum_{i \notin I}\ket{i}_n - \frac{1}{\sqrt{2^n}}\sum_{i \in I}\ket{i}_n.$$

3. $D$: Grover 拡散作用素は $\ket{0}_n$ を変化させず, 他の基底ベクトルを $-1$ 倍するユニタリ作用素である.

Grover アルゴリズムについては以下の記事で簡単に述べた.

https://zenn.dev/quetta_kazunoko/articles/grover_search_differential_geometric_strategy

この記事の鏡映変換 $I_r$ は上記の状態 $\frac{1}{\sqrt{2^n}}\sum_{i \in I}\ket{i}_n$ についての鏡映変換であり, オラクル作用素 $O$ を作用させる変換に対応している. 鏡映変換 $I_a$ は上記の状態 $\frac{1}{\sqrt{2^n}}\sum_{i=0}^{2^n -1}\ket{i}_n$ についての鏡映変換であり, 変換 $ADA^{\dagger}$ に対応している. Grover iteration $G$ は以下のように計算される.

$$G = ADA^{\dagger}O. \tag{1}$$

$A^{\dagger}$ はユニタリ作用素 $A$ のエルミート共役である.

Grover Adaptive Search アルゴリズム

Grover Adaptive Search は Grover Search アルゴリズムを $2$ 値変数 $x \in X=\{ 0,1 \}^n$ の関数 $f:X \longrightarrow \mathbb{R}$ の最小値を見つける問題

$$\min_{x \in X}f(x)$$

に応用したアルゴリズムである^[1:1]. アルゴリズムは以下のような手順で計算される.

1. 閾値 $y \in \mathbb{R}$ に対して条件 $f(x) < y$ を満たす $x \in X$ を探索するユニタリ作用素 $A_y$, $O_y$ を選ぶ.
2. $A_y$, $O_y$ を使った Grover Search を実行して $f(\tilde{x}) < y$ を満たす $\tilde{x} \in X$ を得る.
3. 上記 2. の解 $\tilde{x} \in X$ に対して $y=f(\tilde{x})$ を新しい閾値として 1. の計算から繰り返す.
4. あらかじめ決めておいた上記工程 1, 2 の最大繰り返し回数あるいは計算終了閾値 $y$ に達したところで計算を終了する.

最初の閾値 $y$ を選ぶ方法や工程 2. の Grover iteration の回数などの Grover Adaptive Search アルゴリズムの詳細は論文^[1:2]を参照.

QUBO 問題への応用

QUBO 問題

$2$ 値 $n$ 変数の QUBO 問題は行列 $Q \in \mathbb{R}^{n\times n}$, 定数ベクトル $b \in \mathbb{R}^n$, 定数 $c \in \mathbb{R}$ に対して

$$\min_{x_0,..., x_{n-1} \in \{ 0,1 \}} \left( \sum_{i,j=0}^{n-1}Q_{i,j}x_i x_j + \sum_{i=0}^{n-1}b_i x_i + c \right) \tag{2}$$

と定義される最適化問題である. 上記 Grover Adaptive Search で $f(x) = x^{T}Qx + b^{T}x + c$ としたものである.

本記事では Grover Adaptive Search の QUBO 問題への応用しか紹介しないが, 論文^[1:3]にはより高次の多項式からなる目的関数の問題への応用や CPBO (Constrained Polynomial Binary Optimization) 問題への応用も提案されている.

ユニタリ作用素 $A$ の構成方法

単項式の値の重ね合わせ状態

QUBO 問題を解くために, 以下の $2$ 種類のレジスターを用意する.

• 目的関数の変数の数 $n$ のビット数をもつ入力レジスター $\ket{x}_n$.
• 閾値 $y$ の判定の精度を保証する桁数 $m$ のビット数をもつ出力レジスター $\ket{z}_m$.

通常の Grover Search と同様に, 量子探索の初期状態としてレジスター $\ket{x}_n$, $\ket{z}_m$ はすべての基底ベクトルの重ね合わせ状態をとる.

$$\begin{aligned} \ket{x}_n &= H^{\otimes n}\ket{0}_n, \\ \ket{z}_m &= H^{\otimes m}\ket{0}_m. \tag{3} \end{aligned}$$

目的関数 $f(x) = x^{T}Qx + b^{T}x + c$ を単項式の和として以下のように表現する.

$$P(x) = \sum_{J \subset \{ 0,..., n-1 \}} a_J \prod_{j \in J} x_j. \tag{4}$$

定数項は $a_{\empty}$ と表現している. 後で述べる出力レジスターのエンコードを行うユニタリ作用素 $U_G$ を使って, 多項式 $P(x)$ 中に現れる単項式 $a_J \prod_{j \in J} x_j$ 全てに対して以下の制御ユニタリ作用素を計算する.

$$C^{J}\left(U_{G}\left(\frac{2\pi}{2^m} a_J \right)\right) \ket{x}_n \ket{z}_m = \ket{x}_n U_{G}\left(\frac{2\pi}{2^m} a_J \right)^{\displaystyle \prod_{j \in J} x_j}\ket{z}_m. \tag{5}$$

ユニタリ変換 $U_{G}\left(\frac{2\pi}{2^m} a_J \right)^{\prod_{j \in J} x_j}$ は入力レジスター $\ket{x}_n$ の $J$ に含まれる番号 $j$ の量子ビットたちを制御ビットとして, 出力レジスター $\ket{z}_m$ にユニタリ作用素 $U_{G}\left(\frac{2\pi}{2^m} a_J \right)$ を実行する作用素である.

最適化問題の入力変数 $x$ を $2$ 値変数としたのは, この単項式の値の重ね合わせ状態を用意するときに $x=(x_0, \dots, x_{n-1})$ を制御ビットとして使うためである.

出力レジスターのエンコード

上記の単項式の重ね合わせ状態を用意するため, 出力レジスター $\ket{z}_m$ のエンコードを行う. 角度 $\theta \in [-\pi, \pi)$ に対して, 以下のユニタリ作用素を考える.

$$U_G(\theta): H^{\otimes m}\ket{0}_m \longmapsto \frac{1}{\sqrt{2^m}}\sum_{l=0}^{2^m -1} e^{il\theta}\ket{l}_m. \tag{6}$$

このユニタリ作用素は各番号 $i=0, \dots, m-1$ (虚数ではなく自然数) に対して $m-1-i$ 番目の量子ビットに位相ゲート $R(2^{i}\theta)$ を適用することで得られる^[1:4] ($\theta$ の各位相成分への分解の重ね合わせ状態). 整数 $k$ ($-2^{m-1} \leq k < 2^{m-1}$) に対して, $\theta=2\pi k / 2^m$ ととれば,

$$QFT^{\dagger} \circ U_G \left(\frac{2\pi}{2^m} k \right)(H^{\otimes m}\ket{0}_m) =\ket{k \:\text{(mod } 2^m \text{)}}_m \tag{7}$$

となる. つまり $U_G \left(\frac{2\pi}{2^m} k \right)$ は整数 $k$ を周期 $2^m$ の信号としてエンコードし, 逆量子フーリエ変換 $QFT^{\dagger}$ によってデコードすると, mod $2^m$ の非負整数になる. この表現を $k$ の binary Two's Complement という^[1:5].

• $U_G \left(\frac{2\pi}{2^m} k \right)$ によってエンコードされた信号は式 $(5)$ で単項式の重ね合わせ状態を用意するために使われる.
• 逆量子フーリエ変換 $QFT^{\dagger}$ は多項式 $P(x)$ 中の単項式のすべての重ね合わせ状態に対して適用し (式 $(8)$), そのデコード結果は Grover iteration でオラクル作用素 $O$ が出力レジスター $\ket{z}_m$ の値の正負を判定するのに使われる.
• $k$ が負整数のときは逆量子フーリエ変換 $QFT^{\dagger}$ でデコードされたときに, $2^{m-1}$ より大きな番号に対応する基底ベクトル $\ket{2^m + k}_m$ となる.

$k$ が整数でない時の扱いは論文^[1:6]の Appendix B を参照.

ユニタリ作用素 $A$ の定義

以上のように多項式 $P(x)$ に含まれる単項式の分だけ出力レジスター $\ket{z}_m$ に制御ユニタリ作用素 $U_{G}\left(\frac{2\pi}{2^m} a_J \right)^{\prod_{j \in J} x_j}$ を適用することで, 量子状態を用意するユニタリ作用素 $A$ が以下のように定義できる.

$$A = QFT^{\dagger} \left(\prod_{J \subset \{ 0,...,n-1 \}} C^{J}\left(U_{G}\left(\frac{2\pi}{2^m} a_J \right) \right) \right) (H^{\otimes n}\otimes H^{\otimes m}). \tag{8}$$

先頭の逆量子フーリエ変換は最後に出力レジスター $\ket{z}_m$ をデコードする作用素である. ユニタリ作用素 $A$ の量子回路図や量子ゲート数については論文^[1:7] を参照. Grover Adaptive Search の閾値 $y \in \mathbb{R}$ に対して多項式 $P(x) - y$ を考えることで, 各レジスターでの計算を以下のように表せる.

$$A_y \ket{0}_n \ket{0}_m = \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n -1}\ket{x}_n \ket{P(x) -y}_m. \tag{9}$$

ユニタリ作用素 $O$ の定義

オラクル作用素 $O$ は出力レジスター $\ket{P(x) - y}_m$ の正負を式 $(7)$ に現れる $\ket{k \:\text{(mod } 2^m \text{)}}_m$ として判定する. Grover Adaptive Search の一般的なアルゴリズムの定義としてはオラクル作用素は $O_y$ と閾値 $y$ に応じて決まる書き方をしたが, QUBO への応用をする上では $y$ に関係なく出力レジスターが表す数の符号を判定するだけである.

$$O \ket{x}_n \ket{z}_m = \text{sign}(z) \ket{x}_n \ket{z}_m. \tag{10}$$

$\ket{k \:\text{(mod } 2^m \text{)}}_m$ の正負は $k$ が $2^{m-1}$ より大きいかどうかで判定するので, 出力レジスター $\ket{z}_m$ の最大桁のビットによって判定される.

QUBO 形式のポートフォリオ最適化問題

QUBO 問題の応用として $2$ 値変数のポートフォリオ最適化問題^[3]がある.

$2$ 値変数のポートフォリオ最適化問題

• $n$: 投資対象銘柄とする資産 (asset) の数を入力レジスターの量子ビット数にとる.
• $\mu \in \mathbb{R}^n$: ポートフォリオで期待されるリターン (期待値ベクトル).
• $\Sigma \in \mathbb{R}^{n\times n}$: ポートフォリオのリスク (共分散行列).
• $q \geq 0$: リスクファクター.

このとき, 以下の最適化問題が $2$ 値変数のポートフォリオ最適化問題 (論文^[1:8] の $\sect$ 4.1) である.

$$\min_{x \in \{ 0, 1 \}^n} (q x^{T} \Sigma x - \mu^{T} x). \tag{11}$$

• 上記のポートフォリオ最適化問題は $2$ 値変数の最適化問題であるが, 記事^[4]のように QUBO 形式のポートフォリオ最適化問題を連続変数の最適化に利用する方法もある.
• 連続変数のポートフォリオ最適化問題を量子アルゴリズムで直接解く方法として HHL を使った方法^[5]がある.

Grover Optimizer での実装

論文^[1:9] $\sect$ 4.1.1 の例で Grover Optimizer^[2:1] を使った計算を紹介する. 資産の数 (入力レジスター量子ビット数) を $n=3$, リスクファクターを $q = 0.5$, 期待値 $\mu$ と共分散 $\Sigma$ を,

$$\mu = \begin{pmatrix} 1 \\ -2 \\ 3 \end{pmatrix}, \Sigma = \begin{pmatrix} 2 & 0 & -4 \\ 0 & 4 & -2 \\ -4 & -2 & 10 \end{pmatrix},$$

とすると, 最適化問題 (式 $(11)$) を以下のように書き直せる.

$$\min_{x_1, x_2, x_3 \in \{ 0, 1 \}} (-2x_1 x_3 - x_2 x_3 -x_1 + 2x_2 -3x_3).$$

以下はチュートリアルを参考にした実装. まず行列形式で目的関数を指定し, その目的関数を多項式で出力する.

import numpy as np
from qiskit_optimization.algorithms import GroverOptimizer
from qiskit_optimization.translators import from_docplex_mp
from qiskit import BasicAer
from docplex.mp.model import Model

backend = BasicAer.get_backend("statevector_simulator")

model = Model()
x1 = model.binary_var(name="x1")
x2 = model.binary_var(name="x2")
x3 = model.binary_var(name="x3")

# 行列形式で目的関数の係数を指定する
sigma = np.array([[2.0,   0.0,  -4.0],
                  [0.0,   4.0,  -2.0],
                  [-4.0, -2.0,  10.0]])
q = 0.5
sigma = sigma - np.tril(sigma) # 非対角係数だけを残す
sigma *= q
mu = np.array([1, -2, 3])

qp = from_docplex_mp(model)
qp.minimize(constant = 0, linear = -mu, quadratic = sigma)
print(qp.prettyprint())

以下のような出力となる.

Problem name: docplex_model1

Minimize
  -2*x1*x3 - x2*x3 - x1 + 2*x2 - 3*x3

Subject to
  No constraints

  Binary variables (3)
    x1 x2 x3

出力レジスターの量子ビット数を $6$, Grover Adaptive Search の繰り返し回数を $10$ 回で打ち切るという設定で計算を実行する.

grover_optimizer = GroverOptimizer(6, num_iterations=10, quantum_instance=backend)
results = grover_optimizer.solve(qp)
print(results.prettyprint())

以下のような出力となる.

objective function value: -6.0
variable values: x1=1.0, x2=0.0, x3=1.0
status: SUCCESS

• このコードを PC で実行しても高速計算はできない.
• 量子ゲートマシン実機での計算については論文^[1:10] $\sect$ 4.2 を参照.

まとめ

• Grover Adaptive Search は Grover Search アルゴリズムを $2$ 値変数関数の最適化問題に応用したアルゴリズムである.
• $2$ 値変数関数の最適化問題の例として QUBO 問題がある.
• QUBO 問題に Grover Search アルゴリズムを適用するために, 目的関数の多項式の値を量子状態の重ね合わせとして表現した.
• QUBO 問題の例として $2$ 値変数のポートフォリオ最適化問題がある.
• QUBO 問題はポートフォリオ最適化に限らず広く応用されているものなので, Grover Adaptive Search も様々な問題の解法として使える.

脚注

1. Austin Gilliam, Stefan Woerner, and Constantin Gonciulea. Grover adaptive search for constrained polynomial binary optimization. Quantum, 5:428, 2021.

2. Qiskit Optimization documentation: Grover Optimizer.

3. 金融のための量子アニーリングことはじめ2 | ポートフォリオ最適化 1.

4. 金融のための量子アニーリングことはじめ3 | ポートフォリオ最適化 2.

5. qulacs を使う実装や qiskit を使う実装がある.