2019-02-24

Pytorch：テキストの前処理（torchtext）④

PyTorch torchtext

はじめに

torchtextの使い方に関するメモ。

入力の素性をカスタマイズしたい場合について

例えば、各系列に付与したカラムを入力に含めたい場合
（0,1は系列の何らかの情報）

これをどうハンドリングするかについて少し考えた。

簡単な方法

多分、一番簡単な方法は以下のように系列カラムと同じように処理するだけ。

TEXT = data.Field()
LABEL = data.Field()
train = datasets.SequenceTaggingDataset(
        path=data_dir,  
        fields=[('text', TEXT),
                ('label', LABEL)])

TEXT.build_vocab(train)
LABEL.build_vocab(train)

これだけで基本的には問題ない。
ただ、思ったのは素性が数値の場合、わざわざ「辞書作り」する必要があるのかってこと。
別に辞書作りをしても全く問題ないが、個人的にはそのまま入力できないかなって思った。
で、どうするのが良いか考えてみた。

そのまま処理する方法

1.入力の数値変換

結論としては、そのまま処理できないので、工夫が必要。
まず、初めにデータをDatasetで読み込んだ場合、全ての入力は「文字列」として扱われる。

train.examples[0].label
# ['1', '0', '1']

そして、辞書を作らずにそのまま処理した場合、バッチ化のところで「文字は扱えないよ」って怒られる。

なので、まずはここで入力を数値変換する必要がある。
Datasetでハンドリングできるかもしれないが、ここでは Fieldのpreprocessingで変換させる。

def to_int(x):
    return list(map(int, x))
LABEL = data.Field(use_vocab=False, preprocessing=to_int)

これで、入力情報を数値変換できる。

train.examples[0].label
# [1, 0, 1]

2.pad値の指定

続いて、バッチのpadding時の値を指定する。
これはデフォルトの場合、<unk>となっているので、このままではまたも怒られる。
そこで、paddingの値をpad_tokenで明示する。

LABEL = data.Field(sequential=True, use_vocab=False, preprocessing=to_int, pad_token=-1)

※指定する値は入力で使用しない値をしておくのが無難。

バッチ処理までをまとめて書くとこんな感じ。

from torchtext import data, datasets

data_dir = 'test_ner.txt'

def to_int(x):
    return list(map(int, x))

TEXT = data.Field()
LABEL = data.Field(sequential=True, use_vocab=False, preprocessing=to_int, pad_token=-1)
train = datasets.SequenceTaggingDataset(
        path=data_dir,  
        fields=[('text', TEXT),
                ('label', LABEL)])

TEXT.build_vocab(train)
rain_iter = data.BucketIterator(dataset=train, batch_size=2, device=-1, 
                                 repeat=False)

for i, train_batch in enumerate(train_iter):
    print('text  : \n', train_batch.text)
    print('label : \n', train_batch.label)

"""
text  : 
 tensor([[7]])
label : 
 tensor([[0]])
text  : 
 tensor([[5, 2],
        [6, 3],
        [1, 4]])
label : 
 tensor([[0, 1],
        [1, 0],
        [9, 1]])
"""

全く持って自己満足に過ぎないが、備忘録として残しておく。

2018-11-29

Stan：LDA

LDA Stan

はじめに

自然言語処理の領域では広く知られいるLDA(Latent Dirichlet Allocation)について復習する。

LDAはトピックモデルの１種であり、文書がどのようなトピックから構成されているかを推論するモデル。
推論するパラメータは以下の２つ。

トピック分布：文書ごとのトピック構成比率
単語分布：トピックごとの単語比率

トピックモデルに関する理解はこの１枚に尽きると思う。
f:id:kento1109:20180627110029p:plain
Fast and Scalable Algorithms for Topic Modeling | Center for Big Data Analyticsより引用

後、日本語でのLDAの説明としては視覚的にも以下が分かりやすかった。
LDA for Pokemon analysis | haripo.com

モデリング

数式によるLDAはググれば色々出てくるのでここでは割愛する。今回は「尤度計算」の具体例のみを数式で紹介する。これが分かれば、Stanコードも理解できると思う。

$\begin{eqnarray} L &=& \prod_{d=1}^D\prod_{w=1}^W \sum_{k=1}^K　p(\theta_{d,k})\times p(\phi_{k,w}) \\ \log L &=& \sum_{d=1}^D\sum_{w=1}^W \log\left(\sum_{k=1}^K \{p(\theta_{d,k})\times p(\phi_{k,w})\}\right) \end{eqnarray}$
各記号の説明は以下の通り。

トピック分布（ $\Theta\in\mathbb{R}^{d\times k}$ ）
単語分布（ $\Phi\in\mathbb{R}^{k\times w}$ ）

$d$ は文書数、 $k$ はトピック数、 $w$ は単語数を表す。
これは「各文書内に出現する各単語の尤度（出現確率）」を計算している。

ちなみに単語の生成確率は文書 $d$ の $n$ 番目の単語のトピック割り当てを意味する潜在変数 $z_{d,n}\in\{0,1\}^K$ を用いると以下のように表せた。
$\begin{eqnarray}p(w_{d,n})=\prod_{k=1}^K p(\phi_k)^{z_{d,n,k}} \end{eqnarray}$
しかし、Stanでは、整数型の潜在変数が使えないので、上式のようにしてsumming outしている。

対数尤度の計算

まず、以下の文書のトピックを考える。

play soccer

全ての文書が「音楽」か「スポーツ」の何れかに関するものであるとき、この文書はどちらに分類されるか。
常識で考えると「スポーツ」トピックから生成された（「スポーツ」トピックの方が混合比が大きい）と考えるのが自然だろう。

ただ、そんな常識を機械は知らない。
では、どう考えればよいか。

それを考えるため、この文書のトピック分布（ $\theta_{d,k}$ ）及び単語分布（ $\phi_{k,w}$ ）が事前にディレクレ分布から以下のように生成されたとする。

トピック分布（ $\theta_{d,k}$ ）

music	0.5
sports	0.5

単語分布（ $\phi_{k,w}$ ）

	play	soccer
music	0.01	0.01
sports	0.01	0.01
guitar	0.01	0.01
baseball	0.01	0.01
the	0.01	0.01

※guitar,baseball,theに関しては後で言及する。

この文書（ $d$ ）における尤度は以下のように計算できる。
$\begin{eqnarray} L_1 &=& \sum_{w=1}^W \log\left(\sum_{k=1}^K \{p(\theta_{d,k})\times p(\phi_{k,w})\}\right)\\ &=& \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,play} \\ + \\ \theta_{d,sports}\times \phi_{sports,play}\\ \end{array} \right)+ \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,soccer} \\ + \\ \theta_{d,sports}\times \phi_{sports,soccer}\\ \end{array} \right)\\ &=& \log\left( \begin{array}{ccc} 0.5\times 0.01 \\ + \\ 0.5\times 0.01\\ \end{array} \right)+ \log\left( \begin{array}{ccc} 0.5\times 0.01 \\ + \\ 0.5\times 0.01\\ \end{array} \right)\\ &=& \log(0.01)+\log(0.01)=-4.0 \end{eqnarray}$

次に尤度が大きくなるよう以下のように単語分布をサンプリングされたとする。（トピック分布は固定）

単語分布（ $\phi_{k,w}$ ）

	play	soccer
music	0.01	0.01
sports	0.01	0.1
guitar	0.01	0.01
baseball	0.01	0.01
the	0.01	0.01

この場合、尤度は以下のように更新される。
$\begin{eqnarray} L_1 &=& \sum_{w=1}^W \log\left(\sum_{k=1}^K \{p(\theta_{d,k})\times p(\phi_{k,w})\}\right)\\ &=& \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,play} \\ + \\ \theta_{d,sports}\times \phi_{sports,play}\\ \end{array} \right)+ \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,soccer} \\ + \\ \theta_{d,sports}\times \phi_{sports,soccer}\\ \end{array} \right)\\ &=& \log\left( \begin{array}{ccc} 0.5\times 0.01 \\ + \\ 0.5\times 0.01\\ \end{array} \right)+ \log\left( \begin{array}{ccc} 0.5\times 0.01 \\ + \\ 0.5\times 0.1\\ \end{array} \right)\\ &=& \log(0.01)+\log(0.55)=-1.3 \end{eqnarray}$
単語分布を更新したことで尤度が大きくなったことが確認できた。
※当然だが、機械は「スポーツ」トピックかは知らないので、逆になることもある。（計算結果、人間がこれは「スポーツ」と判断するに過ぎない。）

次に、単語分布を固定してトピック分布を以下のようにサンプリングし直す。

トピック分布（ $\theta_{d,k}$ ）

music	0.3
sports	0.7

尤度は以下のように更新される。
$\begin{eqnarray} L_1 &=& \sum_{w=1}^W \log\left(\sum_{k=1}^K \{p(\theta_{d,k})\times p(\phi_{k,w})\}\right)\\ &=& \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,play} \\ + \\ \theta_{d,sports}\times \phi_{sports,play}\\ \end{array} \right)+ \log\left( \begin{array}{ccc} \theta_{d,music}\times \phi_{music,soccer} \\ + \\ \theta_{d,sports}\times \phi_{sports,soccer}\\ \end{array} \right)\\ &=& \log\left( \begin{array}{ccc} 0.3\times 0.01 \\ + \\ 0.7\times 0.01\\ \end{array} \right)+ \log\left( \begin{array}{ccc} 0.3\times 0.01 \\ + \\ 0.7\times 0.1\\ \end{array} \right)\\ &=& \log(0.01)+\log(0.73)=-1.1 \end{eqnarray}$
さらに尤度が大きくなったことが確認できた。
トピック分布・単語分布は合計１の制約があるので、ある単語の出現確率を大きくすると別の単語の確率を小さくする必要がある。
では、どのように単語の割合を調整するのか。

例えば、最初の文書が

play soccer and baseball

だったとする。
この場合、sportsトピックにおけるsoccerの出現確率の増加分をどのように調整するのが良いか。簡単のため、baseball,guitar,theの何れかの出現確率を小さくして調整するとする。
baseballの確率を小さくした場合、この文書の尤度が小さくなるので、baseballは調整しない。一方、guitar,theの確率を小さくしてもこの文書の尤度には何ら影響はない。では、guitar,theであればどちらでもよいのか。結論から言うと、一般的にはguitarの確率を小さくようにサンプリングされる（勿論、これはコーパスによる）。この文書だけを見ると、どちらも尤度に影響を与えない。しかし、コーパス全体で見た時に、guitar,theはどちらの単語の出現確率が高いか。一般的にはtheの方が高い。なので、theの出現確率を下げると、尤度が小さくなる文書の数が多くなる。一方、guitarはそれほど多くの文書に出現するとは考えられないので尤度が下がる文書の数が限定される。文章全体の尤度を大きくすることを考えた場合、guitarの確率を小さくするのが良いと考えられる。
※baseballは調整しないと書いたが、これは一般的は共起性での話であり、他の文書でbaseball以上にguitarがsoccerと共起する場合、baseballの確率が小さくなるよう調整するかもしれない。

サンプリングを文書全体で繰り返しすることで、同じトピックに高確率で出現する単語同士（soccerとbaseballなど）は「共起性」があると考えることが可能となる。
また、同じようなトピック分布をもつ文書同士は類似した文書である可能性が高いと考えることが可能となる。

さいごに

さいごに、Stanコードを載せておく。
（Stan モデリング言語: ユーザーガイド・リファレンスマニュアルとほとんど同じであるが・・）

data {
  int<lower=2> K;               // num topics
  int<lower=2> V;               // num words
  int<lower=1> M;               // num docs
  int<lower=1> N;               // total word instances
  int<lower=1,upper=V> w[N];    // word n
  int<lower=1,upper=M> doc[N];  // doc ID for word n
  vector<lower=0>[K] alpha;     // topic prior
  vector<lower=0>[V] beta;      // word prior
}
parameters {
  simplex[K] theta[M];   // topic dist for doc m
  simplex[V] phi[K];     // word dist for topic k
}
model {
  for (m in 1:M)  
    theta[m] ~ dirichlet(alpha);  // prior
  for (k in 1:K)  
    phi[k] ~ dirichlet(beta);     // prior
  for (n in 1:N) {
    real gamma[K];
    for (k in 1:K) 
      gamma[k] = log(theta[doc[n], k]) + log(phi[k, w[n]]);
    target += log_sum_exp(gamma); // likelihood;
  }
}

以下のようにしてキックできる。

data <- read.csv("lda.csv")
stan.data = list(K=10, V=max(data$w), M=max(data$d), N=nrow(data),
                 w=data$w, doc=data$d, alpha=rep(0.1,10), beta=(0.1,max(data$w))
stan.fit <- stan(file="lda.stan", data=dat)

あまりよい推定結果は得られないかもしれないが、とりあえず動作検証にはなると思う。
テストデータは以下に置いた。
github.com

2018-10-08

Pretrained language modelsを理解する

論文 Transfer Learning Language Model

はじめに

少し前に「A Review of the Neural History of Natural Language Processing」というまとめ記事を見ました。
blog.aylien.com

自然言語処理の発展の歴史みたいなもので、近年のブレイクスルーをまとめてくれています。

2018年の注目はELMoに代表される「Pretrained language models」の発展であり、これによりSOTAを達成したとの論文が挙がっています。

ELMoはこちらで分かりやすく解説されています。
kamujun.hatenablog.com

ここでは、自分なりに「Pretrained language models」をまとめようと思います。

Language Model（言語モデル）

言語モデルは文章（単語の並び）に関する確率分布を評価するものです。
ある単語の並びが自然であればその文章の確率は高くなり、不自然であれば低くなります。
言語モデルが学習するのは「自然な単語の並び」です。

数式を使って言語モデルを書いてみます。ここでは、 $w_1,...,w_n$ という $n$ 個の単語からなる文章について考えます。このとき、 $w_1,...,w_n$ という並びで単語が出現する確率は同時確率 $P(w_1,...,w_n)$ で表されます。これは条件付き確率で書くと、
$\begin{eqnarray} P(w_1,...,w_n)=\prod_{t=1}^nP(w_t|w_1,...,w_{t-1}) \end{eqnarray}$
となります。
これは後述する言語モデルと区別するため、「forward LM」（前向き言語モデル）と呼ばれることが多いです。

最近の言語モデルでは、前向き・後ろ向きの両方の言語モデルを利用することが一般的です。
後ろ向き言語モデルは「backward LM」と呼ばれ、下記のような式で表されます。
$\begin{eqnarray} P(w_1,...,w_n)=\prod_{t=1}^nP(w_t|w_{t+1},...,w_n) \end{eqnarray}$

イメージとしてはこんな感じです。
f:id:kento1109:20181008110529p:plain

2つ合わせて「biLMs」（双方向言語モデル）と呼びます。

「are」という単語の出現確率の計算に「前の単語の並び」だけではなく、「後ろの単語の並び」も利用します。

双方向の並びを学習することで、中間層の重みが言語モデルを最適化するように調整されます。
最終的にはその重みを結合させます。

f:id:kento1109:20181025164901p:plain

※上の２つはThe Bidirectional Language Model – Motoki Wu – Mediumより引用

「word2vec」などの事前学習で獲得した分散表現を入力層の初期パラメータとすることで自然言語処理のタスクの精度が向上することが知られています。
ただし、入力層に渡される文章の分散表現は単語を順番に並べただけであり、当然ですが文脈情報は含まれていません。
一方、双方向の言語モデルで学習された分散表現の場合、文脈情報が考慮されております。

このあたりが「文脈情報を考慮した分散表現」としてSOTAに貢献しているのかなとも思いました。

ELMo

ELMoで学習した単語の分散表現は以下のように目的のタスクに結合させます。
f:id:kento1109:20181025174941p:plain

※http://ai2-website.s3.amazonaws.com/publications/semi-supervised-sequence.pdfより引用

通常の単語の分散表現は、文章中の単語を分散表現に変換したものの集まりです。
一方、biLMsの分散表現は、文章を言語モデルに通して得られたものです。
上記の例だと、単語「York」の分散表現は、forward LMを通して、「New」の単語ベクトルの情報も含んでいます。また、backward LMを通して、「is located ...」の単語ベクトルの情報も含んでいます。

ELMoの応用

まず、タスクの対象言語が日本語の場合は日本語の言語モデルを使う必要があります。そのため、ELMoなどの既存の学習済のモデルをそのまま使うことは出来ません。
また、タスクがドメイン依存する場合、Wikipedia等のコーパスで学習した言語モデルがどの程度有効かは分かりません。
次回はドメインに依存するタスクでの言語モデルについてまとめたいと思います。

参考

以下の解説ページが分かりやすかったです。

2018-10-03

Pytorch：テキストの前処理（torchtext）③

PyTorch torchtext メモ

はじめに

torchtextの使い方メモ第三弾。

前回の内容は下記参照
kento1109.hatenablog.com

今回の内容は１つだけ。
POSやNERなどのTaggingを考える場合、どのようにtorchtextで読み込めばよいか。

前回まではtorchtextでデータをファイルから読む際、想定されているのは（感情分析等の場合）以下のような形式だった。

I love this moive,1

ただし、学習の目的がTaggingの場合、入力データのフォーマットがそもそも異なる。
NERの場合、だいたいこんな感じだ。

Tom B-PER
is O
my O
friend O

New B-LOC
York I-LOC
is O
always O
crowded O

行ごとに単語が並んでおり、それぞれにタグが付与されている。
文章の区切りは、CoNLLのデータセットの場合、空行で行われる。

こんなフォーマットの場合、どのように処理すればよいだろう。
結論から言うと、SequenceTaggingDatasetを使う。以上だ。

SequenceTaggingDataset

先ほどのデータの場合、以下のようにして読み込む。

from torchtext import data, datasets

TEXT = data.Field()
LABEL = data.Field()

train = datasets.SequenceTaggingDataset(
        path='test.tsv',
        fields=[('text', TEXT), ('label', LABEL)])

すると、以下のようにセンテンス単位となる。

print(vars(train[0]))
# {'text': ['Tom', 'is', 'my', 'friend'], 'label': ['B-PER', 'O', 'O', 'O']}

必要に応じて

TEXT = data.Field(init_token='<bos>', eos_token='<eos>')
LABEL = data.Field(init_token='<bos>', eos_token='<eos>')

とすると、文の先頭と末尾に記号をつけることが出来る。

センテンスの単位であるが、github.com
にあるSequenceTaggingDatasetクラスをみれば分かるが、行が空白の場合、それまでの単語・ラベルをセンテンスとしてまとめる処理をしている。（31行目あたり）

datasetsさえ作れれば、後の処理は前回と同じ。
前回と同じ書き方で辞書も作れる。

TEXT.build_vocab(train)
vocab = TEXT.vocab
vocab.itos
['<unk>',
 '<pad>',
 'is',
 'New',
 'Tom',
 'York',
 'always',
 'crowded',
 'friend',
 'my']

大した内容ではないが、少し嵌ったのでメモした。

2018-09-30

Neo4jを触ってみる

Neo4j

はじめに

自然言語処理を勉強していると、オントロジーやシソーラスってのが出てくる。
オントロジーやシソーラスを勉強すると、グラフ理論ってのが出てくる。
グラフ理論を勉強すると、グラフDBとして有名な「Neo4j」が出てくる。

ということで、今回はNeo4jの基本的なことを勉強したのでまとめようと思う。

以下の内容は公式サイトのドキュメントより引用したものが、ほとんどなので詳しく知りたい場合はそっちを見るのが良い。
neo4j.com

Property Graph Model

Neo4jは、Property Graph Modelと呼ばれている。これは、各ノードがプロパティを持ち、ノード情報を格納するモデル。
グラフ構造の基本形はこんな感じ。
f:id:kento1109:20180930115848p:plain

各構成要素を簡単に説明をする。

ノード

しばしば「エンティティ」を表現される。
1つの実体を表す構成単位みたいなもの。
f:id:kento1109:20180930120002p:plain

プロパティ付きで単一をノードを作成するには、

CREATE (n {name: {key1:value1, key2:value2, ...}})

と書く。

関係

ノード間をつなぐ要素。関係を定義してこそのグラフDB。
f:id:kento1109:20180930120525p:plain
Neo4jは有向グラフなので、ノードの向きも合わせて定義する必要がある。

ノード間の関係を定義するには、

CREATE （identifier）-[:name]->（identifier）

と書く。

プロパティ（属性）

ノードや関係などが持つ実際の値。
ノードの例だと、nameやbornが各プロパティ。関係の場合、rolesがプロパティ。

ラベル

グラフ構造体と呼ばれる複数のノードを１つにまとめるために使われる。オントロジーで言うところの「概念」みたいなもので、世界（ノード）をラベルを付けて、概念化する。
ノードの図で言うと、上のノードには「Person」ラベルがついており、下のラベルには「Movie」ラベルがついていることが分かる。オントロジー上ではPersonやMovieはあくまで概念に過ぎないが、それをインスタンス化してあげると、各エンティティ（ノード）が出来あがるようなイメージ。

ノード・ラベル・プロパティの関係を図示するとこんな感じ。
f:id:kento1109:20180930121915p:plain:w160
（Cypher Query Language(QL)-初級編 #neo4j - クリエーションライン株式会社
より引用）

識別子

ノードや関係性、処理結果などの値を示す構造体。
上の構成要素と異なり、処理中に用いるものに過ぎないが、知っておかないとラベルとの違いに少し混乱する。
こんな感じでノードを表してみる。（emilが識別子）

emil:Person {name:'Emil'}

そうすると、emilがPerson {name:'Emil'}というノード情報を保持した構造体として記憶される。

一般的な使用方法としては、ノードと関係を作成するような例。

CREATE (TheMatrix:Movie {title:'The Matrix', released:1999})
CREATE (Keanu:Person {name:'Keanu Reeves', born:1964})
CREATE (Keanu)-[:ACTED_IN {roles:['Neo']}]->(TheMatrix)

作成されたノードに識別子の情報は一切持たないので、そのノードを使って後の文で処理しないのであれば、名前にこだわる必要は全くない。

どういうことかと言うと、

CREATE (TheMatrix:Movie {title:'The Matrix', released:1999})
CREATE (Keanu:Person {name:'Keanu Reeves', born:1964})

を宣言した後で別のクエリで

CREATE (Keanu)-[:ACTED_IN {roles:['Neo']}]->(TheMatrix)

としても上の識別子のノードではなく、新しいノードが作られてしまう。

最低限こんなこととクエリの書き方さえ覚えればOKかと思う。

2018-09-29

Stan：NMF

はじめに

今回はNMF(Non-negative Matrix Factorization）について勉強する。
日本語では「非負値行列因子分解」と言い、線形次元手法の１つ。
その名前の通り「観測データの非負性」を仮定したモデル。
「非負性」の制約を設けることで、これに適用可能な観測データはより良い表現ができるらしい。
推薦アルゴリズムなどでも使われている。

詳しいことは省略して実装するうえで大事なところを押さえていく。MCMCではなく、最適化により求める方法は「非負値行列因子分解 - NTTコミュニケーション科学基礎研究所」がとても丁寧だった。

尚、今回の記事のいくつかは「ベイズ推論による機械学習入門」を参考とした。

NMF

まず、観測データ（各要素が非負値の行列）行列 $\rm{X}\in \mathbb{R}^{D\times N}$ を以下のように表す。
$\rm{X}\approx\rm{WH}$
尚、 $\rm{W}\in \mathbb{R}^{D\times M}$ ・ $\rm{H}\in \mathbb{M}^{M\times N}$
$M$ は観測データに出現しない潜在変数である。
$\rm{X}$ を要素ごとに書くと次のようになる。
$X_{d,n}\approx\sum_{m=1}^M W_{d,m}H_{m,n}$
例えば、行列 $\rm{X}$ がユーザー×アイテムのレーティング値の場合、 $\rm{W}$ はユーザー×潜在変数（ユーザーの嗜好）を表し、行列 $\rm{H}$ は潜在変数×アイテム（アイテムの特徴）のようなものを表している。
この場合、ユーザー $d$ のアイテム $n$ への評価は、ユーザー $d$ の嗜好ベクトルとアイテム $n$ の特徴ベクトルの内積で決定される。２つのベクトルの方向が近い（ユーザー $d$ の嗜好とアイテム $n$ の特徴が近い）ほど、レーティング値は大きくなるようになる。

モデル

今回は観測データの要素が非負の整数を仮定したモデルを考える。
非負の整数なので、データ生成にはポアソン分布が使える。
各要素は以下のように表せたので、
$X_{d,n}\approx\sum_{m=1}^M W_{d,m}H_{m,n}$
これをそのままポアソン分布にあてはめる。
$p(X_{d,n})\sim{\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n})$

今回はこれが尤度となる。（計算では対数を取る。）
NMFで目指すのは、 $\rm{X}$ を $\rm{WH}$ で近似させること。
各要素 $X_{d,n}$ と $sum_{m=1}^M W_{d,m}H_{m,n}$ が近似できているかどうかは「尤度」で測ることができる。
ポアソン分布の場合、観測値が平均値と一致する値が尤度最大となる値である。この場合、 $X_{d,n}=sum_{m=1}^M W_{d,m}H_{m,n}$ である。観測値が平均から離れるほど尤度は小さくなる
簡単であるが、ポアソン分布での尤度を比べてみた。

dpois(1,1,log=TRUE)
[1] -1
dpois(1,5,log=TRUE)
[1] -3.390562

観測値と平均値が近いほど、尤度が大きいことが分かる。
NMFではポアソン分布の尤度を大きくすることが、 $\rm{X}$ を $\rm{WH}$ させることと同じ意味なのである。

行列全体の尤度は以下のように各要素の積で書ける。
$\begin{eqnarray} L &=& \prod_{d=1}^D\prod_{n=1}^N {\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n}) \\ \log L &=& \sum_{d=1}^D\sum_{n=1}^N \log{\rm Poi} (X_{d,n}|\sum_{m=1}^M W_{d,m}H_{m,n}) \end{eqnarray}$

$W,H$ の各要素は、ガンマ分布を取っておくのが良いと思われる。
$p(W_{d,m})\sim{\rm Gam}(W_{d,m}|a,b)$
$p(H_{m,n})\sim{\rm Gam}(H_{d,m}|a,b)$

コード

コードは以下の通り。
（以前、ネットに落ちていたのを少し改良した。）

data{
  int<lower=1> D; // row dimension of input
  int<lower=1> N; // column dimension of input
  int<lower=1> M; // hidden dimension
  
  int X[D,N]; // original data
  
  real<lower=0> alpha; //hyperparameter
  real<lower=0> beta; //hyperparameter
}

parameters{
  matrix<lower=0>[D, M] W;
  matrix<lower=0>[M, N] H;
}


model{
  //estimate W matrix dist
  for(d in 1:D){
    for(m in 1:M){
      W[d,m] ~ gamma(alpha,beta);
    }
  }
  //estimate H matrix dist
  for(m in 1:M){
    for(n in 1:N){
      H[m,n] ~ gamma(alpha,beta);
    }
  }

  for(d in 1:D){
      for(n in 1:N){
        target += poisson_lpmf(X[d,n] | sum(W[d,:])*sum(H[:,n]));
      }
  }
}

以下のようにしてキックする。

stan.dat <- list(D=nrow(data),D=ncol(data),M=2,X=data,alpha=0.1,beta=0.1)
stan.fit <- stan("NMF.stan", data=stan.dat)

ちなみにデータはこんな感じ。

data
  item1 item2 item3 item4 item5 item6 item7 item8 item9 item10
1     0     2     4     1     0     0     0     1     2      5
2     0     2     3     0     0     0     1     2     2      3
3     1     0     1     4     3     0     1     1     1      0
4     0     0     1     4     5     1     0     1     2      1
5     1     1     5     4     1     1     1     0     1      5

収束はするがデータが小さいので面白い発見は出来ない・・

今回は観測データが非負の整数を取ったが、取る範囲が広いほど推論が難しい気がする。観測データが二値でベルヌーイ分布でやってみるのも面白いかもしれない。

2018-08-13

Pytorch：ライブラリの誤差関数の構造を理解する

PyTorch

はじめに

今まで当たり前のように誤差関数を使っていた。
既に用意されたものであればそれで問題ない。

しかし、誤差関数を自作したいと思った場合、
ライブラリの誤差関数の構造を理解している必要がある。

そんなわけでライブラリの誤差関数について調べたのでメモ。

簡単な復習

簡単に使い方を復習する。

ライブラリの誤差関数を利用する場合、以下のような使い方をする。

import torch
import torch.nn as nn
import torch.nn.functional as F

net = Net() 
outputs = net(inputs)

criterion = nn.MSELoss()

loss = criterion(outputs, targets)
loss.backward()

ネットワーク

今回はシンプルな以下のネットワークを考える。

class Net(nn.Module):

    def __init__(self):
        super(Net, self).__init__()
        self.fc1 = nn.Linear(32, 16)
        self.fc2 = nn.Linear(16, 1)

    def forward(self, x):
        x = F.relu(self.fc1(x))
        x = self.fc2(x)
        return x

パラメータは以下の通り

params = list(net.parameters())
print(params[0].size())
torch.Size([16, 32])

print(params[1].size())
torch.Size([16])

print(params[1])
Parameter containing:
tensor([ 0.0982,  0.0495,  0.1656, -0.1646,  0.1014, -0.0163, -0.0873,
         0.0418, -0.0404,  0.1556, -0.1247,  0.0236, -0.0651,  0.0960,
         0.1342,  0.1203])

入力

以下のように入力データを作成する。

torch.manual_seed(1)

inputs = torch.randn(10, 32)
targets = torch.randn(10)
targets = targets.view(1, -1)
print(targets)
tensor([[ 0.1924,  0.7161, -0.8120, -1.4617,  0.2328,  0.1896, -0.2204,
          0.1491,  0.0100, -0.1243]])

出力

出力結果は以下の通り

outputs = net(inputs)
outputs = outputs.view(1, -1)
print(outputs)
tensor([[-0.2504,  0.0378,  0.1485,  0.1909, -0.1019,  0.2417,  0.1653,
          0.1359, -0.1187,  0.0876]])

誤差計算

まず、誤差関数のインスタンスを生成する。

criterion = nn.MSELoss()

次に出力結果と真の値を誤差関数の入力として、誤差を求める。

outputs = outputs.view(1, -1)
loss = criterion(outputs, targets)
print(loss)
tensor(0.4635)

Backprop

求めた誤差からパラメータの勾配を計算する。

net.zero_grad()

print(net.fc2.weight.grad)
> None
loss.backward()
print(net.fc2.weight.grad)
tensor([[ 0.0611,  0.0976, -0.2338,  0.1976,  0.0020, -0.0788, -0.1145,
         -0.1976,  0.0699, -0.3314,  0.2379, -0.0534,  0.0729, -0.0008,
         -0.2456, -0.2657]])

はじめは勾配がNoneとなっているが、loss.backward()により、勾配ベクトルが計算されていることが分かる。

パラメータ更新（最適化）
今回は確率的勾配降下法（SGD）を用いた。

learning_rate = 0.01
for f in net.parameters():
    f.data.sub_(f.grad.data * learning_rate)

sub_でその値から引数の値を引き算する。

更新後の誤差を確認する。

outputs = net(inputs)
outputs = outputs.view(1, -1)

loss = criterion(outputs, targets)
print(loss)
tensor(0.3956)

誤差が減少していることが確認できた。

実際は多くの最適化関数も用意されているのでそちらを利用する方が良い。

import torch.optim as optim

optimizer = optim.SGD(net.parameters(), lr=0.01)
optimizer.step()

参考：
Neural Networks — PyTorch Tutorials 0.4.1 documentation

基本的なことを復習したところで、もう少し深く理解する。

nn.MSELoss()

ソースを覗いてみる。
torch.nn.modules.loss — PyTorch master documentation

class MSELoss(_Loss):
    def __init__(self, size_average=None, reduce=None, reduction='elementwise_mean'):
        super(MSELoss, self).__init__(size_average, reduce, reduction)

    def forward(self, input, target):
        return F.mse_loss(input, target, reduction=self.reduction)

※長かったのでコメントは省略
実装部分は意外とシンプルでコンストラクタとforward関数のみのクラス。
他の誤差関数も同じようなクラス構造。

F.mse_loss自体が実際に平均二乗誤差を計算している関数。

このクラスが継承している_Lossとは何か。

class _Loss(Module):
    def __init__(self, size_average=None, reduce=None, reduction='elementwise_mean'):
        super(_Loss, self).__init__()
        if size_average is not None or reduce is not None:
            self.reduction = _Reduction.legacy_get_string(size_average, reduce)
        else:
            self.reduction = reduction

このクラス自体はModuleを継承している。

てな感じでみていくと誤差関数のクラスは以下の構造でOKと思われる。

class LossFunction(nn.Module):

def __init__(self):
    super(LossFunction, self).__init__()

def forward(self, inputs, targets):
   loss = function(inputs, targets)
   return loss

ただし、必ずクラス構造である必要はなく、誤差関数だけでも良い。
実際に以下でも正しくパラメータ更新ができた。

loss = F.mse_loss(outputs, targets)