Computational Information Geometry for Binary Classification of High-Dimensional Random Tensors ^†

Evaluating the performance limit for the “Gaussian information plus noise” binary classification problem is a challenging research topic, see for instance [1,2,3,4,5,6,7]. Given a binary hypothesis problem, the Bayes’ decision rule is based on the principle of the largest posterior probability. Specifically, the Bayesian detector chooses the alternative hypothesis ${\mathcal{H}}_1$ if $\mathrm{Pr}\left({\mathcal{H}}_1\right|\mathbf{y})>\mathrm{Pr}\left({\mathcal{H}}_0\right|\mathbf{y})$ for a given N-dimensional measurement vector $\mathbf{y}$ or the null hypothesis ${\mathcal{H}}_0$, otherwise. Consequently, the optimal decision rule can often only be derived at the price of a costly numerical computation of the log posterior-odds ratio [3] since an exact calculation of the minimal Bayes’ error probability $P_e^{\left(N\right)}$ is often intractable [3,8]. To circumvent this problem, it is standard to exploit well-known bounds on $P_e^{\left(N\right)}$ based on information theory [9,10,11,12,13]. In particular, the Chernoff information [14,15] is asymptotically (in N) relied on the exponential rate of $P_e^{\left(N\right)}$. It turns out that the Chernoff information is very useful in many practically important problems as for instance, distributed sparse detection [16], sparse support recovery [17], energy detection [18], multi-input and multi-output (MIMO) radar processing [19,20], network secrecy [21], angular resolution limit in array processing [22], detection performance for informed communication systems [23], just to name a few. In addition, the Chernoff information bound can be tight for a minimal s-divergence over parameter $s\in (0,1)$. Generally, this step requires solving numerically an optimization problem [24] and often leads to a complicated and uninformative expression of the optimal value of s. To circumvent this difficulty, a simplified case of $s=1/2$ is often used corresponding to the well-known Bhattacharyya divergence [13] at the price of a less accurate prediction of $P_e^{\left(N\right)}$. In information geometry, parameter s is often called $\alpha$, and the s-divergence is the so-called Chernoff $\alpha$-divergence [24].

The tensor decomposition theory is a timely and prominent research topic [25,26]. Confronting the problem of extracting useful information from a massive and multidimentional volume of measurements, it is shown that tensors are extremely relevant. In the standard literature, two main families of tensor decomposition are prominent, namely the Canonical Polyadic Decomposition (CPD) [26] and the Tucker decomposition (TKD)/HOSVD (High-Order SVD) [27,28]. These approaches are two possible multilinear generalization of the Singular Value Decomposition (SVD). A natural generalization to tensors of the usual concept of rank for matrices is called the CPD. The tensorial/canonical rank of a P-order tensor is equal to the minimal positive integer, say R, of unit rank tensors that must be summed up for perfect recovery. A unit rank tensor is the outer product of P vectors. In addition, the CPD has remarkable uniqueness properties [26] and involves only a reduced number of free parameters due to the constraint of minimality on R. Unfortunately, unlike the matrix case, the set of tensors with fixed (tensorial) rank is not close [29,30]. This singularity implies that the problem of the computation of the CPD is mathematically ill-posed. The consequence is that its numerical computation remains non trivial and is usually done using suboptimal iterative algorithms [31]. Note that this problem can sometimes be avoided by exploiting some natural hidden structures in the physical model [32]. The TKD [28] and the HOSVD [27] are two popular decompositions being an alternative to the CPD. Under this circumstance, alternative definition of rank is required, since the tensorial rank based on CPD scenario is no longer appropriate. In particular, stardard definition of multilinear rank defined as the set of positive integers $\{R_1,\dots ,R_P\}$ where each integer, $R_p$, is the usual rank of the p-th mode. Following the Eckart-Young theorem at each mode level [33], this construction is non-iterative, optimal and practical. In real-time computation [34] or adaptively computation [35], it is shown that this approach is suitable. However, in general, the low (multilinear) rank tensor based on this procedure is suboptimal [27]. More precisely, for tensors of order strictly greater than two, a generalization of the Eckart-Young theorem does not exist.

The classification performance of a multilinear tensor following the CPD and TKD can be derived and studied. It is interesting to note that the classification theory for tensors is very under studied. Based on our knowledge on the topic, only the publication [36] tackles this problem in the context of radar multidimensional data detection. A major difference with this publication is that their analysis is based on the performance of a low rank detection after matched filtering.

More precisely, we consider two cases where the observations are either (1) a noisy rank-R tensor admitting a Q-order CPD with large factors of size $N_q\times R$, i.e., for $1\le q\le Q$, $R,N_q\to \infty$ with $R^{1/q}/N_q$ converging towards a finite constant, or (2) a noisy tensor admitting a TKD of multilinear $(M_1,\dots ,M_Q)$-rank with large factors of size $N_q\times M_q$, i.e., for $1\le q\le Q$, where $N_q,M_q\to \infty$ with $M_q/N_q$ converging towards a finite constant. A standard approach for zero-mean independent Gaussian core and noise tensors, is to define the Signal to Noise Ratio by $\mathrm{SNR}={\sigma }_{\mathbf{s}}^2/{\sigma }^2$ where ${\sigma }_{\mathbf{s}}^2$ and ${\sigma }^2$ are the variances of the vectorized core and noise tensors, respectively. So, the binary classification can be described in the following way:

Under the null hypothesis ${\mathcal{H}}_0$, $\mathrm{SNR}=0$, meaning that the observed tensor contains only noise. Conversely, the alternative hypothesis ${\mathcal{H}}_1$ is based on $\mathrm{SNR}\ne 0$, meaning that there exists a multilinear signal of interest. First note that there exists a lack of contribution dealing with classification performance for tensors. Since the exact derivation of the error probability is intractable, the performance of the classification of the core tensor random entries is hard to evaluate. To circumvent this audible difficulty, based on computational information geometry theory, we consider the Chernoff Upper Bound (CUB), and the Fisher information in the context of massive measurement vectors. The error exponent can be minimized at $s^{\star }$, which corresponds to the reachable tightest CUB. In general, due to the asymmetry of the s-divergence, the Bhattacharyya Upper Bound (BUB)—Chernoff Information calculated at $s^{\star }=1/2$—cannot solve this problem effectively. As a consequence, we rely on a costly numerical optimization strategy to find $s^{\star }$. However, with respect to different Signal to Noise Ratios (SNR), we provide simple analytical expressions of $s^{\star }$, thanks to the so-called Random Matrix Theory (RMT). For low SNR, analytical expressions of the Fisher information are given. Note that the analysis of the Fisher information in the context of the RMT has been only studied in recent contributions [37,38,39] for parameter estimation. For larger SNR, analytic and simple expression of the CUB for the CPD and the TKD are provided.

We note that Random Matrix Theory (RMT) has attracted both mathematicians and physicists since they were first introduced in mathematical statistics by Wishart in 1928 [40]. When Wigner [41] introduced the concept of statistical distribution of nuclear energy levels, the subject has started to earn prominence. However, it took until 1955 before Wigner [42] introduced ensembles of random matrices. Since then, many important results in RMT were developed and analyzed, see for instance [43,44,45,46] and the references therein. In the last two decades, research on RMT has been constantly published.

Finally, let us underline that many arguments of this paper differ from the works presented in [47,48]. In [47], we tackled the problem of detection using Chernoff Upper Bound in data of type matrix in the double asymptotic regime. In [48], we established the detection problem in tensor data by analyzing the Chernoff Upper Bound. In [48], we assumed that the tensor follows the Canonical Polyadic Decomposition (CPD), we gave some analysis of Chernoff Upper Bound when the rank of the tensor is much smaller than the dimensions of the tensor. Since [47,48] are conference papers, some proofs have been omitted due to limited space. Therefore, this full paper may share the ideas in [47,48] on Information Geometry (s-divergence, Chernoff Upper Bound, Fisher Information, etc.), but completes [48] in a more general asymptotic regime. Moreover, in this work, we give new analysis in both scenarios (SNR small and large) whereas [48] did not, and the important and difficult new tensor scenario of the Tucker decomposition is considered. This is in our view the main difference because the CPD is a particular case of the more general decomposition of TucKer. Indeed, in the CPD, the core tensor is assumed to be diagonal.

The organization of the paper is as follows: In the second section, we introduce some definitions, tensor models, and the Marchenko-Pastur distribution from random matrix theory. The third section is devoted to present Chernoff Information for binary hypothesis test. The fourth section gives the main results on Fisher Information and the Chernoff bound. The numerical simulation results are given in the fifth section. We conclude our work by giving some perspectives in the Section 6. Finally, several proofs of the paper can be found in the appendix.

The Marchenko-Pastur distribution was introduced half a century ago [45] in 1967, and plays a key role in a number of high-dimensional signal processing problems. To help the reader, in this section, we introduce some fundamental results concerning large empirical covariance matrices. Let ${\left({\mathbf{v}}_n\right)}_{n=1,\dots ,N}$ a sequence of i.i.d zero mean Gaussian random M-dimensional vectors for which $\mathbb{E}\left({\mathit{v}}_n{\mathit{v}}_n^T\right)={\sigma }^2{\mathit{I}}_M$. We consider the empirical covariance matrix which can be also written as where matrix ${\mathit{W}}_N$ is defined by ${\mathit{W}}_N=\frac{1}{\sqrt{N}}[{\mathit{v}}_1,\dots ,{\mathit{v}}_N]$. ${\mathit{W}}_N$ is thus a Gaussian matrix with independent identically distributed $\mathcal{N}(0,\frac{{\sigma }^2}{N})$ entries. When $N\to +\infty$ while M remains fixed, matrix ${\mathit{W}}_N{\mathit{W}}_N^T$ converges towards ${\sigma }^2{\mathit{I}}_M$ in the spectral norm sense. In the high dimensional asymptotic regime defined by it is well understood that $∥{\mathit{W}}_N{\mathit{W}}_N^T-{\sigma }^2{\mathit{I}}_M∥$ does not converge towards 0. In particular, the empirical distribution ${\widehat{\nu }}_N=\frac{1}{M}{\sum }_{m=1}^M{\delta }_{{\widehat{\lambda }}_{m,N}}$ of the eigenvalues ${\widehat{\lambda }}_{1,N}\ge \dots \ge {\widehat{\lambda }}_{M,N}$ of ${\mathit{W}}_N{\mathit{W}}_N^T$ does not converge towards the Dirac measure at point $\lambda ={\sigma }^2$. More precisely, we denote by ${\nu }_{c,{\sigma }^2}$ the Marchenko-Pastur distribution of parameters $(c,{\sigma }^2)$ defined as the probability measure with ${\lambda }^{-}={\sigma }^2{(1-\sqrt{c})}^2$ and ${\lambda }^{+}={\sigma }^2{(1+\sqrt{c})}^2$. Then, the following result holds.

We also observe that Theorem 1 remains valid if ${\mathbf{W}}_N$ is not necessarily a Gaussian matrix whose i.i.d. elements have a finite fourth order moment (see e.g., [43]). Theorem 1 means that when ratio $\frac{M}{N}$ is not small enough, the eigenvalues of the empirical spatial covariance matrix of a temporally and spatially white noise tend to spread out around the variance of the noise, and that almost surely, for N large enough, all the eigenvalues are located in a neighbourhood of interval $[{\lambda }^{-},{\lambda }^{+}]$. See Figure 2 and Figure 3.

We denote by $\mathrm{SNR}={\sigma }_{\mathbf{s}}^2/{\sigma }^2$ and $p_i(·)=p(·|{\mathcal{H}}_i)$ with $i\in \{0,1\}$. The binary classification of the random signal based on the equi-probable binary hypothesis test, $\mathbf{s}$, is where ${\mathbf{\Sigma }}_0={\sigma }^2{\mathbf{I}}_N$ and ${\mathbf{\Sigma }}_1={\sigma }^2\left(\mathrm{SNR}\times \mathbf{\Phi }{\mathbf{\Phi }}^T+{\mathbf{I}}_N\right)$. The null hypothesis data-space (${\mathcal{H}}_0$) is defined as ${\mathcal{X}}_0=\mathcal{X}\backslash {\mathcal{X}}_1$ where is the alternative hypothesis (${\mathcal{H}}_1$) data-space. Following the above expression, the log-likelihood ratio test $\mathsf{\Lambda }\left({\mathbf{y}}_N\right)$ and the binary classification threshold ${\tau }^{\prime }$ are given by where $\det (·)$ and $log(·)$ are respectively the determinant and the natural logarithm.

We note that the estimated hypothesis $\widehat{\mathcal{H}}$ is associated to $p({\mathbf{y}}_N|\widehat{\mathcal{H}})=\mathcal{N}\left(\mathbf{0},\mathbf{\Sigma }\right)$. Therefore, the expected log-likelihood ratio is defined by where is the Kullback-Leibler Divergence (KLD) [10]. The expected log-likelihood ratio test admits to a simple geometric characterization based on the difference of two KLDs [8]. However, it is often difficult to evaluate the performance of the test via the minimal Bayes’ error probability $P_e^{\left(N\right)}$, since its expression cannot be determined analytically in closed-form [3,8].

The minimal Bayes’ error probability conditionally to vector ${\mathbf{y}}_N$ is defined as where $P_{i,i^{\prime }}=\mathrm{Pr}\left({\mathcal{H}}_i\right|{\mathbf{y}}_N\in {\mathcal{X}}_{i^{\prime }})$.

According to [24], the relation between the Chernoff Upper Bound and the (average) minimal Bayes’ error probability $P_e^{\left(N\right)}=\mathbb{E}\mathrm{Pr}\left(\mathrm{Error}\right|{\mathbf{y}}_N)$ is given by where the (Chernoff) s-divergence for $s\in (0,1)$ is given by in which $M_X\left(t\right)=\mathbb{E}exp[t\times X]$ is the moment generating function (mgf) of variable X. The error exponent, denoted by $\tilde{\mu }\left(s\right)$, is given by the Chernoff information which is an asymptotic characterization on the exponentially decay of the minimal Bayes’ error probability. The error exponent is derived thanks to the Stein’s lemma according to [13]

As parameter $s\in (0,1)$ is free, the CUB can be tightened by minimizing this parameter:

Finally, using Equations (5) and (7), the Chernoff Upper Bound (CUB) is obtained. Instead of solving Equation (7), the Bhattacharyya Upper Bound (BUB) is calculated by Equation (5) and by fixing $s=1/2$ . Therefore we have the following relation of order:

From now on, to simplify the presentation and the numerical results later on, we denote by for all $s\in [0,1]$, the opposites of the log-moment generating function and its limit.

In the small deviation regime, we assume that $\delta \mathrm{SNR}$ is a small deviation of the SNR. The new binary hypothesis test is where $\mathbf{\Sigma }\left(x\right)=x\times \mathbf{\Phi }{\mathbf{\Phi }}^T+\mathbf{I}.$ The s-divergence in the small $\mathrm{SNR}$ deviation scenario is written as

According to Lemma 2, the optimal s-value at low $\mathrm{SNR}$ is $s^{\star }\stackrel{\delta \mathrm{SNR}\ll 1}{=}\frac{1}{2}$. At contrary, the optimal s-value for larger $\mathrm{SNR}$ is given by the following lemma.

The measurement tensor follows a noisy Q-order tensor of size $N_1\times \dots \times N_Q$ can be expressed as where $\mathcal{N}$ is the noise tensor whose entries are assumed to be centered i.i.d. Gaussian, i.e., ${\left[\mathcal{N}\right]}_{n_1,\dots ,n_Q}\sim \mathcal{N}(0,{\sigma }^2)$ and the core tensor $\mathcal{X}$ follows either CPD or TKD given by Section 2.1.2 and Section 2.1.3, respectively. The vectorization of Equation (10) is given by where $\mathbf{n}=\mathrm{vec}\left({\mathbf{N}}_{\left(1\right)}\right)$ and $\mathbf{x}=\mathrm{vec}\left({\mathbf{X}}_{\left(1\right)}\right)$. Note that ${\mathbf{Y}}_{\left(1\right)}$, ${\mathbf{N}}_{\left(1\right)}$ and ${\mathbf{X}}_{\left(1\right)}$ are respectively the first unfolding matrices given by Definition 4 of tensors $\mathcal{Y},\mathcal{N}$ and $\mathcal{X}$,

We recall that in the CPD case, matrix ${\mathbf{\Phi }}^{\odot }={\mathbf{\Phi }}^{\left(Q\right)}\odot \dots \odot {\mathbf{\Phi }}^{\left(1\right)}$ and ${\left({\mathbf{\Phi }}^{\left(q\right)}\right)}_{q=1,\dots ,Q}$ are matrices of size $N_q\times R$. In the following, we assume that matrices ${\mathbf{\Phi }}_{q=1,\dots ,Q}^{\left(q\right)}$ are random matrices with Gaussian $\mathcal{N}(0,\frac{1}{N_q})$ variate entries. We evaluate the behavior of $\frac{{\mu }_N\left(s\right)}{N}$ when ${\left(N_q\right)}_{q=1,\dots ,Q}$ converge towards $+\infty$ at the same rate and that $\frac{R}{N}$ converges towards a non zero limit.

In the TKD case, we recall that matrix ${\mathbf{\Phi }}^{\otimes }={\mathbf{\Phi }}^{\left(Q\right)}\otimes \dots \otimes {\mathbf{\Phi }}^{\left(1\right)}$, with ${\left({\mathit{\varphi }}^{\left(q\right)}\right)}_{1\le q\le Q}$ are $N_q\times M_q$ dimensional matrices. We still assume that matrices ${\mathbf{\Phi }}_{q=1,\dots ,Q}^{\left(q\right)}$ are random matrices with Gaussian $\mathcal{N}(0,\frac{1}{N_q})$ entries.