IBIS2016 企画セッション「統計理論」

低ランクテンソルの学習理論と計算理論
1,2,3

鈴木大慈

共同研究者：1 金川平志郎，4 清水伸幸，4 小林隼人，4 田頭幸浩
1

東京工業大学数理・計算科学系
2
JST さきがけ
3
理研 AIPC
4
Yahoo! Japan

2016 年 11 月 16 日
IBIS2016@京都大学吉田キャンパス

1 / 59Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

2 / 59高次元パラメータ推定

ベクトル
スパース性
手法

Lasso
Sure Screening
応用
特徴選択
遺伝子データ解析
3 / 59高次元パラメータ推定

ベクトル
スパース性
手法

Lasso
Sure Screening
応用

行列
低ランク性
手法

PCA
トレースノルム正則化
応用

特徴選択

次元削減

遺伝子データ解析

推薦システム三層ニューラルネット
3 / 59高次元パラメータ推定

ベクトル
スパース性
手法

Lasso
Sure Screening
応用

行列
低ランク性
手法

PCA

テンソル
低ランク性
本発表
高次の関係性を記述

トレースノルム正則化
応用

特徴選択

次元削減

遺伝子データ解析

推薦システム三層ニューラルネット
3 / 59概要

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

4 / 59

低ランク行列推定問題
スパースベクトルの推定：
スパースな係数の学習(基底は固定)低ランク行列の推定：

基底の学習 ＋ 係数の学習係数がスパースになるような基底を学習
例：主成分分析 (PCA)
Y = [y1 , . . . , yn ] ∈ Rp×n : 観測値
p
n ∑
∑

min

U∈Rp×k ,V ∈Rn×k⊤ 2(Yj,i − Ui,: Vi,:) = min ∥Y − UV ⊤ ∥2F

i=1 j=1

yi

Y

≒U

VT
score

principal component
5 / 59

推薦システム

Movie A

Movie B

Movie C···

Movie X

User 1

4

8

*

···

2

User 2

2

*

2

···

*

User 3

2

4

*

···

*

6 / 59

推薦システムランク 1 (低ランク)
Movie A

Movie B

Movie C···

Movie X

User 1

4

8

4

···

2

User 2

2

4

2

···

1

User 3

2

4

2

···

1

N
Movie

4
2
2

8
4
4

4
2 =M
2

min
U,V

∑

⊤ 2
(Yi,j −Ui,: Vj,:)(i,j)∈I

User(e.g., Srebro et al. (2005); Candès and Tao (2009); Candès and Recht(2009), NetFlix (Bennett & Lanning, 2007))
6 / 59３層ニューラルネットワーク

f (x) = V ⊤ Ux
f (x) = V ⊤ h(Ux)7 / 59

低ランク行列推定の推定精度真の行列 A∗ は M × N 行列でランク d とする (d ≪ M, N: 低ランク).
最小二乗法:∥Â − A∗ ∥22 = Op

(

MN
n

)

低ランク性を利用できていない．
トレースノルム正則化:()
d(M + N) log(MN)∥Â − A∗ ∥22 = Opn
ただし，∥A∥2 :=

√

d

N
d

M

1
MN ∥A∥F .

(Negahban et al., 2012; Negahban & Wainwright, 2012; Rohde & Tsybakov, 2011;
Koltchinskii, 2012)
(最小二乗法)



MN
n

≫

d(M + N) log(MN)n(正則化学習)
8 / 59概要

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

9 / 59行列からテンソルへ

User

User

Movie

行列
∑d
(1) (2)
Xij = r =1 ur urContext−→

Movieテンソル (多次元配列)∑d(1) (2) (3)Xijk = r =1 ur ur ur

複数項目間の高次の関係性 (K 次)
いわゆる座標返還に関して不変な量という意味は薄れている．しかし，多重線形写像という観点は重要．
さらに高次元 M1 × M2 × · · · × MK
→ 低ランク性 を利用テンソルのランクとは？

10 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

11 / 59テンソルの分解・ランク

CP-分解(Hitchcock, 1927; Hitchcock, 1927; Carroll & Chang, 1970; Harshman, 1970)

Tucker-分解 (Tucker, 1966)
テンソルトレイン分解 (Oseledets, 2011)12 / 59

Tensor rank: CP-rank
c1

=

b1

a1

c2

+

b2

a2

cd

+…+
bd

C

=

B
A

ad

CP-分解
Canonical Polyadic 分解 (Hitchcock, 1927; Hitchcock, 1927)
CANDECOMP/PARAFAC (Carroll & Chang, 1970; Harshman, 1970)∑d∑d
Xijk = r =1 air bjr ckr (X =: [[A, B, C ]]). X = r =1 ar ⊗ br ⊗ cr とも書く．
CP-分解で CP-ランクが決まる．
CP-分解の計算は NP-困難 (Hillar & Lim, 2013).
直交分解できるとは限らない．
ランク２のテンソルでランク３のテンソルをいくらでも良く近似できる例
が作れてしまう．
13 / 59

CP-分解の一意性
CP-分解が一意であるための十分条件
U(3)X = [[U (1) , U (2) , . . . , U (K ) ]]各 U (k) が 列 フ ル ラ ン ク な ら ，CP-分 解 はスケーリングと列の入れ替えの自由度を除いて一意．
（行列のような回転の自由度はない）= U(1)

U(2)14 / 59

CP-分解の一意性
CP-分解が一意であるための十分条件
U(3)X = [[U (1) , U (2) , . . . , U (K ) ]]各 U (k) が 列 フ ル ラ ン ク な ら ，CP-分 解 はスケーリングと列の入れ替えの自由度を除いて一意．
（行列のような回転の自由度はない）= U(1)

U(2)k-rank(A)：Kruskal ランク，行列 A の劣を任意に r 個とってきたときに線形独立である r の最小値．

より詳細な条件
X の CP-ランクが d の時，以下が満たされていれば CP-分解は一意:
2d + (K − 1) ≤K∑

k-rank(U (k) ).

k=1
14 / 59

CP-分解の計算

一般的には CP-分解は NP 困難だが．
．

X = [[U (1) , U (2) , . . . , U (K ) ]]で，各 U (k) が列フルランク とする． (このとき，CP-分解は一意)- すると，多項式時間で CP-分解は計算できる．
- 有限回の固有値分解や行列の同時対角化などで解ける．
- 確率的アルゴリズム
Simultaneous Diagonalization (De Lathauwer, 2006)
Simultaneous Generalized Schur Decomposition (De Lathauwer et al., 2004)
Generalized Eigenvalue Decomposition (Leurgans et al., 1993)

15 / 59

Tucker-ランク

=
X

C

G

B
A

Tucker-分解 (Tucker, 1966)
Xijk =∑r1 ∑r2
l=1

m=1

∑r3

n=1 glmn ail bjm ckn . X =: [[G ; A, B, C ]].

G はコアテンソルと呼ばれる．
G が対角なら CP 分解に帰着．
Tucker-ランク = (r1 , r2 , r3 )． ただし rk は G の各辺の大きさ．これは，
モード k 展開 (次ページ) のランクで与えられる．
16 / 59

行列化
モード k 展開: A(k) ∈ RMk ×N/Mk ,

A

(N =

∏K

k=1 Mk ).

A(k)rk = rank(A(k) ) は Tucker-ランクを与える．

17 / 59

応用
推薦システム（ユーザー × アイテム × 店）
マルチタスク学習
信号処理（空間 (2D) × 時間）
自然言語処理（単語のベクトル表現）

Context

4
4 2 3
3 4 Rating
414
1 21 4
1
Item
2 22
2
1
3 3 1 1
3 33 2 1 4 1
2 2
22
Prediction 4
1
3
11 3
2
2
User
テンソル補完18 / 59その他の応用
脳波解析 (De Vos et al., 2007)

時間 × 周波数 × 空間
EEG monitoring: Epileptic seizure onset localization
雑音除去 (Phien et al., 2016)

種々の方法による雑音除去

画像のテンソルトレイン化
19 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

20 / 59テンソル推定の回帰としての定式化


テンソルへの回帰:

Yi = ⟨ Xi , |{z}A∗ ⟩ + Wi
|{z}
|{z}
応答

真
入力
(テンソル) (テンソル)A∗ , Xi ∈ RM1 ×...MK : テンソル.

⟨Xi , A∗ ⟩ :=∑Xi,(j1 ,...,jK ) A∗j1 ,...,jK

j1 ,...,jK

テンソル間の内積

仮定: A∗ は「低ランク」.
E.g., Xi = ej1 ⊗ ej2 ⊗ · · · ⊗ ejK なら，テンソル補完の問題に帰着．
21 / 59

∑⟨Xi , A∗ ⟩ := j1 ,...,jK Xi,(j1 ,...,jK ) A∗j1 ,...,jK
特殊ケース: for uk ∈ RMk ,⟨u1 ⊗ · · · ⊗ uK , A∗ ⟩ =

···

j1 =1

<
特に,

M1
∑

<

u1,j1 u2,j2 . . . uK ,jK A∗j1 ,...,jK .

jK =1

<

,

MK
∑

=

<

⟨ej1 ⊗ ej2 ⊗ · · · ⊗ ejK , A∗ ⟩ = A∗j1 ,...,jK .
i

j

k

,

=

→ テンソル補完(i,j,k)22 / 59

低ランクテンソルの求め方

凸正則化法
低ランクテンソルの集合上での訓練誤差最小化
交互最適化法Riemannian 最適化
ベイズ推定法テンソル推定における統計理論構築の難しさ
高次元性
特異性 (パラメータ表現の一意性のなさ)
非凸性→ 単純な統計的漸近理論が使えない．
これらを考慮したバウンドを構築する必要がある．

23 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

24 / 59

凸正則化による推定
凸正則化アプローチ:

1∑(Yi − ⟨Xi , A⟩)2 + pen(A).
M1 ×M2 ×···×MK n
A∈R
n

min

i=1

テンソルを行列に変形してトレースノルム正則化
Overlapped Schatten 1-norm (Tomioka et al., 2011)
Latent Schatten 1-norm (Tomioka & Suzuki, 2013; Zheng & Tomioka, 2015)
Square deal (Mu et al., 2014)
Scaled latent regularization (Wimalawarne et al., 2014)d1 b∥A − A∗ ∥22 ≤ C
N



√∏


k Mk

n


Overlapped Schatten 1-norm|||A|||S1 /1 :=

A

A

(k)

K
∑

∥A(k) ∥Tr

k=1
25 / 59

制限強凸性


制限強凸性（要約）:δ
1∑∥A∥22 ≤⟨Xi , A⟩2
N
n
n

i=1

がすべてのランク r のテンソル A に対して成り立つ.


線形回帰における各種制限強凸性に類する条件の関係(Bühlmann & Van De Geer, 2011)26 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

27 / 59

ベイズ推定(Suzuki, 2015)

Ab =∫AdΠ(A|X1:n , Y1:n ).

CP-ランク d ′ のテンソルの事前分布:
ランク d ′ テンソル A の分解′
A = [[U (1) , . . . , U (K ) ]] (U (k) ∈ Rd ×Mk )．
π(U (1) , . . . , U (K ) |d ′ )
}
{
K
d′ ∑
(k) ⊤ (k)
Tr[U
U ] ,
∝ exp − 2
2σp

U(3)

= U(1)

U(2)k=1

ランクの事前分布: ある定数 1 > ξ > 0 を用いて
∑d max′π(d ′ ) ∝ ξ d(π(d ′ ) = 1).
′
d =1



∑K
)
(
d( k=1 Mk )1 b
∗ 2
∥A − A ∥2 ≤ C
log KnR K
Nn※制限強凸性条件は不必要．

28 / 59誤差上界の比較

手法

最小二乗法

∏K
誤差上界

k=1 Mkn行列化による
凸正則化
√∏
d
k Mk
n

大

ベイズ

∑
d( Kk=1 Mk ) log(n)n小ミニマックス最適
up to O(log(n))(Suzuki, 2015)

テンソルと行列の間には大きなギャップがある．

29 / 59

数値実験

0.9

8

setting 2: Bayes
setting 2: Convex
setting 5: Bayes
setting 5: Convex

0.8

0.7

6

Out-of-sample acc.

In-sample acc.

0.6

Bayes

0.5

0.4

0.3

0.2

0.1

0

0.4

0.5

0.6

ns

0.7

0.8

setting 2: Bayes
setting 2: Convex
setting 5: Bayes
setting 5: Convex

7

0.9

Convex reg.
via matricization

5

4

3

2

1

0

0.4

0.5

0.6

ns

0.7

0.8

0.9

Figuur: ベイズ推定法と凸正則化 (overlapped approach) の比較.

30 / 59

ミニマックス最適性

状況：テンソル補完
ＣＰ-ランク d かつ max-ノルムが R 以下のテンソルの集合:{}
Hd (R) := A ∈ RM1 ×···×MK | CP-rank d, ∥A∥max,2 ≤ R .

Xi = ei1 ⊗ ei2 ⊗ · · · ⊗ eiK かつ (i1 , . . . , iK ) が一様.

Theorem
Hd (R) 上のミニマックスリスクは下から次のように抑えられる:}{d(M1 + · · · + MK ) 2K
min ∗max E[∥Ab − A∗ ∥2L2 ]≳ min σ 2,R.
b A ∈Hd (R)nA

31 / 59

概要

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

32 / 59

ベクトルから関数へ

f (x (1) , x (2) ) = x (1)⊤ Ax (2) =d∑⟨x (1) , ur(1) ⟩⟨ur(2) , x (2) ⟩

r =1

x

(1)

: ユーザの特徴，x (2) : 映画の特徴．
33 / 59

ベクトルから関数へ

f (x (1) , x (2) ) =d∑

fr(1) (x (1) )fr(2) (x (2) )

r =1

x

(1)

: ユーザの特徴，x (2) : 映画の特徴．
33 / 59

ベクトルから関数へ: 平滑化(De Vos et al., 2007)min(k){ur }r ,k

u1

X−

d
∑

ur(1) ⊗ ur(2) ⊗ ur(3)r =1

u2 u3 u4

(Cao et al., 2016)
2

+

d ∑
3
∑
r =1 k=1u ⊤ Gu =ur(k)⊤ Gur(k)| {z }平滑化∑

2j (uj − uj+1 )u5t
空間的な平滑化 ⇒ カーネル法(Zdunek, 2012; Yokota et al., 2015a; Yokota et al., 2015b)
34 / 59

非線形テンソルモデル
非線形回帰モデル

yi = f (xi ) + ϵi .
Goal: f をデータ Dn = {xi , yi }ni=1 から推定したい.
(k)

非線形テンソルモデル: fr∈ H (再生核ヒルベルト空間)

f (x (1) , . . . , x (K ) ) =d∑

fr(1) (x (1) ) × · · · × fr(K ) (x (K ) )r =1

線形テンソルモデル:

f (x (1) , . . . , x (K ) ) =d∑⟨x (1) , ur(1) ⟩ × · · · × ⟨x (K ) , ur(K ) ⟩r =1

行列の場合:∑dr =1 ⟨x

(1)

d
∑

(1)
(2)
, ur ⟩⟨x (2) , ur ⟩ = (x (1) )⊤ (

|

r =1

⊤

ur(1) ur(2) ) x (2) .
{z}matrix(x (k) )k は必ずしもユークリッド空間の元である必要はない．
35 / 59

推定方法

凸正則化法 (Signoretto et al., 2013), 少ない計算量.
/ （おそらく）ミニマックス最適性を満たない．
交互最適化法 (Suzuki et al., 2016) (本研究), 凸最適化の繰り返しで解ける. 少ない計算量.
/ ミニマックス最適性を満たすには強い条件が必要.
/ 局所最適解の問題.
ベイズ推定法 (Kanagawa et al., 2016) (本研究), 統計的な性質が良い. ミニマックス最適.
/ 重い計算量.

36 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

37 / 59交互最適化法
1
F ({fr(k) }r ,k ) :=
n



n (
∑
i=1

)2

∗

yi −

d ∏
K
∑

(k)

fr(k) (xi )(経験誤差)

r =1 k=1

(k)
ある (r , k) において他の成分を固定したまま fr のみを更新:{}(k ′ )fˆr(k) = arg min F (fr(k) |{fˆr ′ }(r ′ ,k ′ )̸=(r ,k) ) + λ∥fr(k) ∥2H .
(k)

fr

∈H


目的関数は非凸関数.
(k)
各成分 fr に関しては 凸関数.
局所最適解へは収束 (e.g. 座標降下).

38 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

39 / 59

ノンパラメトリックベイズ推定法
非線形テンソルモデル:

f (x (1) , . . . , x (K ) ) =d∑

fr(1) (x (1) ) × · · · × fr(K ) (x (K ) )

r =1
(1)
(K )
yi = f (xi , . . . , xi ) + ϵi(ϵi ∼ N(0, σ 2 )).

再生核ヒルベルト空間: カーネル関数 k(x, x ′ )⇔ ガウス過程: 共分散関数 Cov(f (x), f (x ′ )) = k(x, x ′ )(k)

ガウス過程事前分布を各成分 fr
事後平均がベイズ推定量:

に置く.

∫

fˆ(x) =

f (x) p(df |Dn )p(df |Dn ) ∝| {z }p(Dn |f )| {z }

事後分布p(df ) =∏r ,kp(df )| {z }

データ Dn の尤度 事前分布p(df (k) )| {zr }
ガウス過程事前分布

40 / 59

ノンパラメトリックベイズ推定法
非線形テンソルモデル:

f (x (1) , . . . , x (K ) ) =d∑

fr(1) (x (1) ) × · · · × fr(K ) (x (K ) )

r =1
(1)
(K )
yi = f (xi , . . . , xi ) + ϵi(ϵi ∼ N(0, σ 2 )).

非線形テンソルモデルの事前分布:
各 (r , k) ∈ [d] × [K ] において,λ(k)∼ Exp (1)r

(スケーリングパラメータ)

fr(k) ∼ GPr ,k (·|λ(k)
r )
(k)

(k)

ただし，GPr ,k (·|λr ) はスケーリングされたカーネル関数 kr ,k /λrされるガウス過程.

で定義

ランクの事前分布: π(d) ∝ ξ d (0 < ξ < 1).
計算: Gibbs サンプリング + Nyström 近似
40 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

41 / 59

学習精度（汎化誤差）の評価

42 / 59

RKHS の複雑さ：カバリングナンバー
カバリングナンバー N(B(H), ϵ, L2 (P)): L2 (P)-ノルムに関する半径 ϵ の球で単位球 B(H) を覆うために必要な最小の個数.

B(H)カバリングナンバー条件 (s)ある 0 < s < 1 と定数 C が存在して，以下の不等式がすべての ϵ > 0 で成り立っていると仮定する.
log N(B(H), ϵ, L2 (P)) ≤ C ϵ−2s .

s が小さければ「単純」，大きければ「複雑」
43 / 59

小球確率とカバリングナンバー条件
再生核ヒルベルト空間のカバリングナンバー条件は対応するガウス過程の小球確率の評価を与える．

Proposition (Kuelbs and Li (1993); Li and Shao (2001))− log(GP({f | ∥f ∥L2 (P) ≤ ϵ})) ∼ ϵ−2s/(1−s)⇔ log N(B(H), ϵ, L2 (P)) ∼ ϵ−2s(小球確率).

ε

再生核ヒルベルト空間: カーネル関数 k(x, x ′ )⇔ ガウス過程: 共分散関数 Cov(f (x), f (x ′ )) = k(x, x ′ )44 / 59交互最適化法の収束
仮定

f ∗ は incoherence 条件を満たしている (定義は次のスライド).
P(X ) = P(X1 ) × · · · × P(XK ).
∗(k)
fr
∈H
その他，細かい条件.
fˆ[t] : t-回更新した後の解. s: 再生核ヒルベルト空間 H の複雑さ．

Theoremd ∗ , K に依存した定数 C1 , C2 が存在して，初期解が d(fˆ[0] , f ∗ ) ≤ C1 を満たすな
らば確率 1 − δ で次を満たす:
(
( )t
)
1
3
[t]
∗ 2
∗
∗
− 1+s
ˆ
∥f − f ∥L2 (P) ≤ C2 d K+d| Kn{z } log(1/δ) .
4| {z }
推定誤差
最適化誤差

局所最適解へ線形収束 (log(n) 回の更新で十分
()
)
1
[0]
∗
[log(n)]∗d(fˆ , f ) = Op (1) =⇒ d(fˆ, f ) = Op d ∗ Kn− 1+s .
「ある程度近いところから始めれば，最適な推定誤差を達成」
45 / 59

条件の詳細Incoherence: ∃µ < 1 s.t.
∗(k)|⟨fr∗(k) , fr ′∗(k)⟩| ≤ µ∥fr∗(k) ∥L2 ∥fr ′(r ̸= r ′ ).

∥L2

(k)

fr*

(k)

fr'*

各成分のノルムに関する上界，下界:

0 < vmin ≤ ∥fr∗(k) ∥L2 ≤ vmax(∀r , k).

sup-ノルムに関する条件: 0 < ∃s2 < 1 s.t.
s22∥f ∥∞ ≤ C ∥f ∥1−s
L2 ∥f ∥H

vr = ∥(∀f ∈ H)

∏K
(k)
∗∗(k)
∗(k)
fr∗(k)= ∗(k)∥L2 , v̂r = ∥ k=1 fˆr ∥L2 , fr
k=1 fr
∥f
∥

∏K

r

ˆ(k)
, fˆr =L2fˆr(k)
(k)

∥fˆr

∥L 2

，

ˆr(k) − fr∗∗(k) ∥L }.
d(fˆ, f ∗ ) = max{|v̂r − vr | + vr ∥fˆ
2(r ,k)

46 / 59

収束理論の図解

Pred. Error

Pred. Risk

Emp. Error

Emp. Risk

True

観測データ小

True

Pred. Risk
Emp. Risk

True

観測データ大予想誤差は真の関数の近傍で局所的に凸関数のような形をしている．
経験誤差はサンプルサイズが大きくなるにつれ，予測誤差に近くなってゆく．

47 / 59

証明に用いた道具
Rademacher 複雑度
H: 関数集合[]n
1∑
Ex suph(xi ) − E[h]|{z}h∈H n i=1
| {z } 予測誤差
経験誤差[
≤2Ex,σ1∑sup
σi h(xi )
h∈H n
n

]

i=1

C
≤√
n

Pred. Risk

Emp. Risk

Uniform bound

ただし，{σi }ni=1 は i.i.d. の Rademacher 確率変数
(P(σi = 1) = P(σi = −1)).

局所 Rademacher 複雑度 + peeling device: 二乗損失の強凸性を利用[
Ex,σ sup
h∈H

| n1

]
s
1
σ
h(x
)|
λ− 2
λ− 2
i
i
i=1
≤C √ ∨ − 1∥h∥L2 + λ
n
n 1+s

∑n

真の周辺でよりタイトなバウンド

Pred. Risk

Emp. Risk

Uniform bound48 / 59交互最適化の収束

relative MSE

10 0
10 -1
10 -2
10 -3
10 -4
0

n=400
n=800
n=1200n=1600n=2000n=2400n=2800
5

10

15

20

Number of iterations

25

MSE E[∥fˆ[t] − f ∗ ∥2 ] v.s. 更新回数 t
複数のサンプルサイズ n で比較.
49 / 59

ベイズ推定の理論

50 / 59

事後平均の収束
fˆ: 事後平均fˆ =∫f Π(df |(x1 , y1 ), . . . , (xn , yn )).

∗f : 真の関数f ∗ (x (1) , . . . , x (K ) ) =∑n∥f ∥2n = i=1 f (xi )2 .

∑d ∗∗(1) (1)∗(K ) (K )(x ) × · · · × fr
(x ).
r =1 fr

Theorem(k)∗

もし，真の関数が fr
推定誤差は[
EY1:n |x1:n∥fˆ − f ∗ ∥2n

]

(k)∗

∈ Hr ,k (∀(r , k)) かつ maxr ,k ∥fr{≤ C (R + 1)∗

2(K −1)d ∑K∑
r =1 k=1n∥∞ ≤ R を満たすなら，d∗(r ,k) +logn− 1+s 1

(

1
κξ

)}
,

と抑えられる．ただし，0 < s(r ,k) < 1 はカバリングナンバー条件の指数:

log N(BHr ,k , ϵ, ∥ · ∥n ) = O(ϵ−2s(r ,k) ),ただし，BHr ,k は RKHSHr ,k の単位球で N(BHr ,k , ϵ, ∥ · ∥n ) はそのカバリングナ
ンバー．
51 / 59

例: Matérn prior
Xk = [0, 1]p(k)∗fr∈ W α,2 [0, 1]p : ソボレフ空間 (α-回微分可能で，その微分係数が L2 )ナイーブな方法: テンソル分解の構造を使わず真の関数 f ∗ を空間全体
X = X1 × · · · × XK 上でそのまま推定.

n− 2α+K p .
2α/ 次元の呪い本研究の方法:d ∗ K n− 2α+p .
2α, 次元の呪いを回避
52 / 59

ガウス過程回帰のリスクバウンド
観測ノイズに適当な仮定を置くことで, 次の定理を得る．

Theorem (PAC-ベイズバウンド (Dalalyan & Tsybakov, 2008))
事前分布 Π に対して絶対連続なすべての確率測度 ρ に対して，以下の不等式が
成り立つ．
[
] ∫
βK(ρ, Π)∗ 2
ˆEY1:n |x1:n ∥f − f ∥n ≤ ∥f − f ∗ ∥n2 dρ(f ) +.
n
平行移動した集合の測度 ガウス測度 µ とすべての原点対称なボレル集合 Aについて
1− log(µ(A + h)) ≤ ∥h∥2Σ−1 − log(µ(A)).
2
これらを合わせると,E[∥fˆ − f ∗ ∥2n ] ≲ ϵ2n +∥f ∗ ∥2H
1− log(GP({f | ∥f ∥2n ≤ ϵ2n }))
{z
}
n
n|
小球確率

53 / 59

ミニマックス最適性
導かれた上界は ミニマックス最適性 を（定数倍を除いて）満たす.
ランク d ∗ のテンソルの集合:{
H(d ∗ ,K ) (R) :=

}

∗

f =

d ∏
K
∑

fr(k)

fr(k) ∈ H(r ,k) (R).

r =1 k=1

Theorem (ミニマックスリスク)

inf

sup

E[∥f ∗ − fˆ∥2L2 (PX ) ] ≳ d ∗ Kn− 1+s ,1fˆ f ∗ ∈H(d ∗ ,K ) (R)ただし，E[·] は訓練データの出方に関する期待値，inf はすべての推定量 fˆ についてとる.
ベイズ推定量は定数倍を除いてミニマックス最適性を満たす.
54 / 59

Outline

1

低ランク行列推定2

線形テンソル推定テンソルのランクテンソル推定の定式化
凸正則化による方法
ベイズ推定法3

非線形テンソル推定
交互最適化法
ベイズ推定法
ノンパラメトリックテンソル推定の推定理論
数値実験: マルチタスク学習

55 / 59数値実験 (実データ)
マルチタスク学習 [非線形回帰]
Restaurant data: 138 人の客× 3 項目 (料理, サービス, 全体)
レストランの評価 (3 段階) を予測
レストランは 44 次元の特徴ベクトルで表現 → レストランに対して評価を予測する関数を客と項目ごとに推定．
∑df (客, 項目, xレストラン ) = r =1 α客,r β項目,r fr (xレストラン )56 / 59数値実験 (実データ)
マルチタスク学習 [非線形回帰]
Restaurant data: 138 人の客× 3 項目 (料理, サービス, 全体)
レストランの評価 (3 段階) を予測
レストランは 44 次元の特徴ベクトルで表現 → レストランに対して評価を予測する関数を客と項目ごとに推定．
∑df (客, 項目, xレストラン ) = r =1 α客,r β項目,r fr (xレストラン )
1.8

GRBF
GRBF(2)+lin(1)GRBF(1)+lin(2)linearscaled latent
ALS (Best)ALS (50)

1.6
1.4
MSE

1.2
1
0.8
0.6
0.4
500

1000

1500
2000
Sample size

2500

56 / 59数値実験 (実データ)
マルチタスク学習 [非線形回帰]
Restaurant data: 138 人の客× 3 項目 (料理, サービス, 全体)
レストランの評価 (3 段階) を予測
レストランは 44 次元の特徴ベクトルで表現 → レストランに対して評価を予測する関数を客と項目ごとに推定．
∑df (客, 項目, xレストラン ) = r =1 α客,r β項目,r fr (xレストラン )

交互最適化
線形
線形+非線形
線形+非線形
交互最適化
ベイズ

非線形な手法

ベイズ法が若干良い．交互最適化もそれなりに良い精度． 56 / 59

オンラインショッピング購買予測
Yahoo!ショッピングからデータの購買データを取得
ショップ × 商品 × 客 (508 ショップ, 100 商品)
その店でその客がその商品をいくつ買うかを予測．
客は 65 次元特徴ベクトルで表現．ガウスカーネルを利用．
ショップも特徴ベクトルで表現し，類似度を構成．
16
15

MSE

14
13
12
11
10

GP-MTL(cosdis)
GP-MTL(cossim)
AMP(cosdis)
AMP(cossim)
4000

6000

8000 10000 12000 14000
Sample size

Figuur: Yahoo! shopping の購買予測．店間の異なる類似度で比較．
57 / 59

各種手法の比較

手法
非線形性
計算量
統計的性質

行列型正則化法
線形
◎
×

線形ベイズ
線形
×
◎

交互最適化
非線形
○
○

非線形ベイズ
非線形×◎

計算量と統計的性質にはトレードオフがある．

58 / 59

まとめ
線形テンソルモデルと非線形テンソルモデルで（ほぼ）最適な推定方法を導出．
ミニマックス最適レートを導出.
統計的性質と計算量にはトレードオフの関係があった.
手法
誤差上界

最小二乗法
∏K
k=1 Mkn線形モデル
行列化による凸正則化
√∏
d
k Mk
n

ベイズ
∑
d( Kk=1 Mk ) log(n)nノンパラメトリックベイズ推定量の予測誤差バウンドE[∥fˆ − f ∗ ∥2n ] ≤ Cd ∗ Kn− 1+s .
1

ミニマックス最適レート
min ∗ maxfˆf ∈H(d ∗ ,K ) (R)E[∥fˆ − f ∗ ∥2L2 ] ≳ d ∗ Kn− 1+s .
159 / 59

Barak, B., & Moitra, A. (2015). Tensor prediction, rademacher complexity andrandom 3-xor. arXiv preprint arXiv:1501.06521.
Bennett, J., & Lanning, S. (2007). The netflix prize. Proceedings of KDD Cup
and Workshop 2007.
Bhojanapalli, S., & Sanghavi, S. (2015). A new sampling technique for tensors.
arXiv preprint arXiv:1502.05023.
Bühlmann, P., & Van De Geer, S. (2011). Statistics for high-dimensional data:
methods, theory and applications. Springer Science & Business Media.
Candès, E., & Tao, T. (2009). The power of convex relaxations: Near-optimal
matrix completion. IEEE Transactions on Information Theory, 56, 2053–2080.
Candès, E. J., & Recht, B. (2009). Exact matrix completion via convex
optimization. Foundations of Computational Mathematics, 9, 717–772.
Cao, W., Wang, Y., Sun, J., Meng, D., Yang, C., Cichocki, A., & Xu, Z. (2016).
Total variation regularized tensor RPCA for background subtraction fromcompressive measurements. IEEE Transactions on Image Processing, 25,
4075–4090.
Carroll, J. D., & Chang, J.-J. (1970). Analysis of individual diﬀerences in
multidimensional scaling via an n-way generalization of eckart-young
decomposition. Psychometrika, 35, 283–319.
Dalalyan, A., & Tsybakov, A. B. (2008). Aggregation by exponential weighting
sharp PAC-Bayesian bounds and sparsity. Machine Learning, 72, 39–61.
59 / 59

De Lathauwer, L. (2006). A link between the canonical decomposition in
multilinear algebra and simultaneous matrix diagonalization. SIAM journal on
Matrix Analysis and Applications, 28, 642–666.
De Lathauwer, L., De Moor, B., & Vandewalle, J. (2004). Computation of thecanonical decomposition by means of a simultaneous generalized schur
decomposition. SIAM journal on Matrix Analysis and Applications, 26, 295–327.
De Vos, M., Vergult, A., De Lathauwer, L., De Clercq, W., Van Huﬀel, S.,
Dupont, P., Palmini, A., & Van Paesschen, W. (2007). Canonicaldecomposition of ictal scalp eeg reliably detects the seizure onset zone.
NeuroImage, 37, 844–854.
Hardt, M., & Wootters, M. (2014). Fast matrix completion without the condition
number. Proceedings of The 27th Conference on Learning Theory (pp.
638–678).
Harshman, R. A. (1970). Foundations of the PARAFAC procedure: Models andconditions for an “explanatory” multi-modal factor analysis. UCLA Working
Papers in Phonetics, 16, 1–84.
Hillar, C. J., & Lim, L.-H. (2013). Most tensor problems are np-hard. Journal of
the ACM (JACM), 60, 45.
Hitchcock, F. L. (1927). Multilple invariants and generalized rank of a p-way
matrix or tensor. Journal of Mathematics and Physics, 7, 39–79.
59 / 59

Jain, P., Meka, R., & Dhillon, I. S. (2010). Guaranteed rank minimization via
singular value projection. Advances in Neural Information Processing Systems(pp. 937–945).
Jain, P., & Oh, S. (2014). Provable tensor factorization with missing data. In
Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence and K. Weinberger (Eds.),
Advances in neural information processing systems 27, 1431–1439. Curran
Associates, Inc.
Kanagawa, H., Suzuki, T., Kobayashi, H., Shimizu, N., & Tagami, Y. (2016).
Gaussian process nonparametric tensor estimator and its minimax optimality.
Proceedings of the 33rd International Conference on Machine Learning(ICML2016) (pp. 1632–1641).
Koltchinskii, V. (2012). Sharp oracle inequalities in low rank estimation (Technical
Report). arXiv:1210.1144.
Kuelbs, J., & Li, W. V. (1993). Metric entropy and the small ball problem forgaussian measures. Journal of Functional Analysis, 116, 133–157.
Leurgans, S., Ross, R., & Abel, R. (1993). A decomposition for three-way arrays.
SIAM Journal on Matrix Analysis and Applications, 14, 1064–1083.
Li, W. V., & Shao, Q.-M. (2001). Gaussian processes: inequalities, small ballprobabilities and applications. Stochastic Processes: Theory and Methods, 19,
533–597.
59 / 59

Mu, C., Huang, B., Wright, J., & Goldfarb, D. (2014). Square deal: Lowerbounds and improved relaxations for tensor recovery. Proceedings of the 31th
International Conference on Machine Learning (pp. 73–81).
Negahban, S., Ravikumar, P., Wainwright, M. J., & Yu, B. (2012). A unified
framework for high-dimensional analysis of m-estimators with decomposable
regularizers. Statistical Science, 27, 538–557.
Negahban, S., & Wainwright, M. J. (2012). Restricted strong convexity andweighted matrix completion: Optimal bounds with noise. Journal of Machine
Learning Research, 13, 1665–1697.
Oseledets, I. V. (2011). Tensor-train decomposition. SIAM Journal on Scientific
Computing, 33, 2295–2317.
Phien, H. N., Tuan, H. D., Bengua, J. A., & Do, M. N. (2016). Eﬃcient tensor
completion: Low-rank tensor train. arXiv preprint arXiv:1601.01083.
Rohde, A., & Tsybakov, A. B. (2011). Estimation of high-dimensional low-rankmatrices. The Annals of Statistics, 39, 887–930.
Shah, P., Rao, N., & Tang, G. (2015). Optimal low-rank tensor recovery fromseparable measurements: Four contractions suﬃce. arXiv preprint
arXiv:1505.04085.
Signoretto, M., Lathauwer, L. D., & Suykens, J. A. K. (2013). Learning tensors in
reproducing kernel Hilbert spaces with multilinear spectral penalties. CoRR,
abs/1310.4977.
59 / 59

Srebro, N., Alon, N., & Jaakkola, T. (2005). Generalization error bounds forcollaborative prediction with low-rank matrices. Advances in Neural Information
Processing Systems (NIPS) 17.
Sun, R., & Luo, Z.-Q. (2014). Guaranteed matrix completion via non-convex
factorization. arXiv preprint arXiv:1411.8003.
Sun, W., Wang, Z., Liu, H., & Cheng, G. (2015). Non-convex statisticaloptimization for sparse tensor graphical model. Advances in Neural Information
Processing Systems (pp. 1081–1089).
Suzuki, T. (2015). Convergence rate of Bayesian tensor estimator and its minimax
optimality. Proceedings of the 32nd International Conference on Machine
Learning (ICML2015) (pp. 1273–1282).
Suzuki, T., Kanagawa, H., Kobayashi, H., Shimizu, N., & Tagami, Y. (2016).
Minimax optimal alternating minimization for kernel nonparametric tensor
learning. The 30th Annual Conference on Neural Information Processing
Systems (NIPS2016) (p. accepted).
Tomioka, R., & Suzuki, T. (2013). Convex tensor decomposition via structured
schatten norm regularization. Advances in Neural Information Processing
Systems 26 (pp. 1331–1339). NIPS2013.
Tomioka, R., Suzuki, T., Hayashi, K., & Kashima, H. (2011). Statistical
performance of convex tensor decomposition. Advances in Neural Information
Processing Systems 24 (pp. 972–980). NIPS2011.
59 / 59

Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis.
Psychometrika, 31, 279–311.
Wimalawarne, K., Sugiyama, M., & Tomioka, R. (2014). Multitask learning meets
tensor factorization: task imputation via convex optimization. In
Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence and K. Weinberger (Eds.),
Advances in neural information processing systems 27, 2825–2833. Curran
Associates, Inc.
Yokota, T., Zdunek, R., Cichocki, A., & Yamashita, Y. (2015a). Smoothnonnegative matrix and tensor factorizations for robust multi-way data analysis.
Signal Processing, 113, 234–249.
Yokota, T., Zhao, Q., & Cichocki, A. (2015b). Smooth parafac decomposition fortensor completion. arXiv preprint arXiv:1505.06611.
Zdunek, R. (2012). Approximation of feature vectors in nonnegative matrix
factorization with gaussian radial basis functions. International Conference on
Neural Information Processing (pp. 616–623).
Zhang, Z., & Aeron, S. (2015). Exact tensor completion using t-svd. arXiv
preprint arXiv:1502.04689.
Zhao, T., Wang, Z., & Liu, H. (2015). A nonconvex optimization framework for
low rank matrix estimation. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama,
R. Garnett and R. Garnett (Eds.), Advances in neural information processing
systems 28, 559–567. Curran Associates, Inc.
59 / 59

Zheng, Q., & Tomioka, R. (2015). Interpolating convex and non-convex tensor
decompositions via the subspace norm. Advances in Neural Information
Processing Systems (pp. 3088–3095).

59 / 59