機械学習ゼミ⑤共クラスタリング
2018年7月7日（土）
五百藏夏穂

概要
関係データ学習＊スペクトラルクラスタリング＊確率的ブロックモデル/無限関係モデル＊行列分解

本ゼミでの説明は，“関係データ学習” 石黒勝彦，林浩平（著）に準拠します．

1

関係データ学習

2

関係データとは関係 複数のデータの間に観測，定義される関係にデータ 着目したデータのこと（SNSユーザー間, 購買データ, etc）
1

有向関係
無向関係

起点によって関係が変わるか否か

2
3

単一ドメイン複数ドメインオブジェクトが同じ種類か否か
ex（単一：SNSユーザー関係，複数：購買データ）

グラフ表現
𝑉 = 1,2,30

10

10

1010

行列表現
𝑋 = 𝑥𝑖,𝑗 𝑖, 𝑗 = 1,2,3対称関係データ：単一ドメインかつ無向関係

スペクトラルクラスタリング
非対称関係データ：それ以外

確率的ブロックモデル/無限関係モデル，行列分解

関係データ学習

3

タスク

＊予測

：関係データ学習を通してみ観測のデータの値を推定→アイテム推薦問題＊知識抽出 ：関係データの特性を解析，モデル化し，知見や知識の情報抽出→コミュニティ抽出，クラスタリング

タスク
関係
デー
タの種類

対称二項関係

多項関係

知識抽出

予測

スペクトラル
クラスタリング

行列分解

非対称 確率的ブロックモデル
無限関係モデルテンソル分解スペクトラルクラスタリング
関係データ学習におけるクラスタリング
対象：オブジェクト（ノード）𝑉のクラスタリング
類似性：各オブジェクト𝑖, 𝑗間の“関係”に基づく分類※対象：エッジ（リンク），類似性：オブジェクト自身が持つ特徴量に基づくクラスタリングはここでは本ゼミでは対象外とします．

機械学
習

教師
あり
教師
なし
強化
学
習
etc

コミュニティ検出と密結合グラフ

4・分類（決定木，NN, …etc）

・回帰分析(線形回帰, …etc)

・クラスタリング
・共起分析
・次元圧縮

(k-means, …etc)(協調フィルタリング, …etc)(主成分分析， 特異値分解, …etc)

・Q学習
・TD学習

密結合クラスタグラフ

密結合クラスタ：クラスタ内部で密な結合，クラスタ外部と疎な結合
疎結合クラスタ：クラスタ内部で疎な結合，クラスタ外部と疎な結合

コミュニティ検出では前者
どのクラスタにどのオブジェクトを帰属させるかを計算

疎結合クラスタグラフ

スペクトラルクラスタリング

5グラフカット
疎・密結合の発見や判断に用いる特徴量

1𝑐𝑢𝑡 𝑃 =2𝐾

𝑥𝑖,𝑗𝑘=1 𝑖∈𝑃𝑘 ,𝑗∈𝑃𝑘

グラフカット𝑃𝑘 : 𝑘番目クラスタに所属する頂点集合
𝑃𝑘 : 𝑘番目クラスタに所属しない頂点集合
𝑥𝑖,𝑗 :異なるクラスタに所属する頂点𝑖, 𝑗間の重み頂点間のつながりの弱い部分で分割すると考えてカット最小化問題として密結合クラスタを抽出
スペクトラルクラスタリング
グラフラプラシアンの固有値に着目し，固有値の値によってグラフのノード間の連結性が推定できる固有値分解とk-meansの併用でグラフを構成する頂点のクラスタリングを行う手法

スペクトラルクラスタリング

6①非正規化グラフラプラシアンによるスペクトラルクラスタリング
グラフラプラシアン
入力：対称データ行列𝑋，クラスタ数𝐾

観測データ行列：隣接行列
次数行列：

𝐷 = 𝑑𝑖,𝑗 ,𝑑𝑖,𝑗 =𝑖, 𝑗 = 1,2, ⋯ 𝑁

𝑁
𝑗=1 𝑥𝑖,𝑗 𝑖 = 𝑗0

グラフラプラシアン𝐿の計算

𝑖≠𝑗
𝐿の固有値分解を計算
固有値を昇順に並び替えて最初𝐾個の列固有ベクトルを𝑉 = 𝑣1 , ⋯ 𝑣𝐾 ∈ 𝑅𝑁×𝐾として取り出す

グラフラプラシアン：

𝐿 =𝑋−𝐷
観測データ行列𝑋

次数行列𝐷の計算

次数行列D

グラフラプラシアンL

行列𝑉の各行ベクトルをk-meansでクラスタリング
出力： 𝐾次元空間オブジェクト配置行列V
オブジェクトのクラスタ割り当て𝑍

スペクトラルクラスタリング

7①非正規化グラフラプラシアンによるスペクトラルクラスタリング
固有値分解によるクラスタ抽出
𝐿 = VΛ𝑉 𝑇または

λ𝑘 𝑣𝑘 = 𝐿𝑣𝑘

𝐿は半正定値行列なので
0 ≤ λ1 ≤ λ2 ≤ ⋯ ≤ λ𝐾
この時定数ベクトルの固有値は0なので最小値λ1 =0を取る
定数ベクトルの要素を分割したベクトルを取ることで，𝐾個のクラスタを持つグラフでは最小𝐾個の固有値0（理想状態）※補足：k-means

クラスタ割り当ての計算
列固有ベクトル𝑉 = 𝑣1 , ⋯ 𝑣𝐾 をN個のK次元ベクトルととらえてk-meansを適用
クラスタ割り当て𝑍 = 𝑧𝑖 = 𝑘 が求まる

K個の中心を決める↓
最も近くの中心クラスタに割り当て↓
クラスタ内平均からクラスタ中心を更新↓
中心が変化しなくなったら終了

スペクトラルクラスタリング

8②正規化グラフラプラシアンによるスペクトラルクラスタリング
正規化カット
カット最小化の傾向として各クラスタの大きさに制約をかけるように正規化を行ったカットを用いる

大きなクラスタと，孤立し多数の単頂点クラスタができやすい

レシオカット1
RatioCut =2𝐾𝑖∈𝑃𝑘 ,𝑗∈𝑃𝑘 𝑥𝑖,𝑗𝑘=1

𝑃𝑘=cut(𝑃𝑘 , 𝑃𝑘 )
𝑃𝑘
𝑘=1
𝐾

・各クラスタのサイズ（頂点数）で割り正規化→クラスタが大きいとカットが小さくなる
・非正規化グラフラプラジアンのアルゴリズムで解ける

ノーマライズドカット
1NCut =2𝐾

𝑘=1𝑖∈𝑃𝑘 ,𝑗∈𝑃𝑘 𝑥𝑖,𝑗vol(𝑃𝑘 )=cut(𝑃𝑘 , 𝑃𝑘 )𝑘=1 vol(𝑃𝑘 )𝐾

vol 𝑃𝑘 =・クラスタ内のリンクの重みの総和で割る→クラスタが大きくかつ相互に密に結合しているとカットが小さくなる
・正規化グラフラプラシアンに基づくスペクトラルクラスタリングに帰着

𝑥𝑖,𝑗𝑖,𝑗∈𝑃スペクトラルクラスタリング

9②正規化グラフラプラシアンによるスペクトラルクラスタリング
酔歩正規化ラプラシアン

対称正規化ラプラシアン
1
1−2−2
𝐿𝑠𝑦𝑚 = 𝐷 𝐿𝐷

𝐿𝑟𝑤 = 𝐷 −1 𝐿
比較
＊コミュニティ（密結合）が対象であれば，密結合も考慮したノーマライズドカットの方が望ましい＊𝐿𝑠𝑦𝑚 はk-means前の正規化の分，計算コストが大きいため，𝐿𝑟𝑤 の方が望ましい

ノーマライズドカットをグラフラプラシアンを用いた二次形式に書き換え，それを最小化する
解を制約付きで解くと導出できる

関係データ𝑋の重みに従って酔歩すると想定すると遷移確率行列はT = 𝐷−1 𝑋
ノーマライズドカット最小化＝クラスタ間遷移確率最小化

スペクトラルクラスタリング

10②実データ適用例
Zachary's Karate Club network data
空手クラブのメンバーの交友関係は，レッスン料の違いから2つの集団に分裂

N=34, K=2として各手法を適用
実際の結果と一致しなかったもの
非正規化/酔歩正規化 ： 4/34
対称正規化 ： 7/34
非正規化

対称正規化

酔歩正規化

確率的ブロックモデル

11

確率的ブロックモデルと無限関係モデル
スペクトラルクラスタリングでは，対象関係データの密結合クラスタに対象が限られていた非対称関係データにも適応可能で密結合以外も抽出可能なデータクラスタリング手法＊潜在的なブロック構造を仮定＊確率的アプローチ
関係データに潜むブロック構造を推定し，この時，この構造によりどのような観測データが確率的に得られるのか？

確率的ブロックモデル
Stochastic Block Model（SBM）

無限関係モデルInfinite Relational Model（IRM）

クラスタ数Kが自動的に求まる確率的ブロックモデル
確率的生成モデル

12SBMもIRMもこの一種である系を構成する各変数の値が依存関係に従って順番に生成される過程を確率で表現するモデル
ex)サイコロの目(100回振って出た目を記録)
100回サイコロを振る：変数をサンプリング
1回ごとにサイコロの形状に従って ：パラメータで規定される確率分布
ランダムに目の値が決まる：観測値が生成される

𝑥𝑖 |𝜋~𝑃 𝜋 , 𝑖 = 1,2, … 100

𝑋 = 𝑥1 , ⋯ 𝑥100 :観測値
𝜋 = 𝜋1 , ⋯ 𝜋𝑑 𝜋𝐷 : 面の出やすさパラメータ
𝐷: 面数𝜋𝑑 = 1𝑑

メリット ＊確率分布の形状やパラメータが全て求まれば，実際に値を生成できる＊事前知識を確率分布の形状や種類として反映できる確率モデルのベイズ推定
観測データから確率分布の形状やパラメータ値を推定
本書で説明されるベイズ推定では，ベイズ事後分布を求める

𝑝(𝑍, {𝜃}|𝑋) ∝ 𝑝(𝑋|𝑍, {𝜃}|)𝑝(𝑍, {𝜃}) 周辺化

𝑝(𝑍|𝑋) ∝𝑝 𝑋 𝑍, 𝜃𝑝 𝑍, 𝜃𝑑 𝜃

確率的ブロックモデル

13

確率的ブロックモデルの定式化
𝑥𝑖,𝑗 ∈ {1,0}行（第一ドメイン）インデックス： 𝑖 = 1, ⋯ , 𝑁1列（第二ドメイン）インデックス：

j= 1, ⋯ , 𝑁2

行クラスタインデックス：

𝑘 = 1, ⋯ , 𝐾

列クラスタインデックス：

l= 1, ⋯ , 𝐿

l= 1, ⋯ , 𝐿

𝑋 = 𝑥𝑖,𝑗

j= 1, ⋯ , 𝑁2

観測データ：

𝑋𝑖 = 1, ⋯ , 𝑁1

𝑘 = 1, ⋯ , 𝐾

𝜋1 | 𝛼1 ~Dirichlet(𝛼1 )𝜋2 | 𝛼2 ~Dirichlet(𝛼2 )𝜋：オブジェクトのクラスタへの混合割合パラメータex)K面のサイコロをランダムに生成する分布
𝛼：ハイパラメータ

𝑧1,𝑖 =k|𝜋1 ~Discrete(𝜋1 )𝑧2,𝑗 =l|𝜋2 ~Discrete(𝜋2 )各ドメインのクラスタの混合割合に従ってオブジェクトのクラス
タ割り当てをサンプリングex)K面のサイコロに従って目を出す𝜃𝑘,𝑙 |𝑎0 , 𝑏0 ~Beta(𝑎0 , 𝑏0 )

𝜃:クラスタ間関係強さ（2値関係なのでベータ分布）

𝑥𝑖,𝑗 | 𝜃𝑘,𝑙 , 𝑧1,𝑖 , 𝑧2,𝑗 ~Bernoulli(𝜃𝑧1,𝑖 , 𝑧2,𝑗 )𝑧:割り当てを表す隠れ変数
観測行列の生成確率的ブロックモデル

14

周辺化ギブスサンプラー（CGS）による推論(𝑠)𝑧𝑖 ~𝑝(𝑧𝑖 |𝑋, 𝑍 /(𝑖) ) ∝𝑝 𝑋 (𝑖) 𝑋 /(𝑖) , 𝑧𝑖 , 𝑍 /(𝑖) 𝜃𝑝 𝑧𝑖 , 𝜃 |𝑋 /(𝑖) , 𝑍 /(𝑖) 𝑑 𝜃※確率モデルは隠れ変数𝑧，パラメータ集合 𝜃 ，観測データ𝑋からなるが，結果的に必要なのは隠れ変数の事後分布だけなのでパラメータは周辺化して推論

𝑠：サンプリング回数

𝑋 (𝑖) ：𝑧𝑖 に依存する観測量
𝑋 /(𝑖) ：𝑋 − 𝑋 (𝑖)

事後分布から生成されたとみなされる確率変数のサンプリングを多数行い，サンプルの経験分布をもって事後分布の推定値を計算する
CGSでは，1回のサンプリングである一つの変数だけのサンプリングを順番に行う．
(𝑠)

𝑧𝑖 は𝑠回目の事後分布サンプリング結果(𝑠)𝑧𝑖 ← 𝑧𝑖

: 既知の値として保存𝑧𝑗 𝑖 ≠ 𝑗 のサンプリングに映る確率的ブロックモデル

15

CGSによるSBMの推論式
第一・第二ドメインで対象なので/(𝑖)𝑝(𝑧1,𝑖 = 𝑘|𝑋, 𝑍1 , 𝑍2 ) ∝/(𝑖)𝑝 𝑋 (𝑖) |𝑋 /(𝑖) , 𝑧1𝑖=𝑘 , 𝑍1 , 𝑍2 , 𝜋1 , 𝜋2 , 𝜃𝑘,𝑙/(𝑖)∙ 𝑝 𝑧1𝑖=𝑘 , 𝜋1 , 𝜋2 , 𝜃𝑘,𝑙 |𝑋 /(𝑖) , 𝑍1 , 𝑍2 𝑑𝜋1 𝑑𝜋2 𝑑 𝜃𝑘,𝑙𝑧1𝑖 に依存する関係データ行列の要素に対する尤度

𝑧1𝑖 及びパラメータの𝑧1𝑖 以外の情報が所与の時の事後分布/(𝑖)𝑝(𝑧1,𝑖 = 𝑘|𝑋, 𝑍1 , 𝑍2 ) ∝𝑝 𝑧1,𝑖 = 𝑘 𝜋1 𝑝(𝜋1 |𝑍1𝑖 )𝑑𝜋1 ∙𝜋1 を周辺化した時のクラスタ割り当て𝑧1,𝑖 の確率

𝐿
𝑙=1

𝑁2
Ⅱ 𝑧2,𝑗 =𝑙
𝑗=1[𝑝(𝑥𝑖,𝑗 |𝜃𝑘,𝑙 )]/(𝑖)𝑝(𝜃𝑘,𝑙 |𝑋 /(𝑖) , 𝑍1 , 𝑍2 ) 𝑑 𝜃𝑘,𝑙第二ドメインのクラスタ𝑙ごとにクラスタ間関係強さ𝜃𝑘,𝑙 を周辺化した時の観測データの尤度の影響生成モデルに従って計算
𝑁

𝑁2 𝑥 Ⅱ 𝑧 =𝑙 ) Γ(𝑏 +2Γ(𝑎𝑘,𝑙 + 𝑗=1𝑘,𝑙𝑖,𝑗2,𝑗Γ(𝑎𝑘,𝑙 +𝑏𝑘,𝑙 )/(𝑖)
𝑗=1(1−𝑥𝑖,𝑗 )Ⅱ 𝑧2,𝑗 =𝑙 )𝑝(𝑧1,𝑖 = 𝑘|𝑋, 𝑍1 , 𝑍2 ) ∝ 𝛼1,𝑘×𝑁
Γ(𝑎𝑘,𝑙 )Γ(𝑏𝑘,𝑙 )Γ(𝑎 +𝑏 + 2 Ⅱ 𝑧 =𝑙 )𝑘,𝑙𝑘,𝑙(+)

𝑎𝑘,𝑙 =𝑎0 + 𝑛𝑘,𝑙

𝑗=12,𝑗(−)

𝑏𝑘,𝑙 =𝑏0 + 𝑛𝑘,𝑙確率的ブロックモデル

16

CGSによるSBMの推論式(十分統計量)
定義

𝑁1

𝑚1,𝑘 =Ⅱ 𝑧1,𝑖 = 𝑘(+)

𝑖=1
𝑁1

(−)

𝑖=1
𝑁1

𝑗=1
𝑁2

𝑖=1

𝑗=1𝑛𝑘,𝑙 =𝑛𝑘,𝑙 =𝑁2:クラスタ𝑘に所属するオブジェクト数

𝑥𝑖,𝑗 Ⅱ 𝑧1,𝑖 = 𝑘 Ⅱ 𝑧2,𝑗 = 𝑙(1 − 𝑥𝑖,𝑗 )Ⅱ 𝑧1,𝑖 = 𝑘 Ⅱ 𝑧2,𝑗 = 𝑙𝑧1,𝑖 のクラスタ割り当てを解除してた時の一つ抜き十分統計量𝒛𝟏,𝒊 のサンプリング時に
𝑚1,𝑘 =(+)
𝑛𝑘,𝑙 =(−)
𝑛𝑘,𝑙 =更新

𝑁1𝑖 ′ ≠𝑖,𝑖 ′ =1

Ⅱ 𝑧1,𝑖 ′ = 𝑘 = 𝑚1,𝑘 − Ⅱ 𝑧1,𝑖 = 𝑘

𝑁1

𝑁2𝑖 ′ ≠𝑖,𝑖 ′ =1

𝑗=1

𝑁1

𝑁2𝑖 ′ ≠𝑖,𝑖 ′ =1

𝑗=1

𝑥𝑖 ′ ,𝑗 Ⅱ 𝑧1,𝑖 ′ = 𝑘 Ⅱ 𝑧2,𝑗 = 𝑙(+)

𝑁2(−)

𝑗=1
𝑁2

𝑛𝑘,𝑙 = 𝑛𝑘,𝑙 + Ⅱ 𝑧1,𝑖 = 𝑘(−)(+)= 𝑛𝑘,𝑙 − Ⅱ 𝑧1,𝑖 = 𝑘(1 − 𝑥𝑖 ′ ,𝑗 )Ⅱ 𝑧1,𝑖 ′ = 𝑘 Ⅱ 𝑧2,𝑗 = 𝑙𝑚1,𝑘 = 𝑚1,𝑘 ＋Ⅱ 𝑧1,𝑖 = 𝑘(+):ブロック 𝑘, 𝑙 において𝑥 = 1(0)となる要素数

𝑛𝑘,𝑙 = 𝑛𝑘,𝑙 + Ⅱ 𝑧1,𝑖 = 𝑘𝑥𝑖,𝑗 Ⅱ 𝑧2,𝑗 = 𝑙(1 − 𝑥𝑖,𝑗 ) Ⅱ 𝑧2,𝑗 = 𝑙

𝑗=1𝑁2
𝑗=1(−)= 𝑛𝑘,𝑙 − Ⅱ 𝑧1,𝑖 = 𝑘𝑥𝑖,𝑗 Ⅱ 𝑧2,𝑗 = 𝑙
𝑁2(1 − 𝑥𝑖,𝑗 ) Ⅱ 𝑧2,𝑗 = 𝑙

𝑗=1

対象となる変数の差し引きだ
けで，1回のサンプリングに対応する更新式計算ができる！

無限関係モデル

17

無限関係モデルの定式化
無限個のクラスタの存在を仮定した確率的生成モデル
潜在クラスタ数𝐾, 𝐿を自動的に決定できる
Chinese Restaurant Process(CRP)

Z1 | 𝛼1 ~CRP(𝛼1 )𝑍2 | 𝛼2 ~CRP(𝛼2 )CRPは無次元ディリクレ分布と無次元離散分布を組み合わせたもの

ノンパラメトリックベイズで用いられる確率過程

𝜃𝑘,𝑙 |𝑎0 , 𝑏0 ~Beta(𝑎0 , 𝑏0 )𝑥𝑖,𝑗 | 𝜃𝑘,𝑙 , 𝑧1,𝑖 , 𝑧2,𝑗 ~Bernoulli(𝜃𝑧1,𝑖 , 𝑧2,𝑗 )＊各テーブルに着席できる客数に上限はない＊中華料理店に配置できるテーブル数に上限はない

入店した客をどのテーブルに配置するか決める
𝑛 + 1人目の客のテーブル選択確率

𝑃(𝑧𝑛+1 = 𝑘|𝑧1:𝑛 , 𝛼) ∝𝑛𝑖=1 Ⅱ 𝑧𝑖 = 𝑘𝛼𝑘 ∈ 1,2, ⋯ , 𝐾𝑘 =𝐾+1

無限個のコンポーネントを持つ混合分布を正しく表現
実際のサンプリングではN+1個のテーブルと，有限になるので実装可能

無限関係モデル

18

CGSによるIRMの推論式/(𝑖)𝑝(𝑧1,𝑖 = 𝑘|𝑋, 𝑍1 , 𝑍2 ) ∝/(𝑖)𝑝 𝑋 (𝑖) |𝑋 /(𝑖) , 𝑧1𝑖=𝑘 , 𝑍1 , 𝑍2 , 𝜃𝑘,𝑙/(𝑖)∙ 𝑝 𝑧1𝑖=𝑘 , 𝜃𝑘,𝑙 |𝑋 /(𝑖) , 𝑍1 , 𝑍2 𝑑 𝜃𝑘,𝑙/(𝑖)𝑝(𝑧1,𝑖 = 𝑘|𝑋, 𝑍1 , 𝑍2 ) ∝/(𝑖)𝑝 𝑧1,𝑖 = 𝑘 𝑍1

𝐿
𝑙=1𝑁2Ⅱ 𝑧2,𝑗 =𝑙
𝑗=1[𝑝(𝑥𝑖,𝑗 |𝜃𝑘,𝑙 )]/(𝑖)𝑝(𝜃𝑘,𝑙 |𝑋 /(𝑖) , 𝑍1 , 𝑍2 ) 𝑑 𝜃𝑘,𝑙

場合分けして計算/(𝑖)𝑝(𝑧1,𝑖 = 𝑘 ∈ 1,2, ⋯ , 𝐾 (既存クラスタ)|𝑋, 𝑍1 , 𝑍2 ) ∝𝑁𝑚1,𝑘𝑁2 𝑥 Ⅱ 𝑧 =𝑙 ) Γ(𝑏 +2Γ(𝑎𝑘,𝑙 + 𝑗=1𝑘,𝑙𝑖,𝑗2,𝑗Γ(𝑎𝑘,𝑙 +𝑏𝑘,𝑙 )𝑗=1(1−𝑥𝑖,𝑗 )Ⅱ 𝑧2,𝑗 =𝑙 )𝐿]𝑁2
𝑙=1[Γ(𝑎 )Γ(𝑏 ) ×Γ(𝑎𝑘,𝑙 +𝑏𝑘,𝑙 + 𝑗=1
Ⅱ 𝑧2,𝑗 =𝑙 )𝑘,𝑙𝑘,𝑙/(𝑖)𝑝(𝑧1,𝑖 = 𝑘 + 1 (新規クラスタ)|𝑋, 𝑍1 , 𝑍2 ) ∝𝑁

𝛼1

𝑁2 𝑥 Ⅱ 𝑧 =𝑙 ) Γ(𝑏 +2Γ(𝑎0 + 𝑗=10𝑖,𝑗2,𝑗Γ(𝑎0 +𝑏0 )𝑗=1(1−𝑥𝑖,𝑗 )Ⅱ 𝑧2,𝑗 =𝑙 )𝐿]𝑁2
𝑙=1[Γ(𝑎 )Γ(𝑏 ) ×Γ(𝑎0 +𝑏0 + 𝑗=1
Ⅱ 𝑧2,𝑗 =𝑙 )00

無限関係モデル

19

CGSによるIRMの推論式(十分統計量)
更新

𝑚1,𝑘 = 𝑚1,𝑘 ＋Ⅱ 𝑧1,𝑖 = 𝑘(+)(+)
𝑛𝑘,𝑙 = 𝑛𝑘,𝑙 + Ⅱ 𝑧1,𝑖 = 𝑘(−)(−)

𝑛𝑘,𝑙 = 𝑛𝑘,𝑙 + Ⅱ 𝑧1,𝑖 = 𝑘

𝑁2
𝑗=1𝑁2

もし𝑧1,𝑖 = 𝑘 + 1ならばハットつきの値を全て0にする

𝑥𝑖,𝑗 Ⅱ 𝑧2,𝑗 = 𝑙(1 − 𝑥𝑖,𝑗 ) Ⅱ 𝑧2,𝑗 = 𝑙

𝑗=1

出力方法
サンプリング実現値を用いて事後分布を計算
最初の適当なB期間のサンプルを棄却
𝑆−𝐵
その後M個置きのサンプル実現値だけを抽出（T= 個）
𝑀1
𝑝 𝑧1,𝑖 = 𝑘 𝑋 =𝑇

𝑇
𝑡=1

Ⅱ(𝑡)
𝑧1,𝑖 = 𝑘1
𝑝 𝑧2,𝑗 = 𝑙 𝑋 =𝑇

具体的なクラス割り当て
事後分布最大化(MAP)に基づく

𝑐1,𝑖 = argmax 𝑝 𝑧1,𝑖 = 𝑘 𝑋
𝑐2,𝑗 = argmax 𝑝 𝑧2,𝑗 = 𝑙 𝑋𝑇
𝑡=1(𝑡)

Ⅱ 𝑧2,𝑗 = 𝑙無限関係モデル

20②実データ適用例＊Zachary's Karate Club network data
4つに分裂
現実と分裂結果の一致度はスペクトラルクラスタリングより良い値に．
人間の思い込みに反するがよりよく関係データを表現＊Enron社内のEメールデータ
メール送信者（K），受信者（L）の二値データ

K=6, L=3クラスタに分割
経営層コミュニティ
関係データ行列真のクラスタ構造がわ
からないデータでもクラスタリングを実行可能

IRM適用結果

繰り返し回数ごとのクラスタ数変化

行列分解

21
𝐽

単純行列分解例として映画の推薦を考えます

𝐼: 顧客 𝐽:映画
要素の値が評価値となる行列

𝐼𝑋

𝑅
≅ 𝐼 𝑈

𝐽𝑉𝑇

𝑅𝑥𝑖𝑗 ≅ 𝑢𝑖,𝑎𝑐𝑡𝑖𝑜𝑛 𝑣𝑗,𝑎𝑐𝑡𝑖𝑜𝑛 + 𝑢𝑖,ℎ𝑜𝑟𝑟𝑜𝑟 𝑣𝑗,ℎ𝑜𝑟𝑟𝑜𝑟
より一般に𝑅 ∈ 𝑁個の評価軸があり，すべての顧客，映画について評価値が得られているとする
𝑋のランクR行列分解

𝑋 ≅ 𝑈 𝑉𝑇

因子行列目的関数

𝐸 = 𝑋 − 𝑈 𝑉𝑇𝑢𝑖,: ∈ R :顧客𝑖の嗜好をまとめたもの𝑣𝑗,: ∈ R :映画𝑗の成分をまとめたもの
𝑈 ∈ R𝐼×𝑅 :行方向にまとめた
𝑉 ∈ R 𝐽×𝑅 :列方向にまとめた

𝑋と𝑈 𝑉 𝑇 の近似誤差を最小化する問題ととらえる
2乗誤差最小化問題
1
1min 𝐸 2Fro = min 𝑋 − 𝑈 𝑉 𝑇 2Fro𝑈,𝑉 2𝑈,𝑉 2

行列分解
𝑙 2 正則行列分解22

𝑋に大きいノイズが加わっていた時などは過学習の危険がある→正則化（制約を付加し，解の範囲を狭める）

1
λ𝑇 2
min 𝑋 − 𝑈 𝑉 Fro +𝑈,𝑉 22

𝑈 2Fro + 𝑉 2Fro
λ ≥ 0は正則化強さを決定する係数

非負行列分解テキストの単語数データなど，非負値しかとらないもの

min 𝐽×𝑅

𝐼×𝑅 ,𝑉∈𝑅𝑈∈𝑅++Ex)1
𝑋 − 𝑈 𝑉 𝑇 2Fro
2

同時確率
𝑥𝑖,𝑗 = 𝑝(文書𝑖, 単語𝑗)

＊単語の確率分布は典型的なパターンを持つ＊文書における単語の確率分布はパターンの組み合わせ

行列分解

23

アルゴリズム
目的関数の導関数を0としても，連立方程式を代数的に解くことができない

勾配法による最適化

一次交互勾配法
𝐾 ∈ 𝑁個の変数𝜃1 , ⋯ 𝜃𝐾 を引数にとる関数ℎ(𝜃1 , ⋯ 𝜃𝐾 ) ∈Rでは𝑛𝑒𝑤𝜃𝑘𝑛𝑒𝑤 = 𝜃𝑘 − η𝛻𝜃𝑘 ℎ(𝜃1𝑛𝑒𝑤 , ⋯ 𝜃𝑘−1, 𝜃𝑘 , 𝜃𝑘+1 , ⋯ , 𝜃𝐾 )𝑙2 正則行列分解の交互最適化
目的関数と連鎖率より

𝛻U 𝑓𝑙2 = −𝑋𝑉 + 𝑈 𝑉 𝑇 𝑉 + λ𝐼
𝛻𝑉 𝑓𝑙2 = −𝑋 𝑇 𝑈 + 𝑉 𝑈 𝑇 𝑈 + λ𝐼
𝑋が密でデータが大きい場合計算量の負荷が大きい

η:学習率

入力：行列𝑋，ランクR，正則化係数λ，学習率η，精度ε
𝑈, 𝑉を乱数で初期化繰り返し𝑔 ← 𝑓𝑙2 𝑈, 𝑉𝑈 ← 𝑈-𝛻𝑈 𝑓𝑙2 (U, V)
𝑉 ← 𝑉-𝛻𝑉 𝑓𝑙2 (U, V)
until 𝑔 − 𝑓𝑙2 (U, V) /𝑓𝑙2 (U, V)< ε
出力：因子行列𝑈, 𝑉

行列分解

24

疑似二次交互勾配法
一次交互勾配法は実装は比較的容易であるが，収束が遅いことも．
二階微分の情報を用いて収束を早くする𝜃𝑘𝑛𝑒𝑤 = 𝜃 − η (𝛻𝛻 T 𝑔(𝜃))−1 𝛻𝑔(𝜃)ヘッセ行列

一次勾配法に比べて収束速度は速いが，ヘッセ行列の逆行列を計算しなければならず，計算量は高次元

入力：行列𝑋，ランクR，正則化係数λ，学習率η，精度ε
𝑈, 𝑉を乱数で初期化

𝑙2 正則行列分解の交互最適化

𝑈 𝑛𝑒𝑤 = 1 − η 𝑈 + ηXV(𝑉 𝑇 𝑉 + λ𝐼)−1
𝑉 𝑛𝑒𝑤 = 1 − η 𝑉 + ηXU(𝑈 𝑇 U + λ𝐼)−1繰り返し𝑔 ← 𝑓𝑙2 𝑈, 𝑉𝑈 ← 1 − η 𝑈 + ηXV(𝑉 𝑇 𝑉 + λ𝐼)−1𝑉 ← 1 − η 𝑉 + ηXU(𝑈 𝑇 U + λ𝐼)−1until 𝑔 − 𝑓𝑙2 (U, V) /𝑓𝑙2 (U, V)< ε
出力：因子行列𝑈, 𝑉

行列分解

25

交互勾配法＋射影勾配降下法（制約条件）
射影勾配降下法

非負制約などの制約を取り入れたい場合

u: 制約を満たす集合

𝑈 𝑛𝑒𝑤 =Proj𝑢 [𝑈 − η𝛻𝑈 𝑓𝑙2 𝑈, 𝑉 ]Proj𝑢 [U]=𝑎𝑟𝑔𝑚𝑖𝑛𝑍∈𝑢 𝑈 − 𝑍 2Fro
通常の勾配法の後に，空間𝑢からはみ出た部分を戻す
非負条件では，入力：行列𝑋，ランクR，正則化係数λ，学習率η，精度ε
𝑈, 𝑉を乱数で初期化

𝑈 𝑛𝑒𝑤 =[𝑈 − η𝛻𝑈 𝑓𝑙2 𝑈, 𝑉 ]+[𝑈]+ ：𝑢𝑖,𝑗 > 𝑜の時𝑢𝑖,𝑗 ，𝑢𝑖,𝑗 ≦ 𝑜の時0繰り返し𝑔 ← 𝑓𝑙2 𝑈, 𝑉𝑈 ← [ 1 − η 𝑈 + ηXV(𝑉 𝑇 𝑉 + λ𝐼)−1 ]+𝑉 ← [ 1 − η 𝑉 + ηXU(𝑈𝑇 U + λ𝐼)−1 ]+until 𝑔 − 𝑓𝑙2 (U, V) /𝑓𝑙2 (U, V)< ε
出力：因子行列𝑈, 𝑉行列分解26
入力：行列𝑋，ランクR，正則化係数λ，学習率η，反復回数T

確率勾配降下法
𝑋が密で𝐼, 𝐽がお真理に大きいとメモリに乗りきらないことがある
少数のサンプルで勾配を近似的に計算し最適化を行う方法．
目的関数𝑓がサンプル𝑛のみに依存する関数𝑓𝑛の和で書けるとき，代わりに𝑓𝑛 の勾配を使う

要素ごと[𝐼]上の一様乱数から𝑖をサンプリング[𝐽]上の一様乱数からjをサンプリング(𝑖,𝑗)𝑢𝑖: ← 𝑢𝑖: + η(t) 𝛻𝑢𝑖: 𝑓𝑙2 (𝑢𝑖: , 𝑣𝑗: )λ= −𝑥𝑖,𝑗 𝑣𝑗: + 𝑢𝑖: (𝑣𝑗:𝑇 𝑣𝑗: + )𝐽λ(𝑖,𝑗)𝑇𝛻𝑣𝑗: 𝑓𝑙2 = −𝑥𝑖,𝑗 𝑢𝑖: + 𝑣𝑗: (𝑢𝑖: 𝑢𝑖: + )𝐼(𝑖,𝑗)

𝛻𝑢𝑖: 𝑓𝑙2(𝑖,𝑗)𝑣𝑗: ← 𝑣𝑗: + η(t) 𝛻𝑣𝑗: 𝑓𝑙2

出力：因子行列𝑈, 𝑉共起行列のような𝑋が複数の行列の和で与えられており要素へのアクセスにもコストがかかる場合（ex. 単語の共起行列）
単純行列分解での目的関数を
1
𝑓𝑠𝑡𝑑 𝑋; 𝑈, 𝑉 = 𝑋 − 𝑈 𝑉 𝑇 2Fro と表記すると
22+ 𝑀2 𝜎 2
1𝜎 2 =𝑃𝑚 𝐶 (𝑚) − 𝑃𝑚 [𝐶 (𝑚) ] Fro , 𝑃𝑚 [∙]=𝑀(𝑢𝑖: , 𝑣𝑗: )𝑒𝑛𝑑 𝑓𝑜𝑟

行列ごと

𝑓𝑠𝑡𝑑 𝑋; 𝑈, 𝑉 = 𝑃𝑚 𝑓𝑠𝑡𝑑 M𝐶 𝑚 ; 𝑈, 𝑉

𝑈, 𝑉を乱数で初期化𝑓𝑜𝑟 𝑡 = 1, ⋯ , 𝑇𝑀
𝑚=1(∙)

入力：行列{𝐶1 ,⋯, 𝐶𝑀 , }，ランクR，正則
化係数λ，学習率η，反復回数T

𝑈, 𝑉を乱数で初期化
𝑓𝑜𝑟 𝑡 = 1, ⋯ , 𝑇𝑚 ← 𝑡 mod 𝑀
𝑈 ← 𝑈 + η(t) 𝛻𝑈 𝑓𝑙2 (M𝐶 (𝑚) ; 𝑈, 𝑉)
𝑉 ← 𝑉 + η(t) 𝛻𝑈 𝑓𝑙2 (M𝐶 (𝑚) ; 𝑈, 𝑉)
𝑒𝑛𝑑 𝑓𝑜𝑟

出力：因子行列𝑈, 𝑉

行列分解

27

欠損値がある場合の行列分解
欠損値を除外欠損部分を完全に未知だとし，学習の際に一切使わない方法
𝑋の観測部分のインデックス集合Ω

目的関数：
勾配：

(𝑥𝑖,𝑗 − 𝑢𝑖:𝑇 𝑣𝑗: )2(𝑖,𝑗)∈Ω𝛻U 𝑓𝑙2 = − 𝑋 − 𝑈𝑉 𝑇 ∗ 𝑀 𝑉 + λ𝑈

マスク𝑚𝑖,𝑗 =

欠損値を補完𝑋の欠損値を何らかの形で補完した行列をXとする

目的関数：

2

X − 𝑈 𝑉 𝑇 Fro

補完の仕方はいくつかあるが，交互最適化と組み合わせる場合，𝑈, 𝑉をランダムに初期化した後t+1回目の反復においてt回目の予測値で補完
(𝑡)

(𝑡)

𝑈 𝑡+1 の更新： i, j ∉ Ωにおいて𝑥𝑖,𝑗 = (𝑢𝑖 )𝑇 𝑣𝑗(𝑡+1) 𝑇 (𝑡)) 𝑣𝑗

𝑉 𝑡+1 の更新： i, j ∉ Ωにおいて𝑥𝑖,𝑗 = (𝑢𝑖1 (𝑖, 𝑗) ∈ Ω0 それ以外

参考文献＊「関係データ学習」 石黒勝彦・林浩平 講談社＊「続・わかりやすい パターン認識」 石井健一郎・上田修功 オーム社

28