Julia言語 αダイバージェンスで非負値行列因子分解

こないだのJuliaでNMFしたコードを一般のαダイバージェンスの尺度でNMFできるようにしました． αダイバージェンスの定義はここを参考にしました． 更新式はこの論文の2節に書いてあります．なんと更新式の一部をα倍したり1/α倍したりするだけ，簡単！

using LinearAlgebra
using Random
Random.seed!(123)

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

function D_alpha(A, B, alpha)
    # D_alpha(A, B, 1) = D_KL(A, B)
    # D_alpha(A, B, 0) = D_KL(B, A)
    if alpha == 1
        return D_KL(A, B)
    end
    if alpha == 0
        return D_KL(B, A)
    end

    A = A / sum(A)
    B = B / sum(B)

    n, m = size(A)
    d_alpha = 0.0
    for i = 1:n
        for j = 1:m
            d_alpha += alpha * A[i,j] + (1-alpha) * B[i,j] - X[i,j].^alpha * B[i,j].^(1-alpha)
        end
    end
    return 1.0 / ( alpha - alpha^2 ) * d_alpha
end

function nmf_alpha(X, r, alpha; max_iter=200, tol = 0.000001)
    m, n = size(X)
    W = rand(m, r)
    H = rand(r, n)
    one_mn = ones(m, n)

    error_at_init = D_alpha(X, W*H, alpha)
    previous_error = error_at_init
    for iter in 1:max_iter
        H .=  H .*( ( W' * ( X ./ (W*H) ).^alpha) ./ ( W' * one_mn )).^(1.0/alpha)
        W .=  W .*( ( ( X ./ (W*H) ).^alpha * H') ./ ( one_mn * H' )).^(1.0/alpha)

        if tol > 0 && iter % 10 == 0
            error = D_alpha(X, W*H, alpha)
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return W, H
end

こーゆのってどうやってテストすればいいんでしょうね． sklearnにはαダイバージェンスでNMFは実装されてないみたいなので結果を比較できなかった．更新式の一部をα倍するだけなので，pythonでも比較的容易に実装できると思う． そのうちついでにJuliaでβダイバージェンスでNMFも実装しておきたいな．

そろそろ研究で必要なので非負テンソル因子分解を実装しないといけない．．．