Frank-Wolfe法でポートフォリオ最適化

本記事は[数理最適化 Advent Calendar 2023](https://qiita.com/advent-calendar/2023/mathematical-optimization)の17日目の記事です。

【注意：途中で投資信託のデータセットを使いますが、本記事は現実の投資行動を推奨するものではありません。】

# 0. はじめに

連続最適化が好きです。（皆さんも好きかと思います。）

特に、凸最適化アルゴリズムが好きです。

ところが、日常の問題解決では、お気に入りのアルゴリズムが役立つ場面は案外少ないです。なぜなら、締切、リソース、歴史的経緯、やる気など、数理的な定式化とは関係ないさまざまな制約が存在するために、肝心の最適化部分は、汎用ソルバーや適当なヒューリスティクスで済ませてしまうことも多いからです。

しかし、ここに往年のポケットモンスターの名言があります。

*つよいポケモン　よわいポケモン
そんなの ひとの かって
ほんとうに つよい トレーナーなら
すきな ポケモンで かてるように がんばるべき* [1]

つまり、強い心をもつ最適化民なら、好きなアルゴリズムで問題が解けるように頑張るべきと言っているわけですね。身に染みます。

というわけで、今回は**好きなアルゴリズムで好きな問題を解いてみました。**

## TL;DR

- Frank-Wolfeアルゴリズムが好き
- ポートフォリオ最適化が好き
- Frank-Wolfeアルゴリズムでポートフォリオ最適化を解いてみた
- 実問題では特定の手法や定式化にこだわらず、目的に沿ったものになっているか吟味しよう

# 1. Frank-Wolfeアルゴリズム

突然ですが、[Frank-Wolfeアルゴリズム](https://en.wikipedia.org/wiki/Frank%E2%80%93Wolfe_algorithm
)が好きです。

以下のような制約つき最小化問題があります。

```math
\displaylines{
\min_x  \quad f(x) \\
\text{s.t.} \quad x \in D
}
```

ここで$D$はコンパクト凸集合、$f$は微分可能な凸関数としておきましょう。つまり、最適化アルゴリズムの中で勾配 $\nabla f(x)$ の情報を使ってもよいです。最適化の民の仕事は、$D$や$f$の具体的な形を見て、どのようなアルゴリズムを使ったら良いか検討することです。

さて、Frank-Wolfeアルゴリズムというのは以下のようなアルゴリズムです。（画像は[2]から引用）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/3205e49c-eaaa-7049-743d-57a743447239.png)

つまり、

1. 現在の点 $x^{(k)}$ において線形化した目的関数を最小化する点 $s$ を探す
1. $x^{(k)}$ と $s$ の内分点をとる

という2ステップを収束まで繰り返すアルゴリズムです。重要な要素はステップ1で、**線形化した問題**
```math
\displaylines{
\min_s  \quad \langle s, \nabla f(x^{(k)}) \rangle \\
\text{s.t.} \quad s \in D
}
```
を毎回解かなければならないため、この部分問題が高速である必要があります。これは集合 $D$ に関する性質になります。どのような集合だと速く解けるでしょうか？

実は、日頃よく使われる制約の中にも線形最適化が速く解けるものが多いです。例えば、L1ノルム球、一般のLpノルム球、確率単体、行列のトレースノルムや作用素ノルム球などでは、この線形最適化が容易に解けます [2]。

代表的な応用例は、スパース推定や行列分解、行列補完などです。例えば、制約型のラッソ

```math
\displaylines{
\min_\beta \quad \| y - X \beta \|^2_2 \quad
\text{s.t.} \quad \| \beta \|_1 \leq t
}
```

に対するFrank-Wolfeアルゴリズムを具体的に考えてみると、forward stepwiseとかforward greedyなどと呼ばれる貪欲アルゴリズムとほぼ同じものが出てきます （具体的には [3] などを参照）。

なお、ステップ2におけるステップ幅 $\gamma$ は、上の画像のように決定的な減衰率 $\gamma = 2 / (k + 2)$ を使う版と、毎回ラインサーチで決める版があります。

個人的になぜこのアルゴリズムが好きかというと、統計学や機械学習[^1]、離散最適化[^2][^3]など、いろいろな分野と関係しているからです。ただ、やや本題から逸れそうなので、脚注に書いておきます。

[^1]: **統計学と機械学習:** 元のアルゴリズム自体は1956年の論文[4]で提案されているようですが、機械学習コミュニティでFrank-Wolfeアルゴリズムを有名にした論文は、おそらくJaggi (2013) [2] でしょう。当時は、機械学習の国際会議でも、スパース推定やスパース行列分解のような問題への関心が高かったですが、スパース推定では何らかの形で貪欲法的な要素があるアルゴリズムが色々提案されていました。こういうった応用事例に対して、Frank-Wolfe法として統一的に収束保証を与えたというのが注目ポイントだったのだと思います。
[^2]: **貪欲法:** L1ノルム球、L∞ノルム球、確率単体、少し複雑な例だと[permutohedron](https://en.wikipedia.org/wiki/Permutohedron)など、線形関数の最大化 or 最小化がなぜか容易にできる集合があります。実はこの性質は劣モジュラ関数と深い関係があります[5]。というのも、これらの集合は特定の劣モジュラ関数の基多面体 (base polyhedron) ないしそれらを貼り合わせたものになっており、こういった集合上での線形最適化は $O(n)$ ($n$ は次元数) の貪欲法でできるということが知られています。
[^3]: **離散最適化:** 連続最適化だけでなく、離散最適化にも役立っています。オリジナルのFrank-Wolfeそのものではないのですが、LMOを使った類似のアルゴリズムとしてWolfe (1976) の最小ノルム点アルゴリズムというものがあります。これを利用して、劣モジュラ関数を最小化するアルゴリズムを構成できます (Fujishige-Wolfe)。興味のある人は[この論文など](https://arxiv.org/abs/1411.0095)を参照してください。

# 2. ポートフォリオ最適化

ポートフォリオ最適化も好きです。

ポートフォリオ最適化というのは、金融工学の一分野で、与えられた効用関数やリスク関数を最適化するように資産配分を決定する問題です。数理最適化 (OR) 分野でも、ポートフォリオ最適化は典型的な応用先のひとつになっています。

リスク資産が $J$ 個あり、その上の重み（ポートフォリオ）$w_1, \ldots, w_J$ を決めたいとします。リスク資産については、投資リターンの期待値 $\mu_j$ や、資産どうしの共分散行列 $\Sigma$ などがわかっているものとします。[^4]

[^4]: 現実には将来の期待リターンを見積もるのは大変難しいですが、本記事ではあくまで最適化にスポットライトを当てるということで、そこは気にしないことにします。現実の投資と最適化問題との間にどのくらいギャップがあるかについて真面目に考察し始めると、それだけで無限に研究ができてしまうので、沼です。

例えば、平均・分散ポートフォリオなら、以下のように書けます。
```math
\displaylines{
\min_w \quad \frac{\lambda}{2} w^\top \Sigma w - w^\top \mu\\
\text{s.t.} \quad \sum_{j=1}^J w_j = 1,
\quad w_j \geq 0 \quad (j = 1, \ldots, J)
}
```
つまり、ポートフォリオ $w$ が確率単体に含まれているという制約のもとで、ポートフォリオのリターンの分散 $w^\top \Sigma w$ をなるべく小さく、期待リターン $\mu^\top w$ をなるべく大きくするという問題で、2次計画になります。(説明の都合上、[Wikipedia](https://ja.wikipedia.org/wiki/%E7%8F%BE%E4%BB%A3%E3%83%9D%E3%83%BC%E3%83%88%E3%83%95%E3%82%A9%E3%83%AA%E3%82%AA%E7%90%86%E8%AB%96) などで紹介されている定式化とは多少違うものにしました。)

また、データドリブンにポートフォリオを決める手段として、**フルスケール最適化** (FSO)と呼ばれる定式化もあります。
```math
\displaylines{
\min_w \quad - \sum_{t=1}^T U(r_t^\top w)\\
\text{s.t.} \quad \sum_{j=1}^J w_j = 1,
\quad w_j \geq 0 \quad (j = 1, \ldots, J)
}
```
ここで、$r_t$ は時刻 $t$ における $J$ 個の資産のリターンを並べたものです。$U$ は効用関数です。つまり、FSOは過去リターンのデータ $r_1, r_2, \ldots, r_T$ の上での効用関数をなるべく大きくするようにポートフォリオ $w$ を決定する問題といえます。[^5]

[^5]: FSOの詳しい背景については、手前味噌で恐縮ですが論文 [6] などを参照してください。機械学習などに親しんでいる人には、一見自然な定式化に見えますが、ポートフォリオ最適化の分野では意外とマイナーな手法です。理由としては、累積リターンや富 (wealth) ではなく、一時刻ごとのリターンに対する効用を考える流儀がどうやら珍しいようです。

なお、上の定式化では非負制約 $w \geq 0$ が入っているので、ロングオンリーポートフォリオなどと呼ばれます。これは、資産を買い持ちしかできないという意味です。現実の制約として資産の空売りができない場合、例えば長期的な年金運用や、投資信託を買う際などに当てはまる問題設定といえます。

# 3. Frank-Wolfe法でポートフォリオ最適化を解く

さて、Frank-Wolfeアルゴリズムを現実の問題解決に応用してみたいと思います。上で紹介したFSOの定式化は、確率単体上での最小化になっているので、$-U$ が凸関数ならFrank-Wolfeが使えるだろう、という目論見があります。

簡単のために指数効用関数 $U(c) = - \exp(-\alpha c)$ を使うことにします。$\alpha > 0$ はリスク回避度（ハイパラ）ですので、今回は適当に1から30くらいに設定しておきましょう。まとめると、解きたい問題は以下のようになります。

```math
\displaylines{
\min_w \quad f(w) := \sum_{t=1}^T \exp(- \alpha r_t^\top w)\\
\text{s.t.} \quad \sum_{j=1}^J w_j = 1,
\quad w_j \geq 0 \quad (j = 1, \ldots, J)
}
```

## 3.1 アルゴリズムの細部

Frank-Wolfeアルゴリズムの具体的な更新則を確認しておきます。

$k$ 回目の更新において、目的関数 $f$ の勾配は次のようになります。
```math
\nabla f(w^{(k)}) = \sum_{t=1}^T -\alpha \exp(-\alpha r_t^\top w^{(k)}) r_t
```

次に、線形化した目的関数 $\langle s, \nabla f(w^{(k)}) \rangle$ を確率単体の上で最小化したいです。これはとても簡単で、解は勾配の値が最小となる成分が1となるようなone-hotベクトルになります。証明としてはKKT条件を書いても確認できますが、個人的には貪欲法 [^2] と考えると素直に納得できます。

したがって、アルゴリズムの全体像としては「**効用関数の勾配の値が最も増えるような成分を各更新のたびに1つだけ増やす**」というものになるため、**スパースなポートフォリオが得られるだろう**と考えられます。

## 3.2 実験

### データセット

【再注意: 本記事は現実の投資行動を推奨するものではありません。】

現実でも考えそうな問題設定として、投資信託のデータを利用しましょう。ここでは、2024年度に始まる新NISA制度で積み立てることができる投資信託にどのように資金を配分したらよいか、という問題を考えてみます。ただし、投資信託だけでも数千個あるので、対応する最適化問題は数千次元になってしまいます。そこで、今回は簡単のため100次元程度まで減らしたいと思います。一般投資家の心情としてありそうな規準として、

- (1) 純資産総額10,000百万円以上
- (2) 信託報酬0.5%未満
- (3) 2023年12月時点の基準価格10,000円以上
- (4) 2年分以上データが取れる

であるものを適当に選び、$J = 102$ まで絞りました。基準価格のデータは、証券会社のウェブサイト等でCSVデータをダウンロードできる場合がありますので、それを利用します。今回は手作業でひとつずつダウンロードしています。

データ期間としては、2021年1月1日〜2023年12月1日までの約3年弱の日次データを使いました。実は、データ期間の選択が結果に対して大変センシティブに影響するのですが、それについては記事の最後に言及します。

まずは資産リターンの相関係数を見てみます。

![correlation.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/4b88ebd6-0e15-cfd2-e212-5a3fe51d3dd0.png)

あらやだ真っ赤。

102銘柄のほとんどが正の相関関係にあり、相関係数が大きい (0.7以上) の銘柄ペアもたくさんあります。原因としては、抽出条件として純資産総額が多く、かつ信託報酬が低いものを選んでいるため、米国や日本の株式指数への連動を目指す似たようなインデックス投信がたくさん釣れてしまったからですね。実際、コンセプトが重複している投資信託が多数含まれており、例えばS&P500連動銘柄が6つあったりします。

**スパース最適化としては良い問題かもしれません。**

ポートフォリオ最適化は本来「どのように分散投資すればよいか」を定式化したものです。しかし、S&P500に連動する投資信託を複数購入するのは現実的に無駄ですから、このグループからは高々1つ選べば十分のはずです。

### 実装

単純なアルゴリズムですので、パフォーマンスなどは気にせず、Pythonでさっと実装してみます。

```python
from typing import Any, Dict, Optional
import numpy as np
from scipy.optimize import minimize_scalar


def loss(
    w: np.ndarray,
    ret: np.ndarray,
    risk_aversion: float = 2.0,
) -> float:
    """目的関数"""
    return float(np.exp(-risk_aversion * (ret @ w)).sum())


def grad(
    w: np.ndarray,
    ret: np.ndarray,
    risk_aversion: float = 2.0,
) -> np.ndarray:
    """目的関数の勾配"""
    return -risk_aversion * (
        np.exp(-risk_aversion * (ret @ w)).reshape(-1, 1) * ret
    ).sum(axis=0)


def line_search_obj(
    gamma: float,
    w_curr: np.ndarray,
    w_prev: np.ndarray,
    ret: np.ndarray,
    risk_aversion: float,
) -> float:
    """ラインサーチの目的関数"""
    return loss(
        (1.0 - gamma) * w_prev + gamma * w_curr,
        ret=ret,
        risk_aversion=risk_aversion,
    )


def line_search(
    w_curr: np.ndarray,
    w_prev: np.ndarray,
    ret: np.ndarray,
    risk_aversion: float,
) -> float:
    """ラインサーチ"""
    result = minimize_scalar(
        line_search_obj, bounds=(0.0, 1.0), args=(w_curr, w_prev, ret, risk_aversion)
    )
    return result.x


def fw_exponential(
    ret: np.ndarray,
    w0: Optional[np.ndarray] = None,
    max_iter: int = 1000,
    risk_aversion: float = 2.0,
    use_line_search: bool = False,
) -> Dict[str, Any]:
    """Frank-Wolfeアルゴリズムによって指数効用関数を最適化する。

    Args:
        ret: 日次リターンの配列 (T * J)
        w0: 初期値 (J,)
        max_iter: 反復回数
        risk_aversion: リスク回避度 (alpha)
        use_line_search: Frank-WolfeアルゴリズムのStep2でラインサーチをする
    """
    J = ret.shape[1]

    if w0 is None:
        w0 = np.ones(J, dtype=ret.dtype) / J

    history = {
        "w": [],
        "loss": [],
        "grad": [],
    }
    w = w0
    history["w"].append(w.copy())
    history["loss"].append(loss(w=w, ret=ret, risk_aversion=risk_aversion))
    for i in range(max_iter):
        # 勾配を計算する。
        g = grad(w, ret, risk_aversion)

        # Step 1: 確率単体上の線形最適化
        # ここでは、勾配の成分のargminを計算するだけ！
        idx = np.argmin(g, keepdims=True)[0]

        s = np.zeros(J, dtype=ret.dtype)
        s[idx] = 1.0


        # Step 2: 固定ステップ幅での更新、またはラインサーチ
        if not use_line_search:
            gamma = 2.0 / (2.0 + i)
        else:
            gamma = line_search(s, w, ret, risk_aversion)
        w = (1.0 - gamma) * w + gamma * s
        
        history["grad"].append(g.copy())
        history["w"].append(w.copy())
        history["loss"].append(loss(w, ret, risk_aversion))
    return history
```
上のコードで`fw_exponential`がコア部分で、今回は記事の可視化のためにヒストリーも保存するようにしてます。ちなみにステップサイズ $\gamma$ をラインサーチする版も一応作りましたが、実際に動かすと定性的な結果はほとんど変わらなかったので、記事内での報告は省略します。

試しに $\alpha = 1$ で動かしてみましょう。

![loss.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/68600065-d3ef-2ec9-f2f5-353be9865adc.png)

目的関数 $f(w)$ はかなり高速に収束して動かなくなっています。得られた解 $w$ をみてみましょう。確率単体上の最適化ですので、わかりやすくするために円グラフで表示してみます。

![a01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/50fe0548-f379-019b-3e3a-8a7a60f650ee.png)

**おお〜、スパース！！**

最適化問題としては102次元あるのですが、解の非ゼロ成分は3次元（実質2銘柄）しかありません。ちゃんとFrank-Wolfeに期待すべき挙動になっている気がします。

他の $\alpha$ でもやってみましょう。

![a10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/bfabe5a6-edbd-ee15-1008-2b7c05e93357.png)
![a20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/373311bd-4dcc-6fca-0930-09578dff1de1.png)
![a30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/166350/1f2fd005-fe3d-825b-349a-8011b9427de2.png)

$\alpha$ を増やしていくと非ゼロ成分が増えていく様子がわかります。以下、考察です。

- $\alpha$ は「リスク回避度」として解釈できる。$\alpha$ が小さいときは、データ期間中のリターンが大きい銘柄を貪欲に選択する。$\alpha$ が大きいときは、相関が低い複数の銘柄に分散投資することでリスクを抑えることを考えるようになる。特に $\alpha$ が大きいときには、国内株式、米国株式、債券、ゴールドなど異なる資産クラスに投資するポートフォリオが作られている。
- 似ている銘柄は1つしか選ばない。例えば「S&P500」「NASDAQ」「TOPIX」「ゴールド」といったコンセプトのインデックス投信は102銘柄中に重複して存在するが、Frank-Wolfeアルゴリズムから得られる解ではその中から1つずつのみ選択されている。

# 4. おわりに

好きなアルゴリズム（Frank-Wolfe）で好きな問題（ポートフォリオ最適化）を解いてみました。

幸い、今回の実験では、「アルゴリズムがこう動くだろう」という期待から外れない結果が得られたように思います。

しかし、現実の意思決定に数理最適化アルゴリズムを使う際には、**実際にやりたいことの要件が最適化問題の定式化から漏れていたケース**が往々にしてありますので、注意が必要です。今回扱ったのはポートフォリオ最適化問題なので、最終的に役立つかどうかの基準は、「得られた最適解に対して、自分が実際に投資したいかどうか」です。そのような観点で改めて考えると、以下のようなことから、今回の実験結果だけではまだ不十分ではないかと思います。

- **定式化やアルゴリズムの性質が自分のニーズと合致しているか？:** データドリブンの意思決定は、訓練データ期間で偶然観察されなかった性質や、データに含まれない見通し（＝やりたいこと）を見落としてしまう危険性があります。 例えば、本来は意味が異なるはずの資産（株と債券など）が、訓練データ中でたまたま非常に高い相関を持ってしまったとします。このとき、Frank-Wolfeアルゴリズムはそのうち1つしか選択してくれませんが、意思決定者の気持ち的には、スパースではなく50:50で持ちたかった、という場合もあるでしょう。こういった要件がある場合は、Frank-Wolfeではなく、内点法[^6]など、よりdenseな解を得やすい最適化アルゴリズムの方が好ましいでしょう。
- **入力データや定式化の変化に対して結論が安定しているか？:** 最適化アルゴリズムを現実の意思決定に用いる際には、そこから得られる結論がデータや定式化の入れ替えに対してどれだけ揺るぎないものなのか、ということも吟味する必要があるでしょう。ポートフォリオ最適化の出力は訓練データ期間に対して非常にセンシティブです。上では恣意的に2021年〜2023年のデータを使っていますが、ここで2020年を含めるとほとんど株（NASDAQ）になったり、2022年以降のみに絞るとほとんどコモディティ（ゴールド）になったりします。したがって、今回の定式化で得られるポートフォリオは、あくまで過去のスナップショットであり、訓練データ期間の市況を反映しただけのものとも言えます。[^7]　また、今回の実験でも観察されたように、効用関数やハイパラ ($\alpha$) の選択によっても大きく異なる解が得られます。

[^6]: 具体的には、Rでよく使われている[Rsolnp](https://cran.r-project.org/web/packages/Rsolnp/index.html)などです。
[^7]: ポートフォリオ最適化の典型的な定式化がデータの時系列性を無視しているため、とも考えられます。なお、この問題意識は論文[6]でも取り扱われています。

以上です。
「好きなアルゴリズムで問題が解けるように頑張るべき」、というよりは、アルゴリズムも適材適所、というスタンスでいるのが大人な最適化民になるコツかなと感じています。

## 参考文献

[1] ポケットモンスター 金・銀 (1999)
[2] Jaggi (2013). Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization. Proceedings of the 30th International Conference on Machine Learning, PMLR 28(1):427-435.
[3] Shalev-Shwartz, Srebro, and Zhang (2010). Trading Accuracy for Sparsity in Optimization Problems with Sparsity Constraints. SIAM Journal on Optimization, 20(6):2807-2832.
[4] Frank and Wolfe (1956). An algorithm for quadratic programming. Naval Research Logistics Quarterly. 3(1–2):95–110.
[5] Edmonds (1970). Submodular Functions, Matroids, and Certain Polyhedra.
[6] 南 賢太郎, 今城 健太郎, 中川 慧, 今長谷 拓 (2022). 予測型フルスケール最適化による資産配分. https://www.jstage.jst.go.jp/article/pjsai/JSAI2022/0/JSAI2022_2J4GS1002/_article/-char/ja/