Julia言語 Non-negative Multiple Matrix Factorization (NMMF)の実装

複合データ分析技術についてこのpdfで勉強してたら，NMMFというものを知った．

ちゃんとした議論は，IJCAIでしている． www.ijcai.org

人工知能学会の予稿や論文，絵本のパターン抽出や，和文論文が参考になる．

私も研究で使うことになったのでささっと実装した．

とりあえず，いつも使っているコスト関数

using LinearAlgebra

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

function cost_function(cost, X, Xr)
    if cost == "LS"
        return norm(X-Xr)

    elseif cost == "KL"
        dkl = get_dkl(X, Xr)
        return dkl
    end
end

としておく． 本体は以下．これで実行すれば，ロスが下がっていく様子がわかる． 初期化は超適当に連続一様分布から行っている．

"""
Non-Negative Multiple Matrix Factorization
proposed by Koh Takeuchi, Katsuhiko Ishiguro, Akisato Kimura and Hiroshi Sawada
See [original paper](https://www.ijcai.org/Proceedings/13/Papers/254.pdf)
You can see the implementation in [this Japanese documents]
(http://www.kecl.ntt.co.jp/icl/ls/members/tatsushi/PDF/IEICE_vol99_no6_543-550.pdf)

# Aruguments
- `X` input I \times J matrix
- `Y` input I \times K matrix
- `R` target rank
"""
function NMMF_KL(X, Y, R; max_iter=200, tol = 1.0E-3, verbose = false)
    (I, J) = size(X)
    (I, K) = size(Y)
    A = rand(I, R)
    B = rand(J, R)
    C = rand(K, R)

    Xhat = A * B'
    Yhat = A * C'

    cost_at_init  = cost_function("KL", X, Xhat) + cost_function("KL", Y, Yhat)
    previous_cost = cost_at_init
    for iter = 1 : max_iter
        A .= A .* ( ( (X ./ Xhat ) * B + (Y ./ Yhat) * C) ./ ( ones(I,J) * B + ones(I,K) * C ) )

        Xhat = A * B'
        Yhat = A * C'

        B .= B .* (( (X ./ Xhat)' * A ) ./ (ones(J,I)*A))
        Xhat = A * B'

        C .= C .* (( (Y ./ Yhat)' * A ) ./ (ones(K,I)*A))
        Yhat = A * C'

        if tol > 0 && (iter % 10 == 0)
            cost = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
            if verbose
                println("step:$iter loss:$cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end
    return A, B, C
end

function NMMF_LS(X, Y, R; max_iter=200, tol = 1.0E-3, verbose = false)
    (I, J) = size(X)
    (I, K) = size(Y)
    A = rand(I, R)
    B = rand(J, R)
    C = rand(K, R)

    Xhat = A * B'
    Yhat = A * C'

    cost_at_init  = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
    previous_cost = cost_at_init
    for iter = 1 : max_iter
        A .= A .* ( (X * B + Y * C) ./ ( Xhat * B + Yhat * C ))
        Xhat = A * B'
        Yhat = A * C'

        B' .= B' .* ( (A' * X) ./ ( A' * Xhat ) )
        Xhat = A * B'

        C' .= C' .* ( (A' * Y) ./ ( A' * Yhat ) )
        Yhat = A * C'

        if tol > 0 && (iter % 10 == 0)
            cost = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
            if verbose
                println("step:$iter loss:$cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end
    return A, B, C
end

X = rand(20,10)
Y = rand(20,8)
A,B,C = NMMF_LS(X, Y, 2, verbose=true)
A,B,C = NMMF_KL(X, Y, 4, verbose=true)