DEIM Forum 2016 H7-2

一般化 KL ダイバージェンスを用いた
非負テンソル補完と交通流解析への応用
竹内

孝†

太†

納谷

上田 修功†† NTT コミュニケーション科学基礎研究所 〒 619–0237 京都府相楽郡精華町光台 2-4
E-mail: †{takeuchi.koh,naya.futoshi,ueda.naonori}@lab.ntt.co.jp
あらまし

本研究では, 一般化 KL ダイバージェンスを用いた非負テンソル補完問題 (NTC: Non-negative Tensor

Completion) を提案する. N 一般化 KL ダイバージェンスを用いることで, 観測値に自然数が含まれる場合に適した
NTC のモデリングが可能になる. 提案法の学習法として, 乗法的更新則を用いた因子の推定法を導出する. 交通流
データの欠損値補完実験を行い, 提案法と既存の欠損値補完法の欠損値補完精度を比較する. 実験から欠損値精度が提案法によって改善されることが確認された. さらに推定された因子の可視化を行い、提案法によって得られた因子の定性的な評価を行う.
キーワード

非負テンソル因子分解法，欠損値推定，乗法的更新則, 交通流データ
たテンソル補完とその推定法を提案する. 一般化 KL ダイバー1. イントロダクション

ジェンスは, ポアソン分布の負の対数尤度の定義域を実数に拡交通流データ, 評価点データ, 画像データなど, N 個の属張した距離規準で, 非負の整数がポアソン分布に従う場合に用性に紐付けられた観測値からなる関係データが多く存在す

いると, 高い汎化性能を示すことが知られている [4]. 提案法でる [3, 16, 12]. これらのデータには, センサーの故障や, デー

因子を推定するために、テンソル補完問題をモード毎に因子をタ転送の失敗, データの未観測などによって欠損値が含まれる

推定する非凸最適化問題として扱うと, 複数の非負行列補完問場合がある. このような場合にデータの欠損値の補完が出来れ

題 [14, 15] となることを用いて, 乗法的更新則を導出する.

ば, ある時刻に欠損したセンサーの観測の擬似的な観測, ユー

提案法を交通流データの欠損値推定問題に適応し, 複数の評

ザのあるアイテムに対する未観測の評価点の推定, 一部のピク

価基準から既存の欠損値推定技術との汎化性能の比較実験を行

セルが欠損した画像データの復元などが可能になる.

う. さらにデータセットに含まれる, 駐車場の駐車数データやテンソル [9] とは, 行列を一般化し, 関係データの N 個の属

都市の気象データを利用した場合の実験を行い, 複数のデータ性を N 個モードに対応させることで, 関係データを高次元配を混ぜあわせた場合の欠損値推定実験を行う.

実験から得ら

列として扱うモデルである (N = 2 のとき, テンソルは行列とれた因子の可視化を行い、提案法の解釈容易性を定性的に確認一致する). テンソル分解 [9] は, テンソルの観測値は少数のパする.

ターンで表現可能であると仮定し, パターンを因子として推定

以下, 本稿の構成を述べる. 2. 章では, テンソル補完の導入する技術である. テンソル補完 [5] は, テンソルに欠損値が含と一般化 KL ダイバージェンスを用いた定式化を行う. 3. 章

まれる場合に合わせてテンソル分解を拡張した技術である。テでは, 乗法的更新則を用いた因子推定法を導出する.

ンソル補完は, 少数の因子のみを用いるため, 少数の観測値から

は, 交通流データを用いた実験の設定と実験結果を示す. 5. 章

高精度で欠損値の推定を行えることが知られている.

では, 考察と今後の課題を述べる.

先述した関係データの例では, 観測値が 2 地点を通過した車

4. 章で

2. テンソル補完

両の台数や, ユーザがアイテムに付与した 5 段階の評価点, ピ
クセル毎の RGB の輝度値など, 非負の整数となっている. 観

テンソルおよび一般のテンソル分解については [9, 3] を参照

測値が非負の整数値からなる関係データは, これら以外にも多されたい. モード数が N で, 各モードの特徴数が I1 , . . . , IN の数存在する [7]. 非負の整数値からなる関係データの解析法とN 階非負テンソルを,して, 非負値行列因子分解法 [10, 11, 6] や, 非負値テンソル因I1 ×···×IN
X ∈ R+,

子分解法 [13] などが提案されている. これらの技術は, データ(1)の非負性に合わせて因子に非負制約を課すことで, 解釈性の高とし, テンソルのインデクスとその集合を i = (i1 , . . . , iN ) ∈いパターンを因子として抽出することが知られている.

I,

しかし, テンソル補完法では、観測値が非負の整数値である場合を考慮したモデルは提案されていない.
本研究では, 観測値の定義域が非負の整数の場合のテンソル補完の問題を扱うために, 一般化 KL ダイバージェンスを用いR

テ ン ソ ル X の n モ ー ド ア ン フォー ル ド を X(n) ∈Im ×(I1 ···In−1 In+1 ···IM )とする.

テンソル X の欠損を示すマスクテンソルを M ∈{0, 1}I1 ×···×IN とし, マスクテンソルの要素 mi を,

(3)

'

X

K
X

k=1

3. 推 定 法

ak

(2)

ak
k

る. この際, 非負テンソル補完における因子行列毎の最小化問
題は, 非負行列補完問題に変形できることを用いる [18, 17],

(1)

(2)

(3)

ak , ak , akで近似した時の誤差を最小化するように因子を推定する.

まず, n 番目モードに対応する因子行列 A(n) を変数, その他の因子行列を定数とみなす. 次に対角成分が 1, それ以外の要
R×···×R
素が 0 の N 階対角テンソルを C ∈ R+とすると, 推定値

X̂ は,
1 (if xi is observed),mi =0 otherwiseX̂ = C ×1 A(1) × · · · ×N A(N ) ,
,

(8)

(2)
となる. このとき, ×n はコアテンソル C と因子行列 A(n) の n
モードのマルチプリケーションである.

とする.
本研究では, テンソル X の観測値を K 個の因子の線形和によって近似する CP 分解をテンソル補完のモデルに用いる [1, 8]. テンソルの n 番目のモードに対応する k 番目の因子(n)ベクトル ak(n)((n′ ) )
X̂(n) = A(n) ⊙N,
n′ =1,n′ =n
| A

(9)

番目の因子行列を除いた, 因子行列のクラーチ・ラオ積,

(n)

(3)

と定め, 因子行列の集合を A = {A(1) , A(2) , . . . , A(N ) } とす
る. CP 分解では, テンソルの X の推定値を, 因子の線形和, す
なわち,
K
∑

推定値 X̂ の n モードアンフォールドは,となる. この時, ⊙ はクラーチ・ラオ積とし, ⊙Nはnn′ =1,n′ =n|∈ RI+n を列ベクトルにもつ因子行列を,A(n) = (a1 , . . . , aK ) ∈ RI+n ×K ,x̂i =適化問題である. 本研究では, 因子行列 A(n) (n = 1, . . . , N )を順番に最適化するブロック最適化を用いて因子行列を推定す

(1)

ak

非 負 テ ン ソ ル 補 完 の 模 式 図. テ ン ソ ル X を 因 子

図 1

非負テンソル補完は, 因子行列の集合 A を最小化する非凸最⊙N= A(1) ⊙ · · · A(n−1) ⊙ A(n+1) · · · A(N ) . (10)n′ =1,n′ =n|とした. 以上から,X(n) ∈ RIn ×J , M(n) ∈ RIn ×J ,
(1)

(2)

(N )

ai1 ,k ai2 ,k · · · aiN ,k ,

(4)

J×K
Y = ⊙N,n′ =1,n′ =n| A(n′ ) ∈ R

(11)

k=1

と定める. コスト関数 f を X と X̂ の要素毎の近似誤差の線形とすると, 因子行列 A(n) に関する最小化問題は,

和, すなわち,f (A) =∑minmi d(xi ∥x̂i ),

(5)

i∈I

小化問題,
A

(
m(n),i,j d x(n),i,j

i=1 j=1K∑)(n)ai,k zj,k

.

(12)

k=1

となる. なお, この問題は因子行列 A(n) 関する凸である.

とする. 以上から, CP 分解を用いた非負テンソル補完は次の最min f (A) subject to A

A(n) >
=0

In ∑
J
∑

上記の最小化問題の最適解を乗法的更新則によって求める. 乗
法的更新則は, Expectation-Minimization アルゴリズムを特殊(1) >(N ) >= 0, . . . , A

= 0,

(6)

例として含む, Majorization-Minimization アルゴリズムの一種である. Majorize-Minimization アルゴリズムは, 凸なコスと定められる. N = 3 の場合の非負テンソル補完の模式図を図

ト関数を直接最小化する代わりに, コスト関数の上限を用いて1 に示す.

パラメータの最適化を行うもので, アルゴリズムによるコスト

近似誤差を評価する尺度には, ユークリッド距離や一般化 KL
ダイバージェンス, 板倉斎藤距離などが存在する [3].

本稿では, 近似誤差に一般化 KL ダイバージェンスを用いる. 一般化

KL ダイバージェンスは,d(xi ∥x̂i )
−xi log(x̂i ) + x̂i + const= x̂i + const




+∞

ゴリズムの収束は保証されていない.)
非負行列補完を行うため乗法的更新則は, [14, 15] によって提案されている. この更新則は, NTC の近似誤差に 2 階微分なブ
レグマン・ダイバージェンスを場合に用いることが出来る. ブ
レグマン・ダイバージェンスのヘシアンを要素に持つ行列を Z(if xi > 0 and x̂i > 0),

とおくと, A(n) の乗法的更新則は,(if xi = 0 and x̂i > 0),
A(n) ← A(n) ∗otherwise.
(7)とする [2].

関数の単調減少が保証されている. (ただし, 特殊例を除きアル(Z ∗ X(n) ∗ M(n) )Y T
(Z ∗ X̂(n) ∗ M(n) )Y T

(13)となる. このとき行列の要素ごとの積を ∗, 要素ごとの商を / とした.

	

Algorithm 1 一般化 KL ダイバージェンスを用いた NTC

6

Input: data tensor X , mask tensor M, factor size K, loop

10

count T

10

Initialize A

10

5

4

頻度

for n = 1 to N do
X(n) = unfold(X , n)M(n) = unfold(M, n)

3

10

2

10

end for1

10

for t = 1 to T do

0

10 0
10

for n = 1 to N do(
) ()A(n) ← A(n) ∗ M(n) ∗ X(n) /A(n) Z T / M(n) Z T

1

2

10

3

10

10

[ /30 ]

end forend for

図2

交通流データのヒストグラム
(3)

(n)

一般化 KL ダイバージェンスは 2 階微分が可能で, その ai,kに関するヘシアンは,(∇2a(n) d xi
i,k

K
∑

)
	

(n)

ai,k zj,kk=1

−xi /x̂i + 1


= 1




0

(if xi > 0 and x̂i > 0),(if xi = 0 and x̂i > 0),

図3

(14)

X1
X2
X3

'

K
X

k=1

ak
tk

(2)

ak
(1)

ak

設定 B における非負テンソル補完の模式図. 駐車場と気象のテンソルを交通流のテンソルに結合させて利用する.

otherwise.
R441×48×61 を作成した. 駐車データからは, 駐車場 (8 箇所), 時である. したがって, 因子行列 A(n) の乗法的更新則は,刻, 日付をモードに持つ 3 階のテンソル X2 ∈ R8×48×61 , 気

(

) (
)
A(n) ← A(n) ∗ M(n) ∗ X(n) /A(n) Z T / M(n) Z T

象データからは, センサ (温度, 気温, 湿度, 風速の 4 種類), 時(15)

刻, 日付をモードに持つ 3 階のテンソル X3 ∈ R4×48×61 を作成した.

となる. この更新則を用いて各因子行列を順番に更新する. 計
算手順をアルゴリズム 1 を示す.

4. 実

実験では 2 つの設定を用いる. 設定 A では, テンソル X1 を用いて交通流の欠損値を推定する. 設定 B では, 交通流データに加えて駐車場データ X2 , 気象データ X3 を用いて, 交通流の欠験

損値を推定する. これは駐車場データと気象データが, 交通流

実験では, Citypulse（注 1） という, デンマークのアーヘン市内

データとの間に何かしらの関係を持つならば, 追加のデータにで計測されたスマートシティに関するデータセットを利用す

よって欠損値の補完精度が高まると考えられるためである. 駐る.

車場データと気象データを追加のテンソルとして扱い, 交通流

Citypulse データセットに含まれる, アーヘン市内の 419地点で観測された交通流データに対し欠損値推定を行う.

な

お, Citypulse データセットには 8 個の駐車場の駐車車両の台
数, 市内の気象データ (温度, 気温, 湿度, 風速）も含む.

のテンソルに結合させて X として利用する. 設定 B で用いる
テンソル補完の模式図を図 3 に示す.
データの欠損割合に対する, 欠損値推定の精度の変化を確

データの観測期間は 2014 年 8 月 1 日から 9 月 31 日まで認するために, いずれの設定でも, 交通流データの実際の欠の 61 日間で, この期間中はデンマークの暦に祝日は存在しな損値に加えて観測値をランダムに欠損させた.

い. データは 30 分毎に観測されるため, データが完全であれ10%, 50%, 70%, 90%, 99% を用いた. 各設定で 5 回の試行を行ば, センサー毎に最大で 2, 928 点の観測値が存在する.

い, 各試行で毎回ランダムに観測値を欠損させ, 残った観測値のセン欠損の割合は

サー毎にデータの観測値を調べ, 24 時間データが完全に観測さ90% をトレーニングデータ, 10% をテストデータとした. 因

れていない場合には欠損とした. 交通流データのヒストグラム

子の数は 5, 10, 20 とし, 5-fold cross validation によって決定を図 2 に示す. 図から, データは通貨車両台数が増えるにつれした.

て頻度が減少する, ポアソン分布様の形状をしていることがわ

欠損値の推定精度は, 提案法の距離規準である一般化 KL ダ

かる. 気象データは, 各データの観測値の最大値が 1 となるよ

イバージェンス (gKL), 時空間データの欠損値推定の精度に用

うに前処理を行った.

いられる Relative Mean Absolute Error (RAME),∑i∈I ∥xi − x̂i ∥1∑
RAME =,i∈I ∥xi ∥1

交通流データからは, 観測地点 (441 点), 時刻 (30 分毎の

24 時間), 日付 (61 日) をモードに持つ 3 階のテンソル X1 ∈（注 1）：http://www.ict-citypulse.eu/(16)

一般の欠損値推定の精度に用いられる Root Mean Squared

Error (RMSE),√
RMSE =

1 ∑(xi − x̂i )2 ,|I| i∈I(17)を採用した.
比較手法として, ユークリッド距離を用いた非負テンソル補完(Non-negative Tensor Completion: NTC ), テンソルをセン
サーのモードに関してアンフォールドした行列に対する, 一般化

KL ダイバージェンスを用いた非負行列補完 (Non-negative Matrix Completion: NMC), 欠損値を 0 として扱う非負値テンソ
ル因子分解法 (Non-negative Tensor Factorization with gKL:NTF ), 交通流データの平均値を推定値としたもの (Mean) を採用した.

4. 1 実 験 結 果
設定 A の実験の結果として, 5 回の試行から得られた評価指標の平均値と標準偏差を表 1, 2, 3 に示す. 最も低い平均値を示した手法の数値を太字にした. 欠損率が 99% のとき, NMCと NTF は計算結果が不安定になったため, 結果から除外した.

図4

設定 B で提案手法から抽出された因子の例.

欠損率が 10% の場合には, NMC が他の手法と比べて最も良い汎化性能を全指標で示した. 一方, 欠損率が 30, 50, 70, 90, 99%の場合には NMC の汎化性能が劣化し, 提案手法の汎化性能が最良となった. NTC と NTF は, 平均値を推定値とした Mean
よりも高い汎化性能を示すが, 全設定で提案手法と NMC よりも汎化性能が低かった.
まず, 提案法と NMC が、NTC よりも高い汎化性能を示した
理由として, 一般化 KL ダイバージェンスがデータの分布に近いためと考えられる. NTF が一般化 KL ダイバージェンスを用いたにも関わらず, 汎化性能が低い原因には, 欠損値を適切に扱っていないためである.
欠損率によって提案法と NMC の順位が変化する原因は, パ
ラメータ数の大小による過学習が考えられる. つまり, 因子数∑を K としたとき, 提案法ではパラメータ数が N
n=1 In K であるのに対し, NMC ではテンソルをアンフォールドしたため, パ
ラメータ数が In (Πn′ =n| In′ )K となり, 提案法は NMC よりもパラメータ数が少なくなる. これはテンソル分解では, 日付が異なっていても交通流データのパターンは同一と仮定するためである. 実験では, 因子数 K を 5, 10, 20 から Cross-validation
によって決定した. したがって, データの欠損率が大きくなる
ほど NMC は過学習を起こし汎化性能が劣化したが, 提案手法

気象のデータによる汎化性能の改善は限定的であった.
設定 B で得られた因子の可視化結果を図 4, 5 に示す. 図にはアーヘン市内の地図が描かれており, 地図上のオレンジ色の丸は交通流の観測地点, 緑色の丸は駐車場の地点に対応する. 円の半径は因子における各地点のパラメータの強さを示す.

図

右下部のグラフは上から, 因子における時刻のパラメータの強さ, 日付のパラメータの強さを示している. 横軸は, 24 時間と61 日間に対応している. 最後の棒グラフは, 因子における気象データに対するパラメータの強さを示している.
図 4 で, 時刻グラフは 9 時から値が上昇し, 正午にかけてピー
クを取り減衰している. 日付グラフでは, 5 日おきに値が上昇
し, 2 日後に値が下がる周期性を示している. また, 気象グラフ
では気温が高い値を持っている. したがって, 休日の昼間帯の晴天における交通流と駐車場の混雑パターンが抽出されていると考えられる. 一方, 図 5 で, 時刻グラフは早朝から午前中にピークを取る. 日付グラフでは, 値の高い日が 5 日続き, その後の 2 日間に値が下がる, 先ほどの図 4 とは逆の周期性を示している. また, 気象データでは湿度が高い値を持っている. したがって, 平日の午前中の雨天における交通流と駐車場の混雑パターンが抽出されていると考えられる.

は NMC よりもパラメータ数が少ないために, 過学習を起こしにくかったと考えられる. テンソルを利用することで, 欠損率の高いデータからも多数の因子を学習できるようになったとも解釈できる.
設定 B での実験の結果を表 4, 5, 6 に示す. 設定 B においても, 設定 A と同様の傾向が結果に見られ, 欠損率が低い場合には NMC が, 欠損率が高い場合には提案法が高い汎化性能を示した. 提案法の性能を設定 A と設定 B で比較すると, 欠損率
が 90% と 99% の場合には, 設定 B の一般化 KL ダイバージェ
ンスの値が, 設定 A からの改善が確認できる. しかし, RAMEと RMSE では同様の傾向が見られない. 追加された駐車場と

5. ま と め
本研究では, 一般化 KL ダイバージェンスを用いた非負テンソル補間法と乗法的更新則の提案を行った. 提案手法を, Citipulse
データセットのの交通流データの欠損値補完問題に適応することで, 既存法との欠損値補完精度の比較を行った. データの欠損率が高いな場合において, 提案法の優位性が確認できた. また, 推定された因子の可視化を行い, 提案法の定性的な性質を確認した.
今後の課題としては, まず設定 B において推定精度が設定 A

表 1 実験設定 A: テストデータの推定値に対する gKL.
欠損率

提案法

NTC

NMC

NTF

Mean

10% −49.72 ± 0.01 −49.58 ± 0.01 −49.86 ± 0.03 −46.50 ± 0.01 −37.91 ± 5.89
50% −49.11 ± 0.01 −49.00 ± 0.02 −49.02 ± 0.04 −32.01 ± 0.02 −32.54 ± 1.48
70% −48.66 ± 0.04 −48.49 ± 0.02 −47.89 ± 0.17 −18.80 ± 0.05 −32.07 ± 0.21
90% −46.48 ± 0.04 −46.02 ± 0.15 −31.45 ± 7.84

67.07 ± 0.01

−32.13 ± 0.12

99% −39.01 ± 0.27 −35.63 ± 0.48

2353 ± 96−32.27 ± 0.28

Nan ± Nan

表 2 実験設定 A: テストデータの推定値に対する RAME.
欠損率

提案法

NTC

NMC

NTF

Mean

10%

0.26 ± 0.01

0.26 ± 0.01

0.24 ± 0.01

0.51 ± 0.01

1.27 ± 0.1450%

0.27 ± 0.01

0.28 ± 0.01

0.27 ± 0.01

0.82 ± 0.01

1.33 ± 0.03

70%

0.29 ± 0.01

0.30 ± 0.02

0.34 ± 0.01

0.92 ± 0.05

1.33 ± 0.01

90%

0.43 ± 0.01

0.47 ± 0.01

1.27 ± 0.45

0.99 ± 0.01

1.33 ± 0.0199%

0.99 ± 0.02

1.03 ± 0.05

Nan ± Nan

Nan ± Nan

1.33 ± 0.01

表3
欠損率

実験設定 A: テストデータの推定値に対する RMSE.

提案法

NTC

NMC

10%

8.96 ± 0.04

8.91 ± 0.05

50%

9.38 ± 0.11

9.38 ± 0.08

70%

10.40 ± 0.11

90%

NTF

Mean

8.52 ± 0.12

18.86 ± 0.03

34.24 ± 4.97

9.65 ± 0.19

29.22 ± 0.01

32.37 ± 0.45

10.59 ± 0.12

15.00 ± 3.08

33.13 ± 0.01

32.08 ± 0.13

17.10 ± 1.06

18.93 ± 1.11

392.95 ± 289.05

34.53 ± 0.01

31.96 ± 0.03

99%

28.97 ± 1.03

30.55 ± 4.76

Nan ± Nan

Nan ± Nan

31.92 ± 0.13

欠損率

提案法

表 4 実験設定 B: テストデータの推定値に対する gKL.
NTC

NMC

NTF

Mean

10% −49.66 ± 0.01 −49.38 ± 0.01 −49.88 ± 0.03 −46.50 ± 0.02 −30.11 ± 0.62
50% −49.64 ± 0.01 −49.32 ± 0.01 −49.60 ± 0.02 −34.40 ± 0.02 −28.69 ± 0.11
70% −49.55 ± 0.30 −49.22 ± 0.09 −49.12 ± 0.18 −19.69 ± 0.03 −26.25 ± 0.39
90% −49.37 ± 0.18 −48.20 ± 0.58 −31.94 ± 4.60

30.48 ± 3.86

−16.65 ± 1.21

99% −42.20 ± 1.10

1588 ± 111

29.56 ± 19.52

10.28 ± 7.59

Nan ± Nan

表 5 実験設定 B: テストデータの推定値に対する RAME.
欠損率

提案法

NTC

NMC

NTF

Mean

10%

0.27 ± 0.01

0.29 ± 0.01

0.24 ± 0.01

0.51 ± 0.01

1.53 ± 0.01

50%

0.28 ± 0.01

0.30 ± 0.01

0.28 ± 0.01

0.82 ± 0.01

1.64 ± 0.01

70%

0.30 ± 0.01

0.32 ± 0.01

0.32 ± 0.01

0.92 ± 0.01

1.81 ± 0.01

90%

0.40 ± 0.02

0.48 ± 0.03

1.36 ± 0.26

0.99 ± 0.01

2.52 ± 0.04

99%

1.18 ± 0.02

1.82 ± 0.44

Nan ± Nan

0.99 ± 0.01

5.64 ± 1.28よりも劣化した問題の解決がある. 原因としては, 気象データ(NICT) の委託研究「ソーシャル・ビッグデータ利活用・基盤

が整数値を取らず, ヒストグラムがガウス分布状の形状をして技術の研究開発」により得られたものです．

いることが挙げられる. つまり, ポアソン分布に従うデータとガウス分布に従うデータが混在した場合にあわせた, より柔
軟なモデリングと, その場合の因子の推定法が必要と考えられ

文

献る. また, 設定 B で可視化した結果, 時刻のグラフにおいて値[1] J. D. Carroll and J.-J. Chang. Analysis of individ-にノイズが乗り, パターンがスムーズで無くなっていた. これ

ual diﬀerences in multidimensional scaling via an n-はテンソルは, 特徴量をすべて独立に扱っており, 時刻の前後関

way generalization of

係を考慮していないためと考えられる. 特徴量のもつ構造を用

Psychometrika, 35(3):283–319, 1970.

Eckart-Youngdecomposition.

いた正則化を課したモデリングの提案を, 今後のさらなる課題[2] C. Chaux, P. L. Combettes, J.-C. Pesquet, and V. R.

としたい.

Wajs. A variational formulation for frame-based in-

6. 謝

辞

本研究成果の一部は，国立研究開発法人情報通信研究機構

verse problems. Inverse Problems, 23(4):1495, 2007.
[3] A. Cichocki, R. Zdunek, A. H. Phan, and S. Amari.

表6
欠損率

実験設定 B: テストデータの推定値に対する RMSE.

提案法

NTC

NMC

NTF

Mean

10%

9.41 ± 0.11

10.12 ± 0.12

8.56 ± 0.02

19.19 ± 0.04

33.83 ± 0.2150%

9.63 ± 0.10

10.44 ± 0.14

9.95 ± 0.08

29.67 ± 0.03

35.10 ± 0.05

70% 10.69 ± 0.16

11.39 ± 0.02

13.45 ± 1.57

33.52 ± 0.01

36.89 ± 0.21

90% 15.19 ± 1.10

17.67 ± 1.65

410.63 ± 298.90 35.64 ± 0.01

48.42 ± 0.52

99% 41.77 ± 4.33 111.38 ± 52.38

Nan ± Nan

35.45 ± 0.01 103.48 ± 10.55[10] D. D. Lee and H. S. Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401:788–791, 1999.
[11] D. D. Lee and H. S. Seung. Algorithms for non-negative
matrix factorization. In NIPS, 2000.
[12] Yasushi Sakurai, Yasuko Matsubara, and Christos
Faloutsos. Mining and forecasting of big time-series
data. In Proc. SIGMOD, pages 919–922, 2015.
[13] A. Shashua and T. Hazan. Non-negative tensor factorization with applications to statistics and computer
vision. In Proc. ICML, 2005.
[14] S. Sra and I. S. Dhillon. Generalized nonnegative matrix approximations with Bregman divergences.

In

NIPS, 2005.
[15] R. Tandon and S. Sra. Sparse nonnegative matrix ap図5

設定 B で提案手法から抽出された因子の例.

proximation: new formulations and algorithms. Max
Planck Institute for Biological Cybernetics, 2010.

Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blindsource separation. Wiley, 2009.

[16] Yilun Wang, Yu Zheng, and Yexiang Xue. Travel time
estimation of a path using sparse trajectories. In Proc.
SIGKDD, pages=25–34, year=2014,.

[4] P. L. Combettes and J.-C. Pesquet. Proximal splitting
methods in signal processing. In Fixed-point algorithms[17] Y. Xu.

Alternating proximal gradient method forfor inverse problems in science and engineering, pages

sparse nonnegative Tucker decomposition. Mathemat-

185–212. Springer, 2011.

ical Programming Computation, pages 1–32, 2013.

[5] S. Gandy, B. Recht, and I. Yamada. Tensor comple-[18] Y. Xu and W. Yin. A block coordinate descent methodtion and low-n-rank tensor recovery via convex opti-

for regularized multiconvex optimization with applica-

mization. Inverse Problems, 27(2):025010, 2011.

tions to nonnegative tensor factorization and comple-[6] E. Gaussier and C. Goutte. Relation between PLSA
and NMF and implications. In SIGIR, 2005.
[7] N. Gillis. The why and how of nonnegative matrix
factorization. arXiv preprint arXiv:1401.5226, 2014.
[8] R. A. Harshman. Foundations of the PARAFAC procedure: models and conditions for an ”explanatory”
multimodal factor analysis. UCLA Working Papers in
Phonetics, 1970.
[9] T. G. Kolda and B. W. Bader. Tensor decompositions
and applications. SIAM Review, 51(3):455–500, 2009.

tion. SIAM Journal on Imaging Sciences, 6(3):1758–
1789, 2013.