Let us consider the problem of approximating a divergence D between two probability distributions P and P' on ？^d from two sets of independent and identically distributed samples

following P and P'.

A divergence approximator can be used for various purposes, such as two-sample testing [1,2], change detection in time-series [3], class-prior estimation under classbalance change [4], salient object detection in images [5], and event detection from movies [6] and Twitter [7]. Furthermore, an approximator of the divergence between the joint distribution and the product of marginal distributions can be used for solving a wide range of machine learning problems [8], including independence testing [9], feature selection [10,11], feature extraction [12,13], canonical dependency analysis [14], object matching [15], independent component analysis [16], clustering [17,18], and causal direction learning [19]. For this reason, accurately approximating a divergence between two probability distributions from their samples has been one of the challenging research topics in the statistics, information theory, and machine learning communities.

A naive way to approximate the divergence from P to P', denoted by D(P||P'), is to first obtain estimators

and

of the distributions P and P' separately from their samples X and X', and then compute a plug-in approximator

However, this naive two-step approach violates Vapnik’s principle [20]:

If you possess a restricted amount of information for solving some problem, try to solve the problem directly and never solve a more general problem as an intermediate step. It is possible that the available information is sufficient for a direct solution but is insufficient for solving a more general intermediate problem.

More specifically, if we know the distributions P and P', we can immediately know their divergence D(P||P'). However, knowing the divergence D(P||P') does not necessarily imply knowing the distributions P and P', because different pairs of distributions can yield the same divergence value. Thus, estimating the distributions P and P' is more general than estimating the divergence D(P||P'). Following Vapnik’s principle, direct divergence approximators

that do not involve the estimation of distributions P and P' have been developed recently [21-25].

The purpose of this article is to give an overview of the development of such direct divergence approximators. In Section II, we review the definitions of the Kullback- Leibler divergence, the Pearson divergence, the relative Pearson divergence, and the L²-distance, and discuss their pros and cons. Then, in Section III, we review direct approximators of these divergences that do not involve the estimation of probability distributions. In Section IV, we show practical usage of divergence approximators in unsupervised change-detection in time-series, semi-supervised class-prior estimation under class-balance change, salient object detection in an image, and evaluation of statistical independence between random variables. Finally, we conclude in Section V.

A function d(·, ·) is called a distance if and only if the following four conditions are satisfied:

Non-negativity: ∀x, y, d(x, y) ≥ 0

Non-degeneracy: d(x, y) = 0 ⇔ x = y

Symmetry: ∀x, y, d(x, y) = d(y, x)

Triangle inequality: ∀x, y, z d(x, z) ≤ d(x, y) + d(y, z)

A divergence is a pseudo-distance that still acts like a distance, but it may violate some of the above conditions. In this section, we introduce useful divergence and distance measures between probability distributions.

The most popular divergence measure in statistics and machine learning is the KL divergence [26], defined as

where p(x) and p'(x) are probability density functions of P and P', respectively.

Advantages of the KL divergence are that it is compatible with maximum likelihood estimation, that it is invariant under input metric change, that its Riemannian geometric structure is well studied [27], and that it can be approximated accurately via direct density-ratio estimation [21,22,28]. However, it is not symmetric, it does not satisfy the triangle inequality, its approximation is computationally expensive due to the log function, and it is sensitive to outliers and numerically unstable, because of the strong non-linearity of the log function, and possible unboundedness of the density-ratio function p/p' [24,29].

The PE divergence [30] is a squared-loss variant of the KL divergence defined as

Because both the PE and KL divergences belong to the class of Ali-Silvey-Csiszar divergences (which is also known as f-divergences) [31,32], they share similar theoretical properties such as the invariance under input metric change.

The PE divergence can also be accurately approxi-mated via direct density-ratio estimation, in the same way as the KL divergence [23,28]. However, its approximator can be obtained analytically in a computationally much more efficient manner than the KL divergence, because the quadratic function the PE divergence adopts is compatible with least-squares estimation. Furthermore, the PE divergence tends to be more robust against outliers than the KL divergence [33]. However, other weaknesses of the KL divergence, such as asymmetry, violation of the triangle inequality, and possible unboundedness of the density-ratio function p/p', remain unsolved in the PE divergence.

To overcome the possible unboundedness of the density- ratio function p/p', the rPE divergence was recently introduced [24]. The rPE divergence is defined as

where, for 0 ≤ α < 1, q_α is defined as the α-mixture of p and p':

When α = 0, the rPE divergence is reduced to the plain PE divergence. The quantity p/q_α is called the relative density-ratio, which is always upper-bounded by 1/α for α > 0, because

Thus, it can overcome the unboundedness problem of the PE divergence, while the invariance under input metric change is still maintained.

The rPE divergence is still compatible with least-squares estimation, and it can be approximated in almost the same way as the PE divergence via direct relative density- ratio estimation [24]. Indeed, an rPE-divergence approximator can still be obtained analytically in an accurate and computationally efficient manner. However, it still violates symmetry and the triangle inequality, in the same way as the KL and PE divergence. Furthermore, the choice of α may not be straightforward, in some applications.

The L²-distance is another standard distance measure between probability distributions, defined as

The L²-distance is a proper distance measure, and thus it is symmetric and satisfies the triangle inequality. Furthermore, the density difference p(x) ？ p'(x) is always bounded as long as each density is bounded. Therefore, the L²-distance is stable without the need of tuning any control parameter such as α in the rPE divergence.

The L²-distance is also compatible with least-squares estimation, and it can be accurately and analytically approximated in a computationally efficient and numerically stable manner via direct density-difference estimation [25]. However, the L²-distance is not invariant under input metric change, which is a unique property inherent to ratio-based divergences.

In this section, we review recent advances in direct divergence approximation.

Suppose that we are given two sets of independent and identically distributed samples

from probability distributions on ？^d, with densities p(x) and p'(x), respectively:

Our goal is to approximate a divergence between from p to p' from samples X and X'.

The key idea of direct KL divergence approximation is to estimate the density ratio p/p' without estimating the densities p and p' [21]. More specifically, a density-ratio estimator is obtained by minimizing the KL divergence from p to r·p' with respect to a density-ratio model r under the constraints that the density-ratio function is non-negative and r·p' is integrated to one:

Its empirical optimization problem, where an irrelevant constant is ignored and the expectations are approximated by the sample averages, is given by

Let us consider the following Gaussian density-ratio model:

where ||·|| denotes the l₂-norm. We define the vector of parameters

as

where, ^T denotes the transpose. In this model, the Gaussian kernels are located on numerator samples

, because the density ratio p/p' tends to take large values in the regions where the numerator samples

exist. Alternatively, Gaussian kernels may be located on both numerator and denominator samples, but this seems not to further improve the accuracy [21]. When n is very large, a (random) subset of numerator samples

may be chosen as Gaussian centers, which can reduce the computational cost.

For the Gaussian density-ratio model (3), the above optimization problem is expressed as

This is a convex optimization problem, and thus the global optimal solution can be obtained easily, e.g., by gradient-projection iterations. Furthermore, the global optimal solution tends to be sparse (i.e., many parameter values become exactly zero), which can be utilized for reducing the computational cost.

The Gaussian width σ is a tuning parameter in this algorithm, and it can be systematically optimized by cross-validation with respect to the objective function. More specifically, the numerator samples X :=

are divided into T disjoint subsets

of (approximately) the same size. Then, a density-ratio estimator

is obtained using X\X_t and

(i.e., all numerator samples without X_t and all denominator samples), and its objective value for the hold-out numerator samples X_t is computed:

where |X_t| denotes the number of elements in the set X_t. This procedure is repeated for t = 1,...,T, and the σ value that maximizes the average of the above hold-out objective values is chosen as the best one.

Given a density-ratio estimator

, a KL-divergence approximator

can be constructed as

A MATLABR^® (MathWorks, Natick, MA, USA) implementation of the above KL divergence approximator (called the KL importance estimation procedure; KLIEP) is available from

“http://sugiyama-www.cs.titech.ac.jp/ ~sugi/software/KLIEP/”.

Variations of this procedure for other density-ratio models have been developed including the log-linear model [34], the Gaussian mixture model [35], and the mixture of probabilistic principal component analyzers [36]. Also, an unconstrained variant, which corresponds to approximately maximizing the Legendre-Fenchel lower bound of the KL divergence [37], was proposed [22]:

The PE divergence can also be directly approximated without estimating the densities p and p' via direct estimation of the density ratio p/p' [23]. More specifically, a density-ratio estimator is obtained by minimizing the p'- weighted squared difference between a density-ratio model r and the true density-ratio function p/p':

Its empirical criterion, where an irrelevant constant is ignored, and the expectations are approximated by the sample averages, is given by

For the Gaussian density-ratio model (3) together with the L₂-regularizer, the above optimization problem is expressed as

where λ ≥ 0 denotes the regularization parameter, and

is the n×n matrix with the (l,l')-th element defined by

and

is the n-dimensional vector with the l-th element defined by

This is a convex optimization problem, and the global optimal solution can be computed analytically as

where I denotes the identity matrix.

The Gaussian width σ and the regularization parameter λ are the tuning parameters in this algorithm, and they can be systematically optimized by cross-validation, with respect to the objective function, as follows: First, the numerator and denominator samples X =

and

are divided into T disjoint subsets

and

respectively. Then, a density-ratio estimator

is obtained using X\X_t and X'\X'_t (i.e., all samples without X_t and X'_t ), and its objective value for the holdout samples X_t and X'_t is computed:

This procedure is repeated for t = 1,...,T, and the σ and λ values that maximize the average of the above hold-out objective values are chosen as the best ones.

By expanding the squared term

in Equation (1), the PE divergence can be expressed as

Note that Equation (7) can also be obtained via Legendre- Fenchel convex duality of the divergence functional [38]. Based on these expressions, PE divergence approximators are obtained using a density-ratio estimator

as

Equation (8) is suitable for algorithmic development because this would be the simplest expression, while Equation (9) is suitable for theoretical analysis because this corresponds to the negative of the objective function in Equation (4).

A MATLAB implementation of the above method (called unconstrained least-squares importance fitting; uLSIF) is available from

“http://sugiyama-www.cs.titech.ac.jp/ ~sugi/software/uLSIF/”.

If the L₂-regularizer

in Equation (4) is replaced with the L₁-regularizer

the solution tends to be sparse [39]. Then the solution can be obtained in a computationally more efficient way [40], and furthermore, a regularization path tracking algorithm [41] is available for efficiently computing solutions with different regularization parameter values.

The rPE divergence can be directly estimated in the same way as the PE divergence [24]:

Its empirical criterion, where an irrelevant constant is ignored and the expectations are approximated by sample averages, is given by

For the Gaussian density-ratio model (3), together with the L₂-regularizer, the above optimization problem is expressed as

where

is the n×n matrix, with the (l, l')-th element defined by

This is a convex optimization problem, and the global optimal solution can be computed analytically as

Cross-validation for tuning the Gaussian width σ and the regularization parameter λ can be carried out in the same way as the PE-divergence case, with Equation (5) replaced by

By expanding the squared term

in Equation (2), the rPE divergence can be expressed as

Based on these expressions, rPE divergence approximators are given, using the relative density-ratio estimator

as

A MATLAB implementation of this algorithm (called relative uLSIF; RuLSIF) is available from

“http://sugiyama-www.cs.titech.ac.jp/ ~yamada/RuLSIF.html”.

The key idea is to directly estimate the density difference p ？ p' without estimating each density [25]. More specifically, a density-difference estimator is obtained by minimizing the squared difference between a density-difference model f and the true density-difference function p ？ p':

Its empirical criterion, where an irrelevant constant is ignored and the expectation is approximated by the sample average, is given by

Let us consider the following Gaussian density-difference model:

where

are Gaussian centers. Then the above optimization problem is expressed as

where the L₂-regularizer

is included, U is the (n + n') × (n + n') matrix, with the (l, l')-th element defined by

d denotes the dimensionality of x, and

is the (n + n')- dimensional vector, with the l-th element defined by

This is a convex optimization problem, and the global optimal solution can be computed analytically as

The above optimization problem is essentially the same form as least-squares density-ratio approximation for the PE divergence, and therefore least-squares density-difference approximation can enjoy all the computational properties of least-squares density-ratio approximation.

Cross-validation for tuning the Gaussian width σ and the regularization parameter λ can be carried out as follows: First, the numerator and denominator samples X =

and

are divided into T disjoint subsets

and

respectively. Then a density-difference estimator f_t(x) is obtained using X\X_t and X'\X'_t (i.e., all samples without X_t and X'_t ), and its objective value for the hold-out samples X_t and X'_tis computed:

Note that the first term can be computed analytically for the Gaussian density-difference model (14):

where

is the parameter vector learned from X\X_t and X'\X'_t.

For an equivalent expression of the L²-distance,

if f is replaced with a density-difference estimator

, and the expectations are approximated by empirical averages, the following L²-distance approximator can be obtained:

Similarly, for another expression

replacing f with a density-difference estimator

gives another L²-distance approximator:

Equations (15) and (16) themselves give valid approximations to L²(p,p'), but their linear combination

was shown to have a smaller bias than Equations (15) and (16).

A MATLAB implementation of the above algorithm (called least-squares density difference; LSDD) is available from

“http://sugiyama-www.cs.titech.ac.jp/ ~sugi/software/LSDD/”.

In this section, we show applications of divergence approximators in machine learning.

In real-world pattern recognition tasks, changes in class balance are often observed between the training and test phases. In such cases, naive classifier training produces significant estimation bias, because the class balance in the training dataset does not properly reflect that in the test dataset. Here, let us consider a binary pattern recognition task of classifying pattern x to class y ∈ {+1, -1}. The goal is to learn the class balance in a test dataset under a semi-supervised learning setup, where unlabeled test samples are provided, in addition to labeled training samples [42].

The class balance in the test set can be estimated by matching a π-mixture of class-wise training input densities,

to the test input density p_test(x) under some divergence measure [4]. Here, π ∈ [0, 1] is a mixing coefficient to be learned, to minimize the divergence (Fig. 1).

We use four UCI benchmark datasets (http://archive. ics.uci.edu/ml/) for numerical experiments, where we randomly choose 10 labeled training samples from each class, and 50 unlabeled test samples, following true classprior

π* = 0.1, 0.2, ... , 0.9.

The graphs on the left-hand side of Fig. 2 plot the mean and standard error of the squared difference between true and estimated class-balances π. These graphs show that LSDD tends to provide better class-balance estimates than the two-step procedure of first estimating probability densities by kernel density estimation (KDE) and then learning π.

The graphs on the right-hand side of Fig. 2 plot the test misclassification error obtained by a weighted L₂-regularized kernel least-squares classifier [43], with weighted cross-validation [44]. The results show that the LSDDbased method provides lower classification errors, which would be brought by good estimates of test class-balances.

The goal is to discover abrupt property changes behind time-series data. Let y(t) ∈ ？^m be an m-dimensional timeseries sample at time t, and let

be a subsequence of time series at time t with length k. Instead of a single point y(t), the subsequence Y(t) is treated as a sample here, because time-dependent information can be naturally incorporated by this trick [3]. Let

be a set of r retrospective subsequence samples starting at time t. Then a divergence between the probability distributions of Y(t) and Y(t + r) may be used as the plausibility of change points (Fig. 3).

In Fig. 4, we illustrate results of unsupervised change detection for the Information Processing Society of Japan Special Interest Group of Spoken Language Processing (IPSJ SIG-SLP) Corpora and Environments for Noisy Speech Recognition (CENSREC) dataset (http://research. nii.ac.jp/src/en/CENSREC-1-C.html) that records human voice in noisy environments, such as a restaurant, and the Human Activity Sensing Consortium (HASC) challenge 2011 dataset (http://hasc.jp/hc2011/), that provides human activity information collected by portable three-axis accelerometers. These graphs show that the KL-based method is excessively sensitive to noise, and thus change points are not clearly detected. On the other hand, the L²- based method more clearly indicates the existence of change points.

It was also demonstrated that divergence-based changedetection methods are useful in event detection from movies [6], and Twitter [7].

The goal is to find salient objects in an image. This can be achieved by computing a divergence between the prob-

ability distributions of image features (such as brightness, edges, and colors) in the center window and its surroundings [5]. This divergence computation is swept over the entire image, with changing scales (Fig. 5).

The object detection results on the MSRA salient object database [45] by the rPE divergence with α = 0.1 are described in Fig. 6, where pixels in gray-scale saliency maps take brighter color if the estimated divergence value

is large. The results show that visually salient objects can be successfully detected by the divergence-based

approach.

The goal is to measure how strongly two random variables U and V are statistically dependent, from paired samples

drawn independently from the joint distribution, with density p_U,V(u,v). Let us consider a divergence between the joint density p_U,V, and the product of marginal densities p_U · p_V. This actually serves as a measure of statistical independence, because U and V are independent if and only if the divergence is zero (i.e., p_U,V = p_U · p_V), and the dependence between U and V is stronger, if the divergence is larger.

Such a dependence measure can be approximated in the same way as ordinary divergences by using the two datasets formed as

The dependence measure based on the KL divergence is called mutual information (MI) [46]:

MI plays a central role in information theory [47].

On the other hand, its PE-divergence variant is called the squared-loss mutual information (SMI):

SMI is useful for solving various machine learning tasks [8], including independence testing [9], feature selection [10,11], feature extraction [12,13], canonical dependency analysis [14], object matching [15], independent component analysis [16], clustering [17,18], and causal direction estimation [19].

An L²-distance variant of the dependence measure is called quadratic mutual information (QMI) [48]:

QMI is also a useful dependence measure in practice [49].

In this article, we reviewed recent advances in direct divergence approximation. Direct divergence approximators theoretically achieve optimal convergence rates, both in parametric and non-parametric cases, and experimentally compare favorably with the naive density-estimation counterparts [21-25].

However, direct divergence approximators still suffer from the curse of dimensionality. A possible cure for this problem is to combine them with dimensionality reduction, based on the hope that two probability distributions share some commonality [51,52]. Further investigating this line would be a promising future direction.