gmm_distance_kld.m

高斯滤波器 matlab toolbox For GMMs and Gaussian kernels

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