データサイエンティスト(仮)

元素粒子論博士。今はデータサイエンティスト(仮)。

Memo:BoostのVersionを確認

BoostのVersionを確認したいときに実行するコード

会社のサーバでとあるソフトウェアをPythonとbindingさせようとしたときに、BoostのVersionが知りたくなったので調べてみました。
以下のようなコードを作成して実行するとVersionがわかるようです。

#include <iostream>
#include <boost/version.hpp>

int main()
{
    int major = BOOST_VERSION / 100000;
    int minor = BOOST_VERSION / 100 % 1000;
    int patch = BOOST_VERSION % 100;
    std::cout << "boost version " << major << "." << minor << "." << patch
        << " or " << BOOST_LIB_VERSION << std::endl;
    return 0;
}


ただググっただけですが、ソース元はこちら。
boost/version.hpp - 1.64.0
e-kwsm.hatenablog.com

トポロジカルデータアナリシス:単体の境界

導入

不定期でトポロジカルデータアナリシス(TDA)に関する紹介をします。今回は単体の境界を数学的にどう表現するかを紹介します。

振り返り

前回、図形の穴は「境界」に着目することで特徴づけができそうだということを紹介しました。

tekenuko.hatenablog.com

今回は、境界をざっくり紹介していきます。その際に、これまで表記を曖昧にしてきたものを幾つか定義しておきます。

単体および向き

単体

単体とは、点(0次元)や辺(1次元)、三角形(2次元)、四面体(3次元)といった図形を一般化したものでした。これらは、N次元空間(以下、\mathbb{R}^Nと表記します)にあるk次元三角形と拡張されます。このようなk次元三角形をk単体といいます。今、k単体の頂点をv_0, v_1, \cdots v_kとした場合、k単体は|v_0 \cdots v_k|と表します。本当はいくつか条件がありますが、ここではざっくりした言い方にとどめておきます。

単体の向き

先程の記法は頂点の順序、つまり向きの情報がありませんでした。向きを明示的に表す記法として、ここでは\langle \rangleを採用します。例えば、2単体|v_0v_1v_2|で頂点の順序がv_0 \rightarrow v_1 \rightarrow v_2と決まっている場合、向きを加味した単体は\langle v_0 v_1 v_2 \rangleと書きます。さらに、ルールとして、ある頂点を隣の頂点と入れ替えた場合はマイナスがつく、というものを設定しておきます。例えば、向きがある1単体 \langle v_0 v_1 \rangleだと
\[ \langle v_0 v_1 \rangle = - \langle v_1 v_0 \rangle \]
といった具合です。上記の向きづけは、厳密には置換という数学を使って定式化をするのですが、ここでは立ち入らないことにします。

境界作用素

境界は、ざっくり言うと境界作用素という、k単体に対するある決まった演算として定義されます。もう少し厳密には、k単体のある集まりに対して定義されるものなのですが、ここではこれ以上立ち入らないことにします*1。今、向きづけられたk単体を\langle v_0 \cdots v_k \rangleとします。k単体に対する境界作用素\partial _kとすると、その作用は
\[\partial _k \langle v_0 \cdots v_k \rangle = \sum _{i = 0}^k (-1)^i \langle v_0 \cdots \hat{v}_i \cdots v_k \rangle \]
というように定義されます。ここで、\langle v_0 \cdots \hat{v}_i \cdots v_k \rangle k-1単体で、\hat{v_i}は頂点v_iが除かれていることを意味しています。つまり、境界作用素k単体をk-1単体へ変える働きを持っています。しかも、ただ単に変えるだけではなく、「向き」も考慮しています。それが式中にある(-1)^iに対応しています。

境界作用素を作用させた例

1単体  \langle v_0 v_1  \rangle

\[ \partial _1 \langle v_0 v_1 \rangle = \langle v_1 \rangle - \langle v_0 \rangle\]

0単体でいうと向きづけは入ってくる方向がプラス、出て行く方向がマイナスとなっており、もともと、1単体としてはv_0 \rightarrow v_1という向きがあったことを鑑みると、境界作用素は向きを反映しつつ端っこ(0単体)を抽出していることがわかります。

2単体  \langle v_0 v_1  v_2 \rangle

\[ \partial _2 \langle v_0 v_1 v_2 \rangle = \langle v_1 v_2 \rangle - \langle v_0 v_2 \rangle + \langle v_0 v_1 \rangle = \langle v_0 v_1 \rangle + \langle v_1 v_2 \rangle + \langle v_2 v_0 \rangle\]

わかりやすくするために少し変形をしました。作用させるもとの2単体は、「反時計まわり」的な向きを持っていますが、境界作用素を作用させた結果は、その向きを反映したような端っこ(1単体の集まり)になっていることがわかります。以下に図的に表現したものを載せておきます。

f:id:tekenuko:20161213231451j:plain

境界の境界は無い

最後に、境界作用素の少しおもしろい性質を紹介します。それは、ある単体に境界作用素を2回作用させると必ず0になってしまう、というものです。2単体  \langle v_0 v_1  v_2 \rangleを例に見てみましょう。  \langle v_0 v_1  v_2 \rangle\partial _2を作用させたものにさらに\partial _1を作用させると
\[ \partial _1 \partial _2 \langle v_0 v_1 v_2 \rangle = \partial _1 \left( \langle v_1 v_2 \rangle - \langle v_0 v_2 \rangle + \langle v_0 v_1 \rangle\right) \\ = \langle v_2 \rangle - \langle v_1 \rangle - \langle v_2 \rangle + \langle v_0 \rangle + \langle v_1 \rangle - \langle v_0 \rangle = 0\]
のように、うまく打ち消し合いが起きて0になります。これは一般のk単体でも成り立つ性質です。

境界から穴の議論をさらにするために

今回、ちょっとごまかしつつも境界を表現する境界作用素という道具を紹介しました。ただし、ここから「図形に穴が空いている」ことを議論するためには

  • 今考えている図形の大元に対応する単体的複体を明示的に考える
  • 単体を「基底」とみなすようなベクトル空間的なものを考え、その空間に対してきちんと境界作用素を考える

といったことを詰めていかなければなりません。ネタを小出しにしていっている感が満載ですが、次回以降はこれらに関係する話を紹介して行こうと考えています。

*1:ちゃんと単体的複体を定式化したり、群などの特定の演算が入っている数学的対象を考えなければならなくなるからです。

トポロジカルデータアナリシス:ホモロジー群の紹介のための準備

導入

不定期でトポロジカルデータアナリシス(TDA)に関する紹介をします。今回は、図形の「穴」を数学的に表現するための準備をします。具体的に表現していくのは次回以降になります。

振り返り

これまで、データから図形を見立てる方法を紹介してきました。

tekenuko.hatenablog.com
tekenuko.hatenablog.com

これまでの話では、見立てた図形から情報を抽出する方法を紹介していませんでした。ここでは、図形の「穴」に着目するという考えを採用します。この「穴」に着目するという考えを数学的に表現するもの(の一つ)がホモロジー群です。

図形の穴を表現する方法は他にもあり、例えば、ホモトピー群というものがあります。これも図形の形や穴の数を反映するような量になっています。理論物理学ではホモトピー群のほうがよく応用に使われ、例えば液晶や渦、宇宙ひもやドメインウォール*1の分類や、非線形微分方程式の特解(ソリトン解)の分類などの応用があります。ただし、ホモトピー群ホモロジー群に比べて具体的に計算することは非常に難しく、計算機などに載せた応用には向いていません。そのため、上記のような具体的な図形の分類の問題ではホモロジー群がよく用いられます。

以下では、「穴」を数学的に表現するための準備的な考え方を紹介していきます。

穴は「境界」で特徴づけられる

今回、1単体(= 辺)の集まりで穴を作る例と作らない例を見ることで、穴を数学的に表現する手がかりを得られないかを考えていきます。ここで、単体とは、点(0次元)や辺(1次元)、三角形(2次元)、四面体(3次元)といった図形を一般化したものです。これらをそれぞれ0次元三角形、1次元三角形、...と呼ぶことにします。これは、N次元空間内のk次元三角形という拡張ができます。このk次元三角形をk単体といいます。単体、およびそれを集めた単体的複体に関しては、ざっくり以下の記事でも紹介しています。

tekenuko.hatenablog.com

以下、1単体の集まりで穴を作る例(左)と作らない例(右)です。
f:id:tekenuko:20161212220610j:plain

まず、左図は直感的にも穴が空いていそうですが、それはおそらく「端っこがない」ように見えるからだと思います。この「端っこ」は、数学的には「境界」と呼ばれます。例えば、頂点1と2からなる1単体の境界(端っこ)は頂点1と2です。しかしながら、頂点1と2には別の1単体(頂点1と3、および2と3からなる1単体)がくっついています。これは頂点1と3、および2と3からなる1単体に関しても同様です。つまり、これら3つの1単体の集まりには境界がありません。

一方で、右図は直感的には穴が空いていなそうです。ただし、1単体の集まりだけに着目すると、左図と構造は同じです。では、なぜ穴が空いていなそうと思うのでしょうか。

それは、1単体の中身を埋めるように2単体(三角形)が存在しているからです。これは、2単体の「端っこ」が1単体の集まりになっているとも言いかえられます。つまり、1単体の集まりに着目すると端っこを持っていませんが、この1単体の集まりは2単体の端っこになっています。

以上のことから、1単体に関しては穴を持つ条件は以下になりそうです。

境界(端っこ)を持たない1単体の集まりで、2単体の集まりの境界になっていないもの」

つまり、ある次元の単体が与えられたとき、「境界」というものに着目すると穴による特徴づけが行なえそうです。実際、上の1単体の例を数学的に表現したものが1次ホモロジー群となっています。また、ここでの類推から、1単体をk単体、2単体をk+1単体とすればkホモロジー群という拡張も行えるのではないかと期待されます。

さらに数学的に定式化するために

今回は、図形の穴の分類において、「境界」というものに着目していけば良さそうであることを説明しました。次回以降、境界の数学的な表現方法を紹介していきます。ただし、そのためには単体などに関する表記や向きなどをきちんと定義した上で進める必要がある(というよりこれらを数式で表現しないと逆にわかりづらくなる)ので、以前の記事よりも若干かっちり定式化しつつ進めていきたいと思っています。

*1:これらはざっくりいうと、初期宇宙で対称性の破れによってできた一次元、二次元的な欠陥のようなものです。

時系列データでお絵かきする

はじめに

最近、趣味の時間があまり取れなかったので更新が滞ってました(汗)。

周期性をもつ(だろう)時系列データをお遊びで可視化してみました。やってみるとなかなか面白かったので、備忘録として残しておくことにしました。

どうやって可視化するのか

ある時点tでの観測値をy_tとするような時系列データ\{y_t, t = 1, \cdots, T \}を考えます。y_tと1期前の観測値y_{t-1}を成分にもつ2次元ベクトル
{
\displaystyle
\begin{eqnarray}
\xi _t= \begin{pmatrix}
              y_t \\
              y_{t-1} \\
           \end{pmatrix} 
\tag{1}
\end{eqnarray}
}
を考え、これをtを変えてプロットします。今回、1期前だけを考えましたが、並べるのは何期前のものでもよいですし、さらに過去の観測値を並べて高次元のベクトルを考えても良いです。とにかく現時点と過去の観測データを幾つかチョイスし、それらを成分として構成したベクトルをプロットする、ということが具体的にやることです。

可視化

サンプルデータ1 : 単一の周期性をもつノイズのない系列

単純な場合として、観測値がy_t = \cos(t)となっている場合を考えます。

# ライブラリ
if(!require(data.table)){
  install.packages("data.table", dependencies = TRUE)
}
require(data.table)
if(!require(dplyr)){
  install.packages("dplyr", dependencies = TRUE)
}
require(dplyr)

# サンプルデータを100期分作成
t <- seq(1, 100, 1)
x1 <- cos(t)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

観測値をobsとし、obs_lagはtを1期さかのぼったデータが格納されています。

観測値をプロットすると、たしかに周期性をもつような振る舞いのデータになっています。

# ggplot
if(!require(ggplot2)){
  install.packages("ggplot2", dependencies = TRUE)
}
require(ggplot2)

# 観測値をプロット
ggplot(data = df, 
       aes(x = time, y = obs)
) + 
  geom_line()

f:id:tekenuko:20161110225451p:plain

次に、obsとobs_lagの散布図を描いてみます。

# プロット
ggplot(data = df, 
       aes(x = obs, y = obs_lag)
       ) + 
  geom_point(color = "hotpink")

結果は以下になります。
f:id:tekenuko:20161110231730p:plain

傾いた楕円のような形に見えるものが生成されました。このような形になることは、\sin^2(t) + \cos^2(t) = 1y_t y_{t-1}で表すことで式でも理解できますが、ここでは省略します。言いたいことは、周期性のある時系列データに関して、時点tと過去の時点の観測値を並べたものを見ると、なにか構造が見えることがあるということです。

傾いた楕円のようなものが生成されるのは、もとの時系列が単一の周期を持っていたからです。この振る舞いは、比較する観測値を変えても変わりません。比較として、y_tと2期前の観測値y_{t-2}を成分とするベクトルをプロットした場合も考えます。中心や傾き方は変わりますが、傾いた楕円のようなものが生成されます。

# 2期前の観測値を並べる
df[, obs_lag2 := lag(obs, 2)]
# プロット
ggplot(data = df, 
            aes(x = obs, y = obs_lag2)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110231928p:plain

サンプルデータ2 : 単一の周期性をもつノイズがある系列

さきほどはきれいな形が見れましたが、ノイズが付与されたデータだとどうでしょうか。今度は、\cos(t)にホワイトノイズN(0, 0.1)が付与されたデータを考えてみます。こちらの場合にobsとobs_lagの散布図を描くと

t <- seq(1, 100, 1)
x1 <- cos(t) + rnorm(100, 0, 0.1)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

# プロット
ggplot(data = df, 
            aes(x = obs, y = obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110233410p:plain

となります。目を凝らすと楕円に見えなくもないです。このように、ノイズがある場合は、なんらかの構造にノイズがのったものが生成されます。

注意 : トレンドがある場合

トレンド(長期的な上昇/下降傾向)があるような時系列データの場合は、差分系列にしてからベクトル化します。以下のような上昇トレンドをもつようなデータを考えます。

t <- seq(1, 100, 1)
x1 <- cos(t) + t + rnorm(100, 0, 0.1) 
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

# プロット
g <- ggplot(data = df,
            aes(x = time, y = obs)
) +
  geom_line()
ggsave("fig_cau.png", g)

f:id:tekenuko:20161110234038p:plain

この場合、obsとobs_lagの関係は以下になります。

ggplot(data = df, 
            aes(x = obs, y = obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161110234152p:plain

楕円とは似ても似つかないものが生成されました。これは、y_tが周期関数以外の部分でtにあらわに依存しているためです。そのため、\sin^2(t) + \cos^2(t) = 1y_t y_{t-1}で表すと、tの寄与が残ってしまい、y_t y_{t-1}の関係式((y_t, y_{t-1})平面での軌道に対応)がtに依存する形になってしまうのです。よって、軌道が閉じなくなり、上のような直線っぽい振る舞いになります。

周期性そのものを可視化したい場合は、このトレンドは調整した後にベクトル化する必要があります。差分系列を取った後、ベクトル化したものを可視化すると以下のようになります。

# 差分系列を作成
df[, diff_obs := ifelse(t == 1, 0, obs - obs_lag)]
df[, diff_obs_lag := lag(diff_obs)]

# プロット
ggplot(data = df, 
            aes(x = diff_obs, y = diff_obs_lag)
) + 
  geom_point(color = "hotpink")

f:id:tekenuko:20161111000219p:plain

今度は傾いた楕円にノイズがのったようなものが生成されました。

サンプルデータ3: 多重周期性をもつノイズのない系列

もう少し複雑な例として、複数の周期性をもつ系列を考えます。観測値がy_t = 3 \sin(t) + 2 \cos(t / 3)となっている場合を仮定します。可視化したときの見栄えの都合上、1000期分のデータを用意します。

t <- seq(1, 1000, 1)
x1 <- 3 * sin(t) + 2 * cos(t / 3) + rnorm(100, 0, 0.1) 
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

100期までの観測値をプロットすると、以下のようになります。
f:id:tekenuko:20161111001519p:plain
2つの周期的な変動が重なったような振る舞いになっています。

こちらもobsとobs_lagを可視化すると、以下のようになります。
f:id:tekenuko:20161111002017p:plain
単一の周期性をもつ場合に比べると、かなり複雑な構造になりました。線が何回重なるか、何回小さな円っぽい軌道を描くか、などは周期が変わると大幅に変わります。今回は\sin(t)\cos(t/3)を比較すると、周期が3倍異なるものが2つある、ということで3回回ってもとのところに戻ってくる、という振る舞いになってそうです。

サンプルデータ4: 多重周期性と周期的なスパイクをもつノイズのない系列

さらに複雑な例として、さきほどの系列に周期的なスパイクを加えます。12期ごとに一回、大きな正の値の何かが印加される状況を考えます。

t <- seq(1, 1000, 1)
x1 <- 3 * sin(t) + 2 * cos(t / 3) + ifelse(t %% 12 == 0, 10, 0)
df <- data.table(time = t, obs = x1)
df[, obs_lag := lag(obs)]

100期までの観測値をプロットすると、以下のようになります。
f:id:tekenuko:20161111003001p:plain
たまに大きな値に観測値が振れるような感じになっています。

obsとobs_lagを可視化すると、以下のようになります。
f:id:tekenuko:20161111003023p:plain
何か両端に変なものが湧いてきました。。ちょっとこのままだとわかりづらいので、100000期のデータに増やしてプロットし直してみます。
f:id:tekenuko:20161111003402p:plain
どうやらさきほどの構造のコピーらしきものが両端にできるようです。大きな正の値が印加される、という状況は、ベクトル\xi_tの言葉でいうと、成分y_tのみが大きな値になる場合と成分y_{t-1}のみが大きな値になる場合の2パターンが新たに生まれることを意味しています。加えて、印加されるものも周期性を持っているため、新たに生まれた2パターンもさきほどの構造が反映された軌道を描く、といった理解になるんでしょうか。また、こちらの場合、データ点を増やさないと詳細が見えてこなかったということを鑑みると、実データの場合もそうとう蓄積しないとこのような構造が見えない可能性があるな、と思いました。

おわりに

観測値をうまくベクトルにして可視化すると、構造が見える場合があって面白いです。他にも色々なパターンを作れるので、やってみると楽しくて無限に時間を溶かせる感じになります。ここで載せたもの以外にも色々試しており、他にも面白いものもあったりする*1んですが、夜な夜なやりこみすぎると怒られたりするんで、このへんでやめときます。

*1:日中も空いた時間にやって一人でテンション上がってたりしてました(笑)

トポロジカルデータアナリシス:ノイズのあるデータへ応用する

導入

不定期でトポロジカルデータアナリシス(TDA)に関する紹介をします。今回は、ノイズをもったデータへのTDAの応用方法を紹介します。

振り返り

データが与えられたとき、それを高次元空間の点の集まりとみなし、それらから図形を見立てる方法を紹介しました。
tekenuko.hatenablog.com

このとき、点の周りにある半径の球を考えたのですが、その半径はよしなにとっていました。半径を変えると、見立てられる図形も変化してしまいます。すると、前回見立てた図形が本当にデータの意味のある情報を取り出していたかが疑問としてでてきます。
f:id:tekenuko:20161013072544j:plain

上の図は大きいリングと小さいリングが1個ずつあるのが真の構造だった場合の例です。実際のデータは、この構造にノイズが上乗せされた状態で観測されます。ここから図形を見立てようとした場合、球の半径を小さくとってしまうと小さいリングの構造は見えますが、大きいリングの構造が見えない、逆に球の半径を大きく取りすぎると小さいリングの構造が見えなくなる、といったことが起こります。つまり、ノイズをもったデータから安定した「構造」を取り出すのは一筋縄ではいかないということです。

この問題に対してTDAはどう答えるのでしょうか。

基本的なアイデア

TDAの考え方の一つに、パーシステントホモロジーというものがあります。これは、大ざっぱに言うと、点の周りの球の半径を連続的に変えていったときの「図形の移り変わり」をまるっと見る考え方です。特定の半径を持つ球から見立てる図形だけでなく、半径を変えていったときのパターンを総合的にみて、データの持つ情報を引き出していきます。パーシステントホモロジーを定式化するためには、様々な前提が必要ですが、ここでは「図形の移り変わりをまるっと見るための量」と思っておきましょう。

半径を変えていったときの移り変わりの様について、図を見ながら考えていきます。
f:id:tekenuko:20161013073656j:plain

先ほどと同じように、大きいリングと小さいリングが1個ずつあるのが真の構造だった場合の例を考えます。データ点の周りに球を考えますが、その半径を一斉に、かつ連続的に大きくしていきます。すると、最初は細かい構造だったものが

  • 小さいリングの構造が見え始める
  • 大きいリングと小さいリングが1個ずつある構造が見える
  • 大きいリングのみの構造が見える
  • 穴がなくなる

と移り変わります。ただし、それぞれの構造は、出現している区間の大きさに違いがあります。TDAでは、半径を大きくしていったときにすぐ構造が変わってしまうものをノイズ的なもの、構造が安定的(長い区間同じ構造を保っている)なものをデータの本質的な情報、といった見方をします。上の例だと、半径を大きくしていったときに、大きいリングと小さいリングが1個ずつある構造が比較的安定して存在していることが見えた場合、この構造が本質であったと判断するわけです。

パーシステント図

構造の移り変わりを構造の「生成」と「消滅」という見方で可視化したものを、パーシステント図といいます。ここでいう構造とは、例えば一次元ホモロジー群の生成元のことを指しますが、これはざっくり「何個かの穴をもったものを表現する量」のことです。

今、球の半径を一斉に、かつ連続的に大きくしていったとします。大きくしていく途中で、新しい構造がでてきた半径(Birth)と、その構造が消えた(つまり別の新しい構造が生まれた)半径(Death)を記録していきます。それらをプロットすると、下の図のようになります。
f:id:tekenuko:20161013075300j:plain

すぐ生成されてすぐ消滅する量は、それぞれの半径の大きさの差が非常に小さくなります。そのため、Death - Birth(寿命)が0に近い、つまり対角線上に近くなります。一方で、長い区間で安定的な構造は、Death - Birthが大きくなるため、対角線上から離れたところに位置します。先ほどの例だと、大きいリングと小さいリングが1個ずつある構造が対角線上から離れた場所にプロットされた場合、この構造が本質であったと判断します。

おわりに

非常にざっくりですが、TDAの考え方の一つであるパーシステントホモロジーの紹介をここではしました。このような考え方の応用先ですが

  • クラスタリング : 最適なクラスタ数を安定的な構造と結びつける
  • 構造解析・画像認識 : ノイズのあるデータから本質的な情報を抽出

といったものが素朴にはありそうです。ただし、実際の応用の事例はまだ多くありません。しかしながら、TDAはデータに対する新しい見方を与えるものであり、これが既存手法と比較してメリットがあることがわかれば、一気に世の中に広まるポテンシャルを秘めていると期待しています。

今回紹介したパーシステントホモロジーは、TDAの考え方の一つにすぎません。例えば、AYASDIというTDAベンチャー企業は、「Mapper」という別の手法を使用しているようです。幸い、パーシステントホモロジーとMapperの両者ともにソフトウェアが存在しています。TDAに関して更に知りたいという方はソフトウェアを使う、Mapperについて論文やAYASDI社の事例を調べてみる、といったことをやってみると良いのかな、と思います。