Julia言語重み付き非負行列因子分解

重み付き非負行列分解というタスクがある．欠損ありのNMFとかに使える．研究に必要だったので，オーソドックスな次の手法を実装した．

www.semanticscholar.org

実装は簡単で普通のNMFの更新式をちょっと変えるだけ．フロベニウスノルム尺度のものとKL尺度のもの．Wが重み行列．W=ones(n,m) で通常のNMFになる．

using LinearAlgebra
using Random
using Arpack

function D_KL(A, B)
    n, m = size(A)
    kl = 0.0
    for i = 1:n
        for j = 1:m
            kl += A[i,j] * log( A[i,j] / B[i,j] ) - A[i,j] + B[i,j]
        end
    end
    return kl
end

function wnmf_euc(A, W, r ; max_iter=200, tol = 1.0E-4, verbose = false)
    m, n = size(A)
    U = rand(m, r)
    V = rand(r, n)

    error_at_init = norm(A - U*V)
    previous_error = error_at_init
    for iter = 1:max_iter
        V .=  V .* ( U' * (W .* A) ) ./ ( U' * (W .* (U * V) ) )
        U .=  U .* ( (W .* A)  * V') ./ ( (W .* (U  * V)) * V' )

        if tol > 0 && iter % 1 == 0
            error = norm(A - U*V)
            if verbose
                println("iter: $iter cost: $error")
            end
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return U, V
end

function wnmf_kl(A, W, r ; max_iter=200, tol = 1.0E-3, verbose = false)
    m, n = size(A)
    U = rand(m, r)
    V = rand(r, n)
    one_mn = ones(m, n)

    error_at_init = D_KL(A, U*V)
    previous_error = error_at_init
    for iter in 1:max_iter
        V .= (V ./ (U'*W)) .* (U' * (( W .* A ) ./ (U*V)))
        U .= (U ./ (W*V')) .* ((( W .* A ) ./ (U*V)) * V' )

        if tol > 0 && iter % 10 == 0
            error = D_KL(A, U*V)
            if verbose
                println("iter: $iter cost: $error")
            end
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return U, V
end