www.pudn.com > gmm_utilities.zip > gmm_distance_KLD.m

function K = gmm_distance_KLD(g1, g2, N) 
 
if nargin == 3 
    K = gmm_KLD_montecarlo(g1, g2, N); 
else 
    K = gmm_KLD_unscented(g1, g2); 
end 
 
% 
% 
 
function K = gmm_KLD_montecarlo(g1, g2, N) 
% Monte Carlo approximation of KLD for Gaussian mixtures. 
s = gmm_samples(g1, N); 
w1 = gmm_evaluate(g1, s); 
w2 = gmm_evaluate(g2, s); 
 
%K = sum(log(w1) - log(w2)) / N; 
 
ii = find(w1 ~= 0 & w2 ~= 0); % warning: this step biases the result; we need a more precise way to compute w1 and w2 
K = sum(log(w1(ii)) - log(w2(ii))) / length(ii); 
 
% 
% 
 
function K = gmm_KLD_unscented(g1, g2) 
% Unscented gmm KLD  
% Reference: Goldberger, An Efficient Image Similarity Measure Based on  
% Approximations of KL-Divergence Between Two Gaussian Mixtures, 2003. 
[D,N] = size(g1.x);  
Ds = sqrt(D); 
 
K = 0; 
for i=1:N 
    Ps = Ds * matrix_square_root(g1.P(:,:,i), 1); 
    x = repvec(g1.x(:,i), D); 
    s = [x+Ps, x-Ps]; % unscented samples for i-th component of g1 
 
    if 0 % Golberger: int f log g 
        w = gmm_evaluate(g2, s); 
        K = K + g1.w(i)*sum(log(w));     
    else % KLD: int f (log f - log g)  
        w1 = gmm_evaluate(g1, s); 
        w2 = gmm_evaluate(g2, s); 
        K = K + g1.w(i)*sum(log(w1)-log(w2));  
    end 
end 
K = K/(2*D); 
 
% 
% 
 
function R = matrix_square_root(P, type) 
switch type 
    case 1 % svd decomposition, P = U*D*U' = R*R' (UDU form is also called modified Cholesky decomposition) 
        [U,D,V] = svd(P); 
        R = U*sqrt(D); 
        %R = reprow(sqrt(diag(D))', size(U,1)) .* U; 
    case 2 % cholesky decomposition (triangular), P = R*R' 
        R = chol(P)'; 
    case 3 % principal square root (symmetric), P = R*R 
        R = sqrtm(P); 
end