Noise-contrastive estimation の雑なメモ

word2vec のモデルの説明で、以下の尤度関数の最大化を行いたいが、分母にある正規化項の計算が重すぎてこのままでは計算できないという話があった。

{ \displaystyle
    \mathcal{L} = \prod_t \prod_{c \in c_t} p( c \mid w_t ) = \prod_t \prod_{c \in c_t} \frac{\exp( \vec{v}_c \cdot \vec{v}_{w_t} )}{\sum_{w' \in V} \exp( \vec{v}_{w'} \cdot \vec{v}_{w_t} )}
}

word2vec では、この尤度関数を直接最大化するのではなく、正例とノイズとを区別する2クラス分類問題の新たに考え、この問題の誤差関数を最小化することでパラメータを求めていた。 しかし普通に考えると、位置tのコンテキストの条件付き確率分布の尤度最大化と、観測データとノイズの2クラス分類問題は全く別の問題のように見える。

Chainer本では前者の問題を解くのと後者の問題を解くのとは「本質的に同じ」と書いてあった。 これらの本質的に同じとはどういうことなのか疑問に思ったので、negative sampling の元ネタらしい Noise-contrastive estimation の論文を読んでみた。

Noise-contrastive estimation: A new estimation principle for unnormalized statistical models

問題

ベクトル  {\bf x} \in \mathbb{R}^n が未知の確率分布  p_d({\bf x}) に従うとする。 この分布を、パラメータ  \alpha を持つ確率分布  p_m({\bf x}; \alpha) でモデル化し、尤度を最大化することでパラメータ  \alpha を決定することを考える。

確率分布は正規化されていなければならない。 つまり、定義域全体に渡る積分が1にならなければならない。

正規化されていない分布  p^0_m({\bf x}; \alpha) があるとき、正規化制約を満たすためには以下のように  p_m({\bf x}; \alpha) を定義すればよい。

{ \displaystyle
    p_m({\bf x}; \alpha) = \frac{p^0_m({\bf x}; \alpha)}{Z(\alpha)}, \ \ \ Z(\alpha) = \int p^0_m({\bf u}; \alpha) d{\bf u}
}

しかし、 Z(\alpha) が解析的に解けることはめったにないし、(word2vec の場合のように)数値的に計算するのもコストが高すぎて困難である場合が多い。 したがって、 Z(\alpha) を計算することなく、尤度を最大化する  \alpha を求めたい。

Noise-contrastive estimation

 Z(\alpha) を計算するのではなく、新たなモデルパラメータとしてモデルに組み込むことを考える。 つまり  p_m を以下のように再定義する。

{ \displaystyle
    p_m({\bf x}; \alpha, Z) = \frac{p^0_m({\bf x}; \alpha)}{Z}
}

 Z をパラメータにしてしまったので、この分布は  Z の値をうまく決めないと正規化されない。 このため、この分布の尤度を目的関数として最大化を行うことはできない。 なぜなら、 Z を小さくすれば尤度をいくらでも大きくできるからだ。

Noise-contrastive estimation では、観測されたデータである  {\bf x} と、人工的に作ったノイズ  {\bf y} \sim p_n をロジスティック回帰で識別する問題を代わりに考える。

 p_nノイズ分布 と呼ばれる分布であり、サンプリングすることができる分布であれば自由に選んでよいらしい。 ただし、 p_d が非ゼロである点では  p_n も非ゼロであることが望ましい。 さらに言うと、 p_d にできるだけ近い分布を選ぶのがよいらしい。

 X = ({\bf x}_1, \dots, {\bf x}_T)を観測データとし、 Y = ({\bf y}_1, \dots, {\bf y}_T) をノイズ分布  p_n からサンプルしたデータとする。 このとき、NCE において最大化すべき目的関数(ロジスティック回帰の尤度関数)は以下のように与えられる。

{ \displaystyle
    \mathcal{L}(\alpha, Z) = \sum_{t=1}^T \left\{
        \log \left[ h({\bf x}_t; \alpha, Z) \right] +
        \log \left[ 1 - h({\bf y}_t; \alpha, Z) \right]
    \right\}
}

ただし、 h({\bf u}; \alpha, Z) は、  \sigma(v)シグモイド関数として、以下のように書ける関数。

{ \displaystyle
    h({\bf u}; \alpha, Z) = \sigma\left( \log p_m({\bf u}; \alpha, Z) - \log p_n({\bf u}) \right)
}

つまり、 p_m p_n の対数尤度比をシグモイドに突っ込んだのが  h

いろいろ条件を仮定すると、この目的関数を最大化するパラメータを  (\alpha^{*}, Z^{*}) としたとき、 \alpha^{*} は元々の問題の尤度関数を最大化する  \alpha と一致し、 Z p^0_m を正しく正規化する値になる ことが示せるらしい。すごい。

まとめ

ということで、「本質的に同じ」というのは「どっちの問題を解いても同じモデルパラメータが得られる」という意味らしい。 ただし、negative sampling は NCE をベースとしているけど異なる手法なので、negative sampling で本当に同じパラメータが得られるのかはよくわからない。 後で調べたい。

Word2Vec メモ その2

昨日の記事の続き。

Word2Vec を Chainer で実装していく。 完全なコードは以下の URL にある。

workspace/learning-chainer/word2vec at master · nojima/workspace · GitHub

model

誤差関数をネットワークとして表現すると下図のようになる。

f:id:nojima718:20170821234058p:plain

普通のニューラルネットワークと違って活性化関数はない。

WW' が入力である単語の hot vector を低次元のベクトル表現に変換するための行列で、Chainer 本の図だと  W = W' になっている。今回の実験では、 W W'が同じ行列の場合と違う行列の場合を両方試してみた。

 W = W' の方のモデルを Chainer のコードに落とすとこんな感じになる。

class Word2Vec(Chain):
    def __init__(self, n_vocabulary: int, n_units: int) -> None:
        super().__init__()
        with self.init_scope():
            self._embed = L.EmbedID(n_vocabulary, n_units)

    def __call__(self, x1: Variable, x2: Variable, t: Variable) -> Variable:
        output = self.forward(x1, x2)
        return F.sigmoid_cross_entropy(output, t)

    def forward(self, x1: Variable, x2: Variable) -> Variable:
        h1 = self._embed(x1)
        h2 = self._embed(x2)
        return F.sum(h1 * h2, axis=1)

(変数名をこれまでの説明と揃えるべきなんだけど、時間がなくて全然揃えれてないので、すみませんが雰囲気で読み替えてください)

L.EmbedID は one-hot ベクトルに最適化された線形レイヤーで、L.Linear よりも forward も backward も効率的に実装されている。L.EmbedID.__call__ は単語ベクトルじゃなくて単語IDを取ることに注意。

 W W' を異なる行列にする場合は L.EmbedID を2つ用意して使い分ければいい。

train

モデルを学習するためのコードを書いていく。

まずハイパーパラメータを保持しておくための型を定義する。

HyperParameters = namedtuple('HyperParameters', [
    'vocabulary_size',      # 語彙数
    'vector_dimension',     # 分散表現ベクトルの次元
    'window_size',          # ウィンドウの半径
    'n_negative_samples',   # 1単語につきネガティブサンプルする単語の数
    'batch_size',           # 1回のミニバッチで処理する単語数
    'n_epoch',              # 総エポック数
])

次に negative sampling で使うためのサンプラを作る。

Chainer には WalkerAlias という便利なサンプラがある。 WalkerAlias のコンストラクタに頻度の列を渡すと、その分布に従うサンプラが作れる。

def _make_sampler(dataset: DataSet) -> walker_alias.WalkerAlias:
    _, counts = np.unique(dataset.data, return_counts=True)
    counts = np.power(counts, 0.75)
    return walker_alias.WalkerAlias(counts)

DataSet は、単語IDの列 data と語彙集合 vocabulary を保持する型)

次に、train 関数を定義する。普通にミニバッチで学習する。ここで出てくる _make_batch_set は訓練データを作る関数で、後で定義する。

def train(dataset: DataSet, params: HyperParameters) -> Word2Vec:

    model = Word2Vec(params.vocabulary_size, params.vector_dimension)

    optimizer = optimizers.Adam()
    optimizer.setup(model)

    sampler = _make_sampler(dataset)

    batch_size = params.batch_size
    n_epoch = params.n_epoch

    for epoch in range(n_epoch):
        indices = np.random.permutation(dataset.size)

        for i in range(0, dataset.size, batch_size):
            model.cleargrads()

            x1, x2, t = _make_batch_set(
                indices[i:i+batch_size], sampler, dataset, params)
            loss = model(x1, x2, t)
            loss.backward()
            optimizer.update()

    return model

後回しにしていた _make_batch_set の定義はこれ。 昨日の記事で説明した方法にしたがって正例と負例を作って返す。

def _make_batch_set(
        indices: np.ndarray,
        sampler: walker_alias.WalkerAlias,
        dataset: DataSet,
        params: HyperParameters) -> Tuple[Variable, Variable, Variable]:

    window_size = params.window_size
    n_negative_samples = params.n_negative_samples

    x1, x2, t = [], [], []

    for index in indices:
        id1 = dataset.data[index]

        for i in range(-window_size, window_size+1):
            p = index + i
            if i == 0 or p < 0 or p >= dataset.size:
                continue
            id2 = dataset.data[p]
            x1.append(id1)
            x2.append(id2)
            t.append(1)
            for nid in sampler.sample(n_negative_samples):
                x1.append(id1)
                x2.append(nid)
                t.append(0)

    return (Variable(np.array(x1, dtype=np.int32)),
            Variable(np.array(x2, dtype=np.int32)),
            Variable(np.array(t,  dtype=np.int32)))

serialize/deserialize

枝葉だけど、モデルの保存と読み込みに結構嵌ったのでここにメモしておく。

モデルの保存は save_npz を使えば簡単にできる。 しかし save_npz ではモデルパラメータのみが保存され、語彙数、ベクトルの次元、ウィンドウサイズなどのハイパーパラメータを保存することができない。 これらのパラメータも一緒に保存したい場合は serializers.DictionarySerializer を使ってハイパーパラメータを追加でシリアライズすればいい。

def save_model(filename: str, model: Word2Vec, params: HyperParameters) -> None:
    serializer = serializers.DictionarySerializer()

    pickled_params = np.frombuffer(pickle.dumps(params), dtype=np.uint8)
    serializer("hyper_parameters", pickled_params)

    serializer["model"].save(model)

    np.savez_compressed(filename, **serializer.target)

モデルの読み込みも、保存のときと同様、 load_npz で一発でできるが、これだとモデルパラメータしかロードできない。 よって、serializers.NpzDeserializer を使ってまずハイパーパラメータを復元し、それを使ってモデルを作り、そのモデルを使ってモデルパラメータを読み込む。

def load_model(filename: str) -> Tuple[Word2Vec, HyperParameters]:
    with np.load(filename) as f:
        deserializer = serializers.NpzDeserializer(f)

        pickled_params = deserializer("hyper_parameters", None)
        params = pickle.loads(pickled_params.tobytes())  # type: HyperParameters

        model = Word2Vec(params.vocabulary_size, params.vector_dimension)
        deserializer["model"].load(model)

        return model, params

実験

Chainer 本で使っていたデータセットである ptb.train.txt を使ってモデルを学習してみた。

ハイパーパラメータがいくつかあるので、以下のすべての組み合わせ(12通り)を使って学習を行った。

  • モデル: SingleMatrix(WW'が等しいモデル), DoubleMatrix(WW'が等しくないモデル)
  • ベクトルの次元: 50, 100, 200
  • ウィンドウサイズ: 3, 5

類義語

ベクトルのコサイン類似度を使って"意味"が近い単語をアドホックに調べてみる。 (Search はコサイン類似度を使って似ている単語を探してくるヘルパークラス。名前が適当なのでもっといい名前がほしい。)

In [7]: model, params = load_model("./word2vec/results/word2vec_singlematrix_vd200_w5_ns5_bs100_epoch29.npz")

In [8]: s = Search(dataset.vocabulary, model)

In [9]: s.find_similar_words("ibm")
Out[9]: 
['ibm',
 'mainframes',
 'computer',
 'digital',
 'machine',
 'mainframe',
 'machines',
 'armonk',
 'software',
 'chips']

In [10]: s.find_similar_words("monday")
Out[10]: 
['monday',
 'late',
 'friday',
 'tuesday',
 'thursday',
 'wednesday',
 'sell-off',
 'opened',
 'close',
 'points']

In [11]: model, params = load_model("./word2vec/results/word2vec_doublematrix_vd200_w5_ns5_bs100_epoch29.npz")

In [12]: s = Search(dataset.vocabulary, model)

In [13]: s.find_similar_words("ibm")
Out[13]: 
['ibm',
 'digital',
 'mainframe',
 'armonk',
 'customers',
 'computer',
 'machine',
 'machines',
 'hewlett-packard',
 'computers']

In [14]: s.find_similar_words("monday")
Out[14]: 
['monday',
 'friday',
 'tuesday',
 'wednesday',
 'thursday',
 'advancers',
 'october',
 'quoted',
 'week',
 'trading']

一応それっぽい結果になったけど、これを見ただけだと上手く学習できているのかはいまいちわからない。

頻出語のプロット

f:id:nojima718:20170822022555p:plain

頻度が高い300単語のベクトルを t-SNE を使って2次元に落としてプロットした。 モデルは SingleMatrix, ベクトル次元200, ウィンドウサイズ5 のものを使った。

ラベルが重なりまくっていてかなり微妙な図だが、よく見ると would, should, might などの助動詞が近くに固まってたり、after-before や up-down、buy-sell などの対義語のペアがいたりしておもしろい。

DoubleMatrix, ベクトル次元200, ウィンドウサイズ5だとこんな図になった。こっちはさらにラベルの重なりが激しくて見づらい。

f:id:nojima718:20170822104836p:plain

国と首都の関係をプロット

各モデルにおいて、国と首都のベクトルがどのように位置しているかを主成分分析を用いて2次元平面にプロットした。 国と首都の相対位置がどのペアにおいてもだいたい似たベクトルになることを期待したが、残念ながらあまりそういう傾向は見られなかった。

f:id:nojima718:20170822010708p:plain

loss のプロット

エポックごとに各モデルのロスをプロットした。 ロスの計算には ptb.test.txt (訓練に用いていないデータ)を用いた。

上の図がウィンドウサイズが3であるときの図で、下の図がウィンドウサイズが5であるときの図。

f:id:nojima718:20170822011025p:plain

SingleMatrix のモデルでは、ベクトルの次元にかかわらずほぼ同じロスに収束している。 DoubleMatrix のモデルでは、ベクトルの次元が小さいほど小さいロスに収束しているが、SingleMatrix よりもロスが大きい。

つまり、このデータセットにおいては SingleMatrix のほうがよいモデルと考えてよさそう。

学習の過程

国と首都のベクトルの関係が、ひとつのモデルの学習の過程でどのように変化していくかをプロットした。

f:id:nojima718:20170822011729p:plain

ロスのグラフでは epoch 9 ぐらいになると全然ロスが減ってなかったけど、この図だと epoch 12 ぐらいまではわりと変化が見られる。 だから何だって言われたら何も言えないけど。

まとめ

Word2Vec を Chainer を使って実装してみた。

実験した範囲では本に書いてあったとおり、W = W' なモデルが良さそうだった。

Reference

Word2Vec のメモ その1

Chainer 本 を読みながら Word2Vec を Chainer で実装してみたので、その過程でわかったことをメモしておく。

注意: 素人なので完全に間違っているかもしれない。

Word2Vec

  • Word2Vec の目的は、各単語の 分散表現 を求めること。
    • 単語の分散表現とは、単語の意味を低次元(100次元とか)の密な実ベクトルで表したもの。
  • Word2Vec では、分散表現を求めるために Continuous Bag-of-Words (CBOW) か、skip-gram が利用できる。Chainer 本では skip-gram のみが解説されていたため、この記事でも skip-gram のみを扱う。

skip-gram

skip-gram とは、雑に言うと、「ある単語が与えられたときに、その周囲に現れる単語を予測する」という問題を考え、これを精度よく答えられるように各単語の分散表現ベクトルを学習してくモデル。

入力されたコーパス上の位置 { t } からオフセット { b } 以内にある単語群 { w_{t-b}, \dots, w_{t-1}, w_{t+1}, \dots, w_{t+b} } を位置 { t } における 文脈 { c_t } と定義する。 文脈という言葉を使って上の問題を言い換えると「位置 { t } にある単語 { w_t } が与えられたとき、その文脈 { c_t } を予測する問題」となる。

文脈内の各単語の条件付き独立性を仮定して、{ p( c_t \mid w_t ) } を以下のようにモデル化する。

{ \displaystyle
    p( c_t \mid w_t ) = \prod_{ c \in c_t } p( c \mid w_t )
}

図で書くとこんな感じ( b=2 のとき)

f:id:nojima718:20170821012826p:plain

さらに、{ p( c \mid w_t ) } を、{ c, w_t } の分散表現 { \vec{v}_c, \vec{v}_{w_t} } を使って、以下のように分散表現の内積のソフトマックスでモデル化する。(各単語の分散表現はパラメータであり、学習によって求める)

{ \displaystyle
    p( c \mid w_t ) = \frac{\exp( \vec{v}_c \cdot \vec{v}_{w_t} )}{\sum_{w' \in V} \exp( \vec{v}_{w'} \cdot \vec{v}_{w_t} )}
}

ただし、 V は語彙集合。

……と本には書いてあるんだけど、自分が元論文を読んだ感じでは単語のベクトル表現は2種類 (入力用、出力用) があり、{ w_t } は入力用のベクトル表現(=分散表現)でベクトル化し、{ c } は出力用のベクトル表現でベクトル化しているように見えた。 つまり、2種類のベクトル表現 { \vec{u}, \vec{v} } を使って以下のようにモデル化する。

{\displaystyle
    p( c \mid w_t ) = \frac{\exp( \vec{u}_c \cdot \vec{v}_{w_t} )}{\sum_{w' \in V} \exp( \vec{v}_{w'} \cdot \vec{v}_{w_t} )}
}

どっちがいいのかよくわからなかったので、両方実装して実験してみたので結果については後述する。

いずれにしても、これで文脈の条件付き確率が定義できたので、コーパスの対数尤度が以下の式で計算できる。

{
\begin{align*}
    \mathcal{L} &= \sum_{t=1}^T \log p( {\bf c}_t \mid w_t ) \\
                &= \sum_{t=1}^T \sum_{c \in {\bf c}_t} \log p( c \mid w_t ) \\
                &= \sum_{t=1}^T \sum_{c \in {\bf c}_t} \left( (\vec{v}_c \cdot \vec{v}_{w_t}) - \log \sum_{w' \in V} \exp( \vec{v}_{w'} \cdot \vec{v}_{w_t} ) \right)
\end{align*}
}

ただし、 T は訓練データのサイズ。

しかし、この尤度は { \log \sum_{w' \in V} \exp( \vec{v}_{w'} \cdot \vec{v}_{w_t} ) } の項があるため、学習のための計算コストが大きい。 元論文では語彙のサイズが70万ぐらいで、分散表現の次元が300であるため、70万×300 = 2億1000万回の演算が一つの入力ごとに必要になる。 入力データセットは、元論文の実験では300億単語あり、context size が 5 なので、1500億個の入力があることになる。 よって、 \mathcal{L} を計算するには 2億1000万×1500億=3150京回の演算が必要になる。 スーパーコンピューターなら計算できるかもしれないけど、普通はこのサイズの計算は無理。

よって、モデルをちょっと変えて、もっと計算しやすいモデルを作る。

次の問題を考える。

「位置  t の単語  w_t と、以下のいずれかの単語  c が与えられる。

  • 正例: 位置  t の文脈  c_t に属する単語
  • 負例: ランダムな単語

このとき、 c が正例なのか負例なのかを判別せよ1

つまり、もともとの問題は与えられた単語の近傍の単語の確率分布を求めていたが、新しい問題は近傍の単語とノイズとを区別する2クラス分類問題になっている。

まず、確率変数  D を導入して、答えが正例であるとき  D=1 とし、負例であるとき  D=0 とする。

 c w_t が与えられたときの  D の確率を以下の式でモデル化する。

{
    p(D=1 \mid c, w_t) = \sigma(\vec{v}_c \cdot \vec{v}_{w_t})
}

ただし、 \sigma(x)シグモイド関数 { 1 / (1 + \exp(-x)) }

要するに、分散表現の内積が大きい単語同士は近傍に出現しやすいというモデルになっている。 (ベクトルを2つ使うモデルだとちょっと違う解釈をしないといけない)

次にパラメータ(=分散表現)の学習のために誤差関数を作る。 Word2Vec ではクロスエントロピー誤差  H(p, q) = -p \log q - (1-p) \log (1-q) を用いる。

訓練データ  (w_t, c, D) に対するこのモデルのクロスエントロピー誤差は以下の式で与えられる。

{
\begin{align*}
    E &= H( D, p(D=1 \mid c, w_t) ) \\
      &= -D \log p(D=1 \mid c, w_t) - (1 - D) \log (1 - p(D=1 \mid c, w_t)) \\
      &= -D \log \sigma(\vec{v}_c \cdot \vec{v}_{w_t}) - (1 - D) \log ( 1 - \sigma(\vec{v}_c \cdot \vec{v}_{w_t}) )
\end{align*}
}

もともとの問題の対数尤度関数  \mathcal{L} と違って語彙集合の上を走る和がないことに注目してほしい。

 E の形をよく見ると、ロジスティック回帰の対数尤度関数とほぼ同じ式になっている。 なので、この手法は2クラス分類問題をロジスティック回帰を用いて解いていると言ってもいいのかもしれない。 (ロジスティック回帰の定義をよく知らないので間違ってるかも)

あとは訓練データを作れば学習できる。 訓練データは後述する ノイズ分布 とハイパーパラメータ  k \in \mathbb{N} を用いて以下のように作る。

  • 学習データ内の各位置  t について
    • 位置  t の文脈内の各単語  c について
      • 正例  (w_t, c, D=1) を訓練データに加える。
      • ノイズ分布から  k 個単語  c' をサンプルし、それぞれの  c' に対して負例  (w_t, c', D=0) を訓練データに加える。

ノイズ分布 は負例をサンプリングするための単語の確率分布で、Word2Vec では以下の分布を使うのが実験的によい結果を残しているらしい。

{ \displaystyle
p(w) = \frac{\mathrm{freq}(w)^{0.75}}{\sum_{w' \in V} \mathrm{freq}(w')^{0.75} }
}

ただし、 \mathrm{freq}(w) は、コーパスにおける単語  w の頻度。

このように負例をサンプリングして学習する手法を negative sampling と呼ぶらしい。

実装と実験

これを Chainer で実装していろいろ実験してみたけど、今日は疲れたので明日の記事で


  1. Chainer 本の説明では「文脈上の単語  c と (1)  w_t または (2) ランダムな単語 が与えられて、(1) か (2) かを区別する」問題を解いており、ここで紹介した問題は解いていない。元の論文では、(自分の理解が正しいとすると)この記事で紹介したほうの問題を解いていると思う。

ICFPC 2017 に参加した

ICFPC 2017 にチーム「ゲームセンターYAGI」として参加しました。メンバーは以下の4人です。

リポジトリhttps://github.com/seikichi/icfpc2017 です。

今年の問題

無向グラフが与えられます。 このグラフの頂点は都市鉱山のいずれかを表していて、辺はを表しています。

このゲームの目的は、川の領有権を主張(claim)して、各鉱山をできるだけ多くの都市と(川によって)繋げることです。

各プレイヤーは自分のターンに以下のいずれかの行動を行います。

  1. 他人によって claim されていない川をひとつ claim する。
  2. パスして「claimする権利」を1貯める。
  3. 溜まっている「claimする権利」を全て消費して、溜まっている権利数+1本の連続する川を claim する。

所定のターン数が経過するとゲームが終了し、各プレイヤーのスコアが以下のルールによって計算されます。スコアが最大であるプレイヤーが勝ちです。

  • 全ての都市 S について、以下を計算して合計する。
    • S からそのプレイヤーが claim した辺によって到達可能な各鉱山 M について、以下を計算して合計する。
      • M から S までの距離 (claim した辺以外も利用できるときの距離) の二乗

これ以外にも細かいルールが存在しますが、詳しくは問題文を参照してください。

タイムライン

1日目

  • 21:00 コンテストスタート。コンテストサーバが落ちる。

  • 21:08 問題文を入手する。

  • 21:40 問題文を大体読み終えたので、実装に入る。

    • 1人がシミュレータを書いて、1人が CI を整備し、残りの2人が AI を書くという分担になった。
  • 23:00 タイムアウト付きの spawn の実装に嵌まる。

    • fread(buffer, 1, sizeof(buffer), file) と書くべきところを fread(buffer, sizeof(buffer), 1, file) と書いていたのが原因。悲しい。
  • 00:04 ようやく spawn のテストが通る。ここからシミュレータの実装に入る。

    • JSON をパースしたりダンプするのが面倒くさい。心を無視にして実装していく。
    • シェルスクリプトで AI を作るとタイムアウトしたときに親プロセスしか殺されず、子プロセスが残る問題を見つける。プロセスグループを使って kill するように spawn を修正する。
  • 03:30 セットアップフェイズまで実装完了。ライブラリの整備も進む。

  • 03:40 寝る。

  • 10:30 起きる。

  • 11:20 シミュレータに GamePlay フェイズを実装。

  • 11:30 昼飯を食べに行く。

  • 15:00 Scoring フェイズを実装。これで一応シミュレータが完成。しかし、シミュレータを作っている間にプロトコルが破壊的に変更されていたので、これに対応しないといけない。

    • 今までは入力を一個入れると出力が一個あるようなインターフェイスだったので、普通の spawn のインターフェイスで何とかなったんだけど、新しいインターフェイスでは、ハンドシェイクを先にやる仕様になったのでこれが使えなくなった。
    • 破壊的に変更を "small change" って表現するのやめてほしい。
  • 17:00 新しいプロトコルに対応するために spawn とシミュレータを書き直した。

  • 20:00 貪欲AI (毎ターン、合法手の中からスコアが最大になるような手を選ぶ) を書いた。

  • 21:00 貪欲AI に枝刈りを追加して高速化。lightning にはこれを提出した。

2日目

  • 21:30 飯を食べに行く。鯖が美味しかった。

  • 01:00 寝た。

  • 11:00 起きた。寝すぎ。

  • 12:00 昼飯を食べに行く。ラーメン。

  • 14:00 鉱山に近づいていくAIを作った。

    • 最も近い鉱山への距離 (他人が claim した辺は通れないとして考える) を最小化する手を選ぶ。
    • 何をやっても鉱山に近づけない場合は貪欲AIにフォールバックする。
  • 16:30 タイムアウトしそうになったら弱いAIにフォールバックする仕組みを作った。

    • しかしこれはちゃんと動かないことが後になって判明する。
  • 18:00 夕食を食べに行った気がする。

3日目

  • 20:00 初手の選び方を改善した。

    • できるだけ強い連結成分を選び、連結成分の中心に近い鉱山を選ぶ、というヒューリスティック
    • 連結成分の強さは、頂点数×鉱山数で表すことにした。
    • 「連結成分の中心への近さ」は「その鉱山から他の都市への距離の和」で測ることにした。この和が小さいほど中心性が高い。
  • 23:30 貪欲に次の手を選ぶ際に、到達したときにスコアが最大となる未達都市への距離を最小にするような手を優先するようにした。

    • 今までは次の1手のスコアしか考慮していなかったので、将来的に遠くまで行けるパスがあっても軽視されることがあった。
  • 仕事のため、自分の ICFPC ここで終了。

所感

今年は仕事のため月曜日は参加できませんでした。1日少ないとかなり短く感じます。 (普通なら有給休暇を取ればいいんだけど、今年は4週間ぐらい入院したので有給休暇をほとんど使い切ってしまっていたので)

1日目の最初の方はずっと spawn に嵌っていて、2日目もタイムアウトの実装で時間を使ってしまったので、ここらへんをライブラリで済ませられればもっと AI の実装に時間を使えた気がします。 シミュレータは別に提出するプログラムじゃないので、もっと手を抜いて実装すべきだったのかも。

AI の出来はいつもどおりっていう感じです。そろそろもっと高度な探索とかをしたいなー。

あと、そろそろ C++ にも飽きてきたら別の言語が使いたい。 来年は Rust で。

Ubuntu 16.04 で pystan を動かす

MCMC サンプラー Stan を Python から呼び出すためのライブラリ pystan を使ってみようとしたら ABI 問題のせいでちょっと嵌ったのでここにメモしておく。

Anaconda3-4.3.0-Linux-x86_64 をインストールした。 そして pystan を pip でインストールした。

pip install pystan

この状態で pystan を使ってみたところ、モデルをコンパイルするところで以下のようなエラーが出た。

>>> stanmodel = StanModel(file='model/model5-3.stan')
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-5-4435d7d1dfc2> in <module>()
----> 1 stanmodel = StanModel(file='model/model5-3.stan')
      2 stanmodel

/home/nojima/anaconda3/lib/python3.6/site-packages/pystan/model.py in __init__(self, file, charset, model_name, model_code, stanc_ret, boost_lib, eigen_lib, verbose, obfuscate_model_name, extra_compile_args)
    313                 os.dup2(orig_stderr, sys.stderr.fileno())
    314 
--> 315         self.module = load_module(self.module_name, lib_dir)
    316         self.module_filename = os.path.basename(self.module.__file__)
    317         # once the module is in memory, we no longer need the file on disk

/home/nojima/anaconda3/lib/python3.6/site-packages/pystan/model.py in load_module(module_name, module_path)
     48         pyximport.install()
     49         sys.path.append(module_path)
---> 50         return __import__(module_name)
     51     else:
     52         import imp

ImportError: /tmp/tmppaa8t0sa/stanfit4anon_model_d6577da659c87ba489087107d1a725f1_1731976423115490565.cpython-36m-x86_64-linux-gnu.so: undefined symbol: _ZTVNSt7__cxx1118basic_stringstreamIcSt11char_traitsIcESaIcEEE

これは、最近 libstdc++ の ABI の変更があり、Ubuntu 16.04 デフォルトの gcc に同梱されている libstdc++ は新しい ABI のものだが、Anaconda に同梱されているライブラリは古い ABI を使っており、うまくリンクできないのが理由っぽい。

ということで、Anaconda が利用しているバージョンである gcc-4.8 をインストールする。 幸いなことに Ubuntu 16.04 には gcc-4.8 パッケージが存在しているため、それを apt-get するだけでインストールできる。

sudo apt-get install gcc-4.8 g++-4.8

この方法でインストールした gccgcc-4.8 という名前で参照できる。gcc コマンドは依然として 5.4 であることに注意。

次に、Cython が gcc-4.8コンパイルに使うようにしなければならない。このためには CC 環境変数gcc-4.8 を指定すればよいらしい。

するとうまくコンパイルできた。

In [1]: from pystan import StanModel

In [2]: import os

In [3]: os.environ['CC'] = 'gcc-4.8'

In [4]: stanmodel = StanModel(file='model/model5-3.stan')
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_d6577da659c87ba489087107d1a725f1 NOW.

In [5]:

記念にトレースプロットを貼っておく。

f:id:nojima718:20170304173332p:plain

明日使えない Linux の capabilities の話

(この記事は KMC アドベントカレンダー 2016 の3日目の記事です)

はじめに

みなさん以下のようなことで困ったことはないでしょうか?

ポート80を listen したいけど特権ポートなので、一般ユーザの権限で動くデーモンでは bind できない。

1024未満のポートは特権ポートと呼ばれ、一般ユーザの権限では bind することはできません。 この問題の解決策を考えてみます。

(なお、長々と説明を書いていますが、結論だけ知りたい人は一番下だけ読んで下さい)

root で起動

まず、root であれば特権ポートを自由に bind できるので、root で対象デーモンを起動すれば、特権ポートを bind できます。 しかし、デーモンを root として動作させるのは一般にリスクが大きいです。 もしそのデーモンに脆弱性があった場合、root 権限を悪用される可能性があるわけです。

したがって、このやり方は却下です。

プロセス分割

プロセスを分割することでこの問題に対処する場合もあります。 nginx は、まず root として起動し、特権ポートを bind します。 そして、子プロセスを fork し、それらを一般ユーザとして実行します。 この際に、listening socket を open したままの状態にしておくことで、子プロセスでもその socket を利用することができます。

しかし、プロセス分割は結構大掛かりです。 Java や Go みたいにランタイムの都合上 fork できない場合もあります。 できればプログラム本体には手を入れないで対処できる方法が欲しいわけです。

capability

今までの話をまとめると

  • root は使いたくないけど
  • 特権ポートを bind したい

ということになります。 今回は Linux の capabilities という機能を用いてこれを実現していきます。

capabilities とは、root が持っている特権をいくつかのグループに分割し、それぞれを独立に enable, disable できるようにしたものです。

capabilities には例えば以下のようなものがあります。

  • CAP_DAC_OVERRIDE
  • CAP_DAC_READ_SEARCH
  • CAP_KILL
    • シグナルを送るときの権限チェックをバイパスする。
  • CAP_NET_BIND_SERVICE
    • 特権ポートにソケットをバインドできる。
  • CAP_SYS_TIME
    • システムの時刻を設定できる。

他にもいっぱい capability はあるので、興味のある人は man capabilities してみましょう。 今回の用途では CAP_NET_BIND_SERVICE を利用すればよさそうです。

File capabilities

プロセスに capabilities を持たせるひとつの方法として、実行ファイルに対して capabilities を設定する方法があります。 この方法を使うと、そのファイルを 誰が実行しても 指定した capabilities を持ってプロセスが起動するようになります。

ここでは例として、パーミッションを無視して任意のファイルを見れる超脆弱 cat を作ってみます。

# まず /bin/cat を自分用にコピー
$ cp /bin/cat ./mycat

# CAP_DAC_READ_SEARCH を付与する
$ sudo setcap cap_dac_read_search=eip ./mycat

# 確認
$ getcap ./mycat
mycat = cap_dac_read_search+eip

# まず普通の cat で /etc/shadow を見ようとしてみる (当然見れない)
$ cat /etc/shadow
cat: /etc/shadow: Permission denied

# 次に mycat で /etc/shadow を表示!!
$ ./mycat /etc/shadow
root:!:17042:0:99999:7:::
daemon:*:17001:0:99999:7:::
bin:*:17001:0:99999:7:::
...

このように catCAP_DAC_READ_SEARCH を付与することでパーミッションを無視して任意のファイルを表示できるようになりました。

sudo setcap cap_dac_read_search=eip ./mycatmycat に file capabilities を付加している部分です。 eip は立てるフラグの種類をしていしています、と言ってもこの時点では意味不明だと思いますが説明すると長くなるので省略します。

さて、冒頭の問題に戻りましょう。

特権ポートを bind する権限を与えたいので、CAP_NET_BIND_SERVICEsetcap を使って対象サービスの実行ファイルに付与すれば良さそうです。

しかし、この方法には問題があります。

file capabilities を使った方法では 誰が実行しても その特権を持ってしまいます。 一般ユーザは通常何の特権の持っていないので、一般ユーザにとっては setcap されたファイルを実行すると特権が増えることになります。 一般に、特権が増えるような系を作ってしまうと、脆弱性のリスクが増えるので、できるだけ避けたいです。

ということで、file capabilities を使わずに、デーモンに CAP_NET_BIND_SERVICE を付与する方法を考えてみます。

Ambient capabilities

ここまでをまとめると、

  • root として動いている親プロセス(今回の場合は init)が存在し、
  • この親プロセスから一部の capabilities (今回の場合は CAP_NET_BIND_SERVICE)だけを継承して、非rootユーザの子プロセスを spawn したい。
  • File capabilities は利用しない。

という要件になります。

実は、最近まで capabilities の継承ルールがクソなせいで、file capabilities を使わないと非rootユーザに capabilities を付与することはできませんでした(多分)。 しかし、Linux 4.3 で Ambient capabilities という機能が追加され、これが簡単にできるようになりました。

Ambient capabilities の説明に行く前に、簡単にプロセスの capabilities を説明します。

各プロセス(厳密にはスレッド)は、Permitted, Effective, Inheritable の3つのフラグセットを持っています。

  • Permitted, Effective, Inheritable はそれぞれ capability の集合を表すフラグセットです。
  • Effectiveカーネルによってパーミッションチェックをされるときに使われる capability の集合を表します。
  • PermittedEffectiveInheritable の superset を表します。Permitted に含まれていない capability を EffectiveInheritable に加えることはできません。

そして、Inheritableexecve(2) で他の実行ファイルを起動するときに継承される capability の集合を表し……ていたら何も問題なかったのですが、Inheritable に含まれる capabilities は条件が揃わないと子プロセスの PermittedEffective に継承されません。 特に、今回のように「非rootユーザ」で「file capabilitiesが設定されていないファイルを実行」する場合、PermittedEffective は常に空になります。

そこで、ambient capabilities を利用します。

ambient capabilities とは、簡単に言うと「子プロセスに継承される capabilities」を表すフラグセットです。 より厳密に言うと、これから実行するファイルに setuid ビットや setgid ビットがセットされておらず、かつ file capabilities も設定されていない場合、ambient capabilities が子プロセスに継承されます。

つまり、子プロセスに継承したい capability があるときは、ambient capabilities にそれをセットしてから exec すればいいわけです。

どうして最近まで存在しなかったのか疑問に思うぐらい、シンプルで自然な機能ですね。

systemd

以上より、冒頭の問題を解決するためには、「init が対象のデーモンを起動するときに ambient capabilities に CAP_NET_BIND_SERVICE を加えておけばよい」ということになります。

systemd はまさにその機能を持っているので、以下の指定をユニットファイルに書くことで簡単に CAP_NET_BIND_SERVICE を持った状態でデーモンを起動することができます。

AmbientCapabilities=CAP_NET_BIND_SERVICE

おわりに

一方ロシアは8080を使った。

Protocol Buffers が本当に遅いのか実際に確かめてみた

Protocol Buffers で検索すると Protocol Buffersは遅い という MessagePack 作者による2008年の記事が未だに上位に来る。 一方で、Protocol Buffersは遅いのか という反論記事も見つかる。 一体遅いのか速いのかどっちなんだ!!ということで、ベンチマークを取ってみた。

2016年8月現在では、Protocol Buffer の最新バージョンは 3.0.0 であり、MessagePack の C++ バインディングの最新バージョンは 2.0.0 なので、これらのバージョンを使ってベンチマークを取ることにした。

ベンチマーク

元の記事を踏襲して、以下の4つのメッセージを使ってベンチマークを行った。

  • Test1: 2つの符号無し整数
  • Test2: 2つの符号付き整数
  • Test3: 8KBのバイト列
  • Test4: Test1 + Test2 + Test3

Test1 と Test2 は 223 個、Test3 と Test4 は 216 個のメッセージをシリアライズして所要時間を計測し、再びそれをデシリアライズして所要時間を計測した。 外乱の影響を抑えるため、それぞれのワークロードを10回ずつ実行して所要時間の最小値を最終結果とした。

環境:

Test1: 2つの符号無し整数

Serialize Deserialize Size
Protobuf 107 msec 145 msec 32 MB
MessagePack 100 msec 1929 msec 24 MB

シリアライズにかかる時間はほとんど同じだったが、MessagePack のデシリアライズが Protobuf の10倍以上遅い。

Test2: 2つの符号付き整数

Serialize Deserialize Size
Protobuf 112 msec 147 msec 32 MB
MessagePack 102 msec 1954 msec 24 MB

Test1 とほぼ同じ結果になった。 上で挙げた記事では Protobuf の符号付き整数のシリアライズ後のサイズが非常に大きいと書いてあるが、sint を使えば MessagePack と同じサイズになる。 (Protobuf と MessagePack のサイズ差はタグの有無によるもの)

Test3: 8KBのバイト列

Serialize Deserialize Size
Protobuf 162 msec 42 msec 512 MB
MessagePack 121 msec 79 msec 512 MB

バイト列のシリアライズは MessagePack の方が速かった。逆にデシリアライズは Protobuf の方が速い。

Test4: Test1 + Test2 + Test3

Serialize Deserialize Size
Protobuf 165 msec 66 msec 513 MB
MessagePack 123 msec 83 msec 513 MB

Test3 と同じくシリアライズは MessagePack の方が速く、デシリアライズは Protobuf の方が速いという結果になった。

考察

ベンチマークの結果から、バイト列のシリアライズは MessagePack の方が速いが、デシリアライズや整数のシリアライズは Protobuf の方が速いことがわかった。 したがって、2016年の現在においては Protobuf は別に遅くないと思ってよさそう。

ベンチマークプログラム

https://gist.github.com/nojima/005cf04bfa35a1fb971adc43b54abbef