A novel formulation of independence detection based on the sample characteristic function

Accepted version of the article published in the 26th European Signal Processing Conference, 2018. DOI: 10.23919/EUSIPCO.2018.8553062. Available at: http://ieeexplore.ieee.org A novel formulation of Independence Detection based on the Sample Characteristic Function Ferran de Cabrera, Student Member, IEEE and Jaume Riba, Senior Member, IEEE Signal Theory and Communications Department, Technical University of Catalonia (UPC) {ferran.de.cabrera, jaume.riba}@upc.edu Abstract—A novel independence test for continuous random sequences is proposed in this paper. The test is based on seeking for coherence in a particular fixed-dimension feature space based on a uniform sampling of the sample characteristic function of the data, providing significant computational advantages over kernel methods. This feature space relates uncorrelation and independence, allowing to analyze the second order statistics as it is encountered in traditional signal processing. As a result, the possibility of utilizing well known correlation tools arises, motivating the usage of Canonical Correlation Analysis as the main tool for detecting independence. Comparative simulation results are provided using a model based on fading AWGN channels. Index Terms—Independence Detection, Characteristic Function, Mutual Information, ITL, Kernel, HSIC, CCA. I. INTRODUCTION Methods for measuring and/or detecting statistical dependence, a problem that has a long history, has recently found applications in a wide variety of areas ([10][19] and references therein). Some important examples from the broad literature are the Maximal Information Coefficient (MIC) [11], the distance Correlation (dCor) [16], the Mutual Information Dimension (MID) [15] and the Hilbert-Schmidt Independence Criterion (HSIC) [5]. Detecting dependence has a close link with machine learning, which, in general terms, tries to discover structure in data. In both cases the objective is to construct a non-linear mapping from the data onto a higherdimensional space (referred to as feature space) in which, in some sense, the original problem becomes simpler, resorting the Cover’s theorem [2]. For the specific problem of detecting statistical dependence, the objective of the mapping is to move the original problem into a simpler problem of measuring second order statistics. Kernel methods ([17][6][1]), which are based on the rigorous mathematical framework of reproducing kernel Hilbert spaces (RKHS), provide an elegant solution for this mapping known as kernel trick: replacing inner products by positivedefinite kernels is equivalent to implicitly processing the data on an extremely high dimensional feature space. However, apart from the need of regularization to control the smoothness of the results, an important drawback of kernel methods is their compromised computational complexity for large data sets, which requires the use of some special techniques to This work has been supported by project TEC2016-76409-C2-1-R (WINTER) (Ministerio de Economia y Competividad, Spanish National Research Plan) and project 2017-SGR-578 (AGAUR) (Catalan Government). avoid the need of accessing to all pairs of input data points, such as the Incomplete Cholesky decomposition [13]. This paper proposes an alternative to the kernel trick based on a particular feature space. Following the rationale of the Information Theoretic Learning (ITL) issued in [14], the objective is to define a detector based on quadratic measures on the space of the Characteristic Function (CF), and we will show that we are capable of defining an independence test based solely on first and second order statistics in this space. The mapping is based on uniformly sampling the CF whose fixed-dimension is independent of the data size, allowing to ameliorate the complexity issues raised by the kernels. This feature space resembles an analog memory-less frequencymodulation of the data, giving rise to second-order autocorrelation matrices, which, jointly with the link between uncorrelation and independence, allows to expand the study to the Canonical Correlation Analysis (CCA), which is shown to avoid the well-known overfitting problem that arises in Kernel CCA (K-CCA) [7]. The rationale is similar to the study of independent component analysis under a CF approach present in [3], in which the empirical joint CF is compared to the marginal ones. However, while their work is based on expanding the analysis of the CF through its cumulants from the Taylor expansion around zero, we are more interested on studying the sampling of the feature space in order to handle the problem with classical signal processing. Finally, we provide a comparative analysis of the method with respect to other known methods, in addition with a genie-aided detector based on the sample Mutual Information (MI) estimation with known Probability Density Function (PDF). To this end, three data models based on Gaussian mixtures are defined at different Signal to Noise Ratio (SNR) regimes, which are inspired on fading AWGN communications channels. II. SAMPLE CHARACTERISTIC FUNCTION AS A FEATURE SPACE Let fx,y(x, y)be the joint PDF of the random variables x, y ∈R. The joint CF is φx,y(ω1;ω2) = Z ZR fx,y(x, y)ejω1x+jω2ydxdy =Eejω1x+jω2y=Eejω1xejω2y(1) for any ω1, ω2∈R. A necessary and sufficient condition for independence is the separability of the CF φx,y(ω1;ω2) = c 2018 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. φx(ω1)φy(ω2) = Eejω1xEejω2y. Considering the random variables z1=ejω1xand z2=ejω2ythe condition for independence can be written as E(z1−E[z1]) (z2−E[z2])∗= 0 (2) which constitutes the requisite for uncorrelation. This property succeeds in moving the independence condition to an equivalent uncorrelation one, which will be used as the basis of independence criteria. A. Sample characteristic function analysis Consider {xi}i=1...N aN-length sequence composed by i.i.d. data observations drawn from the random variable xwith PDF fx(x)and associated CF φx(ω) = R∞ −∞ fx(x)ejxωdx. Since the CF exists at all values of ω, we are interested in studying under which conditions a finite dimensional feature space is sufficient for detecting independence. Taking a specific sample xi, the windowed sample CF is φxi(ω) = ejωxiG(ω)(3) being G(ω)an even and positive window function with G(0) = 1 and R|G(ω)|dω < ∞. The function G(ω) assumes that most of the informative part of the CF is concentrated around the window bandwidth σ. Let us consider an additional mapping based on the Fourier Transform (FT): Φi(β) = ∞ Z −∞ φxi(ω)e−jβωdω = ∞ Z −∞ G(ω)e−jω(β−xi)dω (4) Defining g(β) = RG(ω)e−jωβdω as the FT of G(ω), then we have Φi(β) = g(β−xi)(5) corresponding to a real function g(β)centered around the data value xi. Let us now define the first and second order statistics. The expectation over the data samples is µΦ(β) = E[Φi(β)] = ∞ Z −∞ g(β−x)fx(x)dx = (fx∗g) (β)(6) and its empirical version is given by ˆµΦ(β) = 1 N N X i=1 g(β−xi)(7) which is in fact a kernel density estimate [9] given g(β)is a unitary area function. This result ensures that the expected value of the FT of the sample CF is a smoothed estimate of the original PDF. We define the second order statistic as: rΦ(β1, β2) = E[Φi(β1) Φ∗ i(β2)] = ∞ Z −∞ g(β1−x)g∗(β2−x)fx(x)dx (8) It is worth noting that for g(β)narrow enough, a data association phenomenon occurs such that rΦ(β1, β2)≈0∀|β1−β2|> ε (9) for any small positive ε > 0, whose insight is to measure closeness between data, and resembling the pairwise sample interactions present in [10]. If we restrict Eq. (8) to lay in the same β, then we can express it as rΦ(β) = ∞ Z −∞ h(β−x)fx(x)dx = (fx∗h) (β)(10) and its empirical estimate as ˆrΦ(β) = 1 N N X i=1 h(β−xi)(11) with h(β) = |g(β)|2. Since the window G(ω)is expected to be wide enough in order to cover most of the informative part of the CF, then the support of its squared FT h(β)tends to decrease. If we approximate h(β)as a delta-like function such that h(β)≃δ(β)with Rδ(β)dβ =k1, then rΦ(β)≃(fx∗δ) (β) = kfx(β)(12) which is a scaled version of the original PDF. Let us assume a new set of i.i.d. data {yi}i=1...N extracted from the random variable ywith PDF fy(y). The crosscorrelation between the two data sets evaluated in the feature space can now be expressed in terms of (6) and (8). To this end, consider d(xi) = g(β1−xi)and d(yi) = g(β2−yi) and formulate the squared Pearson coefficient as follows: ρ2=(E[d(xi)d(yi)] −E[d(xi)] E[d(yi)])2 (E[d2(xi)] −E2[d(xi)]) (E[d2(yi)] −E2[d(yi)]) (13) where E[d(xi)d(yi)] is an approximate of the joint PDF fxy (x, y): E[d(xi)d(yi)] = ∞ Z −∞ ∞ Z −∞ g(β1−x)g(β2−y) ×fx,y (x, y)dxdy ≃fx,y (β1, β2)(14) and the expectation E[d(xi)] is the same as in Eq. (6), which is approximately fx(β1), and analogously E[d(yi)] ≃ fy(β2). By analyzing Eq. (13) with the same rationale and using Eq. (12) we get: ρ2→ ∞ Z −∞ ∞ Z −∞ (fx,y (β1, β2)−fx(β1)fy(β2))2 (kfx(β1)−f2 x(β1)) (kfx(β2)−f2 x(β2))dβ1dβ2 (15) Interestingly, if we let kto increase without limit, then the squared Pearson coefficient tends to ρ2−→ k→∞ ∞ Z −∞ ∞ Z −∞ (fx,y (β1, β2)−fx(β1)fy(β2))2 k2fx(β1)fy(β2)dβ1dβ2 =1 k2 ∞ Z −∞ ∞ Z −∞ fx,y (x, y) fx(x)fy(y)−12 fx(x)fy(y)dxdy (16) which results in the kernel-free expression of the Normalized Cross-Covariance Operator (NOCCO) developed by Fukumizu et al. ([4]) that studies the dependence structure based on normalized cross-covariance operators. Eq. (16) constitutes an indicator of independence between random variables through the joint PDF and the product of the marginal ones, and compiles it as a correlation measure based on the squared Pearson coefficient. The possibility of drawing correlation metrics from the first and second order statistics in the proposed feature space, as well as the capability of linking uncorrelation and independence in terms of the CF, motivates us to explore the CCA from the uniform sampling of the CF, as we will see in Section III. B. Uniform sampling In order to operate with the CF from data samples, we will proceed by applying a uniform sampling of the sample feature space φxi(ω1), and equivalently for φyi(ω2), by evaluating ω1and ω2for nα being α∈Rthe sampling period and n=−(Nz−1) /2, ..., (Nz−1) /2with Nzany real number. Then, we construct the mapped data matrices Zx= [zx(1),zx(2), ..., zx(N)] ∈CNz×Nand Zy= [zy(1),zy(2), ..., zy(N)] ∈CNz×Ncontaining in their i-th column the mapping of the CF for a given sample1: [zx(i)]n=ejnαxi [zy(i)]n=ejnαyi(17) This sampling assumes that the relevant CF space is contained in a support of Nzα, and so the feature space dimension Nzis sufficiently large for the detection task but still lower than the data dimensionality (NzN) in order to reduce the computational complexity from ON2present in kernel methods to O(NzN). Moreover, this assumption also minimizes the aliasing induced by the sampling and the Inverse Fourier Transform intrinsic in the CF. Considering the space between replicas to be 1/α, the inverse of the sampling period, then we must ensure that α < 1/σ, with σcorresponding to the window bandwidth. Note also that it is possible to sample asymmetrically with n= 1, ..., Nzsince the CF is Hermitian. The sample mean vectors evaluated in the feature map are then ¯ zx=Zx1N/N and ¯ zy=Zy1N/N, auto-covariance matrices ˆ Cx=ZxP ZH x/N and ˆ Cy=ZyP ZH y/N, and the cross-covariance matrices ˆ Cxy =ZxP ZH y/N and ˆ Cyx = ˆ CH xy being 1NaN×1column vector of ones, P=IN− 1N1T N/N the projection matrix and INthe N×Nidentity matrix. III. CANONICAL CORRELATION ANALYSIS Given two data matrices Zx,Zy∈CNz×N, the CCA problem is based on finding the canonical variables uand v 1Both even and odd Nzare permitted. Although an odd real value induces bias given by [z(i)]n=0 = 1, it does not impact the performance of the detector since the bias is known and can be fixed. Algorithm 1 Inversion algorithm pseudocode. Input: Zmatrix of dimension Nz×N Initialize: R−1(0) = δINzwith δ1 for i= 0,1, ..., N −1do g= 1 + zH(i)R−1(i)z(i) R−1(i+ 1) ←R−1(i)z(i)zH(i)R−1(i)/g end that maximize the correlation between the canonical variates wx=uTZxand wy=vTZy, i.e. arg max u,v ρ2= uTˆ Cxyv ruTˆ CxuvTˆ Cyv  2 (18) where ρ2corresponds to the squared maximum canonical correlation. The solution to this problem, and thereby the canonical variates, is given by the Singular Value Decomposition (SVD) of the Coherence matrix, a whitened version of the cross-covariance matrix, defined as ∆ = ˆ C−1/2 xˆ Cxy ˆ C−1/2 y(19) The eigenvalues of ∆are the non-negative roots of the eigenvalues λm=1,...,Nzof the squared Coherence matrix K= ∆∆Hwith λ1≥λ2≥ ··· ≥ λNz≥0([12]), corresponding to all the canonical correlations ρ2 mwith the highest being ρ2. Thus, the maximum squared canonical correlation ρ2is given by the largest eigenvalue of K, and so we can evaluate the main solution as: λ1=λmax ˆ C−1/2 xˆ Cxy ˆ C−1 yˆ CH xy ˆ C−1/2 x =λmax ˆ C−1 xˆ Cxy ˆ C−1 yˆ Cyx(20) However, the auto-covariance matrices ˆ Cxand ˆ Cyare illconditioned, making the inversion task difficult under real simulations. In this case, the inversion of the auto-covariance matrices can be solved through Algorithm 1 by means of the Woodbury matrix identity ([8]). Finally, we can build the Characteristic Canonical Independence Test (CCIT) as the maximum canonical correlation: CCIT (x,y, Nz, α) = λmax ˆ C−1 xˆ Cxy ˆ C−1 yˆ Cyx(21) Then, the proposed detector is CCIT (x,y, Nz, α) H1 ≷ H0 γ(22) being γthe detector threshold, H0the independence hypothesis and H1the dependence hypothesis. Note that we can also evaluate the trace of Kas a detector such as CCIT2(x,y, Nz, α) = tr (K)since it contains all the subsequent canonical correlations ρ2 mfrom the CCA and permitting to avoid the SVD of K. But, this test contributes with a bias given the addition of the lesser eigenvalues. However, since the canonical correlations are concentrated at the largest eigenvalues, the impact to the test is minimal as we will see in the numerical results. Figure 1. Contour plots of the joint probability functions with SNR = 8 and p= 0.95 in the third model. IV. NUMERICAL RESULTS In order to evaluate the detector, the data is modeled as: x(i) = hx(i)−√p+p1−pzx(i)+wx(i) y(i) = hy(i)√p+p1−pzy(i)+wy(i)(23) corresponding to fast fading AWGN channels, with zx, zy∼ N(0,1) independent from each other under independence hypothesis H0and equal under dependence hypothesis H1, and the noise wx, wy∼ N 0, σ2 wmutually independent for both hypothesis. The intention is to generate data with different kinds of associations through the parameter pand the channels hxand hy, following the examples in [11], while maintaining a close relation with a communication channel. Thus, we are able to relate the difficulty of detecting independence with the SNR. Additionally, each one of the three models, observed in Figure 1, holds the same marginal PDFs for both independence and dependence hypothesis. The basic and first model consists on the presence of correlation and dependence, represented by a linear structure of the data. For this model, the parameters are p= 0, hx(i) = hy(i)=1and its MI is given by I(x, y) = 1 2ln 1+4SNR2/(1 + 4SNR)with SNR = 1/2σ2 w, which for high SNR is asymptotically the MI of an AWGN channel I(x, y) = 1 2ln (1 + SNR). The second model is based on the dependence of variables while being uncorrelated. The intention is to evaluate if the CCIT can properly detect dependence under non-linear associations of the data. For this model, p= 0 and hx(i), hy(i)are binary distributed with equal probability [−1,1] and independent to each other, which act as fast fading channels. Lastly, the third model provides uncorrelation but with an increased difficulty of detection. For this model hx(i)and hy(i)are distributed as the second model but the parameter 0<p<1determines the difficulty of detection, being p= 0 akin to the second model and p= 1 nullifies the dependence. The generated data is shaped as a Gaussian Mixture Model (GMM), allowing an explicit expression of the PDFs and 10-1 100101 SNR 0 0.5 1 1.5 2 Mutual information Second model Third model p=0.95 First model 0.5 In(1+SNR) Figure 2. True MI vs. SNR of the three models obtained by the unbiased genie-aided estimator given in Eq. (24) with N= 106. ensuring a consistent estimate of the MI based on the sample mean as ˆ I=1 N N X i=1 ln fx,y (xi, yi) fx(xi)fx(yi)(24) We will compare the performance of the detector with a benchmark based on a threshold from the genie-aided statistic in Eq. (24), whose values for each model as can be seen in Figure 2. Note also that the proposed models provide an increased difficulty for detecting dependence, as appreciated by a higher threshold effect of MI vs. SNR. Figure 3 provides a comparative analysis of the Probability of Detection PDfor multiple detectors for a fixed false alarm probability of 5%. The results were obtained with 400 Monte Carlo simulations and the parameters optimized for each case. For the HSIC, a Gaussian kernel was used with the bandwidth recommended in [5], which is the heuristic for these cases. Regarding the parameters of the CCIT, Nz= 21 and the sampling period α= 0.001 were chosen, providing a trade-off between computational complexity and leakage in the form of aliasing in the FT. Figure 3a shows the performance of the first model, focusing on the low SNR scenario. It can be seen that for linear associations of the data, HSIC and dCor slightly outperform the CCIT with all of them presenting a high robustness to noise, while MIC and MID are very sensitive to it. On the other hand, Figure 3b and 3c show that CCIT improves the detection rate with respect to other detectors for uncorrelated data. Specifically, CCIT provides a higher convergence for low SNR values, whose SNR threshold for a PDclose to 100% is near to that of the genie-aided estimate. The underperformance of the CCIT for linear dependencies is caused by an increase of the condition number of the autocovariance matrices given precisely by this linear relation, since the data mapped in the feature space is confined in a small dimensional subspace. Under non-linear dependencies, the feature space spreads the subspace dimension and the inverse of the matrices can be better estimated. This affectation can be solved by reducing the sampling period α, thus limiting the mapping into the desired subspace. Then, the CCIT is Figure 3. Probability of Detection for PF A = 5% vs SNR from 0.125 to 64 of the three models. In dashed blue, the PDof Eq. (24) estimate. (N= 300, p = 0.95, Nz= 21, α = 0.001) more aware of the linear dependence and the PDis improved. However, in the figures a fixed value of αis chosen due to the seek of universality of the detector ([19]), and so to be fair in the comparison with the case of non-linear dependencies. Alternatively, the detection rate could be improved by avoiding the inverse process and measuring the canonical correlation through adaptive methods instead of CCA [18]. V. CONCLUSIONS In this paper we have presented a feature space based on the uniform sampling of the marginal sample CFs. This mapping allowed us to take advantage of the relation between uncorrelation and independence from the chosen feature space, and we have used a correlation measure to build the detector. We have shown that by applying a specific condition on the sample characteristic function we can relate, at least asymptotically, the proposed CCIT with other quadratic measures of information. The numerical simulations proved that this method is capable of improving the detection rate with respect to others in the literature, specifically under non-linear associations, and with a clear computational advantage over kernel signal processing. However, an exhaustive analysis on adaptive methods is needed in order to solve the numerical problems when inverting the ill-conditioned matrices. Moreover, the study of dependence under the proposed feature space for non i.i.d. data provides an interesting future research line. REFERENCES [1] F.R. Bach and M.I. Jordan. Kernel independent component analysis. Journal of Machine Learning, 3, pages 1-48, 2002. [2] T.M. Cover. Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition. IEEE Transactions on Electronic Computers, Volume: EC-14, Issue: 3, pages 326-334, 1965. [3] J. Eriksson and V. Koivunen. Characteristic-function-based independent component analysis. Signal Processing, 83(10):2195 – 2208, 2003. [4] K. Fukumizu, A. Gretton, X. Sun, and B. Schölkopf. Kernel measures of conditional dependence. In J. C. Platt, D. Koller, Y. Singer, and S. T. Roweis, editors, Advances in Neural Information Processing Systems 20, pages 489–496. Curran Associates, Inc., 2008. [5] A. Gretton, A. Smola, O. Bousquet, and B. Schölkopf. Measuring statistical dependence with hilbert-schmidt norms. Technical report, Max Planck Institute for Biological Cybernetics, Tübingen, Germany, 2005. [6] A. Gretton, A. Smola, R. Herbrich, O. Bousquet, and B. Schölkopf. Kernel methods for measuring independence. Journal of Machine Learning Research, 2005. [7] D. R. Hardoon, S. Szedmak, and J. Shawe-Taylor. Canonical correlation analysis: An overview with application to learning methods. Neural Comput., 16(12):2639–2664, 2004. [8] Monson H. Hayes. Statistical Digital Signal Processing and Modeling. Number 541-553. John Wiley & Sons, Inc., New York, NY, USA, 1st edition, 1996, pp 541-553. [9] E. Parzen. On a Estimation of a Probability Density Function and Mode, volume 33, 1065-1067. Ann. Math. Statist., 1962. [10] J.C. Príncipe. Information Theoretic Leraning. Springer, 2010. [11] D.N. Reshef, Y.A. Reshef, H.K. Finucane, S.R. Grossman, G. McVean, P.J. Turnbaugh, E.S. Lander, M. Mitzenmacher, and P.C. Sabeti. Detecting novel associations in large data sets. Science, 334(6062):1518–1524, December 2011. [12] P.J. Schreier. A unifying discussion of correlation analysis for complex random vectors. IEEE Transactions on Signal Processing, 56(4), 2008. [13] S. Seth and J.C. Príncipe. On speeding up computation in information theoretic learning. In International Joint Conference on Neural Networks, IJCNN, 2009. [14] S. Seth, M. Rao, I. Park, and J.C. Príncpie. A unified framework for quadratic measures of independence. IEEE Transactions on Signal Processing, 59(8), 2011. [15] M. Sugiyama and K.M Borgwardt. Measuring statistical dependence via the mutual information dimension. In Proceedings of the 23rd International Joint Conference on Artificial Intelligence (IJCAI 2013), 2013. [16] G.J Székely and M.L Rizzo. Brownian distance covariance. Annals of Statistics, 2010. [17] S.V. Vaerenbergh. Kernel Methods for Nonlinear Identification, Equalization and Separation of Signals. PhD thesis, University of Cantabria, Department of Communications Engineering, 2009. [18] J. Vía, I. Santamaría, and J. Pérez. A robust algorithm for adaptive canonical correlation analysis. IEEE International Conference on Acoustics, Speech and Signal Processing, 2005. [19] Q. Wang, S.R.Kulkarni, and S.Verdú. Universal Estimation of Information Measures for Analog Sources. Now publishers, 2009.