Julia 言語 TensorToolbox で CP 分解して，再構成する．

Julia の TensorToolbox を使うと気軽に CP 分解ができます．

github.com

例えば，ランダムな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分解って，こんなにターゲットランクが大きくないと再構成うまくいかないのか？？？