Measuring the Discrepancy between Conditional Distributions: Methods, Properties and Applications

1 Introduction

Measuring the discrepancy or divergence between two conditional distribution functions plays a leading role in numerous real-world machine learning problems. One vivid example is the modeling of the seasonal effects on consumer preferences, in which the statistical analyst needs to distinguish the changes on the distributions of the merchandise sales conditioning on the explanatory variables such as the amount of money spent on advertising, the promotions being run, etc.

Despite substantial efforts have been made on specifying the discrepancy for unconditional distribution (density) functions (see [Anderson et al., 1994; Gretton et al., 2012b; Pardo, 2005] and the references therein), methods on quantifying the discrepancy of regression models or statistical tests

Prior art falls into two categories. The first relies heavily on the precise estimation of the underlying distribution functions using different density estimators, such as the k-Nearest Neighbor (kNN) estimator [Wang et al., 2009] and the kernel density estimator (KDE) [Lee and Park, 2006]. However, density estimation is notoriously difficult for high-dimensional data. Moreover, existing conditional tests (e.g., [Zheng, 2000; Fan et al., 2006]) are always onesample based, which means that they are designed to test if the observations are generated by a conditional distribution p(y|x) in a particular parametric family with parameter θ, rather than distinguishing p1(y|x) from p2(y|x). Another category defines a distance metric through the embedding of probability measures in another space (typically the reproducing kernel Hilbert space or RKHS). A notable example is the Maximum Mean Discrepancy (MMD) [Gretton et al., 2012b] which has attracted much attention in recent years due to its solid mathematical foundation. However, MMD always requires high computational burden and carefully hyper-parameter (e.g., the kernel width) tuning [Gretton et al., 2012a]. Moreover, computing the distance of the embeddings of two conditional distributions in RKHS still remains a challenging problem [Ren et al., 2016].

∗Contact Author

Different from previous efforts, we propose a simple statistic to quantify the discrepancy between two conditional distributions. It directly operates on the cone of symmetric positive semidefinite (SPS) matrix to avoid the estimation of the underlying distributions. To strengthen the discriminative power of our statistic, we make use of the correntropy function [Santamar´ıa et al., 2006], which has demonstrated its effectiveness in non-Gaussian signal processing [Liu et al., 2007], to explicitly incorporate higher order information in the data. We demonstrate the power of our statistic and establish its connections to prior art. Three solid examples of machine learning applications are presented to demonstrate the effectiveness and the superiority of our statistic.

2 Background Knowledge

2.1 BregmanMatrix Divergence and Its Computation

A symmetric matrix is positive semidefinite (SPS) if all its eigenvalues are non-negative. We denote S n+ the set of all n × n SPS matrices, i.e., S n+ = {A ∈ Rn×n |A = AT , A ≽ 0}. To measure the nearness between two SPS matrices, a reliable choice is the Bregman matrix divergence [Kulis et al., 2009]. Specifically, given a strictly convex, differentiable function ϕ that maps matrices to the extended real numbers, the Bregman divergence from the matrix ρ to the matrix σ is defined as:

where tr(A) denotes the trace of matrix A.

When ϕ(σ) = tr(σ log σ - σ), where log σ is the matrix logarithm, the resulting Bregman divergence is:

which is also referred to von Neumann divergence in quantum information theory [Nielsen and Chuang, 2011]. Another important matrix divergence arises by taking ϕ(σ) =

-logdet σ, in which the resulting Bregman divergence reduces to:

and is commonly called the LogDet divergence.

2.2 Correntropy Function: AGeneralized Correlation Measure

The correntropy function of two random variables x and y is defined as [Santamar´ıa et al., 2006]:

where E denotes mathematical expectation, κ is a positive definite kernel function, and FX,Y (x, y) is the joint distribution of (X,Y). One widely used kernel function is the Gaussian kernel given by:

Taking Taylor series expansion of the Gaussian kernel, we have:

Therefore, correntropy involves all the even moments 1 of random variable e = x - y. Furthermore, increasing the kernel size σ makes correntropy tends to the correlation of x and y [Santamar´ıa et al., 2006].

A similar quantity to the correntropy is the centered correntropy:

where FX(x) and FY (y) are the marginal distributions of X and Y , respectively.

The centered correntropy can be interprested as a nonlinear counterpart of the covariance in RKHS [Chen et al., 2017]. Moreover, the following property serves as the basis of our test statistic that will be introduced later.

1 A different kernel would yield a different expansion, for instance the sigmoid kernel κ(x, y) = tanh(〈x, y〉 + θ) admits an expansion in terms of the odd moments of its argument.

Property 2.1. Given n random variables x1, x2, · · · , xn and any set of real numbers α1, α2, · · · , αn, for any symmetric positive definite kernel κ(x, y), the centered correntropy matrix C defined as C(i, j) = U(xi, xj) is always positive semidefinite, i.e., C ∈ S n+.

Proof. By Property 2 in [Rao et al., 2011], U(x, y) is a symmetric positive semidefinite function, from which it follows that:

C is obviously symmetric. This concludes our proof.

In practice, data distributions are unknown and only a finite number of observations {xi, yi}Ni=1 are available, which leads to the sample estimator of centered correntropy 2 [Rao et al., 2011]:

3 Methods

3.1 ProblemFormulation

We have two groups of observations {(x1i , y1 i )}N1 i=1 and {(x2i , y2 i )}N2 i=1 that are assumed to be independently and identically distributed (i.i.d.) with density functions p1(x, y) and p2(x, y), respectively. y is a dependent variable that takes values in R, and x is a vector of explanatory variables that takes values in Rp . Typically, the conditional distribution p(y|x) are unknown and unspecified. The aim of this paper is to suggest a test statistic to measure the nearness between p1(y|x) and p2(y|x).H0 : Pr (p1(y|x) = p2(y|x)) = 1

3.2 Our Statistic and the Conditional Test

We define the divergence from p1(y|x) to p2(y|x) as:

where Cxy ∈ S p+1 + denotes the centered correntropy matrix of the random vector concatenated by x and y, and Cx ∈ S p+ denotes the centered correntropy matrix of x. Obviously, Cx is a submatrix of Cxy by removing the row and column associated with y.

Although Eq. (10) is assymetic itself, one can easily achieve symmetry by taking the form:

2 Throughout this work, we determine kernel width σ with the Silverman's rule of thumb [Silverman, 1986].

Algorithm 1 Test the conditional distribution divergence (CDD) based on the matrix Bregman divergence

Input: Two groups of observations S1 = {(x1i , y1 i )}N1 i=1 and S2 = {(x2i , y2 i )}N2 i=1 , xi ∈ Rp , yi ∈ R; Dϕ,B; Permutation number P; Significant rate η.

Output: Test decision (Is H0 : Pr (p1(y|x) = p2(y|x)) = 1 T rue or F alse?).

• 1: Measure CDD d0 on S1 and S2 with Eq. (11).
• 2: for t = 1 to P do
• 3: (S1 t , S2 t ) ← random split of S1 ⋃ S2 .
• 4: Measure CDD dt on S1 t and S2 t with Eq. (11).
• 5: end for
• 6: if 1+∑Pt=1 1[d0≤dt] 1+P ≤ η then
• 7: decision←False
• 8: else
• 9: decision←True
• 10: end if
• 11: return decision

Given Eq. (10) and Eq. (11), we design a simple permutation test to distinguish p1(y|x) from p2(y|x). Our test methodology is shown in Algorithm 1, where 1 indicates an indicator function. The intuition behind this scheme is that if there is no difference on the underlying conditional distributions, the test statistic on the ordered split (i.e., d0) should not deviate too much from that of the shuffled splits (i.e., {dt}Pt=1 ).

4 Conditional Divergence Properties

We present useful properties of our statistic (i.e., Eq. (10) or Eq. (11)) based on the LogDet divergence. In particular we show it is non-negative (when evaluated on covariance matrices) which reduces to zero when two data sets share the same linear regression function. We also perform Monte Carlo simulations to investigate its detection power and establish its connection to prior art.

Proof. The proof follows immediately from the fact that Dℓd(Σ1yx ||Σ2yx ) - Dℓd(Σ1x ||Σ2x ) is the KL divergence between p1(y|x) and p2(y|x) for Gaussian distributions with zero mean and of covariances Σ1yx ,Σ2yx (see Eq. (15)) and DKL(p1(y|x)||p2(y|x)) ≥ 0.

Property 4.2. Let x1 ∼ N(µ1x ,Σ1x ) and x2 ∼ N(µ2x ,Σ2x ) be two input processes of full rank covariance matrices of size p × p. For a common linear system, defined by a full rank matrix W of size r × p such that y = Wx, we have:

Proof. Denote M = ∣∣∣∣ IpW ∣∣∣∣, where Ip denotes an identity matrix, we have Σixy = WΣixW T , from which we have:

where we used Property 12 and Lemma 5 of [Kulis et al., 2009], since range(Σixy ) = p ≤ r + p.

4.1 Power Test

Our aim here is to examine if our statistic is really suitable for quantifying the discrepancy between two conditional distributions. Motivated by [Zheng, 2000; Fan et al., 2006], we generate four groups of data that have distinct conditional distributions. Specifically, in model (a), the dependent variable y is generated by y = 1 + ∑pi=1 xi + ϵ, where ϵ denotes standard normal distribution. In model (b), y = 1+ ∑pi=1 xi +ψ, where ψ denotes the standard Logistic distribution. In model (c), y = 1+ ∑pi=1 logxi+ϵ. In model (d), y = 1 + ∑pi=1 logxi + ψ. For each model, the input distribution is an isotropic Gaussian.

To evaluate the detection power of our statistic on any two models, we randomly generate 500 samples from each model which has the same dimensionality m on explanatory variable x. We then use Algorithm 1 (P = 500, η = 0.1) to test if our statistic can distinguish these two data sets. We repeat this procedure with 100 independent runs and use the percentage of success as the detection power. The conditional KL divergence (see Eq. (14)) estimated with an adaptive kNN estimator [Wang et al., 2009] is implemented as a baseline competitor. Table 1 summarizes the power test results when p = 3 and p = 30. Although all methods perform good for p = 3, our test statistic is significantly more powerful than conditional KL divergence in high-dimensional space.

We also depict in Fig. 1 the power of our statistic in case I: model (a) against model (b) and case II: model (c) against model (d), with respect to different kernel widths. For case I (the first row), larger width is preferred. This is because the second-order information dominates in the linear model. For case II (the second row), we need smaller width to capture more higher order information in highly nonlinear model.

4.2 Relation to Previous Efforts Gaussian Data

As mentioned earlier, the centered correntropy matrix C evaluated with a Gaussian kernel encloses all even higher order statistics of pairwise dimensions of data. If we replace C with its second-order counterpart (i.e., the covariance matrix Σ), we have:

Taking Eq. (3) into Eq. (12), we obtain:

On the other hand, the KL divergence between two conditional distributions on a pair of random variables satisfies the following chain rule (see proof in Chapter 2 of [Cover and Thomas, 2012]):

Suppose the data is Gaussian distributed, i.e., p1(x) ∼ N(µ1x ,Σ1x ), p2(x) ∼ N(µ2x ,Σ2x ), p1(x, y) ∼ N(µ1xy ,Σ1xy ), p2(x, y) ∼ N(µ2xy , Σ2xy ), Eq. (14) has a closed-form expression (see proof in [Duchi, 2007]):

Comparing Eq. (13) with Eq. (15), it is easy to find that our baseline variant reduces to the conditional KL divergence under Gaussian assumption. The only difference is that the conditional KL divergence contains a Mahalanobis Distance term on the mean, which can be interpreted as the first-order information of data.

Beyond Gaussian Data

Gaussian assumption is always over-optimistic. If we stick to the conditional KL divergence, the probability estimation becomes inevitable, which is notoriously difficult in highdimensional space [Nagler and Czado, 2016]. By making use the correntropy function, our statistic avoids the estimation of the underlying distribution, but it explicitly incorporates the higher order information which was lost in Eq. (12).

5 Machine Learning Applications

We present three solid examples on machine learning applications to demonstrate the performance improvement in the state-of-the-art (SOTA) methodologies gained by our conditional divergence statistic.

5.1 Multitask Learning

Consider an input set X and an output set Y and for simplicity that X ∈ Rp , Y ∈ R. Tasks can be viewed as T functions ft, t = 1, · · · , T, to be learned from given data {(xti, yti) : i = 1, · · · ,Nt, t = 1, · · · , T} ⊆ X ×Y, where Nt is the number of samples in the t-th task. These tasks may be viewed as drawn from an unknown joint distribution of tasks, which is the source of the bias that relates the tasks. Multitask learning is the paradigm that aims at improving the generalization performance of multiple prediction problems (tasks) by exploiting potential useful information between related tasks.

Visualizing Task-relatedness

We first test if our proposed statistic is able to reveal the relatedness among multiple tasks. To this end, we select data from 29 tasks that are collected from various landmine fields 3 . Each object in a given data set is represented by a

3 http://www.ee.duke.edu/∼lcarin/LandmineData.zip.

9-dimensional feature vector and the corresponding binary label (1 for landmine and 0 for clutter). The landmine detection problem is thus modeled as a binary classification problem.

Among these 29 data sets, 1-15 correspond to regions that are relatively highly foliated and 16-29 correspond to regions that are bare earth or desert. We measure the task-relatedness with the conditional divergence Dϕ,B(p1(y|x) : p2(y|x)) and demonstrate their pairwise relationships in a task-relatedness matrix. Thus we expect that there are approximately two clusters in the task-relatedness matrix corresponding to two classes of ground surface condition. Fig. 2 demonstrates the visualization results generated by our proposed statistic and the conditional KL divergence estimated with an adaptive kNN estimator [Wang et al., 2009]. We also compare our method with MTRL [Zhang and Yeung, 2010], a widely used objective to learn multiple tasks and their relationships in a convex manner. Our vN (C) and vN (Σ) clearly demonstrate two clusters of tasks. By contrast, both MTRL and conditional KL divergence indicate a few misleading relations (i.e., tasks in the same group surface condition have high divergence values).

Multitask Learning over Graph Structure

In the second example, we demonstrate how our statistic improves the performance of the Convex Clustering Multi-Task Learning (CCMTL) [He et al., 2019], a SOTA method for learning multiple regression tasks. The learning objective of CCMTL is given by:

where W = [wT 1 ,wT 2 , · · · ,wT T ] ∈ RT×p is the weight matrix constitutes of the learned linear regression coefficients in each task, G is the graph of relations over all tasks (if two tasks are related, then there is an edge to connect them), and λ is a regularization parameter.

CCMTL requires an initialization of the graph structure G0. In the original paper, the authors set G0 as a kNN graph on the prediction models learned independently for each task. In this sense, the task-relatedness between two tasks is modeled as the ℓ2 distance of their independently learned linear regression coefficients.

We replace the na¨ıve ℓ2 distance with our proposed statistic to reconstruct the initial kNN graph for CCMTL and test its performance on a real-world Parkinson's disease data set [Tsanas et al., 2009]. We have 42 tasks and 5, 875 observations, where each task and observation correspond to a prediction of the symptom score (motor UPDRS) for a patient and a record of a patient, respectively. Fig. 3 depicts the prediction root mean square errors (RMSE) under different train/test ratios. To highlight the superiority of CCMTL and our improvement, we also compare it with MSSL [Gonc¸alves et al., 2016], a recently proposed alternative objective to MTRL that learns multiple tasks and their relationships with a Gaussian graphical model. The performance of MTRL is relatively poor, and thus omitted here. Obviously, our statistic improves the performance of CCMTL with a large margin, especially when training samples are limited. Note that, we did not observe performance difference between centered correntropy matrix C and covariance matrix Σ. This is probably because the higher order information is weak or because the Silverman's rule of thumb is not optimal to tune kernel width here.

5.2 Concept Drift Detection

One important assumption underlying common classification models is the stationarity of the training data. However, in real-world data stream applications, the joint distribution p(x, y) between the predictor x and response variable y is not stationary but drifting over time. Concept drift detection approaches aim to detect such drifts and adapt the model so as to mitigate the deteriorating classification performance over time.

Formally, the concept drift between time instants t0 and t1 can be defined as pt1 (x, y) ̸= pt2 (x, y) [Gama et al., 2014]. From a Bayesian perspective, concept drifts can manifest two fundamental forms of changes: 1) a change in the marginal probability pt(x) or pt(y); and 2) a change in the posterior probability pt(y|x). Although any two-sample test (e.g., [Gretton et al., 2012b]) on pt(x) or pt(y) is an option to detect concept drift, existing studies tend to prioritize detecting posterior distribution change, because it clearly indicates the optimal decision rule [Dries and R¨uckert, 2009].

Error-based methods constitute a main category of existing concept drift detectors. These methods keep track of the online classification error or error-related metrics of a baseline classifier. A significant increase or decrease of the monitored statistics may suggests a possible concept drift. Unfortunately, the performance of existing manually designed errorrelated statistics depends on the baseline classifier, which makes them either perform poorly across different data distributions or difficult to be extended to the multi-class classification scenario.

Different from prevailing error-based methods, our statistic operates directly on the conditional divergence Dϕ,B(pt1 (y|x) : pt2 (y|x)), which makes it possible to fundamentally solve the problem of concept drift detection (without any classifier). To this end, we test the conditional distribution divergence in two consecutive sliding windows (using Algorithm 1) at each time instant t. A reject of the null hypothesis indicates the existence of a concept drift. Note that the same permutation test procedure has been theoretically and practically investigated in PERM [Harel et al., 2014], in which the authors use classification error as the test statistic. Interested readers can refer to [Harel et al., 2014; Yu and Abraham, 2017] for more details on concept drift detection with permutation test.

We evaluate the performance of our method against four SOTA error-based concept drift detectors (i.e., DDM [Gama et al., 2004], EDDM [Baena-Garcıa et al., 2006], HDDM [Fr´ıas-Blanco et al., 2014], and PERM) on two realworld data streams, namely the Digits08 [Sethi and Kantardzic, 2017] and the AbruptInsects [dos Reis et al., 2016]. Among the selected competitors, HDDM represents one of the best-performing detectors, whereas PERM is the most similar one to ours. The concept drift detection results and the stream classification results over 30 independent runs are summarized in Table 2. Our method always enjoys the highest recall

(b) Abrupt Insect

Table 2: Quantitative metrics on real-world data sets. The Precision, Recall and Delay denote the concept drift detection precision value, recall value and detection delay, whereas the Accuracy denotes the classification accuracy in the testing set (%).

values and the shortest detection delay. Note that, the classification accuracy is not as high as we expected. One possible reason is that the short detection delay makes us do not have sufficient number of samples to retrain the classifier.

5.3 Feature Selection

Our final application is information-theoretic feature selection. Given a set of variables S = {X1,X2, · · · ,Xn}, feature selection refers to seeking a small subset of informative variables S⋆ ⊂ S, such that S⋆ contains the most relevant yet least redundant information about a desired variable Y . From an information-theoretic perspective, this amounts to maximize the mutual information term I(y; S⋆ ). Suppose we are now given a set of "useless" features ˜ S that has the same size as S⋆ but has no predictive power to y, Eq. (17) suggests that instead of maximizing I(y; S⋆ ), one can resort to maximize DKL(P (y|S⋆ )‖P (y| ˜ S)) as an alternative.

the last line is by our assumption that ˜ S has no predictive power to y such that P(y| ˜ S) = P (y).

Motivated by the generic equivalence theorem between KL divergence and the Bregman divergence [Banerjee et al., 2005], we optimize our statistic in a greedy forward search manner (i.e., adding the best feature at each round) and generate "useless" feature set ˜ S by randomly permutating S⋆ as conducted in [Franc¸ois et al., 2007].

We perform feature selection on two benchmark data sets [Brown et al., 2012]. For both data sets, we select 10 features and use the linear Support Vector Machine (SVM) as the baseline classifier. We compare our method with three popular information-theoretic feature selection methods that target maximizing I(y; S⋆ ), namely MIFS [Battiti, 1994], FOU [Brown et al., 2012], and MIM [Lewis, 1992]. The classification accuracies with respect to different number of selected features (averaged over 10 fold cross-validation) are presented in Fig. 4. As can be seen, methods on maximizing conditional divergence always achieve comparable performances to those mutual information based competitors. This result confirms Eq. (17). If we look deeper, it is easy to see that our statistic implemented with von Neumann divergence achieves the best performance on both data sets. Moreover, by incorporating higher order information, centered correntropy performs better than its covariance based counterpart.

6 Conclusions

We propose a simple statistic to quantify conditional divergence. We also establish its connections to prior art and illustrate some of its fundamental properties. A natural extension to incorporate higher order information of data is also presented. Three solid examples suggest that our statistic can offer a remarkable performance gain to SOTA learning algorithms (e.g., CCMTL and PERM). Moreover, our statistic enables the development of alternative solutions to classical machine learning problems (e.g., classifier-free concept drift detection and feature selection by maximizing conditional divergence) in a fundamental manner.

Future work is twofold. We will investigate more fundamental properties of our statistic. We will also apply our statistic to more challenging problems, such as the continual learning [Kirkpatrick et al., 2017] which also requires the knowledge of task-relatedness.