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

αダイバージェンスやβダイバージェンスでのNMFを一般化するために，αβダイバージェンスというものを定義します．入力行列と出力低ランク行列間のαβダイバージェンスが小さくなるようなNMFアルゴリズムをJuliaで実装しました．αβダイバージェンスやアルゴリズムは以下の論文の通りに実装しました．

Entropy | Free Full-Text | Generalized Alpha-Beta Divergences and Their Application to Robust Nonnegative Matrix Factorization

あんまりちゃんとテストしてないので，使うときは気を付けてね． とりあえず，αβダイバージェンスD_abを定義に従って計算できるようにD_ab関数を用意しておく．α=β=1で通常のフロベニウスノルムになる．

using LinearAlgebra

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

function D_IS(P, Q)
    #Itakura-Saito distance

    n, m = size(P)
    d_is = 0.0
    for i = 1:n
        for j = 1:m
            d_is += log( Q[i,j] / P[i,j] ) + P[i,j] / Q[i,j] + - 1.0
        end
    end
    return d_is
end

function D_logeuc(P, Q)
    n, m = size(P)
    log_P = log.(P)
    log_Q = log.(Q)
    return sum( (log_P - log_Q).^2 )
end

function D_ab(P, Q, alpha, beta)
    # defined by Andrzej Cichocki in
    # Entropy2011,13, 134-170; doi:10.3390/e13010134OPEN
    # See also equation (23) in
    # core.ac.uk/download/pdf/51394763.pdf
    # D_ab = D_euc when alpha = 1 and beta = 1

    if alpha != 0 && beta == 0
        return 1.0 / alpha^2 * D_KL(P.^alpha, Q.^beta)
    end
    if alpha + beta == 0 && alpha !=0 && beta !=0
        return 1.0 / alpha^2 * D_IS(P.^alpha, Q.^alpha)
    end
    if alpha == 0 && beta != 0
        return 1.0 / beta^2 * D_KL(Q.^beta, P.^beta)
    end
    if alpha == 0 && beta == 0
        return 1.0 / 2.0 * D_logeuc(P, Q)
    end

    n, m = size(P)
    d_ab = 0.0
    for i = 1:n
        for t = 1:m
            d_ab += P[i,t]^alpha * Q[i,t]^beta
            d_ab += -alpha / (alpha + beta) * P[i,t]^(alpha+beta)
            d_ab +=  -beta / (alpha + beta) * Q[i,t]^(alpha+beta)
        end
    end

    return -1.0 * d_ab / (alpha*beta)

end

あとは論文の更新ルールをそのまま実装する．

function nmf(P, r ; alpha=1, beta=1, max_iter=200, tol = 1.0E-3, verbose = false)
    # proposed by Andrzej Cichocki in
    # Entropy2011,13, 134-170; doi:10.3390/e13010134OPEN
    # See also equation (69) and (70) in
    # core.ac.uk/download/pdf/51394763.pdf
    #
    # (alpha, beta) = (1.0, 1.0) is usual euc-NMF
    # (alpha, beta) = (1.0, 0.0) is KL-NMF
    # (alpha, beta) = (0.0, 1.0) is IS-NMF

    n, m = size(P)
    A = rand(n, r)
    X = rand(r, m)
    one_mn = ones(n, m)

    error_at_init = D_ab(X, A*X, alpha, beta)
    previous_error = error_at_init
    for iter in 1:max_iter
        # update rule eq(69)
        Q = A*X
        Z = P.^alpha .* Q.^( beta-1 )
        X .= X .* ( (A'*Z) ./ (A'*Q.^(alpha+beta-1))).^(1.0/alpha)


        # update rule eq(70)
        Q = A*X
        Z = P.^alpha .* Q.^( beta-1 )
        A .= A .* ( (Z*X') ./ (Q.^(alpha+beta-1) * X')).^(1.0/alpha)


        if tol > 0 && iter % 10 == 0
            error = D_ab(P, A*X, alpha, beta)
            if verbose
                println("iter: $iter cost: $error")
            end
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return A, X
end

A*Xが入力行列Pを近似するrランク行列となる．

追記

私の記憶が正しければKL-NMFは行和列和を保持する．しかし，このアルゴリズムでα=1,β=0にしても行和列和が保持されない．何か間違えていると思う．そもそも論文の式(71)なんかおかしい気がするのは気のせいだろうか．

追記2

以前アップロードしたNMFの記事のコードと，上のコードとで同じ行列を入れても別の答えが返ってくる．やはりこの実装あやしい．．．（さすがに論文が間違っている可能性は低いと思うので，私の実装がおかしいのだろう．．．）

genkaiphd.hatenablog.com

追記3

α=1, β=0で更新式が以前の記事で引用した文献にある式にならない．Normalize Factor分ずれている気がする．めんどくさいが，今度，自分で更新式を導いて，どこが怪しいか調べる．