確率モデルを用いたテンソル因子化法の拡張に関するサーベイ

確率モデルを用いたテンソル因子化法の拡張に関するサーベイ

Survey of Probabilitic Tensor Factorization Methods

林 浩平

1∗

_{池田 和司}

1

Kohei Hayashi

1

_{Kazushi Ikeda}

1

1

_{奈良先端科学技術大学院大学}

1

_{Nara Institute of Science and Technology}

Abstract: This survey focuses on tensor factorization methods based on stochastic models. They are natural extensions of the conventional methods to Bayesian statistics and are easier to introduce prior knowledge. Some of their applications are also given.

1

はじめに

テンソル因子化法は多次元配列の構造を持つデータ の特徴抽出，次元縮約，あるいはデータ補完に使用さ れる手法であり，fMRI 解析，顔認識，社会ネットワー ク解析，蛍光スペクトル分析など幅広い分野に応用さ れている [1]．またテンソル因子化法は高次元データの 可視化やクラスタリングなどデータマイニング分野に おいても有用な手法である． 近年テンソル因子化法は確率モデルとしても解釈さ れ，いくつかの拡張が行われている．本サーベイはこ こ二年間の内に発表された 4 つの確率モデルに焦点を 当て，それぞれの特徴について述べる．本稿での構成 を以下に示す．2 章では本稿で用いる記法や代表的な テンソル因子化法である Tucker と PARAFAC につい て説明する．3 章では Tucker や PARAFAC を拡張し た 4 つの確率モデルについて記述する．5 章にて本稿 をまとめる．

2

テンソル因子化法

主成分分析 (Principal Component Analysis, PCA) や特異値分解 (Singular Value Decomposition, SVD) を初めとする行列因子化法は，グラフ構造など行列形 式で表現されるデータの特徴抽出や次元縮約を行う手 法として良く知られている．テンソル因子化法は行列 因子化法の一般化であり，テンソルデータをより低次 元なパラメータで近似する．本章では代表的なテンソ ル因子化法である Tucker と PARAFAC について簡単 な紹介を行う．ただし，本稿はテンソル因子化法の確 率モデルに関するサーベイであり，通常のテンソル因 ∗_{連絡先：奈良先端科学技術大学院大学情報科学研究科}       〒 630-0192 奈良県生駒市高山町 8916 番地の 5        E-mail: [email protected] 子化法については簡単に紹介するにとどめた．テンソ ル因子化法全般に関しては [9] や [1] が詳しい．また Faloutsos らによる ICML2007 のチュートリアル1 _は 図が豊富でありテンソル因子化法の直観的な理解に役 立つ．

2.1

準備

M 次元配列を本稿では M 次のテンソルと呼ぶ．テ ンソルの次元はモードという別名を持つ．特別な場合 として，1 次，2 次のテンソルはそれぞれベクトル，行 列となる．本稿では簡単のため 3 次のテンソルのみを 扱う．それ以上の高次テンソルに関しても 3 次の場合 と同様の議論が成り立つ場合がほとんどである．また 本稿では [1] の表記を踏襲し，テンソルを大文字+ボー ルド体+下線 (X) と表記する．同様に行列は大文字+ ボールド体 (X)，ベクトルは小文字+ボールド体 (x) と する．また xijkは 3 次のテンソル X の第 1，2，3 モー ドの添字がそれぞれ i, j, k 番目となる要素を指す．

2.2

Tucker

Tucker は，縮約された 1 つのテンソルとそれぞれの モードに対応する複数の行列によるテンソルの分解法 をいくつか提案した [15]．中でも Tucker3 として知ら れる手法は，例えば 3 次のテンソル X∈ RI_×J×Kを次 式の形で分解する． xijk= Q ∑ q=1 R ∑ r=1 S ∑ s=1 tqrsuiqvjrzks+ εijk (1) 1 http://www.cs.cmu.edu/˜christos/TALKS/ICML-07-tutorial/ 人工知能学会研究会資料 SIG-DMSM-A902-06 (10/18)

ここで，3 次のテンソル T∈ RQ×S×R_{はコアテンソル} と呼ばれる．U ∈ RI×Q_{, V∈ R}J×R_{, Z ∈ R}K×S_は因 子行列と呼ばれ，それぞれのモードの相関情報が抽出 される．U, V, Z はコアテンソル T を通して相互作用 し X を復元する．また E = [εijk] は残差項である．一 般に Q ¿ I, R ¿ J, S ¿ K であり，T は X の低次 元表現として捉えることができる．Tucker は後述する PARAFAC に比べパラメータの自由度が高く，対象の テンソルをよりコンパクトに表現できる．しかしなが ら非常に高い自由度を持つゆえに分解の一意性は保証 されない． Tucker3 はテンソル X が唯一得られた観測だとして 分解する手法であった．しかしながら，応用としてはテ ンソルが複数個得られるような状況も考えられる．こ のようなデータを扱う 1 つの方法としては，これら複 数のテンソル X1, . . . , XN を 1 つにまとめた，Xnよ り 1 次元高次なテンソル ˜X = [X₁, . . . , X_N] に対して Tucker 分解を行なうという方法が考えられる．しかし ながら，こうした場合サンプルに対応するモードの独 立性は陽に仮定されない．すなわちサンプル間に何ら かの相関性を仮定しパラメータ推定を行うことになる． 一般にサンプルは独立同分布 (i.i.d) で得られると仮定 することが多いが，その場合この仮定は不適切となる． サンプルの独立性 (あるいはあるモードに対する独立 性) を考慮した手法としては Tensor PCA [7] が提案さ れている．これは因子行列は全サンプルで共通とし，コ アテンソル_{{Tn} がサンプル {Xn} ごとに用意された} 分解となっている (式 2)． xijkn= Q ∑ q=1 R ∑ r=1 S ∑ s=1 tqrsnuiqvjrzks+ εijkn. (2) Tucker の不定性を減らすため制約を入れた手法は， コアテンソルや因子行列に対して直交性の制約を入れ た高次 SVD [8] をはじめとしていくつか提案されてい る．いくつか提案されている．詳しくは [16] を参照さ れたい．

2.3

PARAFAC

PARAFAC は対象データをランク 1 のテンソルの足 し合せで表現する．ランク 1 のテンソルとは，xijk = uivjzkのように対象テンソルの次数分のベクトルの積で 表現可能なテンソルのことである．ここで X∈ RI×J×K に対する PAFARAC 分解は以下とかける． xijk = Q ∑ q=1 uiqvjqzkq+ εijk. (3) 行列 U∈ RI×Q_{, V∈ R}J×Q_{, Z∈ R}K×Q_{は Tucker3 と} 同じくそれぞれ X の 1，2，3 番目のモードに関する低 次元表現と解釈できる．また，X が 3 次のテンソルの時 に限り，式 (3) を次式のように行列形式で書き直せる． Xk = UDkVT+ Ek, k = 1, . . . , K. (4) ただし，Xk ={xijk|i = 1, . . . , I, j = 1, . . . , J}, D = diag(zk)，zkは Z の k 番目の行ベクトルである．式 (1) から分かる通り，PARAFAC は Tucker の特殊ケースと みなせる．すなわち Tucker 分解において，Q = R = S とし T の成分を i = j = k のとき tijk = 1，それ以外 tijk= 0 とすると PARAFAC と一致する． 式 (3) にはスケール変換による不定性が残っている が，U, V, Z を正規化した行列 ˜U, ˜V, ˜Z に置き換える ことでこの問題を解消することができる．すなわち， xijk= Q ∑ q=1 λqu˜iq˜vjqz˜kq+ εijk. (5) ここで U, V, Z の各行ベクトルのノルムは非負ベクト ル λ∈ R+Qに吸収されることに注意されたい． 2 次のテンソル (すなわち行列)X∈ RI×J_{の PARAFAC} 分解は以下と書ける． X = UDVT+ E. (6) ただし，D = diag(λ) とする．式 (7) に加えて U, V に 正規直交の条件を入れると SVD と等価となる．Tucker とは異なり，3 次以上のテンソルに対して PARAFAC 分解が一意に決まるための十分条件は良く知られてい る．詳しくは [9] を参照されたい．

2.4

パラメータの推定アルゴリズム

通常，テンソル因子化法は因子行列の次元 Q, R, S を 固定した上で二乗誤差最小化によりパラメータを推定 する．例えば Tucker の一般的な目的関数は以下となる． argmin T,U,V,Z I ∑ i=1 J ∑ j=1 K ∑ k=1 ||xijk− yijk||2 (7) where yijk= Q ∑ q=1 R ∑ r=1 S ∑ s=1 tqrsuiqvjrzks. 式 (7) を見れば分かる通り，テンソル因子化法の目的関 数はパラメータ同士が複雑に絡み合っているため，基本 的には解は閉形式で求まらない．代わりに他のパラメー タを固定した上である 1 つのパラメータを最小二乗誤 差に基づいて更新し，これを他のパラメータに対して も繰り返し計算する Alternating Least Square (ALS) が良く用いられる．ALS は実装が容易であるためテン

ソル因子化法を解くアルゴリズムとして一般的に良く 用いられるが，解が収束するまでに通常多くの反復回 数を必要とし，また最適解が得られる保証はない． 式 (7) は誤差項 εijkを平均 0 分散 1 のガウスノイズ とした時の対数尤度関数と等価となる．よって式 (7) の最小化によって得られたパラメータ T, U, V, Z は最 尤推定値と解釈できる．

3

確率モデルへの拡張

従来のテンソル因子化法を，特にベイズ統計の枠組 みに基づいた確率モデルへ拡張することは，以下のよ うな利点がある [3]： 拡張性 有用な事前分布を容易に導入できる． 欠損値の推定 観測に欠損値が含まれていても，適切な 観測モデルと事前分布の導入によりそれらを推定 できる． ベイズ推論 パラメータの事後分布を推定することで， よりロバストな解が得られることが期待できる． モデル選択 因子行列の次元 (X における Q, R, S) を決 定するための 1 つの基準であるエビデンスを計算 できる． 本章ではこれらの利点を活かしたテンソル因子化法の 確率モデルを 4 つ紹介し，モデルの特徴，計算アルゴ リズム，応用についてそれぞれ記述する．表 1 に本稿 で紹介するモデルを簡単にまとめた．

3.1

pTucker

Tucker3 は非常に自由度の高い手法であるが，パラ メータの不定性が高く過学習の危険性があることを前 章で述べた．過学習を防ぐ方法の 1 つとしてパラメー タの正則化が挙げられるが，Chu らは式 (1) における コアテンソル T を隠れ変数とし，T に平均 0 分散 1 の ガウス事前分布を与えた pTucker を提案した [4]． この研究の最大の特徴は T そのものは推定せず，周 辺化した点にある．このモデルにおいてコアテンソル に関する周辺化計算は厳密に解け，周辺尤度はガウス 分布で得られる．この周辺化によってパラメータ数が 多く不定性の高い T を陽に推定する必要がないため， 過学習をうまく抑えることができかつ計算コストの削 減にも役立っている．これは確率モデルならではの発 想であり，確率モデルが持つ利点をうまく活かした拡 張といえる．また周辺尤度は Probabilistic PCA [14] 図 1: pTucker のグラフィカルモデル． と非常に良く似た形となる．pTucker のグラフィカル モデルを図 1 に示す2_． U, V, Z は，平均 0 分散 1 のガウス事前分布のもと で MAP 推定 (=二乗正則化項つき最尤推定) により得 られる． [4] では人工データを用いて通常の Tucker と pTucker の比較実験を行っている．その結果，コアテンソルの 次元 Q, R, S を真の値以上に設定しても pTucker は常 に真の解に近い推定をしていたのに対し，Tucker は汎 化誤差が増大した．また映画レーティングや様々な実 データの欠損予測も行い，いずれも pTucker が優れた 汎化性能を示すことを確認した．

3.2

Probabilistic 2D PCA (p2DPCA)

Tensor PCA を拡張したモデルとして p2DPCA [18] が提案されている．p2DPCA は pTucker と同じくコア テンソル_{{Tn|n = 1, . . . , N} や因子行列 U, V, Z に平} 均 0 分散 1 のガウス事前分布を仮定することで Tucker が持つ不定性を減らした手法となっている． またパラメータの推定に EM アルゴリズムを用いて いる点も pTucker と同じであるが，p2DPCA において は E ステップが厳密に解けないため，事後分布の近似 法の一種である変分法を用いている． 実験では応用として心臓の超音波画像の二値分類を 行っており，実験結果より p2DPCA の優位性を確認し ている．

3.3

Bayesian Tensor Analysis (BTA)

p2DPCA をさらに拡張した確率モデルとして，BTA [12] が提案されている．これは p2DPCA と同じく{Tn} に 平均 0 分散 1 のガウス事前分布を仮定しているが，加え て U, V, Z にそれぞれの行ベクトルごとに独立な精度 (逆分散) パラメータ αU, αV, αZ を持つガウス事前分 布を与えている．この事前分布を用いる枠組みは Au-tomatic Relevance Determination (ARD) [10] と呼ば

2_{グラフィカルモデルとは確率変数の依存関係を記述するダイア}

表 1: 本稿で紹介する手法まとめ．

モデル名称 因子化方法 特徴 推定アルゴリズム 応用

pTucker [4] Tucker コアテンソルを周辺化 MAP, EM レーティング予測 p2DPCA [18] Tensor PCA フルベイズ 変分 EM 画像識別 BTA [12] Tensor PCA ARD によるスパース制約 変分 EM 3 次元顔モデリング DETF [6] PARAFAC 状態空間モデル EP EM 時変ネットワーク予測 れ，行列の余分な次元を削減する効果があることで知 られている．行列因子化法のモデル選択に ARD を用 いた例は Bayesian PCA [2] を初めとして良く知られて いるが，テンソル因子化法に ARD を適用した手法は 我々の知る限り BTA がはじめてである．クロスバリ デーションや赤池情報量基準 (AIC) といった情報量基 準によるモデル選択は，1 つ 1 つのモデルについてそ れぞれ個別に計算を行う必要があるため膨大な計算量 となるが，ARD の枠組みを用いることでモデル選択を 自動かつ効率的に行うことが出来る．BTA のグラフィ カルモデルを図 2 に示す． 図 2: BTA のグラフィカルモデル． [12] では BTA の応用として 3 次元顔データのモデリ ングを行っている．具体的には，モードがそれぞれ個 人，表情，特徴点からなる 3 次のテンソルデータに対 し，BTA による特徴抽出を行っている．著者は実験結 果より ARD によって通常は適切な次元が選択される と結論付けている．

3.4

Dynamic Exponential family

Ten-sor Factorization (DETF)

通常のテンソル因子化法は，各モードにおける添字 の順番に意味はない．すなわち X の添字 i, j, k につい て，2 つの添字を固定した上での残りの添字の並び換 えは対応する因子行列の行ベクトルの並び換えで吸収 される．しかしながら，実世界には時系列データのよ うに添字の順番に意味のあるデータも多く存在する． DETF [6] はあるモードの添字の順番に対して連続性を 仮定した，PARAFAC ベースの確率モデルである．具 体的には，時刻 k 毎に行列 Xk が観測されるような状 況を考え，状態 ztの線形変換 UDkVTが Xkの期待 値パラメータとなる状態空間モデルを仮定する．これ は式 (4) からのアナロジーとして捉えることができる． この時系列への拡張に加えて，DETF はその観測モデ ルを通常のテンソル因子化法が仮定するガウス分布か ら，より広い枠組みである指数型分布族への一般化を 行っている．これにより，2 値や頻度で表現される関係 データ (例えば E メール送受信関係) などをより柔軟 かつ統一的に扱うことができる．DETF のグラフィカ ルモデルを図 3 に示す． 図 3: DETF のグラフィカルモデル． DETF はパラメータ推定のため EM アルゴリズムを 導出しているが，E ステップの近似として Expectation Propagation (EP) 法 [11] を用いている． [6] では DETF の応用としてメール送受信データの 欠損値の予測を行っている．これは，時刻 k における メールネットワークの隣接行列 Xk(xijkは時刻 k にお けるユーザ i から j へのメール回数) を k = 1, . . . , K ま で観測した時，_{{Xk} に含まれるいくつかの欠損値を推} 定する問題である．実験結果より，通常の PARAFAC よりも時系列構造を導入した DETF の方が欠損予測に おいて優れていることが確認された．

4

終わりに

本稿では，Tucker や PARAFAC に代表されるテンソ ル因子化法の確率モデルに関する研究を 4 つ紹介した が，この他にもいくつかの手法が提案されている．BTA の著者 Tao はテンソル因子化法に対して情報量基準に よるモデル選択を適用し，その効果を実データにて検

証している [13]．また [5] は厳密には確率モデルでは ないが，非負テンソル因子化のコスト関数を二乗誤差， KL ダイバージェンス，板倉・斉藤距離など様々な距離 尺度を含む α, β ダイバージェンスに一般化した研究で ある． 今回紹介した手法にはベイズ予測分布を計算したも のは無かった．予測分布が厳密に解ける場合は限られる ため，実際にテンソル因子化法の予測分布を導出する にはなんらかの近似法が必要となるだろう．さらにテ ンソル因子化法は膨大なパラメータを持つため，周辺 化の計算量を抑える工夫も兼ね備えたアルゴリズムの 考案が望まれる．また通常テンソル因子化法では観測 ノイズは等方ガウス分布を仮定することが多いが，実 際はノイズの分散にもなんらかの構造が含まれている 場合もある．行列因子化法においては，ノイズを行列 変量 (matrix-variate) なガウス確率変数に拡張したモ デルが近年活発に提案されている [17] [19]．テンソル 因子化法においては同種の拡張は我々の知る限りまだ 行われていないが，行列因子化法と同様に予測精度の 向上が期待できるだろう．

参考文献

[1] E. Acar and B. Yener. Unsupervised multiway data analysis: A literature survey. IEEE Transactions on

Knowledge and Data Engineering, 21(1):6–20, 2009.

[2] C. M. Bishop. Variational principal components. In

Artificial Neural Networks, 1999. ICANN 99. Ninth International Conference on (Conf. Publ. No. 470),

Vol. 1, 1999.

[3] C. M. Bishop. Pattern Recognition and Ma-chine Learning (Information Science and Statistics).

Springer, October 2007.

[4] W. Chu and Z. Ghahramani. Probabilistic models for incomplete multi-dimensional arrays. In Proceedings

of the 12th International Conference on Artificial In-telligence and Statistics, 2009.

[5] A. Cichocki, R. Zdunek, S. Choi, R. Plemmons, and S. Amari. Non-negative tensor factorization using al-pha and beta divergences. In Proceedings of IEEE

In-ternational Conference on Acoustics, Speech and Sig-nal Processing, 2007., Vol. 3, pp. III–1393–III–1396,

June 2007.

[6] 林, 平山, 石井. 指数族行列因子化の状態空間モデル

への拡張と時系列関係データ解析への応用. 信学技報,

NC2008-175, 2009.

[7] L. De Lathauwer. Signal Processing Based on

Mul-tilinear Algebra. PhD thesis, Katholike Universiteit

Leuven, 1997.

[8] L. De Lathauwer, B. De Moor, and J. Vande-walle. A multilinear singular value decomposition.

SIAM Journal on Matrix Analysis and Applications,

21(4):1253–1278, 2000.

[9] T. G. Kolda and B. W. Bader. Tensor decompositions and applications. Technical report, Sandia National Laboratories, Albuquerque, NM and Livermore, CA, December 2007.

[10] D. J. C. MacKay. Bayesian interpolation. Neural Computation, 4(3):415–447, May 1992.

[11] T. P. Minka. Expectation propagation for approx-imate Bayesian inference. Uncertainty in Artificial

Intelligence, 17:362–369, 2001.

[12] D. Tao, M. Song, X. Li, J. Shen, J. Sun, X. Wu, C. Faloutsos, and S. J. Maybank. Bayesian Tensor Approach for 3-D Face Modeling. IEEE Transac-tions on Circuits and Systems for Video Technology,

18(10):1397–1410, 2008.

[13] D. Tao, J. Sun, X. Wu, X. Li, J. Shen, S. Maybank, and C. Faloutsos. Probabilistic tensor analysis with Akaike and Bayesian information criteria. Neural

In-formation Processing, pp. 791–801, 2008.

[14] M. Tipping and C. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical

Society. Series B, statistical methodology, 61(3):611–

622, 1999.

[15] L. Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279–311, September 1966.

[16] J. Ye. Generalized low rank approximations of matri-ces. Machine Learning, 61(1-3):167–191, November 2005.

[17] K. Yu, W. Chu, S. Yu, V. Tresp, and Z. Xu. Stochas-tic relational models for discriminative link predic-tion. In B. Sch¨olkopf, J. Platt, and T. Hoffman eds.,

Advances in Neural Information Processing Systems 19, pp. 1553–1560. MIT Press, Cambridge, MA, 2007.

[18] S. Yu, J. Bi, and J. Ye. Probabilistic interpretations and extensions for a family of 2D PCA-style algo-rithms. In the KDD’2008 Workshop on Data Mining

using Matrices and Tensors, 2008.

[19] S. Zhu, K. Yu, and Y. Gong. Predictive matrix-variate t models. In J. C. Platt, D. Koller, Y. Singer, and S. Roweis eds., Advances in Neural Information

Processing Systems 20, pp. 1721–1728. MIT Press,