{'fraction_non_alphanumeric': 0.06071546352385746, 'fraction_numerical': 0.023484283061114677, 'mean_word_length': 3.8529312668463613, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 12, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deep Sufficient Representation Learning via Mutual Information

We propose a mutual information-based sufficient representation learning (MSRL) approach, which uses the variational formulation of the mutual information and leverages the approximation power of deep neural networks. MSRL learns a sufficient representation with the maximum mutual information with the response and a user-selected distribution. It can easily handle multi-dimensional continuous or categorical response variables. MSRL is shown to be consistent in the sense that the conditional probability density function of the response variable given the learned representation converges to the conditional probability density function of the response variable given the predictor. Non-asymptotic error bounds for MSRL are also established under suitable conditions. To establish the error bounds, we derive a generalized Dudley's inequality for an order-two U-process indexed by deep neural networks, which may be of independent interest. We discuss how to determine the intrinsic dimension of the underlying data distribution. Moreover, we evaluate the performance of MSRL via extensive numerical experiments and real data analysis and demonstrate that MSRL outperforms some existing nonlinear sufficient dimension reduction methods.

1. Introduction. Learning a good representation of a predictor X for predicting a response Y is a fundamental problem in statistics and machine learning. An effective data representation is crucial to the success for downstream modeling and analysis tasks. In this paper, we propose a nonparametric method for learning a sufficient representation of X. We formulate an objective function based on a variational form of the mutual information and use deep neural networks to approximate the representation function. We shall refer to our proposed method as Mutual-information-based Sufficient Representation Learning (MSRL).
MSRL is inspired by the idea underlying sufficient dimension reduction (SDR, Li, 1991;Cook, 1998). The goal of SDR is to learn a linear representation of X such that the response Y is conditionally independent of X given this representation, that is, the representation is sufficient. SDR is a semiparametric method in the sense that only the representation to be estimated is assumed to take a linear form, but no parametric form is assumed for the conditional distribution of Y given X. Li (1991) introduced the sliced inverse regression (SIR) for estimating the sufficient direction under the assumption that the predictor satisfies a linear expectation property. Subsequently, various related methods have been proposed under similar assumptions, including sliced average variance estimation (SAVE, Cook and Weisberg, 1991), principal Hessian directions (Li, 1992), minimum average variance estimation (Xia et al., 2002), and cumulative slicing estimation (Zhu, Zhu and Feng, 2010), among others. We refer to Cook (2007Cook ( , 2018 and Ma and Zhu (2013) and the references therein for thorough reviews of the SDR methods.
In recent years, several studies have developed nonlinear SDR methods using the kernel tricks. Lee, Li and Chiaromonte (2013) proposed the generalized sliced inverse regression (GSIR) and generalized sliced average variance (GSAVE) methods based on the reproducing kernel Hilbert space (RKHS) theory. Other kernel-based nonlinear SDR methods have been developed by Li, Artemiou and Li (2011), Li and Song (2017), and Shin et al. (2017), among others. For these kernel-based methods, an n × n Gram matrix is involved in the estimation procedure, where n is the sample size. Hence, these methods are computationally expensive when n is relatively large.
MSRL uses an objective function based on the mutual information for learning a sufficient representation R(X) so that Y and X are conditionally independent given R(X). It takes advantage of the strong approximation power of deep neural networks and model the representation R(X) using deep neural networks to extract nonlinear features in a nonparametric fashion. To facilitate the implementation of the proposed method, we consider a dual form of the mutual information criterion. Then the estimation of R(X) can be done using the adversarial training procedure (Goodfellow et al., 2014). MSRL can easily deal with the case when Y is a multi-dimensional response or a categorical variable. MSRL is easy-to-implement and computationally efficient, even for very large n.
The information bottleneck (IB) method (Tishby, Pereira and Bialek, 1999;Tishby and Zaslavsky, 2015) is also based on the mutual information. The goal of IB is to find a maximally compressed mapping of the predictor that can preserve as much information on the response as possible. To achieve this, IB maximizes the mutual information between Y and a parametrically modeled representation R(X), under a constraint that the mutual information between X and R(X) is smaller than a pre-specified nonnegative constant. Different from the IB method and motivated by a result from the optimal transport theory, we introduce a new constraint to enforce the representation to follow a user-specified distribution. Making use of the variational formulation of the f -divergence, we formulate MSRL as a minimax optimization problem and compute the estimated representation using the adversarial learning techniques (Goodfellow et al., 2014;Nowozin, Cseke and Tomioka, 2016).
Our contributions are as follows.
• We propose MSRL, a new nonparametric sufficient representation learning method based on f -mutual information and deep neural networks. We formulate the problem as that of maximizing the mutual information between the response and the representation of the predictor. Computationally, we use a dual form of the mutual information and parameterize the representation using deep neural networks. MSRL can directly handle multi-dimensional continuous response and can be easily modified to deal with categorical response variables. • We show that the learned representation is consistent in the sense that the conditional probability density function (p.d.f) of Y given the learned representation converges to the conditional p.d.f of Y given X. We also study the convergence properties by establishing the non-asymptotic error bounds under suitable conditions. To establish the non-asymptotic error bounds, we derive a generalized Dudley's inequality for an order-two U-process indexed by deep neural network functions using Rademacher chaos. This inequality may be of independent interest in other problems where U-processes arise. • We evaluate the performance of MSRL via extensive numerical experiments and real data analysis and demonstrate that MSRL outperforms some existing nonlinear sufficient dimension reduction methods. We also discuss the extension to accommodate the cases involving categorical variables and provide a feasible algorithm to determine the intrinsic dimension of the sufficient representation.
The rest of the paper is organized as follows. We first give some preliminary informationtheoretic background and describe our proposed MSRL method in Section ??. In Section 3, we establish the consistency property of MSRL and study its non-asymptotic error bounds under reasonable conditions. We provide a high-level description of the proofs for the theoretical results in Section 4 and briefly discuss the extension to handle categorical responses or predictor in Section 5. We describe a training algorithm and a way to determine the intrinsic dimension of the representation in Section 6. In Section 7 we conduct numerical experiments to evaluate the performance of MSRL. Concluding remarks are given in Section 8.
2. f -mutual information and sufficient representations. To learn a sufficient representation R, we use mutual information to characterize conditional independence. Let f : R + → R be a differentiable and convex function with f (1) = 0. Let P 1 and P 2 be two probability measures supported on a set Z and with p.d.f's p 1 (·) and p 2 (·), respectively. The f -divergence (Ali and Silvey, 1966) between P 1 and P 2 is defined as It can be easily verified that D f (P 1 || P 2 ) ≥ 0 and the equality holds if p 1 ≡ p 2 . Commonlyused f -functions include f (t) = t log t, f (t) = −(t + 1) log(t + 1/2) + t log t, and f (t) = (t − 1) 2 , corresponding to the well-known Kullback-Leibler (KL), Jensen-Shanon (JS), and χ 2 divergences. Throughout the paper, we let (X, Y ) be a pair of d X -vector of covariate and d Y -vector of response, where X ∈ X = [0, 1] dX and Y ∈ Y = [0, 1] dY . The probability measures induced by (X, Y ), X and Y are respectively denoted by P XY , P X and P Y , and their corresponding p.d.f's are p XY , p X and p Y . In addition, P X P Y denotes the product measure. Using these notation, the f -mutual information between X and Y is defined as and I f (X; Y ) = 0 if and only if X and Y are independent, when f (t) is strictly convex at t = 1 (Liese and Vajda, 2006;Esposito, Gastpar and Issa, 2020). Thus, the f -mutual information I f (X; Y ) characterises the independence between X and Y . We now state a fundamental data-processing inequality for f -mutual information in Lemma 2.1 below, which serves as the key theoretical basis of the proposed method.
LEMMA 2.1. For any measurable function g : X → R d , we have If all the relevant p.d.f's and conditional p.d.f's are positive, and f is strictly convex, the equality in (2.1) holds if and only if (iff) X and Y are conditionally independent given g(X), that is, The data-processing inequality in (2.1) can be also found in Polyanskiy and Wu (2017) and Sason (2019). We provide a proof of Lemma 2.1 in Section A in the supplementary material. In what follows, if X and Y are independent given a representation R(X), i.e., Y ⊥ ⊥ X|R(X), the representation R is said to be sufficient. By Lemma 2.1, we can characterize a sufficient representation V : X → R dV with d V ∈ N + and d V ≤ d X as a minimizer For a sufficient representation V and any injective transformation T , T • V is also a sufficient representation by the basic property of conditional probability. Based on a result from the optimal transport theory given in Lemma C.1, under some regularity conditions, there exists an injective map T * : R dV → [0, 1] dV depending on V such that which implies that R V (X) is a sufficient representation. We refer to it as a d V -dimensional uniformized sufficient representation (USR).
In addition, define We call d 0 the intrinsic dimension and denote a d 0 -dimensional USR by R 0 . Namely, Our goal is to learn such a d 0 -dimensional USR R 0 . A sufficient and necessary condition to ensure the existence of such a R 0 is that there exists a sufficient representation R with dim(R) = d 0 having the second order moment and absolute continuous distribution. This condition is mild as the sufficient representation condition is also required for all the SDR methods. Moment condition and absolute continuity condition are also found in the SDR literature, such as the conditions for Theorem 4 in Li, Artemiou and Li (2011) and Assumptions 1 and 2 in Shin et al. (2017). In what follows, we assume that a d 0 -dimensional USR R 0 exists.
Remark 1. The reference distribution γ U is user-specified. We choose the uniform distribution for simplicity. Other distributions can also be considered, such as the multivariate Gaussian distribution and the uniform distribution on a sphere.
2.1. Population objective function based on mutual informaiton. In the rest of this paper, we focus on the case when f (t) = t log t, that is, D f is the KL-divergence and I f is the mutual information. The KL-divergence and mutual information, denoted by D KL and I KL respectively, have wide applications in statistics (Nguyen, Wainwright and Jordan, 2010;Berrett and Samworth, 2019). To estimate R 0 in (2.3), by the definition of R 0 , we note that where P T denotes the probability measure induced by T (X) for any transformation T of X. These properties plus the nonnegativity of KL divergence motivate us to consider minimizing the population-level objective function: where dim(R) = d 0 and λ is a nonnegative constant. Heuristically, minimizing −I KL (Y ; R(X)) enforces the conditional independence of X and Y given R(X), and the regularization term λD KL (P R || γ U ) with λ > 0 encourages R(X) to be distributed as the reference distribution Uniform[0, 1] d0 . The following theorem shows that the target R 0 is a minimizer of L(R; λ) over all measurable functions R satisfying dim(R) = d 0 .
A proof of Theorem 2.1 is given in Section B in the supplementary material. To facilitate the estimation of R 0 based on L(R; λ), we use a variational formulation of f -divergence (Lemma 1 of Nguyen, Wainwright and Jordan, 2010).
LEMMA 2.2. Suppose that two distributions P 1 and P 2 are supported on the same set Z ⊆ R d . For any differentiable and convex function f : R + → R with f (1) = 0, the fdivergence between P 1 and P 2 . Here, f denotes the first derivative of f .
When f (t) = t log t, its Fenchel conjugate f * (t) = e t−1 . It then follows from Lemma 2.2 that The sample version of the variational expression in (2.6) can be easily obtained by replacing the expectations by their empirical analogues.
By (2.6), we can write (2.4) as Hence, using the expression in (2.7), we can obtain the sample version of L(R; λ) by approximating the function classes of R, D, and Q, and using the empirical analogues for the expectations. We introduce the neural network function classes used to approximate the function classes of R, D, and Q in the next subsection.
2.2. Deep neural networks. The class of feedforward neural networks (FNNs) F consists of functions F θ : R din → R dout with input dimensionality dim in (F) = d in , output dimensionality dim out (F) = d out , weight and bias parameters θ, depth L, width W, size S, number of neurons U , and F θ satisfying F θ ∞ ≤ B for some 0 ≤ B ≤ ∞, where · ∞ is the sup-norm on some specific domain. Specifically, the architecture of F θ can be expressed as: and bias vector b i ∈ R ki+1 , i = 0, 1, . . . , L, and σ is the componentwise linear rectified linear unit (ReLU) activation function. For i = 1, 2, . . . , L, k i is the number of neurons (width) in the ith hidden layer; then k 0 = d in and k L+1 = d out . For this network, we say that it has L hidden layers when considering the layers containing the activation function, or (L + 1) layers in total when considering the layers involving affine transformations. We write the parameters of the network as θ = (Θ 0 , Θ 1 , . . . , Θ L , b 1 , b 2 , . . . , b L ). The width parameter W = max{k i , i = 1, . . . , L} is the maximum width of hidden layers; the number of neurons U is defined as the number of neurons in f θ , i.e., U = L i=1 k i ; the size S is the total number of parameters in the network, i.e., S = L i=0 k i+1 × (k i + 1). The multi-layer structure implies that S ≤ W(d in + 1) + (W 2 + W)(L − 1) + (W + 1)d out .
Note that the number of hidden layers L may depend on n, d in , and W, U, S may depend on n, d in and d out , but we suppress the dependence for notational simplicity. We simply refer to all ReLU-activated FNNs using the ReLU activation function as ReLU networks. We use three FNNs R, D, and Q with parameters given below to approximate the functions R, D, and Q in (2.7), respectively: In what follows, R is called the representer network, D the MI discriminator network, and Q the push-forward network.
2.3. Empirical objective function. In practice, the distribution of (X, Y ) is usually unknown and only a random sample {(X i , Y i ), i = 1, . . . , n} that are i.i.d. copies of (X, Y ) is available. Let U 1 , . . . , U n be i.i.d. observations sampled from Uniform[0, 1] d0 . The empirical objective function for the proposed MSRL estimator is an empirical version of L(R; λ), defined as Then the MSRL representation iŝ R λ n ∈ argmin R∈R L net n (R; λ). Below, we use an example to illustrate how MSRL differs from two existing nonlinear SDR methods: generalized sliced inverse regression (GSIR) and generalized sliced average variance estimation (GSAVE). We consider the model where β 1 , β 2 are two parameter vectors generated from the uniform distribution on the tendimensional sphere, X ∼ N (0, I 10 ), 1 , 2 ∼ N (0, 0.25) independently, and c = 5. Intuitively, a better representation should be able to better predict the sign of the response. We learn the two-dimensional representations, apply them to 1000 testing data points, and then plot these transformed data along with their response sign. We use two reference distributions for this experiment, γ U ∼ Uniform[0,1] or γ S ∼ sin(Z) with Z ∼ N (0, 1). The training details are given in Section D.1 in the supplementary material. and the results are displayed in Figures  1 and 2. From Figure 1, compared to GSIR and GSAVE, the learned MSRL features are more predictive for the response sign as the transformed data based on the learned MSRL features behave more similarly to the linearly separable data. This partially implies the sufficiency of the learned representation through MSRL. The results in Figure 2 show that the learned representation through MSRL has the desired reference distribution. 3. Theoretical results. In this section, we establish consistency results and nonasymptotic error bounds for the proposed MSRL method. The intrinsic dimension d 0 is assumed known for simplicity. Our theoretical analysis focus on continuous type random variables X and Y . Extension to the categorical case is discussed in Section 5.
3.1. Consistency. For any d ∈ N + , we denote the class of continuous functions and the class of Lipschitz functions by For any transformation T of X, let p T and p Y |T denote the p.d.f of T (X) and the conditional p.d.f of Y given T (X). We need the following regularity conditions in establishing a consistency result.
For ease of notation, we define three functions W (·), L(·), and S(·) on N + that are needed to specify the required neural network structures: for any d ∈ N + , for a positive constant C 1 , where a for any a ∈ R denotes the smallest integer no less than a. We specify the network structures as follows. NS 1. The representer network R has parameters: depth . And for every R ∈ R, 0 ≤ R(x) ≤ 1 for any x ∈ [0, 1] dX , where the inequality holds in a componentwise sense. NS 2. The discriminator network D has parameters: depth NS 3. The discriminator network Q has parameters: depth THEOREM 3.1. Suppose Assumptions 1 & 2 hold and the network parameters of R, D, and Q satisfy (NS1), (NS2) and (NS3) respectively, then The proof of Theorem 3.1 is deferred to Section B in the supplementary material.
COROLLARY 1 (Consistency). Suppose Assumptions 1 & 2 hold and the network parameters of R, D, and Q satisfy (NS1), (NS2) and (NS3) respectively, then Corollary 1 shows that the conditional p.d.f of Y given the the learned representationR λ n converges to the conditional p.d.f of Y given X, which implies the asymptotic sufficiency ofR λ n . The slow convergence or consistency in Theorem 1 is mainly due to the fact that the class of continuous function is much larger than a bounded and smooth function class. The difficulties in the approximation of a very general continuous function class are responsible for the arbitrarily slow convergence rate.
3.2. Error bound. To derive error bounds with faster convergence rate, we impose some slightly stronger smoothness assumptions on the relevant function classes. In this subsection, we establish non-asymptotic error bound and faster convergence rates under the Hölder continuity assumptions. We first give the definition of a Hölder class below.
DEFINITION 1 (Hölder class). A Hölder class H β ([0, 1] d , B) with β = k + a, k ∈ N + , a ∈ (0, 1] and a finite constant B > 0, is a function class consisting of function f : In this subsection, other than Assumptions 1 & 2, we make the following Hölder continuity assumptions for the target and related density functions. We use the same smoothness parameter β for all function classes for simplicity and when different function classes have various smoothness parameters, we can use the smallest smoothness parameter in the analysis. In the following Assumptions 3 & 4, the smoothness parameter β > 0 and B 2 is a constant satisfying 1 ≤ B 2 ≤ ∞. ASSUMPTION 4. For every R ∈ R, p Y |R (y|r) ∈ H β ([0, 1] dY +d0 , B 2 ), and 1/B 2 ≤ inf y,r p Y |R (y|r). In addition, p R (r) ∈ H β ([0, 1] d0 , B 2 ) and 1/B 2 ≤ inf r p R (r).
A straightforward but tedious calculation yields that where c 1 (s) = 2 s−1 , c 2 (s) = c 1 (s)[c 2 (s − 1) + 1] for any s ∈ N + and c 2 (0) = 1, and a for any a ∈ R is the greatest integer smaller than a. Similar results also hold for log p Y |R (y|r), log p Y (y) and log p R (r).
To describe the network structures needed to achieve a faster convergence rate, we define the following three functionsW (·),L(·), andS(·) on N + : for any d ∈ N + , for a positive constant C 2 .
NS 4. Representer network R has parameters: depth . And for every R ∈ R, 0 ≤ R(x) ≤ 1 for any x ∈ [0, 1] dX , where the inequality holds in a componentwise sense.
NS 5. The discriminator network D has parameters: depth THEOREM 3.2 (Non-asymptotic risk bounds). Suppose Assumptions 3 & 4 hold and the network parameters of R, D, and Q satisfy (NS4), (NS5) and (NS6) respectively. Then, COROLLARY 2 (Non-asymptotic error bounds). Under the conditions of Theorem 3.2, we have Here C is a constant only depending on the constants B 2 and β in Assumptions 3 and 4.
The excess error bound (3.5) follows from the error decomposition inequality in Lemma 4.1 and where a b means a is smaller than b up to a constant factor independent of n. The error bounds (3.6) and (3.7) follow from the excess error bound and Pinsker's inequality. We defer the proofs of Theorem 2 and Corollary 2 to Section B in the supplementary material. An appealing feature in Theorem 2 and Corollary 2 is that the prefactors in the error bounds depend on the relevant dimensions polynomially.
Remark 2. The cubic support assumption X × Y = [0, 1] dX × [0, 1] dY is made for technical convenience. It can be relaxed to the unbounded cases. For the unbounded case, under the assumption that P ( Z ∞ ≥ c 1 log n) ≤ c 2 /n for some positive constants c 1 , c 2 and Z = (Y, X), the consistency result in Theorem 1 and error bounds in Theorem 2 can be derived, at a small price of an additional polynomial log n factor in the convergence rate.
Remark 3. By Corollary 2, we see that if λ > 2, C 1λ = Cλ 1/2 and C 2λ = C, where C is defined in Corollary 2. This implies that λ with a large value would slow the convergence rate for the asymptotic sufficiency ofR λ n in (3.6) and does not impact the pushing-forward convergence rate in (3.7). In contrast, if λ ≤ 2, C 1λ = √ 2C and C 2λ = C(2/λ) 1/2 . Thus λ with a small value would not impact the convergence rate for the asymptotic sufficiency but deteriorate the pushing-forward convergence rate. Based on these observations, we can use a constant value for λ. In particular, we use λ = 2 in our numerical experiments below, which leads to satisfactory performance of MSRL.

4.
A high-level description of the proofs. In this section, we provide a description of the main ideas and technical results needed for proving the error analysis results stated in Section 3. 4.1. Error decomposition. We first give two basic lemmas that provide a basis for the error analysis. The first lemma gives an error decomposition for bounding the excess risk and the second lemma bounds the L 1 distance between two density functions in terms of the excess risk.
For any measurable functions R, D and Q, let n and L Push n are defined in (2.8) and (2.9) respectively.
LEMMA 4.1 (Error Decomposition). LetR λ n ∈ argmin R∈R L net n (R; λ), and let R * ∈ argmin R∈R L(R; λ) be the best representation in the neural network class R. Then, for any fixed λ ≥ 0, The proof of Lemma 4.1 is given in Section B in the supplementary material. The error decomposition in Lemma 4.1 is a key step in the proof of Theorem 1, as we only need to bound each term on the right-hand side of the error decomposition inequality in Lemma 4.1.
The terms ∆ 1 , ∆ 2 and ∆ 5 are the approximation errors, which can be controlled using the approximation results with deep neural networks (Shen, Yang and Zhang, 2020;Jiao et al., 2021). The terms ∆ 3 and ∆ 4 are stochastic errors. Specifically, ∆ 3 involves empirical processes and can be bounded through some classic empirical process results such as Dudley's inequality; see Lemma C.3 in Section C in the supplementary material. However, the term ∆ 4 involves an order-two U-process, and the classic Dudley's inequality is no longer applicable. We will derive an inequality to bound this term below.
The second basic lemma bounds the L 1 distances of density functions in terms of the excess risk.
LEMMA 4.2. It holds that Lemma 4.2 can be easily verified as follows. By the definition of the objective function L given in (2.4), for any R ∈ R, where (4.2) follows from the chain rule for conditional mutual information (Theorem 2.5.2 in Cover and Thomas, 2006), i.e., I KL (Y ; (X, R(X))) = I KL (Y ; R(X)) + I KL (Y ; X|R(X)). Therefore, Lemma 4.2 follows directly from Pinsker's inequality (Tsybakov, 2008).
In view of (4.1), to prove (3.1), it suffices to show

4.2.
A maximal inequality for U-process of order 2. To bound ∆ 4 , we derive a generalized Dudley's inequality for an order-two U-process using Rademacher chaos. This inequality stated in the next lemma may be of independent interest. LEMMA 4.3 (Maximal inequality for U-Process of Order 2). Let H be a class of symmetric and real-valued functions on Z 2 and assume that sup whereh(z) = Eh(z, Z) for any z ∈ Z, then there exists a universal constant C such that and N 2n(n−1) (u, · ∞ , H) is the uniform covering number of H (Anthony and Bartlett, 1999). See also Section B in the supplementary material.
We note that some similar entropy bounds but with different metrics have been obtained in the existing literature, see Lemma 8 of Fukumizu et al. (2009) and Theorem 2 of Ying and Campbell (2010). These entropy bounds are similar to a special case for the inequality (4.3), which is where C is a universal constant. Obviously, the integral upper bound in (4.3) can be tighter than that in (4.4). A distinctive feature in (4.3) is the uniform covering number we use, which is naturally related to neural networks . By Theorem 12.2 of Anthony and Bartlett (1999), we have for n > Pdim(H), where Pdim(H) is the pseudo dimension of H ( Bartlett et al., 2019), see also Section B in the supplementary material. By Theorem 6 in Bartlett et al. (2019), if H is a ReLU network with depth L and size S, there exists a universal constant C 1 , such that Pdim(H) ≤ C 1 SL log S. Using these results and our inequality ( 5. Extensions to categorical response variables. We focus on continuous response and predictors in the previous sections and the theoretical analysis also relies on the continuousdistributed assumption. When there are categorical variables, we can easily modify our method to accommodate this case. We first assume only the response Y is categorical variable with finite categories. Namely, Y ∈ {1, 2, . . . , K} for some K ∈ N + . Note that where p k = P (Y = k) for k = 1, 2, . . . , K. In view of this, we consider the following empirical objective function e Q(Ui) , (5.1) and its minimizerR λ n ∈ argmin R∈R L net n,c (R; λ) is defined as a USR estimate, wherep k = n k /n, n k = n i=1 I{Y i = k} for k = 1, 2, . . . , K, and R, Q, and D k , k = 1, 2, . . . , K, are some neural network function classes. Letp = (p 1 ,p 2 , . . . ,p K ) and letỸ i ∈ R K be the one-hot vector with the kth entry being 1 when Y i = k. And D is a K-output neural network function class but each output does not share any weights or bias. Namely, D = (D 1 , D 2 , . . . , D K ) for any D ∈ D and D 1 , D 2 , . . . , D K do not share any parameter. Then, (5.1) can be recast as where e D(·) = (e D1(·) , e D2(·) , . . . , e DK (·) ) . This formulation is in terms of one single neural network and allows direct application of the MSRL training algorithm 1. Note that we can further extend the above approach to allow more general categorical cases. For instance, we assume Y = (Y D , Y C ) and X = (X D , X C ), where the subscripts D and C indicate the categorical and continuous components of a random vector, respectively. Moreover, it is assumed that Y D ∈ {1, 2, . . . , K} and X D ∈ {1, 2, . . . , J} for some J ∈ N + . Our goal is to learn a representation R such that Y ⊥ ⊥ X|(X D , R(X C )). Analogous to (5.1), we consider the following empirical objective function wherep jk = n jk /n, n jk = n i=1 I{Y D,i = k, X D,i = j}, j = 1, 2, . . . , J , k = 1, 2, . . . , K, and R, Q, and D jk for j = 1, 2, . . . , J , and k = 1, 2, . . . , K, are some neural network function classes.
Our theoretical analysis in section 3 can be applied to study the convergence properties of the MSRL for the aforementioned categorical cases without further difficulties. We omit the details.
6. Implementation. In this section, we present a MSRL training algorithm and discuss some practical issues. We rewrite R, D, Q as R θ , D φ , Q ψ to underline their weight and bias parameters θ, φ and ψ, respectively. The training procedure is summarized in Algorithm 1.
When the batch size m is relatively large, such as m = 256 or 512, to reduce the computational burden, we approximate the U-process term in Algorithm 1 using permuted samples. The early-stopping (ES) criterion is applied to avoid the over-fitting phenomena, and the distance covariance or distance correlation (Székely, Rizzo and Bakirov, 2007) is used as the early-stopping criterion.
6.1. Determination of the intrinsic dimension. The determination of the intrinsic dimension d 0 is crucial in real implementation. In practice, if taking different downstream tasks into account, different determination methods may be adopted. For example, when the predictive performance is of more concern, we may determine d 0 by the cross-validation in terms of APE. In sections 7.2-7.3, we determine d 0 as 5 and d 0 as 10 respectively by this approach.
In this section, we propose a dichotomy-based cross-validation method to determine the intrinsic dimension d 0 . For k = 1, 2, . . . , d X , let for k = 1, 2, . . . , d X −1, where the first inequality follows from the data-processing inequality (2.1) and the second inequality is by the definition in (6.1). Notably, when k ≥ d 0 , ). This is because ) for any k, by the definition of USR for R 0 , and Applying MSRL with intrinsic dimension input k, we obtain the minimizer Under some mild conditions, I KL (Y ; R (k) 0 (X)) can be consistently estimated by This motivates a dichotomy-based cross-validation selection method as follows. Given the data {(Y i , X i )} n i=1 and tolerance level η, we partition the sample into folds I t (t = 1, . . . , T ) and denote the complement of I t by I c t . For any nonnegative integer k ≤ d X , we obtain the mutual information estimateÎ k,n whereR λ k,n,t andD λ k,n,t are trained based on I c t as we have done in (6.4) and (6.5). Let UD = d X , LD = 1, MI =Î UD,n KL , and u = (UD + LD)/2 . If |Î u,n KL − MI|/ MI ≤ η, we update UD = u, MI =Î u,n KL and otherwise update LD = u + 1. We continue this procedure until UD = LD and the estimate of the intrinsic dimension is the ultimate UD.
In practice, instead of starting from d X , one may start from certain upper bound d U for d 0 . The proposed method significantly reduces the computational complexity since it only screens log 2 d U alternative values for d 0 . We provide a simple simulation example to illustrate the proposed method in Subsection 7.4.
Fix R θ , update D φ by ascending its stochastic gradient: and update Q ψ by ascending its stochastic gradient:
7. Experiments. We conduct several simulation studies to examine the performance of the proposed method. For comparison, we compute two commonly-used linear SDR methods: sliced inverse regression (SIR) and sliced average variance estimation (SAVE), and two nonlinear SDR methods: generalized sliced inverse regression and generalized sliced average variance estimation (GSIR and GSAVE), as the competitors. In Subsection 7.1, we first adopt the conventional practice in many SDR existing literature to assume the intrinsic dimension is known, see such as Zhu, Miao and Peng (2006) and Lee, Li and Chiaromonte (2013).
After learning a representationR on the training data {(Y j , X j )} n j=1 , we apply it to the testing data {(Y n+j , X n+j )} m j=1 and obtain the transformed data {(Y n+j ,R(X n+j ))} m j=1 . Based on {(Y n+j ,R(X n+j ))} m j=1 , we then calculate the empirical distance correlation (DC) of Y and R(X). To assess the predictive power of the learned representation, we fit a linear regression between Y and R(X) by using the transformed training data {(Y j ,R(X j ))} n j=1 and subsequently employ the fitted linear model LR to the testing data {(Y n+j ,R(X n+j ))} m j=1 . The average prediction error (APE) is calculated by In the simulation studies, we generate 6000 data points from each model below: 1 + X 2 2 ) 1/2 log((X 2 1 + X 2 2 ) 1/2 ) + ; III : Y = (X 1 + X 2 ) 2 /(1 + exp(X 1 )) + ; IV : Y = sin(π(X 1 + X 2 )/10) + X 2 1 + , where ⊥ ⊥ X = (X 1 , X 2 , . . . , X p ) , ∼ N (0, 0.25). Four distributional settings for the predictor X are considered: Notably, model complexity increases from model I to model IV, and from scenario (i) to scenario (iv). Full details of neural network architectures and hyperparameter choices for all the methods are given in Section D.3if the supplementary material.
The simulations are implemented in Python and the results for p = 10 are summarized in Table 1 from which we have several observations. First, our proposed MSRL outperforms all competitors in terms of larger distance correlation in almost all cases, implying that the features extracted by our method are more relevant to the outcome Y . Second, as expected, SIR and SAVE perform the best for model I which is relative simple model, but suffer from serious DC and APE loss for nonlinear models. GSIR and GSAVE show the best performance in terms of APE when the underlying distribution of the predictor X is simple. Among all these cases, our MSRL performs comparably well. Third, the performance of our MSRL is superior to other competitors in terms of APE when the underlying distributions of X are sophisticated. Additional simulations results for p = 30 are provided in Section D in the supplementary material. Overall, the proposed MSRL could be good supplement to other existing competitors in learning sufficient representation. 7.2. The superconductivity dataset. In this subsection, we apply MSDL to the superconductivity dataset, which is available at https://archive.ics.uci.edu/ml/datasets/Superconductivty+Data. In this dataset, there are 81 features of 21, 263 superconductors and the response variable is the critical temperature; see a relevant paper Hamidieh (2018) for the feature explanations. The goal of this analysis is to predict the critical temperature based on these features. The MSRL training procedure in Algorithm 1 is implemented with d 0 = 5, 10 and 20, respectively, and the reported results are based 5-fold validation. Additional implementation details are given in Section D.4 in the supplementary material. The distance correlation, average prediction errors, and their standard errors for all methods are reported in Table 2. It can be seen that MSRL outperforms other methods in terms of both empirical DC and prediction accuracy. In this example, SIR performs better than GSIR and GSAVE, which may be due to the curse of dimensionality for the kernel-based approaches and the relatively small effective sample size. The implementation details for GSIR and GSAVE are given in Section D.4 in the supplementary material.
In Figure 3, we plot the critical temperature against each component of the learned representation with d 0 = 5 for MSRL, GSIR and GSAVE. The plots for SIR and SAVE are given in Section D.5 in the supplementary material. We fit a quadratic model to check how the components impact the response. From Figure 3, we see that the critical temperature is nearly linear in three components of the MSRL estimate and among other SDR methods, only SIR and GSIR show that the critical temperature is nearly linear in their corresponding first dimension reduction direction. This may account for why MSRL has the high empirical DC and low prediction error when using a linear model. 7.3. The Pole-Telecommunication dataset. We next analyze the Pole-Telecommunication dataset (Weiss and Indurkhya, 1995), which is available at https://www.dcc.fc.up.pt/ ltorgo/Regression/DataSets.html. In this dataset, there are 15,000 instances with 48 predictors and a continuous distance-based response. Among the 48 predictors, 22 of them are non-informative, meaning that they have the same value for all subjects. So we delete all these non-informative features. Other implementation details are similar to those in section 7.2, see Section D.4 in the supplementary material, except that we take d 0 = 2, 5, 10, 20 in this example. The distance correlation, average prediction errors, and their standard errors for all methods are reported in Table 3. One can see that MSRL produces better results. The high empirical DC is supporting evidence of the high quality of the learned representation. Different from the previous example, in this example, GSIR and GSAVE outperform the linear SDR methods and GSIR performs the best among the four SDR competitors.
To visualize the predictive power of the learned representations, we classify the data into two classes: one has the response less than 50 and the rest constitutes another class. In Figure  4, we plot the transformed data based on each pair of the components of the learned representation with d 0 = 5. From Figure 4, we see that each pair of features from the MSRL estimate are informative for this classification task, and among other SDR methods, only the first dimension reduction directions of GSIR and GSAVE are predictive for the classification. The plots for SIR and SAVE are included in the Section D.5 in the supplementary material. This indicates that in this example the representation learned by MSRL has more predictive power. 7.4. Intrinsic dimension determination. To evaluate our proposed intrinsic dimension determination method, we consider a simple model where ⊥ ⊥ X = (X 1 , X 2 , . . . , X 10 ) , ∼ N (0, 1) and X ∼ N (0, I 10 ). In this example, the true d 0 = 2. The implementation details are given in Section D.6 in the supplementary material and we calculate the proportions of different estimated values for d 0 over 100 replications. The simulation results for two tolerance levels are provided in Table 4. It shows that the proposed intrinsic dimension determination method works well for this toy example since it did not estimate d 0 as 1, asd 0 = 1 would produce an insufficient representation.
8. Conclusion. In this paper, we have proposed MSRL, a mutual information based deep sufficient representation learning approach. Under reasonable conditions on the underlying target representation and the neural network classes, we have established the consistency and derived the non-asymptotic error bounds for MSRL. Extensive numerical experiments demonstrate that our method outperforms several existing SDR in the simulation models considered. We have also discussed the extensions to the case of categorical response variables and proposed a dichotomy-based cross-validation algorithm to determine the intrinsic dimension. Based on the results from our numerical experiments, MSRL outperforms the existing nonlinear SDR methods GSIR and GSAVE when the sample size is relatively large. With relatively small sample sizes such as n = 300 or 500, MSRL performs similarly to GSIR and GSAVE. Moreover, with big sample sizes, the kernel-based nonlinear SDR methods are computationally prohibitive due to an n × n Gram matrix involved in the estimation procedure. MSRL is still computationally feasible and performs well.
We have only considered the mutual information criterion corresponding to the Kullback-Leibler divergence in this work. In the future, it would be interesting to study the properties of using other mutual information criterions corresponding to divergencies such as Jensen-Shanon and χ 2 -divergences. It would also be interesting to extend the proposed method to the multi-dataset setting and apply it to problems arising from transfer learning and domain generalization.

Supplementary Material
In this supplementary material, we provide technical details and proofs of the lemmas and the theorems in the paper. We also include additional numerical results and details for the experiments.
APPENDIX A: PROOF OF LEMMA 2.1 PROOF OF LEMMA 2.1. By the fact that Y → X → g(X) is a Markov chain, Lemma 2.1 is a direct corollary of Lemma A.3 below.
LEMMA A.1. Suppose that distributions P XY and Q XY have p.d.f's p X (x)p Y |X (y|x) and q X (x)p Y |X (y|x), respectively, then This completes the proof of Lemma A.1.
LEMMA A.2. For two distributions P XY and Q XY , then If all relevant p.d.f's and conditional p.d.f's are positive, when f is strictly convex, the above equality holds if and only if where the inequality (S.3) follows from Jensen's inequality. This completes the proof of (S.1). When all relevant densities and conditional densities are positive, by the strict convexity of f , the inequality (S.3) becomes an equality if and only if This implies the result of (S.2).
If all relevant densities and conditional densities are positive, when f is strictly convex, the above equality holds if and only if PROOF OF LEMMA A.3. Note that (U, X, Y ) has the joint distribution density p U,X,Y satisfying p U,X,Y (u, x, y) = p U |X,Y (u|x, y)p X,Y (x, y) and p U |X,Y (u|x, y) = p U |X (u|x) by the Markov property. Let Q U,X,Y be the distribution with a joint density Then by the result of Lemma A.1, we have By the result of (S.1) in Lemma A.2, where the last equality follows from As p X|U,Y (x|u, y) = q X|U,Y (x|u, y)

By (S.2) in Lemma
where (S.5) follows from the construction of q U,X,Y (u, x, y) in (S.4), and (S.6) is due to the Markov property, which means p U |X,Y (u|x, y) = p U |X (u|x), we have under the assumed conditions that all relevant densities and conditional densities are positive, and f is strictly convex, This completes the proof of Lemma A.3.

APPENDIX B: PROOF OF MAIN RESULTS
In Appendices B and C, let Z 1 , Z 2 , . . . , Z n be i.i.d random vectors in a measurable space (Z, P), for example, , and let be an independent copy of . We use E Z [·] and E [·] to denote the expectation w.
For any R with dim(R) = d 0 , by the data processing inequality (2.1) and the nonnegativity of KL divergence, we know Hence, for any λ ≥ 0, we have

It follows that
This completes the proof of Theorem 2.1. PROOF OF LEMMA 4.1. For any R ∈ R, we have The proof is completed by choosing R = R * ∈ argmin R∈R L(R; λ).
PROOFS OF THEOREM 3.1 AND COROLLARY 1. By Lemma 4.1, we have where R * ∈ argmin R∈R L(R; λ). Using part (a) of Lemmas B.1-B.5, we obtain (S.2) E{L(R λ n ; λ) − L(R 0 ; λ)} → 0, as n → ∞. This completes the proof of Theorem 3.1. By Pinsker's inequality, we have where C is a constant only depending on B 2 and β, by bounding the terms on the right-hand side of the error decomposition inequality of Lemma 4.1.
where the last inequality holds due to the following facts: and QR λ n (r) = log pR λ n (r) on [0, 1] d0 .
Proof of Part (a): By Assumption 2, we have that QR λ n is a Lipschitz continuous function with a Lipschitz constant B 2 1 . By the network structure condition (NS3), and applying Theorem 2.1 in Shen, Yang and Zhang (2020) with N = n d 0 2(2+d 0 ) / log n and M = log n , there exists a ψ R ∈ Q such that By (S.5) and the arbitrariness of δ, we have where the last inequality follows from Q R ∞ ≤ B Q and ψ R ∞ ≤ B Q by Assumption 2 and (NS3).
Proof of Part (b): By Assumption 4 and the same argument as the one for (3.2), ).
By the network structure condition (NS6), and applying Theorem 3.3 in Jiao et al. (2021) with N = 1 and M = n d 0 2(2+d 0 ) , we get a ψ R ∈ Q, such that Due to the arbitrariness of δ, we have where the first inequality in (S.8) is due to Q R ∞ ≤ B Q and ψ R ∞ ≤ B Q by Assumption 4 and (NS6). Combining (S.5), (S.7), and (S.8), we obtain This completes the proof of Lemma B.1.
LEMMA B.2. (a) Suppose Assumption 2 holds and the network parameters of D satisfy (NS2), we have (b) Suppose Assumption 4 holds and the network parameters of D satisfy (NS5), we have PROOF OF LEMMA B.2. Similar to Lemma B.1, we get where the last inequality follows from the fact that, for any R ∈ R, and D R (y, r) = log pY |R (y|r) pY (y) . By Assumptions 2 and 4, respectively, a similar argument to (S.6) and (S.9) can be applied here and we obtain (b) Suppose the network parameters of R and D satisfy (NS4) and (NS5), then there exists a universal constant C, such that for n > Pdim(D ⊗ R) ∨ 2, where C 2,β (d, a) = (1 + a)( β + 1) 9/2 d β +2 .
where the last inequality follows from Lemma C.3. Using Theorem 12.2 in Anthony and Bartlett (1999) again, we have for n > Pdim(D ⊗ R), Similarly, for (S.12), we can get Again, using Theorem 12.2 in Anthony and Bartlett (1999), we have for n > Pdim(D ⊗ R), .
Putting together these integrals and covering number upper bounds, we can find a universal constant C such that for n > Pdim(D ⊗ R) ∨ 2, Pdim(D ⊗ R) log n n .
By part (a) of Lemma C.9, D ⊗ R can be implemented by a ReLU network . By Theorem 6 in Bartlett et al. (2019), for any ReLU network F with depth L and size S, there exists a universal constant C 1 , such that Based on these results, we shall proceed to prove part (a) and (b) separately.
Proof of Part (a): Given the settings in (NS1) and (NS2), there exists a universal constant Hence we can easily conclude that there exists a universal constant C such that Proof of Part (b): Using the settings in (NS4) and (NS5), we can find a universal constant C 3 , such that Applying this upper bound to (S.13), we conclude that there exists a universal constant C such that Then the rest of the proof for this lemma is similar to those for (S.14) and (S.15). This completes the proof of Lemma B.4.
By Assumption 2, we know Q R is a Lipschitz continuous and bounded function with a Lipschitz constant B 2 1 . Using (S.18), we have E|QR n (R n (X)) − QR n (R 0 (X))| ≤ B 2 1 E R n (X) − R 0 (X) → 0, as n → ∞. Invoking these results, the definition of R * , that is R * ∈ argmin R∈R L(R; λ), and (S.16), we have Proof of Part (b): For part (b), by (S.16), we further have for any R ∈ R, The first part of (S.19) is derived as below. When β < 1, for any (y i , r i ) ∈ [0, 1] dY × [0, 1] d0 , i = 1, 2, we have where c 1 lies between p Y |R0 (y 2 |r 2 ) and p Y |R0 (y 1 |r 1 ), and c 2 lies between p Y (y 2 ) and p Y (y 1 ). (S.21) is due to the mean value theorem, (S.22) is because p Y |R0 (y|r) ∈ H β ([0, 1] dY +d0 , B 2 ), and (S.23) follows from inf y,r p Y |R0 (y, r) ≥ 1/B 2 and inf y p Y ≥ 1/B 2 by Assumption 3. Similarly, by Assumption 4, we have Using these facts and that D R0 ∞ ≤ 2 log B 2 by Assumption 3, we can easily obtain the first part of (S.19). The proof for the second part of (S.19) follows the same line of argument as that for the first part of (S.19), with β in (S.23) and (S.24), B 2 in (S.23) and B 2 in (S.24) replaced by 1, , and the same line of argument as that for (S.17), when the network parameters of R satisfy (NS1), we find aR n ∈ R such that Due to the arbitrariness of δ, we have Invoking the definition of R * , that is R * ∈ argmin R∈R L(R; λ), and (S.20), we have This completes the proof of Lemma B.5.
APPENDIX C: SUPPORTING LEMMAS LEMMA C.1 (Example 11.2 in Villani (2009)). Let µ be a probability measure on R d and suppose it has finite second moment and is absolutely continuous with respect to the standard Gaussian measure γ d on R d . Then there exists a unique optimal transportation map T : R d → R d , which is also injective µ-almost everywhere, such that T # µ = γ d ≡ N (0, I d ), where T # µ(A) = µ(T −1 (A)) for any Borel measurable set A in R d .
Given Z 1 , Z 2 , . . . , Z n and a class H of measurable real-valued functions on Z, for any h ∈ H, we define , and the empirical Rademacher complexity for H w.r.t Z 1 , Z 2 , . . . , Z n , that iŝ LEMMA C.2 ( Symmetrization Inequality, Lemma 6.3.2 in Vershynin (2018)). Assume H is a function class of measurable functions, then  (2021)). Let N (u, H, · n ) denote the covering number for class H with covering radius u and w.r.t metric · n and assume sup h∈H h n ≤ b, then LEMMA C.4. Let A be a finite calss of n × n diagonal-free matrices, and 1 , 2 , . . . , n are i.i.d. Rademacher variables, then there exists a universal constant C such that where = ( 1 , 2 , . . . , n ) and A 2 F = i,j a 2 ij .
PROOF OF LEMMA C.4. When |A| = 1, the inequality is trivial. Hence, in what follows, we assume |A| > 1. Note that for any λ > 0, where the last inequality follows from Jensen's inequality. For any n × n, diagonal-free matrix A, by applying Lemma C.5 with f (x) = exp(λx), it holds where is an independent copy of . As i ψ2 = i ψ2 = 1/ log 2 for any i ∈ [n], by Lemma C.7, there exists a universal and positive constant C 1 such that Then by Lemma C.6 and C.8, we have for all λ satisfying |λ| ≤ c/ A F and some universal and positive constants c, C, C 2 .
LEMMA C.5 (Theorem 6.1.1 in Vershynin (2018)). Let A be an n × n, diagonal-free matrix. Let X = (X 1 , X 2 , . . . , X n ) be a random vector with independent mean zero coordinates X i . Then, for every convex function f : R → R, one has where X is an independent copy of X.
LEMMA C.6 (Lemma 6.2.3 in Vershynin (2018)). Consider independent mean-zero subgaussian random vectors X, X in R n with X ψ2 ≤ K and X ψ2 ≤ K, where X ψ2 = sup x∈S n−1 X, x ψ2 and W ψ2 = inf{t > 0 : E exp(W 2 /t 2 ) ≤ 2} for a scalar random variable W . Consider also independent random vectors g, g ∼ N (0, I n ). Let A be an n × n matrix. Then there exists a universal constant and positive C such that for any λ ∈ R.
LEMMA C.7 (Lemma 3.4.2 in Vershynin (2018)). Let X = (X 1 , X 2 , . . . , X n ) be a random vector with independent, mean zero, sub-gaussian coordinates X i . Then X is a subgaussian random vector, and there exists a universal constant and positive C such that LEMMA C.8 (Lemma 6.2.2 in Vershynin (2018)). Let X, X ∼ N (0, I n ) be independent and let A = (a ij ) be an n × n matrix. Then there exists two universal and positive constants c, C such that where the inequality follows from Jensen's inequality and the last equality is due to Fubini's Theorem and a direct corollary of the symmetrization inequality in Sherman (1994).
We prove the second inequality in (4.3) through the chaining technique. Let α 0 = b and for any t ∈ N + , let α t = 2 −t b. For each t, let C t be a α t -cover of H w.r.t. d U,n (·, ·) such that |C t | = N (α t , H, d U,n ). For each h ∈ H, there exists a functionĥ t ∈ C t such that d U,n (h,ĥ t ) ≤ α t . Letĥ 0 ≡ 0 and for any T ∈ N + , we have the following chaining expression for h: Hence for any T ∈ N + , we have where the last inequality is due to Now the second term in (S.4) is the summation of empirical Rademacher chaos w.r.t. the function classes C t − C t−1 = {h t − h t−1 : h t ∈ C t , h t−1 ∈ C t−1 }, t = 1, 2, . . . , T . For any measurable function h ∈ H, let Applying Lemma C.4 to C t − C t−1 , t = 1, 2, . . . , T , we have there exists a universal constant C such that S Deg n (H) ≤ α T + 6C T t=1 α t log(|C t |) n(n − 1) = α T + 12C T t=1 (α t − α t+1 ) log(N (α t , H, d U,n )) n(n − 1) ≤ α T + 12C α0 αT +1 log(N (u, H, d U,n ))du n(n − 1) ≤ α T + 24C α0 αT +1 log(N (u, H, d U,n ))du n .
For any 0 < δ < b, let T = min{K ∈ N + : 2 −K b ≤ δ}, then S Deg n (H) ≤ δ + 24C b δ/2 log(N (u, H, d U,n ))du n ≤ δ + 24C b δ/4 log(N 2n(n−1) (u, · ∞ , H))du n where the last inequality follows from (S.1). This completes the proof of Lemma 4.3. LEMMA C.9. (a) Let D ⊗ R = {h(z) = D(y, R(x)) : D : R dY +d0 → R ∈ D, R : R dX → R d0 ∈ R}. If D is a class of functions implemented by a ReLU network with width W D and depth L D , and R is a class of functions implemented by a ReLU network with width W R and depth L R , then each element of D ⊗ R can be implemented by a ReLU network with width max{W D , W R + 2d Y , 2(d Y + d 0 )} and depth L D + L R + 1.
(b) Let Q • R = {h(z) = Q(R(x)) : Q : R d0 → R dout ∈ Q, R : R dX → R d0 ∈ R}. If Q is a class of functions implemented by a ReLU network with width W Q and depth L Q , and R is a class of functions implemented by a ReLU network with width W R and depth L R , then each element of Q • R can be implemented by a ReLU network with width max{W Q , W R , 2d 0 } and depth L Q + L R + 1.
PROOF OF LEMMA C.9. The proof follows from the method in Subsection B.1.1 in Nakada and Imaizumi (2020), which shows how to construct a neural network from subneural networks through concatenation, parallelization, and realization of identity function.

APPENDIX D: ADDITIONAL NUMERICAL STUDIES AND EXPERIMENT DETAILS
This Appendix contains Table S.1 which includes some additional simulation results. D.1. Experimental Details. Python codes will be made publicly available on a GitHub repository, we still Here we provide the details for the numerical experimental details. In all the experiments, we set λ = 2 and use the Adam optimizer (Kingma and Ba, 2014) in Pytorch with learning rate lr and weight decay parameter wd specified in each experiment. The activation function for the neural networks is the LeakyReLU. The batch size in all experiments is 512.
D.2. Simulation for Model (2.10). The training data sample size is 4000 and no ES criterion is applied. The maximum number of epoches is 3000. When the reference distribution is γ U , R has 2 hidden layer with width (64, 64). When the reference distribution is γ S , R has 3 hidden layer with width (64, 64, 64). Both D and Q has 2 hidden layers with widths (32, 32). lr is 0.003 and wd is 0. For GSIR and GSAVE, we use some default tuning parameters. And the kernel estimates are based on 10000 new samples.
D.3. Simulation Study in Subsection 7.1. In this experiment, the total 6000 data points are first divided into 5000 training-validation data points and 1000 testing data points and then 5000 training-validation data points are divided into 4000 training data points and 1000 validation data points.
In MSRL training, ES criterion is applied and it is the empirical distance correlation (DC) calculated based on the validation data. The maximum number of epochs is 1000 and wd is 0.0001. The patience is 200, where patience is the number of epochs until termination if no improvement is made on the validation dataset. The learning procedure is replicated 10 times with different initial weights for the parameters of neural networks and we choose the one with the largest empirical DC on the validation dataset among the ten representation estimate as the ultimate representation estimate. R has 3 hidden layers with widths (32,16,8). Both D and Q has 2 hidden layers with widths ( For SIR and SAVE, the tuning parameter is the slice number sn. For each sn ∈ {5, 10, 15, 20, 25, 30}, we use the validation data to calculate the empirical DC between the response and the representation learned from training data. We use the slice number with the largest empirical DC as the ultimate choice for the tuning parameter sn. For GSIR and GSAVE, we apply the tuning parameter selection approach provided in Li (2018) to the validation data.
D.4. Real Datasets. For each data split, 20% of data is used as the testing data and then from the other 80% of the data, we randomly select 3000 data points as the validation data and the rest part is the training data.
For GSIR and GSAVE, limited by our computational capability, we only subsample 4000, 2000, 3000 data points from the above training, validation, and testing data as the new training, validation, and testing data, and conduct the estimation procedure we have stated in subsection D.3 In the superconductivty dataset, for our MSRL, R has 3 hidden layers with widths (128, 128, 128). Both D and Q has 2 hidden layers with widths (64, 64). The maximum number of epochs is 2000, lr = 0.0001, wd = 0.001 and the patience is 400. The number of different initial weights is still 10. For the goal of data visualization, we randomly sample 800 data points from the testing data for one fold to conduct plotting.
In the pole-telecommunication dataset, for MSRL, R has 2 hidden layers with widths (30, 25). Both D and Q has 2 hidden layers with widths (16, 8). The maximum number of epochs is 1000, lr = 0.001, wd = 0.0001 and the patience is 200. The number of different initial weights is still 10. For the goal of data visualization, we randomly sample 400 data points from the testing data for one fold to conduct plotting.
D.5. The extra plots for the superconductivity and pole-telecommunication datasets. In Figure S.1, we plot the critical temperature against each component of the learned representation with d 0 = 5 for SIR and SAVE. We fit a quadratic model to check how the components impact the response.