Python αダイバージェンスを小さくするNMF

KLダイバージェンスを最小化するNMFはよく知られていますし，sklearnで'beta_loss=1'にすればすぐに利用できます． 今回はKLダイバージェンスの拡張であるαダイバージェンスを小さくするNMFを実装しました．

この論文の更新式をさくっと実装しました． www.sciencedirect.com

収束条件とかはsklearnの普通のNMFに揃えようとしたけど，ちょっと違うかも． （もしかしたら，sklearnの方はコスト関数のルートをとっているかもしれない）

def D_KL(P,Q):
    epsilon = 1.0e-10
    P = P+epsilon
    Q = Q+epsilon
    return np.sum(P*np.log(P/Q)) - np.sum(P) + np.sum(Q)

def alpha_KL(X,AS,alpha):
    if alpha == 1.0:
        return D_KL(X, AS)
    if alpha == 0.0:
        return D_KL(AS, X)
    
    alpha_kl = 0.0
    n, m = np.shape(X)
    term = 0.0
    for i in range(n):
        for j in range(m):
            term += alpha * X[i,j] + (1-alpha) * AS[i,j] - math.pow(X[i,j],alpha) * math.pow((AS[i,j]),1-alpha)
    
    alpha_kl = term / (alpha*(1-alpha))
    
    return alpha_kl
    
def nmf_alpha(X, r, alpha, cnt_max=100, tol=1.0e-3, verbose=False):
    n, m = np.shape(X)
    S = np.random.rand(r, m) 
    A = np.random.rand(n, r)
    
    if alpha == 0.0:
        alpha += 1.0e-5
    
    error_at_init = alpha_KL(X, np.dot(A,S),alpha)
    previous_error = error_at_init
    for cnt in range(cnt_max):
        # update S
        for i in range(r):
            for j in range(m):
                sumA = np.sum(A, axis=0)
                term = 0.0
                AS = np.dot(A,S)
                for k in range(n):
                    term += A[k,i] * math.pow(( X[k,j] / AS[k,j] ),alpha)
                S[i,j] *= pow( term / sumA[i], 1.0/alpha )

        # update A
        for i in range(n):
            for j in range(r):
                sumS = np.sum(S, axis=1)
                term = 0.0
                AS = np.dot(A,S)
                for k in range(m):
                    term += S[j,k] * math.pow(( X[i,k] / AS[i,k] ),alpha)
                A[i,j] *= pow( term / sumS[j], 1.0/alpha )
                
        if tol>0 and cnt % 10 == 0:
            error = alpha_KL(X, np.dot(A,S), alpha)
            if verbose:
                print(f"iter:{cnt} cost:{error}")
            if (previous_error - error) / error_at_init < tol:
                break
            previous_error = error
    return A, S

とりあえず，こんな感じで実行． 分解後のターゲットランクは2にして，αは1.0にしてみる．

X = np.random.rand(30,30)
A,S = nmf_alpha(X, 2, 1.0, verbose=True, cnt_max=200)

実行結果はこんな感じ．

iter:0 cost:88.0270095179497
iter:10 cost:77.71017624081105
iter:20 cost:75.81128050686453
iter:30 cost:73.81279529499801
iter:40 cost:72.67141167025983
iter:50 cost:72.32208831953886
iter:60 cost:72.20333908949885

ちゃんとロスがさがっているので正しいっぽい． α=1の極限では，普通のKL-NMFと一致するので，sklearnの結果と一致するか検証． beta_loss=1はコスト関数をKLにするということ．

from sklearn.decomposition import NMF
sklearn_NMF = NMF(n_components=2, init='random', random_state=0, max_iter=200, beta_loss=1, solver='mu')
W = sklearn_NMF.fit_transform(X)
H = sklearn_NMF.components_
R = np.dot(W,H)
cost = alpha_KL(X, R,1.0)
print(cost)

すると

72.21912349939544

まぁ，ロスが同じくらいなのであってそう．