Fire Engine

消防士→ITエンジニア→研究者

遺伝的アルゴリズムをGoで実装してみた

こんにちは!つるべーです!
最近は、進化計算と呼ばれるバイオミメティックな計算技法に興味を持っており、実装しながら勉強しています。
前回の記事では粒子群最適化(Particle Swarm Optimization: PSO)を実装しました
今回は、遺伝的アルゴリズムという生物の進化の過程を模倣して作られた最適解探索アルゴリズムをGoでスクラッチで書いてみました。

遺伝的アルゴリズム(Genetic Algorithm: GA)

GAとは何か?という問いに対する答えとして、以下の一文が端的にその特徴を表していると思います。

生物進化における遺伝と適者生存による自然淘汰の仕組みをソフトウェア的に模すことで複雑な問題に対する最適解を探索する手法
引用:「遺伝的アルゴリズム(Genetic Algorithm)を始めよう!」

GAでは、問題に対する解を個体の染色体、解の構成要素を遺伝子として表現し、複数の個体を進化させて、より環境に適応できる個体を見つけることで最適解を導いていきます。その過程をフローチャートとしてまとめたものが以下になります。

f:id:hirotsuru314:20190503224038p:plain

コード

github.com

eago(イーゴ)という進化計算パッケージを鋭意開発中です!
ちなみに、Pythonで書かれたdeapという素晴らしいライブラリがあり、遺伝的アルゴリズムを始めとした様々な進化計算の手法が実装されています。

最適化問題を解いてみる

今回はAckley functionと呼ばれる関数の最小化問題を解いてみたいと思います。数式およびグラフは以下の通りです。グラフは左側が3Dに描いたもので、右側が2次元に写像したもの(左のグラフを上から見たもの)です。

f:id:hirotsuru314:20190504145941p:plain

今回は n=2の条件を採用しています。
Ackley functionは多峰性関数であり、多くの局所解をもつため、最適化アルゴリズムベンチマーク関数として用いられることもあるようです。下のサイトに様々なベンチマーク関数が載っていて眺めるだけでも面白かったです。

Test functions for optimization - Wikipedia

グラフからわかるように、この関数は (x_1,x_2)=(0,0)で最小値 0を取るため、今回の探索の目的は、 (x_1,x_2)=(0,0)を見つけることになります。

サンプルコードは以下の通りです。

package main

import (
    "github.com/tsurubee/eago"
    "log"
    "math"
)

type Variables []float64

func (V Variables) Initialization() eago.Genome {
    return Variables(eago.InitFloatVector(2, 32, -32))
}

func (V Variables) Fitness() float64 {
    return Ackley(V)
}

func (V Variables) Mutation() {
    eago.AddNormalFloat(V, 0.5)
}

func (V Variables) Crossover(X eago.Genome) eago.Genome {
    return Variables(eago.BLXalpha(V, X.(Variables), 0.3))
}

func Ackley(X []float64) float64 {
    a, b, c, d := 20.0, 0.2, 2*math.Pi, float64(len(X))
    var s1, s2 float64
    for _, x := range X {
        s1 += x * x
        s2 += math.Cos(c * x)
    }
    return -a*math.Exp(-b*math.Sqrt(s1/d)) - math.Exp(s2/d) + a + math.Exp(1)
}

func main() {
    var v Variables
    ga := eago.NewGA(eago.GAConfig{
        PopulationSize: 30,
        NGenerations:   20,
        CrossoverRate:  0.8,
        MutationRate:   0.01,
    })
    if err := ga.Minimize(v); err != nil {
        log.Fatal(err)
    }
}

コードの詳細については最後に補足します。
上のコードをmain.goとして以下のように実行します。

$ go run main.go
Generation   1: Fitness=11.208 Solution=[-4.501 0.087]
Generation   2: Fitness=10.001 Solution=[-3.592 -0.833]
Generation   3: Fitness=9.165 Solution=[-3.841 -0.960]
Generation   4: Fitness=5.747 Solution=[-1.426 -0.235]
Generation   5: Fitness=4.258 Solution=[-0.824 0.380]
Generation   6: Fitness=2.797 Solution=[-0.843 -0.011]
Generation   7: Fitness=2.605 Solution=[-0.955 -0.030]
Generation   8: Fitness=2.605 Solution=[-0.955 -0.030]
Generation   9: Fitness=2.604 Solution=[-0.925 0.016]
Generation  10: Fitness=2.604 Solution=[-0.925 0.016]
Generation  11: Fitness=2.479 Solution=[-0.346 -0.054]
Generation  12: Fitness=1.216 Solution=[-0.184 -0.006]
Generation  13: Fitness=0.615 Solution=[-0.111 -0.001]
Generation  14: Fitness=0.501 Solution=[-0.096 -0.007]
Generation  15: Fitness=0.267 Solution=[-0.055 0.024]
Generation  16: Fitness=0.056 Solution=[0.010 -0.014]
Generation  17: Fitness=0.056 Solution=[0.010 -0.014]
Generation  18: Fitness=0.053 Solution=[-0.005 -0.015]
Generation  19: Fitness=0.031 Solution=[0.010 0.003]
Generation  20: Fitness=0.013 Solution=[-0.000 0.004]

出力結果を見てみると、世代数が大きくなるにつれて、Fitnessが 0、パラメータの解が (x_1,x_2)=(0,0)に収束しようとしているため、うまく探索できていることがわかります。

最後に、今回用いたコードとアルゴリズムの補足をしておきます。

実装とアルゴリズムの補足

Genomeインターフェース

実装では、下記のようなGenomeインターフェースを定義しており、上記のサンプルコードのVariablesはGenomeインターフェースを満たすように各メソッドを実装しています。

type Genome interface {
    Initialization() Genome
    Fitness() float64
    Mutation()
    Crossover(Genome) Genome
}

このように実装した理由は、eagoを使う側で自由に独自に実装した選択や交叉などの手法を組み込めるようにするためです。できればメジャーな手法については全てeago側に実装して、特殊なケース以外は全てeagoに実装されたメソッドを呼び出すだけで事足りるようにしていきたいと思っています。

実数値遺伝的アルゴリズム(Real coded GA)

今回の最適化では、解表現に浮動小数点数を用いており、実数値GAと呼ばれる手法を実装しました。
GAにおいては、各遺伝子を0か1で表現したバイナリエンコーディングが最も一般的な解表現だと思います。このような解の表現の仕方により、探索過程の各ステップ(選択・交叉・突然変異)で用いる手法も大きく異なってきます。
それぞれのステップにおいて、今回用いた手法を説明します。

選択

今回用いている選択の手法は、トーナメント法と呼ばれるものです。
任意のn個の個体を組にして、それらの中で最も適応度が高いものを選択するという手法になります。
他の手法として、ルーレット法やランキング法などがあります。

交叉

交叉にはブレンド交叉(BLX-α)という実数値GA特有の手法を用いました。
BLX-αについては下のサイトがわかりやすく解説されていました。

実数型GAに於ける交叉法の改良

親となる2個体によってできる長方形を一定数倍引き伸ばした長方形で囲まれる範囲内でランダムに次世代の個体が生成されます。

突然変異

突然変異として、個体の各遺伝子に任意の確率で一様乱数を加算するようなものを実装しました。一般的に突然変異は、局所解に陥ることを防ぐために取り入れられます。

さいごに

GAは、解の表現の仕方や、選択・交叉・突然変異など各遺伝子操作に様々な手法が提案されていて、何を使うかは問題を解く側に委ねられており、以前実装したPSOと比べて、非常に自由度が高いという印象を受けました。なので、GAの理論や実装をまだまだ深掘りしながら楽しめそうです!
今後は、処理の並列化や多目的最適化への拡張などを中心に勉強していきたいと思っています。

粒子群最適化(Particle Swarm Optimization: PSO)をGoで実装してみた

粒子群最適化とは、群知能による最適化手法の一種です。この手のバイオミメティクス(生物模倣)によるアプローチは、私の学生時代の専門である材料工学でも非常に盛んに研究されていましたが、データサイエンスでも応用されているのを知って興味を持ったので、実装しながら勉強しました。

粒子群最適化(Particle Swarm Optimazation: PSO

粒子群最適化とは、鳥や魚などの群れに見られる社会的行動のシミュレーションを基にモデル化されたヒューリスティックな最適解探索アルゴリズムです。

全体的な最適化の流れを説明します。 まず最初に探索空間内で、決められた個数の粒子を全て初期化します。初期化とは空間内の位置をランダムに決めることを指します。

それぞれの粒子は、目的関数により適合度(どれだけその解が優れているか)を算出できます。全ての粒子をある移動ルールを元に動かしながら、適合度が高い位置を探していくのがPSOの探索方法です。

ここでPSOにおいて重要な移動ルールについて説明します。

粒子 P_{i}の時刻 tから t+1の移動速度 \vec{v}_{i}(t+1)は以下の式で表されます。

f:id:hirotsuru314:20190409205526p:plain:w500

ここで重要なのが、グローバルベスト \vec{g}(t)とパーソナルベスト \vec{p}_{i}(t)です。
グローバルベストとは全ての粒子群で過去の探索も含めて最も適合度が高い位置、パーソナルベストとは特定の1粒子におけるこれまでの探索の中で最も適合度が高い位置を示します。これは繰り返しの探索の中で随時更新されていく変数です。
 Iは慣性係数と呼ばれる定数で、前の移動速度がどれだけ維持されるかを表します。 A_{g}および A_{p}は加速係数と呼ばれる定数です。

この数式からわかるように粒子群の中ではグローバルベストの情報が共有されていて、その方向に向かうように移動速度を調整しながら、より良い解の探索を繰り返します。これはあたかも餌を見つけた鳥が仲間にその情報を伝えて、仲間が餌のある位置に集まってくるようです。

作ったもの

github.com

『Evolutionary Algorithm implemented in Go』の略で名前を付けました。 進化計算全般に興味があるので、ここに実装しながら勉強したいと思っています。

最適化問題を解いてみる

今回は下のシンプルな目的関数の最小値を最適解として、PSOでパラメータを探索してみました。

f:id:hirotsuru314:20190409212421p:plain:w120

 N=2としたときの関数のグラフを描くと以下のようになります。

f:id:hirotsuru314:20190409215708p:plain

左側が3Dに描いたもので、右側が2次元に写像したもの(左のグラフを上から見たもの)です。
グラフや数式から直感的にわかるように、この関数は (x_1,x_2)=(0,0)で最小値 0を取るため、今回の探索の目的は、 (x_1,x_2)=(0,0)を見つけることになります。

サンプルコードは以下の通りです。

package main

import (
    "github.com/tsurubee/eago"
    "log"
)

func objectiveFunc(x []float64) float64 {
    return x[0] * x[0] + x[1] * x[1]
}

func main() {
    pso := eago.NewDefaultPSO()
    pso.NParticle =  5
    pso.NStep = 20
    pso.Min = -20
    pso.Max = 10

    if err := pso.Minimize(objectiveFunc, 2); err != nil {
        log.Fatal(err)
    }
}

上のコードをmain.goとして以下のように実行します。

$ go run main.go
Step   0: Fitness=8.124 Position=[-2.168 1.851]
Step   1: Fitness=8.124 Position=[-2.168 1.851]
Step   2: Fitness=7.394 Position=[-1.618 2.186]
Step   3: Fitness=6.958 Position=[1.068 2.412]
Step   4: Fitness=6.958 Position=[1.068 2.412]
Step   5: Fitness=4.809 Position=[-2.106 0.611]
Step   6: Fitness=2.215 Position=[-1.052 1.053]
Step   7: Fitness=1.535 Position=[0.580 1.095]
Step   8: Fitness=0.163 Position=[0.023 0.403]
Step   9: Fitness=0.163 Position=[0.023 0.403]
Step  10: Fitness=0.055 Position=[-0.175 0.156]
Step  11: Fitness=0.055 Position=[-0.175 0.156]
Step  12: Fitness=0.055 Position=[-0.175 0.156]
Step  13: Fitness=0.055 Position=[-0.175 0.156]
Step  14: Fitness=0.055 Position=[-0.175 0.156]
Step  15: Fitness=0.055 Position=[-0.175 0.156]
Step  16: Fitness=0.055 Position=[-0.180 0.151]
Step  17: Fitness=0.044 Position=[0.019 0.210]
Step  18: Fitness=0.013 Position=[-0.111 0.023]
Step  19: Fitness=0.008 Position=[-0.083 -0.031]

出力結果を見てみると、Stepが進むにつれて、Fitnessが 0、パラメータの解が (x_1,x_2)=(0,0)に収束しようとしていることがわかります。
粒子が速度を調整しながら最適解に収束して行く様子は下のサイトのアニメーションが非常にわかりやすかったです。

Particle swarm optimization - Wikiwand

さいごに

群知能とか生物模倣的なアプローチはもう少し勉強を続けて深掘りしていきたい。
あと、もっと実践的なデータに対して適用してみて、自分の今の領域における応用についても考えていきたい。

参考

shop.ohmsha.co.jp

遺伝的アルゴリズムを用いてコンテナの配置を最適化する論文:「Genetic Algorithm for Multi-Objective Optimization of Container Allocation in Cloud Architecture」

ふとタイトルが目に入ってきて気になって読んでみた。私自身そんなにコンテナ技術を触ってないけど、前から遺伝的アルゴリズムが気になっていた。 ちょいちょい知識が足りずに理解しきれないところ出てきたけど、とりあえず最後まで読んだのでざっくりまとめを書いてみる。

読んだ論文

pdfはダウンロードできなかったけど、ここからOnlineで読めた。

所感

先に論文を読んだ所感を書いておく。

  • コンテナのリソース割り当ての最適化問題遺伝的アルゴリズム(GA)を使って解くというコンセプトが面白かった。Related Workとして多くの論文が挙げられていた(Table 1およびTable2)ので、この辺りは追っていくと面白そう。

  • GAはパラメータの全検索が不可能と思われるくらい膨大な組み合わせがあって、かつ有効な探索方法が定かでないパターンに有効っぽい。さらに、論文内で用いているNon-dominated Sorting Genetic Algorithm Ⅱ(NSGA-Ⅱ)などの手法は、多目的最適化に有効な手法だから、けっこう応用の幅が広そうだなぁと思った。

  • いかに最適化したいものを目的関数として定式化するかっていうところが大変で重要なところなんだろうけど、それができるならば、目的は複数あっても、目的同士が相関関係にあっても、近似解をヒューリスティックに探索してくれるから、GAという手法は学んでいくといろんな問題に適用できそう。

  • 遺伝的アルゴリズムを始めとした進化的計算という分野に興味が湧いたのでやっていきたい。

概要

クラスター内の物理ノード群に対するコンテナの配置は、システム全体のパフォーマンスや信頼性、スケーラビリティに大きな影響を与えるため、システム設計において重要な要素である。Kubernetesのような既存のコンテナクラスター管理は、コンテナのオートスケールや物理マシンへのコンテナの配置に対して単純な実装しかされておらず、それらは物理的なリソースの使用率や閾値にのみ焦点を当てている。

このようなリソース最適化はNP完全問題であり、メタヒューリスティックなアプローチが有効である。中でも、GAのような進化的アプローチはリソース最適化問題に非常に有効な解法の一つである。本論文では、コンテナのリソース割り当てに対してNon-dominated Sorting Genetic Algorithm Ⅱ(NSGA-Ⅱ)というアルゴリズムで最適化している。NSGA-Ⅱという手法はMulti-Objective(多目的)な最適化に有効であるため、今回のように最適化したい目的関数が複数あって、それぞれの目的関数同士が独立でない場合でも適用可能である。

提案手法を用いて最適化したコンテナ配置は、Kubernetesによるコンテナ管理のアプローチよりも、論文中に定義した目的(後述)において、優位な結果が出ることが実験的に示されている。

遺伝的アルゴリズム(GA)

GAとは何か?という問いに対する答えとして「遺伝的アルゴリズム(Genetic Algorithm)を始めよう!」の7枚目のスライドがわかりやすかった。

生物進化における遺伝と適者生存による自然淘汰の仕組みをソフトウェア的に模すことで複雑な問題に対する最適解を探索する手法

私の印象としては、まじめに解こうとしても解けないくらい解の組み合わせの数が膨大なときに、実用的な時間スケールで近似的な解を探索できる点が大きな特徴だと思った。
また、実際に論文内の最適化で用いているNSGA-Ⅱについてもいくつか解析記事があった。

http://mikilab.doshisha.ac.jp/dia/research/mop_ga/moga/3/3-5-5.html

この辺りはこれから学んで数学的な背景からしっかり押さえていくので、また詳しくブログに書こうと思う。

システムのモデル化

最適化問題として解くためには、当然ながら何が最適化を評価できなければならない。さらに出てきた解(GA的には染色体)の候補同士に優劣をつけるために、定量的に評価できる形、つまり定式化しておく必要がある。

本論文ではコンテナのリソース割り当てに対して、以下の4つの目的を置いて定式化している。式中の変数についてはTable3を、定式化の詳細については「3.1 Optimization Objectives」を参照。

これらの定式化したものは、個体(解の候補)の適応度(fitness)を定量的に評価するために用いられる。

1. ボトルネックを回避するためにクラスタ全体でワークロードを均等に分散させる(Threshold Distance)

f:id:hirotsuru314:20190328222537p:plain

2. 新しいアプリケーションの追加を考慮した上で、クラスタのリソースをバランスよく利用する(Balanced Cluster)

f:id:hirotsuru314:20190328222715p:plain

3. クラスタ内のノードに適切にコンテナを分散させることによってマイクロサービスの信頼性を高める(System Failure)

f:id:hirotsuru314:20190328222655p:plain

4. 関連するマイクロサービスのコンテナを短いネットワーク距離の物理マシンに配置することによって、マイクロサービス間の通信のオーバーヘッドを小さくする(Network Distance)

f:id:hirotsuru314:20190328222846p:plain

実験と結果

実験については、GAで最適化したコンテナの配置とKubernetesアロケーションポリシーによるアプローチのそれぞれをサンプルアプリケーションを使ってベンチマークを取って比較している。

比較対象のKubernetesについては以下のように説明されている(私はKubernetesを触ったことないのでピンときていない)
KubernetesはコンテナやPodの配置を管理するために2つのポリシーを採用している。

  1. PodFitResources:物理マシンへのリソース要求の合計がマシンのキャパシティを超えないようにする

  2. LeastRequestPriority:コンテナは一番リソース消費が少ないマシンに配置される

Kubernetesの実装についてはよく知らないが、これを見るだけでも、論文中の目的関数4つの方が多くのことを考慮して最適化しようとしているのがわかる。
論文中では、GAの世代ごとの成長過程についても色々と議論されている。結論としてはGenerations 100、Population size 200というリーズナブルな値で最適化された解が得られたと述べられている。

Kubenetesとの比較の結果については、前述の4つの目的関数全てにおいて、提案手法のアプローチの方が優位な結果が出たと述べられている。例えば、一つの目的関数をピックアップすると以下のような感じ

f:id:hirotsuru314:20190328225640p:plain

変数として、物理マシンの台数を250・300・350・400台と降っている。こんな規模で実験やるんだなぁ。

さいごに

GAのこともKubernetesのこともよく知らないから、へーって感じだったけど、こういったインフラ周りで機械学習を適用している事例は面白いなと思った。GAについては数学的な理論からしっかり学んでいく。

Dynamic Time Warping(動的時間伸縮法)で時系列データをクラスタリングする

最近時系列データのクラスタリングに興味を持ち始めて、いくつか論文読んだり、アルゴリズムについて調べていたら、実装してみたくなったので勉強のために作ってみました。
実装の言語にはGolangを用いていて、クラスタリングアルゴリズムは、Dynamic Time Warping(以下、DTW)とk-medoids法を組み合わせたものです。

作ったもの

github.com

tsclusterはtime series clusteringの略です。 今回は作ったのは、特定のアルゴリズムのみですが、今後興味があるアルゴリズムがあれば、ここに実装していきます。

使い方

func main() {
    var dataset [][]float64
    dataset = append(dataset, []float64{1, 1, 1, 1, 1, 1})
    dataset = append(dataset, []float64{1, 1, 1, 1, 1, 1, 1})
    dataset = append(dataset, []float64{2, 2, 2, 2, 2})
    dataset = append(dataset, []float64{2, 2, 2, 2, 2, 2, 2})

    tc := tscluster.NewTscluster(tscluster.DTW) // 引数:距離関数
    labels, err := tc.Kmedoids(dataset, 2, 20)  // 引数:データ, クラスタ数, 最大反復回数
    if err != nil {
    log.Fatal(err)
    }
    fmt.Println(labels)
}

#=>
[0 0 2 2]

いくつか補足すると、与えるデータセットの型は[][]float64であり、1次元時系列データ([]float64)のスライスです。
また、出力結果labelsは、与えたデータセットのスライスに対してクラスタリングした結果のラベルを返しています。ラベルの値として、データが所属するクラスタのmedoids(後述するがクラスタの代表点のようなもの)のindex番号を返しています。
すなわち、[0 0 2 2] とは与えたデータセットの0番目と1番目は同じクラスタに属しており、datasetの0番目が代表点であるということです。

検証

それでは、もう少し実践的なデータで試してみます。
用いたデータとプログラムは、tsclusterのリポジトリexamplesディレクトにおいています。
1次元の時系列データを30個用意してグラフを描いたものが、以下になります。このデータは意図的に、上昇トレンド・停滞トレンド・下降トレンドそれぞれを持つものを10個ずつプログラムで用意しました。

f:id:hirotsuru314:20190217160343p:plain

これらのデータセットをDTW + k-medoidsでクラスタリングして、クラスタごとに色分けして描いたグラフが以下になります。

f:id:hirotsuru314:20190217160405p:plain

グラフを見るとそれぞれのトレンドごとにきれいにクラスタリングできているように見えます。
今回のような単純な線形のトレンドごとにクラスターに分けるにはもっと単純な手法もありそうです。本来DTWという距離尺度を用いていると、形状は似ているけど時間軸に対してずれているといったデータの類似性が高くなるような特徴を持つため、そういったデータをクラスタリングしてみるのも良いかもしれません。

実装したアルゴリズム

時系列クラスタリングの概要を把握するには、下のサーベイ論文が非常に情報がまとまっていました。

Time-series clustering – A decade review

Dynamic Time Warping(DTW)

DTWの特徴としては、

  • 時系列のデータ長がそろっていなくても計算可能
  • 時間・位相のずれていても形状が似ていれば類似度が高くなる

などが挙げられます。
DTWについてはこちらのブログに直感的に大変わかりやすいアニメーションがあったので、見てみるとイメージが湧きやすいと思います。

データ長がmとnの時系列データのDTW距離を求めるアルゴリズムは以下のステップで構成されます。

  1. 各データ点に対して総当たりで距離を計算し、要素がデータ点同士の距離となる行列D(m × n)を作成
  2. D[0][0] からD[m][n]までに通るパスのうち、行列の要素の合計が最小となるパスを探索する。その際のパスの要素の合計がDTW距離となる

k-medoids

k-medoidsはクラスタリングの手法で大変メジャーなk-meansと非常に似た手法です。
両者の大きな違いはクラスタの代表点をcentroid(重心)ではなくmedoidを選択するという点です。
medoidとは、クラスタ内のデータ点で、同一クラスタ内の他の全ての点までの距離の総和が最小になる点のことを言います。
Wikipediaには、k-meansと比較したk-medoidsの特徴が以下のように表現されていました。

It is more robust to noise and outliers as compared to k-means because it minimizes a sum of pairwise dissimilarities instead of a sum of squared Euclidean distances. (https://en.wikipedia.org/wiki/K-medoids より引用)

k-meansよりノイズや外れ値に強いとのことでした。また、k-medoidsはクラスタの代表点を実際のデータの中から選ぶため、一度データ同士の距離を全て計算してしまえば(距離行列を作成すれば)、k-meansのような重心の計算や重心からのデータ点までの距離の計算などといった処理がないため、計算量も少なくなると思います。

アルゴリズムとしては以下の通りです。

  1. データセットの中からk個のデータをランダムに選択する(modoid)
  2. 全てのデータを一番近いmedoidのクラスタに割り当てる
  3. それぞれのクラスタ内において、クラスタ内の他のすべての点との距離の合計が最小になる点を新たにmedoidとする
  4. 2と3を繰り返す。ただし、3においてmedoidに変化がない、もしくは指定した最大の反復回数に到達すれば終了する

さいごに

今回のようなデータのクラスタリング教師なし学習で、データがあればすぐに適用できるというメリットはあるものの、当然ながら出てきたクラスター群に対する解釈は人間が与えないといけないですし、そもそも何個のクラスターに分けるかはある程度事前知識が必要だと思います(クラスタ数を最適化する方法はあるものの)
また何か気になるものがあれば実装してみたいと思います!

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

現場ですぐ使える時系列データ分析 ~データサイエンティストのための基礎知識~

k-Shapeによる時系列クラスタリングの論文:「k-Shape: Efficient and Accurate Clustering of Time Series」を読んだ

最近、時系列データのクラスタリングに興味を持っているので、k-Shapeというクラスタリング手法に関する論文を読みました。
なぜ興味を持っているかというと、サーバの各種メトリクス(CPU使用率・メモリ使用率など)を使って、似たような特徴を持っているサーバ群をクラスタリングできないかと考えいるためです。例えば、負荷の高いサーバ群と負荷の低いサーバ群などにグルーピングできると、面白いのではないかと考えています。

元論文はこちらです。

k-Shape: Efficient and Accurate Clustering of Time Series

この論文では、提案するk-Shapeのアルゴリズムの説明だけでなく、時系列データをクラスタリングする上での理論的な背景から説明されていて、時系列クラスタリング手法の全体像を把握するのに大変良い論文でした。

概要

本論文では、k-Shapeという時系列クラスタリング手法を提案しています。

k-Shapeの特徴は、

  • 時系列データの形状に着目したshape-basedクラスタリングである
  • データ間の距離尺度として、規格化した相互相関を用いている
  • 高効率かつ高精度で、幅広いデータに適用できる

などが挙げられます。
以下、本論文の内容をまとめていきますが、理論の詳細に関しては、論文をご参照ください。

Time-Series Invariances

時系列データをクラスタリングする際には、何に着目してクラスタリングするか(何をもって似ているとするのか)というのが重要になります。
例えば、形状は似ているけど、時間軸に対してずれている2つの時系列を同一クラスと捉えるかどうかは、そのときの用途によって異なると思います。

k-Shapeの場合、ScalingShiftingに対する不変性 (Invariances)に着目しています。
簡単に言うと、Scalingは軸に対してスケールした際にデータ同士の性質が似ているかどうか、Shiftingは時間軸(位相)に対してずらした場合にデータ同士の性質が似ているかどうかに着目しており、この特徴からk-Shapeはデータの形状に着目したshape-basedなクラスタリングであるとされています。

Time-Series Distance Measures

クラスタリングを行う際に、データ間の類似性は距離として表現されます。すなわち、距離が近いほど類似性が高く、距離が遠いほど類似性は低いということです。

このことから、前述の「何を持って似ているとするか」は、データ間の距離をどのように定義するかに強く依存するため、距離尺度の定義が極めて重要となってきます。論文中でも距離をどうやって決めるかが、クラスタリングアルゴリズム自体より重要だと述べらていました。

時系列データの「形状」に着目したクラスタリングをするには、振幅と位相のずれをうまく扱う距離尺度が必要となります。
k-Shapeでは、一般的によく用いられるユークリッド距離や動的時間伸縮法ではなく、独自のShape-based distance(SBD)という距離尺度を用いています。

Shape-based distance(SBD)

SBDは、規格化した相互相関を用いています。
相互相関とは、信号処理や画像処理に広く使われている信号間の類似性の尺度です。
ざっくりいうとデータをずらしながら内積を計算していき、内積が一番大きくなった位置を探す感じです。

2つの時系列データxyにおけるSBDは以下の式で表されます。

f:id:hirotsuru314:20190206220202p:plain

ただし、

f:id:hirotsuru314:20190206220243p:plain

Shape-based Time-Series Clustering

距離尺度が決まってしまえば、クラスタの重心が計算できます。(重心の計算方法は論文中の「3.2 Time-Series Shape Extraction」を参照)
あとのクラスタリングアルゴリズムはk-meansとほぼ同じような感じで、以下のように反復的な処理を行います。

  1. 時系列データを各クラスタの重心と比較し、最も近い重心をもつクラスタにデータを割り当てる
  2. クラスタに属するメンバーから重心を計算して更新する

これを繰り返したのち、クラスタのメンバーに変更がないか、反復回数の最大値に達するまで上の手順を繰り返します。

感想

手持ちのデータをクラスタリングしたいときには、どの不変性に着目して、何を似ているとしたいのかを考えて、それに適したアルゴリズムを選定するかが重要だと思いました。
今回のアルゴリズムは、時間軸のずれや、伸縮に対して頑強なアルゴリズムである印象を受けました。サーバメトリクスの場合、scalingは重要な違いを表す要素だと思うので、サーバメトリクスのクラスタリングに対してはもっと適したアルゴリズムがあるかもなーと思いました。もっといろんな手法をサーベイしていきます。

k-Shapeに関しては、論文中に擬似コードが載っていたのでスクラッチで書けそうだし、tslearnというPythonライブラリもあるようなので試したい。

Krylov部分空間を導入して特異スペクトル変換による異常検知の処理を高速化した

1年くらい前に特異スペクトル変換法による異常検知ライブラリを作ったんですが、作ったっきり放置していたので、開発当初からやりたかった計算の高速化処理を書きました。
ずっと放置してた割にはちょいちょいGitHubのスターを押してもらえてて、データサイエンスの流行を感じた。自分ももう一回ちゃんと学び直していこうという気になったので、まずは昔書いたやつの拡張からやっていく。

【目次】

特異スペクトル変換とは?

特異スペクトル変換法の特徴については以前のブログに書いているので、ぜひそちらも読んでください。

特異スペクトル変換法の全体像は以下のようになっています。

f:id:hirotsuru314:20171011214314p:plain
出典:上の図は井手剛氏の著書「入門 機械学習による異常検知―Rによる実践ガイド」のP200 図7.4を元に作成しました。

図のように過去と今のパターンを行列として取り出し、それぞれの特徴抽出をした行列(上図の主部分空間の部分)を生成します。生成した二つの行列の食い違いの大きさを変化度として異常検知を行います。

特異スペクトル変換は、特定の確率分布を仮定していないため、多様な変化に頑強に対応できる点やパラメータチューニングが容易などといった特徴が挙げられます。

一方、特異スペクトル変換の一つの課題として計算コストが高いことが挙げられます。先ほど述べた「行列からの特徴抽出」の部分で、特異値分解と呼ばれる行列計算を用いており、これがかなり計算コストが高いです。したがって、検知を行いたい時間間隔や、パラメータの最適値によっては、リアルタイム計算が厳しくなってきます。

今回は、Krylov部分空間という概念を導入することにより、特異スペクトル変換による異常検知の処理を高速化することが目的です。

Krylov部分空間の導入

理論的な詳細は下記の元論文をご参照ください。この論文は今回実装した異常検知のアルゴリズムについてはもちろんのこと、部分空間同士の距離について丁寧に書かれていて数学的にも大変勉強になりました。

Change-Point Detection using Krylov Subspace Learning

Banpeiの既存実装だと、データから取り出した密行列に特異値分解かけるので、その行列のサイズが大きくなると処理はどんどん重くなりますが、論文のアルゴリズムを用いると、次元数k(通常2~5程度)のKrylov部分空間を導入して、最終的にはk×kの固有値問題に帰着します。

上の論文や他の参考情報(最後に載せている)を元に高速化処理を書きました。

github.com

検証結果

それでは、特異値分解(SVD)を使った既存実装と、今回実装したLanczos法を使った新規実装を比較していきます。(Lanczos法はKrylov部分空間法の計算手法の一つ)

用いるデータは、Banpeiのリポジトリに置いている周波数異常データです。
まず以下のようにしてデータを読み込みます。(ライブラリのインストール方法についてはREADME参照)

import banpei
import pandas as pd

model   = banpei.SST(w=100)
raw_data = pd.read_csv('./tests/test_data/periodic_wave.csv')
data = raw_data['y']

Jupyter Notebookを使うと以下のように%timeというマジックコマンドでプログラムの実行時間を簡単に計測できます。

# 既存実装(SVD)
%time results_svd = model.detect(data, is_lanczos=False)
#出力結果=> 
CPU times: user 4.22 s, sys: 80.4 ms, total: 4.3 s
Wall time: 2.6 s

# 新規実装(Lanczos)
%time results_lanczos = model.detect(data, is_lanczos=True)
#出力結果=> 
CPU times: user 944 ms, sys: 51.2 ms, total: 995 ms
Wall time: 872 ms

用いたサンプルデータ、既存実装により計算した変化度スコア、新規実装により計算した変化度スコアを描くと以下のようになりました。

f:id:hirotsuru314:20190203093607p:plain

実行時間に関しては、4.3sから995msとなり、約4.3倍早くなりました。
この結果は完全に取り出した行列のサイズに依存しています。今回の例で用いた行列のサイズは50×100ですが、このサイズが大きくなればなるほど、新規実装の方が顕著に処理が早くなります。

次に変化度スコアのグラフについてですが、一番下の新規実装(Lanczos)は既存実装(SVD)に比べてグラフがノイジーで粗くなっているように見えます。この理由は詳細には調べていませんが、異常検知の観点でいくと、しっかり変化度スコアが立ち上がってくれているなら、多少のグラフの粗さは問題ないかと思っています。

さいごに

難解な数学的な処理をプログラムに落とすと謎の達成感があるんだけど、書いたものは使ってなんぼなので、その辺りは今後の課題としてやっていきます。
データサイエンス周りは1年間くらい追ってなかったので、基礎的なインプットも再開してやっていく!

参考

ConsulのACL Bootstrapをリセットしてもう一度やり直す方法

こんにちは、つるべーです!
先日、Consul ACLの記事を書きましたが、今回もACLのちょっとした小ネタについて書きます。
内容としては、ACLをBootstrapしたあとにMater TokenのSecret IDをなくしてしまい、Consul関連の操作が何もできなって、発狂しそうになったときの対応方法です。

前回の記事

blog.tsurubee.tech

環境

DockerでConsul ACLの検証環境を用意しています。使い方はREADMEをご覧ください。

github.com

Consulのバージョンは1.4.0です。

/ # consul version
Consul v1.4.0
Protocol 2 spoken by default, understands 2 to 3 (agent will automatically use protocol >2 when speaking to compatible agents)

ACLのBootstrap

ACLを使い始めるには、最初にBootstrap、いわゆる初期化をしてやる必要があります。
公式ドキュメントではこのあたりに説明があります。

Consulの設定ファイルでACLを有効にしたのち、serverモードのサーバでBootstrapをすると、下のようにBootstrap Token(Master Token)が払い出されます。

/ # consul acl bootstrap
AccessorID:   3cff8208-af8e-feb9-9ebb-fc75a5ddb76e
SecretID:     91118e17-445b-942d-f090-f3b34defa4bf
Description:  Bootstrap Token (Global Management)
Local:        false
Create Time:  2018-12-21 14:26:58.574906326 +0000 UTC
Policies:
   00000000-0000-0000-0000-000000000001 - global-management

ACLのBootstrapは一度しかできない

上記のACLのBootstrapは一度しかできず、二回目を叩くと下のようにエラーが出ます。

/ # consul acl bootstrap -token=91118e17-445b-942d-f090-f3b34defa4bf
Failed ACL bootstrapping: Unexpected response code: 403 (Permission denied: ACL bootstrap no longer allowed (reset index: 18))

通常は、Bootstrapをし直す必要などないのですが、例えば、Bootstrapしたあとに、Bootstrap時に出力されるMater TokenのSecret IDをなくしてしまったらどうなるでしょう?
そうです、Bootstrapした直後の時点ではMaster TokenなしではConsulコマンドが何も実行できません。(Bootstrap時にAnonymous Tokenには何もPolicyがアタッチしていないため)
したがって、Bootstrap後にMater TokenのSecret IDがわからなくなるとConsul関連の操作が何もできなくなり、にっちもさっちもいかなくなります。

Bootstrapをリセットしてやり直す方法

Consulの公式ドキュメントを探しても、Bootstrapをもう一回やり直す方法が見当たりませんでした。(2018年12月22日現在)
しかし、同じくhashicorpが開発しているNomadのドキュメントにResetting ACL Bootstrapの説明がありました。このNomadのドキュメントと同じ手順でConsul ACLのBootstrapをリセットすることができました。

実際に手順を示します。まず二回目のBootstrapを叩いた時に出力されるエラーメッセージのreset index: 18の部分に着目してください。

ACL bootstrap no longer allowed (reset index: 18)

このindexの値は固定値でなく、個々で異なる値です。(クラスタ内では共有の値)
この値を埋め込んだファイルを、acl-bootstrap-resetという名前で、Leaderノードのデータディレクトリ配下に作成します。データディレクトリとは設定ファイルのdata_dirの部分です。(ここでは仮に/tmp/dataとする)

$ echo 18 >> /tmp/data/acl-bootstrap-reset

このようにreset indexを埋め込んだacl-bootstrap-resetファイルを作成すると、Bootstrapコマンドを再度実行できるようになります。

/ # consul acl bootstrap
AccessorID:   180d0d38-424f-69be-1897-d98b41abed7e
SecretID:     32ff58d9-df47-05d0-1279-f63ae28af609
Description:  Bootstrap Token (Global Management)
Local:        false
Create Time:  2018-12-22 07:09:03.8573353 +0000 UTC
Policies:
   00000000-0000-0000-0000-000000000001 - global-management

ここで全ての権限が与えられたTokenができるため、このSecretIDでConsulの操作ができるようになります。
ちなみに作成したacl-bootstrap-resetファイルは勝手に削除されています。

/ # cat /tmp/data/acl-bootstrap-reset
cat: can't open '/tmp/data/acl-bootstrap-reset': No such file or directory

また、実は再Bootstrapをしたからといって、完全に初期化されたのでなく、1回目のBootstrap時に作成されたMater TokenとAnonymous Tokenは残ったまま、新しいMaster Tokenが作られた形になります。

/ # consul acl token list -token=32ff58d9-df47-05d0-1279-f63ae28af609
AccessorID:   180d0d38-424f-69be-1897-d98b41abed7e
Description:  Bootstrap Token (Global Management)
Local:        false
Create Time:  2018-12-22 07:09:03.8573353 +0000 UTC
Legacy:       false
Policies:
   00000000-0000-0000-0000-000000000001 - global-management

AccessorID:   3cff8208-af8e-feb9-9ebb-fc75a5ddb76e
Description:  Bootstrap Token (Global Management)
Local:        false
Create Time:  2018-12-21 14:26:58.574906326 +0000 UTC
Legacy:       false
Policies:
   00000000-0000-0000-0000-000000000001 - global-management

AccessorID:   00000000-0000-0000-0000-000000000002
Description:  Anonymous Token
Local:        false
Create Time:  2018-12-21 14:18:38.848568562 +0000 UTC
Legacy:       false
Policies:

そのため、古いMaster Tokenは消しておいてもいいですし、新しいMaster Tokenが発行された今であれば、古いMaster TokenのSecretIDを確認することもできます。

/ # consul acl token read -id=3cff8208-af8e-feb9-9ebb-fc75a5ddb76e -token=32ff58d9-df47-05d0-1279-f63ae28af609
AccessorID:   3cff8208-af8e-feb9-9ebb-fc75a5ddb76e
SecretID:     91118e17-445b-942d-f090-f3b34defa4bf
Description:  Bootstrap Token (Global Management)
Local:        false
Create Time:  2018-12-21 14:26:58.574906326 +0000 UTC
Policies:
   00000000-0000-0000-0000-000000000001 - global-management

さいごに

最初Bootstrapをもう一回やりたいとなった時、bootstrapコマンドに--forceみたいなオプションないかなぁーと思ったのですが、ありませんでした。
今回見たように再Bootstrapできたとしても、古いTokenが消えるわけではないためすぐには事故に繋がらなさそうですが、Consulでは一意なreset indexを指定のファイル名・パスに埋め込むというオペミスが起こり得ない安全な設計になっており、自分で何かしらツールを作る時にもこういった設計は参考になりそうだなぁと思いました。

consul image