implement by a linear layer. Then, any $f\in\mathcal{F}_{\mathcal{X}}^{\,k}$
can be expressed in this way by setting corresponding weights in the linear
layer. Similarly, we can express $g\in\mathcal{F}_{\mathcal{Y}}^{\,k}$ using
another linear layer.

(a) $f_{i}^{*}(x)$ for each $i=1,2,3$ (b) $g_{i}^{*}(y)$ for each $i=1,2,3$

(c) $\sigma_{i}$

Figure 13: Top three dependence modes learned from a discrete dataset, in
consistent with theoretical results.

In the experiment, we set $|\mathcal{X}|=8$, $|\mathcal{Y}|=6$, and randomly
generate a $P_{X,Y}\in\mathcal{P}^{\mathcal{X}\times\mathcal{Y}}$. We generate
$N=30\,000$ training samples from $P_{X,Y}$, and learn $k=3$ dimensional
features $f,g$ by maximizing
$\mathscr{H}(f,g;\\{(1)^{k};\,\mathcal{F}_{\mathcal{X}};\,\mathcal{F}_{\mathcal{Y}}\\})$.
Then, we normalize each $f_{i},g_{i}$ to obtain corresponding estimations of
$f_{i}^{*}$, $g_{i}^{*}$, and $\sigma_{i}$ by applying (37). We show the
estimated features and singular values in Figure 13, which are consistent with
the corresponding theoretical values computed from $P_{X,Y}$.

####  7.1.2 Continuous Data

We proceed to consider a continuous dataset with degenerate dependence modes,
i.e., the singular values $\sigma_{i}$’s are not all distinct.

In particular, we consider $X,Y$ taking values from
$\mathcal{X}=\mathcal{Y}=[-1,1]$, where the joint probability density function
$p_{X,Y}$ takes a raised cosine form:

$\displaystyle
p_{X,Y}(x,y)=\frac{1}{4}\cdot\left[1+\cos(\pi(x-y))\right],\quad(x,y)\in[-1,1]^{2}.$
(97)

Then, it can be verified that the corresponding marginal distributions of
$X,Y$ are uniform distributions $p_{X}=p_{Y}=\mathrm{Unif}([-1,1])$. In
addition, the resulting CDK function is

$\displaystyle\mathfrak{i}_{X;Y}(x,y)=\frac{p_{X,Y}(x,y)-p_{X}(x)p_{Y}(y)}{p_{X}(x)p_{Y}(y)}$
$\displaystyle=\cos(\pi(x-y)).$ (98)

Note that we have $\cos(\pi(x-y))=\cos(\pi x+\theta_{0})\cdot\cos(\pi
y+\theta_{0})+\sin(\pi x+\theta_{0})\cdot\sin(\pi y+\theta_{0})$, for any
$\theta_{0}\in[-\pi,\pi)$. Therefore, we have
$\operatorname{rank}(\mathfrak{i}_{X;Y})=2$, and the associated dependence
modes are given by $\sigma_{1}=\sigma_{2}=1/2$ and the maximal correlation
functions

$\displaystyle f_{1}^{*}(x)=\sqrt{2}\cdot\cos(\pi x+\theta_{0}),\qquad
f_{2}^{*}(x)=\sqrt{2}\cdot\sin(\pi x+\theta_{0}),$ (99a) $\displaystyle
g_{1}^{*}(y)=\sqrt{2}\cdot\cos(\pi y+\theta_{0}),\qquad
g_{2}^{*}(y)=\sqrt{2}\cdot\sin(\pi y+\theta_{0}),$ (99b)

for any $\theta_{0}\in[-\pi,\pi)$.

(a) Histogram of Training Data

(b) Learned Dependence Modes

Figure 14: Dependence modes learned from continuous data, in consistent with
theoretical results.

During this experiment, we first generate $N=50\,000$ sample pairs of $X,Y$
for training, with histogram showing in Figure 14a. Then, we learn $k=2$
dimensional features $f_{1},f_{2}$ of $\mathcal{X}$ and $g_{1},g_{2}$ of
$\mathcal{Y}$ by maximizing the nested H-score (35), where $f$ and $g$ are
parameterized neural feature extractors detailed in Appendix D.1.2. Figure 14b
shows the learned functions after normalization (37). The learned results well
match the theoretical results (99): (i) The learned $f_{1}^{*}$ and
$f_{2}^{*}$ are sinusoids differ in phase by $\pi/2$, and (ii) $g_{i}^{*}$
coincides with $f_{i}^{*}$, for each $i=1,2$. It is also worth mentioning that
due to the degeneracy $\sigma_{1}=\sigma_{2}$, the initial phase $\theta_{0}$
in learned sinusoids (99) can be different during each run of the training
algorithm.

Figure 15: MMSE estimators ${\mathbb{E}}\left[\psi_{i}(Y)|X=x\right]$ obtained
from learning dependence modes, in comparison with theoretical results.

Based on the learned dependence modes, we then demonstrate estimating
functions of $Y$ based on observed $X=x$. Here, we consider the functions
$\psi_{1}(y)=y,\psi_{2}(y)=y^{2},\psi_{3}(y)=e^{y}$. From Proposition 11, we
can compute the learned MMSE estimator
${\mathbb{E}}\left[\psi_{i}(Y)|X=x\right]$ for each $i$, by estimating
${\mathbb{E}}\left[\psi_{i}(Y)\right]$ and $\Lambda_{\psi_{i},g}$ from the
training set and then applying (24).

For comparison, we compute the theoretical values

$\displaystyle{\mathbb{E}}\left[\psi_{i}(Y)|X=x\right]=\int_{-1}^{1}p_{Y|X}(y|x)\psi_{i}(y)\,dy,$

with $p_{Y|X}(y|x)=\frac{1}{2}\cdot\left[1+\cos(\pi(y-x))\right]$, which gives

$\displaystyle{\mathbb{E}}\left[Y\middle|X=x\right]=\frac{1}{\pi}\cdot\sin(\pi
x),\qquad\displaystyle{\mathbb{E}}\left[Y^{2}\middle|X=x\right]=\frac{1}{3}-\frac{2}{\pi^{2}}\cos(\pi
x),$ (100a)
$\displaystyle{\mathbb{E}}\left[e^{Y}\middle|X=x\right]=\frac{e^{2}-1}{2e(1+\pi^{2})}\cdot(\pi\sin(\pi
x)-\cos(\pi x)+\pi^{2}+1).$ (100b)

Figure 15 shows the learned estimators, which are consistent with the
theoretically optimal estimators given by (100).

####  7.1.3 Sequential Data

We proceed with an example of learning dependence modes among sequence pairs.
For simplicity, we consider binary sequences $\underline{X}$ and
$\underline{Y}$, of lengths $l$ and $m$, respectively. Suppose we have the
Markov relation $\underline{X}-U-V-\underline{Y}$ for some unobserved binary
factors $U,V\in\mathcal{U}=\mathcal{V}=\\{0,1\\}$. In addition, we
assume777For convenience, we adopt the vector notation to represent sequences.
$\underline{X}=(X_{1},\dots,X_{l})^{\mathrm{T}},\underline{Y}=(Y_{1},\dots,Y_{m})^{\mathrm{T}}$
satisfy

$\displaystyle\underline{X}|U=i\sim\mathrm{BMS}(l,q_{i}),\quad\underline{Y}|V=i\sim\mathrm{BMS}(m,q_{i}),\quad\text{for
$i=0,1$},$ (101)

where $\mathrm{BMS}(l,q)$ denotes the distribution of a binary first-order
Markov sequence of length $l$ and state flipping probability $q$. The
corresponding state transition diagram is shown in Figure 16a. Therefore, if
$\underline{Z}\sim\mathrm{BMS}(l,q)$, then $Z_{1}\sim\mathrm{Unif}(\\{0,1\\})$
and $(Z_{1},\dots,Z_{l})$ forms a first order Markov chain over binary states
$\\{0,1\\}$, and flipping probability $\mathbb{P}\left\\{Z_{i+1}\neq
Z_{i}|Z_{i}=z\right\\}=q$ for both $z=0,1$. Formally,
$\underline{Z}\sim\mathrm{BMS}(l,q)$ if and only if

$\displaystyle
P_{\underline{Z}}(\underline{z})=\frac{1}{2}\cdot\prod_{i=1}^{l-1}\left[(1-q)^{\delta_{z_{i}z_{i+1}}}\cdot
q^{1-\delta_{z_{i}z_{i+1}}}\right]\quad\text{for all
$(z_{1},\dots,z_{l})\in\\{0,1\\}^{l}$}.$

$0$$1$$1-q$$q$$q$$1-q$

(a) State Transitions of $\mathrm{BMS}(l,q)$

(b) Generated Samples of $\underline{X},\underline{Y}$ Pair

Figure 16: Sequential Data: Model and Generated $\underline{X},\underline{Y}$
Samples Pairs

As a consequence, the resulting alphabets are
$\mathcal{X}=\\{0,1\\}^{l},\mathcal{Y}=\\{0,1\\}^{m}$, with sizes
$|\mathcal{X}|=2^{l},|\mathcal{Y}|=2^{m}$. In our experiment, we set
$l=40,m=30,q_{0}=0.1,q_{1}=0.9$, and use the following joint distribution
$P_{U,V}$:

Prob. | $U=0$ | $U=1$
---|---|---
$V=0$ | 0.1 | 0.2
$V=1$ | 0.4 | 0.3

We generate $N=50\,000$ training sample pairs of
$\underline{X},\underline{Y}$, with instances shown in Figure 16b. We also
generate $N^{\prime}=10\,000$ sample pairs in the same manner, as the testing
dataset.

Figure 17: Features Learned from Sequential Data $\underline{X}$ and
$\underline{Y}$

Then, we learn $k=1$ dimensional features $f$ and $g$ by maximizing
$\mathscr{H}(f,g)$ over the training set. We plot the extracted features in
Figure 17. In particular, each point represents an
$(f(\underline{x}),g(\underline{y}))$ pair evaluated on an instance from
testing set, with corresponding values of binary factors $(U,V)$ shown for
comparison. For ease of demonstration, here we plot only $1,000$ sample pairs
randomly chosen from the testing set. As we can see in the figure, the learned
features are clustered according to the underlying factors. This essentially
reconstructs the hidden factors $U,V$. For example, one can apply a standard
clustering algorithm on the features, e.g., $k$-means (Hastie et al., 2009),
then count the proportion of each cluster, to obtain an estimation of
$P_{U,V}$ up to permutation of symbols.

For a closer inspection, we can compare the learned features with the
theoretical results, formalized as follows. A proof is provided in Appendix
C.21.

###### Proposition 34

Suppose $\underline{X},\underline{Y}$ satisfy the Markov relation
$\underline{X}-U-V-\underline{Y}$ with $U,V\in\\{0,1\\}$ and the conditional
distributions (101). Then, we have
$\operatorname{rank}(\mathfrak{i}_{X;Y})\leq 1$, and the corresponding maximal
correlation functions $f_{1}^{*},g_{1}^{*}$ satisfy

$\displaystyle f_{1}^{*}(\underline{x})$
$\displaystyle=c\cdot\left[\tanh\left(2w\cdot\varphi(\underline{x})+b_{U}\right)-\tanh(b_{U})\right],$
(102a) $\displaystyle g_{1}^{*}(\underline{y})$
$\displaystyle=c^{\prime}\cdot\left[\tanh\left(2w\cdot\varphi(\underline{y})+b_{V}\right)-\tanh(b_{V})\right],$
(102b)

for some $c,c^{\prime}\in\mathbb{R}$, where
$w\triangleq\frac{1}{2}\log\frac{q_{1}}{q_{0}},b_{U}\triangleq\frac{1}{2}\log\frac{P_{U}(1)}{P_{U}(0)},b_{V}\triangleq\frac{1}{2}\log\frac{P_{V}(1)}{P_{V}(0)}$,
and where we have defined $\varphi\colon\\{0,1\\}^{*}\to\mathbb{R}$, such that
for each $\underline{z}=(z_{1},\dots,z_{l})^{\mathrm{T}}\in\\{0,1\\}^{l}$, we
have
$\displaystyle\varphi(\underline{z})\triangleq\frac{l-1}{2}-\sum_{i=1}^{l-1}\delta_{z_{i}z_{i+1}}.$

Then, we compute the correlation coefficients between $f(\underline{X})$ and
$f_{1}^{*}(\underline{X})$, and between $g(\underline{Y})$,
$g_{1}^{*}(\underline{Y})$, respectively, using sample pairs in the testing
set. The absolute values of both correlation coefficients are greater than
$0.99$, demonstrating the effectiveness of the learning algorithm.

###  7.2 Learning With Orthogonality Constraints

We verify the feature learning with orthogonality constraints presented in
Section 4.4 on the same dataset used for Section 7.1.2. Here, we consider the
settings ${\bar{k}}=k=1$, i.e., we learn one-dimensional feature $f(x)$
uncorrelated to given one dimensional feature
$\phi\in\mathcal{F}_{\mathcal{X}}$.

Note that without the orthogonality constraint [cf. (99)], the optimal feature
will be sinusoids with any initial phase, i.e.,
$f_{1}^{*}(x)=\sqrt{2}\cdot\cos(\pi x+\theta_{0})$ for any
$\theta_{0}\in[-\pi,\pi)$. Here, we consider the following two choices of
$\phi$, $(x\mapsto x)$ and $(x\mapsto x^{2})$, which are even and odd
functions, respectively. Since the underlying $p_{X}$ is uniform on $[-1,1]$,
we can verify the optimal features under the two constraints are
$f^{*}_{1}(x)=\sqrt{2}\cos(\pi x)$ for $\phi(x)=x$, and
$f^{*}_{1}(x)=\sqrt{2}\sin(\pi x)$ for $\phi(x)=x^{2}$, respectively.

(a) $\phi(x)=x$ (b) $\phi(x)=x^{2}$

Figure 18: Learning feature uncorrelated to given $\phi$ under different
settings of $\phi$. The learned results are compared with theoretical results.

By maximizing the nested H-score restricted to $\bar{f}=\phi$ [cf. (49)], we
can learn the optimal feature $f_{1}^{*}$, as shown in Figure 18. The learned
features are in consistent with the theoretical ones, validating the
effectiveness of the learning algorithm.

###  7.3 Learning With Side Information

We design an experiment to verify the connection between our learning
algorithm and the multitask classification DNN, as demonstrated in Theorem 25.

(a) Learned feature $f(x)$ (b) Learned feature $g(s,y)$

Figure 19: Experimental verification of the connection between learning with
side information and training a multihead neural network.

In particular, we consider the discrete $X,Y,S$ with
$|\mathcal{X}|=8,|\mathcal{S}|=|\mathcal{Y}|=3$, and randomly choose a joint
distribution $P_{X,S,Y}$ on $\mathcal{X}\times\mathcal{S}\times\mathcal{Y}$.
Then, we generate $N=50\,000$ training samples of $(X,S,Y)$ triples. In our
implementation, we set ${\bar{k}}=|\mathcal{S}|-1=2,k=1$ and maximize the
nested H-score configured by $\mathcal{C}_{\mathsf{MC}}$ [cf. (56)] on the
training set.

For comparison, we also train the multihead network shown in Figure 10 by
maximizing the log-likelihood function (67) to learn the corresponding feature
$f$ and weight matrices $G_{s}$ for all $s\in\mathcal{S}$. Then, we convert
the weights to $g\in\mathcal{F}_{\mathcal{S}\times\mathcal{Y}}$ via the
correspondence [cf. (66)] $g(s,y)=G_{s}(1,y)$. The features learned by our
algorithm (labeled as “SideInfo”) and the multihead neural network are shown
in Figure 19, where the results are consistent.

###  7.4 Multimodal Learning With Missing Modalities

To verify the multimodal learning algorithms presented in Section 6, we
consider multimodal classification problems in two different settings. Suppose
$X_{1},X_{2}$ are multimodal data variables, and $Y\in\mathcal{Y}=\\{-1,1\\}$
denotes the binary label to predict. In the first setting, we consider the
training set with complete $(X_{1},X_{2},Y)$ samples. In the second setting,
only the pairwise observations of $(X_{1},X_{2})$, $(X_{1},Y)$, and
$(X_{2},Y)$ are available, presented in three different datasets.

In both settings, we set $\mathcal{X}_{1}=\mathcal{X}_{2}=[-1,1]$ with

$\displaystyle
P_{X_{1},X_{2}}(x_{1},x_{2})=\frac{1}{4}\cdot\left[1+\cos(2\pi(x_{1}-x_{2}))\right].$
(103)

We consider predicting $Y$ based on the learned features, where some modality
in $X_{1}$, $X_{2}$ might be missing during the prediction.

####  7.4.1 Learning from Complete Observations

We consider the $X_{1},X_{2},Y$ dependence specified by (103) and the
conditional distribution

$\displaystyle
P_{Y|X_{1},X_{2}}(y|x_{1},x_{2})=\frac{1}{2}+\frac{y}{4}\cdot\left[\cos(\pi
x_{1})+\cos(\pi x_{2})+\cos(\pi(x_{1}+x_{2}))\right]$ (104)

for $x_{1},x_{2}\in[-1,1]$ and $y=\pm 1$. It can be verified that $P_{Y}$
satisfies $P_{Y}(1)=P_{Y}(-1)=\frac{1}{2}$. The corresponding CDK function and
dependence components [cf. (71)] are given by

$\displaystyle\mathfrak{i}_{X_{1},X_{2};Y}(x_{1},x_{2},y)$
$\displaystyle=\frac{y}{2}\cdot\left[\cos(\pi x_{1})+\cos(\pi
x_{2})+\cos(\pi(x_{1}+x_{2}))\right]$ (105a)
$\displaystyle[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)$
$\displaystyle=\frac{y}{2}\cdot\left[\cos(\pi x_{1})+\cos(\pi x_{2})\right],$
(105b)
$\displaystyle[\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)$
$\displaystyle=\frac{y}{2}\cdot\cos(\pi(x_{1}+x_{2})).$ (105c)

Therefore, we have
$\operatorname{rank}(\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}))=\operatorname{rank}(\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y}))=1$,
and the functions obtained from modal decompositions [cf. (88)] are
$\bar{g}_{1}^{*}(y)=g_{1}^{*}(y)=y$ and

$\displaystyle\bar{f}_{1}^{*}(x_{1},x_{2})=\cos(\pi x_{1})+\cos(\pi
x_{2}),\qquad f_{1}^{*}(x_{1},x_{2})=\sqrt{2}\cdot\cos(\pi(x_{1}+x_{2})).$
(106)

(a) $X_{1},X_{2}$ Histogram

(b) $\bar{f}_{1}^{*}(x_{1},x_{2})$ (c) $f_{1}^{*}(x_{1},x_{2})$

(d) $P_{Y|X_{1},X_{2}}(1|x_{1},x_{2})$

Figure 20: Features and posterior probability learned from multimodal data
$X_{1},X_{2},Y$, in comparison with theoretical results.

In the experiment, we first generate $N=50\,000$ triples of $(X_{1},X_{2},Y)$
for training. The histogram of $(X_{1},X_{2})$ pair is shown in Figure 20a.
Then, we set ${\bar{k}}=k=1$, and learn the features $\bar{f},f,\bar{g},g$ by
maximizing the nested H-score configured by $\mathcal{C}_{\sf BI}$ [cf. (75)].
We then normalize $\bar{f}$, $f$ to obtain the estimated $\bar{f}_{1}^{*}$ and
$f_{1}^{*}$, and compute the posterior probability
$P_{Y|X_{1},X_{2}}(1|x_{1},x_{2})$ based on (80a). The results of learned
features $\bar{f}_{1}^{*}$, $f_{1}^{*}$ and posterior
$P_{Y|X_{1},X_{2}}(1|x_{1},x_{2})$ are shown in Figure 20, which are
consistent with the theoretical results.

We then consider the prediction problem with missing modality, i.e., predict
label $Y$ based on unimodal data $X_{1}$ or $X_{2}$. In particular, based on
learned $\bar{f}=\bar{f}^{(1)}+\bar{f}^{(2)}$, we train two separate networks
to operate as $\tau_{{1}}$ and $\tau_{{2}}$, then apply (80b) and (80c) to
estimate the posteriors $P_{Y|X_{1}}$ and $P_{Y|X_{2}}$. Then, for each
$i=1,2$, the MAP prediction of $Y$ based on observed $X_{i}=x_{i}$ can be
obtained by comparing $P_{Y|X_{i}}(1|x_{i})$ with the threshold $1/2$, via

$\displaystyle\operatorname*{arg\,max}_{y\in\mathcal{Y}}P_{Y|X_{i}}(y|x_{i})=\begin{cases}1&\text{if
$P_{Y|X_{i}}(1|x_{i})>1/2$},\\\ -1&\text{if $P_{Y|X_{i}}(1|x_{i})\leq
1/2$}.\end{cases}$

(a) $P_{Y|X_{1}}(1|x_{1})$ (b) $P_{Y|X_{2}}(1|x_{2})$

Figure 21: Label prediction from unimodal data using learned features

We plot the estimated results in Figure 21, in comparison with the threshold
$1/2$ and the theoretical values

$\displaystyle P_{Y|X_{i}}(y|x)=\frac{1}{2}+\frac{y}{4}\cdot\cos(\pi
x),\quad\text{for $i=1,2$.}$ (107)

From the figure, the estimated posteriors $P_{Y|X_{1}},P_{Y|X_{2}}$ have
consistent trends with the ground truth posteriors, and the induced $Y$
predictions are well aligned.

####  7.4.2 Learning from Pairwise Observations

We proceed to consider the multimodal learning with only pairwise
observations. Specifically, we consider the joint distribution of
$(X_{1},X_{2},Y)$, specified by (103) and

$\displaystyle
P_{Y|X_{1},X_{2}}(y|x_{1},x_{2})=\frac{1}{2}+\frac{y}{4}\cdot\left[\cos(\pi
x_{1})+\cos(\pi x_{2})\right].$ (108)

for $x_{1},x_{2}\in[-1,1]$ and $y=\pm 1$. It can be verified that
$P_{Y}(1)=P_{Y}(-1)=\frac{1}{2}$, and the associated CDK function satisfies

$\displaystyle\mathfrak{i}_{X_{1},X_{2};Y}(x_{1},x_{2},y)=\frac{y}{2}\cdot\left[\cos(\pi
x_{1})+\cos(\pi x_{2})\right].$ (109)

and
$\mathfrak{i}_{X_{1},X_{2};Y}=\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})$.
Therefore, the interaction dependence component
$\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})=0$, and the joint dependence
can be learned from all pairwise samples, as we discussed in Section 6.5.1.

We then construct an experiment to verify learning joint dependence from all
pairwise observations. Specifically, we generate $N=50\,000$ triples of
$(X_{1},X_{2},Y)$ from (103) and (108). Then, we adopt the decomposition (86)
on each triple, to obtain three pairwise datasets with samples of
$(X_{1},X_{2})$, $(X_{1},Y)$, $(X_{2},Y)$, where each dataset has $N$ sample
pairs.

(a) $\bar{f}_{1}^{*}(x_{1},x_{2})$ (b) $P_{Y|X_{1},X_{2}}(1|x_{1},x_{2})$

Figure 22: Features and posterior probability learned from pairwise datasets
of $(X_{1},X_{2}),(X_{1},Y)$, and $(X_{2},Y)$, in comparison with theoretical
results.

We use these three pairwise datasets for training to learn one dimensional
$\bar{f}\in\mathcal{F}_{\mathcal{X}},\bar{g}\in\mathcal{F}_{\mathcal{Y}}$ that
maximize $\mathscr{H}(\bar{f},\bar{g})$. Here, we compute
$\mathscr{H}(\bar{f},\bar{g})$ based on the minibatches from the three
pairwise datasets, according to (90). Based on learned $\bar{f},\bar{g}$, we
then compute the normalized $\bar{f}_{1}^{*}$ and posterior distribution
$P_{Y|X_{1},X_{2}}(1|x_{1},x_{2})$, as shown in Figure 22, where the learned
results match the theoretical values.

(a) $P_{Y|X_{1}}(1|x_{1})$ (b) $P_{Y|X_{2}}(1|x_{2})$

Figure 23: Label prediction from unimodal data, with features learned from
pairwise datasets.

Similar to the previous setting, we consider the unimodal prediction problem,
and show the estimated results in Figure 23. It is worth noting that from
(105b), (109), the joint distributions $P_{X_{1},X_{2},Y}$ in both settings
contain the same bivariate dependence component. Therefore, the theoretical
results for $\bar{f}_{1}^{*}$ and $P_{Y|X_{1}},P_{Y|X_{2}}$ are the same.

##  8 Related Works

##### HGR Maximal Correlation and Learning Algorithms

The Hirschfeld–Gebelein–Rényi (HGR) maximal correlation (Hirschfeld, 1935;
Gebelein, 1941; Rényi, 1959) provides an important connection between
statistical dependence and functional space. The same concept has been studied
with various formulations, and often in different terminologies, such as
correspondence analysis (Greenacre, 2017), functional canonical variates
(Buja, 1990), and principal inertial components (du Pin Calmon et al., 2017).
See, e.g., (Huang et al., 2019, Section II), for a detailed discussion. The
first practical algorithm for learning such functions is the alternating
conditional expectations (ACE) algorithm (Breiman and Friedman, 1985; Buja,
1990), which learns the top dependence mode and are mostly used for processing
low-dimensional data. To learn maximal correlation functions by neural feature
extractors, several different training objectives have been proposed (Wang et
al., 2019; Hsu et al., 2021). Specifically, the H-score was derived from the
low-rank approximation, referred as the Soft-HGR objective (Wang et al.,
2019). It was also introduced as the local approximation of the log-loss
function (Xu et al., 2022). The application of H-score for classification
tasks was exploited in Xu and Huang (2020). In an independent line of work
(HaoChen et al., 2021), a special form of H-score was proposed for self-
supervised learning tasks, referred to as the _spectral contrastive loss_
therein.

##### Informative Features and Local Analysis

Huang et al. (2019) provided an in-depth characterization of informative
features by applying information-theoretic tools, with a focus on and
bivariate learning problems and the local analysis regime. In particular, it
was shown in Huang et al. (2019) that for bivariate settings, there are a
series of statistical and learning formulations that lead to the same optimal
features, characterized as the HGR maximal correlation functions. For
multivariate problems, the features introduced in (88a) were also studied by
Xu and Huang (2021) in characterizing distributed hypothesis testing problems.

##### Decomposition of Probability Distributions

In addition to the modal decomposition, the decomposition of probability
distributions has also been studied, particularly in the context of
information geometry (Amari and Nagaoka, 2000). For example, Amari (2001)
established a decomposition in distribution space, which also investigated the
maximum entropy formulation (84), cf. (Amari, 2001, Theorem 7).

##  9 Conclusions and Discussions

We have presented a framework for designing learning systems with neural
feature extractors, which allows us to learn informative features and assemble
them to build different inference models. Based on the feature geometry, we
convert learning problems to corresponding geometric operations in feature
spaces. We then introduce the nesting technique for implementing such
geometric operations in feature spaces, which provides a systematic design of
both feature learning and feature assembling modules. We demonstrate the
applications by considering conditioned inference and multimodal learning
problems.

The established framework provides abundant opportunities for further
explorations. For the ease of presentation, our discussions focused on basic
and typically simplified settings, which can be extended to more general
analyses and applications. For example, we can extend the analysis to
decompose the dependence of a random processes and develop the corresponding
feature learning algorithms. Specifically, we can get refined dependence
components by iteratively applying the decomposition (40). Our
characterizations on learning factors, such as expressive power of feature
extractors (cf. Section 3.3) and sample size of training data (cf. Section
6.5.2), can also be extended to more in-depth discussions. In addition, we can
apply the feature geometry to analyze existing algorithm designs, e.g.,
evaluating kernel choices in kernel methods (Xu and Zheng, 2023).

Acknowledgments and Disclosure of Funding

This work was supported in part by the National Science Foundation (NSF) under
Award CNS-2002908 and the Office of Naval Research (ONR) under grant
N00014-19-1-2621.

##  Appendix Appendix A Feature Space: Canonical Basis and Vector Notations

For completeness, we briefly summarize the information vector and
corresponding matrix conventions introduced in Huang et al. (2019) and related
works.

To begin, we assume the random variable $Z$ takes finite many possible values,
i.e., $|\mathcal{Z}|<\infty$, then the resulting feature space
$\mathcal{F}_{\mathcal{Z}}$ is a finite-dimensional vector space. It is
sometimes more convenient to represent features using vector and matrix
notations. Specifically, each $f\in\mathcal{F}_{\mathcal{Z}}$ can be
equivalently expressed as one vector in $\mathbb{R}^{|\mathcal{Z}|}$ as
follows. Suppose $R_{Z}$ is the metric distribution, then we can construct an
orthonormal basis $\left\\{\mathsf{b}_{Z}^{(z)}\colon
z\in\mathcal{Z}\right\\}$ of $\mathcal{F}_{\mathcal{Z}}$, where

$\displaystyle\mathsf{b}_{Z}^{(z)}(z^{\prime})\triangleq\frac{\delta_{zz^{\prime}}}{\sqrt{R_{Z}(z)}}\quad\text{for
all $z,z^{\prime}\in\mathcal{Z}$}.$ (110)

We refer to this basis as the canonical basis of $\mathcal{F}_{\mathcal{Z}}$.
For all $f\in\mathcal{F}_{\mathcal{Z}}$, we can represent $f$ as a linear
combination of these basis functions, i.e.,

$\displaystyle
f=\sum_{z^{\prime}\in\mathcal{Z}}\xi(z^{\prime})\cdot\mathsf{b}_{Z}^{(z^{\prime})},$

where the coefficient $\xi(z)$ for each $z\in\mathcal{Z}$ is given by

$\displaystyle\xi(z)=\left\langle
f,\mathsf{b}_{Z}^{(z)}\right\rangle=f(z)\sqrt{R_{Z}(z)}.$ (111)

In particular, when $f$ is the density ratio $\tilde{\ell}_{P_{Z}}$ for some
$P_{Z}\in\mathcal{F}_{\mathcal{Z}}$, the corresponding coefficient $\xi(z)$
for each $z\in\mathcal{Z}$ is

$\displaystyle\xi(z)=\left\langle\tilde{\ell}_{P_{Z}},\mathsf{b}_{Z}^{(z)}\right\rangle=\frac{P_{Z}(z)-R_{Z}(z)}{\sqrt{R_{Z}(z)}}.$
(112)

This establishes a one-to-one correspondence between $\tilde{\ell}_{P_{Z}}$
(or $P_{Z}$) and the vector
$\underline{\xi}\triangleq[\xi(z),z\in\mathcal{Z}]^{\mathrm{T}}\in\mathbb{R}^{|\mathcal{Z}|}$,
which is referred to as the information vector associated with
$\tilde{\ell}_{P_{Z}}$ (or $P_{Z}$).

Similarly, for $X$ and $Y$ with $|\mathcal{X}|<\infty,|\mathcal{Y}|<\infty$,
we can represent the CDK function
$\mathfrak{i}_{X;Y}\in\mathcal{F}_{\mathcal{X}\times\mathcal{Y}}$ as an
$|\mathcal{X}|\times|\mathcal{Y}|$ matrix $\tilde{B}_{X;Y}$:

$\displaystyle\tilde{B}_{X;Y}(x;y)$
$\displaystyle\triangleq\frac{P_{X,Y}(x,y)-P_{X}(x)P_{Y}(y)}{\sqrt{P_{X}(x)}\sqrt{P_{Y}(y)}},$
(113)

which is referred to as the canonical dependence matrix (CDM) of $X$ and $Y$.
With the metric distribution $P_{X}P_{Y}$ on
$\mathcal{F}_{\mathcal{X}\times\mathcal{Y}}$, each $\tilde{B}_{X;Y}(x;y)$ is
the coefficient associated with the basis function $\mathsf{b}_{X,Y}^{(x,y)}$
[cf. (110)].

In addition, we have
$\operatorname{rank}(\tilde{B}_{X;Y})=\operatorname{rank}(\mathfrak{i}_{X;Y})$.
Suppose $\sigma_{i}(f_{i}^{*}\otimes
g_{i}^{*})=\operatorname{\zeta}_{i}(\mathfrak{i}_{X;Y})$, for each $1\leq
i\leq\operatorname{rank}(\mathfrak{i}_{X;Y})$, then $\sigma_{i}$ is the $i$-th
singular vector of $\tilde{B}_{X;Y}$, and the corresponding $i$-th left and
right singular vector pair are given by
$\underline{\xi}_{i}^{X}\in\mathbb{R}^{|\mathcal{X}|},\underline{\xi}_{i}^{Y}\in\mathbb{R}^{|\mathcal{Y}|}$,
with [cf. (111)]

$\displaystyle\xi_{i}^{X}(x)=\sqrt{P_{X}(x)}\cdot
f_{i}^{*}(x),\quad\xi_{i}^{Y}(y)=\sqrt{P_{Y}(y)}\cdot g_{i}^{*}(y).$ (114)

Therefore, for small-scale discrete data, we can use the connection (114) to
obtain the modal decomposition by solving the SVD of the corresponding CDM.

##  Appendix Appendix B Characterization of Common Optimal Solutions

We consider optimization problems defined on a common domain $\mathcal{D}$.
Specifically, given $k$ functions $h_{i}\colon\mathcal{D}\to\mathbb{R}$,
$i=1,\dots,k$, let us consider optimization problems

$\displaystyle\operatorname*{maximize}_{u\in\mathcal{D}}h_{i}(u),\quad
i=1,\dots,k.$ (115)

For each $i=1,\dots,k$, we denote the optimal solution and optimal value for
$i$-th problem by

$\displaystyle\mathcal{D}_{i}^{*}\triangleq\operatorname*{arg\,max}_{u\in\mathcal{D}}h_{i}(u),\quad
t_{i}^{*}\triangleq\max_{u\in\mathcal{D}}h_{i}(u),$ (116)

respectively, and suppose $\mathcal{D}_{i}^{*}\neq\varnothing,i=1,\dots,k$.

Then, the set $\mathcal{D}^{*}\triangleq\cap_{i=1}^{k}\mathcal{D}_{i}^{*}$
represents the collection of common optimal solutions for all $k$ optimization
problems (115). When such common solutions exist, i.e., $\mathcal{D}^{*}$ is
nonempty, we can obtain $\mathcal{D}^{*}$ by a single optimization program, by
using an objective that aggregates original $k$ objectives. We formalize the
result as follows.

###### Proposition 35

If $\mathcal{D}^{*}\neq\varnothing$, we have
$\displaystyle\mathcal{D}^{*}=\operatorname*{arg\,max}_{u\in\mathcal{D}}\,\Gamma(h_{1}(u),h_{2}(u),\dots,h_{k}(u))$
for every $\Gamma\colon\mathbb{R}^{k}\to\mathbb{R}$ that is strictly
increasing in each argument.

Proof  Let
$\mathcal{D}^{**}\triangleq\operatorname*{arg\,max}_{u\in\mathcal{D}}\,h(u)$
with $h(u)\triangleq\Gamma(h_{1}(u),h_{2}(u),\dots,h_{k}(u))$. Then the
proposition is equivalent to $\mathcal{D}^{*}=\mathcal{D}^{**}$. We then
establish $\mathcal{D}^{*}\subset\mathcal{D}^{**}$ and
$\mathcal{D}^{**}\subset\mathcal{D}^{*}$, respectively.

To prove $\mathcal{D}^{*}\subset\mathcal{D}^{**}$, take any
$u^{*}\in\mathcal{D}^{*}$. Then, for all $u\in\mathcal{D}$, we have
$h_{i}(u)\leq h_{i}(u^{*})$, $i=1,\dots,k$, which implies that

$\displaystyle
h(u)=\Gamma(h_{1}(u),\dots,h_{k}(u))\leq\Gamma(h_{1}(u^{*}),\dots,h_{k}(u^{*}))=h(u^{*}).$

Therefore, we have $u^{*}\in\mathcal{D}^{**}$. Since $u^{*}$ can be
arbitrarily chosen from $\mathcal{D}^{*}$, we have
$\mathcal{D}^{*}\subset\mathcal{D}^{**}$.

We then establish $\mathcal{D}^{**}\subset\mathcal{D}^{*}$, which is
equvialent to

$\displaystyle(\mathcal{D}\setminus\mathcal{D}^{*})\subset(\mathcal{D}\setminus\mathcal{D}^{**}).$
(117)

Note that (117) is trivially true if
$\mathcal{D}\setminus\mathcal{D}^{*}=\varnothing$. Otherwise, take any
$u^{\prime}\in\mathcal{D}\setminus\mathcal{D}^{*}$. Then, we have
$h_{i}(u^{\prime})\leq h_{i}(u^{*})$ for all $i\in[k]$, and the strict
inequality holds for at least one $i\in[k]$. This implies that

$\displaystyle
h(u^{\prime})=\Gamma(h_{1}(u^{\prime}),\dots,<h_{k}(u^{\prime}))<\Gamma(h_{1}(u^{*}),\dots,h_{k}(u^{*}))=h(u^{*}).$

Hence, $u^{\prime}\notin\mathcal{D}^{**}$, and thus
$u^{\prime}\in(\mathcal{D}\setminus\mathcal{D}^{**})$, which establishes
(117).

To apply Proposition 35, the first step is to test the existence of common
optimal solutions. A naive test is to solve all optimization problems (115)
and then check the definition, which can be difficult in practice. Instead, we
can consider a related multilevel optimization problem as follows. Let
$\mathcal{D}_{0}\triangleq\mathcal{D}$, and for each $i=1,\dots,k$ as, we
define $\mathcal{D}_{i}$ and $t_{i}\in\mathbb{R}$ as

$\displaystyle\mathcal{D}_{i}\triangleq\operatorname*{arg\,max}_{u\in\mathcal{D}_{\\!i-1}}\,h_{i}(u),\quad
t_{i}\triangleq\max_{u\in\mathcal{D}_{\\!i-1}}\,h_{i}(u).$ (118)

Note that for each $i$, the $\mathcal{D}_{i}$ solved by level $i$ optimization
problem gives the elements in $\mathcal{D}_{i}$ that maximize $h_{i}$, and
$t_{i}$ denotes the corresponding optimal value. Therefore, $\mathcal{D}_{k}$
can be obtained by sequentially solving $k$ optimization problems as defined
in (118).

Then, the following result provides an approach to test the existence of
common optimal solutions.

###### Proposition 36

The following three statements are equivalent:

  1. 1.

The optimization problems (115) have common optimal solutions, i.e.,
$\mathcal{D}^{*}\neq\varnothing$;

  2. 2.

$t_{i}=t_{i}^{*}$ for all $i=1,\dots,k$;

  3. 3.

$\mathcal{D}_{i-1}\cap\mathcal{D}_{i}^{*}\neq\varnothing$ for all
$i=1,\dots,k$.

In addition, if one of these statements holds, then we have
$\mathcal{D}_{k}=\mathcal{D}^{*}$.

Proof  We establish the equivalence of the statements 1 to 3, by showing that
“1” $\implies$ “2”, “2” $\implies$ “3”, and “3” $\implies$ “1”.

#####  “1” $\implies$ “2”

Suppose “1” holds, and take any
$u^{*}\in\mathcal{D}^{*}=\cap_{i=1}^{k}\mathcal{D}_{i}^{*}$. We then establish
“2” by induction. First, note that $u^{*}\in\mathcal{D}_{0}$. For the
induction step, we can show that for each $i=1,\dots,k$, if
$u^{*}\in\mathcal{D}_{i-1}$, then $u^{*}\in\mathcal{D}_{i}$ and
$t_{i}^{*}=t_{i}$. Indeed, we have

$\displaystyle t_{i}^{*}\geq t_{i}=\max_{u\in\mathcal{D}_{i-1}}h_{i}(u)\geq
h_{i}(u^{*})=t_{i}^{*},$

where the first inequality follows from the fact that
$\mathcal{D}=\mathcal{D}_{0}\supset\mathcal{D}_{1}\supset\dots\supset\mathcal{D}_{k}$,
where the second inequality follows from the inductive assumption
$u^{*}\in\mathcal{D}_{i-1}$, and where the last equality follows from that
$u^{*}\in\mathcal{D}^{*}\subset\mathcal{D}_{i}^{*}$.

##### “2” $\implies$ “3”

For each $i=2,\dots,k$, $t_{i}=t_{i}^{*}$ implies that there exists some
$u_{i}\in\mathcal{D}_{i-1}$, such that
$h_{i}(u_{i})=t_{i}=t_{i}^{*}=\max_{u\in\mathcal{D}}h_{i}(u)$, and thus
$u_{i}\in\mathcal{D}_{i}^{*}$. Therefore,
$u_{i}\in\mathcal{D}_{i-1}\cap\mathcal{D}_{i}^{*}$, which establishes “3”.

##### “3” $\implies$ “1”

For each $i=2,\dots,k$, from
$\mathcal{D}_{i-1}\cap\mathcal{D}_{i}^{*}\neq\varnothing$ and the definitions
(116) and (118), we have
$\mathcal{D}_{i}=\mathcal{D}_{i-1}\cap\mathcal{D}_{i}^{*}$. It can also be
verified that $\mathcal{D}_{i}=\mathcal{D}_{i-1}\cap\mathcal{D}_{i}^{*}$ holds
for $i=1$. Therefore, we obtain

$\displaystyle\mathcal{D}_{k}=\mathcal{D}_{k-1}\cap\mathcal{D}_{k}^{*}=\mathcal{D}_{k-2}\cap\mathcal{D}_{k-1}^{*}\cap\mathcal{D}_{k}^{*}=\dots=\mathcal{D}_{0}\cap\left(\bigcap_{i=1}^{k}\mathcal{D}_{i}^{*}\right)=\mathcal{D}\cap\mathcal{D}^{*}=\mathcal{D}^{*}.$
(119)

This implies that
$\mathcal{D}^{*}=\mathcal{D}_{k-1}\cap\mathcal{D}_{k}^{*}\neq\varnothing$.

Finally, from (119) we know that statement 3 implies
$\mathcal{D}_{k}=\mathcal{D}^{*}$. Since all three statements are equivalent,
each statement implies $\mathcal{D}_{k}=\mathcal{D}^{*}$.

##  Appendix Appendix C Proofs

###  C.1 Proof of Proposition 5

Let
$\hat{\gamma}\triangleq\Pi\left(\gamma;{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}}\right)$.
We first consider the second equality of (15), i.e.,

$\displaystyle\operatorname{\zeta}_{\leq
k}(\hat{\gamma})=\operatorname*{arg\,min}_{\begin{subarray}{c}\gamma^{\prime}\colon\gamma^{\prime}=f\otimes
g,\\\
f\in\mathcal{G}_{\mathcal{X}}^{\,k},\,g\in\mathcal{G}_{\mathcal{Y}}^{\,k}\end{subarray}}\|\gamma-\gamma^{\prime}\|$
(120)

For each $k\leq\operatorname{rank}(\hat{\gamma})$, consider
$\gamma^{\prime}=f\otimes g$ with
$f\in\mathcal{G}_{\mathcal{X}}^{\,k},g\in\mathcal{G}_{\mathcal{Y}}^{\,k}$.
Then, we have

$\displaystyle\|\gamma-\gamma^{\prime}\|^{2}=\|\gamma-\hat{\gamma}+\hat{\gamma}-\gamma^{\prime}\|^{2}=\|\gamma-\hat{\gamma}\|^{2}+\|\hat{\gamma}-\gamma^{\prime}\|^{2}\geq\|\gamma-\hat{\gamma}\|^{2}+\|r_{k}(\hat{\gamma})\|^{2},$
(121)

where to obtain the second equality we have used the orthogonality principle
with fact that
$(\gamma-\hat{\gamma})\perp\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}$
and
$(\hat{\gamma}-\gamma^{\prime})\in\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}$.
Note that the inequality in (121) holds with equality if and only if
$\gamma^{\prime}=\operatorname{\zeta}_{\leq k}(\hat{\gamma})$, which
establishes (120).

We then establish
$\operatorname{\zeta}_{k}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})=\operatorname{\zeta}_{k}(\hat{\gamma})$
by induction. To begin, set $k=1$ in (120), and the right hand side becomes
$\operatorname{\zeta}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})$,
which implies
$\operatorname{\zeta}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})=\operatorname{\zeta}(\hat{\gamma})$,
i.e.,
$\operatorname{\zeta}_{1}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})=\operatorname{\zeta}_{1}(\hat{\gamma})$.
As the inductive hypothesis, suppose we have

$\displaystyle\operatorname{\zeta}_{i}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})=\operatorname{\zeta}_{i}(\hat{\gamma}),\quad
i=1,\dots,m.$ (122)

From (14), we have

$\displaystyle\operatorname{\zeta}_{m+1}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})$
$\displaystyle=\operatorname{\zeta}\left(\gamma-\sum_{i=1}^{m}\operatorname{\zeta}_{i}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})\middle|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}}\right)$
$\displaystyle=\operatorname{\zeta}\left(\gamma-\sum_{i=1}^{m}\operatorname{\zeta}_{i}(\hat{\gamma})\middle|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}}\right)$
(123)
$\displaystyle=\operatorname{\zeta}\left(\Pi\left(\gamma-\sum_{i=1}^{m}\operatorname{\zeta}_{i}(\hat{\gamma});{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}}\right)\right)$
(124)
$\displaystyle=\operatorname{\zeta}\left(\hat{\gamma}-\sum_{i=1}^{m}\Pi\left(\operatorname{\zeta}_{i}(\hat{\gamma});{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}}\right)\right)$
(125)
$\displaystyle=\operatorname{\zeta}\left(\hat{\gamma}-\sum_{i=1}^{m}\operatorname{\zeta}_{i}(\hat{\gamma})\right)=\operatorname{\zeta}_{m+1}(\hat{\gamma}),$
(126)

where (123)–(124) follow from the inductive assumption (122), where (125)
follows from the linearity of projection operator. To obtain the first
equality of (125), we have again applied the assumption (122): for
$i=1,\dots,m$,
$\operatorname{\zeta}_{i}(\hat{\gamma})=\operatorname{\zeta}_{i}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})\in\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}$.

Finally, $\operatorname{\zeta}_{\leq
k}(\gamma|\mathcal{G}_{\mathcal{X}},\mathcal{G}_{\mathcal{Y}})=\operatorname{\zeta}_{\leq
k}(\hat{\gamma})$ can be readily obtained by definition. $\blacksquare$

###  C.2 Proof of Proposition 7

It is easy to verify that
$\operatorname{cov}(\hat{f}_{i}^{*},\hat{g}_{i}^{*})=\bigl{\langle}\mathfrak{i}_{X;Y},\hat{f}_{i}^{*}\otimes\hat{g}_{i}^{*}\bigr{\rangle}=\|\operatorname{\zeta}_{i}(\mathfrak{i}^{\prime}_{X;Y})\|=\hat{\sigma}_{i}^{*}$,
where
$\mathfrak{i}^{\prime}_{X;Y}=\Pi\left(\mathfrak{i}_{X;Y};{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}}\right)$.

From Fact 4, we have
$(\hat{f}_{i}^{*},\hat{g}_{i}^{*})=\operatorname*{arg\,max}_{f_{i},g_{i}}\,\bigl{\langle}\mathfrak{i}^{\prime}_{X;Y},f_{i}\otimes
g_{i}\bigr{\rangle}$ where the maximization is taken over all
$f_{i}\in\mathcal{F}_{\mathcal{X}}$ and $g_{i}\in\mathcal{F}_{\mathcal{Y}}$
that satisfy (19).

Therefore, for each $i$ and
$f_{i}\in\mathcal{G}_{\mathcal{X}},g_{i}\in\mathcal{G}_{\mathcal{Y}}$ that
satisfy (19), we have

$\displaystyle\operatorname{cov}(f_{i},g_{i})=\langle\mathfrak{i}_{X;Y},f_{i}\otimes
g_{i}\rangle=\bigl{\langle}\mathfrak{i}^{\prime}_{X;Y},f_{i}\otimes
g_{i}\bigr{\rangle}\leq\hat{\sigma}_{i}^{*}=\operatorname{cov}(\hat{f}_{i}^{*},\hat{g}_{i}^{*}),$

where the second equality follows from that $f_{i}\otimes
g_{i}\in\mathcal{G}_{\mathcal{X}}\otimes\mathcal{G}_{\mathcal{Y}}$. Hence, we
obtain
$(\hat{f}_{i}^{*},\hat{g}_{i}^{*})=\operatorname*{arg\,max}_{f_{i},g_{i}}\operatorname{cov}(f_{i},g_{i})$
where the maximization is taken over
$f_{i}\in\mathcal{G}_{\mathcal{X}},g_{i}\in\mathcal{G}_{\mathcal{Y}}$ with the
constraint (19). $\blacksquare$

###  C.3 Proof of Proposition 11

The result of $\|\mathfrak{i}_{X;Y}\|$ directly follows from Property 1. In
addition,

$\displaystyle
f^{\mathrm{T}}(x)g(Y)=\mathfrak{i}_{X;Y}(x,y)=\frac{P_{X,Y}(x,y)-P_{X}(x)P_{Y}(y)}{P_{X}(x)P_{Y}(y)},\quad\text{for
all }(x,y)\in\mathcal{X}\times\mathcal{Y},$

which implies $P_{Y|X}(y|x)=P_{Y}(y)\left(1+f^{\mathrm{T}}(x)g(y)\right)$,
i.e., (23).

Therefore, for all $\psi\in\mathcal{F}_{\mathcal{Y}}^{\,d}$, we have

$\displaystyle{\mathbb{E}}\left[\psi(Y)|X=x\right]=\sum_{y\in\mathcal{Y}}P_{Y|X}(y|x)\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)(1+f^{\mathrm{T}}(x)g(y))\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)\psi(y)+\sum_{y\in\mathcal{Y}}P_{Y}(y)(\psi(y)g^{\mathrm{T}}(y))f^{\mathrm{T}}(x)$
$\displaystyle={\mathbb{E}}\left[\psi(Y)\right]+\Lambda_{\psi,g}\cdot f(x),$

which gives (24). $\blacksquare$

###  C.4 Proof of Proposition 12

We have

$\displaystyle{\mathbb{E}}\left[\psi(Y)|X=x\right]$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y|X}(y|x)\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)(1+\mathfrak{i}_{X;Y}(x,y))\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)[1+\operatorname{\zeta}_{\leq
k}(\mathfrak{i}_{X;Y})](x,y)\psi(y)+\sum_{y\in\mathcal{Y}}P_{Y}(y)\sum_{i>k}\sigma_{i}^{*}f_{i}^{*}(x)g_{i}^{*}(y)\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)(1+f^{\mathrm{T}}(x)g(y))\psi(y)$
$\displaystyle=\sum_{y\in\mathcal{Y}}P_{Y}(y)\psi(y)+\sum_{y\in\mathcal{Y}}P_{Y}(y)(\psi(y)g^{\mathrm{T}}(y))f^{\mathrm{T}}(x)$
$\displaystyle={\mathbb{E}}\left[\psi(Y)\right]+\Lambda_{\psi,g}\cdot f(x),$

where the fourth equality follows from the fact that

$\displaystyle\sum_{y\in\mathcal{Y}}P_{Y}(y)\sum_{i>k}\sigma_{i}^{*}f_{i}^{*}(x)g_{i}^{*}(y)\psi(y)=\sum_{i>k}\sigma_{i}^{*}f_{i}^{*}(x){\mathbb{E}}\left[g_{i}^{*}(Y)\psi(Y)\right]=0,$

due to $\psi_{j}\in\\{1,g_{1}^{*},\dots,g_{k}^{*}\\}$ for each $j\in[d]$.
$\blacksquare$

###  C.5 Proof of Property 3

The property is equivalent to
$\mathcal{L}(f,g)=\mathcal{L}(f+\underline{u},g+\underline{v})$ for all
$\underline{u},\underline{v}\in\mathbb{R}^{k}$. Hence, it suffices to prove
that

$\displaystyle\mathcal{L}(f+\underline{u},g+\underline{v})\leq\mathcal{L}(f,g),\qquad\text{for
all $\underline{u},\underline{v}\in\mathbb{R}^{k}$}.$ (127)

Note that for all $b\in\mathcal{F}_{\mathcal{Y}}$, since

$\displaystyle\left(f(x)+\underline{u}\right)^{\mathrm{T}}\left(g(y)+\underline{v}\right)+b(y)=f(x)\cdot
g(y)+\underline{v}^{\mathrm{T}}f(x)+\underline{u}^{\mathrm{T}}g(y)+b(y)+\underline{u}^{\mathrm{T}}\underline{v},$
(128)

we have [cf. (28)]
$\tilde{P}_{Y|X}^{(f+\underline{u},g+\underline{v},b)}=\tilde{P}_{Y|X}^{(f,g,b+\underline{u}^{\mathrm{T}}g)}$,
which implies that
$\mathcal{L}(f+\underline{u},g+\underline{v},b)=\mathcal{L}(f,g,b+\underline{u}^{\mathrm{T}}g)$.

Therefore, we obtain

$\displaystyle\mathcal{L}(f+\underline{u},g+\underline{v})=\max_{b\in\mathcal{F}_{\mathcal{Y}}}\mathcal{L}(f+\underline{u},g+\underline{v},b)\leq\max_{b\in\mathcal{F}_{\mathcal{Y}}}\mathcal{L}(f,g,b+\underline{u}^{\mathrm{T}}g)\leq\mathcal{L}(f,g).$
(129)

$\blacksquare$

###  C.6 Proof of Proposition 13

We first prove a useful lemma.

###### Lemma 37

Suppose $p>0,q>0$, $p+q=1$, then we have

$\displaystyle\log\left(p\cdot\exp\left(\frac{u}{p}\right)+q\cdot\exp\left(-\frac{u}{q}\right)\right)\geq\min\left\\{u^{2},u_{0}|u|\right\\},\quad\text{for
all $u\in\mathbb{R}$},$

where $u_{0}\triangleq\frac{\ln 2}{3}\cdot\min\\{p,q\\}$.

Proof [Proof of Lemma 37] Let $p_{\min}\triangleq\min\\{p,q\\}$.
$h(u)\triangleq\log\left(p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)\right).$
Then, we have

$\displaystyle
h^{\prime}(u)=\frac{\exp\left(p^{-1}u\right)-\exp\left(-q^{-1}u\right)}{p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)}$
$\displaystyle h^{\prime\prime}(u)$
$\displaystyle=\left[p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)\right]^{-2}$
$\displaystyle\quad\cdot\left[\left(\frac{1}{p}\exp\left(\frac{u}{p}\right)+\frac{1}{q}\exp\left(-\frac{u}{q}\right)\right)\cdot\left(p\exp\left(\frac{u}{p}\right)+q\exp\left(-\frac{u}{q}\right)\right)-\left[\exp\left(\frac{u}{p}\right)-\exp\left(-\frac{u}{q}\right)\right]^{2}\right]$
$\displaystyle\geq\frac{\left[\exp\left(p^{-1}u\right)+\exp\left(-q^{-1}u\right)\right]^{2}-\left[\exp\left(p^{-1}u\right)-\exp\left(-q^{-1}u\right)\right]^{2}}{\left[p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)\right]^{2}}$
$\displaystyle=4\cdot\frac{\exp\left((p^{-1}-q^{-1})\cdot
u\right)}{\left[p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)\right]^{2}}.$

Moreover, for all $|u|\leq u_{0}$, we have

$\displaystyle\exp\left((p^{-1}-q^{-1})\cdot
u\right)\geq\exp\left(-|p^{-1}-q^{-1}|\cdot
u_{0}\right)\geq\exp\left(-\frac{u_{0}}{p_{\min}}\right)$ $\displaystyle
p\cdot\exp\left(p^{-1}u\right)+q\cdot\exp\left(-q^{-1}u\right)\leq\exp\left(\frac{u_{0}}{p_{\min}}\right).$

As a result, for all $|u|\leq u_{0}$, we have $h^{\prime\prime}(u)\geq
4\exp\left(-\frac{3u_{0}}{p_{\min}}\right)=2$.

Therefore, $h^{\prime}(u_{0})\geq h^{\prime}(0)+2\cdot(u_{0}-0)=2u_{0}>u_{0}$,
and similarly, $-h^{\prime}(-u_{0})=h^{\prime}(0)-h^{\prime}(-u_{0})\geq
2u_{0}$, i.e., $h^{\prime}(-u_{0})\leq-2u_{0}\leq-u_{0}$. Moreover, for all
$|u|\leq u_{0}$, $h(u)\geq h(0)+h^{\prime}(0)\cdot u+\frac{1}{2}u^{2}\cdot
2=u^{2}.$ Therefore, for all $u>u_{0}$, $h^{\prime}(u)\geq
h^{\prime}(u_{0})>u_{0}$, which implies that $h(u)\geq
h(u_{0})+u_{0}(u-u_{0})\geq u_{0}^{2}+u_{0}u-u_{0}^{2}=u_{0}u$. Similarly, we
have $h(u)\geq-u_{0}u$ for all $u<-u_{0}$.
Proceeding to the proof of Proposition 13, we consider zero-mean
$k$-dimensional $f$, $g$. Without loss of generality, we assume
$b\in\mathcal{F}_{\mathcal{Y}}$ satisfies
${\mathbb{E}}\left[b(Y)\right]=-H(Y)$. Then, let
$a\in\mathcal{F}_{\mathcal{Y}|\varnothing}$ be $a(y)\triangleq b(y)-\log
P_{Y}(y)$, and define $\gamma\in\mathcal{F}_{\mathcal{X}\times\mathcal{Y}}$ as
$\gamma(x,y)\triangleq f(x)\cdot g(y)+a(y)$. Note that since

$\exp(f(x)\cdot g(y)+b(y))=P_{Y}(y)\exp(f(x)\cdot
g(y)+a(y))=P_{Y}(y)\exp(\gamma(x,y)),$

we have

$\displaystyle\tilde{P}^{(f,g,b)}_{Y|X}(y|x)=\frac{\exp(f(x)\cdot
g(y)+b(y))}{\sum_{y^{\prime}\in\mathcal{Y}}\exp(f(x)\cdot
g(y^{\prime})+b(y^{\prime}))}$
$\displaystyle=\frac{P_{Y}(y)\exp(\gamma(x,y))}{\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(x,y^{\prime}))}.$

Therefore,

$\displaystyle\mathcal{L}(f,g,b)$
$\displaystyle\quad={\mathbb{E}}_{(\hat{X},{\hat{Y}})\sim
P_{X}P_{Y}}\left[(\mathfrak{i}_{X;Y}(\hat{X},{\hat{Y}})+1)\cdot\log\tilde{P}^{(f,g,b)}_{Y|X}({\hat{Y}}|\hat{X})\right]$
$\displaystyle\quad={\mathbb{E}}_{(\hat{X},{\hat{Y}})\sim
P_{X}P_{Y}}\left[(\mathfrak{i}_{X;Y}(\hat{X},{\hat{Y}})+1)\cdot\left(\log
P_{Y}({\hat{Y}})+\gamma(\hat{X},{\hat{Y}})-\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right)\right]$
$\displaystyle\quad=-H(Y)+\langle\mathfrak{i}_{X;Y},\gamma\rangle-{\mathbb{E}}_{\hat{X}\sim
P_{X}}\left[\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right].$
(130)

As a result, for all $(f,g,b)$ that satisfies $\mathcal{L}(f,g,b)\geq-H(Y)$,
we have

$\displaystyle\langle\mathfrak{i}_{X;Y},\gamma\rangle\geq{\mathbb{E}}_{\hat{X}\sim
P_{X}}\left[\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right].$
(131)

In addition, note that for all $x\in\mathcal{X}$ and $y\in\mathcal{Y}$, we
have

$\displaystyle\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(x,y^{\prime}))$
$\displaystyle=P_{Y}(y)\exp(\gamma(x,y))+(1-P_{Y}(y))\sum_{y^{\prime}\in\mathcal{Y}\setminus\\{y\\}}\frac{P_{Y}(y^{\prime})}{1-P_{Y}(y)}\exp(\gamma(x,y^{\prime}))$
$\displaystyle\geq
P_{Y}(y)\exp(\gamma(x,y))+(1-P_{Y}(y))\exp\left(-\frac{P_{Y}(y)}{1-P_{Y}(y)}\cdot\gamma(x,y)\right).$

where the inequality follows from Jensen’s inequality and
$\sum_{y\in\mathcal{Y}}P_{Y}(y)\gamma(x,y)=0$.

Let us define

$\displaystyle q_{X}\triangleq\min_{x\in\mathcal{X}}P_{X}(x)>0,\quad
q_{Y}\triangleq\min_{y\in\mathcal{Y}}P_{Y}(y)>0.$

Then, from Lemma 37 we have

$\displaystyle\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(x,y^{\prime}))$
$\displaystyle\geq\log\left[P_{Y}(y)\exp(\gamma(x,y))+(1-P_{Y}(y))\exp\left(-\frac{P_{Y}(y)}{1-P_{Y}(y)}\cdot\gamma(x,y)\right)\right]$
$\displaystyle\geq\min\left\\{(P_{Y}(y)\gamma(x,y))^{2},\frac{\ln 2}{3}\cdot
q_{Y}|P_{Y}(y)\gamma(x,y)|\right\\}$ $\displaystyle\geq\frac{\ln 2\cdot
q_{Y}^{2}}{3}\cdot\min\left\\{(\gamma(x,y))^{2},|\gamma(x,y)|\right\\},$

which implies that

$\displaystyle{\mathbb{E}}_{\hat{X}\sim
P_{X}}\left[\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right]$
$\displaystyle\geq
P_{X}(x)\cdot\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(x,y^{\prime}))$
$\displaystyle\geq\frac{\ln 2\cdot
q_{X}q_{Y}^{2}}{3}\cdot\min\left\\{(\gamma(x,y))^{2},|\gamma(x,y)|\right\\}.$
(132)

On the other hand, since $\|\mathfrak{i}_{X;Y}\|=O(\epsilon)$, there exists a
constant $C>0$ such that $\|\mathfrak{i}_{X;Y}\|\leq C\epsilon$. Therefore,

$\displaystyle\langle\mathfrak{i}_{X;Y},\gamma\rangle\leq\|\mathfrak{i}_{X;Y}\|\cdot\|\gamma\|\leq\gamma_{\max}\cdot\epsilon,$
(133)

where
$\gamma_{\max}\triangleq\max_{x\in\mathcal{X},y\in\mathcal{Y}}|\gamma(x,y)|$.

Hence, by combining (131), (132) and (133), we obtain

$\displaystyle\frac{\ln 2\cdot
q_{X}q_{Y}^{2}}{3}\cdot\min\left\\{\gamma_{\max}^{2},\gamma_{\max}\right\\}\leq
C\cdot\gamma_{\max}\cdot\epsilon,$ (134)

where we have taken
$(x,y)=\operatorname*{arg\,max}_{x^{\prime},y^{\prime}}|\gamma(x^{\prime},y^{\prime})|$
in (132).

From (134), if $\epsilon<\frac{\ln 2}{3}\cdot\frac{q_{X}q_{Y}^{2}}{3C}$ we
have $\gamma_{\max}<\epsilon$. This implies that

$\displaystyle|\gamma(x,y)|<\epsilon,\quad\text{for all
}x\in\mathcal{X},y\in\mathcal{Y},$ (135)

and we obtain $\|\gamma\|<\epsilon$. In addition, note that since
$\|\gamma\|^{2}=\|f\otimes g+a\|^{2}=\|f\otimes g\|^{2}+\|a\|^{2}$, we obtain
$\|f\otimes g\|=O(\epsilon)$.

From (135), we have

$\displaystyle\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(x,y^{\prime}))$
$\displaystyle=\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\left(1+\gamma(x,y^{\prime})+\frac{(\gamma(x,y^{\prime}))^{2}}{2}+o(\epsilon^{2})\right)$
$\displaystyle=1+\frac{1}{2}\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})(\gamma(x,y^{\prime}))^{2}+o(\epsilon^{2})$
$\displaystyle=1+\frac{1}{2}\cdot{\mathbb{E}}_{{\hat{Y}}\sim
P_{Y}}\left[\gamma(x,{\hat{Y}})\right]^{2}+o(\epsilon^{2}).$

Therefore,

$\displaystyle{\mathbb{E}}_{\hat{X}\sim
P_{X}}\left[\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right]=\frac{1}{2}\cdot{\mathbb{E}}_{(\hat{X},{\hat{Y}})\sim
P_{X}P_{Y}}\left[\gamma(\hat{X},{\hat{Y}})\right]^{2}+o(\epsilon^{2})=\frac{1}{2}\cdot\|\gamma\|^{2}+o(\epsilon^{2}),$

and the likelihood (130) becomes

$\displaystyle\mathcal{L}(f,g,b)$
$\displaystyle=-H(Y)+\langle\mathfrak{i}_{X;Y},\gamma\rangle-{\mathbb{E}}_{\hat{X}\sim
P_{X}}\left[\log\sum_{y^{\prime}\in\mathcal{Y}}P_{Y}(y^{\prime})\exp(\gamma(\hat{X},y^{\prime}))\right]$
$\displaystyle=-H(Y)+\langle\mathfrak{i}_{X;Y},\gamma\rangle-\frac{1}{2}\cdot\|\gamma\|^{2}+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-\|\mathfrak{i}_{X;Y}-\gamma\|^{2}\right)-H(Y)+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-\|\mathfrak{i}_{X;Y}-f\otimes
g-a\|^{2}\right)-H(Y)+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-\|\mathfrak{i}_{X;Y}-f\otimes
g\|^{2}-\|a\|^{2}\right)-H(Y)+o(\epsilon^{2}),$ (136)

where the last equality follows from the fact that

$\displaystyle\|\mathfrak{i}_{X;Y}-f\otimes g-a\|^{2}$
$\displaystyle=\|\mathfrak{i}_{X;Y}-f\otimes g\|^{2}+\|a\|^{2},$

due to the orthogonality $(\mathfrak{i}_{X;Y}-f\otimes
g)\perp\mathcal{F}_{\mathcal{Y}}\ni a$.

Finally, from (136), for given $f,g$, $\mathcal{L}(f,g,b)$ is maximized when
$\|a\|=0$. Therefore, we have

$\displaystyle\mathcal{L}(f,g)=\max_{b\in\mathcal{F}_{\mathcal{Y}}}\mathcal{L}(f,g,b)$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-\|\mathfrak{i}_{X;Y}-f\otimes
g\|^{2}\right)-H(Y)+o(\epsilon^{2}),$

which gives (30). $\blacksquare$

###  C.7 Proof of Theorem 17

Throughout this proof, we consider

$\displaystyle
L_{1}(\bar{f},\bar{g},f,g)\triangleq\bigl{\|}\mathfrak{i}_{X;Y}-\bar{f}\otimes\bar{g}\bigr{\|}^{2},\quad
L_{2}(\bar{f},\bar{g},f,g)\triangleq\bigl{\|}\mathfrak{i}_{X;Y}-\bar{f}\otimes\bar{g}-f\otimes
g\bigr{\|}^{2}$

defined on the domain
$\mathcal{D}\triangleq\\{(\bar{f},\bar{g},f,g)\colon\bar{f}\in\mathcal{G}_{\mathcal{X}}^{\,{\bar{k}}},\bar{g}\in\mathcal{F}_{\mathcal{Y}}^{\,{\bar{k}}},f\in\mathcal{F}_{\mathcal{X}}^{\,k},g\in\mathcal{F}_{\mathcal{Y}}^{\,k}\\}$,
and let

$\displaystyle\mathcal{D}_{1}^{*}\triangleq\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{1}(\bar{f},\bar{g},f,g),\quad\mathcal{D}_{2}^{*}\triangleq\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{2}(\bar{f},\bar{g},f,g).$

We also define
$\bar{\gamma},\gamma\in\mathcal{F}_{\mathcal{X}\times\mathcal{Y}}$, as

$\displaystyle\bar{\gamma}\triangleq\Pi\left(\mathfrak{i}_{X;Y};{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{F}_{\mathcal{Y}}}\right),\quad\gamma\triangleq\mathfrak{i}_{X;Y}-\bar{\gamma}=\Pi\left(\mathfrak{i}_{X;Y};{(\mathcal{F}_{\mathcal{X}}\mathbin{\boxminus}\mathcal{G}_{\mathcal{X}})\otimes\mathcal{F}_{\mathcal{Y}}}\right).$
(137)

Then, it suffices to establish that

$\displaystyle\mathcal{D}^{*}\triangleq\mathcal{D}_{1}^{*}\cap\mathcal{D}_{2}^{*}=\\{(\bar{f},\bar{g},f,g)\in\mathcal{D}\colon\bar{f}\otimes\bar{g}=\bar{\gamma},f\otimes
g$ $\displaystyle=\operatorname{\zeta}_{\leq k}(\gamma)\\}.$ (138)

To see this, note that the common solution $\mathcal{D}^{*}$ (138) coincides
with the optimal solution (44) to be established. From Proposition 35, since
$\mathcal{D}^{*}\neq\varnothing$, we have

$\displaystyle\mathcal{D}^{*}=\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{1}(\bar{f},\bar{g},f,g)+L_{2}(\bar{f},\bar{g},f,g).$
(139)

Furthermore, from the definition of H-score (cf. Definition 10), we have

$\displaystyle\mathscr{H}\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g\end{bmatrix};\mathcal{C}_{\pi}\right)$
$\displaystyle=\mathscr{H}(\bar{f},\bar{g})+\mathscr{H}\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\ g\end{bmatrix}\right)$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-L_{1}(\bar{f},\bar{g},f,g)\right)+\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y}\|^{2}-L_{2}(\bar{f},\bar{g},f,g)\right)$
$\displaystyle=\|\mathfrak{i}_{X;Y}\|^{2}-\frac{1}{2}\left(L_{1}(\bar{f},\bar{g},f,g)+L_{2}(\bar{f},\bar{g},f,g)\right),$

which implies that

$\displaystyle\operatorname*{arg\,max}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}\mathscr{H}\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g\end{bmatrix};\mathcal{C}_{\pi}\right)=\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{1}(\bar{f},\bar{g},f,g)+L_{2}(\bar{f},\bar{g},f,g)=\mathcal{D}^{*}.$

It remains only to establish (138). From Proposition 36, it suffices to
establish that $t_{2}=t_{2}^{*}$ and (cf. (138))

$\displaystyle\mathcal{D}_{2}=\\{(\bar{f},\bar{g},f,g)\in\mathcal{D}\colon\bar{f}\otimes\bar{g}=\bar{\gamma},f\otimes
g=\gamma\\},$ (140)

where we have defined

$\displaystyle
t_{2}^{*}\triangleq\min_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{2}(\bar{f},\bar{g},f,g),\quad
t_{2}\triangleq\min_{(\bar{f},\bar{g},f,g)\in\mathcal{D}_{1}^{*}}L_{2}(\bar{f},\bar{g},f,g),\quad\mathcal{D}_{2}\triangleq\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}_{1}^{*}}L_{2}(\bar{f},\bar{g},f,g).$

To this end, let us define $\bar{h}\in\mathcal{G}_{\mathcal{X}}^{\,k}$ and
$\bar{h}\in(\mathcal{F}_{\mathcal{X}}\mathbin{\boxminus}\mathcal{G}_{\mathcal{X}})^{k}$
as $\bar{h}\triangleq\Pi\left(f;{\mathcal{G}_{\mathcal{X}}}\right)$, and
$h\triangleq
f-\bar{h}=\Pi\left(f;{\mathcal{F}_{\mathcal{X}}\mathbin{\boxminus}\mathcal{G}_{\mathcal{X}}}\right)$.
Then, we have

$\displaystyle L_{2}(\bar{f},\bar{g},f,g)$
$\displaystyle=\bigl{\|}\mathfrak{i}_{X;Y}-\bar{f}\otimes\bar{g}-f\otimes
g\bigr{\|}^{2}$ (141)
$\displaystyle=\bigl{\|}\bar{\gamma}+\gamma-\bar{f}\otimes\bar{g}-(\bar{h}+h)\otimes
g\bigr{\|}^{2}$ (142)
$\displaystyle=\bigl{\|}\bar{\gamma}-\bar{f}\otimes\bar{g}-\bar{h}\otimes
g\bigr{\|}^{2}+\bigl{\|}\gamma-h\otimes g\bigr{\|}^{2}$ (143)
$\displaystyle\geq\bigl{\|}\gamma-h\otimes g\bigr{\|}^{2}$ (144)
$\displaystyle\geq\bigl{\|}r_{k}(\gamma)\bigr{\|}^{2}$ (145)

where (143) follows from the orthogonality principle, since
$\left(\bar{\gamma}-\bar{f}\otimes\bar{g}-\bar{h}\otimes
g\right)\in\mathcal{G}_{\mathcal{X}}\otimes\mathcal{F}_{\mathcal{Y}}$ and
$\left(\gamma-h\otimes
g\right)\perp\mathcal{G}_{\mathcal{X}}\otimes\mathcal{F}_{\mathcal{Y}}$.
Moreover, it is easily verified that the lower bound (145) is tight: all
inequalities hold with equality for $(\bar{f},\bar{g},f,g)$ satisfying (44).

Therefore, we have

$\displaystyle
t_{2}^{*}=\min_{(\bar{f},\bar{g},f,g)\in\mathcal{D}}L_{2}(\bar{f},\bar{g},f,g)=\bigl{\|}r_{k}(\gamma)\bigr{\|}^{2}$
(146)

On the other hand, since ${\bar{k}}\geq\operatorname{rank}(\bar{\gamma})$,
from Proposition 5, we have
$\mathcal{D}_{1}^{*}=\\{(\bar{f},\bar{g},f,g)\colon\bar{f}\otimes\bar{g}=\bar{\gamma}\\}$.
Hence, for all $(\bar{f},\bar{g},f,g)\in\mathcal{D}_{1}^{*}$, we have

$\displaystyle
L_{2}(\bar{f},\bar{g},f,g)=\bigl{\|}\mathfrak{i}_{X;Y}-\bar{f}\otimes\bar{g}-f\otimes
g\bigr{\|}^{2}=\bigl{\|}\mathfrak{i}_{X;Y}-\bar{\gamma}-f\otimes
g\bigr{\|}^{2}=\bigl{\|}\gamma-f\otimes
g\bigr{\|}^{2}\geq\bigl{\|}r_{k}(\gamma)\bigr{\|}^{2},$

where the inequality holds with equality if and only if $f\otimes
g=\operatorname{\zeta}_{\leq k}(\gamma)$.

As a result,
$\mathcal{D}_{2}=\operatorname*{arg\,min}_{(\bar{f},\bar{g},f,g)\in\mathcal{D}_{1}^{*}}L_{2}(\bar{f},\bar{g},f,g)$
is given by (140), and we have

$\displaystyle
t_{2}=\min_{(\bar{f},\bar{g},f,g)\in\mathcal{D}_{1}^{*}}L_{2}(\bar{f},\bar{g},f,g)=\bigl{\|}r_{k}(\gamma)\bigr{\|}^{2}=t_{2}^{*}.$

$\blacksquare$

###  C.8 Proof of Theorem 18

To begin, we have

$\displaystyle\mathscr{H}\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g\end{bmatrix};\mathcal{C}_{\pi}^{\star}\right)$
$\displaystyle=\sum_{i=1}^{{\bar{k}}}\mathscr{H}(\bar{f}_{[i]},\bar{g}_{[i]})+\sum_{i=1}^{k}\mathscr{H}\left(\begin{bmatrix}\bar{f}\\\
f_{[i]}\end{bmatrix},\begin{bmatrix}\bar{g}\\\ g_{[i]}\end{bmatrix}\right)$
$\displaystyle=\sum_{i=1}^{{\bar{k}}-1}\mathscr{H}(\bar{f}_{[i]},\bar{g}_{[i]})+\sum_{i=1}^{k}\left(k^{-1}\cdot\mathscr{H}(\bar{f},\bar{g})+\mathscr{H}\left(\begin{bmatrix}\bar{f}\\\
f_{[i]}\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g_{[i]}\end{bmatrix}\right)\right).$

Note that for $\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g\end{bmatrix}\right)\in\operatorname{dom}(\mathcal{C}_{\pi}^{\star})$ and
each $i\in[{\bar{k}}]$, $\mathscr{H}(\bar{f}_{[i]},\bar{g}_{[i]})$ is
maximized if and only if

$\displaystyle\bar{f}_{[i]}\otimes\bar{g}_{[i]}=\operatorname{\zeta}_{\leq
i}(\mathfrak{i}_{X;Y}|\mathcal{G}_{\mathcal{X}},\mathcal{F}_{\mathcal{Y}}).$
(147)

In addition, from Theorem 17, for each $i\in[k]$, the term
$\left(k^{-1}\cdot\mathscr{H}(\bar{f},\bar{g})+\mathscr{H}\left(\begin{bmatrix}\bar{f}\\\
f_{[i]}\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g_{[i]}\end{bmatrix}\right)\right)$ is maximized if and only if

$\displaystyle\bar{f}\otimes\bar{g}=\Pi\left(\mathfrak{i}_{X;Y};{\mathcal{G}_{\mathcal{X}}\otimes\mathcal{F}_{\mathcal{Y}}}\right),\qquad
f_{[i]}\otimes g_{[i]}=\operatorname{\zeta}_{\leq
i}(\Pi\left(\mathfrak{i}_{X;Y};{(\mathcal{F}_{\mathcal{X}}\mathbin{\boxminus}\mathcal{G}_{\mathcal{X}})\otimes\mathcal{F}_{\mathcal{Y}}}\right)).$
(148)

Note that (45) is the common solution of (147) and (148). Hence, the proof is
completed by applying Proposition 35. $\blacksquare$

###  C.9 Proof of Proposition 20

The relation $\mathfrak{i}_{X;S}=\tilde{\ell}_{P^{\mathsf{M}}_{X,S,Y}}$ can be
directly verified from definition. To establish
$\pi_{\mathsf{M}}(\mathfrak{i}_{X;S,Y})=\Pi\left(\mathfrak{i}_{X;S,Y};{\mathcal{F}_{\mathcal{X}\times\mathcal{S}}}\right)=\mathfrak{i}_{X;S}$,
from Fact 2, it suffices to show that
$(\mathfrak{i}_{X;S,Y}-\mathfrak{i}_{X;S})\perp\mathcal{F}_{\mathcal{X}\times\mathcal{S}}$.

To this end, note that

$\displaystyle\mathfrak{i}_{X;S,Y}(x,s,y)-\mathfrak{i}_{X;S}(x,s)$
$\displaystyle=\frac{P_{X,S,Y}(x,s,y)-P^{\mathsf{M}}_{X,S,Y}(x,s,y)}{R_{X,S,Y}(x,s,y)}$
$\displaystyle=\frac{P_{X,S,Y}(x,s,y)-P_{X|S}(x|s)P_{S}(s)P_{Y|S}(y|s)}{R_{X,S,Y}(x,s,y)}.$

Therefore, for all $f\in\mathcal{F}_{\mathcal{X}\times\mathcal{S}}$, we have

$\displaystyle\bigl{\langle}\mathfrak{i}_{X;S,Y}-\mathfrak{i}_{X;S},f\bigr{\rangle}$
$\displaystyle=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{X,S,Y}(x,s,y)f(x,s)-\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{X|S}(x|s)P_{S}(s)P_{Y|S}(y|s)\cdot
f(x,s)$
$\displaystyle={\mathbb{E}}_{P_{X,S}}\left[f(X,S)\right]-{\mathbb{E}}_{P_{X,S}}\left[f(X,S)\right]$
$\displaystyle=0.$

$\blacksquare$

###  C.10 Proof of Proposition 21

Let $\gamma=\mathfrak{i}_{X;S,Y}^{(Q)}$, then the statement is equivalent to

$\displaystyle\gamma\in\mathcal{F}_{\mathcal{X}\times\mathcal{S}}\quad\iff\quad
Q_{X,S,Y}=Q_{X|S}Q_{S}Q_{Y|S}.$

##### “$\implies$”

If $\gamma\in\mathcal{F}_{\mathcal{X}\times\mathcal{S}}$, then we have

$\displaystyle Q_{X,S,Y}(x,s,y)$
$\displaystyle=Q_{X}(x)Q_{S,Y}(s,y)\left(1+\gamma(x,s)\right)$
$\displaystyle=P_{X}(x)P_{S,Y}(s,y)\left(1+\gamma(x,s)\right)$
$\displaystyle=P_{X}(x)P_{S}(s)\left(1+\gamma(x,s)\right)\cdot
P_{Y|S}(y|s),\quad\text{for all $x,s,y$}.$ (149)

Therefore,

$\displaystyle
Q_{X,S}(x,s)=\sum_{y\in\mathcal{Y}}Q_{X,S,Y}(x,s,y)=P_{X}(x)P_{S}(s)\left(1+\gamma(x,s)\right).$
(150)

From (149) and (150), we obtain

$\displaystyle
Q_{X,S,Y}(x,s,y)=Q_{X,S}(x,s)P_{Y|S}(y|s)=Q_{X,S}(x,s)Q_{Y|S}(y|s)=Q_{X|S}(x|s)Q_{S}(s)Q_{Y|S}(y|s).$
(151)

##### “$\impliedby$”

It suffices to note that

$\displaystyle\mathfrak{i}_{X;S,Y}^{(Q)}(x,s,y)=\frac{Q_{X,S,Y}(x,s,y)}{Q_{X}(x)Q_{S,Y}(s,y)}-1=\frac{Q_{X,S}(x,s)Q_{Y|S}(y|s)}{Q_{X}(x)Q_{S}(s)Q_{Y|S}(y|s)}-1=\frac{Q_{X,S}(x,s)}{Q_{X}(x)Q_{S}(s)}-1.$

$\blacksquare$

###  C.11 Proof of Proposition 23

The results on $\|\mathfrak{i}_{X;S}\|$ and $\|\mathfrak{i}_{X;Y|S}\|$
directly follows from Property 1. To establish (60), note that from
Proposition 20, we have

$\displaystyle
P_{X,S,Y}(x,s,y)-P_{X|S}(x|s)P_{S}(s)P_{Y|S}(y|s)=P_{X}(x)P_{S,Y}(s,y)\cdot\mathfrak{i}_{X;Y|S}(x,s,y),$

which implies that

$\displaystyle P_{Y|X,S}(y|x,s)$
$\displaystyle=P_{Y|S}(y|s)\cdot\left(1+\frac{P_{X}(x)P_{S}(s)}{P_{X,S}(x,s)}\cdot\mathfrak{i}_{X;Y|S}(x,s,y)\right)$
$\displaystyle=P_{Y|S}(y|s)\cdot\left(1+\frac{1}{1+\mathfrak{i}_{X;S}(x,s)}\cdot\mathfrak{i}_{X;Y|S}(x,s,y)\right)$
(152)
$\displaystyle=P_{Y|S}(y|s)\cdot\left(1+\frac{f^{\mathrm{T}}(x)g(s,y)}{1+\bar{f}^{\,\mathrm{T}}(x)\bar{g}(s)}\right),$
(153)

which further implies (61). $\blacksquare$

###  C.12 Proof of Proposition 24

When $X$ and $(S,Y)$ are $\epsilon$-dependent, we have
$\|\mathfrak{i}_{X;S,Y}\|=O(\epsilon)$, and thus

$\displaystyle\mathfrak{i}_{X;S}(x,s)=O(\epsilon),\quad\mathfrak{i}_{X;Y|S}(x,s,y)=O(\epsilon).$

Therefore, we have

$\displaystyle\frac{1}{1+\mathfrak{i}_{X;S}(x,s)}\cdot\mathfrak{i}_{X;Y|S}(x,s,y)=\left(1-\mathfrak{i}_{X;S}(x,s)\right)\cdot\mathfrak{i}_{X;Y|S}(x,s,y)+o(\epsilon)$
$\displaystyle=\mathfrak{i}_{X;Y|S}(x,s,y)+o(\epsilon).$

Then, it follows from (152) that

$\displaystyle
P_{Y|X,S}(y|x,s)=P_{Y|S}(y|s)\cdot\left(1+\mathfrak{i}_{X;Y|S}(x,s,y)\right)+o(\epsilon).$

Finally, the proof is completed by using the decomposition (63).
$\blacksquare$

###  C.13 Proof of Theorem 25

For each $s,x,y$, let us define $R_{X,Y}^{(s)}(x,y)\triangleq
P_{X|S=s}(x)P_{Y|S=s}(y)$ and

$\displaystyle\mathfrak{i}^{(s)}_{X;Y}(x,y)\triangleq\frac{P_{X,Y|S=s}(x,y)-P_{X|S=s}(x)P_{Y|S=s}(y)}{P_{X|S=s}(x)P_{Y|S=s}(y)}.$
(154)

In addition, we define $\underline{\mu}_{s}\in\mathbb{R}^{k}$ as
$\underline{\mu}_{s}\triangleq{\mathbb{E}}\left[f(X)|S=s\right]$ for each
$s\in\mathcal{S}$. Also, let
$\beta\in\mathcal{F}_{\mathcal{S}\times\mathcal{Y}}$ be defined as
$\beta(s,y)\triangleq\underline{\mu}_{s}^{\mathrm{T}}g(s,y)$.

Then, for each $s\in\mathcal{S}$, from (69) and Proposition 13 we have

$\displaystyle|f^{\mathrm{T}}(x)g(s,y)-\beta(s,y)|=|(f(x)-\underline{\mu}_{s})^{\mathrm{T}}g(s,y)|=O(\epsilon),\quad\text{for
all $x,s,y$},$ (155)

which implies that

$\displaystyle\left|f^{\mathrm{T}}(x)g(s,y)\right|=O(\epsilon),\quad\text{for
all $x,s,y$}.$ (156)

In addition, from Property 3, we have

$\displaystyle\mathcal{L}_{S}^{(s)}(f,g^{(s)})$
$\displaystyle\qquad=\mathcal{L}_{S}^{(s)}(f-\underline{\mu}_{s},g^{(s)})$
$\displaystyle\qquad=\frac{1}{2}\cdot\left(\bigl{\|}\mathfrak{i}^{(s)}_{X;Y}\bigr{\|}_{R_{X,Y}^{(s)}}^{2}-\bigl{\|}\mathfrak{i}^{(s)}_{X;Y}-(f-\underline{\mu}_{s})\otimes
g^{(s)}\bigr{\|}_{R_{X,Y}^{(s)}}^{2}\right)-H(Y|S=s)+o(\epsilon^{2})$
$\displaystyle\qquad=\bigl{\langle}\mathfrak{i}^{(s)}_{X;Y},f\otimes
g^{(s)}-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}\bigr{\rangle}_{R_{X,Y}^{(s)}}-\frac{1}{2}\cdot\bigl{\|}f\otimes
g^{(s)}-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}\bigr{\|}_{R_{X,Y}^{(s)}}^{2}-H(Y|S=s)+o(\epsilon^{2}),$
(157)

where $\langle\,\cdot\,,\,\cdot\,\rangle_{R}$ and $\|\,\cdot\,\|_{R}$ denote
the inner product and corresponding induced norm on the functional space, with
respect to the metric distribution $R$.

For the first two terms in (157), we compute their expectations over $P_{S}$.
For the first term,

$\displaystyle\sum_{s\in\mathcal{S}}P_{S}(s)\bigl{\langle}\mathfrak{i}^{(s)}_{X;Y},f\otimes
g^{(s)}-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}\bigr{\rangle}_{R_{X,Y}^{(s)}}$
$\displaystyle\qquad=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{S}(s)R_{X,Y}^{(s)}(x,y)\mathfrak{i}^{(s)}_{X;Y}(x,y)\left(f^{\mathrm{T}}(x)g^{(s)}(y)-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}(y)\right)$
$\displaystyle\qquad=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{X}(x)P_{S,Y}(s,y)\mathfrak{i}_{X;Y|S}(x,s,y)\cdot\left(f^{\mathrm{T}}(x)g(s,y)-\beta(s,y)\right)$
$\displaystyle\qquad=\bigl{\langle}\mathfrak{i}_{X;Y|S},f\otimes
g\bigr{\rangle}-\bigl{\langle}\mathfrak{i}_{X;Y|S},\beta\bigr{\rangle}$
$\displaystyle\qquad=\bigl{\langle}\mathfrak{i}_{X;Y|S},f\otimes
g\bigr{\rangle}$ (158)

where to obtain the second equality we have used the facts that

$\displaystyle\mathfrak{i}^{(s)}_{X;Y}(x,y)$
$\displaystyle=\frac{P_{X,Y|S=s}(x,y)-P_{X|S=s}(x)P_{Y|S=s}(y)}{P_{X|S=s}(x)P_{Y|S=s}(y)}$
$\displaystyle=\frac{P_{X,S,Y}(x,s,y)-P^{\mathsf{M}}_{X,S,Y}(x,s,y)}{P_{X}(x)P_{S,Y}(s,y)}\cdot\frac{P_{X}(x)}{P_{X|S=s}(x)}$
$\displaystyle=\mathfrak{i}_{X;Y|S}(x,s,y)\cdot\frac{1}{1+\mathfrak{i}_{X;S}(x,s)}$
(159)

and

$\displaystyle P_{S}(s)R_{X,Y}^{(s)}(x,y)=P^{\mathsf{M}}_{X,S,Y}(x,s,y)$
$\displaystyle=P_{X}(x)P_{S,Y}(s,y)\cdot(1+\mathfrak{i}_{X;S}(x,s)),$ (160)

and where to obtain (158) we have used the fact that
$\mathfrak{i}_{X;Y|S}\perp\mathcal{F}_{\mathcal{S}\times\mathcal{Y}}\ni\beta$.

For the second term of (157), we have

$\displaystyle\sum_{s\in\mathcal{S}}P_{S}(s)\bigl{\|}f\otimes
g^{(s)}-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}\bigr{\|}_{R_{X,Y}^{(s)}}^{2}$
$\displaystyle\qquad=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P^{\mathsf{M}}_{X,S,Y}(x,s,y)\cdot\left[(f(x)-\underline{\mu}_{s})^{\mathrm{T}}g(s,y)\right]^{2}$
(161)
$\displaystyle\qquad=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{X}(x)P_{S,Y}(s,y)\cdot\left[(f(x)-\underline{\mu}_{s})^{\mathrm{T}}g(s,y)\right]^{2}+o(\epsilon^{2})$
(162)
$\displaystyle\qquad=\sum_{x\in\mathcal{X},s\in\mathcal{S},y\in\mathcal{Y}}P_{X}(x)P_{S,Y}(s,y)\cdot\left[f^{\mathrm{T}}(x)g(s,y)-\beta(s,y)\right]^{2}+o(\epsilon^{2})$
(163) $\displaystyle\qquad=\left\|f\otimes
g-\beta\right\|^{2}+o(\epsilon^{2})$ (164)
$\displaystyle\qquad=\left\|f\otimes
g\right\|^{2}+\left\|\beta\right\|^{2}+o(\epsilon^{2}),$ (165)

where to obtain the second equality we have used (155), (160), and the fact
that $\|\mathfrak{i}_{X;S}\|=O(\epsilon)$. Furthermore, we can show that
$\|\beta\|=o(\epsilon)$. To see this, note that

$\displaystyle\beta(s,y)=\underline{\mu}_{s}^{\mathrm{T}}g(s,y)$
$\displaystyle=\sum_{x\in\mathcal{X}}P_{X|S}(x|s)f^{\mathrm{T}}(x)g(s,y)$
$\displaystyle=\sum_{x\in\mathcal{X}}P_{X}(x)\cdot[1+\mathfrak{i}_{X;S}(x,s)]\cdot
f^{\mathrm{T}}(x)g(s,y)$
$\displaystyle=\sum_{x\in\mathcal{X}}P_{X}(x)\cdot\mathfrak{i}_{X;S}(x,s)\cdot
f^{\mathrm{T}}(x)g(s,y).$

Then, from $\|\mathfrak{i}_{X;S}\|=O(\epsilon)$ and (156), we obtain
$|\beta(s,y)|=O(\epsilon^{2})$ and $\|\beta\|=O(\epsilon^{2})=o(\epsilon)$.

Therefore, we can refine (165) as

$\displaystyle\sum_{s\in\mathcal{S}}P_{S}(s)\bigl{\|}f\otimes
g^{(s)}-\underline{\mu}_{s}^{\mathrm{T}}g^{(s)}\bigr{\|}_{R_{X,Y}^{(s)}}^{2}=\left\|f\otimes
g\right\|^{2}+o(\epsilon^{2}).$ (166)

Combining (157), (158), and (166), we have

$\displaystyle\mathcal{L}_{S}(f,g)=\sum_{s\in\mathcal{S}}P_{S}(s)\cdot\mathcal{L}_{S}^{(s)}(f,g^{(s)})$
$\displaystyle=\left\langle\mathfrak{i}_{X;Y|S},f\otimes
g\right\rangle-\frac{1}{2}\cdot\bigl{\|}f\otimes
g\bigr{\|}^{2}-H(Y|S)+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\left(\|\mathfrak{i}_{X;Y|S}\|^{2}-\bigl{\|}\mathfrak{i}_{X;Y|S}-f\otimes
g\bigr{\|}^{2}\right)-H(Y|S)+o(\epsilon^{2}),$

which is maximized if and only if $f\otimes g=\operatorname{\zeta}_{\leq
k}(\mathfrak{i}_{X;Y|S})+o(\epsilon)$. $\blacksquare$

###  C.14 Proof of Proposition 26

Given any
$h\in\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}+\mathcal{F}_{\mathcal{X}_{2}\times\mathcal{Y}}$,
we can represent $h=h^{(1)}+h^{(2)}$ with
$h^{(i)}\in\mathcal{F}_{\mathcal{X}_{i}\times\mathcal{Y}},i=1,2$, i.e.,
$h(x_{1},x_{2},y)=h^{(1)}(x_{1},y)+h^{(2)}(x_{2},y)$.

Then, for any $Q_{X_{1},X_{2},Y}\in\mathcal{Q}_{\mathsf{B}}$, we have

$\displaystyle\langle\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y},h^{(i)}\rangle$
$\displaystyle=\sum_{x_{1}\in\mathcal{X}_{1},x_{2}\in\mathcal{X}_{2},y\in\mathcal{Y}}P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y}(x_{1},x_{2},y)\cdot
h^{(i)}(x_{i},y)$
$\displaystyle=\sum_{x_{1}\in\mathcal{X}_{1},x_{2}\in\mathcal{X}_{2},y\in\mathcal{Y}}[Q_{X_{1},X_{2},Y}(x_{1},x_{2},y)-P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)]\cdot
h^{(i)}(x_{i},y)$
$\displaystyle=\sum_{x_{i}\in\mathcal{X}_{i},y\in\mathcal{Y}}[P_{X_{i},Y}(x_{i},y)-P_{X_{i}}(x_{i})P_{Y}(y)]\cdot
h^{(i)}(x_{i},y)$
$\displaystyle=\langle\mathfrak{i}_{X_{i};Y},h^{(i)}\rangle\qquad\qquad\text{for
$i=1,2$.}$

As a result,

$\displaystyle\langle\pi_{\mathsf{B}}(\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y}),h\rangle$
$\displaystyle=\langle\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y}-\pi_{\mathsf{I}}(\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y}),h\rangle$
$\displaystyle=\langle\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y},h\rangle$
$\displaystyle=\langle\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y},h^{(1)}\rangle+\langle\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y},h^{(2)}\rangle\equiv\langle\mathfrak{i}_{X_{1};Y},h^{(1)}\rangle+\langle\mathfrak{i}_{X_{2};Y},h^{(2)}\rangle,$

where the second equality follows from the fact that
$\langle\pi_{\mathsf{I}}(\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y}),h\rangle=0$.

Hence, we obtain
$\langle\pi_{\mathsf{B}}(\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y})-\pi_{\mathsf{B}}(\tilde{\ell}_{P_{X_{1},X_{2},Y}}),h\rangle=0$
for all
$h\in\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}+\mathcal{F}_{\mathcal{X}_{2}\times\mathcal{Y}}$,
which implies that
$\pi_{\mathsf{B}}(\mathfrak{i}^{(Q)}_{X_{1},X_{2};Y})-\pi_{\mathsf{B}}(\tilde{\ell}_{P_{X_{1},X_{2},Y}})=0$.
$\blacksquare$

###  C.15 Proof of Proposition 28

From the definition, we have

$\displaystyle P_{X_{1},X_{2},Y}(x_{1},x_{2},y)$
$\displaystyle=P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)\left(1+[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)+[\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)\right)$
$\displaystyle=P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)\left[1+\bar{f}^{\,\mathrm{T}}(x_{1},x_{2})\bar{g}(y)+f^{\mathrm{T}}(x_{1},x_{2})g(y)\right],$

which implies (80a).

To obtain (80b), note that from (167), we have

$\displaystyle P_{X_{1},Y}(x_{1},y)$
$\displaystyle=\sum_{x_{2}\in\mathcal{X}_{2}}P^{\sf
B}_{X_{1},X_{2},Y}(x_{1},x_{2},y)$
$\displaystyle=P_{X_{1}}(x_{1})P_{Y}(y)\left[1+{\mathbb{E}}\left[\bar{f}^{\,\mathrm{T}}(x_{1},X_{2})\bar{g}(y)|X_{1}=x_{1}\right]\right]$
$\displaystyle=P_{X_{1}}(x_{1})P_{Y}(y)\left[1+\left(\bar{f}^{(1)}(x_{1})+[\tau_{{1}}(\bar{f}^{(2)})](x_{1})\right)^{\mathrm{T}}\bar{g}(y)\right],$

where to obtain the last equality we have used the fact that
$\tau_{{1}}(\bar{f})=\bar{f}^{(1)}+\tau_{{1}}(\bar{f}^{(2)})$. Similarly, we
can obtain (80b).

Finally, (81) can be readily obtained from (80). $\blacksquare$

###  C.16 Proof of Proposition 29

Note that since
$\mathfrak{i}_{X_{1},X_{2};Y}\in\mathcal{F}_{(\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y})|\varnothing}$,
we have
$\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\in\mathcal{F}_{(\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y})|\varnothing}$,
which implies that

$\displaystyle\sum_{(x_{1},x_{2},y)\in\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}}P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)\cdot[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)=\langle\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}),1\rangle=0.$

Therefore, it follows from the definition of $P^{\sf B}_{X_{1},X_{2},Y}$ that
$\displaystyle\sum_{(x_{1},x_{2},y)\in\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}}P^{\sf
B}_{X_{1},X_{2},Y}(x_{1},x_{2},y)=1$. Similarly, we have
$\displaystyle\sum_{(x_{1},x_{2},y)\in\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}}P^{\sf
I}_{X_{1},X_{2},Y}(x_{1},x_{2},y)=1$.

Moreover, from (82), we have $P^{\sf B}_{X_{1},X_{2},Y}(x_{1},x_{2},y)\geq
0,P^{\sf I}_{X_{1},X_{2},Y}(x_{1},x_{2},y)\geq 0$ for all $(x_{1},x_{2},y)$.
Therefore, we obtain $P^{\sf B}_{X_{1},X_{2},Y},P^{\sf
I}_{X_{1},X_{2},Y}\in\mathcal{P}^{\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}}$.

Since
$\mathfrak{i}_{X_{1},X_{2};Y}\perp\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{X}_{2}}$,
we have
$\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\perp\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{X}_{2}}$.
Therefore, for any
$(\hat{x}_{1},\hat{x}_{2})\in\mathcal{X}_{1}\times\mathcal{X}_{2}$, let
$f(x_{1},x_{2})=\delta_{x_{1}\hat{x}_{1}}\delta_{x_{2}\hat{x}_{2}}$, then we
have

$\displaystyle\sum_{y\in\mathcal{Y}}[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](\hat{x}_{1},\hat{x}_{2},y)$
$\displaystyle=\sum_{(x_{1},x_{2},y)\in\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}}[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)\cdot\delta_{x_{1}\hat{x}_{1}}\delta_{x_{2}\hat{x}_{2}}$
$\displaystyle=\langle\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}),f\rangle=0.$

This implies that $P^{\sf B}_{X_{1},X_{2}}=P_{X_{1},X_{2}}$. Similarly, we
have $P^{\sf I}_{X_{1},X_{2}}=P_{X_{1},X_{2}}$.

Finally, note that since

$\displaystyle
P_{X_{1},X_{2},Y}(x_{1},x_{2},y)=P_{X_{1},X_{2}}(x_{1},x_{2})P_{Y}(y)\left[1+[\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)+[\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2},y)\right],$

we have

$\displaystyle P_{X_{1},Y}(x_{1},y)-P^{\sf B}_{X_{1},Y}(x_{1},y)$
$\displaystyle=P^{\sf I}_{X_{1},Y}(x_{1},y)-P_{X_{1}}(x_{1})P_{Y}(y)$
$\displaystyle=\sum_{x_{2}^{\prime}\in\mathcal{X}_{2}}P_{X_{1},X_{2}}(x_{1},x_{2}^{\prime})P_{Y}(y)[\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})](x_{1},x_{2}^{\prime},y)$
$\displaystyle=0.$

To obtain the last equality, note that for any
$\hat{x}_{1},{\hat{y}}\in\mathcal{X}_{1}\times\mathcal{Y}$, let us define
$\gamma\in\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}$ as
$\gamma(x_{1},y)=\delta_{x_{1}\hat{x}_{1}}\delta_{y{\hat{y}}}$. Then, from
$\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})\perp\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}$,
we have

$\displaystyle
0=\langle\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y}),\gamma\rangle=\sum_{x_{2}^{\prime}\in\mathcal{X}_{2}}P_{X_{1},X_{2}}(\hat{x}_{1},x_{2}^{\prime})P_{Y}({\hat{y}})[\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})](\hat{x}_{1},x_{2}^{\prime},{\hat{y}}).$
(167)

Similarly, we can show that

$\displaystyle P_{X_{2},Y}(x_{2},y)-P^{\sf B}_{X_{2},Y}(x_{2},y)$
$\displaystyle=P^{\sf I}_{X_{2},Y}(x_{2},y)-P_{X_{2}}(x_{2})P_{Y}(y)=0.$

$\blacksquare$

###  C.17 Proof of Proposition 30

Note that since

$\displaystyle H(Q_{X_{1},X_{2},Y})$
$\displaystyle=H(P_{X_{1},X_{2}})+H(P_{Y})-I_{Q}(X_{1},X_{2};Y),$

we have

$\displaystyle
P^{{\mathsf{ent}}}_{X_{1},X_{2},Y}=\operatorname*{arg\,min}_{Q_{X_{1},X_{2},Y}\in\mathcal{Q}_{\mathsf{B}}}I_{Q}(X_{1},X_{2};Y),$
(168)

where $I_{Q}(X_{1},X_{2};Y)$ denotes the mutual information between
$(X_{1},X_{2})$ and $Y$ with respect to $Q_{X_{1},X_{2},Y}$.

Specifically, when we take $P_{X_{1},X_{2},Y}$ as the $Q_{X_{1},X_{2},Y}$, we
have $\displaystyle
I_{P}(X_{1},X_{2};Y)=\frac{1}{2}\cdot\|\mathfrak{i}_{X_{1},X_{2};Y}\|^{2}+o(\epsilon^{2}).$
Therefore, to solve (168), it suffices to consider $Q_{X_{1},X_{2},Y}$ with
$I_{Q}(X_{1},X_{2};Y)<I_{P}(X_{1},X_{2};Y)=O(\epsilon^{2})$. Since
$Q_{X_{1},X_{2}}Q_{Y}=P_{X_{1},X_{2}}P_{Y}\in\operatorname{relint}(\mathcal{P}^{\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{Y}})$,
we have $\|\tilde{\ell}_{Q_{X_{1},X_{2},Y}}\|^{2}=O(\epsilon^{2})$ [cf. (Sason
and Verdú, 2016, Eq. (338))].

Therefore,

$\displaystyle I_{Q}(X_{1},X_{2};Y)$
$\displaystyle=\frac{1}{2}\cdot\|\tilde{\ell}_{Q_{X_{1},X_{2},Y}}\|^{2}+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\|\pi_{\mathsf{B}}(\tilde{\ell}_{Q_{X_{1},X_{2},Y}})\|^{2}+\frac{1}{2}\cdot\|\pi_{\mathsf{I}}(\tilde{\ell}_{Q_{X_{1},X_{2},Y}})\|^{2}+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\|^{2}+\frac{1}{2}\cdot\|\pi_{\mathsf{I}}(\tilde{\ell}_{Q_{X_{1},X_{2},Y}})\|^{2}+o(\epsilon^{2})$
$\displaystyle\geq\frac{1}{2}\cdot\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\|^{2}+o(\epsilon^{2})$

where the inequality holds with equality when
$\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{Q_{X_{1},X_{2},Y}})\bigr{\|}=o(\epsilon)$.
Hence, we obtain
$\mathfrak{i}^{({\mathsf{ent}})}_{X_{1},X_{2};Y}=\tilde{\ell}_{P^{{\mathsf{ent}}}_{X_{1},X_{2},Y}}=\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})+o(\epsilon)$,
which gives (85). $\blacksquare$

###  C.18 Proof of Proposition 31

For all
$\phi\triangleq\phi^{(1)}+\phi^{(2)}\in\mathcal{F}_{\mathcal{X}_{1}}+\mathcal{F}_{\mathcal{X}_{2}}$
and $\psi\in\mathcal{F}_{\mathcal{Y}}$ with $\|\psi\|=1$, we have
${\mathbb{E}}_{P_{X_{1},X_{2}}P_{Y}}\left[\psi(Y)\phi(X_{1},X_{2})\right]=0$.
Hence,

$\displaystyle{\mathbb{E}}\left[\psi(Y)\phi(X_{1},X_{2})\right]$
$\displaystyle={\mathbb{E}}_{P_{X_{1},X_{2}}P_{Y}}\left[(1+\mathfrak{i}_{X_{1},X_{2};Y}(X_{1},X_{2},Y))\cdot\psi(Y)\phi(X_{1},X_{2})\right]$
$\displaystyle=\langle\mathfrak{i}_{X_{1},X_{2};Y},\phi\otimes\psi\rangle$
$\displaystyle=\langle\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}),\phi\otimes\psi\rangle+\langle\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y}),\phi\otimes\psi\rangle$
$\displaystyle=\langle\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}),\phi\otimes\psi\rangle,$

where the first equality follows from the definition of
$\mathfrak{i}_{X_{1},X_{2};Y}$, and where the last equality follow from the
fact that
$\pi_{\mathsf{I}}(\mathfrak{i}_{X_{1},X_{2};Y})\perp(\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}+\mathcal{F}_{\mathcal{X}_{2}\times\mathcal{Y}})\ni\phi\otimes\psi$.

Therefore, we can rewrite the objective (89) as

$\displaystyle{\mathbb{E}}\left[\left(\psi(Y)-\phi^{(1)}(X_{1})-\phi^{(2)}(X_{2})\right)^{2}\right]$
$\displaystyle={\mathbb{E}}\left[\left(\psi(Y)-\phi(X_{1},X_{2})\right)^{2}\right]$
$\displaystyle=\|\psi\|^{2}+\|\phi\|^{2}-2\cdot{\mathbb{E}}\left[\psi(Y)\phi(X_{1},X_{2})\right]$
$\displaystyle=1+\|\phi\otimes\psi\|^{2}-2\langle\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}),\phi\otimes\psi\rangle$
$\displaystyle=1+\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})-\phi\otimes\psi\|^{2}-\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\|^{2}.$

Finally, the proof is completed by noting that
$\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})\|^{2}=\sum_{i=1}^{{\bar{K}}}\bar{\sigma}_{i}^{2}$,
and we have

$\displaystyle\|\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y})-\phi\otimes\psi\|^{2}\geq\|{r_{1}(\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}))}\|^{2}=\sum_{i=2}^{{\bar{K}}}\bar{\sigma}_{i}^{2},$

where the inequality becomes equality if and only if
$\phi\otimes\psi=\operatorname{\zeta}(\pi_{\mathsf{B}}(\mathfrak{i}_{X_{1},X_{2};Y}))=\bar{\sigma}_{1}(\bar{f}_{1}^{*}\otimes\bar{g}_{1}^{*})$.

$\blacksquare$

###  C.19 Proof of Theorem 32

We first introduce a useful lemma.

###### Lemma 38

For all $k\geq
1,f\in\mathcal{F}_{\mathcal{X}}^{\,k},g\in\mathcal{F}_{\mathcal{Y}}^{\,k}$,
let
$\bar{h}\triangleq\Pi\left(f;{\mathcal{F}_{\mathcal{X}_{1}}+\mathcal{F}_{\mathcal{X}_{2}}}\right)$,
$h=f-\bar{h}$. Then, we have

$\displaystyle\mathscr{H}_{\rm{m}}(f,g)$
$\displaystyle=\frac{1}{2}\cdot\left[L(R_{X_{1},X_{2},Y})-\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}-L_{\sf B}(\bar{h}\otimes g)\right],$ (169)

where we have defined

$\displaystyle L_{\sf
B}(\gamma)\triangleq\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\pi_{\mathsf{B}}(\gamma)\bigr{\|}^{2}+\eta_{1}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(1)}_{X_{1},Y}}-\pi_{\mathsf{M}_{1}}(\gamma)\bigr{\|}^{2}+\eta_{2}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(2)}_{X_{2},Y}}-\pi_{\mathsf{M}_{2}}(\gamma)\bigr{\|}^{2}.$
(170)

Proof By definition, we have

$\displaystyle\mathscr{H}(f,g;\hat{P}^{(0)}_{X_{1},X_{2},Y})$
$\displaystyle=\frac{1}{2}\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\bigr{\|}^{2}-\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}-f\otimes
g\bigr{\|}^{2}\right)$
$\displaystyle=\frac{1}{2}\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\bigr{\|}^{2}-\bigl{\|}(\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g)+(\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\bar{h}\otimes
g)\bigr{\|}^{2}\right)$
$\displaystyle=\frac{1}{2}\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\bigr{\|}^{2}-\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}-\bigl{\|}\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\bar{h}\otimes
g\bigr{\|}^{2}\right),$ (171)

where to obtain the last equality we have used the orthogonality between
$\left(\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\right)\perp(\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}+\mathcal{F}_{\mathcal{X}_{2}\times\mathcal{Y}})$
and
$\left(\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\bar{h}\otimes
g\right)\in(\mathcal{F}_{\mathcal{X}_{1}\times\mathcal{Y}}+\mathcal{F}_{\mathcal{X}_{2}\times\mathcal{Y}})$.

In addition, for each $i=1,2$, from $\tau_{{i}}(f)=\tau_{{i}}(\bar{h})$ we
have

$\displaystyle\mathscr{H}(\tau_{{i}}(f),g;\hat{P}^{(i)}_{X_{i},Y})=\mathscr{H}(\tau_{{i}}(\bar{h}),g;\hat{P}^{(i)}_{X_{i},Y})$
$\displaystyle=\frac{1}{2}\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{i},Y}}\bigr{\|}^{2}-\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{i},Y}}-\tau_{{i}}(\bar{h})\otimes
g\bigr{\|}^{2}\right),$
$\displaystyle=\frac{1}{2}\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{i},Y}}\bigr{\|}^{2}-\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{i},Y}}-\pi_{\mathsf{M}_{i}}(\bar{h}\otimes
g)\bigr{\|}^{2}\right),$ (172)

where the last equality follows from that $\pi_{\mathsf{M}_{i}}(\bar{h}\otimes
g)=\tau_{{i}}(\bar{h})\otimes g$.

From (171) and (172), we can rewrite (92) as

$\displaystyle\mathscr{H}_{\rm{m}}(f,g)$
$\displaystyle=\eta_{0}\cdot\mathscr{H}(f,g;\hat{P}^{(0)}_{X_{1},X_{2},Y})+\eta_{1}\cdot\mathscr{H}\left(\tau_{{1}}(f),g;\hat{P}^{(1)}_{X_{1},Y}\right)+\eta_{2}\cdot\mathscr{H}\left(\tau_{{2}}(f),g;\hat{P}^{(2)}_{X_{2},Y}\right)$
$\displaystyle=\frac{\eta_{0}}{2}\cdot\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\bigr{\|}^{2}-\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}-\bigl{\|}\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\bar{h}\otimes
g\bigr{\|}^{2}\right)$
$\displaystyle\qquad+\sum_{i=1}^{2}\frac{\eta_{i}}{2}\cdot\left(\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{i},Y}}\bigr{\|}^{2}-\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{1},Y}}-\pi_{\mathsf{M}_{i}}(\bar{h}\otimes
g)\bigr{\|}^{2}\right)$
$\displaystyle=\frac{1}{2}\cdot\left[L(R_{X_{1},X_{2},Y})-\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}\right]$
$\displaystyle\qquad-\frac{1}{2}\left[\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{B}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-\bar{h}\otimes
g\bigr{\|}^{2}+\sum_{i=1}^{2}\eta_{i}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(i)}_{X_{1},Y}}-\pi_{\mathsf{M}_{i}}(\bar{h}\otimes
g)\bigr{\|}^{2}\right]$
$\displaystyle=\frac{1}{2}\cdot\left[L(R_{X_{1},X_{2},Y})-\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}-L_{\sf B}(\bar{h}\otimes g)\right],$

where the last equality follows from the fact that $\bar{h}\otimes
g=\pi_{\mathsf{B}}(\bar{h}\otimes g)$.

Let use define
$\bar{h}\triangleq\Pi\left(f;{\mathcal{F}_{\mathcal{X}_{1}}+\mathcal{F}_{\mathcal{X}_{2}}}\right)$
and $h\triangleq f-\bar{h}$. Then, from Lemma 38, we have

$\displaystyle\mathscr{H}_{\rm{m}}(\bar{f},\bar{g})=\frac{1}{2}\cdot\left[L(R_{X_{1},X_{2},Y})-\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})\bigr{\|}^{2}-L_{\sf
B}(\bar{f}\otimes\bar{g})\right],$ (173)

and, similarly,

$\displaystyle\mathscr{H}_{\rm{m}}\left(\begin{bmatrix}\bar{f}\,\\\
f\,\end{bmatrix},\begin{bmatrix}\bar{g}\\\
g\end{bmatrix}\right)=\frac{1}{2}\cdot\left[L(R_{X_{1},X_{2},Y})-\eta_{0}\cdot\bigl{\|}\pi_{\mathsf{I}}(\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}})-h\otimes
g\bigr{\|}^{2}-L_{\sf B}(\bar{f}\otimes\bar{g}+\bar{h}\otimes g)\right].$
(174)

Therefore, (173) is minimized if and only if

$\displaystyle\bar{f}\otimes\bar{g}=\pi_{\mathsf{B}}\Bigl{(}\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\Bigr{)},$
(175)

and (174) is minimized if and only if

$\displaystyle\bar{f}\otimes\bar{g}+\bar{h}\otimes
g=\pi_{\mathsf{B}}\Bigl{(}\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\Bigr{)},\qquad
h\otimes g=\operatorname{\zeta}_{\leq
k}\left(\pi_{\mathsf{I}}\Bigl{(}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\Bigr{)}\right).$
(176)

Hence, the common solution of (175) and (176) is

$\displaystyle\bar{f}\otimes\bar{g}=\pi_{\mathsf{B}}\Bigl{(}\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\Bigr{)},\qquad\bar{h}\otimes
g=0,\qquad h\otimes g=\operatorname{\zeta}_{\leq
k}\left(\pi_{\mathsf{I}}\Bigl{(}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\Bigr{)}\right),$
(177)

which is equivalent to (94). Finally, from Proposition 35, this is also the
solution that maximizes the nested H-score (93). $\blacksquare$

###  C.20 Proof of Theorem 33

We first introduce two useful facts.

######  Fact 6 (Cover and Thomas 2006, Theorem 11.1.2)

Given $m$ samples $Z_{1},\dots,Z_{m}$ i.i.d. generated from
$Q_{Z}\in\mathcal{P}^{\mathcal{Z}}$, th probability of observing
$\\{Z_{i}\\}_{i=1}^{m}=\\{z_{i}\\}_{i=1}^{m}$, denoted by
$\mathbb{P}\left\\{\\{z_{i}\\}_{i=1}^{m};Q_{Z}\right\\}$, satisfies
$\displaystyle-\log\mathbb{P}\left\\{\\{z_{i}\\}_{i=1}^{m};Q_{Z}\right\\}=m\left[H(\hat{P}_{Z})+D(\hat{P}_{Z}\|Q_{Z})\right]$,
where $\hat{P}_{Z}$ is the empirical distribution of $\\{z_{i}\\}_{i=1}^{m}$.

######  Fact 7 (Huang et al. 2019, Lemma 10)

Given a reference distribution
$R\in\operatorname{relint}(\mathcal{P}^{\mathcal{Z}})$. Then, for
$P,Q\in\mathcal{P}^{\mathcal{Z}}$ with $\|\tilde{\ell}_{P;R}\|=O(\epsilon)$,
$\|\tilde{\ell}_{Q;R}\|=O(\epsilon)$, we have

$\displaystyle
D(P\|Q)=\frac{1}{2}\cdot\|\tilde{\ell}_{P;R}-\tilde{\ell}_{Q;R}\|^{2}+o(\epsilon^{2}).$

Note that since the data from three different datasets are generated
independently, we have

$\displaystyle\mathbb{P}\left\\{\mathscr{D}_{0},\mathscr{D}_{1},\mathscr{D}_{2};Q_{X_{1},X_{2},Y}\right\\}$
$\displaystyle=\mathbb{P}\left\\{\mathscr{D}_{0};Q_{X_{1},X_{2},Y}\right\\}\cdot\mathbb{P}\left\\{\mathscr{D}_{1};Q_{X_{1},X_{2},Y}\right\\}\cdot\mathbb{P}\left\\{\mathscr{D}_{2};Q_{X_{1},X_{2},Y}\right\\}$
$\displaystyle=\mathbb{P}\left\\{\mathscr{D}_{0};Q_{X_{1},X_{2},Y}\right\\}\cdot\mathbb{P}\left\\{\mathscr{D}_{1};Q_{X_{1},Y}\right\\}\cdot\mathbb{P}\left\\{\mathscr{D}_{2};Q_{X_{2},Y}\right\\}.$

Therefore, from Fact 6,

$\displaystyle-\log\mathbb{P}\left\\{\mathscr{D}_{0},\mathscr{D}_{1},\mathscr{D}_{2};Q_{X_{1},X_{2},Y}\right\\}$
$\displaystyle\qquad=-\log\mathbb{P}\left\\{\mathscr{D}_{0};Q_{X_{1},X_{2},Y}\right\\}-\log\mathbb{P}\left\\{\mathscr{D}_{1};Q_{X_{1},Y}\right\\}-\log\mathbb{P}\left\\{\mathscr{D}_{2};Q_{X_{2},Y}\right\\}$
$\displaystyle\qquad=n_{0}\cdot
D\bigl{(}\hat{P}^{(0)}_{X_{1},X_{2},Y}\big{\|}Q_{X_{1},X_{2},Y}\bigr{)}+n_{1}\cdot
D\bigl{(}\hat{P}_{X_{1},Y}^{(1)}\big{\|}Q_{X_{1},Y}\bigr{)}+n_{2}\cdot
D\bigl{(}\hat{P}_{X_{2},Y}^{(2)}\big{\|}Q_{X_{2},Y}\bigr{)}$
$\displaystyle\qquad\qquad+n_{0}\cdot
H(\hat{P}^{(0)}_{X_{1},X_{2},Y})+n_{1}\cdot
H(\hat{P}_{X_{1},Y}^{(1)})+n_{2}\cdot H(\hat{P}_{X_{2},Y}^{(2)})$
$\displaystyle\qquad=n\cdot L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})+n_{0}\cdot
H(\hat{P}^{(0)}_{X_{1},X_{2},Y})+n_{1}\cdot
H(\hat{P}_{X_{1},Y}^{(1)})+n_{2}\cdot H(\hat{P}_{X_{2},Y}^{(2)})$

where the second equality follows from (6), and where we have defined

$\displaystyle L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})\triangleq\eta_{0}\cdot
D\bigl{(}\hat{P}^{(0)}_{X_{1},X_{2},Y}\big{\|}Q_{X_{1},X_{2},Y}\bigr{)}+\eta_{1}\cdot
D\bigl{(}\hat{P}_{X_{1},Y}^{(1)}\big{\|}Q_{X_{1},Y}\bigr{)}+\eta_{2}\cdot
D\bigl{(}\hat{P}_{X_{2},Y}^{(2)}\big{\|}Q_{X_{2},Y}\bigr{)}.$

Hence, we can rewrite the maximum likelihood solution
$P^{(\mathrm{ML})}_{X_{1},X_{2},Y}$ as

$\displaystyle P^{(\mathrm{ML})}_{X_{1},X_{2},Y}$
$\displaystyle=\operatorname*{arg\,max}_{Q_{X_{1},X_{2},Y}}\leavevmode\nobreak\
\mathbb{P}\left\\{\mathscr{D}_{0},\mathscr{D}_{1},\mathscr{D}_{2};Q_{X_{1},X_{2},Y}\right\\}=\operatorname*{arg\,min}_{Q_{X_{1},X_{2},Y}}L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y}).$

From $L(R_{X_{1},X_{2},Y})=O(\epsilon^{2})$, we obtain
$\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}\bigr{\|}=O(\epsilon),\bigl{\|}\tilde{\ell}_{\hat{P}^{(1)}_{X_{1},Y}}\bigr{\|}=O(\epsilon),\bigl{\|}\tilde{\ell}_{\hat{P}^{(2)}_{X_{2},Y}}\bigr{\|}=O(\epsilon)$.
By definition, we have
$L({P^{(\mathrm{est})}_{X_{1},X_{2},Y}})<L(R_{X_{1},X_{2},Y})=O(\epsilon^{2})$,
which implies that
$\bigl{\|}\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}=O(\epsilon)$.

We first consider $Q_{X_{1},X_{2},Y}$ with

$\displaystyle\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}\leq\epsilon.$
(178)

Then, we have

$\displaystyle\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}\bigr{\|}\leq\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}+\bigl{\|}\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}=O(\epsilon),$
(179)

and it follows from Fact 7 that

$\displaystyle L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})$
$\displaystyle=\eta_{0}\cdot
D\bigl{(}\hat{P}^{(0)}_{X_{1},X_{2},Y}\big{\|}Q_{X_{1},X_{2},Y}\bigr{)}+\eta_{1}\cdot
D\bigl{(}\hat{P}_{X_{1},Y}^{(1)}\big{\|}Q_{X_{1},Y}\bigr{)}+\eta_{2}\cdot
D\bigl{(}\hat{P}_{X_{2},Y}^{(2)}\big{\|}Q_{X_{2},Y}\bigr{)}$
$\displaystyle=\frac{1}{2}\Bigl{(}\eta_{0}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(0)}_{X_{1},X_{2},Y}}-\tilde{\ell}_{Q_{X_{1},X_{2},Y}}\bigr{\|}^{2}+\eta_{1}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(1)}_{X_{1},Y}}-\tilde{\ell}_{Q_{X_{1},Y}}\bigr{\|}^{2}$
$\displaystyle\qquad\quad+\eta_{2}\cdot\bigl{\|}\tilde{\ell}_{\hat{P}^{(2)}_{X_{2},Y}}-\tilde{\ell}_{Q_{X_{2},Y}}\bigr{\|}^{2}\Bigr{)}+o(\epsilon^{2})$
$\displaystyle=\frac{1}{2}\cdot L({Q_{X_{1},X_{2},Y}})+o(\epsilon^{2}).$ (180)

Therefore, for $Q_{X_{1},X_{2},Y}$ that satisfies (178), the minimum value of
$L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})$ is achieved by
$Q_{X_{1},X_{2},Y}=P^{(\mathrm{est})}_{X_{1},X_{2},Y}+o(\epsilon)$.

Now we consider the case of $Q_{X_{1},X_{2},Y}$ with
$\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}>\epsilon$.
Let
$\epsilon^{\prime}=\epsilon/\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}<1$
and define

$\displaystyle\bar{Q}_{X_{1},X_{2},Y}\triangleq\epsilon^{\prime}\cdot
Q_{X_{1},X_{2},Y}+(1-\epsilon^{\prime})\cdot
P^{(\mathrm{est})}_{X_{1},X_{2},Y}.$ (181)

Then, we have
$\tilde{\ell}_{\bar{Q}_{X_{1},X_{2},Y}}=\epsilon^{\prime}\cdot\tilde{\ell}_{Q_{X_{1},X_{2},Y}}+(1-\epsilon^{\prime})\cdot\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}$,
which implies

$\displaystyle\bigl{\|}\tilde{\ell}_{\bar{Q}_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}=\epsilon^{\prime}\cdot\bigl{\|}\tilde{\ell}_{Q_{X_{1},X_{2},Y}}-\tilde{\ell}_{P^{(\mathrm{est})}_{X_{1},X_{2},Y}}\bigr{\|}=\epsilon.$
(182)

As a result, we can apply the same analysis on $\bar{Q}_{X_{1},X_{2},Y}$ and
obtain [cf. (180)]

$\displaystyle L^{(\mathrm{ML})}(\bar{Q}_{X_{1},X_{2},Y})=\frac{1}{2}\cdot
L({\bar{Q}_{X_{1},X_{2},Y}})+o(\epsilon^{2}).$ (183)

Hence, we obtain from (182) that

$\displaystyle
L^{(\mathrm{ML})}(\bar{Q}_{X_{1},X_{2},Y})>L^{(\mathrm{ML})}(P^{(\mathrm{est})}_{X_{1},X_{2},Y})=\frac{1}{2}\cdot
L({P^{(\mathrm{est})}_{X_{1},X_{2},Y}})+o(\epsilon^{2})$ (184)

for $\epsilon$ sufficiently small.

In addition, since $L^{(\mathrm{ML})}$ is convex, it follows from Jensen’s
inequality that (cf. (181))

$\displaystyle
L^{(\mathrm{ML})}(\bar{Q}_{X_{1},X_{2},Y})\leq\epsilon^{\prime}\cdot
L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})+(1-\epsilon^{\prime})\cdot
L^{(\mathrm{ML})}(P^{(\mathrm{est})}_{X_{1},X_{2},Y}).$ (185)

As a result, from (184) and (185) we have
$L^{(\mathrm{ML})}(Q_{X_{1},X_{2},Y})>L^{(\mathrm{ML})}(P^{(\mathrm{est})}_{X_{1},X_{2},Y})$.

Combining both cases of $Q_{X_{1},X_{2},Y}$, we obtain (96). $\blacksquare$

###  C.21 Proof of Proposition 34

It suffices to establish that there exist
$\alpha\colon\mathcal{U}\to\mathbb{R}$, $\beta\colon\mathcal{V}\to\mathbb{R}$,
such that

$\displaystyle\mathfrak{i}_{\underline{X};U}(\underline{x},u)$
$\displaystyle=\alpha(u)\cdot\left[\tanh\left(2w\cdot\varphi(\underline{x})+b_{U}\right)-\tanh(b_{U})\right],$
(186a) $\displaystyle\mathfrak{i}_{\underline{Y};V}(\underline{y},v)$
$\displaystyle=\beta(v)\cdot\left[\tanh\left(2w\cdot\varphi(\underline{y})+b_{V}\right)-\tanh(b_{V})\right].$
(186b)

To see this, note that from the Markov relation
$\underline{X}-U-V-\underline{Y}$, we have

$\displaystyle P_{\underline{X},\underline{Y}}(\underline{x},\underline{y})$
$\displaystyle=\sum_{u\in\mathcal{U},v\in\mathcal{V}}P_{\underline{X}|U}(\underline{x}|u)\cdot
P_{\underline{Y}|V}(\underline{y}|v)\cdot P_{UV}(u,v)$
$\displaystyle=P_{\underline{X}}(\underline{x})P_{\underline{Y}}(\underline{y})\sum_{u\in\mathcal{U},v\in\mathcal{V}}P_{U,V}(u,v)\cdot(1+\mathfrak{i}_{\underline{X};U}(\underline{x},u))\cdot(1+\mathfrak{i}_{\underline{Y};V}(\underline{y},v))$
$\displaystyle=P_{\underline{X}}(\underline{x})P_{\underline{Y}}(\underline{y})\left(1+\sum_{u\in\mathcal{U},v\in\mathcal{V}}P_{U,V}(u,v)\cdot\mathfrak{i}_{\underline{X};U}(\underline{x},u)\cdot\mathfrak{i}_{\underline{Y};V}(\underline{y},v)\right),$

where to obtain the last equality we have used the fact that

$\displaystyle\sum_{u\in\mathcal{U}}P_{U}(u)\cdot\mathfrak{i}_{\underline{X};U}(\underline{x},u)=\sum_{v\in\mathcal{V}}P_{V}(v)\cdot\mathfrak{i}_{\underline{Y};V}(\underline{y},v)=0.$

Therefore, we obtain

$\displaystyle\mathfrak{i}_{\underline{X};\underline{Y}}(\underline{x},\underline{y})$
$\displaystyle=\frac{P_{\underline{X},\underline{Y}}(\underline{x},\underline{y})}{P_{\underline{X}}(\underline{x})P_{\underline{Y}}(\underline{y})}-1$
$\displaystyle=\sum_{u\in\mathcal{U},v\in\mathcal{V}}P_{U,V}(u,v)\cdot\mathfrak{i}_{\underline{X};U}(\underline{x},u)\cdot\mathfrak{i}_{\underline{Y};V}(\underline{y},v)$
$\displaystyle={\mathbb{E}}\left[\alpha(U)\beta(V)\right]\cdot\left[\tanh\left(2w\cdot\varphi(\underline{x})+b_{U}\right)-\tanh(b_{U})\right]\cdot\left[\tanh\left(2w\cdot\varphi(\underline{y})+b_{V}\right)-\tanh(b_{V})\right],$

which gives (102).

It remains only to establish (186). For symmetry, we consider only (186a). To
begin, by definition, we have

$\displaystyle P_{\underline{X}|U=u}(\underline{x})$
$\displaystyle=\frac{1}{2}\cdot\prod_{i=1}^{l-1}\left[q_{u}^{(1-\delta_{x_{i}x_{i+1}})}(1-q_{u})^{\delta_{x_{i}x_{i+1}}}\right],$

which implies

$\displaystyle\log P_{\underline{X}|U=u}(\underline{x})=\log\frac{1}{2}+\log
q_{u}\cdot\sum_{i=1}^{l-1}(1-\delta_{x_{i}x_{i+1}})+\log(1-q_{u})\cdot\sum_{i=1}^{l-1}\delta_{x_{i}x_{i+1}}.$

Therefore, we obtain

$\displaystyle\frac{1}{2}\log\frac{P_{\underline{X}|U=1}(\underline{x})}{P_{\underline{X}|U=0}(\underline{x})}$
$\displaystyle=\frac{1}{2}\left[\log{P_{\underline{X}|U=1}(\underline{x})}-\log{P_{\underline{X}|U=0}(\underline{x})}\right]$
$\displaystyle=\frac{1}{2}\left[\log\frac{q_{1}}{q_{0}}\cdot\sum_{i=1}^{l-1}(1-\delta_{x_{i}x_{i+1}})+\log\frac{1-q_{1}}{1-q_{0}}\cdot\sum_{i=1}^{l-1}\delta_{x_{i}x_{i+1}}\right]$
$\displaystyle=\log\frac{q_{1}}{q_{0}}\cdot\left(\frac{l-1}{2}-\sum_{i=1}^{l-1}\delta_{x_{i}x_{i+1}}\right)$
$\displaystyle=2w\cdot\varphi(\underline{x}).$

As a consequence,

$\displaystyle\frac{P_{\underline{X}|U=1}(\underline{x})-P_{\underline{X}|U=0}(\underline{x})}{P_{\underline{X}}(\underline{x})}$
$\displaystyle=\frac{P_{\underline{X}|U=1}(\underline{x})-P_{\underline{X}|U=0}(\underline{x})}{P_{U}(1)\cdot
P_{\underline{X}|U=1}(\underline{x})+P_{U}(0)\cdot
P_{\underline{X}|U=0}(\underline{x})}$
$\displaystyle=\frac{1}{P_{U}(0)}\cdot\frac{\frac{P_{\underline{X}|U=1}(\underline{x})}{P_{\underline{X}|U=0}(\underline{x})}-1}{\frac{P_{U}(1)}{P_{U}(0)}\cdot\frac{P_{\underline{X}|U=1}(\underline{x})}{P_{\underline{X}|U=0}(\underline{x})}+1}$
$\displaystyle=\frac{1}{2P_{U}(0)P_{U}(1)}\cdot\left[\frac{\frac{P_{U}(1)}{P_{U}(0)}\cdot\frac{P_{\underline{X}|U=1}(\underline{x})}{P_{\underline{X}|U=0}(\underline{x})}-1}{\frac{P_{U}(1)}{P_{U}(0)}\cdot\frac{P_{\underline{X}|U=1}(\underline{x})}{P_{\underline{X}|U=0}(\underline{x})}+1}-\frac{\frac{P_{U}(1)}{P_{U}(0)}-1}{\frac{P_{U}(1)}{P_{U}(0)}+1}\right]$
$\displaystyle=\frac{1}{2P_{U}(0)P_{U}(1)}\cdot\left[\tanh\left(2w\cdot\varphi(\underline{x})+b_{U}\right)-\tanh(b_{U})\right],$

where to obtain the last equality we have used the fact that
$\frac{t-1}{t+1}=\tanh\left(\frac{1}{2}\log t\right)$.

Hence, with $u^{\prime}=1-u$, we have

$\displaystyle\mathfrak{i}_{\underline{X};U}(\underline{x},u)$
$\displaystyle=\frac{P_{\underline{X},U}(\underline{x},u)-P_{\underline{X}}(\underline{x})P_{U}(u)}{P_{\underline{X}}(\underline{x})P_{U}(u)}$
$\displaystyle=\frac{P_{\underline{X}|U=u}(\underline{x})-P_{\underline{X}}(\underline{x})}{P_{\underline{X}}(\underline{x})}$
$\displaystyle=\frac{P_{\underline{X}|U=u}(\underline{x})-P_{U}(u)P_{\underline{X}|U=u}(\underline{x})-P_{U}(u^{\prime})P_{\underline{X}|U=u^{\prime}}(\underline{x})}{P_{\underline{X}}(\underline{x})}$
$\displaystyle=P_{U}(u^{\prime})\cdot\frac{P_{\underline{X}|U=u}(\underline{x})-P_{\underline{X}|U=u^{\prime}}(\underline{x})}{P_{\underline{X}}(\underline{x})}$
$\displaystyle=P_{U}(u^{\prime})\cdot(-1)^{u+1}\cdot\frac{P_{\underline{X}|U=1}(\underline{x})-P_{\underline{X}|U=0}(\underline{x})}{P_{\underline{X}}(\underline{x})}$
$\displaystyle=\frac{(-1)^{u+1}}{2P_{U}(u)}\cdot\left[\tanh\left(2w\cdot\varphi(\underline{x})+b_{U}\right)-\tanh(b_{U})\right],$

which gives (186a) as desired, with $\alpha(u)=\frac{(-1)^{u+1}}{2P_{U}(u)}$.
$\blacksquare$

##  Appendix Appendix D Implementation Details of Experiments

We implement our experiments in Python 3 (Van Rossum and Drake, 2009), where
we use the PyTorch (Paszke et al., 2019) library for neural network training
and use the Matplotlib (Hunter, 2007) library for plotting. We also make use
of NumPy (Harris et al., 2020) and SciPy (Virtanen et al., 2020) for the
computation.

In the experiments, we apply Adam (Kingma and Ba, 2015) as the optimizer with
the default parameters: a learning rate of $10^{-3}$,
$\beta_{1}=0.9,\beta_{2}=0.999$, and $\epsilon=10^{-8}$. For each MLP
(multilayer perceptron) used in the experiments, we set the activation
function to be the softplus function $x\mapsto\log(1+e^{x})$, which are
applied to all layers except the output layer.

It is worth mentioning that our choices of network architectures, optimizers
and hyperparameters are not tuned, and thus might not be optimal. It is
possible to further optimize such choices to improve the performance or
convergence.

###  D.1 Learning Maximal Correlation Functions

We first introduce the implementation details for Section 7.1, where the goal
is to learn maximal correlation functions for different data. The
corresponding learning objective is the nested H-score (35), which are
maximized during the training.

####  D.1.1 Implementation of Section 7.1.1

We set $|\mathcal{X}|=8,|\mathcal{Y}|=6$, and feature dimension $k=3$. To
generate the discrete distributions $P_{X,Y}$, we draw
$(|\mathcal{X}|\cdot|\mathcal{Y}|)$ i.i.d. numbers from $\mathrm{Unif}[0,1]$
and divide each number by their sum. We then use the resulting
$|\mathcal{X}|\times|\mathcal{Y}|$ table as the values for the probability
mass function $P_{X,Y}$. To ensure reproducible results, we set the random
seed of NumPy to $20\,230\,606$ in the generating process.

Then, we generate $N=30\,000$ training sample pairs of $(X,Y)$ from $P_{X,Y}$,
then apply one-hot encoding such that the inputs are represented as
$|\mathcal{X}|$ and $|\mathcal{Y}|$ dimensional vectors. Then, we use two one-
layer linear networks as the feature extractors $f$ and $g$.

We train the networks with a minibatch size of 128 for 100 epochs. Then, we
obtain the estimated $f_{i}^{*}$, $g_{i}^{*}$, and $\sigma_{i}$ by applying
(37) and compare them with corresponding theoretical values, which we compute
from the SVD of corresponding CDM matrix [cf. (113)], with the results shown
in Figure 13. Note that since $f_{i}^{*}\otimes
g_{i}^{*}=(-f_{i}^{*})\otimes(-g_{i}^{*})$, both $(f_{i}^{*},g_{i}^{*})$ and
$(-f_{i}^{*},-g_{i}^{*})$ are the optimal feature pairs. For the sake of
presentation, we applied a sign modification before the plotting.

####  D.1.2 Implementation of Section 7.1.2

In this experiment, we first generate $N=50\,000$ samples of
$(X,Y)\in\mathbb{R}^{2}$ for training, to learn $k=2$ dimensional features $f$
and $g$. We use two MLPs of the same architecture as the feature extractors
for $f$ and $g$. Specifically, each MLP is with three layers, where the
dimensions for all intermediate features, from input to output, are:
$\text{input}=1$ – 32 – 32 – $2=\text{output}$. We then train the networks
with a minibatch size of 256 for 100 epochs and use the learned features for
estimation tasks, as demonstrated in Section 7.1.2.

####  D.1.3 Implementation of Section 7.1.3

In this experiment, we set $k=1$. To extract $f$ and $g$ from input sequences
$\underline{X}$ and $\underline{Y}$, we use one-dimensional convolutional
neural networks as the feature extractors, which are used in sentence
classifications (Kim, 2014; Zhang and Wallace, 2017). In particular, $f$ and
$g$ are of the same architecture, composed of an embedding (linear) layer, a 1
dimensional convolutional layer, an average pooling layer, and a fully
connected (linear) layer.

We use feature extractor $f$ as an example to illustrate the processing of
sequential data. First, we represent $\underline{x}$ sequence as a one-hot
encoded list, i.e., each
$x_{i}\in\\{(1,0)^{\mathrm{T}},(0,1)^{\mathrm{T}}\\}$. Then, the embedding
layer maps each $x_{i}$ to a 4-dimensional vector. The one-dimensional
convolutional layer then processes the length-$l$ list of embedded
4-dimensional vectors, by 32 convolutional kernels of size 4. We then activate
the convolution results by the ReLU function $x\mapsto\max\\{x,0\\}$. The
output from each convolutional kernel is further averaged by the average
pooling layer, leading to a 32 dimensional feature, with each dimension
corresponding to a convolutional kernel. Finally, we feed the 32 dimensional
feature to the fully connected layer and generate $k=1$ dimensional output.

Then, we train the feature extractors $f$ and $g$ with a minibatch size of 128
for 100 epochs. The learned features are shown in Section 7.1.3.

###  D.2 Learning With Orthogonality Constraints

In this experiment, we use the same dataset generated in Appendix D.1.2. We
set ${\bar{k}}=k=1$, i.e., we learn one-dimensional feature $f$ from $X$
orthogonal to given one-dimensional $\bar{f}$. To this end, we use three MLPs
of the same architecture as the feature extractors for $\bar{g},f,g$, with
dimensions $\text{input}=1$ – 32 – 32 – $1=\text{output}$. We then train the
networks with a minibatch size of 256 for 100 epochs to maximize the nested
H-score restricted to $\bar{f}=\phi$ [cf. (49)], for $\phi(x)=x$ and
$\phi(x)=x^{2}$, respectively.

###  D.3 Learning With Side Information

We set $|\mathcal{X}|=8,|\mathcal{S}|=|\mathcal{Y}|=3$, and generate
$P_{X,S,Y}$ in a manner similar to Appendix D.1.1, with the same random seed.
Then, we generate $N=50\,000$ training samples of $(X,S,Y)$ triples. In our
implementation, we set ${\bar{k}}=|\mathcal{S}|-1=2$ and $k=1$, and set
feature extractors
$\bar{f}\in\mathcal{F}_{\mathcal{X}}^{\,2},\bar{g}\in\mathcal{F}_{\mathcal{S}}^{\,2}$,
$f\in\mathcal{F}_{\mathcal{X}}$,
$g\in\mathcal{F}_{\mathcal{S}\times\mathcal{Y}}$ as corresponding one-layer
linear networks with one-hot encoded inputs. In particular, we convert each
$(s,y)$ to one unique one-hot vector of in
$\mathbb{R}^{|\mathcal{S}|\cdot|\mathcal{Y}|}$ as the input to the network
$g$. Then, we train these feature extractors on the training set with a
minibatch size of 256 for 100 epochs, to maximize the nested H-score
configured by $\mathcal{C}_{\mathsf{MC}}$ [cf. (56)].

For comparison, we train a multihead network shown in Figure 10 on the same
dataset, with the same minibatch size and epochs. The feature $f$ is again
implemented by a one-layer linear network. In particular, we maximize the log-
likelihood function (67) to learn the corresponding feature and weights. Then,
we convert the weights to $g\in\mathcal{F}_{\mathcal{S}\times\mathcal{Y}}$ via
the correspondence [cf. (66)] $g(s,y)=G_{s}(1,y)$. The comparison between
features learned from two approaches are shown in Figure 19. For the sake of
presentation, we have normalized the features before plotting, such that $f$
and each $g(s,\cdot)$ are zero-mean, and unit variance with respect to
$\mathrm{Unif}(\mathcal{X})$ and $\mathrm{Unif}(\mathcal{Y})$, respectively.

###  D.4 Multimodal Learning With Missing Modalities

####  D.4.1 Implementation of Section 7.4.1

We first generate $N=50\,000$ triples of $(X_{1},X_{2},Y)$ for training. In
implementing the algorithm, we ${\bar{k}}=k=1$. To represent
$\bar{f}\in\mathcal{F}_{\mathcal{X}}$, we set each $\bar{f}^{(i)}$,
$i\in\\{1,2\\}$ as an MLP with dimensions $\text{input}=1$ – 32 – 32 –
$1=\text{output}$. To represent $f$, we use an MLP with dimensions
$\text{input}=2$ – 32 – 32 – $1=\text{output}$ with the input set to
$X_{1}\mathbin{+\mkern-10.0mu+}X_{2}=(X_{1},X_{2})^{\mathrm{T}}$. Since $Y$ is
discrete, we use one-layer linear network as $\bar{g}$ and $g$, where the
inputs are one-hot encoded $Y$. Therefore, both $\bar{g}$ and $g$ are with one
linear layer, taking $|\mathcal{Y}|=2$ dimensional input and outputting one-
dimensional feature.

We then train these feature extractors on the training set with a minibatch
size of 256 for 100 epochs, to maximize the nested H-score configured by
$\mathcal{C}_{\sf BI}$ [cf. (75)]. The learned features are then shown in
Figure 20.

For the prediction problem with missing modality, we use two MLPs to represent
conditional expectation operators
$\phi_{1}\triangleq\tau_{{1}}(\bar{f}^{(2)})$ and
$\phi_{2}\triangleq\tau_{{2}}(\bar{f}^{(1)})$, each with dimensions
$\text{input}=1$ – 32 – 32 – $1=\text{output}$. These two networks are
optimized by minimizing the corresponding mean square error. For example, we
train $\phi_{1}$ by minimizing
${\mathbb{E}}\left[\left(\bar{f}^{(2)}(X_{2})-\phi_{1}(X_{1})\right)^{2}\right]$
over all parameters in $\phi_{1}$ network. The training of $\phi_{1},\phi_{2}$
is with minibatch size of 256, and 100 epochs.

####  D.4.2 Implementation of Section 7.4.2

We generate $N=50\,000$ triples of $(X_{1},X_{2},Y)$ from (103) and (108), and
adopt the decomposition (86) on each triple. This gives three pairwise
datasets with samples of $(X_{1},X_{2})$, $(X_{1},Y)$, $(X_{2},Y)$, where each
dataset has $N$ sample pairs. Then, we adopt the same setting of networks to
represent one-dimensional $\bar{f}^{(1)},\bar{f}^{(2)}$ and $g$.

We then train these feature extractors on the three datasets for 100 epochs
with a minibatch size of 256, to maximize $\mathscr{H}(\bar{f},\bar{g})$.
Here, we compute $\mathscr{H}(\bar{f},\bar{g})$ based on the minibatches from
the three pairwise datasets according to (90).

To solve the prediction problem with missing modality, we use the same network
architectures and training settings to learn $\tau_{{1}}(\bar{f}^{(2)})$ and
$\tau_{{2}}(\bar{f}^{(1)})$, as introduced in Appendix D.4.1.

## References

  * Amari (2001) S-I Amari.  Information geometry on hierarchy of probability distributions.  _IEEE transactions on information theory_ , 47(5):1701–1711, 2001.
  * Amari and Nagaoka (2000) Shun-ichi Amari and Hiroshi Nagaoka.  _Methods of information geometry_ , volume 191.  American Mathematical Soc., 2000.
  * Belghazi et al. (2018) Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm.  Mutual information neural estimation.  In _International conference on machine learning_ , pages 531–540. PMLR, 2018.
  * Breiman and Friedman (1985) Leo Breiman and Jerome H Friedman.  Estimating optimal transformations for multiple regression and correlation.  _Journal of the American statistical Association_ , 80(391):580–598, 1985.
  * Buja (1990) Andreas Buja.  Remarks on functional canonical variates, alternating least squares methods and ace.  _The Annals of Statistics_ , pages 1032–1069, 1990.
  * Caruana (1993) R Caruana.  Multitask learning: A knowledge-based source of inductive bias1.  In _Proceedings of the Tenth International Conference on Machine Learning_ , pages 41–48. Citeseer, 1993.
  * Chan et al. (2022) Kwan Ho Ryan Chan, Yaodong Yu, Chong You, Haozhi Qi, John Wright, and Yi Ma.  Redunet: A white-box deep network from the principle of maximizing rate reduction.  _The Journal of Machine Learning Research_ , 23(1):4907–5009, 2022.
  * Cortes and Vapnik (1995) Corinna Cortes and Vladimir Vapnik.  Support-vector networks.  _Machine learning_ , 20(3):273–297, 1995.
  * Cover and Thomas (2006) Thomas M Cover and Joy A Thomas.  Elements of information theory (wiley series in telecommunications and signal processing), 2006.
  * du Pin Calmon et al. (2017) Flavio du Pin Calmon, Ali Makhdoumi, Muriel Médard, Mayank Varia, Mark Christiansen, and Ken R Duffy.  Principal inertia components and applications.  _IEEE Transactions on Information Theory_ , 63(8):5011–5038, 2017.
  * Eckart and Young (1936) Carl Eckart and Gale Young.  The approximation of one matrix by another of lower rank.  _Psychometrika_ , 1(3):211–218, 1936.
  * Ekert and Knight (1995) Artur Ekert and Peter L Knight.  Entangled quantum systems and the schmidt decomposition.  _American Journal of Physics_ , 63(5):415–423, 1995.
  * Fisher (1936) Ronald A Fisher.  The use of multiple measurements in taxonomic problems.  _Annals of eugenics_ , 7(2):179–188, 1936.
  * Gebelein (1941) Hans Gebelein.  Das statistische problem der korrelation als variations-und eigenwertproblem und sein zusammenhang mit der ausgleichsrechnung.  _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik_ , 21(6):364–379, 1941.
  * Greenacre (2017) Michael Greenacre.  _Correspondence analysis in practice_.  CRC press, 2017.
  * HaoChen et al. (2021) Jeff Z HaoChen, Colin Wei, Adrien Gaidon, and Tengyu Ma.  Provable guarantees for self-supervised deep learning with spectral contrastive loss.  _Advances in Neural Information Processing Systems_ , 34:5000–5011, 2021.
  * Harris et al. (2020) Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Hastie et al. (2009) Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer, 2009.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Hirschfeld (1935) Hermann O Hirschfeld.  A connection between correlation and contingency.  In _Proceedings of the Cambridge Philosophical Society_ , volume 31, pages 520–524, 1935.
  * Hofmann et al. (2008) Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola.  Kernel methods in machine learning.  2008\.
  * Hotelling (1936) Harold Hotelling.  Relations between two sets of variates.  _Biometrika_ , 28(3/4):321–377, 1936.
  * Hsu et al. (2021) Hsiang Hsu, Salman Salamatian, and Flavio P Calmon.  Generalizing correspondence analysis for applications in machine learning.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 44(12):9347–9362, 2021.
  * Huang et al. (2019) Shao-Lun Huang, Anuran Makur, Gregory W Wornell, and Lizhong Zheng.  On universal features for high-dimensional learning and inference.  _arXiv preprint arXiv:1911.09105_ , 2019.
  * Hunter (2007) John D Hunter.  Matplotlib: A 2D graphics environment.  _Computing in science & engineering_, 9(03):90–95, 2007.
  * Jaynes (1957a) Edwin T Jaynes.  Information theory and statistical mechanics.  _Physical review_ , 106(4):620, 1957a.
  * Jaynes (1957b) Edwin T Jaynes.  Information theory and statistical mechanics. ii.  _Physical review_ , 108(2):171, 1957b.
  * Kim (2014) Yoon Kim.  Convolutional neural networks for sentence classification.  In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 1746–1751, Doha, Qatar, October 2014. Association for Computational Linguistics.  doi: 10.3115/v1/D14-1181.  URL https://aclanthology.org/D14-1181.
  * Kingma and Ba (2015) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Yoshua Bengio and Yann LeCun, editors, _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.  URL http://arxiv.org/abs/1412.6980.
  * Krizhevsky et al. (2017) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  _Communications of the ACM_ , 60(6):84–90, 2017\.
  * LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Ngiam et al. (2011) Jiquan Ngiam, Aditya Khosla, Mingyu Kim, Juhan Nam, Honglak Lee, and Andrew Y Ng.  Multimodal deep learning.  In _ICML_ , 2011.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32, 2019.
  * Pearson (1901) Karl Pearson.  Liii. on lines and planes of closest fit to systems of points in space.  _The London, Edinburgh, and Dublin philosophical magazine and journal of science_ , 2(11):559–572, 1901.
  * Rényi (1959) Alfréd Rényi.  On measures of dependence.  _Acta Mathematica Academiae Scientiarum Hungarica_ , 10(3-4):441–451, 1959.
  * Ruder (2017) Sebastian Ruder.  An overview of multi-task learning in deep neural networks.  _arXiv preprint arXiv:1706.05098_ , 2017.
  * Sason and Verdú (2016) Igal Sason and Sergio Verdú.  $f$-divergence inequalities.  _IEEE Transactions on Information Theory_ , 62(11):5973–6006, 2016.  doi: 10.1109/TIT.2016.2603151.
  * Schmidt (1907) Erhard Schmidt.  Zur theorie der linearen und nichtlinearen integralgleichungen. i. teil: Entwicklung willkürlicher funktionen nach systemen vorgeschriebener.  _Mathematische Annalen_ , 63:433–476, 1907.
  * Shannon (1948) Claude E Shannon.  A mathematical theory of communication.  _The Bell system technical journal_ , 27(3):379–423, 1948.
  * Tishby and Zaslavsky (2015) Naftali Tishby and Noga Zaslavsky.  Deep learning and the information bottleneck principle.  In _2015 ieee information theory workshop (itw)_ , pages 1–5. IEEE, 2015.
  * Tishby et al. (2000) Naftali Tishby, Fernando C Pereira, and William Bialek.  The information bottleneck method.  _arXiv preprint physics/0004057_ , 2000.
  * Van Rossum and Drake (2009) Guido Van Rossum and Fred L. Drake.  _Python 3 Reference Manual_.  CreateSpace, Scotts Valley, CA, 2009.  ISBN 1441412697.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Wang et al. (2019) Lichen Wang, Jiaxiang Wu, Shao-Lun Huang, Lizhong Zheng, Xiangxiang Xu, Lin Zhang, and Junzhou Huang.  An efficient approach to informative feature extraction from multimodal data.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 5281–5288, 2019.
  * Weidmann (2012) Joachim Weidmann.  _Linear operators in Hilbert spaces_ , volume 68.  Springer Science & Business Media, 2012.
  * Xu and Huang (2020) Xiangxiang Xu and Shao-Lun Huang.  Maximal correlation regression.  _IEEE Access_ , 8:26591–26601, 2020.
  * Xu and Huang (2021) Xiangxiang Xu and Shao-Lun Huang.  An information theoretic framework for distributed learning algorithms.  In _2021 IEEE International Symposium on Information Theory (ISIT)_ , pages 314–319. IEEE, 2021.
  * Xu and Zheng (2022) Xiangxiang Xu and Lizhong Zheng.  Multivariate feature extraction.  In _2022 58th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 1–8. IEEE, 2022.
  * Xu and Zheng (2023) Xiangxiang Xu and Lizhong Zheng.  Kernel subspace and feature extraction.  In _2023 IEEE International Symposium on Information Theory (ISIT) (ISIT’2023)_ , pages 1032–1037, Taipei, Taiwan, June 2023.
  * Xu et al. (2022) Xiangxiang Xu, Shao-Lun Huang, Lizhong Zheng, and Gregory W Wornell.  An information theoretic interpretation to deep neural networks.  _Entropy_ , 24(1):135, 2022.
  * Zhang and Wallace (2017) Ye Zhang and Byron C Wallace.  A sensitivity analysis of (and practitioners’ guide to) convolutional neural networks for sentence classification.  In _Proceedings of the Eighth International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 253–263, 2017.