wildcatでポートフォリオ最適化問題を試してみた

#はじめに
wildcatの使い方は https://github.com/mdrft/wildcat を見てください。
こちらも参考になるかもしれません。wildcatで4x4の数独を解いてみた
また、以前別の問題を試した記事の参考になるかもしれません、。wildcatで巡回セールスマン問題を試してみた
また、この投稿はこの勉強会の内容がネタになっています。【関東・中級】wildcatで量子アニーリング金融系アルゴリズムの実装
また、次の論文の内容に基いています。Solving the Optimal Trading Trajectory Problem Using a Quantum Annealer

#ポートフォリオ最適化問題
今回解きたい問題のハミルトニアンは次のように書けます。

H=\sum_{t}\left\{-\sum_{i}\mu_{t,i}\omega_{t,i}+M\left(\sum_{i}\omega_{t,i}-K\right)^2\right\}

第一項の不符号を除いた部分（$\sum_{t,i}\mu_{t,i}\omega_{t,i}$）が最大化したい利益になります。$\mu_{t,i}$はリターンを表し$\omega_{t,i}$は各銘柄の保有数になります。第二項は各期に$K$個まで保有できるという制約条件です。また、今回は簡単のため、リスク項と取引手数料項は除いていますが、計算の本筋は変わりません。
#QUBOに変形
ハミルトニアンを下記の形式に変形します。

H = \sum_{i,j}Q_{i,j}q_{i}q_{j}

量子ビット$q$は$(0,1)$しかとらないため、0,1しか表せません。保有数$\omega$はそうではないので少し工夫が必要です。$\omega$を複数の量子ビットを用いて次のように表します。

\omega = \sum_{a=0}^S2^{a}q_{a}

これは2進数から10進数に変換しているだけで、$S$は問題に必要な桁数を適当に選びます。また、[2]の記事でやったように次の性質を用います。

q_{i} = q_{i}^2

以上を用いることでハミルトニアンは次の形に変形できます。

H=\sum_{t}\sum_{i,a}\sum_{j,b}J_{t,i,a,j,b}q_{t,i,a}q_{t,j,b}+TMK^2

J_{t,i,a,j,b}=-2^{a}\mu_{t,i}\delta_{i,j}\delta_{a,b}+M\left(2^{a+b}-2^{1+a}K\delta_{i,j}\delta_{a,b}\right)

あとは$J$を計算した後、変形して$Q$を導きます。以上で数学の計算は終わりです。
#コード
それではコードで書いていきます。
まず、ライブラリを読み込みます。そして断熱計算の設定を行います。今回はローカルで計算をします。

from wildcat.solver.qubo_solver import QuboSolver
from wildcat.network.local_endpoint import LocalEndpoint
from wildcat.annealer.simulated.simulated_annealer import SimulatedAnnealer
from wildcat.annealer.simulated.single_spin_flip_strategy import SingleSpinFlipStrategy
from wildcat.annealer.simulated.temperature_schedule import TemperatureSchedule
from wildcat.util.matrix import hamiltonian_energy
    
schedule = TemperatureSchedule(initial_temperature=1000, last_temperature=0.1, scale=0.99)
strategy = SingleSpinFlipStrategy(repetition=1000)
annealer = SimulatedAnnealer(schedule=schedule, strategy=strategy)
local_endpoint = LocalEndpoint(annealer=annealer)

次に問題設定をします。

#行列muを作成
mu = np.array([
    [0.0,0.2,0.3],
    [0.2,0.2,0.1],
    [0.1,0.3,0.3],
    [0.4,0.1,0.5],
    [0.5,0.3,0.1]
])

# T : 期間
T = mu.shape[0]
# N : 銘柄の数
N = mu.shape[1]
# K : 選ぶ銘柄の数
K = 10
# S : ２進数の桁数
S = 5
# M : 定数
M = 1

それでは$J$を計算します。

#クロネッカーデルタを定義
def delta(i,j):
    if(i == j):
        return 1
    else:
        return 0

# 行列Jを計算
J = np.empty((T,N,S,N,S),dtype = np.float64)
for t in range(T):
    for i in range(N):
        for a in range(S):
            for j in range(N):
                for b in range(S):
                    J[t,i,a,j,b] = -2**a * mu[t,i] * delta(i,j) * delta(a,b) + M * (2**(a + b) - 2**(1 + a) * K * delta(i,j) * delta(a,b))


#定数項の計算
C = T * M * K**2

$J$を変形して$Q$を作ります。

#行列Qを計算
Q = np.zeros((T*N*S, T*N*S), dtype = np.float64)
for t in range(T):
    for i in range(N):
        for a in range(S):
            for j in range(N):
                for b in range(S):
                    Q[a + i*S + t*S*N,b + j*S + t*S*N] = J[t,i,a,j,b]

それでは$Q$が出来たので、最適化計算を行います。今回は最適化計算を繰り返し、最小のエネルギーが得られた結果を採用することにします。（エネルギーを計算する関数はwildcatのソースコードを参考にしています。）

#最適化計算の定義
solver = QuboSolver(Q)
def callback(arrangement):
    e = solver.hamiltonian_energy(arrangement)

#(-1,1)->(0,1)に変える関数
def chgangeExpression(vector):
    D = vector.shape[0]
    for x in range(D):
        if(vector[x] == -1):
            vector[x] = 0
    return vector

#エネルギーを計算する関数
def calcEnergy(vector, matrix, const):
    return vector.T.dot(matrix).dot(vector) + const

#繰り返し計算し、最小エネルギーを探す
MinEnergy = 0
MinArrangement = np.zeros(N*T, dtype = np.int)
for i in range(1000):
    arrangement = chgangeExpression(solver.solve(callback, endpoint=local_endpoint).result())
    if(i == 0):
        MinEnergy = calcEnergy(arrangement, Q, C)
        MinArrangement = arrangement
    else:
        if(calcEnergy(arrangement, Q, C) < MinEnergy):
            MinEnergy = calcEnergy(arrangement, Q, C)
            MinArrangement = arrangement

結果は次のようになります。

print("Energy: ", calcEnergy(MinArrangement, Q, C))
print("Spins: ", MinArrangement)

Energy:  -12.399999999999977
Spins:  [0 0 0 0 0 1 1 1 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1
 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 1 0 1 0
 0]

これを見やすい形に変形すると次のようになります。

result = np.zeros((T, N), dtype = np.int)
i = 0
for t in range(T):
    for i in range(N):
        for a in range(S):
            result[t, i] += 2**a * MinArrangement[a + i*S + t*S*N]
print("Result: \n", result)

Result: 
 [[0 7 3]
 [2 3 5]
 [8 2 0]
 [4 1 5]
 [5 0 5]]

これは縦軸が時間軸で、横軸が銘柄です。時間は上から下へ進みます。これからわかるように、各期での保有数は10個になっており、問題設定の制約条件を満たしていることをが分かります。また、問題設定の$\mu$の値からわかるように、リターンが大きいものをなるべく多く保有している傾向が見えなくもないです。
興味がある方はぜひ試してみてください。コードは上げておきます。https://github.com/Yokomakura/wildcat-portfolio-optimisation
急ぎ足でやったので、ミスってたらごめんなさい。
#参考

1. wildcatで4x4の数独を解いてみた
2. wildcatで巡回セールスマン問題を試してみた
3. Solving the Optimal Trading Trajectory Problem Using a Quantum Annealer