TCA
テンソル分解は様々な種類があるが、そのうち比較的単純とされているCP分解を紹介する。 3階テンソルは、 \[\begin{equation} X_{i,j,k} \end{equation}\] として添字が3つある行列を奥行方向に並べたようなデータを意味する。添字を \(N\) 個にして、N階テンソルを考えることもある (そちらのほうが一般的)。
CP分解もまた、特別な制約なしには直交しない成分を抽出する。そして、データテンソル \(\boldsymbol{X} \in \mathbb{R}^{S \times T \times K}\) を、ベクトル \(\boldsymbol{s}_r \in \mathbb{R}^{S \times 1}\)、\(\boldsymbol{t}_r \in \mathbb{R}^{T \times 1}\)、\(\boldsymbol{u}_r = (u_{1r}, u_{2r}, ..., u_{Kr}) \in \mathbb{R}^{K \times 1}\)を利用して、 \[\begin{equation} \boldsymbol{X}_{:,:,k} = \sum_{r=1}^R \lambda_r \boldsymbol{s}_r\boldsymbol{t}_r^T u_{kr} \end{equation}\] として分解することである。ここで、\(R\) は、いくつのベクトル要素に分解するかを表す指標である。ここでのポイントは、\(\boldsymbol{s}_r\boldsymbol{t}_r^T\)の部分は ３要因目 \(k\) に依存しないことである。つまり、行列分解と同様に、例えば関節角度・筋電ならば空間情報、時間情報に分けて表現する。この空間情報、時間情報は３要因目 \(k\) に依存せず、すべてに共通しているという想定を置く。そして、\(\boldsymbol{s}_r\boldsymbol{t}_r^T\)が \(k\)番目の要因で用いられた大きさを \(u_{kr}\) で表現している。
例えば、\(\boldsymbol{X}_{:,:,k}\) が \(k\)番目の速度にて歩行/走行していたときの関節角度データの場合、\(\boldsymbol{s}_r\) は関節のグループ(=空間情報)、\(\boldsymbol{t}_r\) はその時間変動 (=時間変動) となり、これらはすべての \(k\) で共通して用いる。そして、これらの時空間成分がどれだけ\(k\)番目の速度にて関わっていたかを \(u_{kr}\) で推定しようというものである。
実際にCP分解を実行してみよう。
library(rTensor)
X = cbind(matrix(runif(100, 0, 1), 5, 20), matrix(runif(100, 1, 2), 5, 20), matrix(runif(100, 2, 3), 5, 20))
# 正規化 (各空間情報の最小値が0、最大値が1になるようにする)
dimX = dim(X)
mX = apply(X, 1, mean)
sX = apply(X, 1, sd)
X = (X - mX%*%matrix(1, 1, dimX[2]))/(sX%*%matrix(1, 1, dimX[2]))
X = array(X, c(5, 20, 3))
dat = X
X = as.tensor(dat)
R = 3
P = cp(X, num_components = R)##
|
| | 0%
|
|=== | 4%
|
|====== | 8%
|
|======== | 12%
|
|=========== | 16%
|
|============== | 20%
|
|================= | 24%
|
|==================== | 28%
|
|====================== | 32%
|
|========================= | 36%
|
|============================ | 40%
|
|=============================== | 44%
|
|================================== | 48%
|
|==================================== | 52%
|
|======================================= | 56%
|
|========================================== | 60%
|
|============================================= | 64%
|
|================================================ | 68%
|
|================================================== | 72%
|
|===================================================== | 76%
|
|======================================================== | 80%
|
|=========================================================== | 84%
|
|============================================================== | 88%
|
|================================================================ | 92%
|
|=================================================================== | 96%
|
|======================================================================| 100%Spatial_modules = as.matrix(unlist(P$U[[1]]))
Temporal_modules = as.matrix(unlist(P$U[[2]]))
Task_components = as.matrix(unlist(P$U[[3]]))
lambda = P$lambdas
par(mfrow=c(3,3))
# 以下はnon-normalized
for(dam in 1:3){
barplot(Spatial_modules[,dam])
plot(Temporal_modules[,dam], type = "l")
plot(Task_components[,dam])
}
(結構重要な) 余談: 様々 \(R\) を決める指標は提案されているが、多くは元データを再現する割合を元に議論される。正直、様々提案されているが、これといって \(R\) を決める方法はない。例えばAICを使う方法もあるが、例えば\(\lambda_r = 2\)のとき、\(\boldsymbol{s}_r\)を2倍するのか \(\boldsymbol{t}_r\) を2倍するのかどちらでも問題ない。本来AICが想定する状況は、冗長性のない状況であり、NNMFにAICを適用することは妥当ではない。また、Cross-validationという方法もあるが、基本は\(R\)は大きいほどいい。とはいえ、大きすぎる \(R\) はノイズも拾い、解釈も難しい。どうしようもないのが現実だが、元データを80%再現、90%再現あたりが無難な閾値として用いられる。