Julia の TensorToolbox を使うと気軽に CP 分解ができます.
例えば,ランダムな4×4×4×4テンソルX
をCP分解するなら
using TensorToolbox X = rand(4,4,4,4) F = cp_als(X, 3)
でOKです.ここで,ターゲットランクは3にしました.これで,F.lambda
に固有値もどき(テンソル分解なので行列の意味での固有値ではない)が格納され,F.fmat[1]
,F.fmat[2]
,F.fmat[3]
,F.fmat[4]
には因子が格納されます.これらの因子から元のテンソルを復元すると,(ターゲットランクが十分に大きければ)それは入力テンソルX
の近似になっているはずです.
ということで,F.mat
とF.lambda
から元のテンソルを再構成するコードをかきました.もっとよい方法があったら教えて下さい.
using LinearAlgebra using TensorToolbox function reconst(F) D = ndims(F) R = length(F.lambda) F1 = F.lambda[1] * ttt(F.fmat[1][:,1], F.fmat[2][:,1]) for d = 3:D F1 = ttt(F1, F.fmat[d][:,1]) end for r = 2:R Fr = F.lambda[r] * ttt(F.fmat[1][:,r], F.fmat[2][:,r]) for d = 3:D Fr = ttt(Fr, F.fmat[d][:,r]) end F1 .= F1 .+ Fr end return F1 end
これで reconst(F) は X を近似するはず.
julia> X = rand(3,3,3,3); julia> for r=1:5:30 @show r,norm( reconst(cp_als(X,r)) - X ) / norm(X) end (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (1, 0.4363688254611265) (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (6, 0.29372126298825085) (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (11, 0.19977496289421565) (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (16, 0.06927809613284029) (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (21, 0.022644681361821475) (r, norm(reconst(cp_als(X, r)) - X) / norm(X)) = (26, 7.873258786342794e-6)
確かにターゲットランク r
を大きくするとどんどん再構成誤差(を入力テンソルの総和で割った値)が小さくないって行く様子がわかる.一応,テストとして,
mrank( reconst( cp_als(X,r) ) ) == (r,r,r,r)
になっている.
しかし,CP分解って,こんなにターゲットランクが大きくないと再構成うまくいかないのか???