Deep Principal Support Vector Machines for Nonlinear Sufficient Dimension Reduction
Authors: Yinfeng Chen
Jin Liu
Rui Qiu
Venue: ICML 2025 poster (2025)
Keywords: principal support vector machines
neural networks
nonlinear sufficient dimension reduction
nonasymptotic convergence rate
Primary area: deep_learning->other_representation_learning
TL;DR: Sufficient Dimension Reduction; Deep Learning

Abstract: The normal vectors obtained from the support vector machine (SVM) method offer the potential to achieve sufficient dimension reduction in both classification and regression scenarios. Motivated by it, we in this paper introduce a unified framework for nonlinear sufficient dimension reduction based on classification ensemble. Kernel principal SVM, which leverages the reproducing kernel Hilbert space, can almost be regarded as a special case of this framework, and we generalize it by using a neural network function class for more flexible deep nonlinear reduction. We theoretically prove its unbiasedness with respect to the central $\sigma$-field and provide a nonasymptotic upper bound for the estimation error. Simulations and real data analysis demonstrate the considerable competitiveness of the proposed method, especially under heavy data contamination, large sample sizes, and complex inputs.

Deep Principal Support Vector Machines for Nonlinear Sufficient Dimension Reduction
YinFeng Chen"
Jin Liu?
Rui Qiu
Abstract
following conditional independence holds:
Y JL X
(1)
The normal vectors obtained from the support vector machine (SVM) method offer the potential to achieve sufficient dimension reduction in both classification and regression scenarios. Motivated by it, we in this paper introduce a unified frame- work for nonlinear sufficient dimension reduction based on classification ensemble. Kernel princi- pal SVM, which leverages the reproducing kernel Hilbert space, can almost be regarded as special case of this framework, and we generalize it by us- ing a neural network function class for more flex- ible deep nonlinear reduction. We theoretically prove its unbiasedness with respect to the central o- o-fieed and provide nonasymptotic upper bound for the estimation error. Simulations and real data analysis demonstrate the considerable competi- tiveness of the proposed method, especially under heavy data contamination, large sample sizes, and complex inputs.
where JL stands for statistical independence. The iden- tifiable parameter in (1) is the subspace spanned by the columns of 1, rather than itself, referred to as the suffi- cient dimension reduction subspace. Under mild conditions (see Cook (1998); Yin et al. (2008)), the intersection of all such spaces is also a sufficient dimension reduction sub- space and is called the central subspace, denoted by Fyxx- Representative methods of linear SDR include sliced inverse regression (SIR) (Li, 1991), sliced average variance estimate (SAVE) (Cook & Weisberg, 1991), directional regression (DR) (Li & Wang, 2007), minimum average variance esti- mation (MAVE) (Xia et al., 2002), and so on. We refer to Ma & Zhu (2013) for an extensive survey on linear SDR.
Sometimes due to complex regression relationships, linear combinations of predictor variables often fail to achieve effective SDR. To address this limitation, a nonlinear exten- sion of the linear SDR, as formulated in Cook (2007), is considered to seek any function f(X)
IR + Rd such that
Y Y X
(2)
1. Introduction
The kernel trick stands out as the most popular approach to extend linear SDR to nonlinear SDR, and related works in- clude Wu (2008), Yeh al. (2009), among others. Recently, Lee al. (2013) proposed rigorous SDR framework based on a-fileds, which integrates both linear and nonlinear SDR into a unified system. They define the SDR problem as seeking a sub o-field & of (XX such that
With the advent of the big data era, high-dimensional data have become increasingly common across diverse scientific domains. To this end, sufficient dimension reduction (SDR) techniques developed to extract liw-dimesiional rep- resentations of data while preserving much information as possible. These representations serve to facilitate both visualization and downstream tasks of data analysis. Let X E R2 be the predictor variables of dimension p, and Y E IR be the response variable. The classical linear SDR method aims to find a p x d(d < p) matrix 7) such that the
Y L 1XS.
(3)
Similar to the linear scenario, the intersection of all sub offields satisfying relation (3) still conforms to (3) under mild conditions. This minimal sub o-fiedd is then called the central o o-iiedd and denoted by Syjx. Definitions of (1) and (2) can be seamlessly integrated into this framework by considering & as ann X) and (f)()),, respectively. Compared to the linear case, the development of nonlin- ear SDR methods is relatively limited. As the nonlinear counterparts of SIR and SAVE, Lee et al. (2013) introduced
'Equal contribution 'School of Mathematics and Statistics, Xi' an Jiaotong University, Xi'an China "School of Data Science, The Chinese University of Hong Kong (Shenzhen). Shenzhen, China "Scoool of Mathematical Sciences, Peking University, Bei- jing. China. Correspondence to: Rui Qiu <rqiu@ @ pku.edu.cn>.
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Deep principal support vector machines
the generalized sliced inverse regression (GSIR) and gener- alized sliced average variance estimator (GSAVE). These methods may struggle with large samples and complex input structures.
is given. Encouragingly, the assumption can be satis- fied by Vapnik-Chervonent (VC) classes, with neural networks being a representative example. As result, we derive a specific convergence rate tailored to the neural network class, showcasing its optimality up to a polynomial order of log n when compared against the minimax rate of nonparametric regression.
The rapid development of machine learning has sparked interest in incorporating machine learning techniques into traditional dimension reduction, such as k-means inverse regression (Setodji & Cook, 2004), principal support vector machines (Li et al., 30, and dimension reduction forests (Loyal et al., 2022). In particular, deep neural networks have attracted the attention of many researchers due to their powerful feature extraction capabilities. A prominent repre- sentative of deep unsupervised dimension reduction is the auto-encoder models with bottleneck structures in represen- tation learning (Hinton & Salakhutdinov, 2006). Inspired by the mechanism of generative adversarial networks (GANs), Zheng et al. (2022) trained two deep neural networks ad- versarially and proposed mutual information-basee suf- ficient representation learning approach. Further, to avoid estimating the ratio of the joint density and the marginal densities, Huang et al. (2024) replaced mutual information with distance covariance and established the consistency of the proposal but not the rate of convergence. However, their strategy for limiting the distribution of sufficient representa- tion distinguished itself from the central o-filed framework used in our paper. Recently, Chen et al. (2024) proposed GMDDNet for nonlinear SDR, which combines the gener- alized martingale difference divergence criteria with neural networks. In addition, Kapla et al. (2022) introduced two- stage SDR estimation utilizing neural network instead of local linear smoothing typically employed in MAVE (Xia et al., 2002). However, it remains a linear method.
3.
The structure order determination method from classi- cal sufficient dimension reduction can be easily applied to our method, which is not the case for other neural network-based methods.
For any vector x, || and || represent the I, and supremum norm of æ, respectively. Let L2 (Px) be the space of square-integrable measurable functions with respect to the measure P'x of X. For any function
Notations:
I E L2 (Px), we define [[ = (FI 1/2 and Il = sup, lxx) as the L2(xx) norm and supre- mum norm of f. respectively. The e-covering number of a set with respect to the metric d is denoted by N A, d).
2. Preliminaries
To better elucidate our nonlinear SDR method and related theories, we first provide a detailed introduction to the frame- work of nonlinear SDR. Deep nonlinear SDR is the focus of this paper, for which we will briefly describe the structure of neural networks.
2.1. Nonlinear SDR
Let (2, denote a probability space, and let (2x, 3x), (Ny, and (Ixx, B xx) represent measurable spaces where Qx C R., Ry C R" and Ωxx C RR+0, By default, Zxy Ωx x Gy and ExY = 3x x By. Let X, Y, and (X, Y) be random vectors taking values in 2x Gy, and Ωxy, respectively, with distributions Px. Py, and PxY. The o-filed generated by X is oXX) = X-1 (xx). The conditional distribution of X given Y is denoted by (-I) : 2x x Qy + R. We consider the following nonlinear SDR problem
This paper will introduce a new nonlinear SDR method based on classification ensembles within the standard statis- tical framework established by (3). Our contribution lies in the following aspects:
1. Our nonlinear SDR methods can be applied a broad range of function classes. Among them, the neural network approximator shows significant promise. This deep method demonstrates flexibility in handling large samples and diverse input types including images and text. Importantly, it is robust for regression and clas- sification tasks, particularly suited for data with many outliers. Our deep method does not require additional special design of network structures such as bottleneck and generative adversarial structures.
Y L x(f(()
(4)
where fo(:) is a -ad-ud function (d < p). That is to say, at least functions of p-dimensional X are needed to fully capture the relationship between Y and X. Here, we assume the structural dimension d is known for simplic- ity. Obviously, fo is unidentifiable since we can substitute fo with any one-to-one function of it, without altering the conditional independence (4). In fact, this conditional inde- pendence relies solely on the o-iildd generated by (X), which is unique under very mild assumptions (Lee et al., 2013). Without loss of generality, we may always assume
2. We rigorously establish the unbiasedness of optimal solutions to the proposed objective function, indicat- ing that they belong to the central class. Under the assumption concerning the covering number of func- tion classes, the upper bound of the estimation error
2

---

Deep principal support vector machines
Assumption. The family of probability measures (+v)) : y [t) is dominated by a-finite measure. Then there exists a unique minimal sufficient o-field (or central o-field) Byjx such that
o (x) max 0} is the ReLU activation function ap- plied element-wise. Here, kn and kc+1 are the input and output dimensions, respectively. And kg i = 1, 1,2.... C stands for the number of neurons (width) in the ith hid- den layer. The depth of a network, denoted by L, refers to the number of hidden layers, while its width denotes the maximum width, A = max(ki, , kc). among the hidden layers. Moreover, the size of a network can be fur- ther characterized by two quantities: the total number of neurons U1 = -E kü and the total number of parameters S - Ei k++1 - (k; + 1). As the most widely used func- tion class, neural networks will serve as the main estimator for the nonlinear SDR functions in our methods.
Y IL X
Moreover, fo in (4) satisfies
F(fo(X)) 0, Var Var(Fo(X)) =Ta:
and
offo = By1x-
We introduce some examples to build intuition in Appendix A.I. Even though the central o-field is identifiable, it can be quite abstract to handle. Following Lee et al. (2013), let My be the central class corresponding the central a-field a (X)). This class is the collection of all square- integrable and lently, 1xx is a subset of L2(Px) containing all elements that can be expressed as a function of fo- To introduce the concept of unbiasedness for nonlinear SDR, we first re- view the unbiasedness in linear SDR. The central subspace FYix Span(1) can be viewed member of the param- eter space that encompasses all d-dimensional subspaces of IRP Let Pn be the empirical distribution of an i.i.d. sample {(X,,Y))...... Y.).... (X,, } from the distribution Pxv. T is a mapping from all possible distributions of (X, Y) to Rexd If (paan[((Pxx) xI ) S .YIx, we call T(P,) is an unbiased estimator of Syix- This equivalent to the inner product of X and any column vector of T(xx)) being measurable with respect to ann")) From this perspective, we can easily generalize the concept of unbiasedness to nonlinear SDR:
3. Motivation and Methodology
In this section, we introduce a new nonlinear SDR method based on the principal support vector machines (PSVM) proposed by Li et al. (2011). Then we give a reliable strat- egy for determining the structure dimension d in practical applications.
functions. Equiva-
3.1. Sample estimation for nonlinear SDR
When the response Y E {-1, 1] is binary, Li et al. (2011) observed that the normal vectors of the hyperplane learned by support vector machines (SVM) align with the directions forming the SDR subspace. Therefore these normal vectors can be naturally leveraged to construct the SDR subspace. Specifically, given training data 2n = {(X,, Y;)) Li et al. (2011) demonstrated the validity of the following objective function for linear SDR
min (V.D)RPXRR
Σ [1
Y. Ex)))) - -t]]"
Definition 2.1. A function f(X) from [2(P is unbiased for the central o o-field Byjx if f(X) E Myyxx. Such function is called an SDR function.
Var ( Tx),
(5)
where À is a hyperparameter, (.) and Var (-) denotes
the empirical mean and variance, and [x] = max (ar 0} is the hinge loss function. In particular, Li et al. (2011) proposed combining these normal vectors through principal component analysis form the basis of the SDR subspace. Zhou & Zhu (2016) then introduced the principal minimax support vector machines formulation
Unbiasedness is the most fundamental requirement for SDR estimators. It ensures that the information extracted by each SDR function solely originates from the regression relationship between X and Y.
2.2. Deep neural network
We now briefly introduce the feedforward neural networks with the rectified linear unit (ReLU) activation function. Specifically, the feedforward neural network can be repre- sented as
Al|ab11 + c² "T (v.1)ER" XR
min
2p + Var, Tx)
n
f(x) = WooooWl Wcooooo-oooooooooo
W1 o Wo(x),
where X, C are hyperparameters and et is special number. Additionally, Shin et al. (2017) proposed
where x E RP, Wx = W,F + be denotes an affine trans- formation with the weight matrix Wi E R xk, and inter- cept vector by E Rl+t for 0, 0, ; L. The function
A Σ
min (v.DERFXR XR
3 T
wx (Yt) (1 Y, T x, xxt)))
i=
3

---

Deep principal support vector machines
I. (OVA) In addition to the above {g-},=1 let 20 min min(Yi,..... (Yi,..., Yn) and Th = max max(YV,... 1,.,... Yn). For
where À is a hyperparameter and wx (Y)) - if
each (r, 6) satistying 1 rr<< - h, set
and π. otherwise.
The kernel trick or reproducing kernel Hilbert space (RKHS) provides an effective approach for extending linear SDR methods their nonlinear counterparts. However, finding suitable kernel for non-vector input remains challenging. For example, when dealing with tensor-valued input such as images, the conventional kernel trick often requires flat- tening tensors, which inevitably results in information loss. In contrast, neural networks show remarkable flexibility in handling diverse data types. Moreover, the computational burden of neural networks is typically lower than that of the kernel trick when handling large samples. Given that the kernel trick frequently involves inverting an n x n kernel matrix, the computational cost ( becomes formidable for large n.
8: I( < S..S0 q.) T(g-1 <,, S (5). -[g--1<
Over Fn x R. (Fres 4,) is the minimizer of L (J,) with Y, replaced by Y
Let Vi...., Va be the d leading eigenvectors of either of the matrices M with entries
Mij = Σ ((X))))()))
where f is the minimizer from different dividing
schemes (1 I i.j S h - for LVR and 1 S i,j S h(h 1)/2 for OVA). Then the ith (1 S S d) non-
Motivated by the superiority of neural networks over the kernel trick and the validity of normal vectors from SVM in linear SDR, we combine (5) with general functional class Fn (e.g., neural networks) and general loss p to achieve nonlinear SDR. Specifically, let Y be the binarization of Y, then we define an objective function
linear SDR function is
Eenf..
where Uyy is the jth entry of I-dimensional v, (1 = h for LVR and = h(h - 1)/2 for OVA).
tfr =X Var, Var,,f(X)) + Xt2
Li et al. (2011) empirically demonstrated that LVR works best when the response is continuous variable and OVA works best when the response is categorical.
12 p (-Y En(((X)) -t]
(6)
21
(f(x,) )(x,) 2t2
3.2. Estimation of structure dimension.
n(n-
1)
ISi<jSn
Determining the structural dimension d' is critical aspect of SDR. Here employ the ladle method proposed by Luo & Li (2016). This method utilizes not only the eigenvalue information of the matrix M but also its eigenvectors to determine the structure dimension. Let BN = (VI...., vk) and (Ai...., >x) respectively denote the x matrix com- prising the first eigenvectors of M and the vector compris- ing the first & eigenvalues of M. Utilizing bootstrapping we sample instances with replacement from dataset of size n. With these n samples, we re- re- re-implement B5, denot- ing its ith instance as B This bootstrapping procedure is repeated n times to generate multiple bootstrap samples. To measure the disparity between B, and its bootstrapped counterpart, we define the following function:
1 Σ - [(x.) 11x) - " ((x)) -
where E R. f is a member of Fa and p is any non- decreasing convex loss function. Particularly, p(x) =
max(1+ 11, x, 0} corresponds to the hinge loss in (5). Here we follow the slicing techniques proposed in Li et al. (2011) to discretize Y into the binary Y. One, which we call "left versus right (LVR)", repeatedly divides the data 2, into two groups according to set of cutting points for the response. The other, which we call "one versus another (OVA)", par- titions 2, into several slices and pairs up all possible slices. Below we summarize the the algorithm process:
1. (LVR) Let 9r = 1....,h h - 1. be - 1 cutting points. For example, they can be equally spaced sample percentiles of (Yi...., .., Y,).
= 0,
(2())) =
Σ Di=1 (1 [det(B
BID),
k == l.....l
For < r S h 1, set
Yr- [(Y > 9.) I(Y < gr). Over F,, XR, (f.. f,) is the minimizer of L,, (f. t) with y; replaced by Yr.
We then normalize (k( follows
(k(k) = 1+ Σ (1(1)

---

Deep principal support vector machines
[//lo if > 10, and [a] denotes the largest integer not exceeding The impact of eigenvalues is quantified by the function
where rj 1 if 10. r
at a certain convergence rate. Before stating the results, we need the following assumptions. Let [f",t"] be minimizer of L(f,t) over L2 x R.
Assumption 4.2. There exists an absolute constant B > 1 such that || I B and || I B for any f E Fu. Assumption 4.3. f is 3-Hölder continuous, i.e., Iff (x1)) [*(22) S Laial |25 - =22l for any 21, 32 [0, IP, where Litd is some positive constant.
X On (k) 1+ +1
for l- = 0,1,...,,7
Finally, the ladle method determines the dimension using the formula
Assumption 4.4. Nondecreasing convex loss function
d = argmin (f:. + In (+)))
R + Rt is Lipschitz continuous with Lipschitz constant Lisp L.e. lp(x) p(y) S Lipplæ - yl for any 2, y E R. Assumption 4.5. There exist an universal constant C and a parameter V depending on Fn such that
k 0
4. Theoretical Results
For nonlinear SDR, our primary objective is identify functions within the central class Mylx. The population form of (6) is defined as
log Nc. ll ll(x())) S (1
where Q refers to Px and Pn. the empirical probability measure of X.
t(fit) => Var (f(()) + At?
Assumption 4.2 states that the ground-truth function f and neural networks are bounded. The parameter / in Assump- tion 4.3 characterizes the smoothness or regularity of f". A higher value of / indicates greater smoothness, which generally facilitates more accurate approximation. The con- vergence rate of our estimator is jointly determined by 3 and the input dimension p; see Corollary 4.7 for details. Assumption 4.2 and 4.3 are common requirements in many existing works (Györfi et al., 2002; Schmidt-Hieber, 2020; Farrell et al., 2021). Assumption 4.4, crucial for deriving the convergence result, is satisfied by common loss func- tions such as hinge loss and Huber loss. Regarding the assumption 4.5, this upper bound always holds for any VC class F, with the parameter being its VC dimension (see Theorem 2.6.7 in Vaart & Wellner (1996), such as splines, RKHS with polynomial kemel, and deep neural network class. The parameter V is measure of the size (capacity, complexity, or expressive power) of a hypothesis class Fn. which is relative to the trade-off between statistical error and approximation error of our estimator. Classes with larger V have stronger approximation capabilities but also carry a higher risk of overfitting.
+ ((())(((()-
where& & E IR, f belongs to a function class F, and p repre- sents any nondecreasing convex loss function. The follow- ing theorem shows that it can achieve effective nonlinear SDR.
Theorem 4.1. Suppose that
1. F is a dense subset of L2(x)).
2. Y IL XX)
If (f",'" minimizes (7) over F R. then (X) E My x
Linear SDR methods typically rely on the linear conditional mean assumption concerning E (Xln) T X) to ensure unbi- asedness (Li, 1991; Li & Dong, 2009). Remarkably, for the nonlinear method proposed here, such stringent con- ditions are unnecessary. The first condition of Theorem 4.1 can be satisfied by continuous function class, RKHS with a Gaussian kernel, or Lz(Px) itself. In general, we usually choose F Lz2 (x)) directly for population opti- mization. By Definition 2.1, Theorem 4.1 establishes the validity of our nonlinear SDR method (7) since the optimal solution f is unbiased. In other words, we can conclude that (f (X)) S (X)) The significance of Theorem 4.1 lies in its transformation of an abstract nonlinear SDR problem into a specific and solvable optimization problem (7). It is essential to emphasize that our goal is iden- tify functions that capture as many nonlinear relationships between Y and X as possible, rather than approximating specific predetermined functions, such as fo-
When p(x) = max (1 + 0}. optimizing (6) over RKHS and R is similar to kernel principal SVM (Li et al., 201 I). Flexibility in the optimization range is another advantage over classical SDR methods.
Theorem 4.6. Let (f. f) be the minimizer of Ln t) over F. R and (f*, t be the minimizer of L(f.t) over L2(Px) x R. Let 5 > 0. Under the assumptions 4.2, 4.4 and 4.5, in 2 V, the following inequality holds with
probability at least 1 - 46
L(f,E) - L(f",t))
More generally, given some function class F,, below we prove that the sample optimization (6) over F,, x RR can con- verge to the population optimization (7) over T,2 (x)) x IR.
0 V
log + inf ||f - 1 log(1/8)
n1
V
FEF.
5

---

Deep principal support vector machines
Now we focus on the neural network class. By Theorem 2.6.7 in Vaart & Wellner (1996) and Theorem 7 in Bartlett et al. (2019), the log covering number of the ReLU-activated neural network class Fn with respect to (2))) norm for any probability measure Q can be bounded by
5. Numerical Experiments
We evaluate the finite sample performance of the proposed DPSVM method on synthetic data. For comparison, we consider six other nonlinear dimension reduction methods: kernel principal support vector machine (KPSVM) with Gaussian kernel (Li et al., 201 1), generalized sliced in- verse regression (GSIR) (Lee et al., 2013), kernel sliced inverse regression (KSIR) (Yeh et al., 2009), deep dimen- sion reduction (DDR) (Huang et al., 2024), deep sufficient representation learning (DSRL. (Zheng et al., 2022) and GMDDNet (Chen et al., 2024). KPSVM and DSRL are im- plemented based on the algorithm described in the original paper. GSIR and KSIR can be implemented by R package nsdr and MAVE, respectively. DDR and GMDDNet are implemented using the source code provided by the original authors. More details on the computer configuration and hyperparameters are provided in the Appendix.
log N (e, Il (Iaa())) 1 K1
< K2 - SC (og(()) log(1/e)
where K1, K2 are two universal constants, VC(F,) is the VC dimension of the neural network class, C is the depth of networks, and is the total number of network param- eters. Hence the neural network class satisfies the assump- tion 4.5 with V < (og)), Concerning the approxima- tion error, 3-Hölder continuous function f* can be effec- tively approximated by ReL.U networks. Specifically, con- sidering ReLU neural network class Fr, with width N = 30+3 max (pl [ N +) and depth L = 12L 14 Theorem 1.1 in Shen (2020) tells that
Throughout the simulation, we set the training sample size 500 and the input dimension p 10. We first gener- ate the intermediate response H from the following three models:
If- = -(NL) *). SEF
inf
Here N and L stand for arbitrary positive constants deter- mining the width N' and the depth L of neural networks, respectively. With these preparations, we have the follow- ing result on the function class of width-fixer ReLU neural networks.
A: H -011 x2)) + B: H sin ((X7 + xx3x/110) 10)
X? +
C: H = =lg xog v x3 + X3 +
Corollary 4.7. For any arbitrary N E Nt, consider Fn as the scalar-valued ReLU neural network class with width N = =3 3 max (pl NN J, 1) and depth L = 12n 2tr223) + 14+ + 2p. Let 5 > O. Under the assumption 4.2, 4.3 and 4.4, in 2 SC (o)) the following result holds
where the noise ci independent of X and follows standard t-distribution with degree of freedom 4. Then we generate a random variable B {00, 1) by Bernoulli distribution with success probability pr = 0.1. Then we acquire the final response Y = 10H if B = 1 and Y = H otherwise. This design makes Y an outlier with probability pr, increasing the difficulty of nonlinear dimension reduction. Addition- ally, three different distributional scenarios for the predictor vector X (Xi,.... :xp) are involved:
with probability at least 1 48
L(f.)) LLS = 0 ("
20 (n~
) (1 11+10 log(1/8) log l)
Further,
II - (1) =0 -n p+2z
20
1 + log(1/5)) logn).
I: X ~ NOOpp Ip): II : X (Op, III: X2tt(O0, Sp).
Stone (1982) proved the minimax rate n p+23 for nonpara- metric regression under the Hölder continuous assumption. Remarkably, our nonlinear SDR method, leveraging deep neural networks, achieves this optimal rate up to log n. The- orem 4.1 shows that f is a nonlinear SDR function, while Corollary 4.7 demonstrates the convergence of ff to f*, thereby ensuring the rationality of the first step of nonlinear SDR algorithm process in section 3.1. Regarding the way of aggregating the optimal solutions from all different dividing schemes, the second step of the algorithm process is inspired by the kernel PSVM (Li et al., 2011 1). As in that work, we do not delve into the theoretical justification behind it here. In particular, we call our nonlinear SDR method optimized over neural networks with p(x) = max(1 + xx 0} the deep principal support vector machines VDPS
where Ip is the p x p identity matrix, S,, is a p ma- trix with entry Sty 0.5 and N(A, B) is mul- tivariate normal distribution with mean A and covari- ance B, and t (A, B) is multivariate t-distribution with mean A, covariance B and degree of freedom k. The true nonlinear SDR functions f for the model A, B, C are / + X3), sin ((XP X2)m/10) + X? and log(0.1 + Xi + X? + X2. respectively. The struc- ture dimension is assumed to be known as 1. We use the em- pirical distance correlation DR,, (X) f(x)) developed in Szkely et al. (2007) to reflect the merit of the estimated dimension reduction function. The detailed formulation of
6

---

Deep principal support vector machines
empirical distance correlation is given in the Appendix. In brief, A higher distance correlation indicates more accuracy in extracting nonlinear information. The average distance correlations (with standard deviations in parentheses) over 100 simulation runs are shown in Table 1. Our method out- performs other methods in seven of the nine settings, This suggests that our approach is more robust to outliers, largely due the process of bisecting response values (LVR or OVA step in our algorithm).
13
DPSVM with LeNet
DDR with LeNet
Figure I. 2D visualization of MNIST using DPSVM and DDR.
Table 1. Average distance correlations (standard deviations) of dif- ferent methods for all settings over 100 simulation runs. Bold-
faced numbers indicate the best performers.
DPSVM
KPSVM
GSIR
KSIR
DDR
DSRL
GMDDNe:
(1000.01)
0.7700011
0.7500.01)
0.7500002
(.330110)
0.76/0.00
0.7660.01)
0.4
0.7300.1)
0.640.04)
0.570.083
0.620.02)
(.770001)
63210.03)
0.7266.11
06623.04)
BAN0009)
0.700.000
06600.02)
0.55005))
0.560022)
0.440.000
0.150.27)
0.850.17)
0.8:0.17)
05900 (4)
B.II
(.110005)
0.59/0.10)
0.5560.04)
0.250.18)
0.55(.11
0.4400.011
0.65(0117)
I
(255005))
06100.12)
0.570.04)
0.2700.14
0.58(011)
0.57(0.8))
0.500.118)
C11
(17608844)
(.710.15)
0.5800.05)
0.200..29
0.54(0.17)
0.69)..99
0.700.02
C-II
(.780004)
0,730006
0,700.00)
0.350.10
0620012
0.78/0.14
0.580.15)
0810005)
0.760.08)
0.7300.0)
DPSVM with LeNet
DDR with LeNet
MNIST database. The MNIST dataset contains 60, 000 training images and 10, 000 testing images of hand-written digits. Here, we take DPSVM and DDR as the representa- tives of deep methods for 21) visualization. LeNet (Leecn et al., 1998) is used as the basic neural network structure for both DPSVM and DDR. The ladle method in section 3.2 gives the structural dimension estimation d = 9. Then train DPSVM and DDR respectively to obtain nine es- timated dimension reduction functions. Finally, we train the logistic regression on training images with correspond- ing nine-dimensional representations of original features and get the classification accuracy result on testing images. The accuracy 0.9862 and 0.9883 for DPSVM and DDR, respectively. As a reference, the accuracy of direct LeNet using the original features is 0.9899. The closeness of the results suggests that both our method and DDR are effec- tive in extracting nonlinear information from the regression relationship. Lastly, we plot the projections of the original features along the first two estimated dimension reduction directions of DPSVM and DDR to compare visualization effects (see Figure 2). appears that DPSVM exhibits stronger discriminant power compared DDR, as the points for each class cluster more distinetly.
Figure 2. 2D visualization of CRIME using DPSVM and DDR.
methods are applied to the training set to extract a low- dimensional representation. Then, we fit either linear regres- sion logistic regression the representation, depending whether it is regression task classification task. The mean square error for regression and classification accuracy for classification are calculated on the testing set. To reduce randomness, we repeat the above splitting-extir eetreeeetoa-nnittttiiiit fitting process 100 times. The KSIR method consistently incurs errors on the OPTDIGITS dataset across all 100 repetitions, its result is excluded from the plot. Note that the y-axis of the first three plots represents accuracy (higher values are better), while the y-axis of the last plot represents the mean square error (lower values are better). Our method shows comparable performance across all datasets, highlighting its competitiveness. Additionally, we present a 2D visualiza- tion comparing our approach with the DDR method on the CRIME dataset. The results indicate that our method yields more visually coherent representations.
6. Conclusions and Limitations
In this paper, we present a new approach to nonlinear suffi- cient dimension reduction that avoids the issues associated with generalized cigendecomposition of n X n matrices in conventional methods. The approach includes kernll principle SVM (Li et al., 2011) and extends it by incorporat- ing neural networks for more flexible deep nonlinear SDR, particularly beneficial for severe data contamination, large samples and complex inputs. The fusion of deep learning with nonlinear SDR brings new energy to the SDR field. From theoretical standpoint, we provide guarantees of un-
Other datasets. The datasets' involved are Communities and Crime (CRIME) for regression tasks; Connectionist Bench Sonar (SONAR), Optical Recognition of Handwritten Digits (OPTDIGITS), Breast Cancer Wisconsin (WDBC) for additional classification tasks, For each dataset, we ran- domly split the data, using 2/3 for the training set and 1/3 for the testing set. Different nonlinear dimension reduction
All the datasets are downloaded from UCT machine learning repository
7

---

Deep principal support vector machines
Impact Statement
This paper presents work whose goal is to advance the field of sufficient dimension reduction. There are many potential societal consequences of our work, none which feel must be specifically highlighted here.
SONAR
OPTDIGITS
References
Bartlett, P. L.., Harvey, N., Liaw, C., and Mehrabian, A. Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks. Journal of Machine Learning Research, 20(63):1-17, 2019.
Blanchard, G., Bousquet, O., and Massart, P. Statistical per- formance of support vector machines. Annals of Statistics, 36(2):489- 531, 2008
WDBC
CRIME
Figure 3. Performance of different methods on different datasets
Chen, Y., Jiao, Y., Qiu, R., and Yu, Z. Deep nonlinear sufficient dimension reduction. Annals of Statistics, 52 (3):1201 1226, 2024.
Clémençon, S., Lugosi, G., and Vayatis, N. Ranking and em- pirical minimization of U-statistics. Annals of Statistics, 36(2):844 874, 2008
biasedness with respect to the central o-field and establish strong convergence rates for estimation error. Numerical simulations and real data analysis demonstrate the excep- tional performance of our method.
Cook, R. D. Principal hessian directions revisited. Journal of the American Statistical Association, 93(441):84-94 1998.
The limitation of our method is its speed. For real-valued Y, we must first discretize Y into binary variables. Typi- cally, we perform this discretization 10 times with different split schemes and fit neural networks on each discretization. Additionally, discretizing a scalar Y is typically based on quantiles, but calculating quantiles for multivariate is challenging and requires more effort.
Cook, R. D. Fisher lecture: Dimension reduction in regres- sion. Statistical Science, 22(1):1-26, 2007.
Cook, R. D. and Weisberg, S. Discussion of sliced inverse regression for dimension reduction. Journal of the Amer- ican Statistical Association, 86(414):328-3 1991.
Acknowledgements
Farrell, M. H., Liang, T., and Misra, S. Deep neural net- works for estimation and inference. Econometrica, 89(1): 181-213. 2021.
We thank all the anonymous reviewers for their help- ful comments and suggestions. Yinfeng Chen is sup- ported by the National Key R&D Program of China (Grant No. 2023YFA 1008700 and 2023YFA 1008703). Rui Qiu is supported by the Postdoctoral Fellowship Program and China Postdoctoral Science Foundation (Grant No. 2024M760060 and BX20250069). Jin Liu is supported by the National Natural Science Foundation of China (Grant No. (2371283). University Development Fund (Grant No. UDF01003033) from The Chinese University of Hong Kong, Shenzhen, the Program for Guangdong Introducing Innova- tive and Entrepreneurial Teams (Grant No. (0232T10X044)) Shenzhen Science and Technology Program (Shenzhen Key Laboratory Grant No. ZDSYS2023062609 1302006),
Györfi, L., Kohler, M., Krzyzak, A., Walk, H., et al. A distribution theory of nonparametric regression. Springer, 2002.
Hinton, G. E. and Salakhutdinov, mensionality of data with neural networks. Science, 313 (5786):504-507 2006.
Reducing the di-
Huang, J., Jiao, Y., Liao, X., Liu, J., and Yu, Z. Deep di-
mension reduction for supervised representation learning. IEEE Transactions on Information Theory, 2024.
Kapla, D., Fertl, L., and Bura, E. Fusing sufficient di-
mension reduction with neural networks. Computational Statistics & Data Analysis, 167:107390, 2022.
the Guangdong Provincial Key Laboratory of Mathemat- ical Foundations for Artificial Intelligence (Grant No. 2023B 0023B121201000 1212010001, Shenzhen Stability Science Program and Shenzhen Fundamental Research Program (General Program) (Grant No. JCYJ2024081 131 13518024).
Lecun, Y., Bottou, L., Bengio, Y., and Haffner, P. Gradient- based learning applied to document recognition. Proceed- ings of the IEFE, 86(1 1) 1-277-22222 1998.
8

---

Deep principal support vector machines
Lee, K.-Y,, Li, B., and Chiaromonte, F. A general theory for nonlinear sufficient dimension reduction: Formulation and estimation. Annals of Statistics, 41(1):221 249, 2013.
Vaart, A. W. and Wellner, J. A. Weak convergence and empirical processes. Springer, 1996.
Wu, H.-M. Kernel sliced inverse regression with applica- tions to classification. Journal of Computational and Graphical Statistics, 17(3):990-610, 2008.
Li, B. and Dong, Y. Dimension reduction for nonelliptically distributed predictors. Annals of Statistics, 3773):1272- 1298, 2009.
Xia, Y., Tong, H., Li, W. K., and Zhu, L.-X. An adaptive estimation of dimension reduction space. Journal of the Royal Statistical Society: Series B (Statistical Methodol- ogy), 64()):63-4110 2002.
Li, and Wang, On directional regression for dimension reduction. Journal of the American Staristical Associa- lion, 102(479);997-1088, 2007.
Yeh, Y.-.., Huang, S.-Y., and Lee, Y.J. Nonlinear dimen- sion reduction with kernel sliced inverse regression. IEEE Transactions on Knowledge and Data Engineering, 21 (11):1590-1603, 2009.
Li, B., Artemiou, A., and Li, L. Principal support vector machines for linear and nonlinear sufficient dimension reduction. Annals of Statistics, 39(6):31 2011.
Li, K.-C. Sliced inverse regression for dimension reduction. Journal of the American Statistical Association, 4):
Yin, X., Li, B., and Cook, R. D. Successive direction ex- traction for estimating the central subspace in a multiple- index regression. Journal of Multivariate Analysis, 99(8): 1733-1757, 2008
316-327, 1991.
Loyal, J. D., Zhu, R., Cui, Y., and Zhang, X. Dimension reduction forests: local variable importance using struc- tured random forests. Journal of Computational and Graphical Statistics, (144):104-111: 2022.
Zheng, S., Lin, Y., and Huang, J. Deep sufficient repre- sentation learning via mutual information, 2022. URL attps://arxis
Luo, W. and Li, B. Combining eigenvalues and variation of eigenvectors for order determination. Biometrika, 103(4): 875-887,2011
Zhou, J. and Zhu, L. Principal minimax support vector machine for sufficient dimension reduction with contami- nated data. Computational statistics & data analysis, 94: 33-48, 2016
Ma, Y. and Zhu, L. A review dimension reduction. International Statistical Review, 81(1):134-150, 2013.
Massart, Concentration inequalities and model selection. Springer Verlag, 2003.
Schmidt-Hiebed J. Nonparametric regression using deep neural networks with relu activation function. Annals of Statistics, 48(4):1875188) 2020.
Setodji, C. M. and Cook, R. D. K-means inverse regression. Technometrics, (144))421-229, 2004.
Shen, Z. Deep network approximation characterized by number of neurons. Communications in Computational Physics, 28(5): 1768)-1768181,, 2020.
Shin, S. J., Wu, Y., Zhang, H. H., and Liu, Y. Principal weighted support vector machines for sufficient dimen- sion reduction in binary classification. Biometrika, 104 (1):67-81, 2017.
Stone, C. Optimal global rates of convergence for nonpara- metric regression, Annals of Statistics, 1040-1001,, 1982.
Szkely, G. J., Rizzo, M. L.,, and Bakirov, N. K. Measur- ing and testing dependence by correlation distances. Annals of Statistics, 35(6):2769 2794, 2007.

---

Deep principal support vector machines
A. Additional materials
A.1. Illustrative examples
To illustrate the motivation behind sufficient dimension reduction (SDR), we present two simple examples.
EXAMPLE I: CLASSIFICATION PROBLEM
Consider the classification problem:
=
x
where X2 ~ ~N(0, 1) and & -N(00 Then, the normalized SDR function is given by
1o(X) Xi + +
2
Here, a linear transformation is applied to ensure B(f((X)) 0, Var(fo(X)) = 1.
EXAMPLE 2: MULTIVARIATE REGRESSION PROBLEM
Now, consider a multivariate regression problem where E R
Y1 sin sinn(1 + X + +1,,
Y/ cos X4) + 4.2,
where X1, X2, X3, X4 NO 1) and 61 0.2). In this case, the normalized SDR function is
(ag (in + x3) + bu a1 cos (X3 + b1 where a, and b = 2, are introduced to ensure ((fo(X)) =0. O. Var(So(X)) 12-
fo(X)=
EXAMPLE 3: GENERAL FORMULATION FOR REGRESSION
More generally, consider the model
Y= (=g(fo(X)) +e,
where fo IRP + Rd and g RRd R Our goal is identify the low-dimensional nonlinear representation fo(X). In this general setting, traditional linear sufficient dimension reduction breaks down.
A.2. Definition of empirical distance correlation
The empirical distance covariance between two random vectors Y, Y is defined as
DC,, (Y,Y")
71
Σ
h ((Y,, Ya,), (Y Yuz)) VES) (Y Y Vza)).
4
where (Y/) (Y:) are the samples of Y and Y", and
1o (Y4, YI)
Σ
(I - Y1l12 || - -Y'll2
24
1 t
2 2| 11ll2 ||Y! YEll, + +1|1; Ysll, || 1/
Then the empirical distance correlation is defined as
DC, (Y,Y")
DR, (Y,Y')
(8)
DCn (Y DCn (Y'Y")
A.3. Simulations for large p
Since neural networks require more data to fit overparameterized models, we set the sample size to n = 5000 for p = 100. The results, shown Table 2, indicate that our method remains competitive under this configuration.
10

---

Deep principal support vector machines
Table Average distance correlations (standard deviations) of different methods for all settings over 100 simulation runs. Bold-faced numbers indicate the best performers.
Setting
DPSVM
KPSVM
GSIR
KSIR
DDR
DSRL
GMDDNet
A-I
0.77(0.01)
0.72(0.01)
0.74(0.01)
0.76(0.01)
0.55(0.05)
0.67(0.05)
0.61(0.06)
A-ll 0.73(0.01)
0.65(0.01)
0.71(0.01)
0.68(0.01)
0.50(0.04)
0.62(0.05)
0.62(0.04)
A-Ill
0.75(0.01)
0.65(0.01)
0.63(0.01)
0.69(0.01)
0.51(0.06)
0.68(0.03)
0.65(0.03)
B-I
0.70(0.03)
0.45(0.03)
0.06(0.01)
nan(nan)
0.42(0.05)
0.69(0.05)
0.47(0.08)
B-II
0.67(0.03)
0.39(0.03)
0.45(0.03)
0.25(0.19)
0.53(0.09)
0.67(0.04)
0.64(0.06)
B-III
0.68(0.04)
0.41(0.02)
0.45(0.02)
0.35(0.25)
0.54(0.11)
0.69(0.02)
0.63(0.06)
C-I
0.71(0.04)
0.58(0.02)
0.08(0.02)
0.48(0.64)
0.38(0.14)
0.55(0.24)
0.53(0.14)
C-II 0.64(0.16)
0.54(0.02)
0.59(0.02)
0.26(0.05)
0.63(0.04)
0.75(0.03)
0.70(0.06)
C-III
0.73(0.02)
0,59(0.01)
0.55(0.01)
0.25(0.05)
0.61(0.09)
0.73(0.03)
0.71(0.02)
A.4. Analysis of the computational complexity
Let n, P, f, and N represent the sample size, input dimension, number of hidden layers, and maximum width of the hidden layers, respectively. Let the batch size and training epochs be b and t, respectively. The computation overhead of each iteration is primarily attributed to calculating the gradient of that batch. This computation overhead is ObL [p, Considering that we need to update n//b batches for each epoch, the total complexity of each binarization is
Ont max (p,
Additionally, using the LVR procedure involves h binarizations. After performing h binarizations, we need to construct the h h matrix based on the inner product of two h matrices and compute its eigendecomposition. These two operations require (hh2n) and (1h)) computations, respectively. Therefore, the final computational cost is
((hntCmax(p,.V) + h²n + h2).
In contrast, deep methods without binarization incur slightly better computational cost,
OntC max(p,N33).
The primary computational complexity of GSIR/PSVM stems from the inversion of an n x n matrix, which is O(n²) for the naive inversion algorithm. When comparing this with the complexity Ohntt max(p, it is evident that the computational advantage neural networks becomes apparent only when n exceeds a certain threshold. Otherwise, the computational overhead of gradient evaluation becomes more expensive.
The above conclusion is validated by empirical evaluation in Figure 4. Here, deepsvm represents proposed method, psvm corresponds to KPSVM, and deepsdr refers to DDR. The other names correspond to their respective original methods as indicated.
A.5. Computer configuration and network architecture
All the experiments are run computer with 80 Intel(R) Xeon(R) Gold 5218R CPU @ 2. CPU and 251 GB memory. The batch size and epoch number are set as 100 and the optimizer is set as adam with default parameters. This network structure follows an expansion-contractio pattern, a reverse pattern of the bottleneck structure in auto-encoders. The initial increase width to map the original features to higher-dimensionall space, similar to the kernel trick used in a reproducing kernel Hilbert space. The subsequent decreasing width is employed to perform dimension reduction, ultimately obtaining a low-dimensional embedding. The default network is feedforward ReLU network with hidden dimensions
21, 2Pi,2Pr+h P1++ 2D, P1-1, 16,
where D, = |log2 pl + and is the input dimension. Our method performs well with a simple feedforward network, so we did not explore more complex architectures. Specifically, we generate a default neural network with the following Python code.

---

Deep principal support vector machines
1400
deepsvm
psvm
1200
gsir
1000
deepsdr
DSRL
timesseconds
GMDDNet
400
200
1000
2000
3000
4000
5000
6000
7000
8000
sample
Figure Computational time comparison of different methods.
def generate_defaull_nnet (input_ding dim, http:hhhdewwwwwhhhhheeeeeee width-None, last hidden width-16, output_dim-l)
if first__iiii gosst_hidden_widtt is None:
2 ** (int (np.1og2 (input_dim)) + 1)
net Sn.geuuential()
net append (nn. (input_dim, (idett)
net append (nn. Rel.
next_width - fiiss_hiidde_wwide
for in range (incrreaesssc 1):
(mn. (next (next_width, next_widtt width 2))
net apt. next_widthh width = next_widtss width *
for i in range (increasenee 1): (tt.appendmm..iineerr (next_width, next_width // 2))
apt. ()).
next next_width neet_ee_wwiee width // 2
if last_hidden_width is not None:
while (net__whhiii // 2 >= @et_hhidde_wiidth) append Lnn.ri (next_width, next_width // 2))
aet.append (nn. (nn. (e).)
next_width // net aet.ap (nn, hiar (next_width, output_dim))
next_widte
return net
A.6. Sensitivity analysis with different hyperparameters
To evaluate the robustness of our method, we conduct sensitivity analysis on various hyperparameters: binarization count, X, learning rate, and neural network structure. The results presented Table 3. DPSVM is not sensitive to the number of binarizations. The learning rate should be chosen suitably small to ensure optimal performance. As shown in Table 3, smaller values of À lead to improved performance, as it shifts the focus toward dimension reduction rather than variance constraints. Specifically, À should range from 0.01 to 0.1. Different configurations of hidden layers and widths of neural networks have the most significant impact on performance in settings B and Four different structures I-IV are generated by the code
anneraate_effullt_nee (10, last_hidden_wwidth-1 increase_count-2 default generate_ deffull__nee._ie (10, heiethhidee__wiii
increase, (
ineerate_eefall-- (10, last hust_hidden siidden_widhhh6,,, increase increase_count-3] count=3)
ineerate_defultt.og (10,
first_ hidden siidden_widdb=32 last_ hidden, 1iiidddn_widdh=13 increase_countt
12

---

Deep principal support vector machines
The optimal configuration the default network (Setting 1), 10 binarizations to balance performance and computational cost, À= 0.01, and learning rate of 0.001, chosen by cross-validation.
Setting
5
20
Setting
0.01
0.1
10
A-I
0.80(0.01)
(.81()
(.81(0.01)
0.82(0.02)
A-I
0.79(0.02)
0.81(0.01)
0.78(0.02)
0.46(0.08)
A-Il
0.74(0.01)
0.73(0.02)
0.73(0.01)
0.74(0.02)
A-Il
0.72(0.02)
0.75(0.02)
0.70(0.02)
0.32(0.08)
A-llI
0.76(0.01)
0.77(0.01)
0.77(0.01)
0.77(.002)
A-IlI 0.75(0.03)
0.76(0.02)
0.72(0.02)
0.34(0.09)
B-1
0.73(0.04)
0.73(0.05)
0.74(0.04)
0.73(0.04)
B-I
0.76(0.03)
0.74(0.05)
0.68(0.07)
0.31(0.06)
B-ll
0.71(0.04)
0.70(0.03)
0.72(0.04)
0.70(0.05)
B-11
0.76(0.03)
0.71(0.04)
0.61(0.08)
0.27(0.07)
B-IlI
0.74(0.05)
0.75(0.06)
0.76(0.05)
0.75(0.05)
B-IIII
0.79(0.03)
0.75(0.03)
0.68(0.05)
0.34(0.07)
C-1
0.75(0.04)
0.76(0.02)
0.79(0.03)
0.78(0.03)
C-I
0.77(0.07)
0.76(0.04)
0.71(0.03)
0.36(0.08)
C-11
0.77(0.02)
0.77(0.03)
0.77(0.05)
0.78(0.04)
C-I1
0.82(0.02)
0.78(0.04)
0.72(0.04)
0.36(0.07)
C-1II
0.80(0.03)
0.82(0.03)
0.81(0.03)
0.81(0.02)
C-III
0.84(0.03)
0.81(0.03)
0.76(0.04)
0.39(0.08)
Sensitivity analysis with a different number of binarizations.
Sensitivity analysis with different X.
Setting
0.001
0.01
0.1
Setting
II
IV
A-I 0.81(0.01)
0.82(0.01)
0.64(0.12)
0.09(0.05)
A-I
0.79(0.01)
0.80(0.02)
0.81(0.02)
0.80(0.03)
A-II
0.74(0.02)
0.78(0.01)
0.07(0.03) 0.59(0.11)
A-Il 0.74(0.02)
0.74(0.03)
0.75(0.01)
0.74(0.02)
A-III
0.76(0.03)
0.80(0.01)
0.61(0.14)
0.09(0.06)
A-IlI 0.78(0.01)
0.77(0.22)
0.79(0.02)
0.76(0.03)
B-I
0.75(0.06)
0.68(0.04)
0.56(0.16)
0.07(0.02)
B-I
0.74(0.03)
0.54(0.04)
0.59(0.05)
0.57(0.05)
B-ll
0.72(0.05)
0.65(0.04)
0.60(0.08)
0.08(0.04)
B-11
0.72(0.03)
0.54(0.05)
0.52(0.05)
0.54(0.05)
B-III
0.74(0.04)
0.67(0.04)
0.60(0.13)
0.11(0.08)
B-IIII
0.77(0.02)
0.63(0.03)
0.58(0.03)
0.60(0.03)
C-I
0.77(0.03)
0.69(0.04)
0.51(0.21)
0.09(0.06)
C-I
0.75(0.04)
0.62(0.06)
0.63(0.04)
0,60(0.06)
C-II
0.78(0.03)
0.74(0.05)
0.60(0.15)
0.09(0.04)
C-II
0.77(0.04)
0.66(0.05)
0.63(0.04)
0.64(0.04)
C-III
0.77(0.04) 0.80(0.03)
0.61(0.17)
0.07(0.04)
C-III
0.82(0.02)
0.69(0.02)
0.64(0.03)
0.67(0.04)
Sensitivity analysis with a different learning rate.
Sensitivity analysis with a different network.
Table Results of sensitivity analysis with different hyperparameters
B. Proofs
Lemma B.1. Under the assumption for each L E R, the function f > L(J,t) is continuous with respect to the L2 (Px) norm.
Proof.
The proof follows the same idea as Lemma 1 in Li et al. (2011). For any fi, E L2(Px), the following inequality holds
-
<X| Var(fi(X)) Var(fa(X))
+ +Elp(- - E(f(()) -t}) - -(-)-(((((()) - E(f((x)) - +
By Schwarz inequality,
Var(Sa(X))
Var(Ji(X) - /2(X) + (x()) - Var(/2(X))
Var(fi(X) - f2(x)) +2 Varlf(X)) f2(x)) Var(5((X)))112
(X) - +21111))))))))))))1))))))) (A0AI000-22////xxxxXXX00110000
+
By Lipschitz continuity of p.
E C1P(-Y(f(X) (x) - -(((()) +))
<Tuip (X) ( - -
<2Tnp (X)) -f2(x))
13

---

Deep principal support vector machines
22p [I((X) (x())))
Then conclude that the function f I(J,() is continuous with respect the (2x)x) norm.
Proof of Theorem 4.1. Without loss of generality, assume B(S(X)) 0. Note that
Ep(-Y - -
By Jensen inequality and conditional independence Y 1 xl)o(x),
((E]((-7(r0)-44))).00 2 ((((E) - -V(11x)-
(9)
= =E(p( - -((([1())16000
Now apply the law of total variance, then we get
Var(f(X)) Var [=(F(X))f((X)) Var Combining (9) and (10), have
Var (f(X))fo(x)) 2 Var (E(f(X))f(()))
(10)
L(S,)) 2 Var [B(I(X]1F((X)) + X²² + +E[o( ----://))))----(---------------- (00[(((((())666200))))) -EU((()))5(X))
The equality in (10) holds when
E Var (F(X)1f((X)) =0 a.s.
which holds if and only if there version of that measurable with respect to offo(X)). Therefore, El Var (f()) > 0 or ff is not measurable with respect to a (ff(())) then
I(S >> ((g,t)
where g(X) (xx( is function in L2(Px). And F is dense in L2(Px), so can find function g such that
- -lllaP)))(Ps)) Se
for all > 0. By Lemma B.1 and the arbitrariness of e, we find gy E F such that
L(f.t) > L(gg,t). t).
which concludes that any function that is not measurable with respect to a (fo(X)) can't be the minimizer.
Proof of Theorem 4.6. To facilitate theoretical analysis, we define an intermediate quantity
(f.tt =A Var, (f(x)) + 2 + Σ
+})
(LI)
21
[s(X)) - f(x))) +62
n(n - 1) 15K4JSn
Σ p (-Y((((x)) Ef(X) t))
Σ hsa(Z,,
n(n - 1)
18icsšn
14

---

Deep principal support vector machines
where h,, Z3) A62 -(XX))X(XX + - (L
symmetric kernel
and Z; = -Xi,, Yi).
Define 9 2, 2z) = h,, Zy) - -- (Z,, 2z), where (f",t") is the optimal solution of L(f.t) over L2(Px) x R.
Now we consider
(f.t) -m((s))- - (,) n(n - 1)
2
Σ
qf, 2z).
Igi<sšn
It's obvious that the expectation of (f,t) the excess risk -(f
Given the conditional kernel
g / (2:1trr2))dP()))dPPP3 = A(fit),
(15a(²1)) 22)
12,11(511, 22)
then the Hoeffding decomposition of (f,
2 n
(n. = +
Σ
(Z,) -0.5+.}
n
+
1
Σ
{92(Ze, Z,) (..,,((i)) -1,,,,(2))) + qo.
(12)
n(n 1)
ISi$jsn
= A1, + A2,
where (f,t fo.f. + Noting that A,, (f, t) is an independent sum, let new loss function
2 Σー - ( (Z.) -o.r. = E - A(f.t) and (2,n(fit) is the third term.
(ff.t.Z) -((
Now we consider new objective function
n n Σ (f,t
on (f.tt) (1,x
whose minimizer over F, x R is denoted by (fe, t&) and centered empirical form is
1
T, (f.
Zf -(f.t)
Step 1: Loss gap.
The sample estimator (f,, i) is obtained by minimizing Ly, while our convergence analysis turns to the independent sum vn(f t). We should derive the loss gap between Uy (f, è) and e, (fe, te) at first. Let the minimizer of L,, be (f. f), then
0 < - (f, f) L.. (f. (In.f.)) (f.(f.))
= D, + D2-
The term D2 can be bounded by
D2 I..[.. 0) Ino(J,t) In((st) (F.))
5 (n.(fit (F,((
15

---

Deep principal support vector machines
L - -.,
En (f-( - -
I Liip Enn -E(X())))
Noting that f is bounded by B, by Hoeffding's inequality, with probability at least 1 - exp(
2nd²
it holds that
8²
Lasp -E(((())))(((( S Lispö,
Or equivalently,
log(1/5)
D2 S Lii
2n
with probability at least - 8. And D1 can be bounded by the same logic. Therefore, with probability at least 25,
log(1/6)
[f,(f))) - Ln(F, & BV
2n
Therefore with probability at least 1 - -2,, we have
log (og(1/))
(n - tt) - 2L BV
< Ln. - In(fttt)
2n
which implies
in, 2Li log(1/5)
S An((ft) S (fer
since (f,t) is also the minimizer of (f.t... With (f,tt) (f, + (2, t), further have
log(1/5)
AI $ BV
(f, (f. 2n
(13)
Additionally, by corollary 6 in Clémençon et al. (2008), with probability least 1 5 it holds that
(log (2) < (log 1/8)/m), fe
sup
(14)
where C' is a constant and V is defined in Assumption 4.5.
Combining (13) and (14), it holds that with probability at least - 35
log log(1/8)
A - ()) < 2Lis BV
+ 2C(V/n (log 1/6)/n).
2n.
Choosing small enough so that log(1/6) 2n
> 1, then we can conclude that with probability at least 1 - 35,
vn (f.)) - (ee s0V/ (log 1/8)/n).
Step Pseudo-metric.
Define the pseudo-metric ((2,))))) - f2ll2:(s)) + 141 - 12)), where the universal constant C₁ will be determined later to satisfy some properties. Since
(z.f+(Z) =x² + Mf(() - E(f(X)) p(- ([((X))E(((X)) -)) t))
2
2
16

---

Deep principal support vector machines
Epl Epl- ( -(((X)) - -1)
+
2
we have
Var [0(S1,11, 7) 22, 711
=Var [201.f1.112 (7) 22 (7))
= VVa (2)) E(f(())) 1² + 2X2 + - - - +1))
22) - F()2(x))) 2 - - F((((X)) - -2))
Noting that f1, f2 is bounded by B, Var [MA( -(A(-))))) -(((())))))))))))- -( can be bounded by
A [1((X) (X)( 12(x) F(2(())) = (x)(()
-
52() F()(()))))
<32B²>² [(fi(X) -ECC)))))) -(((())))
<64B2>2115 f2ll)2((x)
With Lipschitz continuity of p. the remaining term can be bounded similarly by
20i 2l.fi fall2a(xx) + lt1 -t2l3)
Let Ci = -2iip² 2 + 64 Then
Var [( & (.)(ff,.t).(f2.tt)))
The following deduction is still valid if the zero mean assumption doesn't hold. We can just replace f and f by f f f and ft = f" - Ef". respectively. With loss of generality, we assume EX)) = 0 here. Let f be any function such that E(f(X)) = 0 and be any constant, then
T( I(f",t"
=X - --0------000000----------------- +n))
-x()) + -
+Ep(-8([(X)- - -x(-((((()-+t)
-XE(J())-(XX)} - -11--1++A
(XE(J(X) +
=XI]ʃ - 1 + +N/t -t7.
if
A = -
(15)
20. 0.
Now we prove the condition (15) by contradiction. Let's assume A < 0. By convexity of loss function p. the following holds for any 0 < c < 1,
sxxe F(((X) - d ( - P
L(ef +
+ 2Xct (*)
17

---

Deep principal support vector machines
The root - + -
+ cA
224
0
0.
2X
Therefore for some small < c < 1, holds that
L(of + +1- c) at + +t ,t)
which contradicts with the fact that (f*,t") is the minimizer.
Therefore, we get the lower bound of t(f,t) -(f"tt")
C1 L(, w/ -
dp((i, (22,
X
which implies that w<xx) = a: x: a linear function of x.
Step 3: é bound
Next need to find function o: [0, oo) + [0, oo) with the property
(a) VnE
sup (J,)) - un (1) s {(()),
where S argminger, L(S,). Applying Lemma 6.10 in Blanchard et al. (2008), we get
(o)
(.NN 11 )d
((())))((.)))lll - (
Vn
sa//
V - C.
)d
110))
Vn
<a/ l." continuous coniitttionnoooii
( By Loptiuus
H1x0 log(B/8) + Ho) (E(o))2 0² v2vn"" ) By ((x) 2 x/HHo assume V B)
Vn
<C(E(0)+
=
where H.s (8) = log N88, Fn. |l llxo) and E(x) = Jo / (Ho s C(VVx+ VV log(1/6)))))
Denote x. the solution of equation {(x) = Vna² and a. = cr. for a suitable choice of constant e 2 2C > 2, Since x- ==x((=) a decreasing function, holds that
E(o.) = {(cx,) S c((xx)) = c/nx² vnoi/c.
Plugging this into ((r) yields
(l0.) _ C S vnoi
18

---

Deep principal support vector machines
which implies that ((o.)/o3 < Vū. And d and therefore /(x)/a² is non-increasing function, which conclude that cr, the upper bound for the solution of equation ((r)/0² -vn. Let the solution of C(v vV So log(1/6)db) vna? be v., Similarly,
$ vV //" log(1/8)d6)/zz s vū.
which implies 2,, is a upper bound of the solution of E(x) = vnz² with the non non-increasing property of Ex2
Define the non-increasing function g(x) _e /log(1//)d with property 9(1) 0.88 1/2. Therefore.g(z) > 1/2
and vV S 2
vVg(z) for x (0, With the same logic as in the last paragraph the solution of the equation of V is upper bounded by that of 30 30 Vala) x
x
C( (( v
nn. Finally, the solution of equation ((x) v
of xt upper bounded by the solutions
3C vVg(a) Vn.
x
For arbitrary 0 < 1, we have
log 1/ede/0. fet
Vlog 1/ede
C
n
V
x1+o
V
Le V /ogll/cfdtt It's easy to verify that 1o
and it's a bounded by some constant C. Therefore, we can bound x. by
x1-u
r. Vlog 1/ede
1
C²
V
V
1-7
I
x
lta
+1ta
vo
n
12
where y 1+a
is an arbitrary number in 10.1.2).
Therefore, the upper bound of solution of ((w((c)) = ((..) -vnc is
2= = c2 V 1-7
where y is arbitrary number in (0, 1/2).
Step 4: Final result
Combine all the results above and Theorem 8.3 in Massart (2003), we can conclude that with probability at least 45,
L(f,t)
L(F",t) 1-y
V
log(1/8) +2 inf SEFn,TER
{L(f.t) - L(F",))} + 1 (M()) 1-y
log(1/5)
0
+
n
n
n
1-7
{L(ftt") {t)},,} + 1og(1/6)
0
+
2
inf
feF
7
1-7
f 113 + log(1/8)
<0
inf
||f
fefn
n
Note that the above bound holds for all > E 10, and
C - f' <d,
1). (f*,t)) (( -(f",t")
<0
V
1 1
log(1/8)
+ inf [l((z)) - (1))
(16)
feF
19

---

Deep principal support vector machines
Proof of Corollary 4.7. Given any arbitrary N, L N+, assume that Fn is the deep ReLU networks Fn with width
N 3 max (pLNV//J,, J. + 1)
and depth L 127 + 2p. By Theorem 1.1 in Shen (2020), we have the approximation error
inf |lf - f" = (WVD) 3 ))
(17)
feFu
Since the ReLU-activated neural network class is VC class, the parameter V in assumption 4.5 can be its VC dimension (see Theorem of Bartlett et al. (2019))
lo (5).
(18)
Based on (17) and (18), the following holds with probability at least 45 for any 5 > 0 by Theorem 4.6,
42 - = (SLlog (S) + (NC)
SC (S) log log(1/0))
n
Notice that, for any neural network f E Fn, its size satisfies
8 S 1 (C- +
Therefore,
i) (() N2 log (N + (NC) 14 N2 log (N (1/))
n
n
To achieve the balance between the first two terms, consider the width-fixed ReLU neural network class F, with = n 26-22]] and N being any positive integers, which leads to
LÍID) - -( =0 (n p-27 log ni + n
F+23
28
log n log(11)))
28 =0 vt2
{1 log(1/6)} logn)
By equality (16), |1f - -f- has the same convergence rate with L(f,t) - L(f,tt)
20


Reviews:
[{"review_id": "pdSgx7KrA3", "reviewer": "Reviewer_RtJT", "summary": "The paper introduces a general framework of nonlinear sufficient dimension reduction (SDR) based on the previous works of principal support vector machine (PSVM). Classical algorithms such as linear SDR and kernel PSVM are special cases of this new framework. When applying deep neural network to this framework, we obtain the deep principal support vector machine (DPSVM) algorithm. \nThis paper also provides theoretical guarantees for the effectiveness of this new framework and the DPSVM algorithm, including the conditions for the unbiasedness of the optimal solutions, and the minimax convergence rate for the mean square estimation error under Holder assumption when applying deep neural network to this framework.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the insightful feedback. Below are our responses.\n\n> **Q1** The convergence results (Theorem 4.6 and Corollary 4.7) only analyze the properties of the empirical loss minimizer, while does not put the network training into consideration. Since the training speed may be one of the limitations of DPSVM, the analysis of the network optimization is a critical issue for DPSVM.\n\n**A1**. Thank you for your valuable opinion. The training loss or optimization error is often overlooked in the statistical literature. For instance, Schmidt-Hieber (2020) establishes a theoretical bound for deep neural networks in nonparametric regression, where he defines the gap between the estimator derived from training and the true global empirical minimizer as $\\Delta_n$. However, he does not provide any theorems specifically addressing this term, instead treating it as either negligible or dominated by statistical error.\n\nThe analysis of the training process in neural networks is indeed complex and demands considerable effort. Some recent works, such as Jentzen & Welti (2023) and Beck, Jentzen & Kuckuck (2022), incorporate training error into their analysis of least squares regression. However, as noted by the authors, the resulting convergence rates are far from optimal and are hindered by the curse of dimensionality. While we could adopt a similar decomposition—including statistical, optimization, and approximation errors—such an approach may obscure the core theoretical contribution of our work.\n\n* Schmidt-Hieber, J. (2020). Nonparametric regression using deep neural networks with ReLU activation function. The Annals of Statistics, 48(4), 1875–1897.\n\n* Beck, C., Jentzen, A., & Kuckuck, B. (2022). Full error analysis for the training of deep neural networks. Infinite Dimensional Analysis, Quantum Probability and Related Topics, 25(02), 2150020.\n\n* Jentzen, A., & Welti, T. (2023). Overall error analysis for the training of deep neural networks via stochastic gradient descent with random initialization. Applied Mathematics and Computation, 455, 127907.\n\n>**Q2**  The structure dimension estimation discussed in section 3.2 lacks further theoretical discussion.\n\n**A2**. Yes, you are correct. Structure dimension estimation is inspired by traditional linear sufficient dimension reduction (SDR) methods such as SIR (Li, 1991). Here, we follow the ladle method from Luo & Li (2016)  and construct a positive-definite matrix, $E_n (\\widehat{\\mathbf{f}}^\\intercal(X)\\widehat{\\mathbf{f}}(X))$, which converges to $E ({\\mathbf{f}}^\\intercal(X){\\mathbf{f}}(X))$.  Luo & Li (2016) have demonstrated that the ladle method consistently selects the true structure dimension as $n$ approaches infinity.  We have used this method in artificial datasets, and will subsequently add simulation experiments to verify the effectiveness of structure dimension estimation.\n\n* Li, K.-C. (1991). Sliced inverse regression for dimension reduction. Journal of the American Statistical Association, 86(414), 316–327.\n\n* Luo, W., & Li, B. (2016). Combining eigenvalues and variation of eigenvectors for order determination. Biometrika, 103(4), 875–887.\n\n>**Q3**. The DPSVM algorithm requires the binary discretization of the response variable. For the continuous response variable, the error arising from discretization is not discussed in this paper.\n\n**A3**. If we discretize $Y$ into $\\tilde{Y}$, our goal is to find $f$ such that $\\tilde{Y} \\perp X \\mid f(X)$.  It is known that  \n$$Y \\perp X \\mid f(X) \\implies \\tilde{Y} \\perp X \\mid f(X),$$\nbut the reverse does not necessarily hold. However, both functions $f$ belong to the central class, which is fundamental for sufficient dimension reduction (SDR).  In nonlinear SDR, the central class forms an infinite-dimensional space when $X$ is real-valued. Since fully recovering this space is infeasible, we can only approximate it. Therefore, the error introduced by discretization isn't a big issue here.\n\n>**Q4**. This paper does not provide experimental results to compare the speed of DPSVM with those of other algorithms.\n\n**A4**. Theoretical computational costs can be derived, but for brevity, we present only the final result. \nLet $n, p, \\mathcal{L}, \\mathcal{N}$ be the sample size, input dimension, number of hidden layers, and maximum layer width, respectively. With batch size $b$ and epochs $t$, our method has a total cost of\n$$\\mathcal{O}(hnt\\mathcal{L} \\max\\\\{p, \\mathcal{N}\\\\}^2 + h^2n + h^3),$$\nwhere $\\mathcal{O}(nt\\mathcal{L} \\max\\\\{p, \\mathcal{N}\\\\}^2)$ accounts for training per binarization, $\\mathcal{O}(h^2n)$ arises from matrix multiplication in $E_n (\\widehat{\\mathbf{f}}^\\intercal(X)\\widehat{\\mathbf{f}}(X))$, and $\\mathcal{O}(h^3)$ is due to eigendecomposition of an $h \\times h$ matrix.  \nIn contrast, deep methods without binarization incur\n$$\\mathcal{O}(nt\\mathcal{L} \\max\\\\{p, \\mathcal{N}\\\\}^2).$$\nWe will add relevant experiments and theoretical complexity analysis in the revision.", "claims_and_evidence": "The statements of this paper are clear. The new algorithm introduced by this paper is based on a series of previous works, and a number of references and experimental data are provided.", "methods_and_evaluation_criteria": "The nonlinear SDR framework introduced in this paper has the following strengths: \n1. It is very flexible and can be applied to various types of model spaces, including reproducing kernel Hilbert spaces and deep neural network spaces. \n2. In contrast with many previous works on PSVM, the new framework allow the response variable to be multi-dimensional, implying that it is suitable to more various types of regression and classification tasks. \nThere are a number of limitations in this paper: \n1. The convergence results (Theorem 4.6 and Corollary 4.7) only analyze the properties of the empirical loss minimizer, while does not put the network training into consideration. Since the training speed may be one of the limitations of DPSVM, the analysis of the network optimization is a critical issue for DPSVM. \n2. The structure dimension estimation discussed in section 3.2 lacks further theoretical discussion; \n3. The DPSVM algorithm requires the binary discretization of the response variable. For continuous response variable, the error arising from discretization is not discussed in this paper.", "theoretical_claims": "We have checked the proofs and found no evident mistakes.", "experimental_designs_or_analyses": "The paper provides adequate experimental results for the comparison between DPSVM and other nonlinear SDR algorithms on several artificial datasets such as MNIST, demonstrating impressive competitiveness of the new algorithm. \nAs a potential limitation of the experiment part, section 6 mentioned that one of the limitations of DPSVM is the speed, while this paper does not provide experimental results to compare the speed of DPSVM with those of other algorithms.", "relation_to_broader_scientific_literature": "None.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "None.", "other_comments_or_suggestions": "None.", "questions_for_authors": "The DPSVM algorithm discussed in this paper requires the discretization of the response variable, which may cause deceleration of training speed. Is there any nonlinear SDR algorithms that can directly deal with continuous response? If it exists, how is the comparison between it and DPSVM?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Nl4KzaM4HV", "reviewer": "Reviewer_Eu69", "summary": "This paper introduce a unified framework for nonlinear sufficient dimension reduction method based on classification ensemble. The framework proposed in this paper almost includes kernel principal SVM, which is a nonlinear sufficient dimension reduction method using the reproducing kernel Hilbert space and SVM. Here, the neural network function class is considered to replace the reproducing kernel Hilbert space to implement a more flexible deep nonlinear sufficient dimension reduction. The authors demonstrate theoretical unbiasedness of the optimal solution of the population objective function and a non-asymmetric upper bound for the estimation error. All these results demonstrate considerable competitiveness of the newly proposed deep nonlinear sufficient dimension reduction method.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for acknowledging our theoretical contribution. We will address your concerns in the revision and below are our detailed responses.\n\n> **Q1**. In the Preliminary part Section 2.1, when the authors introduce the nonlinear dimension reduction problem, it would be better if there could be some concrete example, as this introduction seems very abstract.\n\n**A1**. To illustrate the motivation behind sufficient dimension reduction (SDR), we present two simple examples.  \n### Example 1: Classification Problem  \nConsider the classification problem:  \n$$  \nY = \\mathbb{1} \\left( X_1^2 + X_2^2 - 1 + \\varepsilon > 0 \\right),  \n$$  \nwhere  $X_1, X_2 \\sim \\mathcal{N}(0,1)$ and $\\varepsilon \\sim \\mathcal{N}(0, 0.2)$.  Then, the normalized SDR function is given by  \n$$  f_0 (X) = \\frac{X_1^2 + X_2^2}{2} - 1.  $$  \nHere, a linear transformation is applied to ensure $E\\big(f_0 (X)\\big)=0，\\operatorname{Var}\\big(f_0 (X)\\big)=1$.  \n### Example 2: Multivariate Regression Problem  \nNow, consider a multivariate regression problem where $Y \\in \\mathbb{R}^2$:  \n$$  Y_1 = \\sin(X_1 + X_2^2) + \\varepsilon_1, \\quad Y_2 = \\cos(X_3 X_4) + \\varepsilon_2,  $$  \nwhere $X_1, X_2, X_3, X_4 \\sim \\mathcal{N}(0,1)$ and $\\varepsilon_1, \\varepsilon_2 \\sim \\mathcal{N}(0, 0.2)$.  In this case, the normalized SDR function is  $$  \\boldsymbol{f}_0(X) = (a_0 \\sin(X_1 + X_2^2) + b_0, a_1\\cos(X_3 X_4) + b_1),  $$  \nwhere $a_i$ and $b_i, i=1,2$, are introduced to ensure $E\\big(\\boldsymbol{f}_0 (X)\\big)=0，\\operatorname{Var}\\big(\\boldsymbol{f}_0(X)\\big)=I_2$.  \n### General Formulation for Regression\nMore generally, we consider the model  \n$$  \nY = \\boldsymbol{g}(\\boldsymbol{f}_0(X)) + \\boldsymbol{\\epsilon},  \n$$  \nwhere $\\boldsymbol{f}_0: \\mathbb{R}^p \\to \\mathbb{R}^d$ and $\\boldsymbol{g}: \\mathbb{R}^d \\to \\mathbb{R}^q$. Our goal is to identify the low-dimensional nonlinear representation $\\boldsymbol{f}_0(X)$.  In this general setting, traditional linear sufficient dimension reduction breaks down.\n\n>**Q2**. I wonder if the authors can present some concrete meanings of Assumptions 4.2-4.5 listed in Section 4 Theoretical Results. They seem to be very theoretical.\n\n**A2**. **Assumption 4.2** is a standard boundedness assumption commonly used in statistical learning theory. It states that both the ground-truth function and the neural network are bounded.  \n\n**Assumption 4.3** assumes that $f^*$ is $\\beta$-Hölder smooth. A function $f: [0,1]^p \\to \\mathbb{R}$ is $\\beta$-Hölder continuous if there exist $\\beta \\in (0,1]$ and $\\zeta > 0$ such that  \n\n$$  \n\\left| f(x_1) - f(x_2) \\right| \\leq \\zeta \\|x_1 - x_2\\|_2^{\\beta}, \\quad \\text{for any} \\quad x_1, x_2 \\in [0,1]^p.  \n$$  \n\nSince $\\|x_1 - x_2\\|_2^{\\beta_1} \\geq \\|x_1 - x_2\\|_2^{\\beta_2}$ for $\\beta_1 \\leq \\beta_2$ and $\\|x_1 - x_2\\|_2 \\leq 1$, it follows that functions with smaller $\\beta$ exhibit more pronounced variations in local regions. In other words, the parameter $\\beta$ quantifies the smoothness or regularity of the function:  \n\n- A higher $\\beta$ corresponds to a smoother function with fewer abrupt changes.  \n- A lower $\\beta$ indicates a less smooth, more irregular function with greater variation.  \n\nThe smoothness level $\\beta$ determines how easily the function can be approximated, with larger $\\beta$ indicating easier approximation. Finally, $\\beta$ and the input dimension $d$ jointly determine the convergence rate of our estimator, please see Corollary 4.7.\n\n**Assumption 4.4** requires that the convex loss function is Lipschitz continuous, a condition satisfied by the hinge loss.\n\n**Assumption 4.5** states that \n\n$$\\log N\\left(\\epsilon, \\mathcal{F_n},\\|\\cdot\\|_{L_2(Q)}\\right) \\leq C V \\left(1+\\log (1 / \\epsilon)\\right). $$\n\nHere, $V$ can be regarded as the VC dimension of a certain VC class $\\mathcal{F}_n$, which always satisfies this assumption. For example, the VC dimension of the neural network class is determined by parameters such as the width and depth of the network. The VC dimension is a measure of the size (capacity, complexity, or expressive power) of a hypothesis class $\\mathcal{F_n}$:  \n\n- A higher VC dimension indicates a more powerful hypothesis class, which can lead to overfitting.  \n- A lower VC dimension reduces flexibility but may lead to underfitting. \n\nHence the performance of our estimator is relative to the trade-off between statistical error $V/n$ and approximation error $\\inf_{f \\in \\mathcal{F_n}} \\|f - f^*\\|_{\\infty}^2$, both affected by the complexity $V$ of $\\mathcal{F_n}$: \n\n- A larger, more complex class $\\mathcal{F}_n$ increases $V$, leading to lower approximation error but higher statistical error.  \n- A smaller, simpler class $\\mathcal{F_n}$ decreases $V$, reducing statistical error but potentially increasing approximation error.  \n\nThe result of the trade-off is formalized in Corollary 4.7.", "claims_and_evidence": "Yes they are.", "methods_and_evaluation_criteria": "Make sense.", "theoretical_claims": "Yes. They are solid.", "experimental_designs_or_analyses": "Experiments are good", "relation_to_broader_scientific_literature": "N/A", "essential_references_not_discussed": "No.", "other_strengths_and_weaknesses": "Strengths: Provides both theoretical and empirical soundness of the proposed method.\n\nWeaknesses: In the Preliminary part Section 2.1, when the authors introduce the nonlinear dimension reduction problem, it would be better if there could be some concrete example, as this introduction seems very abstract.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "I wonder if the authors can present some concrete meanings of the Assumptions 4.2-4.5 listed in the Section 4 Theoretical Results. They seem to be very theoretical.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {"ethical_review_concerns": "N/A"}}, {"review_id": "E0b09zktdw", "reviewer": "Reviewer_dTHf", "summary": "The paper focuses on SDR replacing functions in RKHS with neural networks.\nThis is motivated by the potential advantages of neural networks in handling complex data structures. The authors theoretically demonstrate the unbiasedness of their method and provide a non-asymptotic upper bound for the estimation error, using principles of classification ensembles for nonlinear SDR. While the core idea is sound, the paper lacks a truly comprehensive and balanced approach, especially on the complexity and practical considerations, failing to show the true value of the proposal.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks to the reviewer for helpful comments. We will refine the manuscript.\n\n> **Q1**. The paper doesn't provide a detailed analysis of the computational complexity vs. kernel methods. \n\n**A1**. A brief theoretical analysis of computational costs is provided in **Q4** of Reviewer 1. The cost of our methods is\n$$\\mathcal{O}(hnt\\mathcal{L} \\max\\\\{p, \\mathcal{N}\\\\}^2 + h^2n + h^3),$$  which is linear with $n$.\nOther deep methods avoid binarization, reducing it to\n$$ \\mathcal{O}(nt\\mathcal{L} \\max\\\\{p, \\mathcal{N}\\\\}^2).$$  \n\nKernel-based SDR methods, requiring the eigendecomposition of an $n \\times n$ matrix, incur\n$$\\mathcal{O}(n^3).$$  \nWhile discretization increases computation, it enhances robustness for both regression and classification tasks, particularly for data with many outliers.  \n\n>**Q2**. The experiments are limited in scope.\n\n**A2**. Thanks for your comments. Our focus is to provide a theoretical perspective on a new robust nonlinear SDR framework, especially considering the embedding of neural networks. Importantly, we establish a faster non-asymptotic convergence rate compared to Zheng et al. (2022) and Huang et al. (2024). Our rate nearly reaches the minimax rate of nonparametric regression.\nThe existing experiments show results comparable to other approaches. More experiments like simulation with large $p=100$ is added in the revision.\n\n* Zheng, S., Lin, Y., & Huang, J. (2022). Deep Sufficient Representation Learning via Mutual Information. arXiv preprint.\n* Huang, J., Jiao, Y., Liao, X., Liu, J., & Yu, Z. (2024). Deep dimension reduction for supervised representation learning. IEEE Transactions on Information Theory.\n\n>**Q3**. Literature.\n\n**A3**. We have included the latest deep dimension reduction methods for comparison.\n\n>**Q4**. How was the specific neural network architecture chosen, and what was the rationale behind this choice? Did you explore other architectures?\n\n**A4**. This structure follows an expansion-contraction pattern, similar to a reverse bottleneck in autoencoders: widening layers map $x$ to a higher-dimensional space (akin to RKHS feature mapping), while narrowing layers extract a lower-dimensional representation. The default network is a feedforward ReLU network with hidden dimensions $2^{D_1}, 2^{D_1+1}, 2^{D_1}, 2^{D_1-1}, \\dots, 16, 1$, where $D_1 = \\lfloor \\log_2 p \\rfloor + 1$ and $p$ is the input dimension. Our method performs well with a simple feedforward network, so we did not explore more complex architectures. \n\n>**Q5**. Could you elaborate on the hyperparameter tuning process? What range of values was considered for each hyperparameter, and how was the optimal configuration determined?\n\n**A5**. The hyperparameter tuning process is primarily introduced in section A.3. Some hyperparameters, such as optimizer(Adam), batch size (100) and the number of epochs(100), remain unchanged, as they are less significant compared to the others: binarization count, $\\lambda$, learning rate, and neural network structure. The range of them can be found in Table 2. Table 2 shows that increasing the binarization count improves performance, consistent with classical PSVM results. A smaller $\\lambda$ also leads to better performance, as it shifts the focus toward dimension reduction rather than variance constraints. Notably, $\\lambda$ is not a regularization parameter like in $L_2$-regularized regression. The ablation study on learning rate and network structure suggests that the default setting is optimal. The optimal configuration uses the default network, 10 binarizations to balance performance and computational cost, $\\lambda = 0.01$, and a learning rate of 0.001, chosen by cross-validation.\n\nFollowing the reviewer's suggestion, we will conduct a more in-depth investigation into the sensitivity of the proposed method to different hyperparameters.\n\n>**Q6**. What are the practical limitations of the proposed method in terms of computational resources when applied to very high-dimensional datasets?\n\n**A6**. The theoretical computational cost, addressed in **A1**, scales quadratically with the input dimension $p$, which is consistent with most neural network-based methods. Consequently, the practical limitations of our approach are similar to those of typical deep learning methods.\n\n>**Q7**. How can the learned representations be interpreted in a way that provides meaningful insights into the relationship between the predictors and the response variable? \n\n**A7**. For classification tasks like MNIST, the 2D visualization in Figure 1 offers meaningful insights into the relationship between predictors and the responses. Additionally, we will include scatter plots where colors correspond to Y (e.g., brighter red for larger Y and light green for smaller Y). These visualizations will provide a clearer understanding of the learned representations. Additionally, we appreciate your suggestion—conducting feature importance analysis will further enhance the interpretability of our method.", "claims_and_evidence": "The paper claims that using neural networks offers advantages over traditional kernel methods in SDR. However, it doesn't provide a detailed analysis of the computational complexity vs. kernel methods. While it's generally stated that neural networks can be more efficient for large datasets, the specific computational costs are not comprehensively quantified and compared. Moreover, a comparison of the number of parameters should be done.\n\nThe paper's assertion of superior performance lacks substantial supporting evidence. Little empirical evidence might be the greatest weakness of the paper. The simulations and real data analysis are limited in scope and don't fully validate the claimed advantages, as the scope is very specific. The results presented, while showing some improvements, are not conclusive enough to definitively establish the superiority of the proposed method in a broader context.", "methods_and_evaluation_criteria": "The evaluation is limited by the scope of the datasets used. A wider range of benchmark datasets, especially higher-dimensional ones, and a more rigorous analysis of computational complexity would further strengthen the evaluation. The datasets utilized are relatively low-dimensional, and the performance on these might not generalize to more complex, high-dimensional data encountered in modern scientific applications. The lack of comparison with state of the art on high dimentional or more datasets limits the work. The lack of ablation study in order to understand the relative importance of different components in the network is an issue.", "theoretical_claims": "The theoretical claims regarding the unbiasedness of the optimal solution and the non-asymptotic upper bound for the estimation error seem mathematically OK.", "experimental_designs_or_analyses": "The experimental design is limited in its ability to fully support the claims made in the paper. The claim of outperformance needs qualification. The improvements over other methods in Table 1 are sometimes marginal, and the experiments are limited in scope (low dimensionality, specific datasets).", "relation_to_broader_scientific_literature": "The relevant literature is discussed, placing the proposed method within the context of existing SDR techniques. The paper appropriately cites key works in the field of SDR, including both linear and nonlinear approaches. However, the discussion could be expanded to include a more critical comparison with a wider range of contemporary methods, particularly those that also leverage neural networks or deep learning for dimension reduction.", "essential_references_not_discussed": "N\\A", "other_strengths_and_weaknesses": "A significant weakness is the potential lack of interpretability inherent in using neural networks. Given that they utilize neural networks, interpretability can be reduced, and there should be some discussion around this. While neural networks can offer flexibility and representational power, they often come at the cost of being \"black boxes.\" The paper doesn't address how the learned representations can be interpreted or understood in the context of the original data, which is crucial for many scientific applications where domain knowledge is important.\n\nThe impact statement is missing.", "other_comments_or_suggestions": "The paper could benefit from a more thorough investigation of the sensitivity of the proposed method to different hyperparameter settings. The current sensitivity analysis in the supplementary material is limited. Exploring a wider range of architectures and regularization techniques would provide a more complete picture of the method's robustness. Additionally, a discussion of potential strategies for mitigating the \"black box\" nature of neural networks in this context would be valuable. For instance, exploring techniques like attention mechanisms or feature importance analysis could enhance the interpretability of the learned representations. Finally, it would improve the presentation to include a figure illustrating the proposed neural network architecture.", "questions_for_authors": "How was the specific neural network architecture (number of layers, number of neurons per layer, activation functions) chosen, and what was the rationale behind this choice? Did you explore other architectures, and if so, how did their performance compare?\n\nCould you elaborate on the hyperparameter tuning process? What range of values was considered for each hyperparameter, and how was the optimal configuration determined?\n\nWhat are the practical limitations of the proposed method in terms of computational resources (training time, memory usage) when applied to very high-dimensional datasets?\n\nHow can the learned representations be interpreted in a way that provides meaningful insights into the relationship between the predictors and the response variable? Are there any plans to incorporate techniques to improve interpretability?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduced a unified framework for nonlinear sufficient dimension reduction based on classification ensemble. Under the proposed framework, kernel principal SVM is extended by replacing kernel function with neural network to implement a more flexible deep nonlinear sufficient dimension reduction method. Reviewers consistently praised the soundness of theoretical results and competitive experimental results. However, reviewers also raised some key concerns including the lack of empirical computation complexity analysis and the datasets used are not high-dimensional, so the applications of the proposed method to real datasets can be limited. The authors’ responses try to clarify some of these concerns such as providing theoretical computation analysis of the proposed method instead of empirical comparisons, but not satisfactory especially for the concern of high-dimensional data (p=100 simulation data is not considered as high-dimensional). As a result, weak accept is recommended.
Decision: Accept (poster)