Julia言語 非負テンソルをKL情報量でCP分解．CP-APRによる低ランク近似の実装

研究に必要だったので，Julia言語でこの論文のアルゴリズム3を実装した．「KL情報量でCP分解しましょう」って話だけど，中身はポアソン分布とかいろいろ書いてあってちゃんと読んでないから分からん．とりあえず非負のテンソルCP分解の手法の一つ．

arxiv.org

Matlabではcp_arpコマンドですぐに呼び出せるらしい．

評価用に２つのテンソル間のKLダイバージェンスを求める関数を用意しておく．この関数はCP-APR本体では使わない．0 log 0 = 0 としておこう．

function DKL(P, Q)
    tensor_shape = size(P)
    dkl = 0.0
    for idx in CartesianIndices(tensor_shape)
        if P[idx] < 1.0e-8
            continue
        end
        if Q[idx] < 1.0e-8
            continue
        end
        dkl += P[idx] * log( P[idx] / Q[idx] )
    end
    return dkl - sum(P) + sum(Q)
end

次に初期テンソルMを決める関数と，Sを更新する関数を用意しておく．Mの決め方は6.1.節を参照にしたけど，まぁ謎．

function get_M(I, N, R)
    A = []
    for i = 1 : N
        a = rand(I[i], R)
        for j = 1:R
            coin = rand()
            if coin < (1.0 / R)
                a[j] = 100*rand()
            else
                a[j] = rand()
            end
        end

        for r = 1 : R
            a[:,r] = normalize(a[:,r],1)
        end
        push!(A, a)
    end
    lamb = rand(R)
    return lamb, A
end

function getS(k, In, R, An, Phin, kap, kaptol)
    S = zeros(In, R)
    if k == 1
        return S
    else
        for i = 1 : In
            for r = 1 : R
                if An[i,r] < kaptol && Phin[i,r] > 1.0
                    S[i, r] = kap
                else
                    S[i, r] = 0.0
                end
            end
        end
        return S
    end
end

次に本体をスクラッチ．TensorToolboxとか，TensorDecompositionsとか便利ですね．特にテンソルのKhatri–Rao productがTensorDecompositionsにあるのは便利でした．(LinearAlgebraのkronを繰り返せば同じことできるけど．．．）

入力テンソルXのCPランクをRにする関数．RはIntね．分解するだけが目的なら#reconstract以降は不要．

using LinearAlgebra
using TensorDecompositions
using TensorToolbox
using InvertedIndices

function CPARP(X , R ; kmax = 600, lmax = 10, tau = 1.0e-1, kap = 0.1, kaptol = 1.0e-3, epsilon = 1.0E-10 )
    I = size(X)
    N = ndims(X)

    lamb, A = get_M(I, N, R)
    #A[i] \in R^{I[i] \times R}

    Phi = []
    for n = 1 : N
        push!(Phi, rand(I[n], R))
    end

    S = undef
    for k = 1 : kmax

        isConverged = true
        for n = 1 : N
            S = getS(k, I[n], R, A[n], Phi[n], kap, kaptol)
            B = (A[n]+S) * diagm( vec(lamb) )
            PI = TensorDecompositions.khatrirao( reverse(A[Not(n)]) )'

            for l = 1:lmax
                Phi[n] = ( tenmat(X,n) ./ max.(B*PI,epsilon) ) * PI'
                if sum(abs.( min.(B, ones(size(Phi[n])) - Phi[n]) )) < tau
                    break
                end
                isConverged = false
                B .= B .* Phi[n]
            end
            lamb = ones(1, size(B)[1]) * B
            A[n] = B * diagm( vec( 1 ./ lamb )  )
        end
        if isConverged == true
            break
        end
        if k == kmax
            println("NOT CONVERGE k=$k")
        end
    end

    #reconstract
    G = zeros( fill(R, N)... )
    for i = 1:R
        idx = fill(i, N)
        G[idx...] = lamb[i]
    end

    Y = ttm(G, A[1], 1)
    for n = 2:N
        Y = ttm(Y, A[n], n)
    end
    return Y
end

試しに20×20×20のテンソルを低CPランクで近似してみる．

X = rand(20,20,20)
Rs = [20, 17, 13, 10, 8, 5, 4, 3, 2, 1]
for R in Rs
    Y = CPARP(X, R)
    println("R=$R ", "\t LS=", norm(X-Y), "\tKL=", DKL(X, Y) )
end

順調にLSエラーとKLエラーが上がっていくのが分かる．（CPランクが低いテンソルの表現力が低いという意味）

R=20      LS=22.880251175094255  KL=606.339285365445
R=17     LS=23.234143014881077  KL=624.7517055938652
R=13     LS=23.750102851245273  KL=653.4094801796946
R=10     LS=24.248429669696232  KL=679.6102072187259
R=8      LS=24.540291085977568  KL=695.6061065373501
R=5      LS=25.012954092190817  KL=721.7526514489091
R=4      LS=25.194662960281825  KL=731.3447268180694
R=3      LS=25.3753714333728    KL=740.7357057303489
R=2      LS=25.53270583385743   KL=748.5909542140225
R=1      LS=25.69656543115171   KL=756.886141185732

入力テンソルを大きくしたり，Rを大きくしたりすると，収束しなくなったりするので，うまい具合にハイパラをチューニングする必要がある． tau と kap を調節するとうまくいくことが多い．論文じゃいろいろ考察があるけどlmaxはどう設定したらよいかよくわからん．