【Qiskit】HHLアルゴリズムでポートフォリオ最適化を実装する

金融分野への量子コンピュータの応用先の1つであるポートフォリオ最適化問題をQiskitで実装してみます。

ポートフォリオ最適化とは

投資の世界ではリスク分散のためにポートフォリオを組んだ分散投資が一般的です。

巨額の資金を運用する機関投資家は最小限のリスクでリターンを最大化する最適なポートフォリオを組むことが求められます。

理論的には現代ポートフォリオ理論に基づいてモデル化出来るとされています。

ポートフォリオ最適化の定式化

現代ポートフォリオ理論では投資家が望む期待値リターンを設定したときに、それを達成するポートフォリオの中で収益率の分散が、つまりリスクが最小となるものを求めます。

\begin{aligned}
\min &~\bm{w} \Sigma \bm{w} \\
s.t. &~\bm{R^T} \bm{w} =\mu \\
&~\bm{\Pi^T} \bm{w} = \zeta 
\end{aligned}

各記号の意味は次の通りです。

Quantum Native Dojoでは問題を簡単にするため投資家の予算は考慮せずに、ポートフォリオへの投資する割合が合計で1となるような制限にしています。つまり数式で表すと下のようになります。

\begin{aligned}
\min &~\bm{w} \Sigma \bm{w} \\
s.t. &~\bm{R^T} \bm{w} =\mu \\
&~\bm{I} \bm{w} = 1
\end{aligned}

このような等式制約付きの最大最小問題はラグランジュの未定乗数法を利用することで解くことができます。ラグランジュ乗数$\eta, \theta$を用いると、次のような方程式で表すことができます。

\left( \begin{array}{ccc} 0 & 0 & \bm{R^T} \cr 0 & 0 & \bm{I^T} \cr \bm{R^T} & \bm{I^T} & \Sigma \end{array} \right) \left( \begin{array}{c} \eta \cr \theta  \cr \bm{w^T} \end{array} \right)  = \left( \begin{array}{c} \mu \cr \zeta \cr \bm{0} \end{array} \right)

よって求めたい最適$\bm{w}$がわかればいいので、表現行列$W$の逆行列が求まればよいということになります。そのような問題にはHHLアルゴリズムが有効に働きます。

W = \left( \begin{array}{ccc} 0 & 0 & \bm{R^T} \\ 0 & 0 & \bm{I^T} \\ \bm{R^T} & \bm{I^T} & \Sigma \end{array} \right)

それではQiskitを使って実装してみましょう。

Qiskitで実装してみる

では実際にQiskitを使って実装してみましょう。題材はNative Dojoと変えずに2017年のGAFAの株価時系列データを使ってポートフォリオを最適化してみます。

時系列データを入手する

まずはGAFA4社の時系列データを入手します。どんな方法でもいいですが、一例としてpandas_datareaderを使ってデータを取得します。

# 銘柄選択
codes = ["GOOG", 'AAPL', 'FB', 'AMZN'] # GAFA

# 2017年の1年間のデータを使用
start = datetime.datetime(2017, 1, 1)
end = datetime.datetime(2017, 12, 31)

# Yahoo! Financeから日次の株価データを取得
data = web.DataReader(codes, 'yahoo', start, end)

df = data['Adj Close']

データの前処理

先程のラグランジュ方程式に必要な量を求めていきます。

# 期待値リターン
daily_return = df.pct_change()
expected_return = daily_return.dropna(how='all').mean() * 252 # 年率換算のため年間の営業日数252を掛ける
# 共分散行列
cov = daily_return.dropna(how='all').cov() * 252 # 年率換算のため

これらの量を使って表現行列$W$を定義します。

R = expected_return.values
Pi = np.ones(4)
S = cov.values

row1 = np.append(np.zeros(2), R).reshape(1,-1)
row2 = np.append(np.zeros(2), Pi).reshape(1,-1)
row3 = np.concatenate([R.reshape(-1,1), Pi.reshape(-1,1), S], axis=1)
W = np.concatenate([row1, row2, row3])

np.set_printoptions(linewidth=200)

右辺のベクトルを定義します。ここでハイパーパラメータであるトータル期待リターン$\mu$は10%に設定します。

mu = 0.1 # ポートフォリオのリターン
xi = 1.0 
# 右辺ベクトル
mu_xi_0 = np.append(np.array([mu, xi]), np.zeros_like(R))

ここで今用意したベクトルは6次元のベクトルですが、量子ビットで扱うために行列の次元を$2^n$に拡張する必要があります。QiskitではHHLモジュールの.matrix_resize()メソッドを利用することで簡単に拡張することができます。

matrix, vector, truncate_powerdim, truncate_hermitian = HHL.matrix_resize(W, mu_xi_0)

HHLアルゴリズムの準備

それではHHLアルゴリズムを組み立てていきます。といってもQiskitにはHHLアルゴリズムのもモジュールが用意されているのでこちらを利用します。

処理プロセスである固有値探索と制御回転の条件を設定します。今回固有値探索に用いる量子位相推定には補助量子ビットを7量子ビット利用します。

def create_eigs(matrix, num_ancillae, num_time_slices, negative_evals):
    ne_qfts = [None, None]
    if negative_evals:
        num_ancillae += 1
        # The QFT and IQFT circuits for handling negative eigenvalues
        ne_qfts = [QFT(num_ancillae - 1), QFT(num_ancillae - 1).inverse()]

    """
    Specifically, this class is based on PhaseEstimationCircuit with no measurements 
    and has additional handling of negative eigenvalues
    """
    return EigsQPE(MatrixOperator(matrix=matrix), # The Hamiltonian Operator object
                   QFT(num_ancillae).inverse(),
                   num_time_slices=num_time_slices,
                   num_ancillae=num_ancillae,
                   expansion_mode='trotter', 
                   expansion_order=2, # The suzuki expansion order, has a minimum value of 1
                   evo_time=None,  # This is t, can set to: np.pi*3/4
                   negative_evals=negative_evals, # Set True to indicate negative eigenvalues need to be handled
                   ne_qfts=ne_qfts)

# 扱う問題の次元
orig_size = len(vector)

# 固有値探索の初期化
eigs = create_eigs(matrix, 7, 50, True)
num_q, num_a = eigs.get_register_sizes()

# 制御回転の初期化
c = 0.5*(2 * np.pi * (1. / 2**(num_a) ))
reci = LookupRotation(scale=c, negative_evals=eigs._negative_evals, evo_time=eigs._evo_time)

最後に初期状態、つまり右辺ベクトルをセットしてHHLアルゴリズムのインスタンスを生成します。

# 初期状態の設定
init_state = Custom(num_q, state_vector=vector)

algo = HHL(
    matrix,
    vector, 
    truncate_powerdim, 
    truncate_hermitian, 
    eigs, # The eigenvalue estimation instance
    init_state, # The initial quantum state preparation
    reci, # The eigenvalue reciprocal and controlled rotation instance
    num_q, # Number of qubits required for the matrix Operator instance
    num_a, # Number of ancillary qubits for Eigenvalues instance
    orig_size # The original dimension of the problem (if truncate_powerdim)
    )

HHLアルゴリズムの実行

それでは組み立てたアルゴリズムを実行してみましょう。比較のためNumPyを用いて計算した厳密解と比較してみます。

## 結果
x_HHL = algo.run(QuantumInstance(Aer.get_backend('statevector_simulator')))
## 厳密解
# x_exact = np.linalg.lstsq(matrix, vector, rcond=0)[0]
x_exact = NumPyLSsolver(matrix, vector).run()

結果から最適ポートフォリオの重みベクトルを取り出した結果がこちらのグラフです。青が厳密解、赤がQiskitで計算した結果です。プラスは買い、マイナスは空売りを表していて概ねよく結果を再現できていると言えるでしょう。

ためしに解の信頼度を求めてみると97%の信頼度と求まりました。この信頼度は補助ビットの数を変化させることで上下しますが、私のPC環境では補助ビット8以上では計算が終わりませんでした。。

# 解の信頼度を求める
def fidelity(hhl, ref):
    solution_hhl_normed = hhl / np.linalg.norm(hhl)
    solution_ref_normed = ref / np.linalg.norm(ref)
    fidelity = state_fidelity(solution_hhl_normed, solution_ref_normed)
    print("Fidelity:\t\t %f" % fidelity)

fidelity(x_HHL.solution, x_exact.solution) # => Fidelity: 0.968245

ちなみに今回組み立てた回路の情報はこちらです。使用した総量子ビットは14、量子ゲートは1252個、そのうちCNOTゲートは404個と、この程度の問題でも実機で動かすにはまだまだ大きな回路となってしまいました。

circuit_width:	 14
circuit_depth:	 1252
CNOT gates:	 404