💹

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

2023/02/24に公開

はじめに

本記事では論文[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.
  1. O: オラクル作用素.
    探索対象番号列 I \subset \{ 0, \dots, 2^n -1 \} に対してオラクル作用素 OI に属する番号の基底ベクトルは -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.
  1. 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}

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 に達したところで計算を終了する.

QUBO 問題への応用

QUBO 問題

2n 変数の 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 としたものである.

ユニタリ作用素 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}_nJ に含まれる番号 j の量子ビットたちを制御ビットとして, 出力レジスター \ket{z}_m にユニタリ作用素 U_{G}\left(\frac{2\pi}{2^m} a_J \right) を実行する作用素である.

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

上記の単項式の重ね合わせ状態を用意するため, 出力レジスター \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].

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}

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}

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

まとめ

  • 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 を使う実装がある. ↩︎

Discussion