元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

Transformerを理解したい

Transformerを理解したい記事です。

Transformerが提案された論文「Attention is All you need」をベースにまとめていきます。

以下の知識があることを前提としています。

  • Deep Learningに関する基礎的な知識
  • RNNに関する基礎的な知識
  • 上記に関連する数学的知識

大雑把に言えば、Deepは使っているけどTransformer関連だけは経験がないような人向けになります。

Transformer誕生の背景

Transformerは、自然言語処理など系列モデリングの文脈で登場したモデルです。

入力系列を別の系列に変換するタスクはSeq2Seqと呼ばれ、入力を隠れ状態にエンコードするEncoderと、出力系列へ変換するDecorderからなるため、Encoder-Decoderモデルとも呼ばれます。

Transformer登場以前のDeep LearningによるEncoder-DecoderモデルにはRNNベースのモデルが使われていました。しかし、RNNベースのEncoder-Decoderモデルには以下の問題がありました。

  1. Encoderの出力が固定長ベクトルになる。
  2. Decoderへの入力が、Encoderから出力された最後の隠れ層のみ。

これらの問題を解決したのがTransfomerです。そして、その本質は、Transformerブロック内部にあるAttention機構にあります。以下でまとめていきます。

Transformerの構成要素

Transformerは以下のような構造をしています。

[1] Figure 1 より引用

Nxは、このブロックがN回繰り返されるという意味です。

Transformerを理解するには、この図にあるブロック

  1. Multi-Head Attention
  2. Positional Encoding
  3. Feed Forward
  4. Add & Norm

を理解する必要があります。また、Multi-Head Attentionを理解するためには、Attentionを理解しておく必要があります。

Attention

Multi-Head Attentionの前に、Attentionについて知っておきましょう。Attentionとその派生形であるMulti-Head Attentionは以下のような構造をしています。

[1] Figure 2より引用

計算内容はともかく、Multi-Head Attentionは複数個のAttentionの計算結果を結合したものになっているのがわかると思います。

Attentionの概念

Attentionは入力のどの部分に注目するべきか選択する機構で、以下の式で与えられます。

\begin{align} \text{Attention}(Q, K, V) = \text{softmax}\left(\frac{QK^T}{\sqrt{d_k}}\right)V \label{Eq:attention} \end{align}

 Q K Vはそれぞれクエリ、キー、バリューと呼ばれる行列です。

一旦数学のことは忘れて用語だけ説明すると、クエリは検索キーワード、キーとバリューは辞書に対応します。つまり、入力された検索キーワードにもっとも近しいキーを選び、対応するバリューを得るということを上記の式は表しています。

  • 近しいキーを選ぶ → QとKの内積を取ってsoftmax
  • 対応するバリューを得る → softmaxの結果とVの内積を取る

ということです。キーを選ぶといっても本当に1つの要素だけを選ぶと微分不可能になる不都合があるので、softmaxで代用しています。softmax後に最も値が大きい要素を選ぶと、我々が辞書を引くのと同じような「キーを選ぶ」という操作になります。「 K, V はどう決めるんだ」、「 \sqrt{d_k}ってなんだ」という声が聞こえてきそうですが、この後説明します。

Scaled Dot-Product Attention

さて、式 \eqref{Eq:attention} はScaled Dot-Product Attentionというのが正式名称で、一般的にAttentionと呼ばれるものは大体これです。

Dot-Productというのが 式 \eqref{Eq:attention} における  QK^T のことで、Scaled というのが  \sqrt{d_k} のことですね。ここで、 \sqrt{d_k} はキー  K の次元を表します。なぜこんな計算をするのか理解するために式 \eqref{Eq:attention} を視覚的に確認しましょう。

まず、  QK^T の Dot-Product(内積)では以下の計算を行っています。簡単のためにひとつのクエリで考えます。

 QK^Tの新しくAttention weightというものが登場していますが、要はクエリとキーの類似度を表しています。与えられたクエリに対して類似度を計算し、softmaxを取ることでクエリに対するキーの重みを得ます。そして、得られたAttention weightでバリューの加重平均を取ります。

簡単のため、Attention weightの要素を  w_{i,j} = \boldsymbol{q}_i \cdot \boldsymbol{k}^T_j と置きました。Attention weightのいずれかの要素が  1 それ以外は  0 の場合、クエリで検索したキーに対応するバリューを選んでいることになります。

"Scaled"の謎

 n 次元単位超球を考えましょう。ナニソレという方は、円と球をイメージしたら良いです。半径はいくつになるでしょうか。単位円の半径は  \sqrt{1^2 + 1^2} = \sqrt{2}、単位球の半径は  \sqrt{1^2 + 1^2 + 1^2} = \sqrt{3} になりますね。 n 次元だと  \sqrt{1^2 + 1^2 + \ldots + 1^2} = \sqrt{n} になります。つまり、次元が大きくなるとこの半径はどんどん大きくなってしまいます。

ここで、  QK^T の計算を振り返ってみると、Attention weightの次元はキーの次元で決まっていることがわかると思います。我々がこの計算で知りたいのは、クエリとキーの類似度です。次元の大きさによって計算結果が変化してしまうと困るので、 \sqrt{d_k} でスケーリングするわけですね。

Multi-Head Attention

さて、ここまでAttentionことScaled Dot-Product Attentionの説明をしてきました。複雑に見えて原理は単純な事がわかったかと思います。そしてここまで学習可能なパラメータが一切出てきていないことに気が付いたでしょうか?つまり、単純な内積と加重平均のみでAttention機構は実現されています(すごい)。ただ、このままでは当然モデルの学習はできません。Attentionを学習可能な形にしたものがMulti-Head Attentionというわけです。

ここで、Multi-Head Attentionの図を見直してみましょう。Attention layerの入力前にLinear layerによる線形変換が行われていることがわかります。この部分がAttention機構の学習を担ってくれています。そして、 h 個のAttention(それぞれをheadと呼ぶ)の結果を結合するのでMulti-Head Attentionと呼びます。式で書くと以下のようになります。

\begin{align} \text{MultiHead}(Q, K, V) &= \text{Concat} (\text{head}_1, \ldots, \text{head}_h) W^O \\ \text{where head_i} &= Attention \left(QW_i^Q, KQ_i^K, VW_i^V\right) \end{align}

ここで、 W_i^{\left\{Q,K,V\right\}} が図のLinear layerに対応します。Concatでheadを1行に結合して、 W_O でheadの加重平均を取っています。

Concatするところで出力の次元がheadの数だけ増えてしまいそうですが、Linear layerのところで次元を  1 / h に落としています。

なぜMulti-Headにしたかなのですが、その方が性能が良かったからと著者らは言っています。異なる位置の異なる部分空間表現を学習できるからとのことです(特にこれ以上のことは書かれていませんでした)。

K、Vの決め方

今まで  Q, K, V は暗黙的に与えられるものとして来ましたが、実は基本的に  Q = K =V です。つまり、 Q, K, Vは入力それ自身なのでした。

このような場合のAttentionを特にSelf-Attentionと呼びます。 Q はクエリなので、入力であることに違和感はありませんが  K, V も同じってどういう事なのでしょうか。

 Q=K=V の状況で式 \eqref{Eq:attention} の 計算を考えてみると、各クエリと似ているキーは、入力クエリと同じ値を持つキーになり、対応するバリューは結局入力したクエリになるという意味の分からない式になってしまいます。

Multi-Head AttentionはLinear layerによって  Q, K, V を別々に変換するので、このような状況にならず入力系列間の関係性を学習することができます。

Positional Encoding

Multi-Head Attentionは系列の要素間の関係を学習することができますが、このままでは系列の順序まで考慮してくれません。そこで、Positional Encodingという手法で位置情報を入力へ足しこみます。

Position Encoding  PE は以下の式で与えられます。

\begin{align} PE_{\left(pos, 2i \right)} &= \sin \left(pos / 10000^{2i / d_{model}} \right) \\ PE_{(\left(pos, 2i + 1 \right)} &= \cos \left(pos / 10000^{2i / d_{model}} \right) \end{align}

ここで、 pos は系列における位置(キーで例えるなら何番目のキーか)、 2i 2i + 1 は要素ベクトルの次元(キーで例えるなら  pos 番目のキーベクトルの何番目の要素か)を表します。[1] では、"The wavelength form a geometric progression from  2\pi to  10000 \cdot 2\pi" と言っていることから、 0 \le 2i \le d_{model} であることがわかります。

式で考えるより、グラフを見た方がわかりやすいのでこちらをご覧ください。

グラフを出力するコードは本記事の最後に記載しました。さて、なぜこれで位置をエンコーディングできるのでしょうか。グラフを解剖していきましょう。

位置のエンコード

positionごとに値をいくつか取り出すと以下のようになっています。

位置が変わるごとにくびれている位置が変化していっていますね。くびれの位置を見ればどの位置かわかります。

次元のエンコード

なんと、系列の位置だけでなく次元までエンコードているのがPositional Encodingの凄いところ。次元方向に値をいくつか取り出すと以下のようになっています。

周期を見ればどの次元かわかるようになっています。

Add & Norm

Add & Normalization layerの略です。名前の通りskip-connection由来の入力を足して、Normalizationする層です。Transformerで用いられているようなskip-connectionはresidual connectionと呼ばれており、ResNetの思想を受け継いだ構造になっています。residual connectionを入れることでより深い層での学習を可能にしています。

Residual connection

residual connectionを知っている人は読み飛ばして大丈夫です。

ResNet登場以前はVGG16、VGG19やGoogLeNetなどが深い層を持つモデルとされていたのですが、せいぜい20層程度の深さが限界でした。当時、Batch normalizationの登場により勾配消失・爆発問題は防ぎやすくなっていたものの、これ以上層を深くすると性能が劣化(degradation)する問題があったためです(下図)。

[2] Figure 1より引用

ResNetではこのresidual blockを用いることで、100以上の層を持つモデルでの学習を可能にしました(下図)。

左が単なる直列のblock、右がresidual blockです。Identity mapping(恒等写像)が追加されただけですね。直感的理解としては、追加の層が不要でも  \mathcal{F}(x) が0になれば良いだけなので層を深くしても問題ないということになります。

Layer Normalization

Normalizationにはlayer normalizationが用いられています(下図)。

[3] Figure 2より引用

 N, C, H, W は順に、バッチサイズ、チャネル数、高さ、幅を表しています。画像が前提の説明なので  H, W がありますが、系列モデリングを前提としたTransformer では  H, W が系列の位置に相当します。この図では青のブロックの範囲でnormalizationをすることを示しています。

画像系から入ってきた人にはBatch Normalizationが馴染み深いと思うのですが、バッチサイズが小さいと値が不安定になる問題がありました。

Layer Normalizationはバッチサイズに依存せず、一つのlayer内で正規化を行うことでこの問題を克服していることがわかると思います。

Feed Forward

図ではFeed Forwardとだけ書かれていましたが、[1]ではPosition-wise Feed-Forward Networksと呼んでいます。以下の式で与えられます。

\begin{align} \text{FFN}(\boldsymbol{x}) = \text{ReLU}(\boldsymbol{x}W_1 + \boldsymbol{b}_1)W_2 + \boldsymbol{b}_2 \end{align}

 \boldsymbol{x} が小文字であることに注意してください。position-wiseという名前の通り、系列の位置ごとに同一の関数を適用するということです( W_1, \boldsymbol{b}_1, W_2, \boldsymbol{b}_2 は各位置で共通ということです)。

まとめ

以上がTransformerを構成するブロックの全貌になります。Encoder、Decoder側でブロックを繰り返してやればTransformerの完成です(Decoder側2つ目のMulti-Head AttentionはEncoder側の  K, V を使用していることに注意)。

input/output embeddingやブロックの外のlinearはここでは説明していませんので悪しからず。Masked Multi-Head Attentionもここでは説明しませんでしたが、これは予測時に未来の情報がリークしないようにマスク処理が追加されたMulti-Head Attentionです。

付録

Positional Encodingのコード

import numpy as np
from numpy.typing import NDArray
import matplotlib.pyplot as plt
import seaborn as sns


def get_positional_encoding(pos: int, dmodel: int) -> NDArray[np.floating]:
    """Get positional encoding.

    Parameters
    ----------
    pos : int
        position
    dmodel : int
        dimension of the feature

    Returns
    -------
    pe : (dmodel,) NDArray[np.floating]
        positional encoding
    """
    pe = np.zeros(dmodel)
    pe_odd = np.arange(0, dmodel, 2)
    pe_even = np.arange(1, dmodel, 2)
    pe[pe_odd] = np.cos(pos / (10000 ** (pe_odd / dmodel)))
    pe[pe_even] = np.sin(pos / (10000 ** (pe_even / dmodel)))
    return pe


pe = np.vstack([get_positional_encoding(i, 100) for i in range(100)])
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
sns.heatmap(pe, ax=ax, cmap="viridis")
ax.set_ylabel("position", fontsize=16)
ax.set_xlabel("dimension", fontsize=16)
ax.set_title("position encoding", fontsize=16)
plt.show()

参考文献

  1. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L. & Polosukhin, I. (2017). Attention Is All You Need. arXiv. https://doi.org/10.48550/arxiv.1706.03762
  2. He, K., Zhang, X., Ren, S. & Sun, J. (2015). Deep Residual Learning for Image Recognition. arXiv. https://doi.org/10.48550/arxiv.1512.03385
  3. Wu, Y. & He, K. (2018). Group Normalization. arXiv. https://doi.org/10.48550/arxiv.1803.08494

特徴量重要度を理解したい

特徴量重要度を理解したい記事です。

まさか、どういうものか理解せずに使ってないですよね?

特徴量重要度とは

特徴量重要度は、名前の通り特徴量の予測に対する重要さの指標です。ファイナンス機械学習では以下の3つが紹介されています。

平均不純度変少量(Mean Decrease Impurity; MDI)

ランダムフォレストのようなツリーベース手法特有の重要度です。ジニ重要度(gini importance)とも呼ばれます。式で確認していきます。

決定木における葉を除く任意のノードを  t として、ノード  t の左右の子ノードをそれぞれ  t^{left} t^{right} と置きます。この時の不純度減少量は以下のように定義されます。

\begin{equation} \Delta_\mathcal{I}(t) = \text{Impurity}(t) - \frac{N (t^{left})}{N (t)} \text{Impurity} (t^{left}) - \frac{N (t^{right})}{N (t)} \text{Impurity} (t^{right}) \end{equation}

 \text{Impurity} は不純度、 N(t) はノード  t のサンプルサイズを表しています。データを分割することで不純度がどの程度減少したかを表しているのがわかると思います。

なお、不純度の定義については記事の最後に記載しました。

これを葉を除く全ノードについて計算し、合計が1になるように正規化したものがMDIです。

高速に計算できるのが利点ですが、予測に全く寄与しない場合でも値を持ってしまうというわかりにくさもあります。

平均正解率減少量(Mean Decrease Accuracy; MDA)

MDIは学習済みモデルを用いたアウトオブサンプルに対する重要度です。以下のように計算します。

  1. 学習済みモデルでアウトオブサンプルの予測についてパフォーマンスを算出する(パフォーマンスはAccuracyやRecallなど目的に応じて決める)
  2. 特徴量を1つ選びランダムシャッフルし、アウトオブサンプルの予測についてパフォーマンスを算出する。
  3. 特徴量重要度としてパフォーマンスの低下量を計算する
  4. 2、3をすべての特徴量について完了するまで繰り返し

どんな分類器でも使えるのが利点です。ただし、予測を繰り返すため低速です。

単一特徴量重要度(Single Feature Importance; SFI)

これは名前の通り、1つの特徴量でアウトオブサンプルでのパフォーマンスを計算する方法です。

代替効果(後述します)が生じないのがメリットです。 また、どんな分類器でも使えます。ただし、予測を繰り返すため低速です(単一特徴量を用いるという意味では高速かもしれませんが)。

留意すべきこと

ファイナンス機械学習では色々な留意事項が掛かれていますが、特に重要そうな代替効果についてまとめておきます。

代替効果(Substitution effect)

重回帰における多重共線性(通称マルチコ。Multicollinearity)のような効果です。特徴量間に相関があると問題になります。

相関のある2つの特徴量  a, b を考えましょう。極端な例として、相関1の場合を思い浮かべると簡単かと思います。

 MDI では上記二つの特徴量は同確率でランダムに選択されるので、重要度が分配される(半減する)ということがわかるでしょう。

 MDA でも代替効果は問題で、上記の例で例えば  b がランダムシャッフルされた場合、 a b の予測力を完全に補ってしまいます(相関1なので)。つまり、 b は全く重要でないと判断されます。そして逆もまた然りです。上記の例では、せめて「どちらかが重要」もしくは「どちらも重要」程度のことがわかれば特徴選択の役に立ちますが、MDAでは「どちらも重要でない」 と判断する可能性があるということになります。

SFIは特徴量間の相関を考慮しない。相乗効果も考慮しない。

SFIは特徴量1つずつに着目するため、特徴量間の相関を考慮しません。これは、代替効果のことを考えなくて良いという点で有用です。

一方で、特徴量の組み合わせが予測に重要な効果を持っている場合考慮できないという欠点があります。

代替効果に対する対策

特徴量の相関が問題になるなら直交基底に落とせば良いので、主成分分析(PCA)を行えば良いです。

そんなことかと思うかもしれませんが、PCAはモデルがオーバーフィッティングを確かめる強力な手段になっています。

PCAによって固有値の大きさから、主成分の分散説明力や主成分に対する各特徴の寄与度が得られるわけですが、上でまとめた各特徴量重要度のランキングと、PCAからわかる特徴量のランキングが一致していればオーバーフィッティングをしていないという証拠が得られたことになります。

Appendix

不純度(Impurity)

ジニ不純度  G(k) は以下のように定義されます。

\begin{equation} G(k) = \sum_{i=1}^{n} p(i)\left\{1 - p(i)\right\} \end{equation}

ここで、 nはターゲットラベルの数、 p(i) はノード  k におけるターゲットラベル  i の割合を表しています。

参考文献

アンサンブル法を理解したい #1

まだまだ使う機会の多いアンサンブル法ですが、理論とアルゴリズムからちゃんと追いかけられている人ってどの程度いるのでしょうか?

データサイエンティスト()にならないためにちゃんと理解してきたい記事です。

今回はバギングについてまとめます。

問題設定

以下の関係が成り立つと仮定する。

\begin{align} y_i &= f(x_i) + \varepsilon_i \\ E[\varepsilon_i] &= 0 \\ E[\varepsilon_i^2] &= \sigma_{\varepsilon}^2 \\ \end{align}

ここで  \displaystyle \left\{x_i \right\}_{i=1}^{n} は説明変数、 \displaystyle \left\{y_i \right\}_{i=1}^{n} は目的変数、 \displaystyle \varepsilon はホワイトノイズである。

この  f(x) をうまく近似する関数  \hat f(x) を推定したい。つまり、MSE(Mean Square Error; 平均二乗誤差)の期待値

\begin{align} E\left[\left\{y_i - \hat f(x_i)\right\}^2\right] \end{align}

を最小化するような  \hat f(x) を推定したいということになります。そしてMSEはバイアス、バリアンス、ノイズに分解できることが知られています。

\begin{align} E\left[\left\{y_i - \hat f(x_i)\right\}^2\right]&= \left( f(x_i) - E\left[\hat f(x_i)\right] \right)^2+ E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] + \sigma_i^2 \end{align}

右辺の左から順にバイアス、バリアンス、ノイズです。式変形は本記事の最後に記載しています。バイアスは対象モデルの真のモデルからの乖離度、バリアンスはモデルの安定性、ノイズはモデルで説明できない誤差を表しています。機械学習モデルに限った話ではないですが、理想のモデルに近く(バイアスの最小化)、ノイズをシグナルと勘違いしない(バリアンスの最小化)ようなモデルを作りたいということです。

アンサンブル法とは

アンサンブル法は複数のモデル出力の期待値を取ることで、個別モデルより良い性能を得ようとする手法の総称です。バギングとブースティングが有名です。冒頭で述べた通り、今回はバギングに関するまとめです。

バギングとは

複数のモデルの単純平均をアンサンブルの出力とする方法です。予測がカテゴリカル変数の場合は多数決になります。

バギングはバリアンスを小さくする

バギングはバリアンスを小さくすることが知られています。数式で証明できるので確認してみましょう。

バリアンスは以下の式で与えられました。

\begin{equation} \left( f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \end{equation}

ここで、バギング予測部分(期待値をとっているところ)を、各モデル  \left\{\hat f_{m}\right\}_{m=1}^{M} で明示的に書き下し、真のモデルと各モデルとの誤差をそれぞれ  \varepsilon_m と置きます。

\begin{align} \left( f(x_i) - E\left[\hat f(x_i)\right] \right)^2 &= \left( f(x_i) - \frac{1}{M} \sum_{m=1}^{M} \hat f_m (x_i) \right)^2 \\ &= \left( \frac{1}{M} \sum_{m=1}^{M} \varepsilon (x_i) \right)^2 \end{align}

さて、これでバギングによる平均二乗和誤差が得られたわけですが、個別モデルの平均二乗誤差をすべてのモデルについて平均した

\begin{equation} \frac{1}{M} \sum_{m=1}^{M} E\left[\varepsilon _m (x_i)^2\right] \end{equation}

よりバギング誤差の方が小さくなることが確認できれば、「バギングはバリアンスを小さくする」ということが言えます。これは、Jensenの不等式を利用することで簡単に証明することができ

\begin{align} \left( \frac{1}{M} \sum_{m=1}^{M} \varepsilon (x_i) \right)^2 \le \frac{1}{M} \sum_{m=1}^{M} E\left[\varepsilon _m (x_i)^2\right] \end{align}

が成り立ちます。等号が成り立つのは、誤差間に相関が無いときです。なお、証明は記事の最後に記載しました。

ちなみに、上記ではモデルごとに重複ありのランダムサンプルを用いた場合の議論でしたが、決定木モデルごとに重複なしのランダムサンプリングを行うようにしたのがランダムフォレストです。重複なしにすることで、モデル間の相関を下げる効果が期待されます。

Appendix

バイアス・バリアンス分解

 \hat f(x) の期待値  E \left[ \hat f(x) \right] を導入することで、MSEを分解します。

\begin{align} E\left[\left\{y_i - \hat f(x_i)\right\}^2\right] &= E\left[\left\{y_i - E\left[\hat f(x_i)\right] - \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right) \right\}^2\right] \\ &= E\left[\left( y_i - E\left[\hat f(x_i)\right] \right)^2 \right] \\ &\quad - 2 E\left[ \left( y_i - E\left[\hat f(x_i)\right] \right) \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right) \right] \\ &\quad + E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] \end{align}

ここで、観測時点で値が確定するものについて期待値  E の外に出すと

\begin{align} &= \left( y_i - E\left[\hat f(x_i)\right] \right)^2 \\ &\quad - 2 \left( y_i - E\left[\hat f(x_i)\right] \right) \left(E\left[ \hat f(x_i) \right] - E\left[\hat f(x_i)\right] \right) \\ &\quad + E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] \\ &= \left( y_i - E\left[\hat f(x_i)\right] \right)^2+ E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] \\ &= \left( f(x)_i + \varepsilon_i - E\left[\hat f(x_i)\right] \right)^2+ E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] \label{bb1}\\ &= \left( f(x)_i - E\left[\hat f(x_i)\right] \right)^2+ E\left[ \left(\hat f(x_i) - E\left[\hat f(x_i)\right] \right)^2 \right] + \sigma_i^2 \label{bb2} \end{align}

\eqref{bb1}から\eqref{bb2}のところだけ少し式変形を省略しましたが、ここまで式を追いかけられるなら簡単にわかるかと思います。

バギングによるバリアンス削減

Jensenの不等式は、凸関数  f に対して以下が成り立つというものです。

\begin{equation} f\left(\sum_{i=1}^M \lambda_i x_i \right) \le \sum_{i=1}^M \lambda_i f(x_i) \end{equation}

ここで、 \lambda_i = \frac{1}{M} x_i = \varepsilon_m (x_i) f(x_i) = x_i^2 と置けば

\begin{equation} \left( \frac{1}{M} \sum_{m=1}^{M} \varepsilon (x_i) \right)^2 \le \frac{1}{M} \sum_{m=1}^{M} E\left[\varepsilon _m (x_i)^2\right] \end{equation}

となります。

参考文献

インバランスバーを理解したい

 ファイナンス機械学習で、非定常時系列の前処理・解析の勉強をしているんですが再序盤である第2章で躓いてしまいました。読んで理解するだけなら難しくありませんが、インバランスバーを実装してみるとここで躓く人も多いのではないでしょうか。

 第2章3節ではティックからバーを生成する手法が列挙されているのですが、いくつか実用に耐えない問題があります。今回はその問題と対策についてまとめていきます。

情報ドリブンバー

 チャートでバーと言えば〇分足のような時間で区切ったタイムバーが一般的ですが、タイムバーはリターン分布(各時刻における価格変化分布)の統計的性質が良くないとされています。機械学習モデルではデータが独立同分布(i.i.d)で生成されていることを前提としているため、タイムバーを使用して学習・推論することは適していないといえます。2章ではさまざまなバーを紹介してくれていますが、今回はその中でも、定義と実装の両方に問題があると思われるインバランスバーについてまとめます。

ティックインバランスバー

 ティックインバランスバーは、ティックの不均衡度合いが期待値を超える度にバーを更新するバーです。書籍で以下のように定義されています。

 \displaystyle \left\{(p_t, v_t)\right\}_{t=1,\ldots,T} を時刻  \displaystyle t において、価格  \displaystyle p_t、出来高  v_tのティック系列とする。この時、 \displaystyle b_t \in \left\{-1, 1\right\}である以下の系列を定義する。

$$ \begin{equation} b_t = \left\{ \begin{aligned} b_{t-1} \quad &\text{if } \Delta p_t = 0 \\ \frac{|\Delta p_t|}{\Delta p_t} \quad &\text{if } \Delta p_t \ne 0 \end{aligned} \right. \end{equation} $$

そして、ティックインバランスを以下のように定義する。

$$ \begin{equation} \theta_T =\sum_{t=1}^{T} b_t \end{equation} $$

そして、このティックインバランスが閾値を超えたとき(つまり、ティック不均衡となったとき)に、バーを更新します。

 閾値の決め方は、筆者によると以下のように決定できるとされています。

$$ \begin{equation} T^* = \text{argmin}_T \left\{|\theta_T| \ge\mathbb{E_0}[T]|2\mathrm{P}[b_t=1]-1|\right\} \end{equation} $$

しかしこの式には大きな問題点があります。

閾値が大きくなり続ける

 実は、この式の通りに実装を行うとインバランスバーの閾値が大きくなり続け、最終的にバーが更新されなくなってしまいます。これは、現在の閾値を超えたらバーを更新するということと、定義式の \displaystyle 2P[b_t = 1 ] - 1の最大値が  1 であることに起因します。閾値  \displaystyle E[\theta] を求めるときに  \displaystyle 2P[b_t = 1] -1 の上限が決まっているので  \displaystyle T の方で調整するしかなく、「現在の閾値を超えたら」バーが更新されるので、要求される  \displaystyle T はどんどん大きくなっていきます(細かいことを言えば原因はこれだけじゃないですが...)。

対策

2022-09-04 追記
ツイッターにて、間違いをご指摘いただいたので修正しています [2] 。

 閾値の定義式を嚙み砕いて言ってしまえば、「ラベルの偏り率が \displaystyle 2P[b_t = 1 ] - 1のとき、ティック数 Tが経過したらどの程度偏りが蓄積しているか」ということなので、閾値の定義は偏りの正常とされる許容範囲を示しているに過ぎません。

 「インバランスを検知したい」という基本に立ち返るならば、固定閾値をこえたらバーを更新し、都度 b_tをリセットすれば良さそうです。式で書くと以下のようになります。

$$ \begin{equation} \begin{aligned} \theta_T &=\sum_{t=1}^{T} b_t \\ T^* &= \text{argmin}_T \left\{|\theta_T| \ge \text{const.} \right\} \end{aligned} \end{equation} $$

 \text{const.}をどの程度の値に固定すれば良いのかは決め内でも良いと思いますし、動的に変化させたいなら b_t系列について偏りの移動平均を計算してみても良いと思います。この辺りが妥協点じゃないでしょうか?

実際に本対策法でFTXJPのBTC-PERPのリターン分布を計算してみたところ以下のようになりました。1か月分のデータに対してバーの数が1500~2000程度になるような閾値で計算しています。閾値次第な部分もあると思うのであまり参考にならないかもしれないですが、発散さず計算できたことがわかります。

 閾値によって分布は変わるのでちゃんとした比較にはなりませんが、ボリュームバーやドルバーは、ティックバーよりもやや緩やかなで、より対称な分布になっている気がします。

さいごに

  • これはあくまで私の見解です。間違いやもっと良いアイデアがあるよという方はコメントいただけると嬉しいです。
  • インバランスバーの次にランバーというものも紹介されていますが、定義を見る限り同じ問題が発生しそうです。

参考

  1. https://twitter.com/kodojibtc/status/1566034056542027776?s=20&t=s_JJkLKNgQXaZnXwRYFwUw

強化学習を理解したい #2

前回

hotoke-x.hatenablog.com

の続きで、今回はモデルフリー手法の話です。 モデルベース手法では  \pi (a|s) T(s'|s, a) が定義可能である前提でした。

モデルフリーの問題では直接これらを定義できない為に、蓄えたデータ(経験)から推定する必要があります。 また、経験を活用して行動しなければ、いつまでたっても報酬は得られません。

このようにモデルフリータスクでは経験と活用のバランスを考える必要があります。

また、どの範囲の経験をフィードバックするか、価値評価改善と戦略改善のどちらへフィードバックするかなどの問題もあります。

ここでは、

  • 経験と活用のバランスをとるEpsilon-Greedy法
  • フィードバック範囲を調整する方法
  • 価値評価と戦略のどちらを改善するか

についてまとめます。

経験と活用のバランスをとるEpsilon-Greedy法

Epsilon-Greedy法は至って単純で、確率  \mathcal{\epsilon} でエージェントをランダムに行動させるだけです。

これによって、戦略に従って行動しているだけでは取りにくい行動(選択確率の低い行動)も一定確率で行うようになり、状態の遷移と遷移によって得られる報酬を探索することができます。

フィードバック範囲を調整する方法

次に考えるべきは、経験を使った学習です。

経験を積みつつもエージェントは時々刻々と行動し続けているので、フィードバックすべき経験の範囲を考えるのも、強化学習におけるひとつのポイントとなります。

典型的には以下の方法が有名です。

  • TD法(Temporal Difference)
    • 行動するたびに逐次的にフィードバックします。
  • Monte Carlo法
    • エピソード終了時にフィードバックします。

上記の他に中間的な方法もありますが、結局TD法、Monte Carlo法含め一つの式で書けるので深く考える必要はないです。

では、TD法とMonte Carlo法における価値更新の式を見てみましょう

 \displaystyle
V(s _t) \leftarrow V(s _t) + \alpha (r _{t+1} + \gamma V(s _{t+1}) - V(s _t)) \\
V(s _t) \leftarrow V(s _t) + \alpha \left((r _{t+1} + \gamma r _{t+2} + \gamma^ 2 r _{t+3} + \ldots + \gamma ^ {T-t-1} r _ {T-t}) - V(s _ t)\right)

上段がTD法、下段がMonte Carlo法です。1 step先の価値を予測して誤差の大きさに応じて学習率  \alpha で価値関数を更新するのがTD法、 エピソード終了まで待って価値関数を更新するのがMonte Carlo法というわけですね。

ただ、フィードバック範囲が異なるだけなら1つの式で書けそうだという直感はあっていて、実際以下の式(TD ( \lambda) 法)が知られています。

 
G^\lambda_t = (1 - \lambda)\sum_{n=1}^{T-t-1} \lambda^{n-1}G_t^{(n)} + \lambda^{T-t-1}G_t^{(T-t)}

ここで、 \lambda \in [0, 1] であり、  \lambda = 0 のときはTD法、 \lambda = 1 のときはMonte Carlo法と等価になります。

 \lambda の値次第でフィードバックタイミングを調整可能ということです。

価値評価と戦略のどちらを改善するか

そして最後に、価値評価と戦略のどちらを改善するかですが、言ってしまえば扱いたい問題がValueベースか、Policyベースかという問題です。

内容まではまとめませんが、以下のような手法が知られています。

  • 価値評価の改善(Valueベース)
    • Q-learning
      • Deep learningを利用したQ-learningが、かの有名なDeep Q Network (DQN) です。
    • 戦略関係なしに行動するので、この行動基準をoff-policyと呼びます。
  • 戦略の改善(Policyベース)
    • SARSA (State-Action-Reward-State-Action)
    • 戦略の通りに行動するので、この行動基準をon-policyと呼びます。
    • 価値評価と戦略評価に同じQ-table(価値関数)を使っているのが特徴

SARSAの特徴として「価値評価と戦略評価に同じQ-table(価値関数)を使っているのが特徴」と挙げましたが、価値評価と戦略評価を別々に行うActor-Critic法という手法もあります。

Policy iterationでは、価値の推定と戦略の評価が別々に行われていたことから、価値と戦略を別々に考えることにも違和感はないと思います。

まとめ

というわけで、モデルフリー手法についても、ちょっと詳しい目次程度にまとめました。

モデルフリーとはいえ、まだ扱う状態が離散で、価値がテーブルに収まっているので真のモデルフリーではない気もします。

以降では、状態や価値を関数で近似してしまうことを考えます。

データを入力したら状態を出力する関数、状態を入力したら価値を出力する関数を求めるということですね。

しかしモデルがわかっていませんから当然関数の形も未知です。

その為、表現力の高いDeep learningモデル等を利用しているのが昨今の強化学習になります。

強化学習業界を追いかけるのは容易ではないと思いますが、これくらいまで勉強しておけば業界全体の流れは汲み取れたのではないでしょうか。

後は実装あるのみですね。

参考書籍

Pythonで学ぶ強化学習

強化学習を理解したい #1

仕事でDeep Learningを使う機会がちょいちょい出てきたので、tensorflowやPyTorchの練習も兼ねて強化学習を勉強中。

手法が色々あってよくわからなくなってきたので記事にまとめます。

  • 全体を俯瞰することに努めます。
  • 正確性よりわかりやすさを優先します。
  • 基本的な統計・機械学習の知識があることを前提にまとめます。

問題設定

目的

状態  \displaystyle s に対して行動  \displaystyle a を起こし、得られる報酬の総和  \displaystyle \sum r を最大化する。

方法

価値評価と戦略を学習する

  • 価値評価
    • 価値の算出方法のこと。
    • 価値とは、得られる報酬の総和の見積値のこと(期待報酬*とも呼ぶ)。
  • 戦略
    • 状態を受け取り行動を出力する関数のこと(いわゆるモデルのこと)

*期待報酬とは言うものの、実際は「割引現在価値」の総和 で、いわゆる統計的な期待値計算とは異なります。 現在時点から将来にわたって得られる報酬を見積もるわけですから、 未来の不確実性を考慮して、得られる(と期待されている)報酬に対して割引率を掛けて和をとります。

登場用語

大まかな問題設定がわかったところで、強化学習に登場する基本的な用語をまとめておきます。

  • 環境とエージェント
    • 問題設定を規定するものと、問題設定に従って行動する主体のことです。
    • 自動運転で例えるなら、道路と車の関係になります。
  • 状態  \displaystyle s
  • 行動  \displaystyle a
  • 遷移確率  \displaystyle T
    • transitionのTです。
    • 行動が決まれば遷移先が決まる...というわけでは無いので分けて考える必要があります。
      • 歩こうと思っても躓いたり転んだりすることだってあるということです。
  • 即時報酬  \displaystyle R
    • rewardの頭文字です。
    • 状態が遷移したときに実際に得られる報酬(行動1回分の報酬の)です。

「価値評価」と「戦略」

強化学習では「価値評価」と「戦略」を学習するわけですが、「価値」とは何でしょうか?

ここの理解を疎かにするとこの後混乱するだけなので、価値とは何かについてまとめます。

結局のところ、報酬によって定義されるということを知っておくとこの後がわかりやすくなると思います。

強化学習における「価値」

現時刻を  \displaystyle t として、戦略  \displaystyle \pi に基づいた行動  \displaystyle a_t によって 現状態  \displaystyle s_t から 次の状態  \displaystyle s_{t+1}へ遷移することを考えます。

この時、 \displaystyle s_t から  \displaystyle \pi に基づいて行動して得られる価値は、以下のように定義されます。

 \displaystyle
V_\pi(s_t) = \text{E}_{\pi}\left[r_{t+1} + \gamma V_{\pi}(s_{t+1}) \right]

なぜこの式で良いのかというと、再帰的に定義された式を書き下してみると

 \displaystyle
V_\pi(s_t) = \text{E}_{\pi}\left[r_{t+1} + \gamma \text{E}_{\pi}\left[r_{t+2} + \gamma V_{\pi}(s_{t+2}) \right] \right]

となり、行動  a_t によって得られる即時報酬  r_{t+1} と、その後に得られるであろう価値の期待値を取ったものになっています( \gamma は上で説明した割引率です)。

これが、「 s_t から  \pi に基づいて行動して得られる価値」なわけですが、式を見てわかる通り、 「価値」とは割引率を考慮した報酬の期待値です。つまり、状態が遷移する際の報酬が決まっていれば、自然と定まるものになっています。

ValueベースとPolicyベース

「価値」が何かわかったところで、次はエージェントの行動指針である「戦略」についてまとめます。

戦略には大きく分けてValueベースとPolicyベースの2つが存在します。

  • Valueベース
    • 各状態で常に最大価値の行動を取る手法のこと。
  • Policyベース
    • 各状態で確率的に価値の高い行動を取る手法のこと。

要は決定論的か確率的かという違いですね。 言葉にすると簡単ですが、Policyベースの方は最適化のプロセスが少々複雑です。

Value iteration と Policy iterateion

さて。ValueベースとPolicyベースについては一言でまとめてしまいましたが、それぞれの戦略が機能するのは価値を正しく見積もれる場合のみです。

最初から正しい価値など知るはずもありませんから、ここからが強化学習の本領となります。

では、価値の定義を思い出してみましょう。

 \displaystyle
V_\pi(s_t) = \text{E}_{\pi}\left[r_{t+1} + \gamma V_{\pi}(s_{t+1}) \right]

先ほどは説明のために「期待値」という言葉で  \text{E} の表記をそのままにして置きましたが、展開すると

 \displaystyle
V_{\pi} (s_t) = \sum_{a} \pi(a|s) \sum_{s'} T(s'|s, a)\left\{R(s, s`) + \gamma V_{\pi} (s')\right\}

と書き直すことができます。ここで、新しく登場した表記について書いておくと、

  •  \pi(a|s)
    • 状態  s で行動  a を選択する確率
  •  T(s'|s, a)
    • 状態  s で行動  a を選択したときに状態  s' へ遷移する確率(遷移確率)
  •  R(s, s')
    • 状態  s から  s' へ遷移した際に得られる即時報酬

となっています。ちなみにこの式をBellman equationといいます。この式において、ValueベースとPolicyベースは以下のような関係になっています。

f:id:hotoke-X:20220206041638p:plain:w600

Bellman equationと戦略の関係がわかったところで、それぞれの戦略の最適化方法についてまとめます。

Value iteration

まず、Valueベース手法の最適化についてですが、Bellman equationでは、状態  s から  s' へ遷移したときの価値を定義しているのでした。

つまり  \pi (a|s) T(s'|s, a) さえ定義されていれば全ての遷移について価値の計算が可能ということです。

そこで、Valueベース手法は以下の式を繰り返し計算するだけで最適化が完了します。

 \displaystyle
V_{i+1} (s) \overset{\text{def}}{=} \sum_{s'} T(s'|s, a)\left\{R(s, s') + \gamma V_{i} (s')\right\}

ここで、 i は最適化の繰り返し数を表しています。

つまり、上式で全状態について一通り価値を計算し、それを次の最適化計算へ利用しているだけです。それぞれの状態について価値の更新が行われなくなったら計算終了とします。

価値を繰り返し計算するので、この最適化プロセスをValue iterationと呼びます。

Policy iteration

次に、Policyベース手法の最適化についてです。

同じBellman equationが基礎となっていますから、Value iterationと同様繰り返し計算で最適化が完了しますが、計算のプロセスは少々複雑です。

Policyベースでは、得られる価値を最大化するように戦略を更新する必要がありますがエージェントは戦略に従って確率的に行動します。

確率的に行動するのに最大化?と思う方も居るかもしれませんが、要は戦略に従った場合に得られる価値の期待値を最大化すれば良いということです。

しかし、困ったことにそもそも正確な価値がわかっていないのでした。そこで、Policyベースの最適化 (Policy iteration)では

  1. 価値の推定
  2. 戦略の評価

の2段構えで繰り返しの最適化が行われます。価値と戦略の両方が交互に最適化されるわけですね。

価値の推定はほとんどValue iterationと同じ形の式ですが、最大価値を求めるのではなく期待値を求めるのが異なる点です。具体的には以下のようになります。

 \displaystyle
V_{\pi} (s_t) = \sum_{a} \pi(a|s) \sum_{s'} T(s'|s, a)\left\{R(s, s`) + \gamma V_{\pi} (s')\right\}

つまり、素直にBellman equationの計算をしているだけですね。あとはValue iterationと同じ要領で価値が収束するまで計算を繰り返せば価値の推定は完了します。

価値の推定ができたら、戦略の評価を行います。戦略の評価では、それぞれの状態  s について以下のアルゴリズムで戦略を更新します。

1. 現在の戦略に従った場合最も価値の高い行動を取得する( a_{\text{policy}}とおく)
2. 推定価値に基づいて各行動の価値を計算し、最大価値の行動を取得する( a_{\text{best}}とおく)
3.  a_ {\text {policy}} a_ {\text {best}} が異なっていたら a_ {\text {best}}を選ぶように戦略を更新する

戦略の更新が終わったら再び価値の推定に戻り、戦略の更新が行われなくなるまで以下繰り返しです。

まとめ

以上が強化学習の基礎となる考え方でした。

結局どの手法でも価値の計算が必要で、その価値は報酬によって規定されていることを理解すると全体像が見えてくると思います。

また、今回まとめた方法は  \pi (a|s) T(s'|s, a) が定義可能である前提でした。

これらが定義可能な問題を、強化学習モデルが明らかであるということからモデルベースと呼びます。

これらが定義できないモデルフリーの問題にどう対処するのかについては今回触れませんが、基本となる考え方は同じです。

モデルベースの問題をちゃんと理解していれば理解に困ることは少ないでしょう(多分)。

参考書籍

Pythonで学ぶ強化学習

Big SurでのBUILD FAILED (OS X 11.5.1 using python-build 20180424)を解決したい

macOS Big Sur (Version 11.5.1)にてpyenvで3.7系のpythonをインストールしようとすると

BUILD FAILED (OS X 11.5.1 using python-build 20180424)

のエラーが発生した。

結論から言うと、brewを経由せずにgithubから直接pyenvを取ってきたら解決した。
既にgithubリポジトリから直接インストールしている人は

cd $(pyenv root)
git pull

で最新版のpyenvを取ってくれば良い。

ただ、これだけでは3.6.8が結局入らず、今度は

BUILD FAILED (OS X 11.5.1 using python-build 2.0.6)

のエラーが発生。

pyenvのissueに解決策があったので、コマンドを以下にまとめておいた。


3.6.8のインストール

CFLAGS="-I$(brew --prefix openssl)/include -I$(brew --prefix bzip2)/include -I$(brew --prefix readline)/include -I$(xcrun --show-sdk-path)/usr/include" LDFLAGS="-L$(brew --prefix openssl)/lib -L$(brew --prefix readline)/lib -L$(brew --prefix zlib)/lib -L$(brew --prefix bzip2)/lib" pyenv install --patch 3.6.8 < <(curl -sSL https://github.com/python/cpython/commit/8ea6353.patch\?full_index\=1)

3.6.12のインストール

CFLAGS="-I$(brew --prefix openssl)/include -I$(brew --prefix bzip2)/include -I$(brew --prefix readline)/include -I$(xcrun --show-sdk-path)/usr/include" LDFLAGS="-L$(brew --prefix openssl)/lib -L$(brew --prefix readline)/lib -L$(brew --prefix zlib)/lib -L$(brew --prefix bzip2)/lib" pyenv install --patch 3.6.12 < <(curl -sSL https://github.com/python/cpython/commit/8ea6353.patch\?full_index\=1)

3.7.9のインストール

CFLAGS="-I$(brew --prefix openssl)/include -I$(brew --prefix bzip2)/include -I$(brew --prefix readline)/include -I$(xcrun --show-sdk-path)/usr/include" LDFLAGS="-L$(brew --prefix openssl)/lib -L$(brew --prefix readline)/lib -L$(brew --prefix zlib)/lib -L$(brew --prefix bzip2)/lib" pyenv install 3.7.9

3.8.0のインストール

pyenv install --patch 3.8.0 < (curl -sSL https://github.com/python/cpython/commit/8ea6353.patch\?full_index\=1 | psub)

3.8.6のインストール

LDFLAGS="-L$(xcrun --show-sdk-path)/usr/lib" pyenv install 3.8.6

とりあえずこれで動いているものの、気持ちが悪い。

参考

MLOpsを理解したい #1

会社で「MLOpsについて勉強しよう」となったので、最近Googleより発表されたVertex AIについて勉強していきたい記事です。
初心者すぎるので、MLOpsとは何かというところから始めていきます。

MLOpsとは

機械学習オペレーション (machine learning operations) の略。
これについては色々な方がまとめてくれている。
今回参照させていただいたものはこちら。

speakerdeck.com

大事だと思ったこと抜粋

  • 「機械学習エンジニアは、自分が開発した機械学習アプリケーションについて、デプロイまで責任を持つ」

あまりできていないと思ったこと

  • 実験データのバージョン管理
  • 入力データのバリデーション
  • 複数の機械学習モデル利用における依存関係を考慮した設計

Julia tips #11: REPL起動と同時にプロジェクト環境をアクティベートする

~/.julia/config/startup.jl

using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end

と書いておくだけ。この設定はVSCodeのJulia extensionで自動起動されるREPLでも有効。

参考

K-SVDアルゴリズムを理解したい

辞書学習アルゴリズムであるK-SVDアルゴリズムを理解したい記事です。

また、過去記事

hotoke-x.hatenablog.com

の続きでもあります。

辞書学習の問題設定

許容誤差を \epsilonとして、以下の問題設定で辞書 Aとスパース表現ベクトル \boldsymbol{x}を推定する。 (許容誤差のことを書籍ではモデル誤差と呼んでいるが、誤差が既知な状況はほぼないのでここでは許容誤差と呼ぶことにする。)

$$ \begin{align} \min _{\boldsymbol{A}, \left\{\boldsymbol{x}_{i}\right\}_{i=1}^{M}} \sum_{i=1}^{M}||\boldsymbol{x}_{i}||_{0} \quad \text { subject to } \quad ||\boldsymbol{y}_{i}-\boldsymbol{A} \boldsymbol{x}_{i}||_{2} \leq \epsilon, 1 \leq i \leq M \end{align} $$

また、ペナルティとスパース性を入れ替えて

$$ \begin{align} \min _{\boldsymbol{A}, \left\{\boldsymbol{x}_{i}\right\}_{i=1}^{M}} \sum_{i=1}^{M}||\boldsymbol{y}_{i}-\boldsymbol{Ax}_{i}||_{2}^{2} \quad \text { subject to } \quad||\boldsymbol{x}_{i}||_{0} \leq k_{0}, 1 \leq i \leq M \end{align} $$

という誤差最小化問題を考える。

 \epsilon=0とすると、解の一意性が保証されるらしい(驚き)。

辞書とスパース表現ベクトルを同時に求める問題は行列分解とみることもできる。

K-SVDアルゴリズム

A中の j_0番目以外の列(辞書の列のことをアトムと呼ぶ)を固定し、 j_0番目に対応するアトムとそのアトムにかかる係数(スパース表現ベクトル)を更新する。そこで、誤差計算を以下のように分解する。

$$ \begin{align} ||\mathbf{Y}-\mathbf{A} \mathbf{X}||_{F}^{2} &=\left|\left|\mathbf{Y}-\sum_{j=1}^{m} \mathbf{a}_{j} \mathbf{x}_{j}^{\top}\right|\right|_{F}^{2} \\ &=\left|\left|\left(\mathbf{Y}-\sum_{j \neq j_{0}} \mathbf{a}_{j} \mathbf{x}_{j}^{\top}\right)-\mathbf{a}_{j_{0}} \mathbf{x}_{j_{0}}^{\top} \right|\right|_{F}^{2} \end{align} $$

 \mathbf{x}_j^{\top} \mathbf{X} j番目の行(列じゃないことに注意)とし、更新ステップで \mathbf{a}_{j_0} \mathbf{x}_{j_0}^{\top}の両方を更新する。 j_{0}番目以外は固定しているので、

$$ \begin{align} \mathbf{E}_{j_{0}} = \mathbf{Y}-\sum_{j \neq j_{0}} \mathbf{a}_{j} \mathbf{x}_{j}^{\top} \label{ksvd_error} \end{align} $$

は計算済みとしている。

ここで、\eqref{ksvd_error}を最小化する \mathbf{a}_{j_0} \mathbf{x}_{j_0}^{\top}は、 \mathbf{E}_{j_{0}}の低ランク近似(ランク1近似)とみることができる。

このことは一見するとわかりにくいが、 \mathbf{a}_{j_{0}} \mathbf{x}_{j_{0}}^{\top} \mathbf{Y}と同じサイズになっていることに気がつけば難しくない。仮に \mathbf{Y}_{j_{0}} = \mathbf{Y} - \mathbf{E}_{j_{0}}とすれば、 \mathbf{Y}_{j_{0}}  \mathbf{a}_{j_{0}} \mathbf{x}_{j_{0}}^{\top}で近似しているだけだからである。

これはSVDで求めることができるが、通常 \mathbf{x}_{j_0}^{\top}は密な解が求まってしまう。

そこで、  \mathbf{x}_{j_{0}}^{\top}の非ゼロ要素を用いて計算された \mathbf{E}_{j_{0}}の部分集合を考えることで、非ゼロ要素を固定することができる。

この計算を簡単にするために、制限作用素 \mathbf{P}_{j_{0}} \in \mathbb{R}^{M \times M_{j_{0}}}を定義し、

$$ \begin{align} \mathbf{E}_{j_{0}} \mathbf{P}_{j_{0}} \end{align} $$

として、不要な列を除去する。この制限作用素を用いて、 \mathbf{x}_{j_{0}}の非ゼロ要素だけを取り出したものを

$$ \begin{align} \left(\mathbf{x}_{j_{0}}^{R}\right)^{\mathrm{T}}=\mathbf{x}_{j_{0}}^{\mathrm{T}} \mathbf{P}_{j_{0}} \end{align} $$

と定義する。

この部分行列 \mathbf{E}_{j_{0}} \mathbf{P}_{j_{0}}に対してSVDによるランク1近似を適用すると、アトム \mathbf{a}_{j_{0}}と対応するスパース表現の要素 \mathbf{x}_{j_{0}}^{R}を同時に求めることができる。

すなわち

$$ \begin{align} \mathbf{E}_{j_{0}}^{R} = \mathbf{U} \boldsymbol{\Delta} \mathbf{V}^{\top} \end{align} $$

として、 \mathbf{a}_{j_0} = \mathbf{u}_{1},  \mathbf{x}_{j_{0}}^{R} = \Delta [1, 1] \mathbf{v}_{1}とする。ここで [1, 1]はインデックスを表していることに注意。

最後に以下の式で誤差を評価し、許容誤差未満になるまで反復を繰り返す。

$$ \begin{align} \left|\left| \mathbf{Y}-\mathbf{A}_{(k)} \mathbf{X}_{(k)} \right|\right|_{F}^{2} \le \epsilon \end{align} $$

また、書籍ではSVDではなくブロック座標降下法の考え方を用いた計算方法も紹介されているが、ここでは割愛する。

参考