ひょっこりTucker Decomposition in R

急襲に遭い、トイレに籠城していたお供として村上春樹『職業としての小説家』を読み直していた。
この記事の主題とは全く関係ないが、なんとなく印象に残った箇所を書き残しておく。

アイザック・ディネーセンは「私は希望もなく、絶望もなく、毎日ちょっとずつ書きます」と言っています。
それと同じように、僕は毎日十枚の原稿を書きます。とても淡々と。
「希望もなく、絶望もなく」というのは実に言い得て妙です。
朝早く起きてコーヒーを温め、四時間か五時間、机に向かいます。
一日十枚原稿を書けば、一か月で三百枚書けます。単純計算すれば、半年で千八百枚が書けることになります。
（p.155）

「海辺のカフカ」の第一稿が1,800枚で、プロ野球の開幕と同時に書き始めて、日本シリーズが始まるころに書き終えたという。
毎日XX枚書けばYY日でZZ枚かけるぞ、っていうのは小学生でも考えそうな皮算用だが、それを淡々と積み重ねて現実にするのが村上のしなやかな強さである。
希望も絶望もない、という心的状態が大事なのは、自分も歳を重ねて何となく理解できる気がしてきている。平熱の思考の重要性というか。

Introduction

なんとなく使うかもしれないと多少はインプットしていたが、
忘れて思い出してを繰り返し続けるので、使って覚えるシリーズ。

文系がどうとか理系がどうとかという類型論の不毛さは重々承知なのだが、それでもやっぱり文系には行列計算はきついのである。
2次元の行列で悲鳴を上げているのに、N(3+)次元のテンソルと来た日には眩暈と吐気が同時襲来せざるをえない。
ひとまず使いながら身体に覚えさせるしかないと思った次第である。

参考書籍＆URL。

◆手法についての解説

石黒・林（2016）「関係データ学習」第5-6章
Alternating Least Square (ALS) でCP分解 - でかいチーズをベーグルする
Tucker分解の導出と実装 - でかいチーズをベーグルする

関係データ学習 (機械学習プロフェッショナルシリーズ)

作者: 石黒勝彦,林浩平
出版社/メーカー: 講談社
発売日: 2016/12/07
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

◆実際の応用の仕方とかまでカバーしたやつ（Rコードあり）
Understanding the Tucker decomposition, and compressing tensor-valued data (with R code) – Alexej Gossmann
→ Tucker decomposition推し。とてもわかりやすくてよい。

◆使用するpackage（rTensor）についての参考資料

定義と基本的な操作

定義

as.tensorに配列（array）を引数として与えて作成する。
たとえば2×3×2のテンソルを作成するとき、以下のようにする。
テンソルの次元は、arrayを作成するときに既に与えておく（as.tensorの引数では次元を指定できない）

tnsr1 <- as.tensor( array(1:12 , dim= c(2,3,2) ) 
tnsr1@data
# , , 1
# 
# [,1] [,2] [,3]
# [1,]    1    3    5
# [2,]    2    4    6
# 
# , , 2
# 
# [,1] [,2] [,3]
# [1,]    7    9   11
# [2,]    8   10   12

基本演算

内積とフロベニウスノルムはinnerProd()とfnorm()で。

tnsr_r <- rand_tensor(c(2,3,2) ) #標準正規分布からランダムなtensorを作成
innerProd(tnsr1 , tnsr_r) #内積をとる

fnorm( tnsr1) #フロベニウスノルム
sqrt( innerProd( tnsr1, tnsr1)) #これとやってることは同じ

モードlのi番目のスライスの取得は、[テンソル]@dataを用いて、Arrayオブジェクトと同じようなindexによるアクセスでいける。

tnsr1@data[,3,] #モード2の3番目のスライス
# [,1] [,2]
# [1,]    5   11
# [2,]    6   12
tnsr1@data[2,,] #モード1の2番目のスライス
# [,1] [,2]
# [1,]    2    8
# [2,]    4   10
# [3,]    6   12

モード積

モード積による次元の縮約や特定のモードの次元の伸縮（石黒・林2016: 122-124参照）は、ttm()を使う。

一応モード積とは何ぞや、についてざっくりと説明する。
まず行列について考えるときに、I× Jの行列AとJ次元ベクトルvの積Avは、
$A = \begin{pmatrix} a_1 & a_2 & a_3 \end{pmatrix}$ とAをI次元のベクトル $a_j$ の集まりと考えると、
$Av = \Sigma_j^J a_j v_j$ である。つまり、各列に該当するベクトル（要素は行数Iだけある）の線形結合を得ていると考えられる

I×Jの行列AとJ×Kの行列Bの積ABも同じように考えることができて、 $a_j$ の線形結合であらわされる次元IのベクトルがK個横に並んでいると見なせる。

これをテンソルの場合に拡張する。簡単のため3次元のテンソルX（次元はI × J × K）について話す。
行列の場合は、つぶしたり伸縮したりする対象の次元（＝列）をindexとして得た別次元（＝行）の要素を、線形結合していた。
基本のアイディアは同じで、モードl積を取る場合は、次元lをindexとして1からLまで動かしていったときに得られるスライス $X_i^{\{l\}}$ を線形結合すればよいだけの話である。

だから、XとJ次元ベクトルvのモード2積を考えると
$X \times_2 v = \Sigma_{j=1}^J v_j X_j^{\{J\}}$
となる。
行列とのモード積についても同様に考えることができて、XとQ×J行列Yとのモード2積（次元はI × Q × Kとなる）のq番目のsliceは以下の式で得られる。
$X \times_2 Y = \Sigma_{j=1}^J Q_{qj} X_j^{\{J\}}$

なんというか図らずもテンソルについて勉強しなおしたことで行列計算への理解も上がった気がする。
めでたしめでたし。

rTensor::ttmの話に戻ると、第一引数にテンソル、第二引数に掛け合わせる行列、第三引数にモード積の対象となるモードを指定する
注意点として、モード積の対象としたいモードの次元がJのとき、第二引数の列数をJにせねばならない、ということだ。
行列どうしの積ABの場合、Bの行数をAの列数に合わせる必要があるので、ちょっと直感的には戸惑う。
あと細かい注意点として、ベクトルを第二引数に指定するときも1×Jのmatrix型にしなければならない

#先ほど作成したテンソルをベクトルとのモード2積をとる。
vec <- matrix( 1:3 , nrow=1)
tnsr1_prd_vec  <- ttm( tnsr1  , mat= vec ,m=2) 
tnsr1_prd_vec@data
# , , 1
# 
# [,1]
# [1,]   22
# [2,]   28
# 
# , , 2
# 
# [,1]
# [1,]   58
# [2,]   64 
# モード2が潰れて2×2のテンソル（行列）が得られる

#先ほど作成したテンソルを2×3の行列とのモード2積をとることで、2×2×2のテンソルにする（モード2を縮める）
mat1 <- matrix( 1:6 , byrow=T , ncol=3)
tnsr1_prd_mat <- ttm( tnsr1 , mat = mat1 , m = 2)
tnsr1_prd_mat@data
# , , 1
# 
# [,1] [,2]
# [1,]   22   49
# [2,]   28   64
# 
# , , 2
# 
# [,1] [,2]
# [1,]   58  139
# [2,]   64  154
#モード2が次元2に縮んで、2×2×2のテンソルが得られる。

ちなみに加減乗除は普通にas.tensor()で得られたオブジェクトどうしを+, -, *, / あたりで足したりかけたりすればいい。直感的。

Tucker Decompositionしてみる

Tucker分解自体の説明は、参考URLや参考書籍、その他もろもろの世間に出回っている情報を参照
CP分解と二大メソッドみたいになってるけど、特異値分解の素直な拡張としてとらえられるのはTuck
er分解らしい。
ちなみに、Tucker分解の特殊ケースがCP分解なので、前者は後者を包含する関係にある。

rTensorではtucker()を使うことでTucker decompositionを実行できる。

tnsr3 <- rand_tensor(c(20,15,10))
tucker_res <- tucker( tnsr3 , ranks=c(3,3,2), max_iter=500)

戻り値についてであるが、$Zでコアテンソルに、$Uで因子行列にアクセスできる。
$norm_percentで分散（というかフロベニウスノルム）の説明率に、$fnorm_residで説明しきれていないフロベニウスノルムの残余の値にアクセスできる。
ちなみにもとのテンソルからfnorm(tensor)でもともとのテンソルのフロベニウスノルムは取得できる。