Common Information Components Analysis

Abstract

Wyner’s common information is a measure that quantifies and assesses the commonality between two random variables. Based on this, we introduce a novel two-step procedure to construct features from data, referred to as Common Information Components Analysis (CICA). The first step can be interpreted as an extraction of Wyner’s common information. The second step is a form of back-projection of the common information onto the original variables, leading to the extracted features. A free parameter $\gamma$ controls the complexity of the extracted features. We establish that, in the case of Gaussian statistics, CICA precisely reduces to Canonical Correlation Analysis (CCA), where the parameter $\gamma$ determines the number of CCA components that are extracted. In this sense, we establish a novel rigorous connection between information measures and CCA, and CICA is a strict generalization of the latter. It is shown that CICA has several desirable features, including a natural extension to beyond just two data sets.

1. Introduction

Understanding relations between two (or more) sets of variates is key to many tasks in data analysis and beyond. To approach this problem, it is natural to reduce each of the sets of variates separately in such a way that the reduced descriptions fully capture the commonality between the two sets, while suppressing aspects that are individual to each of the sets. This permits to understand the relation between the two sets without obfuscation. A popular framework to accomplish this task follows the classical viewpoint of dimensionality reduction and is referred to as Canonical Correlation Analysis (CCA) [1]. CCA seeks the best linear extraction, i.e., we consider linear projections of the original variates. Via the so-called Kernel trick, this can be extended to cover arbitrary (fixed) function classes.

Wyner’s common information is a well-known and established measure of the dependence of two random variables. Intuitively, it seeks to extract a third random variable such that the two random variables are conditionally independent given the third, but at the same time that the third variable is as compact as possible. Compactness is measured in terms of the mutual information that the third random variable retains about the original two. The resulting optimization problem is not a convex problem (because the constraint set is not a convex set), and therefore, not surprisingly, closed-form solutions are rare. A natural generalization of Wyner’s common information is obtained by replacing the constraint of conditional independence by a limit on the conditional mutual information. If the limit is set equal to zero, we return precisely to the case of conditional independence. Exactly like mutual information, Wyner’s common information and its generalization are endowed with a clear operational meaning. They characterize the fundamental limits of data compression (in the Shannon sense) for a certain network situation.

1.1. Related Work

Connections between CCA and Wyner’s common information have been explored in the past. It is well known that, for Gaussian vectors (standard, non-relaxed), Wyner’s common information is attained by all of the CCA components together, see [2]. This has been further interpreted, see, e.g., [3]. Needless to say, having all of the CCA components together essentially amounts to a one-to-one transform of the original data into a new basis. It does not yet capture the idea of feature extraction or dimensionality reduction. To put our work into context, it is only the relaxation of Wyner’s common information [4,5] that permits to conceptualize the sequential, one-by-one recovery of the CCA components, and thus, the spirit of dimensionality reduction.

CCA also appears in a number of other problems related to information measures and probabilistic models. For example, in the so-called Gaussian information bottleneck problem, the optimizing solution can be expressed in terms of the CCA components [6], and an interpretation of CCA as a (Gaussian) probabilistic model was presented in [7].

Generalizations of CCA have appeared before in the literature. The most prominent is built around maximal correlation. Here, one seeks arbitrary remappings of the original data in such a way as to maximize their correlation coefficient. This perspective culminates in the well-known alternating conditional expectation (ACE) algorithm [8].

Feature extraction and dimensionality reduction have a vast amount of literature attached to them, and it is beyond the scope of the present article to provide a comprehensive overview. In a part of that literature, information measures play a key role. Prominent examples are independent components analysis (ICA) [9] and the information bottleneck [10,11], amongst others. More recently, feature extraction alternations via information theory are presented in [12,13]. In [12], the estimation of Rényi’s quadratic entropy is studied, whereas, in [13], standard information theoretic measures such as Kullback–Leibler divergence are used for fault diagnosis. Other slightly related feature extraction methods that perform dimensionality reduction on a single dataset include [14,15,16,17,18,19,20]. More concretely, in [14], a sparse Support Vector Machine (SVM) approach is used for feature extraction. In [15], feature extraction is performed via regression by using curvilinearity instead of linearity. In [16], compressed sensing is used to extract features when the data have a sparse representation. In [17], an invariant mapping method is invoked to map the high-dimensional data to low-dimensional data that is based on a neighborhood relation. In [18], feature extraction is performed for a partial learning of the geometry of the manifold. In [19], distance correlation measure (a measure with similar properties as the regular Pearson correlation coefficient) is proposed as a new feature extraction method. In [20], kernel principal component analysis is used to perform feature extraction and allow for the extraction of nonlinearities. In [21], feature extraction is done by a robust regression based approach and, in [22], a linear regression approach is used to extract features.

1.2. Contributions

The contributions of our work are the following:

• We introduce a novel suit of algorithms, referred to as CICA. These algorithms are characterized by a two-step procedure. In the first step, a relaxation of Wyner’s common information is extracted. The second step can be interpreted as a form of projection of the common information back onto the original data so as to obtain the respective features. A free parameter $\gamma$ is introduced to control the complexity of the extracted features.
• We establish that, for the special case where the original data are jointly Gaussian, our algorithms precisely extract the CCA components. In this case, the parameter $\gamma$ determines how many of the CCA components are extracted. In this sense, we establish a new rigorous connection between information measures and CCA.
• We present initial results on how to extend CICA to more than two variates.
• Via a number of paradigmatic examples, we illustrate that, for discrete data, CICA gives intuitively pleasing results while other methods, including CCA, do not. This is most pronounced in a simple example with three sources described in Section 7.1.

1.3. Notation

A bold capital letter such as $\mathbf{X}$ denotes a random vector, and $\mathbf{x}$ its realization. The probability distribution function of random variable X will be denoted by $p_X$ or $p\left(x\right)$ depending on the context. A non-bold capital letter such as K denotes a (fixed) matrix, and $K^H$ its Hermitian transpose. Specifically, $K_{\mathbf{X}}$ denotes the covariance matrix of the random vector $\mathbf{X}.$ $K_{\mathbf{X}\mathbf{Y}}$ denotes the covariance matrix between random vectors $\mathbf{X}$ and $\mathbf{Y}.$ Let $\mathcal{P}$ be the set of all probability distribution, discrete or continuous depending on the context. Let us denote with $I_n$ the identity matrix of dimension $n\times n$ and $0_n$ the zero matrix of dimension $n\times n$. We denote by $\mathcal{L}{\left[f\left(x\right)\right]}_x$ the lower convex envelope of $f\left(x\right)$ with respect to x and for random variables $\mathcal{L}{\left[f\left(X\right)\right]}_{p_X}$ is the lower convex envelope of $f\left(X\right)$ with respect to $p_X$. We denote by $h_b\left(x\right):=-xlogx-(1-x)log1-x$ the binary entropy for $0\le x\le 1$.

1.4. A Simple Example with Synthetic Data

To set the stage and under the guise of an informal problem statement, let us consider a simple example involving synthetic data. Specifically, we consider two-dimensional data, that is, the vectors $\mathbf{X}$ and $\mathbf{Y}$ are of length 2. The goal is to extract, separately from each of the two, a one-dimensional description in such a way as to extract the commonality between $\mathbf{X}$ and $\mathbf{Y}$ while suppressing their individual features. For simplicity, in the present artificial example, we will assume that the entries of the vectors only take value in a small finite set, namely, $\{0,1,2,3\}.$ To illustrate the point, we consider the following special statistical model:

$$\begin{array}{cc}\hfill \mathbf{X}& =\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)=\left(\begin{array}{c}U\oplus X_2\\ X_2\end{array}\right),\hfill \end{array}$$

(1)

and

$$\begin{array}{cc}\hfill \mathbf{Y}& =\left(\begin{array}{c}{Y}_1\\ Y_2\end{array}\right)=\left(\begin{array}{c}U\oplus Y_2\\ Y_2\end{array}\right),\hfill \end{array}$$

(2)

where $U,X_2,$ and $Y_2$ are mutually independent uniform random variables over the set $\{0,1,2,3\}$ and ⊕ denotes addition modulo 4.

The reason for this special statistical structure is such that it is obvious what should be extracted, namely, $\mathbf{X}$ should be reduced to $U,$ and $\mathbf{Y}$ should also be reduced to $U.$ This reduces both $\mathbf{X}$ and $\mathbf{Y}$ to “one-dimensional” descriptions, and these one-dimensional descriptions capture precisely the dependence between $\mathbf{X}$ and $\mathbf{Y}.$ In this simple example, all the commonality between $\mathbf{X}$ and $\mathbf{Y}$ is captured by $U.$ More formally, conditioned on $U,$ the vectors $\mathbf{X}$ and $\mathbf{Y}$ are conditionally independent.

The interesting point of this example is that any pair of components of $\mathbf{X}$ and $\mathbf{Y}$ are independent of each other, such as, for example, $X_1$ and $Y_1.$ Therefore, the joint covariance matrix of the merged vector $(\mathbf{X},\mathbf{Y})$ is a scaled identity matrix. This implies that any method that only uses the covariance matrix as input, including CCA, cannot find any commonalities between $\mathbf{X}$ and $\mathbf{Y}$ in this example.

By contrast, the algorithmic procedure discussed in the present paper will correctly extract the desired answer. In Figure 1, we show numerical simulation outcomes for a couple of approaches. Specifically, in (a), we can see that, in this particular example, CCA fails to extract the common features. This, of course, was done on purpose: For the synthetic data at hand, the global covariance matrix is merely a scaled identity matrix, and since CCA’s only input is the covariance matrix, it does not actually do anything in this example. In (b), we show the performance of the approximate gradient-descent based implementation of the CICA algorithm proposed in this paper, as detailed in Section 6. In this simple example, this precisely coincides with the ideal theoretical performance of CICA as in a Generic Procedure 1, but, in general, the gradient-descent based implementation is not guaranteed to find the ideal solution.

At this point, we should stress that, for such a simple example, many other approaches would also lead to the same, correct answer. One of them is maximal correlation. In that perspective, one seeks to separately reduce $\mathbf{X}$ and $\mathbf{Y}$ by applying possibly nonlinear functions $f(·)$ and $g(·)$ in such a way as to maximize the correlation between $f\left(\mathbf{X}\right)$ and $g\left(\mathbf{Y}\right).$ Clearly, for the simple example at hand, selecting $f\left(\mathbf{X}\right)=X_1\oplus X_2$ and $g\left(\mathbf{Y}\right)=Y_1\oplus Y_2$ leads to correlation one, and is thus a maximizer.

Finally, the present example is also too simplistic to express the finer information-theoretic structure of the problem. One step up is the example presented in Section 5 below, where the commonality between $\mathbf{X}$ and $\mathbf{Y}$ is not merely an equality (the component U above), but rather a probabilistic dependency.

2. Wyner’s Common Information and Its Relaxation

The main framework and underpinning of the proposed algorithm is Wyner’s common information and its extension, which is briefly reviewed in the sequel, along with its key properties.

2.1. Wyner’s Common Information

Wyner’s common information is defined for two random variables X and Y of arbitrary fixed joint distribution $p(x,y).$

Definition 1

([23]). For random variables X and Y with joint distribution $p(x,y),$ Wyner’s common information is defined as

$$\begin{array}{cc}\hfill C(X;Y)& =\underset{p\left(w\right|x,y)}{inf}I(X,Y;W) such that I(X;Y|W)=0.\hfill \end{array}$$

(3)

Wyner’s common information satisfies a number of interesting properties. We state some of them below in Lemma 1 for a generalized definition given in Definition 2.

We note that explicit formulas for Wyner’s common information are known only for a small number of special cases. The case of the doubly symmetric binary source is solved completely in [23] and can be written as

$$\begin{array}{cc}\hfill C(X;Y)& =1+h_b\left(a_0\right)-2h_b\left(\frac{1-\sqrt{1-2a_0}}{2}\right),\hfill \end{array}$$

(4)

where $a_0$ denotes the probability that the two sources are unequal (assuming without loss of generality $a_0\le \frac{1}{2}$). In this case, the optimizing W in Equation (3) can be chosen to be binary. Further special cases of discrete-alphabet sources appear in [24].

Moreover, when X and Y are jointly Gaussian with correlation coefficient $\rho ,$ then $C(X;Y)=\frac{1}{2}log\frac{1+\left|\rho \right|}{1-\left|\rho \right|}.$ Note that, for this example, $I(X;Y)=\frac{1}{2}log\frac{1}{1-{\rho }^2}.$ This case was solved in [25,26] using a parameterization of conditionally independent distributions, and we have recently found an alternative proof that also extends to the generalization of Wyner’s common information discussed in the next subsection [5].

2.2. A Natural Relaxation of Wyner’s Common Information

A natural generalization of Wyner’s common information (Definition 1) is obtained by replacing the constraint of conditional independence with a limit $\gamma$ on the conditional mutual information, leading to the following:

Definition 2

(from [4,5,27]). For random variables X and Y with joint distribution $p(x,y),$ we define

$$\begin{array}{cc}\hfill C_{\gamma }(X;Y)& =\underset{p\left(w\right|x,y)}{inf}I(X,Y;W) such that I(X;Y|W)\le \gamma .\hfill \end{array}$$

(5)

This definition appears in slightly different form in Wyner’s original paper (Section 4.2 in [23]), where an auxiliary quantity $\Gamma ({\delta }_1,{\delta }_2)$ is defined satisfying $C_{\gamma }(X;Y)=H(X,Y)-\Gamma (0,\gamma ).$ The above definition first appears in [4]. Comparing Definitions 1 and 2, we see that, for $\gamma =0,$ we have $C_0(X;Y)=C(X;Y),$ the regular Wyner’s common information. In this sense, one may refer to $C_{\gamma }(X;Y)$ as relaxed Wyner’s common information.

In line with the discussion following Definition 1, it is not surprising that explicit solutions to the optimization problem in Definition 2 are very rare. In fact, the only presently known general solution concerns the case of jointly Gaussian random variables [5]. The corresponding formula is given below in Theorem 1.

By contrast, the case of the doubly symmetric binary source remains open. An upper bound for this case is given by choosing the auxiliary W as

$$\begin{array}{cc}\hfill W& =\left\{\begin{array}{cc}X\oplus V,\hfill & \mathrm{if}X=Y,\hfill \\ U,\hfill & \mathrm{if}X\ne Y,\hfill \end{array}\right.\hfill \end{array}$$

(6)

where V is Bernoulli with probability $\alpha$ and U is Bernoulli with probability $\frac{1}{2}$. Thus, the upper bound is

$$\begin{array}{cc}\hfill C_{\gamma }(X;Y)& \le I(X,Y;W)=1-(1-a_0)h_b\left(\alpha \right)-a_0,\hfill \end{array}$$

(7)

where $\alpha \ge {\alpha }_W$ is chosen such that

$$\begin{array}{cc}\hfill I(X;Y|W)& =2h_b\left((1-\alpha )(1-a_0)+\frac{a_0}{2}\right)-(1-a_0)h_b\left(\alpha \right)-a_0-h_b\left(a_0\right)=\gamma ,\hfill \end{array}$$

(8)

where

$$\begin{array}{c}\hfill {\alpha }_W=\frac{{(1-\sqrt{1-2a_0})}^2}{4(1-a_0)}.\end{array}$$

(9)

Numerical studies (Section 3.5 in [28]) suggest that this upper bound is tight, but no formal proof is available to date.

The following lemma summarizes some basic properties of $C_{\gamma }(X;Y).$

Lemma 1

(partially from [5]). $C_{\gamma }(X;Y)$ satisfies the following basic properties:

1. 

$C_{\gamma }(X;Y)\ge max\{I(X;Y)-\gamma ,0\}.$

2. 

Data processing inequality: If $X-Y-Z$ forms a Markov chain,

then $C_{\gamma }(X;Z)\le min\{C_{\gamma }(X;Y),C_{\gamma }(Y;Z)\}.$

3. 

$C_{\gamma }(X;Y)$ is a convex and continuous function of γ for $\gamma \ge 0.$

4. 

Tensorization: For n independent pairs ${\left\{(X_i,Y_i)\right\}}_{i=1}^n,$ we have that

$$C_{\gamma }(X^n;Y^n)=min\sum _{i=1}^n{C}_{{\gamma }_i}(X_i;Y_i),$$

where the min is over all non-negative ${\left\{{\gamma }_i\right\}}_{i=1}^n$ satisfying ${\sum }_{i=1}^n{\gamma }_i=\gamma .$

5. 

If $Z-X-Y$ forms a Markov chain, then $C_{\gamma }((X,Z);Y)=C_{\gamma }(X;Y).$

6. 

The cardinality of $\mathcal{W}$ may be restricted to $\left|\mathcal{W}\right|\le \left|\mathcal{X}\right|\left|\mathcal{Y}\right|+1.$

7. 

If $f(·)$ and $g(·)$ are one-to-one functions, then $C_{\gamma }(f\left(X\right);g\left(Y\right))=C_{\gamma }(X;Y).$

8. 

For discrete $X,$ we have $C_{\gamma }(X;X)=max\{H\left(X\right)-\gamma ,0\}.$

Proofs of items 1–4 are given in [5], and the proofs of items 5–8 are given in Appendix A.

2.3. The Non-Convexity of the Relaxed Wyner’s Common Information Problem

It is important to observe that the optimization problem of Definition 2 is not a convex problem. First, we observe that $I(X,Y;W)$ is indeed a convex function of $p\left(w\right|x,y),$ which is a well-known fact, see, e.g., (Theorem 2.7.4 in [29]). The issue is with the constraint set. The set of distributions $p\left(w\right|x,y)$ for which $I(X;Y|W)\le \gamma$ is not a convex set. To provide some intuition for the structure of this set, let us consider $I(X;Y|W)$ as a function of $p\left(w\right|x,y)$, and examine its (non-)convexity. The relation between the two is described by the epigraph

$$\begin{array}{c}\hfill \mathrm{epigraph}\left\{I\right(X;Y\left|W\right)\}=\{\left(p\right(w|x,y),\gamma ):p(w|x,y)\in \mathcal{P},\gamma \ge I(X;Y|W\left)\right\}.\end{array}$$

(10)

The function $I(X;Y|W)$ is convex in $p\left(w\right|x,y)$ if and only if its epigraph is a convex set which would imply that the set of distributions $p\left(w\right|x,y)$ for which $I(X;Y|W)\le \gamma$ is also convex. Now, we present an example that $I(X;Y|W)$ is not a convex function of $p\left(w\right|x,y)$.

Example 1.

Let the distributions $p(x,y),p_1\left(w\right|x,y),p_2\left(w\right|x,y)$ be

$$\begin{array}{c}\hfill p(x,y)=\left[\begin{array}{cc}\frac{2}{5}& \frac{1}{10}\\ \frac{1}{10}& \frac{2}{5}\end{array}\right],p_1\left(w\right|x,y)=\left[\begin{array}{cccc}\frac{1}{4}& \frac{1}{4}& \frac{1}{4}& \frac{1}{4}\\ \frac{3}{4}& \frac{3}{4}& \frac{3}{4}& \frac{3}{4}\end{array}\right],p_2\left(w\right|x,y)=\left[\begin{array}{cccc}\frac{1}{2}& \frac{3}{4}& \frac{3}{4}& \frac{3}{4}\\ \frac{1}{2}& \frac{1}{4}& \frac{1}{4}& \frac{1}{4}\end{array}\right],\end{array}$$

(11)

respectively. For this example, one can evaluate numerically that, under $p_1\left(w\right|x,y),$ we have $I_{p_1}(X;Y|W)<0.279$ and under $p_2\left(w\right|x,y),$ we have $I_{p_2}(X;Y|W)<0.262$. By the same token, one can show that, under $(p_1\left(w\right|x,y)+p_2\left(w\right|x,y))/2,$ we have $I_{(p_1+p_2)/2}(X;Y|W)>0.274$. Hence, we conclude that, for this example,

$$\begin{array}{cc}\hfill I_{(p_1+p_2)/2}(X;Y|W)& >\frac{1}{2}\left(I_{p_1}(X;Y|W)+I_{p_2}(X;Y|W)\right),\hfill \end{array}$$

(12)

which proves that $I(X;Y|W)$ cannot be convex.

2.4. The Operational Significance of the Relaxed Wyner’s Common Information Problem

It is important to note that Wyner’s common information has clear and well-defined operational significance. This is perhaps not central to the detailed explanations and examples given in the sequel. However, it has a role in the appreciation of the rigorous connection established in our work. An excellent description of the operational significance is given in (Section I in [23]), where two separate aspects are identified. The first concerns a source coding scenario with a single encoder and two decoders, one interested in X and the other in $Y.$ Three bit streams are constructed by the encoder: One public (to both decoders), and two private streams, one for each decoder. Then, $C(X;Y)$ characterizes the minimum number of bits that must be sent on the public stream such that the total number of bits sent stays at the global minimum, which is well known to be the joint entropy of X and $Y.$ If the rate on the public bit stream dips below $C(X;Y),$ it is no longer possible to keep the total rate at the joint entropy. Rather, there is a strict penalty now, and this penalty can be expressed via $C_{\gamma }(X;Y).$ The second rigorous operational significance concerns the distributed generation of correlated randomness. We have two separate processors, one generating X and the other $Y.$ For a fixed desired resulting probability distribution $p(x,y),$ how many common random bits (shared between both processors) are required? Again, the answer is precisely $C(X;Y).$ A connection between caching and the Gray–Wyner network is developed in [30].

3. The Algorithm

The main technical result of this paper is to establish that the outcome of a specific procedure induced by the relaxed Wyner’s common information is tantamount to CCA whenever the original underlying distribution is Gaussian. In preparation for this, in this section, we present the proposed algorithm. In doing so, we will assume that the distribution of the data are $p(\mathbf{x},\mathbf{y}).$ In many applications involving CCA, the data distributions may not be known, but, rather, a number of samples of $\mathbf{X}$ and $\mathbf{Y}$ are provided, based on which CCA would then estimate the covariance matrix. A similar perspective can be taken on our procedure, but is left for future work. A short discussion can be found in Section 8 below.

3.1. High-Level Description

The proposed algorithm takes as input the distribution $p(\mathbf{x},\mathbf{y})$ of the data, as well as a level $\gamma .$ The level $\gamma$ is a non-negative real number and may be thought of as a resolution level or a measure of coarseness: If $\gamma =0$, then the full commonality (or common information) between $\mathbf{X}$ and $\mathbf{Y}$ is extracted in the sense that it is conditioned on the common information, $\mathbf{X}$ and $\mathbf{Y}$ are conditionally independent. Conversely, if $\gamma$ is large, then only the most important part of the commonality is extracted. Fixing the level $\gamma ,$ the idea of the proposed algorithm is to evaluate the relaxed Wyner’s Common Information of Equation (5) between the information sources (data sets) at the chosen level $\gamma .$ This evaluation will come with an associated conditional distribution $p_{\gamma }\left(w\right|x,y),$ namely, the conditional distribution attaining the minimum in the optimization problem of Equation (5). The second half of the proposed algorithm consists of leveraging the minimizing $p_{\gamma }\left(w\right|x,y)$ in such a way as to separately reduce $\mathbf{X}$ and $\mathbf{Y}$ to those features that best express the commonality. This may be thought of as a type of projection of the minimizing random variable W back onto $\mathbf{X}$ and $\mathbf{Y},$ respectively. For the case of Gaussian statistics, this can be made precise.

3.2. Main Steps of the Algorithm

The algorithm proposed here starts from the joint distribution of the data, $p(\mathbf{x},\mathbf{y}).$ Estimates of this distribution can be obtained from data samples ${\mathbf{X}}^n$ and ${\mathbf{Y}}^n$ via standard techniques. The main steps of the procedure can then be described as follows:

Generic Procedure 1 

(CICA).

1. 

Select a real number $\gamma ,$ where $0\le \gamma \le I(\mathbf{X};\mathbf{Y}).$ This is the compression level: A low value of γ represents low compression, and, thus, many components are retained. A high value of γ represents high compression, and, thus, only a small number of components are retained.

2. 

Solve the relaxed Wyner’s common information problem,

$$\begin{array}{c}\hfill \underset{p\left(w\right|\mathbf{x},\mathbf{y})}{min}I(\mathbf{X},\mathbf{Y};W) such that I(\mathbf{X};\mathbf{Y}|W)\le \gamma ,\end{array}$$

(13)

leading to an associated conditional distribution $p_{\gamma }\left(w\right|\mathbf{x},\mathbf{y}).$

3. 

Using the conditional distribution $p_{\gamma }\left(w\right|\mathbf{x},\mathbf{y})$ found in Step 2), the dimension-reduced data sets can now be found via one of the following three variants:

(a) 

Version 1: MAP (maximum a posteriori):

$$\begin{array}{cc}\hfill u\left(\mathbf{x}\right)& =arg\underset{w}{max}{p}_{\gamma }\left(w\right|\mathbf{x}),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill v\left(\mathbf{y}\right)& =arg\underset{w}{max}{p}_{\gamma }\left(w\right|\mathbf{y}).\hfill \end{array}$$

(15)

(b) 

Version 2: Conditional Expectation:

$$\begin{array}{cc}\hfill u\left(\mathbf{x}\right)& =\mathbb{E}\left[W\right|\mathbf{X}=\mathbf{x}],\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill v\left(\mathbf{y}\right)& =\mathbb{E}\left[W\right|\mathbf{Y}=\mathbf{y}].\hfill \end{array}$$

(17)

(c) 

Version 3: Marginal Integration:

$$\begin{array}{cc}\hfill u\left(\mathbf{x}\right)& ={\int }_{\mathbf{y}}p\left(\mathbf{y}\right)\mathbb{E}\left[W\right|\mathbf{X}=\mathbf{x},\mathbf{Y}=\mathbf{y}]d\mathbf{y},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill v\left(\mathbf{y}\right)& ={\int }_{\mathbf{x}}p\left(\mathbf{x}\right)\mathbb{E}\left[W\right|\mathbf{X}=\mathbf{x},\mathbf{Y}=\mathbf{y}]d\mathbf{x}.\hfill \end{array}$$

(19)

The present paper focuses on the three versions given here because, for these three versions, we can establish Theorem 2, showing that, in the case of Gaussian statistics, all three versions lead exactly to CCA. Second, we note that, for concrete examples, it is often evident which of the versions is preferable. For example, in Section 5, we consider a binary example where the associated W in Step 2 of our algorithm is also binary. In this case, Version 1 will reduce the original binary vector $\mathbf{X}$ to a binary scalar, which is perhaps the most desirable outcome. By contrast, Versions 2 and 3 require an explicit embedding of the binary example in the reals, and will reduce the original binary vector $\mathbf{X}$ to a real-valued scalar, which might not be as insightful.

4. For Gaussian, CICA Is CCA

In this section, we consider the special case where $\mathbf{X}$ and $\mathbf{Y}$ are jointly Gaussian random vectors. Since the mean has no bearing on either CCA or Wyner’s common information, we will assume it to be zero in the sequel, without loss of generality. One key ingredient for this argument is a well-known change of basis, see, for example [2], which we will now introduce in detail. Note that the mean will not change any mutual information term, thus we assume it to be zero without a loss of generality. We first need to introduce notation for CCA. To this end, let us express the covariance matrices, as usual, in terms of their eigendecompositions as

$$\begin{array}{c}\hfill K_{\mathbf{X}}=Q_x\left(\begin{array}{cc}{\Lambda }_{r_X}& 0\\ 0& 0_{n-r_X}\end{array}\right)Q_x^T\end{array}$$

(20)

and

$$\begin{array}{c}\hfill K_{\mathbf{Y}}=Q_y\left(\begin{array}{cc}{\Lambda }_{r_Y}& 0\\ 0& 0_{n-r_Y}\end{array}\right)Q_y^T,\end{array}$$

(21)

where $r_X$ and $r_Y$ denote the rank of $K_{\mathbf{X}}$ and $K_{\mathbf{Y}},$ respectively. Starting from this, we define the matrices

$$\begin{array}{c}\hfill K_{\mathbf{X}}^{-1/2}=Q_x\left(\begin{array}{cc}{\Lambda }_{r_X}^{-1/2}& 0\\ 0& 0_{n-r_X}\end{array}\right)Q_x^T\end{array}$$

(22)

and

$$\begin{array}{c}\hfill K_{\mathbf{Y}}^{-1/2}=Q_y\left(\begin{array}{cc}{\Lambda }_{r_Y}^{-1/2}& 0\\ 0& 0_{n-r_Y}\end{array}\right)Q_y^T,\end{array}$$

(23)

where, for a diagonal matrix $\Lambda$ with strictly positive entries, ${\Lambda }_{r_Y}^{-1/2}$ denotes the diagonal matrix whose diagonal entries are the reciprocals of the square roots of the entries of the matrix $\Lambda .$ Using these matrices, the key step is to apply the change of basis

$$\begin{array}{cc}\hfill \widehat{\mathbf{X}}& =K_{\mathbf{X}}^{-1/2}\mathbf{X}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill \widehat{\mathbf{Y}}& =K_{\mathbf{Y}}^{-1/2}\mathbf{Y}.\hfill \end{array}$$

(25)

In the new coordinates, the covariance matrices of $\widehat{\mathbf{X}}$ and $\widehat{\mathbf{Y}},$ respectively, can be shown to be

$$\begin{array}{c}\hfill K_{\widehat{\mathbf{X}}}=\left(\begin{array}{cc}{I}_{r_X}& 0\\ 0& 0_{n-r_X}\end{array}\right)\end{array}$$

(26)

and

$$\begin{array}{c}\hfill K_{\widehat{\mathbf{Y}}}=\left(\begin{array}{cc}{I}_{r_Y}& 0\\ 0& 0_{n-r_Y}\end{array}\right).\end{array}$$

(27)

Moreover, we have

$$\begin{array}{c}\hfill K_{\widehat{\mathbf{X}}\widehat{\mathbf{Y}}}=K_{\mathbf{X}}^{-1/2}{K}_{\mathbf{X}\mathbf{Y}}{K}_{\mathbf{Y}}^{-1/2}.\end{array}$$

(28)

Let us denote the singular value decomposition of this matrix by

$$\begin{array}{cc}\hfill K_{\widehat{\mathbf{X}}\widehat{\mathbf{Y}}}& =U\Sigma V^H.\hfill \end{array}$$

(29)

where $\Sigma$ contains, on its diagonal, the ordered singular values of this matrix, denoted by ${\rho }_1\ge {\rho }_2\ge \dots \ge {\rho }_n.$ In addition, let us define

$$\begin{array}{cc}\hfill \tilde{\mathbf{X}}& =U^H\widehat{\mathbf{X}}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \tilde{\mathbf{Y}}& =V^H\widehat{\mathbf{Y}},\hfill \end{array}$$

(31)

which implies that $K_{\tilde{\mathbf{X}}}=K_{\widehat{\mathbf{X}}},$$K_{\tilde{\mathbf{Y}}}=K_{\widehat{\mathbf{Y}}},$ and $K_{\tilde{\mathbf{X}}\tilde{\mathbf{Y}}}=\Sigma$.

Next, we will leverage this change of basis to establish Wyner’s common information and its relaxation for the Gaussian vector case, and then to prove the connection between Generic Procedure 1 and CCA.

4.1. Wyner’s Common Information and Its Relaxation in the Gaussian Case

For the case where $\mathbf{X}$ and $\mathbf{Y}$ are jointly Gaussian random vectors, a full and explicit solution to the optimization problem of Equation (5) is found in [5]. To give some high-level intuition, the proof starts by mapping from $\mathbf{X}$ to $\tilde{\mathbf{X}}$ and from $\mathbf{Y}$ to $\tilde{\mathbf{Y}},$ as in Equations (30) and (31). This preserves all mutual information expressions as well as joint Gaussianity. Moreover, due to the structure of the covariance matrices of the vectors $\tilde{\mathbf{X}}$ and $\tilde{\mathbf{Y}},$ we have that ${\left\{({\tilde{X}}_i,{\tilde{Y}}_i)\right\}}_{i=1}^n$ are n independent pairs of Gaussian random variables. Thus, by the tensorization property (see Lemma 1), the vector problem can be reduced to n parallel scalar problems. The solution of the scalar problem is the main technical contribution of [5], and we refer to that paper for the detailed proof. The resulting formula can be expressed as in the following theorem.

Theorem 1

(from [5]). Let $\mathbf{X}$ and $\mathbf{Y}$ be jointly Gaussian random vectors of length n and covariance matrix $K_{(\mathbf{X},\mathbf{Y})}$. Then,

$$\begin{array}{cc}\hfill & C_{\gamma }(\mathbf{X};\mathbf{Y})=\underset{{\gamma }_i:{\sum }_{i=1}^n{\gamma }_i=\gamma }{min}\sum _{i=1}^n{C}_{{\gamma }_i}(X_i;Y_i),\hfill \end{array}$$

(32)

where

$$\begin{array}{c}\hfill C_{{\gamma }_i}(X_i;Y_i)=\frac{1}{2}{log}^{+}\frac{(1+{\rho }_i)(1-\sqrt{1-e^{-2{\gamma }_i}})}{(1-{\rho }_i)(1+\sqrt{1-e^{-2{\gamma }_i}})}\end{array}$$

(33)

and ${\rho }_i$ (for $i=1,\cdots ,n$) are the singular values of $K_{\mathbf{X}}^{-1/2}{K}_{\mathbf{X}\mathbf{Y}}{K}_{\mathbf{Y}}^{-1/2},$ where $K_{\mathbf{X}}^{-1/2}$ and $K_{\mathbf{Y}}^{-1/2}$ are defined to mean that only the positive eigenvalues are inverted.

As pointed out above, we refer to (Theorem 7 in [5]) for a rigorous proof of this theorem.

4.2. CICA in the Gaussian Case and the Exact Connection with CCA

In this section, we consider the proposed CICA algorithm in the special case where the data distribution is $p(\mathbf{x},\mathbf{y})$, a (multivariate) Gaussian distribution. We establish that, in this case, the classic CCA is a solution to all versions of the proposed CICA algorithm. In this sense, CICA is a strict generalization of CCA. CCA is briefly reviewed in Appendix B. Leveraging the matrices U and V defined via the singular value decomposition in Equation (29), CCA performs the dimensonality reduction

$$\begin{array}{cc}\hfill u\left(\mathbf{x}\right)& =U_k^H\widehat{\mathbf{x}}=U_k^H{K}_{\mathbf{X}}^{-1/2}\mathbf{x}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill v\left(\mathbf{y}\right)& =V_k^H\widehat{\mathbf{y}}=V_k^H{K}_{\mathbf{Y}}^{-1/2}\mathbf{y},\hfill \end{array}$$

(35)

where the matrix $U_k$ contains the first k columns of U (that is, the k left singular vectors corresponding to the largest singular values), and the matrix $V_k$ the respective right singular vectors. We refer to these as the “top k CCA components.”

Theorem 2.

Let $\mathbf{X}$ and $\mathbf{Y}$ be jointly Gaussian random vectors. Then:

1. 

The top k CCA components are a solution to all three versions of Generic Procedure 1.

2. 

The parameter γ controls the number k as follows:

$$\begin{array}{cc}\hfill k\left(\gamma \right)& =\left\{\begin{array}{cc}n,\hfill & if 0\le \gamma <ng\left({\rho }_n\right),\hfill \\ n-1,\hfill & if ng\left({\rho }_n\right)\le \gamma <(n-1)g\left({\rho }_{n-1}\right)+g\left({\rho }_n\right),\hfill \\ n-2,\hfill & if (n-1)g\left({\rho }_{n-1}\right)+g\left({\rho }_n\right)\le \gamma \hfill \\ \hfill & <(n-2)g\left({\rho }_{n-2}\right)+g\left({\rho }_{n-1}\right)+g\left({\rho }_n\right),\hfill \\ ⋮,\hfill & ⋮,\hfill \\ \ell \hfill & if (\ell +1)g\left({\rho }_{\ell +1}\right)+\sum _{i=\ell +2}^{n}g\left({\rho }_i\right)\le \gamma \hfill \\ \hfill & <\ell g\left({\rho }_{\ell }\right)+\sum _{i=\ell +1}^{n}g\left({\rho }_i\right),\hfill \\ ⋮,\hfill & ⋮,\hfill \\ 1,\hfill & if 2g\left({\rho }_2\right)+\sum _{i=2}^{n}g\left({\rho }_i\right)\le \gamma <\sum _{i=1}^{n}g\left({\rho }_i\right),\hfill \\ 0,\hfill & if \sum _{i=1}^{n}g\left({\rho }_i\right)\le \gamma ,\hfill \end{array}\right.\hfill \end{array}$$

(36)

where $g\left(\rho \right)=\frac{1}{2}log\frac{1}{1-{\rho }^2}.$

Remark 1.

Note that $k\left(\gamma \right)$ is a decreasing, integer-valued function. An illustration for a special case is given in Figure 2.

Proof. 

The main contribution of the theorem is the first item, i.e., the connection between CCA and Generic Procedure 1 in the case where $\mathbf{X}$ and $\mathbf{Y}$ are jointly Gaussian. The proof follows along the steps of the CICA procedure: We first show that, in Step 2, when $\mathbf{X}$ and $\mathbf{Y}$ are jointly Gaussian, then the minimizing W may be taken jointly Gaussian with $\mathbf{X}$ and $\mathbf{Y}.$ Then, we establish that, in Step 3, with the W from Step 2, we indeed obtain that the dimension-reduced representations $u\left(\mathbf{x}\right)$ and $v\left(\mathbf{y}\right)$ turn into the top k CCA components. In detail:

Step 2 of Generic Procedure 1: The technical heavy lifting for this step in the case where $p\left(\mathbf{x},\mathbf{y}\right)$ is a multivariate Gaussian distribution is presented in [5]. We shall briefly summarize it here. In the case of Gaussian vectors, the solution to the optimization problem in Equation (5) is most easily described in two steps. First, we apply the change of basis indicated in Equations (24) and (25). This is a one-to-one transform, leaving all information expressions in Equation (5) unchanged. In the new basis, we have n independent pairs. By the tensorization property (see Lemma 1), when $\mathbf{X}$ and $\mathbf{Y}$ consist of independent pairs, the solution to the optimization problem in Equation (5) can be reduced to n separate scalar optimizations. The remaining crux then is solving the scalar Gaussian version of the optimization problem in Equation (5). This is done in (Theorem 3 in [5]) via an argument of factorization of convex envelope. The full solution to the optimization problem is given in Equations (32) and (33). The remaining allocation problem over the non-negative numbers ${\gamma }_i$ can be shown to lead to a water-filling solution, given in (Theorem 8 in [5]). More explicitly, to understand this solution, start by setting $\gamma =I(\mathbf{X};\mathbf{Y}).$ Then, the corresponding $C_{\gamma }(\mathbf{X};\mathbf{Y})=0$ and the optimizing distribution $p_{\gamma }\left(w\right|\mathbf{x},\mathbf{y})$ trivializes. Now, as we lower $\gamma ,$ the various terms in the sum in Equation (32) start to become non-zero, starting with the term with the largest correlation coefficient ${\rho }_1.$ Hence, an optimizing distribution $p_{\gamma }\left(w\right|\mathbf{x},\mathbf{y})$ can be expressed as ${\mathbf{W}}_{\gamma }=U_k^H{K}_{\mathbf{X}}^{-1/2}\mathbf{X}+V_k^H{K}_{\mathbf{Y}}^{-1/2}\mathbf{Y}+\mathbf{Z},$ where the matrices $U_k$ and $V_k$ are precisely the top k CCA components (see Equations (34) and (35) and the following discussion), and $\mathbf{Z}$ is additive Gaussian noise with mean zero, independent of $\mathbf{X}$ and $\mathbf{Y}.$

Step 3 of Generic Procedure 1: For the algorithm, we need the corresponding conditional marginals, $p_{\gamma }\left(w\right|\mathbf{x})$ and $p_{\gamma }\left(w\right|\mathbf{y}).$ By symmetry, it suffices to prove one formula. Changing basis as in Equations (24) and (25), we can write

$$\begin{array}{cc}\hfill \mathbb{E}\left[W\right|\mathbf{X}]& =\mathbb{E}[U_k^H\widehat{\mathbf{X}}+V_k^H\widehat{\mathbf{Y}}+\mathbf{Z}|\widehat{\mathbf{X}}]\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & =U_k^H\widehat{\mathbf{X}}+V_k^H\mathbb{E}\left[\widehat{\mathbf{Y}}\right|\widehat{\mathbf{X}}]\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill & =U_k^H\widehat{\mathbf{X}}+V_k^H\left(\mathbb{E}\left[\widehat{\mathbf{Y}}{\widehat{\mathbf{X}}}^H\right]{\left(\mathbb{E}\left[\widehat{\mathbf{X}}{\widehat{\mathbf{X}}}^H\right]\right)}^{-1}\widehat{\mathbf{X}}\right)\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & =U_k^H\widehat{\mathbf{X}}+V_k^H{K}_{\widehat{\mathbf{Y}}\widehat{\mathbf{X}}}\widehat{\mathbf{X}}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & =U_k^H\widehat{\mathbf{X}}+{\left(K_{\widehat{\mathbf{X}}\widehat{\mathbf{Y}}}{V}_k\right)}^H\widehat{\mathbf{X}}.\hfill \end{array}$$

(41)

The first summand contains exactly the top k CCA components extracted from $\mathbf{X},$ which is the claimed result. The second summand requires further scrutiny. To proceed, we observe that, for CCA, the projection vectors obey the relationship (see Equation (A12))

$$\begin{array}{cc}\hfill \mathbf{u}& =\alpha K_{\widehat{\mathbf{X}}\widehat{\mathbf{Y}}}\mathbf{v},\hfill \end{array}$$

(42)

for some real-valued constant $\alpha .$ Thus, combining the top k CCA components, we can write

$$\begin{array}{cc}\hfill U_k& =D{K}_{\widehat{\mathbf{X}}\widehat{\mathbf{Y}}}{V}_k,\hfill \end{array}$$

(43)

where D is a diagonal matrix. Hence,

$$\begin{array}{cc}\hfill \mathbb{E}\left[W\right|\mathbf{X}]& =U_k^H\widehat{\mathbf{X}}+D^{-1}{U}_k^H\widehat{\mathbf{X}}\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & =\tilde{D}{U}_k^H\widehat{\mathbf{X}},\hfill \end{array}$$

(45)

where $\tilde{D}$ is the diagonal matrix

$$\begin{array}{cc}\hfill \tilde{D}& =I+D^{-1}.\hfill \end{array}$$

(46)

This is precisely the top k CCA components (note that the solution to the CCA problem (A7) is only specified up to a scaling). This establishes the theorem for the case of Version 2) of the proposed algorithm. Clearly, it also establishes that $p_{\gamma }\left(w\right|\mathbf{x})$ is a Gaussian distribution with mean given by (45), thus establishing the theorem for Version 1) of the proposed algorithm. The proof for Version 3 follows along similar lines and is thus omitted. □

5. A Binary Example

In this section, we carry through a theoretical study of a somewhat more general case of the example discussed in Section 1.4 that is believed to be within the reach of practical data. In order to do a theoretical study, we need to constrain the data into binary for the reason that computing the Wyner’s common information for doubly binary symmetric source is already known.

Let us illustrate the proposed algorithm via a simple example. Consider the vector $(U,X_2,V,Y_2)$ of binary random variables. Suppose that $(U,V)$ is a doubly symmetric binary source. This means that U is uniform and V is the result of passing U through a binary symmetric (“bit-flipping”) channel with flip probability denoted by $a_0$ to match the notation in (Section 3 in [23]). Without loss of generality, we may assume $a_0\le \frac{1}{2}.$ Meanwhile, $X_2$ and $Y_2$ are independent binary uniform random variables, also independent of the pair $(U,V).$ We will then form the vectors $\mathbf{X}$ and $\mathbf{Y}$ as

$$\begin{array}{cc}\hfill \mathbf{X}& =\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right)=\left(\begin{array}{c}U\oplus X_2\\ X_2\end{array}\right),\hfill \end{array}$$

(47)

and

$$\begin{array}{cc}\hfill \mathbf{Y}& =\left(\begin{array}{c}{Y}_1\\ Y_2\end{array}\right)=\left(\begin{array}{c}V\oplus Y_2\\ Y_2\end{array}\right),\hfill \end{array}$$

(48)

where ⊕ denotes the modulo-reduced addition, as usual. How do various techniques perform for this example?

• Let us first analyze the behavior and outcome of CCA in this particular example. The key observation is that any pair, amongst the four entries in these two vectors, $X_1,X_2,Y_1,$ and $Y_2,$ are (pairwise) independent binary uniform random variables. Hence, the overall covariance matrix of the merged random vector ${({\mathbf{X}}^T,{\mathbf{Y}}^T)}^T$ is merely a scaled identity matrix. This, in turn, implies that CCA as described in Equations (34) and (35) merely boils down to the identity mapping. Concretely, this means that, for CCA, in this example, the best one-dimensional projections are ex aequo any pair of one coordinate of the vector $\mathbf{X}$ with one coordinate of the vector $\mathbf{Y}.$ As we have already explain above, any such pair is merely a pair of independent (and identically distributed) random variables, so CCA does not extract any dependence between $\mathbf{X}$ and $\mathbf{Y}$ at all. Of course, this is the main point of the present example.
• How does CICA perform in this example? We selected this example because it represents one of the only cases for which a closed-form solution to the optimization problem in Equation (13) is known, at least in the case $\gamma =0.$ To see this, let us first observe that, in our example, we have

  $$\begin{array}{c}\hfill p(u,v,x_2,y_2)=p(u,v)p\left(x_2\right)p\left(y_2\right).\end{array}$$

  (49)

  Next, we observe that

  $$\begin{array}{cc}\hfill C_{\gamma }(\mathbf{X};\mathbf{Y})& =C_{\gamma }(U,X_2;V,Y_2)\hfill \end{array}$$

  (50)

  $$\begin{array}{cc}\hfill & =C_{\gamma }(U;V,Y_2)\hfill \end{array}$$

  (51)

  $$\begin{array}{cc}\hfill & =C_{\gamma }(U;V)\hfill \end{array}$$

  (52)

  where (51) follows from Lemma 1, Item 5, and the Markov chain $X_2-U-(V,Y_2)$ that is satisfied from (49). The last Equation (52) follows from Lemma 1, Item 5, and the Markov chain $Y_2-V-U$ that is satisfied from (49). That is, in this simple example, solving the optimization problem of Equation (13) is tantamount to solving the optimization problem in Equation (52). For the latter, the solution is well known, see (Section 3 in [23]). Specifically, we can express the conditional distribution $p_{\gamma }\left(w\right|\mathbf{x},\mathbf{y})$ that solves the optimization problem of Equation (13) and is required for Step 3 of Generic Procedure 1 as follows:

  $$\begin{array}{cc}\hfill p_{\gamma =0}\left(w\right|\mathbf{x},\mathbf{y})& =\left\{\begin{array}{cc}1-\nu ,& \mathrm{if} w=0,x_1\oplus x_2=0,y_1\oplus y_2=0,\hfill \\ \nu ,& \mathrm{if} w=1,x_1\oplus x_2=0,y_1\oplus y_2=0,\hfill \\ \nu ,& \mathrm{if} w=0,x_1\oplus x_2=1,y_1\oplus y_2=1,\hfill \\ 1-\nu ,& \mathrm{if} w=1,x_1\oplus x_2=1,y_1\oplus y_2=1,\hfill \\ \frac{1}{2},& \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

  (53)

  where

  $$\begin{array}{cc}\hfill \nu & =\frac{1}{2}-\frac{\sqrt{1-2a_0}}{2(1-a_0)}.\hfill \end{array}$$

  (54)
  Let us now apply Version 1 (the MAP version) of Generic Procedure 1. To this end, we also need to calculate $p_{\gamma }\left(w\right|\mathbf{x})$ and $p_{\gamma }\left(w\right|\mathbf{y}).$ Again, for $\gamma =0,$ these can be expressed in a closed form as follows:

  $$\begin{array}{cc}\hfill p_{\gamma =0}\left(w\right|\mathbf{x})& =\left\{\begin{array}{cc}1-a_1,& \mathrm{if} w=0,x_1\oplus x_2=0,\hfill \\ a_1,& \mathrm{if} w=1,x_1\oplus x_2=0,\hfill \\ a_1,& \mathrm{if} w=0,x_1\oplus x_2=1,\hfill \\ 1-a_1,& \mathrm{if} w=1,x_1\oplus x_2=1,\hfill \end{array}\right.\hfill \end{array}$$

  (55)

  where

  $$\begin{array}{cc}\hfill a_1& =\frac{1}{2}\left(1-\sqrt{1-2a_0}\right).\hfill \end{array}$$

  (56)

  The formula for $p_{\gamma }\left(w\right|\mathbf{y})$ follows by symmetry and shall be omitted. The final step is to follow Equations (14) and (15) and find $arg{max}_w{p}_{\gamma =0}\left(w\right|\mathbf{x})$ for each $\mathbf{x}$ as well as $arg{max}_w{p}_{\gamma =0}\left(w\right|\mathbf{y})$ for each $\mathbf{y}.$ For the example at hand, these can be compactly expressed as

  $$\begin{array}{cc}\hfill u\left(\mathbf{x}\right)& =arg\underset{w}{max}{p}_{\gamma }\left(w\right|\mathbf{x})=x_1\oplus x_2=u,\hfill \end{array}$$

  (57)

  $$\begin{array}{cc}\hfill v\left(\mathbf{y}\right)& =arg\underset{w}{max}{p}_{\gamma }\left(w\right|\mathbf{y})=y_1\oplus y_2=v,\hfill \end{array}$$

  (58)

  from the fact that $0\le a_0\le \frac{1}{2}$ that implies $0\le a_1\le \frac{1}{2}$. Hence, we find that, for CICA as described in Generic Procedure 1, an optimal solution is to reduce $\mathbf{X}$ to U and $\mathbf{Y}$ to $V.$ This captures all the dependence between the vectors $\mathbf{X}$ and $\mathbf{Y},$ which appears to be the most desirable outcome.
  As a final note, we point out that it is conceptually straightforward to evaluate Versions 2 and 3 (conditional expectation) of Generic Procedure 1 in this example, but this would require embedding the considered binary alphabets into the real numbers. This makes it a less satisfying option for the simple example at hand.

6. A Gradient Descent Based Implementation

As we discussed above, in our problem, the objective is indeed a convex function of the optimization variables (but the constraint set is not convex). Clearly, this gives hope that gradient-based techniques may lead to interesting solutions. In this section, we examine a first tentative implementation and check it against ground truth for some simple examples.

In theory for convex problems, gradient descent will guarantee convergence to the optimal solution; otherwise, it will guarantee only local convergence. Gradient descent runs in iterative steps, where each step does a local linear approximation and the step size depends on a learning parameter that is $\alpha$ for our problem. In our work, we want to minimize the objective $I(W;X,Y)$ when the constraint $I(X;Y|W)$ is held below a $\gamma -$level.

Instead, we apply a variant of gradient descent where we minimize the weighted sum of objective $I(W;X,Y)$ and the constraint $I(X;Y|W)$, which is $I(W;X,Y)+\lambda I(X;Y|W)$. The parameter $\lambda$ will permit some control on the constraint, thus sweeping all its possible values. We present the algorithm where $C\left(p\right(w|x,y))$ will be a function of $p\left(w\right|x,y)$ that will represent $I(W;X,Y)$, and $J\left(p\right(w|x,y))$ will be a function of $p\left(w\right|x,y)$ that will represent $I(X;Y|W)$.

The exact computation of the stated update step is presented in the following lemma.

Lemma 2

(Computation of the update step). Let $p(x,y)$ be a fixed distribution, then the updating steps for the gradient descent are

$$\begin{array}{cc}\hfill \frac{\partial C\left(p\right(w|x,y)}{\partial p\left(w\right|x,y)}& =p(x,y)log\frac{p\left(w\right|x,y)}{{\sum }_{x^{\prime },y^{\prime }}p(x^{\prime },y^{\prime })p\left(w\right|x^{\prime },y^{\prime })},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill \frac{\partial J\left(p\right(w|x,y)}{\partial p\left(w\right|x,y)}& =p(x,y)log\frac{p\left(w\right|x,y){\sum }_{x^{\prime },y^{\prime }}p(x^{\prime },y^{\prime })p\left(w\right|x^{\prime },y^{\prime })}{{\sum }_{x^{″}}p\left(w\right|x^{″},y)p\left(x^{″}\right|y){\sum }_{y^{″}}p\left(w\right|x,y^{″})p\left(y^{″}\right|x)}.\hfill \end{array}$$

(60)

Proof. 

Let the function C be as defined above

$$\begin{array}{cc}\hfill C\left(p\right(w|x,y))& =\sum _{x,y,w}p\left(w\right|x,y)p(x,y)log\frac{p\left(w\right|xy)}{{\sum }_{x^{\prime },y^{\prime }}p\left(w\right|x^{\prime },y^{\prime })p(x^{\prime },y^{\prime })},\hfill \end{array}$$

(61)

and, in terms of information theoretic terms, the function is $C\left(p\right(w|x,y))=I(W;X,Y)$. In addition, $C\left(p\right(w|x,y))$ is a convex function of $p\left(w\right|x,y)$, shown in (Theorem 2.7.4 in [29]). Taking the first derivative, we get

$$\begin{array}{cc}\hfill \frac{\partial C\left(p\right(w|x,y))}{\partial p\left(w\right|x,y)}& =p(x,y)log\frac{p\left(w\right|x,y)}{{\sum }_{x^{\prime },y^{\prime }}p\left(w\right|x^{\prime },y^{\prime })p(x^{\prime },y^{\prime })}+p\left(w\right|x,y)p(x,y)\frac{1}{p\left(w\right|x,y)}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & -\sum _{x^{″},y^{″}}p\left(w\right|x^{″},y^{″})p(x^{″},y^{″})\frac{p(x,y)}{{\sum }_{x^{\prime },y^{\prime }}p\left(w\right|x^{\prime },y^{\prime })p(x^{\prime },y^{\prime })}\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill & =p(x,y)log\frac{p\left(w\right|x,y)}{{\sum }_{x^{\prime },y^{\prime }}p\left(w\right|x^{\prime },y^{\prime })p(x^{\prime },y^{\prime })}.\hfill \end{array}$$

(63)

On the other hand, the term $I(X;Y|W)$ can be expressed as

$$\begin{array}{cc}\hfill I(X;Y|W)& =I(W;X,Y)-I(W;X)-I(W;Y)+I(X;Y)\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill & =C\left(p\right(w|x,y))-C(p\left(w\right|x\left)\right)-C\left(p\right(w\left|y\right))+I(X;Y).\hfill \end{array}$$

(65)

Taking the derivative with respect to $p\left(w\right|x,y)$ becomes easier once $I(X;Y|W)$ is written in terms of function C and we already know the derivative of C from (63). Thus, the derivative would be

$$\begin{array}{cc}\hfill \frac{\partial J\left(p\right(w|x,y)}{\partial p\left(w\right|x,y)}& =\frac{\partial C\left(p\right(w|x,y))}{\partial p\left(w\right|x,y)}-\frac{\partial C\left(p\right(w\left|x\right))}{\partial p\left(w\right|x,y)}-\frac{\partial C\left(p\right(w\left|y\right))}{\partial p\left(w\right|x,y)}\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill & =\frac{\partial C\left(p\right(w|x,y))}{\partial p\left(w\right|x,y)}-\frac{\partial C\left(p\right(w\left|x\right))}{\partial p\left(w\right|x)}\frac{\partial p\left(w\right|x)}{\partial p\left(w\right|x,y)}-\frac{\partial C\left(p\right(w\left|y\right))}{\partial p\left(w\right|y)}\frac{\partial p\left(w\right|y)}{\partial p\left(w\right|x,y)}\hfill \\ \hfill & =p(x,y)log\frac{p\left(w\right|x,y)}{{\sum }_{x^{\prime },y^{\prime }}p\left(w\right|x^{\prime },y^{\prime })p(x^{\prime },y^{\prime })}-p\left(x\right)log\frac{p\left(w\right|x)}{{\sum }_{x^{″}}p\left(w\right|x^{″})p\left(x^{″}\right)}p\left(y\right|x)\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill & -p\left(y\right)log\frac{p\left(w\right|y)}{{\sum }_{y^{″}}p\left(w\right|y^{″})p\left(y^{″}\right)}p\left(x\right|y)\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill & =p(x,y)log\frac{p\left(w\right|x,y){\sum }_{x^{\prime },y^{\prime }}p(x^{\prime },y^{\prime })p\left(w\right|x^{\prime },y^{\prime })}{{\sum }_{x^{″}}p\left(w\right|x^{″},y)p\left(x^{″}\right|y){\sum }_{y^{″}}p\left(w\right|x,y^{″})p\left(y^{″}\right|x)}.\hfill \end{array}$$

(69)

where (67) is an application of the chain rule, and the rest is straightforward computation. □

Remark 2.

In practice, it is useful and computationally cheaper to replace the derivative formulas in Lemma 2 by their standard approximations. That is, the updating step in line 7 of Algorithm 1 is replaced by

$$\begin{array}{cc}\hfill & \frac{\partial C\left(p\right(w|x,y))}{\partial p\left(w\right|x,y)}\approx \frac{C\left(p\right(w|x,y)+\Delta )-C(p\left(w\right|x,y\left)\right)}{\Delta },\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill & \frac{\partial J\left(p\right(w|x,y))}{\partial p\left(w\right|x,y)}\approx \frac{J\left(p\right(w|x,y)+\Delta )-J(p\left(w\right|x,y\left)\right)}{\Delta },\hfill \end{array}$$

(71)

for a judicious choice of $\Delta .$ This is the version that was used for Figure 1b, with $\Delta ={10}^{-3}$. We point out that, in the general case, the error introduced by this approximation is not bounded.

Algorithm 1: Approximate CICA Algorithm via Gradient Descent

7. Extension to More than Two Sources

It is unclear how one would extend CCA to more than two databases. By contrast, for CICA, this extension is conceptually straightforward. For Wyner’s common information, in Definition 1, it suffices to replace the objective in the minimization by $I(X_1,X_2,\dots ,X_M;W)$ and to keep the constraint of conditional independence. To obtain an interesting algorithm, we now need to relax the constraint of conditional independence. The most natural way is via the conditional version of Watanabe’s total correlation [31], leading to the following definition:

Definition 3

(Relaxed Wyner’s Common Information for M variables). For a fixed probability distribution $p(x_1,x_2,\dots ,x_M),$ we define

$$\begin{array}{cc}\hfill C_{\gamma }(X_1;X_2;\dots ;X_M)& =infI(X_1,X_2,\dots ,X_M;W)\hfill \end{array}$$

(72)

such that ${\sum }_{i=1}^{M}H\left(X_i\right|W)-H(X_1,X_2,\dots ,X_M|W)\le \gamma ,$ where the infimum is over all probability distributions $p(w,x_1,x_2,\dots ,x_M)$ with marginal $p(x_1,x_2,\dots ,x_M).$

Not surprisingly, an explicit closed-form solution is difficult to find. One simple case appears below as part of the example presented in Section 7.1, see Lemma 4. By analogy with Lemma 1, we can again state basic properties.

Lemma 3.

$C_{\gamma }(X_1;X_2;\dots ;X_M)$ satisfies the following basic properties:

1. 

$C_{\gamma }(X_1;X_2;\dots ;X_M)\ge \frac{1}{M-1}max\{{\sum }_{i=1}^{M}H\left(X_i\right)-H(X_1,X_2,\dots ,X_M)-\gamma ,0\}.$

2. 

$C_{\gamma }(X_1;X_2;\dots ;X_M)$ is a convex and continuous function of γ for $\gamma \ge 0.$

3. 

If $Z-X_1-(X_2,\dots ,X_M)$ forms a Markov chain,

then $C_{\gamma }((X_1,Z);X_2;\dots ;X_M)=C_{\gamma }(X_1;X_2;\dots ;X_M).$

4. 

The cardinality of $\mathcal{W}$ may be restricted to $\left|\mathcal{W}\right|\le {\prod }_{i=1}^M\left|{\mathcal{X}}_i\right|+1.$

5. 

If $f_i(·)$ are one-to-one functions,

then $C_{\gamma }(f_1\left(X_1\right);f_2\left(X_2\right);\dots ;f_M\left(X_M\right))=C_{\gamma }(X_1;X_2;\dots ;X_M).$

6. 

For discrete X, we have $C_{\gamma }(X;X;\dots ;X)=max\{H\left(X\right)-\frac{\gamma }{M-1},0\}.$

Proofs for these basic properties can be found in Appendix C.

Leveraging Definition 3, it is conceptually straightforward to extend CICA (that is, Generic Procedure 1) to the case of M databases as follows. For completeness, we include an explicit statement of the resulting procedure.

Generic Procedure 2 

(CICA with multiple sources).

1. 

Select a real number $\gamma ,$ where $0\le \gamma \le {\sum }_{i=1}^{M}H\left({\mathbf{X}}_i\right)-H({\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_M).$ This is the compression level: A low value of γ represents low compression, and, thus, many components are retained. A high value of γ represents high compression, and, thus, only a small number of components are retained.

2. 

Solving the relaxed Wyner’s common information problem, -4.6cm0cm

$$\begin{array}{c}\hfill \underset{p\left(w\right|{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_M)}{min}I({\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_M;W) such that \sum _{i=1}^{M}H\left({\mathbf{X}}_i\right|W)-H({\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_M|W)\le \gamma ,\end{array}$$

(73)

leading to an associated conditional distribution $p_{\gamma }\left(w\right|{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_M).$

3. 

Using the conditional distribution $p_{\gamma }\left(w\right|{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_M)$ found in Step 2, the dimension-reduced data sets can now be found via one of the following three variants:

(a) 

Version 1: MAP (maximum a posteriori):

$$\begin{array}{cc}\hfill u_i\left({\mathbf{x}}_i\right)& =arg\underset{w}{max}{p}_{\gamma }\left(w\right|{\mathbf{x}}_i),\hfill \end{array}$$

(74)

for $i=1,2,\dots ,M.$

(b) 

Version 2: Conditional Expectation:

$$\begin{array}{cc}\hfill u_i\left({\mathbf{x}}_i\right)& =\mathbb{E}\left[W\right|{\mathbf{X}}_i={\mathbf{x}}_i]\hfill \end{array}$$

(75)

for $i=1,2,\dots ,M.$

(c) 

Version 3: Marginal Integration:

$$\begin{array}{cc}\hfill u_i\left({\mathbf{x}}_i\right)& ={\int }_{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_{i-1},{\mathbf{x}}_{i+1},\dots ,{\mathbf{x}}_M}p({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{i-1},{\mathbf{x}}_{i+1},\dots ,{\mathbf{x}}_M)\hfill \\ \hfill & \mathbb{E}\left[W\right|{\mathbf{X}}_1={\mathbf{x}}_1,\dots ,{\mathbf{X}}_M={\mathbf{x}}_M]d{\mathbf{x}}_1\cdots d{\mathbf{x}}_{i-1}d{\mathbf{x}}_{i+1}\cdots d{\mathbf{x}}_M\hfill \end{array}$$

(76)

for $i=1,2,\dots ,M.$

Clearly, Generic Procedure 2 closely mirrors Generic Procedure 1. The key difference is that there is no direct analog of Theorem 2. This is no surprise since it is unclear how CCA would be extended to beyond the case of two sources. Nonetheless, it would be very interesting to explore what Generic Procedure 2 boils down to in the special case when all vectors are jointly Gaussian. At the current time, this is unknown. In fact, the explicit solution to the optimization problem in Definition 3 is presently an open problem.

Instead, we illustrate the promise of Generic Procedure 2 via a simple binary example in the next section. The example mirrors some of the basic properties of the example tackled in Section 5.

7.1. A Binary Example with Three Sources

In this section, we develop an example with three sources that borrows some of the ideas from the example discussed in Section 5. In a sense, the present example is even more illustrative because, in it, any two of the original vectors ${\mathbf{X}}_1,{\mathbf{X}}_2,$ and ${\mathbf{X}}_3,$ are (pairwise) independent. Therefore, any method based on pairwise measures, including CCA and maximal correlation, would not identify any commonality at all. Specifically, we consider the following simple statistical model:

$$\begin{array}{c}\hfill {\mathbf{X}}_1=\left(\begin{array}{c}U\\ Z_1\end{array}\right),{\mathbf{X}}_2=\left(\begin{array}{c}V\\ Z_2\end{array}\right),{\mathbf{X}}_3=\left(\begin{array}{c}U\oplus V\\ Z_3\end{array}\right),\end{array}$$

(77)

where $U,V,Z_1,Z_2,Z_3$ are independent uniform binary variables and ⊕ denotes modulo-2 addition. We observe that, amongst these three vectors, any pair is independent. This implies, for example, that any correlation-based technique (including maximal correlation) will not identify any relevant features, since correlation is a pairwise measure. By contrast, we can show that one output of Algorithm 2 is indeed to select $W=(U,V),$ for $\gamma =0.$ Thus, the algorithm would reduce each of the three vectors to their first component, which is the intuitively pleasing answer in this case. By going through the steps of the Generic Procedure 2, for $\gamma =0$, where the the joint distribution satisfies

$$\begin{array}{c}\hfill p(u,v,u\oplus v,z_1,z_2,z_3)=p(u,v,u\oplus v)p\left(z_1\right)p\left(z_2\right)p\left(z_3\right)\end{array}$$

(78)

we have that

$$\begin{array}{cc}\hfill C({\mathbf{X}}_1;{\mathbf{X}}_2;{\mathbf{X}}_3)& =C(U,Z_1;V,Z_2;U\oplus V,Z_3)\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill & =C(U;V,Z_2;U\oplus V,Z_3)\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill & =C(U;V;U\oplus V,Z_3)\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill & =C(U;V;U\oplus V)\hfill \end{array}$$

(82)

where we use Lemma 3, Item 3, together with the Markov chain $Z_1-U-(Z_2,V,Z_3,U\oplus V)$ that follows from (78) to prove step (80). Similarly, the Markov chain $Z_2-V-(U,Z_3,U\oplus V)$ proves step (81) by making use of Lemma 3, Item 3. A similar argument is used for the last step (82). Managing to compute $C(U;V;U\oplus V)$ is equivalent to computing $C({\mathbf{X}}_1;{\mathbf{X}}_2;{\mathbf{X}}_3)$, and we demonstrate how to compute it in the next part.

Lemma 4.

Let $U,V$ be independent uniform binary variables and ⊕ denotes modulo-2 addition. Then, the optimal solution to

$$\begin{array}{c}\hfill C_{\gamma =0}(U;V;U\oplus V)=\underset{W:H\left(U\right|W)+H(V\left|W\right)+H(U\oplus V|W)-H(U,V,U\oplus V\left|W\right)=0}{inf}I(W;U,V,U\oplus V)\end{array}$$

(83)

is $W=(U,V)$, where the expression evaluates to two.

The proof is given in Appendix D. If we apply Version 1 of Step 3 of Generic Procedure 2, we obtain

$$\begin{array}{c}\hfill arg\underset{w}{max}{p}_{\gamma =0}\left(w\right|{\mathbf{x}}_1)=\{(u,0),(u,1)\},\end{array}$$

(84)

that is, in this case, the maximizer is not unique. However, as we observe that the set of maximizers is a deterministic function of u alone, it is natural to reduce as follows:

$$\begin{array}{cc}\hfill u_1\left({\mathbf{x}}_1\right)& =u.\hfill \end{array}$$

(85)

By the same token, we can reduce

$$\begin{array}{cc}\hfill u_2\left({\mathbf{x}}_2\right)& =v,\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill u_3\left({\mathbf{x}}_3\right)& =u\oplus v.\hfill \end{array}$$

(87)

In this example, it is clear that this indeed extracts all of the dependency there is between our three sources, and, thus, is the correct answer.

As pointed out above, in this simple example, any pair of the random vectors ${\mathbf{X}}_1,{\mathbf{X}}_2,$ and ${\mathbf{X}}_3$ are (pairwise) independent, which implies that the classic tools based on pairwise measures (CCA, maximal correlation) cannot identify any commonality between ${\mathbf{X}}_1,{\mathbf{X}}_2,$ and ${\mathbf{X}}_3.$

8. Conclusions and Future Work

We introduce a novel two-step procedure that we refer to as CICA. The first step consists of an information minimization problem related to Wyner’s common information, while the second can be thought of as a type of back-projection. We prove that, in the special case of Gaussian statistics, this two-step procedure precisely extracts the CCA components. A free parameter $\gamma$ in CICA permits selecting the number of CCA components that are being extracted. In this sense, the paper establishes a novel rigorous connection between CCA and information measures. A number of simple examples are presented. It is also shown how to extend the novel algorithm to more than two sources.

Future work includes a more in-depth study and consideration to assess the practical promise of this novel algorithm. This will also require moving beyond the current setting where it was assumed that the probability distribution of the data at hand was provided directly. Instead, this distribution has to be estimated from data, and one needs to understand what limitations this additional constraint will end up imposing.