Random-projection ensemble dimension reduction
Authors: Wenxing Zhou
Timothy Ivor Cannings
Venue: ICLR 2026 Poster (2026)
Keywords: High-dimensional
random projection
sufficient dimension reduction
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We propose a method for dimension reduction, based on aggregating an ensemble of carefully chosen projections.

Abstract: We introduce a new, flexible, and theoretically justified framework for dimension reduction in high-dimensional regression, based on an ensemble of random projections. Specifically, we consider disjoint groups of independent random projections, retain the best projection in each group according to the empirical regression performance on the projected covariates, and then aggregate the selected projections via singular value decomposition. The singular values quantify the relative importance of corresponding projection directions and guide the dimension selection process. We investigate various aspects of our framework, including the choice of projection distribution and the number of projections used. Our theoretical results show that the expected estimation error decreases as the number of groups of projections increases. Finally, we demonstrate that our proposal consistently matches or outperforms state-of-the-art methods through extensive numerical studies on simulated and real data.

Published as a conference paper at ICLR 2026
RANDOM-PROJECTION ENSEMBLE DIMENSION REDUCTION
Wenxing Zhou, Timothy I. Cannings
School Mathematics and Maxwell Institute for Mathematical Sciences University of Edinburgh
W. 1.zhoo2211 21e @mss. ed. ac.uk, timochyycamiingseed.ee
ABSTRACT
We introduce new, flexible, and theoretically justified framework for dimension reduction in nigh-dimensional regression, based on an ensemble of random pro- jections. Specifically,
consider disjoint groups of independent random projec- tions, retain the best projection in each group according to the empirical regression performance on the projected covariates, and then aggregate the selected projec- tions via singular value decomposition. The singular values quantify the relative importance of corresponding projection directions and guide the dimension se- lection process. We investigate various aspects of our framework, including the choice of projection distribution and the number of projections used. Our theo- retical results show that the expected estimation error decreases the number of groups of projections increases. Finally, we demonstrate that our proposal consis- tently matches or outperforms state-of-the-ar methods through extensive numeri- cal studies on simulated and real data.
INTRODUCTION
In regression, we seek to understand the relationship between a p-dimensional vector of predictors X IR" and a response Y R based on n observations. In many modern datasets, p large and may even exceed the sample size n. In these high-dimensional problems, many classical methods suffer from the curse of dimensionality and can become intractable. For instance, ordinary least squares is not applicable when p n because the sample covariance matrix is singular. Moreover, nonparametric methods, which are often based on smoothness assumptions, hinge on having suffi- cient data around the points of interest for prediction. As a result, the sample size required for such methods scales exponentially with p (e.g.. Wainwright 2019 Chapter 1.2).
One general strategy to address the curse of dimensionality is to seek low-dimensional represen- tation of the predictors that preserves all the information in X about Y. In particular, sufficient dimension reduction (SDR) (Cook) 2007) aims to find function f : R -> Rd for some d p. such that X and Y are conditionally independent given f(X). A closely related problem focuses on the so-called central mean subspace (Cook & Li. 2002). where the conditional mean of Y given X can be written in terms of f(X). In practice, f is often restricted to be linear, which yields more interpretable representations.
OUR CONTRIBUTIONS
We propose random projection ensemble dimension reduction (RPEDR) (Algorithm I]. a novel framework that is distinct from existing SDR methods. This is supplemented by two further algorithms: Algorithm[]seeee
an appropriate projection dimension and Algo- further refinement in some cases (see the diagram in Figure 2
2. We systematically study the key tuning choices in our framework (e.g., projectio distribu- tions, number of random projections) and provide default recommendations.
3. We establish theoretical results for our proposal (Theorem T and Theorem 2). and demon- strate its strong performance through extensive simulated and real data experiments.

---

Published as a conference paper at ICLR 2026
I The main idea in Algorithm is to apply many low-dimensional random projections to the data, fit a base regression model after each projection, and choose good projections based on their empiri- cal performance. The selected projections (SVD). I our random projection-basee domly chosen 2-dimensional projection retains little of the relationship between X and Y, whereas a projection carefully selected from group of M 200 random projections yields more promising results. By repeating this process and aggregating the chosen projections, we are able retain large amount of the structure observed after applying the oracle projection.
then aggregated through singular value decomposition
method on a toy regression problem. A ran-
A5s mi
Figure 1: Dimension reduction example (n = 100, = 20). True model: Y Y= A + 0.3e with A (1, 0,...,.. (.....0)/// /V2 TR20, X ~ N20 and ~ (V001) 1). Left to right: against the oracle projection 2t 22 (red line: true regression function); Y against a random projection; Y against best projection among a group of 200 random projections; Y against the projection estimated by our proposal (with I. = 200).
Our proposed framework is highly flexible and allows user-specified choices of the random projec- tion distribution and the base regression method, depending on the problem at hand - these choices, along with other practical aspects, are investigated in Section Appendix [A] In particular, The- orem I in Section 3.2 shows that the expected estimation error converges at rate no slower than L-1/2 1/2 where L is the number of projection groups This rate is observed in empirical investigations.
One key feature of our method is that the singular values in the SVD step provide measure of the relative importance of the corresponding singular vectors. These values can therefore be used to determine the projection dimension if needed. In Section 3.3] propose Algorithm [2] select the projection dimension by comparing the singular values obtained under our method with those ob- tained when the group size is (no selection). We also explore the potential benefits of applying Algorithm [ second time on the projected data from the initial application in Section 3.4 This refinement step. formalized Algorithm 3 can improve performance certain settings by com- bining signals and removing irrelevant directions. A more detailed demonstration of Algorithm 3 together with the associated theoretical result (Theorem D provided in Appendix A.5
Algorithm (doublc RPEDR)
Algonthm on [s., 2h.
dingi D)
Algorithm
Algoniiaas Will (sccond pass)
Output Às
Yes
d prespecitte??
Output (Ui..... Dir)
No
Algunthm on (6..5
diagi D)
Algorithmm
Ourput (Ua..... Ugg
Figure 2: Workflow of Algorithms 1-3.
Section 4 presents the empirical performance of our proposal on simulated data. We compare our method with several other state-of-the-arr approaches to SDR. Our proposal is very competitive across wide range of settings: it achieves the best results in fifteen out of eighteen cases in the simulation study, as well as in six out of ten scenarios the real-data analysis (see Appendix D) ranking second in all remaining cases. Moreover, in some settings, that demonstrates non-trivial performance. Notably, in cases where > n, some comparators become intractable, whereas our method remains effective in such regimes.
method is the only approach
2

---

Published as a conference paper at ICLR 2026
The remainder of our paper is as follows: an overview of the related work is presented in Section In Section 2 we introduce the statistical setting and present our main algorithm. Section focuses on practical considerations and provides guidance on selecting method inputs. The main simulation sults are presented in Section[4 We conclude with a discussion in Section[5] The appendices contain additional details on practical considerations (Appendix (A), the theoretical proof' (Appendix [B), the full results of our simulation study (Appendix C and real-data experiments (Appendix (D).
1.1 RELATED WORK
Work on sufficient dimension reduction dates back to the idea of inverse regression, which targets the expectation of X given Y rather than regressing Y on X. The earliest work is sliced inverse re- gression (STR) (Li 1991b. which uses 'slices' of Y estimate Cov(X Y). the eigenvectors of this estimate are then used as the basis of the dimension reduction subspace. This idea led to long line of work, including sliced average variance estimation (SAVE) (Cook & Weisberg, 1991], principal Hessian directions (PHD) (Li 1992), sparse SIR (Li Nachtsheim 2006) and directional regression (DR) (Li & Wang 2007). In another line of work, the nonparametric minimum average variance estimation (MAVE) (Xia et al., 2002) estimates the dimension reduction subspace by constructing a minimization problem that quantifies how well the model fits the observations. Comprehensive reviews are provided by Ma & Zhu (2013). More recent proposals include gradient-based kernel di-
mension reduction (gKDR) (Fukumizu & Leng 2014), dimension reduction and MARS (drMARS)
(Liu et al. 2023), robust SDR (Huang et al.) 2024a), sparse kernel SDR (Liu & Xue 2024) and
Neural Network based SDR (Xu & Yu 2025) see also the survey by Ghojogh et (2021)
Random projections offer a natural approach to dimension reduction. The Johnson LiiissonnLnddentruuus lemma (Johnson & Lindenstrauss 1984 Dasgupta & Gupta 2003) shows that applying a random
projection can approximately preserve pairwise distances between observations. However, as il- lustrated in Figure [ naively applying single random projection typically fails. Our proposal is motivated by Cannings & Samworth (2017), who introduced a random- projjeecioo ensemble method for classification by aggregating the results of applying a base classifier to the projected data from many carefully selected random projections. Related ideas have been used in sparse principal com- ponent analysis (Gataric al.] 2020), sparse sliced inverse regression (Zhang al. 2025), and semi- supervised learning (Wang et 2025). Other works on (unsupervised) dimension reduction using random projection techniques include [Binghmm & Mannila (2001) and Reeve et al. (2024). Work
on other problems includes precision matrix estimation (Marzetta et al.] (2011). two-sample testing
(oopes al. 2011]. clustering (Fern & Brodley 20033 Anderlucci al. [0022) high-dimensional
regression (Thanei et al. 20177] Slawski, 2018 Dobriban & Liu| [019]] Ahfock et al. 22020), and
differential privacy (Huang et al. 2024b). See also Xie et a2.17 and (annings(2022 (2021).
STATISTICAL SETTING AND METHODOLOGY
Let P denote distribution on R IR and suppose that the covariate-responss pair (X,Y) P. Let n(x)== := E(Y X = x) be the regression function. We seek to find a dimension d E {l,...,p- and corresponding projection matrix A E Apxd {A E Rxxdd AAA Laxd] for which there exists function g : Rd R satisfying
7()) = (ATTx).
(1)
Here A maps X to d-dimensional subspace, and the representation A X contains all the infor- mation available in X about the conditional mean of Y given X. Equation T is of course satisfied by taking d p. A Ipxp. and g = n7 The main interest, however, in finding a solution to equa- tion for a minimal choice of d, which we will denote by dg- The corresponding projection matrix and link function will be denoted by Ag and go, respectively. While dg unique, the associated Ag is not. Indeed, if Ao Apxd, and go satisfy equation] I then so do A₁ = B and gi (-)- go (B)..), for any B E Adgxds-
Since Ao is not uniquely identifiable, our focus is on the space spanned by the columns of Ao. For a projection matrix A E Apxd, we write S(A) = span(A) for the d-dimensional subspace spanned by the columns of A. If A satisfies equation I then S(A) is called a mean dimension reduction subspace (Cook& & Li (2002). The space spanned by the columns of the do dy-dimensional projection An. S(Ao), is called the central mean subspace (CMS).
3

---

Published as a conference paper at ICLR 2026
For the purposes of this paper, assumee the existence and uniqueness of the CMS. This assumption is standard and mild (see, for example, Cook, 1996] Lemma 1), and allows us to focus on our main goal of estimating the space S ((Ao) based on dataset of n independent and identically distributed pairs (X1, Yi),.... (Xa, Yn) ~ P. Our estimate in this problem will be the space spanned by the columns of a projection Au E A pxdg for some du E [p], where do and Ag are (possibly randomised) functions of the data D Yi)......... (X,, Y.)). ) To measure the distance between the spaces S(Ao) and (Aa), we will focus on di (S(Ao), ((A)). where
( : - Ao AoAD
1
2
The scaling in equation 2 ensures that, if do = du. then dp is equivalent to the sin-theta distance 1970; Yu et al. 2015).
Davis & Kahan
2.1
RANDOM PROJECTION ENSEMBLE DIMENSION REDUCTION
procedure for dimension reduction, which is presented in Algo- rithm D We first outline the notation used. Given a dataset of n1 covariate-response pairs ((21, M1 ...., (Zn Bn,)) (IRd R)"1, d-dimensional data-dependent regression method is measurable function 9 R2 x (IR" x [R)"; + R. The function I uses the data (21, (21, (Zru Vns) estimate a regression function g R + R. We often drop all but the first argument of g(: (1), (Z Vr, )) and simply write g(-) for the estimated regression function. We write Id,n, for the set of all such regression methods. In Algorithm [ chosen) projections of subset of the p-dimensional covariate observations æ1,..., Tin-
We now introduce
æ æ Zng are different (randomly
The algorithm also takes as inputs a projection dimension E Ip]: corresponding projection distribution Q Opxd := {A E jæxxd'; diag(A" = (l...., J: number of groups of projections L E N; a group size M E N; a subsample size n1 E [n - 1]: and a base regression method I E Gd,. A full investigation, alongside recommendations on how to choose these inputs in practice is given in Sections and Appendix
Algorithm 1: Handom--rrjjcctin ensemble dimension reduction (RPEDR)
Input: Data ((+1,9) (i)....................... (Itn )n)) (R" x R)", projection dimension d (p). projection distribution Q 2pxa, number of groups L. group size M, sample-split size nı [n 1]. base regression method
9 E Gd,m1-
Let P1.1,..., PL.M L.d. Q.
for f [T.] do
Let Nú he random sample of size ny chosen without replacement from (1,...,n).
for m E [M] do
for i Mf, let he (v:) = ((P) Ti; (P Tm. (i
Let fi,m = z- EiENy {Mr he.m(x;))
Choose the projection that minimizes the mean-square square error within the group, taking the one with the smallest index in the case of ties: Pr.. == Pt.mz where m² == sargminme [M] fic...
Set 1= 1,
E'E P...
10 Calculate the singular value decomposition of π UDUT.
11 Output: The matrix U E _pxP and diag( E IR..
Algorithm starts by generating L- M random projection matrices Pt.m E Rpxd from the distribu- tion Q and partitioning them into disjoint groups of size M. For each projection Pr,ms apply it to all observed covariates, fit the hase regressor Ò on a subsample of the projected data, evaluate the mean-squared error on the complement, and retain the bert-rerforming projection in each group. These selected projections Pt,, are then aggregated by averaging their outer products to obtain II. Finally, the algorithm returns the matrix U and the vector diag(D), provided by the SVD of fī.
The columns of U are the estimated dimension reduction directions, and the singular values in D provide information about the relative importance of each direction. If target dimension dn is
4

---

Published as a conference paper at ICLR 2026
predetermined, then set As (U1....,Uß) estimate of Ag- Otherwise, the singular values provide guidance we propose a data-dependent choice of do in Algorithm in section 3.3
8
Figure 3: Output of Algorithm applied the data used in Figure Inputs: d = 5, 200, M 200, Ti 66, Q: entries rescaled standard Cauchy, J: quadratic least squares. Left: dDag(D). Right: entries of U1 (red filled circles), U2 (blue x), Us (green +), true dimension reduction direction (black filled diamonds).
To demonstrate how our approach works, we applied Algorithm to the data used in Figure 1 In Figure 3 the first two entries of D relatively large, suggesting that the singular vectors U and U2 capture most of the signal. In this case, the true SDR direction Ao is one-dimensional, and nearly all the signal is found in the first singular vector, with |0;1 Aol 0.996. The second singular vector contains only small portion of the signal, with [U, Aol = 0.048. The remaining entries of D are small, suggesting that ..,... Up can be discarded (indeed, max |U] Aul < 0.01).
PRACTICAL CONSIDERATIONS
3.1 RANDOM PROJECTION GENERATING DISTRIBUTION
In this section, we investigate how the choice of distribution for generating random projections af- fects Algorithm [] Projections with independent Gaussian entries are widely used in the literature and have desirable rotational invariance properties (see, for example, Vershynin [018]] Proposi- tion 3.3.2). On the other hand, Cauchy random projections can offer distinet advantages certain applications for example, Li et al. (2007) extended the norm, and Ramiree et al.] (2012) used them for sparse signal reconstruction. As the columns of projection represent directions in the original covariate space, we generate random projections slightly differently from standard practice, Formally, we construct a p x projection matrix P by drawing d i.i.d. unit-norm columns P3 E IR? and concatenating them, i.e., P (Pı,..., Pa). We focus on two choices for the distribution of these columns:
lemma to the f1-
Gaussian columns: We write P, ~ Pf & z/||Z]], where ~ No.
Cauchy columns: We write Pj ~ Qc if P3 d W/|[W/|. where W = (WI,... W,) and Wi...., Wp iid. Cauchy(O, 1).
We conduct a numerical experiment using three different versions of the following model:
Model 1: Let X (0, and 2(A x)² + e, where (0, 1). AT J (1g. Op-g)- for q E [p]. Here 1, (l...., 1) R2 and 0p-q (O...., 0) R--
We set p = 20 and vary q to simulate different sparsity levels: q = (Model la), 9 = 10 (Model 1b) and q = 20 (Model Ic). Three different projection distributions are considered, namely 0d (Gaussian), d (Cauchy) and 1 + 1 (mixture, i.e., with probability we draw P ~ Q and with probability choice affects performance across sparsity levels. Signal strength designed to be the same across settings, and all other inputs to Algorithm [ are held fixed.
draw P 2) The goal is examine how the seaccioonngneratiinn
5

---

Published as a conference paper at ICLR 2026
p-0 (r00
13
05-
08-
11-
06-
17
0.
Figure 4: Boxplots of the sin-theta distance (sce equation 2 between U1 from the output of Algorithm and the true direction Ag J (1,, Op-q) over 100 repeats of the simulation. We present the results for three different instances of Model namely when 4 (left), 4 10 (middle) and 20 (right). We compare the results when the projection distribution Q is Q% (Gaussian), Q8 (Cauchy), and 2 2N + 2 (mixture).
Figure 4 shows that, in sparse models, Cauchy projections outperform Gaussian projections; how- ever, the heavy-tailed Cauchy distribution less suitable in settings without sparsity. Due to the rotational invariance of the problem, Gaussian projections exhibit the same performance regardless of the sparsity level. Using a 50-50 mixture Gaussian and Cauchy projections yields good per- formance in both sparse and dense settings. Therefore, when the sparsity level of Ag is unknown, we recommend using this mixture as the default option. This strategy will be employed in all of the numerical experiments presented later in the paper.
3.2 CHOICE OF L, M AND
Algorithm [ combines the results of evaluating the empirical performance of the base regression method across total of L M d projection directions. More precisely, we aggregate the results of (d-dimensional) projections, where each projection selected as the best from disjoint group of size M. In this section, we explore how varying these parameters affects performance. As we will show, the performance typically improves as L increases, in the sense that the sin-theta distance decreases at rate 1-2 (assuming the other quantities fixed); see Theorem | and our numerical results in Figure 9 Overall, these results lead us to recommend L 200 as the default. The effects of M and d are slightly less straightforward. and their choices are influenced by the ambient dimension p. We recommend M 10p and d min /27 10).
We focus here on the choice of L; details regarding the choices of M and d are provided in Ap- pendix A.2 Theorem [ below compares the sin-theta distance between the first do columns of U from the output of Algorithm run with L groups of projections, with the corresponding output when considering an infinite-simulation version of the algorithm which takes "L oo". More precisely, given data { vi).... (In, yn)} (here treated as fixed pairs in RP x R), d [p]. a projection distribution Q on 2pxd. L E N, M N, nh E [n - 1] and 9 Ga. let denote the first do columns of the output U. We explicitly emphasise the dependence on and assume the true do known. Further define I == E(P,,.TT.) E Spxp to be the expected value of il in line 9 of Algorithm [ the expectation is taken only the ran- domness in the projections and the sample-split in line 4 of the algorithm, but not over the data. Let Ag = (UPo ,..,, U2.), where ",.,,.., U are the singular vectors of 1|
AD := Ao (U....., Ud))
Theorem 1. We have
A04 1lop 2(2m)) 1/2 dp
Edy (S(A)) 2d
11
L1/2
The proof of Theorem is given in Appendix B The second term in the bound suggests that the error of our algorithm decreases at rate L 1/2 as L increases. The first term does not depend on L and can be seen as the infinite-simulation error of our method; it depends on the choices of M and
1/2

---

Published as a conference paper at ICLR 2026
d (as well the base regression method and projection distribution) and would be small when, average, the projections selected are close to Ay. See Appendix A.2 for a numerical demonstration.
3.3 CHOICE OF do
In some scenarios, the user may have a predetermined target dimension, dg in mind. In such cases, one should simply retain the first do columns of U from Algorithm [ and set As = -U1,.... . Ud Indeed, the singular vectors of U can be interpreted as the most frequently selected dimensions (in decreasing order of importance) across projection groups in the algorithm. When the desired dimension is unknown, we propose to use the information contained in the singular values D from the output of Algorithm [ The main idea is to compare these singular values with the component- wise median of those obtained by applying the same procedure to random projections with the same distribution, but without selection within each group (i.e., M 1); see Algorithm Aggorihmm@bblow We use the median rather than the mean, as it provides a more stable summary of the no-selection random projection baseline and is less sensitive to outliers, which can occur under heavy-tailed Q.
Algorithm 2: RPEDR Dimension Estimator
Input:
Projection dimension E [p]. projection distribution Q on 2 number of groups and the vector diag(D) (Di,... ... Dp) of singular values from the output of Algorithm as well a random projection resample size 17 E N.
Let P P........ PL.a
1.1.
Q.
for E [ do
Set Set 11
ΣE P.
Calculate the singular value decomposition of f((r) = ( (17))))))) and let
diag(D(r))
Let = / - (o!
2 1/2, S j). (Here maximum over the Output: The dimension do max E [p]: :
empty set is taken to be zero, interpreted as evidence that there no correlation between and V.)
The motivation is that, for particular direction to be selected, there should be evidence that it carries more signal than would be expected at random. For the leading singular vector U1 (with associated value D1). we require evidence that this direction is more likely to be selected by our algorithm than under random chance. For subsequent components, note that since our projections are rescaled to have trace we have Σー 1 D = trace([i)= = d. This induces a natural penalty for selecting more dimensions when earlier singular values in are relatively large.
Medel 1-5, ww00) p-22
Medde 2-3. price
Molin
3-
2-
3-
3-
Ordend Asles
Comparison between the vector singular values from the output of Algorithm I for Models 1a. and 3, and the median of corresponding (r) in Algorithm 2 for p 20 (left) and 100 (right). We present only the first 10 singular values in each case.
Figure 5:
Figure iillssrrtte the procedure for selecting do in Models 1a, 2 and 3. For Model la, the output suggests retaining the first two singular vectors. This is partly because, seen in Figure the single

---

Published as a conference paper at ICLR 2026
signal 1 (X1 X2) is often split into two. The reason is, in part, that even though the direction v2
va (1, -1,- ................................ ...0.......0 0) is orthogonal to the true signal direction, the response Y is not (marginally) independent of X1 X2. As result, our algorithm suggests selecting (approximately) this direc- tion. We demonstrate in Section 3.4 how applying our algorithm twice in such cases can accurately recover only the true one-dimensional signal. For the other two models, there is clear indication that we should take dg dg i... dy = 1 in Model and dg = 2 in Model 3, and in fact, we recover most of the signal: for Model 2, | esl > 0.999; for Model 3, min [UT eal, | e7 12 0.995.
Further evaluation of dg selection can be found in the full simulation study in Appendix C
3.4 DOUBLE RANDOM PROJECTION ENSEMBLE DIMENSION REDUCTION
In this subsection, we explore how applying our main algorithm second time may help further reduce the dimension in cases where the user is not satisfied with the suggested dimension from an initial application of Algorithms I and 2 This double application strategy particularly effective when we have desired dimension, do say, in mind and the number of covariates contributing to ((A)) exceeds do- In such situations, single application of our algorithm may favour the coordinate axes of the relevant covariates, but fail to combine these effectively into the smallest possible number of dimension reduction directions.
One perspective of this approach, is that the first application of our algorithm acts as a screening step select potentially relevant directions. while the second application then combines these directions to yield a lower-dimes liwer-dimensional projection that still retains a large amount of the signal. The full procedure is presented in Algorithm [3] and a detailed demonstration is provided in Appendix A.5 where we also provide an extension of Theorem | applicabble to this algorithm (see Theorem
Algorithm 3: Double RPEDR
Input: Data ( (ar1, y/), .... (z.. yn)) E [I. R)" and desired dimension du E [p].
Let U and D be the output of Algorithm applied to (21,31),, ... (Zn, Vn)) with L = 200, M = 10p.
d d= minn(pp"/11, 10), n1 = [2n/3], 1 Q = 1/2 e + 1 1Q84 and 9 = ÎMARS.
Let do be the output of Algorithm applied to D, with L = 200, = min( 12,,
Q Q= 10 +
--10 +1000
do > do then
Set Ao == (U,,..., U2,)) and x AJ x. for e [n].
Let ú and D be the output of Algorithm applied to (21,91): v..,,... (zns Un) with L = 200, M = 10do.
d (in([V/d),) 10), 1ı [2n/3], 1084 + Q8" and 9MARS.
Let dn AUr...... ( E Roxdo
else
Let Au (U1, U.o) Roxdo
10 Output: The projection Age
Additional practical considerations including the selection of base regressor, the choice of the split ratio n/nn, and computational notes are presented in Appendix
4 NUMERICAL SIMULATIONS
In this section, we summarize the main findings from a large simulation study comparing the per- formance of our proposal with several existing methods for estimating the central mean subspace S(Ao); the full experimental design and complete results are provided in Appendix C We evaluate three versions of our proposed method (see also Figure 2). The first variant applies Algorithm on its own, and the second is the double dimension reduction approach presented in Algorithm 3 Both variants are suitable when the true projection dimension dg is known. The third variant, which com- bines Algorithm and 2 is designed for situations where do is unknown. For comparison, the default implementations of several existing methods, including SIR (Li 1991), pHd (Li 1992).
MAVE (Xia et al. 2002), DR (Li & Wang 2007, gKDR (Fukumizu & Leng 2014), and drMARS
8

---

Published as a conference paper at ICLR 2026
Liu et al. 2023). The models used in experiments include Model (with 2, introduced in Section 3.1 as well as eight additional models given in Appendix A..1
A.5 and
In the main text focus the slightly simpler scenario where the dimension dg of the central mean subspace S(An) is known. this case, our first proposal, denoted as RPE in Figures 6 and 7 sets is the output of Algorithm with the default inputs
Au = (U1:..., U. where U = (Ul,...,
Up) .',,)..UUp)) recommended. In particular, we set d = min [p 10). I. 200, M 10p,
2n/33]
Q 1º - and 0 .MARS. We also present the performance of Algorithm with do do. which denoted RPE2 in Figures 6 and 7 For the competing methods, we use the R packages dr [Weisberg] 2002) for SIR and pHd, MAVE (Hang & Xial (021) for MAVE, and the relevant code available via GitHub for DR (bttps:
com JSongLab/ HC), gKDR and drMARS (https: / com/ wtom//uuuussttr/r/aaaaa
In each case, we use the corresponding default recommendations to select appropriate tuning parameters.
Model
Model
Model
100
A4
103
ande 17 190
100
Nodel
Model 5
Hodel
30
195
UNISS 15 1912
jide endes 1/6 1h2
195
194
Nodal
Modd 8
Modil
301
245
1955
2014 6NOPS
1912
36
24
005
pile 09005
1955
2011
69006
104
8712
Figure 6: Boxplots of sin-theta distance between the subspaces spanned by the estimated projection Åo and true projection Ao our RPE-based methods the results models 1a-9 over 100 repeats of the experiment, with 200 and p = 50.
competing methods when do known. present
Figures 6 and present boxplots of the sin-theta distance d (So,),,,, ((A)) when 200 and E (50, 500 setting across 100 replicates for the nine models. We observe that both RPE and RPE2 approaches are highly effective. Models la and 4, RPE2 demonstrates clear improvement over RPE, as highlighted in Section 3.4 For the other models, single application of Algorithm performs very well. The drMARS approach is also competitive in Models la, and The other competitors suffer prohibitively due to the curse of dimensionality and are typically ineffective in the 50-dimensional examples shown here; when p == 500, SIR, pHd and DR are intractable because p > n. In the high-dimensional setting, our RPE-based methods show even more pronounced advantages: RPE or RPE2 achieve the best performance across all nine models. There are a few cases where none of the methods achieve sin-theta distance close to zero: In Models 5 and 7, our algorithm tends to find the space spanned by éj and 62, but misses the weak signal in the es direction; Model 6 is particularly challenging, all three covariates X₁ X2 and X3 needed to obtain non-trivial prediction, and the regression problem remains difficult even given the oracle projection. We observe better performance in these models when p 20 and 500- Appendix C

---

Published as a conference paper at ICLR 2026
044
044
92
KFF
06
05
02 02
MARE
BHARS
RPE2
NAVE
BWARS
KFF
0KbK
EPARRS
RPE2
2-
Md0e
gtha:
BHARS
RFF
RP52
SHARS
RPF
RPP2
gkha
BHARS
KFF
RP22
Figure 7: Boxplots of sin-theta distance between the subspaces spanned by the estimated projection Åg and true projection Ag for our RPE-based methods and the three competing methods when da is known. We present the results for models la-9 over 100 repeats of the experiment, with 200 and p 500.
DISCUSSION AND EXTENSIONS
We introduced a new general approach to sufficient dimension reduction, based on aggregating an ensemble of carefully chosen projections. The framework highly flexible, allowing different distributions of random projections and user-specified choice of base regression method. We have provided default recommendations for these aspects and other tuning parameters based on extensive numerical experiments. While our current theory focuses primarily on the effect of L, extending the theory other tuning parameters and developing more automated defaults are promising directions for future work. Our default implementation prioritises statistical performance and can be slower than classical SDR methods. However, the computation is highly parallel across projection groups and candidates and, discussed in Appendix C.3 competitive performance is often achievable with modest runtimes using multi-core CPU parallelism. Therefore, an obvious direction for further improving scalability is to exploit GPU acceleration.
There are several ways in which our framework could extended further. First, a natural assump- tion in high- nigh-dimensional settings is sparsity, specifically the case where the projection Ag is sparse. In this scenario, it may be more effective to consider sparse random projections; perhaps more importantly, employing sparse SVD) (see, for example, Yang al., 2014) to aggregate the chosen projections can encourage sparse estimate Ån- Second, in the context of additive index models (see, for example, |riedman Stuetzle, 1981; Ruan & Yuan 2010), we may seek to sequentially find the signal directions one by one. For the first direction, apply Algorithm and keep only the leading singular vector. To identify the subsequent directions, one could then apply modification of our procedure that always includes the signals found so far and considers random projections that are orthogonal to the previously identified signals. Finally, while this paper has focused on dimension reduction reduction--scccifcalll -pecificcally estimating the projection matrix Ag-it would be interesting to consider al- ternative ways of aggregating the results when the primary objective is to predict the response Y. In particular, improved predictive performance might be achieved by directly aggregating the predic- tions from the chosen projections, opposed to fitting
regression model after projecting the data using the estimated projection Åg. We leave these extensions for future investigation.
10

---

Published as a conference paper at ICLR 2026
REPRODUCIBILITY STATEMENT
Code to implement our algorithms in R is submitted in a .zip file as supplementary material, and also available at https: pettps//githut //gttuub com/ Www.ing /REED The algorithms are fully specified in Algorithm (RPEDR), Algorithm (dimension selection), and Algorithm 3 (double RPEDR) in Sections 2 3.3 and 3.4 respectively. Default settings (base regressor, projection distri- bution, L, M, d. R. /n) and implementation details are summarized in Section and Appendix A The full simulation design (including all nine models, different values of n and p) is reported in Section 4 and Appendix C Real-data preprocessing, splits and evaluation are given in Appendix D Complete proofs of our theoretical results appear Appendix B
ACKNOWLEDGMENTS
We sincerely thank the Area Chair and the three reviewers for their comprehensive and insightful comments, which helped to improve the paper.
REFERENCES
Daniel C. Ahfock, William J. Astle, and Sylvia Richardson. Statistical properties of sketching algorithms. Biometrika, 108(2):283-22 297, 2020.
Laura Anderlucci, Francesca Fortunato, and Angela Montanari. High-dimensional clustering via random projections. Journal of Classification, 39:191-22 2-6,, 2022.
Ella Bingham and Heikki Mannila. Random projection in dimensionality reduction: applications to image and text data. In Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD '01, pp. 245-250. Association for Computing Machinery, 2001.
Timothy I. Cannings, Random projections: Data perturbation for classification problems. WIREs Computational Statistics, 13(1)=e149, 2021.
Timothy I. Camnings and Richard Samworth. Random-projection ensemble classification. Journal of the Royal Statistical Society, Series B (with discussion), 79(4):959-1035, 2017.
R. Dennis Cook. Graphics for regressions with binary response. Journal of the American Statistical Association, 91(435):983-992 1996.
R. Dennis Cook. Fisher Lecture: Dimension Reduction in Regression. Statistical Science, 22(1): 1-26, 2007.
R. Dennis Cook and Bing Li. Dimension reduction for conditional mean regression. The Annals of Statisties, 30(2):455-477 2002.
R. Dennis Cook and Sanford Weisberg. Sliced inverse regression for dimension reduction: Com- ment. Journal of the American Statistical Association, 86(414):328-2) -32, 1991.
Sanjoy Dasgupta and Anupam Gupta. An elementary proof of a theorem of Johnson and Linden- strauss. Random Structures & Algorithmss 22(1):60-655 65, 2003.
Chandler Davis and William Morton Kahan. The rotation of eigenvectors by a perturbation. III.
SIAM Journal on Numerical Analysis, 7(1):1-46, 1970.
Edgar Dobriban and Sifan Liu. Asymptotics for sketching in least squares regression. Advances in Neural Information Processing Systems, 32, 2019.
Xiaoli Z Fern and Carla Brodley, Random projection for high dimensional data clustering: A clus- ter ensemble approach. In Proceedings of the 20th international conference on machine learning (ICML-03). pp. 186-193, 2003.
Jerome H. Friedman. Multivariate Adaptive Regression Splines. The Annals of Statistics, 19(1): 1-67, 1991.
11

---

Published as a conference paper at ICLR 2026
Jerome H. Friedman and Werne Stuetzle. Projection pursuit regression. Journal of the American Statistical Association, 76(376):817-827 1981.
Kenji Fukumizu and Chenlei Leng. Gradient-based kernel dimension reduction for regression. Jour- nal of the American Statistical Association, 109(505):359-37777 2014.
Milana Gataric, Tengyao Wang, and Richard J. Samworth. Sparse principal component analysis via axis-aligned random projections. Journal of the Royal Statistical Society Series B: Statistical Methodology, 82()):22-335, 2020.
Benyamin Ghojogh, Ali Ghodsi, Fakhri Karray, and Mark Crowley. Sufficient dimension reduction for high-dimensional regression and low-dimensional embedding: Tutorial and survey. arXiv preprint 2tt.:::21100092220 2021.
Kam Hamidieh. A data-driven statistical model for predicting the critical temperature of super- conductor. Computational Materials Science, 154:346-320 -344, 2018.
Weigiang Hang and Yingeun Xia. MAVE: Methods for Dimension Reduction, 2021. URL https: jetoj oxq/packkag==MVE package version 1.3.11.
Roger Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 1985.
Hsin Ssiin- Huang, Feng Yu, and Teng Zhang. Robust sufficient dimension reduction via a- distance covariance. Journal of Nonparametric Statistics, pp. 1-16, 2024a.
Zhanliang Huang, Ata Kabán, and Henry Reeve. Efficient learning with projected histograms. Data Mining and Knowledge Discovery, 38(6):3948-4000, 2024b.
William Johnson and Joram Lindenstrauss. Extensions of Lipschitz maps into a Hilbert space. Contemporary Mathematics, 26:189-206, 1984.
Bing Li and Shaoli Wang. On directional regression for dimension reduction. Journal of the Amer- ican Statistical Association, 102(479):997-1008 2007.
Ker-Chau Li. Sliced inverse regression for dimension reduction. Journal of the American Statistical Association, 86(414):316 327, 1991.
Ker-Chau Li. On principal Hessian directions for data visualization and dimension reduction: An- other application of Stein's Lemma. Journal of the American Statistical Association, 87(420): 1025-1039, 1992.
Lexin Li and Christopher J. Nachtsheim. Sparse sliced inverse regression. Technometrics, 48(4): 503-510, 2006.
Ping Li, Trevor J. Hastie, and Kenneth W. Church. Nonlinear estimators and tail bounds for dimen- sion reduction in 11 using Cauchy random projections. Journal of Machine Learning Research, (83):2497-2532. 2007.
Bingyuan Liu and Lingzhou Xue. Sparse kemel sufficient dimension reduction. Journal of Non- parametric Statistics, 0(0): 24, 2024.
Yu Liu, Degui Li, and Yingeun Xia. Dimension reduction and MARS. Journal of Machine Learning Research, 24(309):1-33 2023.
Miles Lopes, Laurent Jacob, and Martin J Wainwright. A more powerful two-sample test in high dimensions using random projection. In Shawe-Taylor, R. Zemel, Bartlett, Pereira, and K.Q. Weinberger (eds.), Advances in Neural Information Processing Systems, volume 24. Curran Associates, Inc., 2.11.
Yanyuan Ma and Liping Zhu. A review on dimension reduction. International Statistical Review, 81(1):134-150, 2013.
Thomas L. Marzetta, Gabriel H. Tucci, and Steven H. Simon. A random matrix-theoreti approach to handling singular covariance estimates. IEEE Transactions on Information Theory, 57(9): 6256-6271, 2011.
12

---

Published as a conference paper at ICLR 2026
Stephen Milborrow. earth: Multivariate Adaptive Regression Splines, 2024. URL https://
R wwt.://epooeecooggggoooooo wwtpppacaaeeearttt....cckkcaccc. R package version 5.3.3.
Mohammad Hossein Rafici and Hojjat Adeli. A novel machine learning model for estimation of sale prices of real estate units. Journal of Construction Engineering and Management, 142(2): 04015066, 2016.
Ana B. Ramirez, Gonzalo R. Arce, Daniel Otero, Jose-Luis Paredes, and Brian M. Sadler. Recon- struction of sparse signals from & dimenniioalliy-reduuee Cauchy random projections. IEEE Transactions on Signal Processing, 60(11):5725-5777 5737, 2012.
Michael Redmond and Alok Baveja. A data-driven software tool for enabling cooperative infor- mation sharing among police departments. European Journal of Operational Research, 141(3): 660-678, 2002.
Henry W. J. Reeve, Ata Kabán, and Jakramate Bootkrajang. Heterogeneous sets in dimensionality reduction and ensemble learning. Machine Learning, 113(4):6633-104, 2024.
Lingyan Ruan and Ming Yuan. Dimension reduction and parameter estimation for additive index models. Statistics and its Interface, 3:493-499, 2010.
Martin Slawski. On principal components regression, random projections, and column subsampling. Electronic Journal of Statisties, 12 12:3673-3712, 2018.
Gian-Andrea Thanei, Christina Heinze, and Nicolai Meinshausen. Random projections for large- scale regression. In S. Ejaz Ahmed (ed.). Big and Complex Data Analysis: Methodologies and Applications, pp. 51-68, Cham, 2017. Springer International Publishing.
Joel A. Tropp. Usser-friendly tail bounds for sums of random matrices. Foundations of Computa- tional Mathematics, 12(4))389 434, 2012.
Roman Vershynin. High-Dinn Biggh-Dmmessonal Probability: An Introduction with Applications in Data Sci- ence. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2018.
Martin J. Wainwright. High- Statistics: A Non-asymptotic Viewpoint, Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2019.
Tengyao Wang, Edgar Dobriban, Milana Gataric, and Richard J. Samworth. Sharp-SSI. Selective
high-dimensions axis-aligned random projections for superviee learning. Journal of the American Statistical Association, 120(549):395-400 2025.
Sanford Weisberg. Dimension reduction regression in R. Journal of Statistiial Software, 7(1):1-22, 2002.
Yingcun Xia, Howell Tong, W. K. Li, and Li-Xing Zhu. An adaptive estimation dimension reduction space. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 64 (3):363 410, 2002.
Haozhe Xie, Jie Li, and Hanqing Xue. A survey of dimensionality reduction techniques based on random projection. arXiv preprint arXiv: /00i:/7706.044 0770000471, 2017.
Shuntuo Xu and Zhou Yu. Neural networks perform sufficient dimension reduction. In Proceedings of the Thirty-Ninth AAAI Conference on Artificial Intelligence and Thirty-Seventh Conference on Innovative Applications of Artificial Intelligence and Fifieenth Symposium on Educational Advances in Artificial Intelligence. AAAI Press, 2025.
Dan Yang, Zongming Ma, and Andreas Buja. sparse singular value decomposition method for high-dimensional data. Journal of Computational and Graphical Statistics, 23(4):923-442, 2014.
Yi Yu, Tengyao Wang, and Richard J. Samworth. A useful variant of the Davis-Kahan theorem for statisticians. Biometrika, 102(2):315-322 2015.
Jia Zhang, Runxiong Wu. and Xin Chen. Optimal sparse sliced inverse regression via random projection. Journal of Computational and Graphical Statistics, pp. 1-14, 2025.
13

---

Published as a conference paper at ICLR 2026
Fang Zhou, Claire, and Ross D. King. Predicting the geographical origin of music. 2014 IEEE International Conference on Data Mining, pp. 1115-1 120, 2014.
A FURTHER PRACTICAL CONSIDERATIONS
A.1 CHOICE BASE REGRESSION METHOD
The choice of base regression method 9 E Id.n, used in Algorithm [ depends on two primary factors: first, the method should be computationally efficient, as will be applied L .M times to different randomly projected datasets; second, the base regression method must capture at least some of the signal after good projection has been applied. We will show that a global polynomial baasedd method is often effective, even when the true signal is not exactly polynomial. However, there are scenarios in which a more flexible nonparametric approach is needed.
We investigate four options in detail. In each case, we describe how Gd.. is defined based dataset (21, yz),..... (2m Ynm) E (RR x R)"" of size m + d(d 3)/2 and test point IR..
1. Global linear least squares (LLS): (2) = à + BT where
(ā, 8) == arg min
(y: a - B 2
(n.f))RRx x
i=1
2. Global quadratic least squares (QLS): gaQ (z) == à bb+b 2T Cz, where
(a. 8 C) =
arg min (e,b,C) xSaxd
(ui a a b)z, 2, Czi)3}
(3)
and Saxd denotes the set of d x d symmetric matrices.
3. Nadaraya-Warte Watson: Fix kernel K and bandwidth 0, then let
ÍNw(2) = Ex= yKK (h 1| =|)
En K(h-1
1 2 =
82 and 0.1.
In experiments below, take K(t)=
/2
4. Multivariate Adaptive Regression Splines (MARS): Let îMARS denote the estimator pro- duced by the forward and backward passes (Algorithms 2 and 3) Friedman (1991), which fits piecewise linear models to the data and automatically selects interactions between vari- ables. We use the implementation of this method from the R-package ea arth (Milborrow 2024). The maximum degree of interaction taken to be 3.
While many alternative base regression methods are available, we focus these four methods for their balance between flexibility and computational feasibility.
To demonstrate how these different options perform in practice, carry out further set of experi- ments using three different models:
Model la: used in Section 3.1 with Ao (12, Op-2)
Model 2: Let X (X1.....
Unif([-1, 1Ip) and Y 2sin(2 + with
c ~ (00 1/4) independent of X.
Model + Let X (X1,
Xp) No and
X6
= 1/2 + (X7 3/2)² +
with (0, 1/4) independent of X. Here Ao (es, E7) is 2-dimensional, where e, denotes the jth canonical Euclidean basis vector in R
ttii sample size requirement is mild and ensures that can apply the global quadratic method below it not satisfied, then either should he reduced, or we are limited using the linear least squares method which only requires m > d. Model 3 used the simulation studies in
1991
Xia al. 2002 Liu al. 2023
14

---

Published as a conference paper at ICLR 2026
In each case, p 20 and 200, and apply Algorithm with d 5, L 200, M 200, 1 2m, and Q _ QN _0 c
PAR
Figure 8: Boxplots of the sin-theta distance (see equation 2 between U1:d0 from the output of Algorithm 1 and the true direction Ag over 100 repeats of the simulation. We present the results Model 1a (left), Model (middle) and Model 3 (right), for p 20 and n = 200. The base method is linear least squares, quadratic least squares, Madarayy-WWaton
MARS.
Figure X exhibits some interesting behaviour. In Model la, the LLS is not effective because no linear signal exists in any d-dimensional projection of the data. For similar reason, the linear and quadratic least squares methods do not perform well in Model 2. However, in Model 1a, the quadratic base method is very effective, as this model is correctly specified after applying the acle projection Ag. MARS also performs well here, although at a higher computational cost. For Model 2, the flexibility of the nouparametric methods required to find the signal, and MARS outperforms Madaraya-Watson The results for Model 3 show that, even when the base model is incorrectly specified, good performance still achievable. Here the true regression function after the oracle projection is not polynomial. The linear base method captures the (approximately) linear signal in X6. but misses the signal in X7. However, the signal is sufficiently well-approximated by quadratic function, allowing us recover much of it using the quadratic base method; the slightly more computationally expensive MARS method performs similarly to the quadratic least squares.
Based on these results. MARS is recommended performs well across all examples presented here and across wide range of settings in our large simulation study (Section 4 and Appendix C
suitable default base regression method. It
A.2 CHOICE OF L., M AND d
To further elucidate the effect L in practice, we conduct a numerical experiment using versions
of Models la. 2 and 3 from the previous subsection with 20 and p 100; see Figure 9 We set n = 200, with the other inputs to Algorithm [ being = 5, . = 10p, nı = 2n /3,
Q = 1 + 10 for the projection distribution, and 9 set to MARS. The number of projection groups L varies from 10 to 1000.
The empirical results in Figure confirm that, in these examples, performance of our method im- proves as increases, line with Theorem I In Model la, observe greater variability across repeats- partly because, in this example, our algorithm captures only part of the true signal within the leading singular vector of U some repeats. As also discussed in Section 3.3 and Appendix C nearly all of the signal is found in the first two singular vectors (which is not fully reflected by the sin-theta distance measure presented here). In fact, the extension Section majority of the signal in a one-dimensional projection using a double application of our procedure.
able find the
We now turn to the choice of M and d, which together determine the total number of projection directions considered within each of the L groups in Algorithm [] In contrast to the choice of L, there is trade-off when choosing and M. For instance, M is taken to small (e.g., less than 20 when p 20), it is unlikely that good projection appears within a block of size M. On the other hand, M is taken to be very large (e.g.. greater than 1000 when = 20), we may start to overfit- although this effect is often negligible due to the data resampling strategy taken in Algorithm I In our experiments, we see that the performance tends to level off rather than deteriorate as M increases. Regarding the projection dimension d, if is too small-especiaa if is smaller than (the
15

---

Published as a conference paper at ICLR 2026
PIO
Wold pnb
Figure 9: The median (solid line) of sin-theta distance (see equation 2 between U1:do from the output of Algorithm| and the true direction Ag over 100 repeats as varies. We present the results for Model la (left), Model 2 (middle) and Model 3 (right), for = 20 (top row) and p 100 (bottom row). The grey shaded region shows area between the 5%-95% quantiles for each value of L.
unknown) do-one or more relevant projection directions will always missed, as the algorithm will consistently choose only more dominant dimension reduction direction. Nevertheless, aim to keep d relatively small for computational considerations and the fact that the base regression methods may suffer from the curse of dimensionality.
To explore these aspects in practice, we repeat the experiments presented Figure 9 but now fix L at the default recommendation of 200, and vary d 21, 2, 10, 20) and M [1000]. The results are shown in Figure IO When p 20, taking d = 1 is not effective across all three models considered, and there clear advantage to taking d do- To elucidate this point, consider an example where do = 1. Setting d = requires one of the M projections in each group to be relatively highly correlated with the true projection direction, which rare unless M very large. With a larger d, it is more likely that we find projection with some signal, which then yields good performance in the final aggregation step. However, d is too large, the base method may begin to overfit-as as seen, for example, Model 1a when p 20 and d = 10.
Regarding M, when p 20 and 5, the performance of our proposal improves rapidly as M increases, and then levels off for M between 100 and 200. When 100 in Model la, M needs to be very large in order to capture of the signal the first singular vector of U: fact, nearly all of the signal lies in the first two singular vectors even for moderately small values of M. For Model 3 with p = 100, the algorithm recovers one the do = directions for small M (for every d), but larger d and M are needed to capture the second direction. This explains the two platcaux seen in the bottom-right plot. Overall, recommend d min min([pp//]] 10} and M 10p as sensible default choices.
A.3 CHOICE OF 71
The purpose this subsection is to present a sensitivity analysis regarding the sample-split ra- tio in Algorithm [] since ni/n trades off selection stability against validation reliability. We use {20, 100} and n = 200, and apply Algo- rithm with d = min{ (2). L= 200, M = 10p, and Q = 1 10 + 10 We vary m1/n 11/2,2/3, 3/4, 4/53. The results presented in Figure 11 where see that the perfor-
Models la, 2 and from Appendix A.l set
16

---

Published as a conference paper at ICLR 2026
thate drannce
Figure 10: The median of the sin tin-thet distance (see equation between U1stp from the output of Algorithm and the true direction Ap over 100 repeats as M varies from 10 to 1000, with = (black), d = 2 (blue), d = 5 (red), d = 10 (green) and d 20 (yellow). We present the results for Model 1a (left), Model 2 (middle) and Model (right), for 20 (top row) and p 100 (bottom row).
mance is very stable varies over this range. Nevertheless, nn ni/n 2/3 throughout the experiments in this paper.
Model 1, 200, 20
Model 2, = 200, p
Model 3, 200, -0
1.0
pedommeeee therr deoveee
b6 0.
0.5
00 0.5
06
0.4
02
0.0
1/2
9,/0
0g/0
Model
Model 100
Model
1.0
metormrcce can SteeN (htreeoo)
50
0.5
to 02 0.4
0.4
U1 0/2
12
23
3/4
4.5
1/2
20
2/4
12
20
01.10
Figure 11: Boxplots of the sin-theta distance between U1:d9 from the output of Algorithm D and the direction Ag over 100 repeats as n2/n varies {1/2, 2/3, 3/4, 4/5). present results for Model la (left), Model 2 (middle) and Model 3 (right), for p = 20 (top row) and = 100 (bottom row).
A.4 COMPUTATIONAL CONSIDERATIONS
The dominant cost of Algorithm I comes from applying the base regression method L x M times; by contrast, projecting the data and computing the final SVD are relatively cheap if the ambient dimension p modest (tens to low hundreds). Notably, these I. x M fits are highly parallelizable,
17

---

Published as a conference paper at ICLR 2026
multi-core execution substantially reduces wall-clock time. In practice, parallelize only the outer loop over the L groups (one group per worker), since these groups are independent of each other; within each worker, the M candidates are evaluated sequentially, and the selected projections are aggregated afterwards.
In Appendix I have shown that global linear least squares (LLS) and global quadratic least squares (QLS) are effective in certain applications. We report below the illustrative runtimes for the setting n 200, p = 50 (Model la). Specifically, we measure the wall-clock time per experiment (using 10 CPU cores) on a 3.20 GHz Intel i9-14900KF computer for MARS (default), QLS, and the results; these numbers are illustrative and complement the full timings in Appendix C.3 (all measured on the same machine).
LLS. Tablee[summann
Base regressor
Time (s)
Notes
MARS (default)
56.133
slow; flexible nonparametric default, robust to misspecification faster; less strong but still often effective
QLS
20.780
LLS
12.794
fastest; effective only when the signal is close to linear
Table I: Runtime (seconds) run for RPE employed with three base regressors on 10 cores; Model la with n 200, p = 50; L=200, M=10p, d= [V). Q = 10 + 200
ad
An appealing direction is to first screen candidate projections within each group using a cheap base regressor (e.g.. QLS), then apply MARS only to the top candidates (for instance, the top 5%) before selecting the best in the group. This two-stage strategy has the potential to preserve good performance at much lower cost; a thorough evaluation is left to future work. Nevertheless, we report runtimes our method with default settings in Appendix
A.5 DEMONSTRATION OF DOUBLE RANDOM PROJECTION ENSEMBLE DIMENSION REDUCTION
We demonstrate the procedure of Algorithm 3 jusing Model la and new model follows:
Model 4: Let X = X,) (N, I) and
Xi X2 X3
=
te
with c ~ NO,, 1/5) independent of X. Here Ao
1 v3 (1,11-1))) is one-
dimensional.
In both Models la and 4, the true projection Ag is one-dimensional and therefore we seek to find the best one-dimensional projection. Figure [2 compares the leading singular vector obtained from
Algorithm D with the output of Algorithm using the recommended default inputs and dg 1. Specifically, we apply Algorithm with L 200, M 10p, min [p' 10). ni = [2n /2n/3], the projection distribution Q N 1000 and ğ = ÔMARS. In both models, second appli- cation yields substantially improved performance. For Model Ia, the first application (Algorithm I typically identifies the first two coordinates as important, but splits the signal across the first two sin- gular vectors. The second application (Algorithm 3 then effectively combines these separate signals into a better estimate of the one-dimensional Ag. similar and even more pronounced improvement is observed in Model 4.
When the initial application Algorithms I and correctly identifies the true dimension dg. ap- plying Algorithm provides no further benefit. However, as long as dg 2 do. there is typically no disadvantage either-see, for example, the results for Model 3 presented in Figure of the main text.
We now provide an extension of Theorem 1, which provides guarantees for the output of Algo- rithm To set the scene for this new result, suppose we input do dg in Algorithm 3, and that Algorithm suggests we should set dg > do, then lines 5-7 of Algorithm apply. Let A Ad xdo be the Q matrix in QR decomposition of AJ Ao- Intuitively, Aă captures exactly the part of
18

---

Published as a conference paper at ICLR 2026
11
Model la
vokee
Model 4
single
souble
Figure 12: Demonstration of the double RPE dimension reduction approach in Algorithm We compare single application of Algorithm| with dg = to double application Algorithm with dg = 1. Left: first 10 entries of the estimated Ag for single run; right: the sin-theta distance over 100 repeats. We report results for E (20, 100}, with Model 1a (top) and Model 4 (hottom).
the true CMS S(Ao) that is preserved inside S(Ao) after the first-stage projection. As in Theo- rem 1, we treat the (/21,81)-. (Zn, }n)} as fixed pairs in R x R. From Algorithm 3, we obtain Ú Und, (U1...., U,,)) and D. Algorithm 3 then outputs Åg Áo(U,..... Ús.) AgU1.do-
Let 11 = EÜÜ
Ao) E dd
Theorem 2 bounds the expected error of the output of Algorithm 3 conditionally on the initial output from the first application of Algorithm The bound includes the false negative error arising from the first application, defined as
dFN Åo An) [ A
dEN measures false negatives and quantifies the amount of the space spanned by Ag that is missed by the projection Åg. Indeed, we have dEN(Ao, Ao) when S(Ao) C S(Ag). This should be in- terpreted as quantifying the amount of the signal missed by the first application of Algorithm I The other terms in the bound in Theorem 2 should be interpreted similarly to those in Theorem I note, however, that since the dimension has been reduced from to do in the first application, the final term in the bound no longer depends on the ambient dimension p.
Theorem 2. We have that
E (S(A., Aol
S de Ao) + /2 A3(A5) Tllpp + 2(2m) d do
12/2
The proof of Theorem 2 is given in Appendix
B
19

---

Published as a conference paper at ICLR 2026
B PROOFS OF OUR TECHNICAL RESULTS
The aim of this section is to prove Theorem 1 and Theorem
Proof of Theorem 4 Recall that the data (11, 81). ..... (Z yn)) are treated as n fixed (i.e. non ran- dom) pairs in R x R. From Algorithm we have that Iī 1 Ei= Pe., P(,,, and Ag E Rexdo denotes the matrix whose columns are the dy leading eigenvectors of πř, ordered according to the corresponding eigenvalues. Recall also that xh = E.Pr..PP. P..) Let II = Ao E RPxP, and for j E [do] write aj for the jth column of Ag. so that Ao (a1, (...... ad.). ). Since Ag Ao = xdo the do largest cigenvalues of II are 1 and the remainder are 0. Indeed, rank (II) = rank (Ao) do. so has do nonzero eigenvalues. Moreover, πa; = 4g a, Ao ( a, = Age, = aj for j E [do]- Then, since IT and IT are symmetric, by the Davis-Kahan Theorem Yu et al. 2015 Theorem 2) we have
dF (((A))) = ||sin (A AO)l - 20/ II
(4)
2d
I Tlo + _1
-[opp))
The second term in equation 4 does not depend L on the randomness in the projections. The remainder of the proof therefore involves controlling the expectation of the first term on the right hand side of equation To that end, we first show that
P(I 11 op 21) P
(5)
for > To see this, for E [L] let 11 = P...P. and write ==
([le- 11.). Then
1 Pr.P 1
Σ (=1
Jr.
Furthermore, B(Jz) = 0, and Jr (and thus also J{) is symmetric. We now show that J & Ixx
in other words d Ipxp - JR is positive semidefinite almost surely, or equivalently that
d² 2
sup U J S (-:|)*]]1)
(6)
First, write Pr.. = (P......... P.,.)). Since diag (PT, Pt..) = (1...., 1)T E Rd, we have 2 = 1 for j [d]. Then, for u R' with |llll2 1, by Cauchy-Schward we have
1)T
v Ih |P 01
PPT
PP. b)
Pe.. 1 || d.
Therefore, u 1l E [0, dj. Moreover, by the linearity of expectation, we also have
u 11 === TE(L))) EEv Ilee)= [0,
whenever |ll||2 = follows that
T (I, IT TT 17 E [-d, dd,].. d),
for all ||0||2 1.
Hence we can write II, - = U,D,UUT where U, E R2xx is orthonormal and D, is diagonal p X matrix with all diagonal entries between -d dl and d (Hom & Johnson. 1985 Corollary 2.5. We deduce that, for all E IRP with ||l||2 = 1, we have
1 L
1 (I, 1 1 SUDDUU
& w11
d² I
J
I I
which establishes equation
Finally, because each of the projections Py,, is chosen from a disjoint group of independently gener- ated random projections, we have that Pt,. for & [7.] are independent, and therefore Jr for & [7.]
20

---

Published as a conference paper at ICLR 2026
also independent. It then follows from the Matrix Hoeffding inequality (see, for example Tropp Theorem 1.3)), that
2012
t²
( I1-
2
S
exp
8|| E
d Exx
llop
which establishes the bound in equation follows:
S
To complete the proof we bound the expectation as
E [[
11
llopp)
P P(|
11 op 2 21)dt
not dt d pv 2T
$ e.e e
VL
as required.
Proof of Theorem
2
First, by triangle inequality we have
dp ((A()))) S de (S(A)U)aa)) S(Ao4)) + (4(55))44))))))))))))))
Now for any B E Rdxxdo
A 1| AA B 43 = AT
=ir[BB |
It follows that
dr (A =
AoU1Ido U1do AD - AAA)(AA)) A lI.
V2
1
|oo(Urdd 10:0 UTddo AD(A5) TA4 (ATT)
V2
|UI.do U..do - Ab(4))
V2
dp (S(
(8)
Further, by an application of Theorem to the data projected by Åo and working conditionally on Åo. we have
B dr (S(r)))(D)))))AAA)))) Ao} s 2 AT(A4) 1po + 2(2m) 1/2 22
d do
L1/2
Now consider the second term on the right hand side of equation 7 For simplicity of exposition,
let P == AAA)(AA) TA Q = AoA Both and Q are symmetric, p2 = P. Q2 = Q. and
tr(PQ) = in). Further
tr(Q) tr(AoA urA Ao) tr(Ia,) do.
and similarly
(r(P) = (r(Ao AT = 1r( ((AJ A1]
= trr((Ab)) An} = vr( = du-
follows that
1 | PP Q12
4((S(A)4))
( ()
2
1((P Q
(10)
2
tr(p2 PQ - QP+Q
2
tr - 2PQ) = do - tP(QQ)
2
21

---

Published as a conference paper at ICLR 2026
Recall that A5 E Adgxd. xdo is given by the QR decomposition of M = Ag Au, Le.
= = A3 Ao = R,
(11)
where R E IRdoxdo is some upper triangular matrix. Then solving equation 11 gives
R = (A5) Ag Ag-
Therefore
tr( = (r4444((45) AJ A.45) = AJ A Ao45) = tr(RRT).
(12)
On the other hand, recall from the definition of false negative measure drx(Au, Ao) we have
( Au)
||(1 AgA
-Ir(A, A0 AA
=Ir(Ar tr(AD (T An AAo}
-rrAAA0 (I Ao
(13)
or(AuAo) A0- AoAD
= do -rr A Ao}
do tr(MM
Finally, we have
tr(M AAARRR (46)
( =t tr RR (A8) (A)) Ir(RRT)
(14)
Combining equation
equation
12
equation
13 and equation 14 we have
de (S()))) deN (Ao.AA.Ap))
(15)
Now taking conditional expectation with respect to Åo on both sides of equation in combination
with equation 9 and equation [5] completes the proof.
C FULL SIMULATION RESULTS
In this section, present the full results of our simulation study, which follows the design of experiments described in Dection4off the main text. The models considered are Model 1 (with q = 2, introduced Section B.Ip. Model and 3 (introduced in Appendix A.II. Model 4 (introduced in Appendix A.5), and five additional models defined as follows: let X (X1,..., X2)) T ~ (0, and (00, 14), with and independent:
Model 5: Y X5+ -01 =E,
Model Y 5X1X2X3 +e
Model 7: Y 4(X, -X2 + Xx) sin (inn (X, + X2)) tc
Model Y (X, + X₂ + 1) +e
Model 9: Y 10 (6.X) + Xa+1 +e
Models 5, 6 and 7 were used in Liu al. 2023). and Model 8 used in Li (1991) ) Xia et al. 2002 Liu et al. (2023).
For each model, consider the settings with {20, 50, 100, 200, 500) and (50, 200, 500). The full description the algorithms considered, along with how tuning parameters were chosen, given in Section 4
22

---

Published as a conference paper at ICLR 2026
C.I DIMENSION do KNOWN
First, we report the complete results for the setting where the dimension do of the central mean subspace S(A)) is known, providing the full picture of this scenario described in the main text. To save space, rather than presenting boxplots here, we focus on the mean and standard error of the sin-theta distance (when do is known) over 100 repeats of the experiments. The results are given in Table 2 (Models 1a-3), Table 3 (Models 4-6), and Table (Models 7-9).
The broad message from the results is that our random projection based algorithms are competitive across all nine models, for different values of n and p. The relative performance is broadly similar to the p E (50, 500). = 200 cases presented in the main text. More specifically, we note the following. First, in the cases where p 2 n, SIR, pHd and DR methods are not applicable. Second, in relatively large p settings (for example, over existing methods. Overall, our approach enjoys the best performance for Models 2, 3, 4. 8, and 9. For the remaining models, our random projection based method and drMARS perform similarly and typically outperform the other competitors.
< 10p), our method shows substantial improvement
Setting
Competing Methods
Our Methods
Model
p
SIR
pHd
MAVE
DR
gKDR
drMARS
RPE
RPE2
50
0.97p.os
0.81o.10
0.870.14
0.97a.06
0.92p.09
0.130.18
0.32g.14
0.12o
20
200
0.97p.os
0.40o.os
0.24g.04
0.330.06
0.29p.os
0.03g.01
0.
0.06p.02
500
0.96p.o7 /
0.246.05
0.10g.02
0.199.03
0.14p.03
0.01o.01
0.13g.0s
0.05p.o2
50
0.93p.08
/
0.976.03
0.34
0.46o.17
0.140 os
200
0.99c.02
0.74o.09
0.930.0s
0.65o.12
0.980.04 04
0.04o.01
0.32n.17
0.06g.02
500
0.98p
0.43o
0.42p.12
0.34
0.57p.20
0.02p.02
0.288.15
0.05g
0.95g.07
0.98p.02
0.33o.2s
0.570.17
0.29o. 20
Ido
1)
100
200
1. 1.00g.om
0.95o
0.970.04
0.99g.01
0.99p.02
0.05g.02
0.46g.17
0.06g.00
500
0.98p.o1
0.720.a7
0.94g.os
0.54a.04
0.99p.oz
0.01o.01
0.38p.18
0,05p.02
200
50
0.97g.06
0.98p.02
0.31o.27
0.64g.19
0.35g
200
0.970.04
0.990.01
0.23o.
0.55e.16
0.06o 02
0.980.05
0.996.01
0.340.31
0.83e.15
0.78a.28
500
200
0.98p.03
0.99p.o1
0.22o.19
0.61o.15
0.08%
30
0,96p os
0.980
0.
0.97a. oa
0.94g.oz os
0.45p.35
0.20p 26
0.26g as
20
200
0.88p.
0.980
0.780.11
0.96g os
0.88p.11
0.09.
0.02o 0o
0.02g
500
0.60p.10 /
0.990.02
0.95% nc 0.61
0.81
0.050.02
0.02g.00 00
0.02g.o0
50
0.91g.07
0.98p.0s
0.
0.210 29
0.43p.41
50
200
0.98o.o2
0.990
0.870.05
0.980.03
0.92p.o6
0.096.01
0.02a.c0
0.02a.om
500
0.92p.o7 /
0.990.01
0. 76.60.
0.980.08 /
0.75p.o7
0.06o.02
0.02..00 00
0.02o.00
50
0.95g.os
0.99p.01
0.77o.51
0.27u
0.54..43
(do = 1)
100
200
0.99p oz
1.000 on
0.91g.05
0.98..0m 02
0.96p.o
0.10p.o1
0.02n no
0.02g
500
0.99
1.00o no
0.84g.05
0.99 02
0.84p.05
0.06o.02
0.02g.00 00
0.02g
50
0.96g.0s
1.00p.oa
0.750.31
0. 0.266.37
0.30o.ss
200
200
0.92p.04
0.98p.02
0.27o.19
0.01o.00
0.01a.00
50
0.97p.0s
1.000.00
0.85o.21
0.33a
0.34g
500
200
0.95g.os
0.99p.01
0.33..25
0.01o.00
0.01o.00
50
1.30g.o7
1.27o.05
1.07o.10
1.29,.0
1.166.11
0.92o.
0.50
0.50
20
200
0. 1.
0.88o.15
0.72g.13
0.990.11
0.89p.12
0.48o.33
0. 0.100.02
0. 0.100.02
500
0.51p.os /
0.350.0s os
0.39p.o7
0.67a.11 /
0.60p.
0.27o.26
0.07o.01
0.07a.
50
1.12p.10
1.32p.05
1.21o.17
0.66p
0.70,.50
50
200
1.276.06
1.290.06
1.070.0s
0.
1.03e.00 os
0.540.27
0.10o.03 03
0.10o 04
500
0.870.10 /
0.930.09
0.80p.10
0.980.00 1
0.790.07
0.32o.26
0,060.01
0.06o.01
50
1.166.11
1.37p.02
1.18p.19
0.71 35
0.80p.31
Ido
2
100
200
1.38
1.380.02
1.090.0s
1.290.04
1.220.04
0.63o.26
0.10o.08
0.110.05
500
1. 1.22p.0s
1.27o.
1.06p.os
1.16g. ne
0.94p.0s
0.36o.27
0.06g.u1
0.07p.01
1.19p.0a
1.39p.o1
1.21o.18
0.83a 32
0.83a
200
200
1.13p.07
1.32p.03
1.11o.19
0.13p 09
0.13g
50
1.220.00
1.41o.00
1.21o.18
0.94o.2 28
0.93o 28
500
200
1.160.07
1.38p.01
1.15o.19
0.25u 22
0.25o 22
Table 2: The average sin-theta distance between Au and Au over 100 repeats for Models 1a, 2 and 3, with {20, 50, 100, 200, 500} and n (50, 200, 500). We compare two versions of our method with six competing approaches. For each setting, we also present 10 times the standard error for each method as subscript. The "Y" entries indicate that the method is not applicable.
23

---

Published as a conference paper at ICLR 2026
Setting
Competing Methods
Our Methods
Model
P
11
SIR
plld
MAVE
DR
gKDR
drMARS
RPE
RPE2
0.900.10
0.93g os
0,51a.12 12
0.850.13
0.65p.13
0,50a 21
0.53a
0.48o.
20
200
0.300..0
0.61g.11
0.29g.05
0.376.07
0.38g
0.28p.
0.34g.14
0.16g
0.170.03
0.38p.08
0.188.03
0.230.04
0.23g
0.140.13
0.29o.12
0.11o.0s
50
0.55o.15
0.86g
0.820.20
0.66g.12
0.56o.10
200
0.67o.12
0.90g os
0. 0.133.08
0. 00.00.77
0.199.os
0.322.10
0.61a.15
0.14u. 04
500
0.30o.03
0.66a 09
0.30g.o3
0.370.04
0.34a
1.180.08
0.53g
0.09o o3
50
0. 616 16
0.95a 04
0. 0.820.20
0.74o.10
0. 0.62o.15
(do
100
0.996..1
0.988 02
0.466.08
0.770.07
0.64g.o4
0. 0.43p.0a
0.70g.10
0.14n.
500
0.52o.o5
0.88g
0.39g.04
0.52o.04
0.41p.03
0.200.10
0.63o.15
0.090
0.690 1o
0.98g.01
0.84o.20
0.80g.09
0.78a.17
200
200
0.488.0s
0.82a
0.85e.18
0.74a.os
0.17u.15
0.770.16
0.996.00 00
0.85
0.84g.os
0.85g.
500
200
0.540.06
0.94a 04
0.80o.15
0.79g 64
0.25o 24
50
1.280.07
1.33p
1.150.00
1.250.00
1.15p
1.060.19
1.06g.02
1.030
1.000.08
1.07p.08
0.95g.12
1.0
0.96a.09
0.730.. 0.730.31
1.00g.00
0.93g.22
500
0.85p.12
0.99g.0s
0.61p.14
0.85p.12
0.7
0.44a
1.00g.00
0.66a.39
1.20p.10
1.30g 05
1.
1.12p.11
1.130.10
50
200
1.240.05
1.38g os
1.150.05
1.180.05
1.09p 05
0.84o 29
1.00g 80
0.97p.1
500
1.040.04
1.11
1.06g.05
1.06..04
0.96g.05
0.500.30
1.00g.00
0.70p.30
1.250.0s
1.36g.0s
1.30o.11
1.15u.13
1.19a.11
(do
2)
100
200
1.400.02
1.40p
1.180.04
1.270.0a
1.20g
1.02a.21
1.00g.00
0.960.18
L.170.03
1.39p.os
1.15g.0s
1.15p.0a
1.05g.04
0.53g 29
1.00a.00
0.74g.34
50
1.29p.07
1.39g
1.290.10
1 1.21n.12
1.21n
200
1.20p.0a
1.29p 03
1.30o.
1 1.00p
1.00g oo
50
0.31
1.41
1.27o.14
1. 1.230.1 18
1.230 18
500
200
1.24p.04
1.36p.01
1.30p.11
1 1.000.00 00
1.000.00
20
1.590.06
1.50p 68
1.
1.590.
1.52p 68
1.100.46
1 14 27
1.20p 26
200
1.58p.07
1.39g 09
1.266.15
1.38p.
1.16g
0.66a
0.73n
0.73p 29
500
1.596.0
1.36g 09
0.41g.05
0.98p.
0.69p.10
0.49p.56
0.36
0.360.19
50
1.470.11
1.64p.03
1.530.18
1. 1.360.23
1.44a.
50
200
1. G8
1,59g
1.55p.0s
1.050.03
1.65p
1.09u.51
1 1.08u.35
1.18g.36
1.676.03
05
1.53p.08
1.496.07
1.47g.
0.66g a5
0.67n.31
0.66n 32
50
1.540.11
1.68g
1.50o.21
1.510.17
1.56g.14
(do- 3)
100
200
1.710.08
1.680 02
1.560.05
1.700.01
1.69p.02 02
1.290 48
1.31o.28
1.41a.21
500
1.706.01
1.64g
1.63p.04
1.670.03
1.69p 02
0.88a as
1.02g 34
1.12g.a5
1.570.10
1.71
1.48o.19
1.52g.19
1.52g.19
200
200
1.59g.0a
1.710.01
1.670.19
1.390.
1.390.28
50
1. 1.644.0a
1.72g 00
1.510.17
1.61g.16
1.61g.
500
1.60p.05
1.72g 00
1.57p.20
1.50p.19
1.50p 19
Table 3: The average sin-theta distance between Åo and Ao over 100 repeats for Models 4-6, with p {20, 50, 100, 200, 500} and E (50, 200, 500). We compare versions of our method with six com- peting approaches. For cach setting, we also present 10 times the standard error for cach method as a subscript. The "Y" entries indicate that the method is not applicable.
24

---

Published as a conference paper at ICLR 2026
Setting
Competing Methods
Our Methods
Model
p
7
SIR
plld
MAVE
DR
gKDR
drMARS
RPE
RPE2
50
1.34p.o5
1.299.c0 06
I. 13
1.330.06
1.30p.ov os
1.090.24
I.130.19
1.15p.18
20
200
1.34p.
1.130.1
1.05p.15
1.07o.08
1.02p.07
0.60g.
0.65p.19
0.610
500
1.33p.06
0.990..1 18
0.54g.22
0.950.00 00
0.780.16
0.550.16
0.55 02
0.540.17
50
1.30g.09
1.30p.os
1.280.18
1.220.
1.20o.17
50
200
1.39p.02
1.31o.0s
1.299.07
1.333.06
1.37p.os
0.670.29
0. 79
0.78p.27
500
1.39p.02 /
1:21o.09
1.200.0a
1. on
1.20p
0.56o.20
0.57g.01
0.53g
50
1.35p.07
1.390.02
1. 1.250.11
1.310.13
1.33p.11
(do
2)
100
200
1.40p.o1
1.39...0 02
1.32g.08
1.400..0 01
1.39p.02
0.87o.a1
0.95g.26
1.01p.26
500
1. 1.40p.0
1.330
1.33g.00
1.30
1.396.00
0.53
0.58g.01
0.56ig.0
1.38g.08
1.40p.01
1.27p.17
1.33g.11
1.33p.14
200
200
1.34g.08
1.40p.o1
1.
1.00u
0.99u
50
1.39p.08
1.41p.00
1.230.18
1.35g.12
1.35p.12
500
200
1.30o.04
1.410.00
1.170.20
1.030 26
1.03g
50
1.30g
1.21o.00
1.02g.11
1.29
1.160.13
0.380
0.21o.15
0.200..1
200
1.03p.15
1.030.07
0.710.12
0.92o.13 13
0.83p.13
0.20o.
0.08o. 02
0.08o.0z
500
0. 68
0.880.
0.38p.07
0.68o.14 1
0.58p.
0.070.05
0.07a.01
0.07a
50
1.11g.11
1.326.o5
0.730.30
0.24n 24
0.27p.24
50
200
1.34g.00
1.200.04
1.020.06
1.16
1.136.00
0.260.16
0.08o
0.08o 01
500
1.06p.11 /
1.77.02
0.82g.05
0.940.10
0.83p.0s
0. 0.140.10
0.070.01
0.07g.01
1.150.10
1.376.02
0.81o.27
0.33o. 31
0.35a 29
Ido
2
100
200
1.40p.o
1.360
1.066.0s
1.320.04
1.280.05
0.360.14
0.07o.0
0.08p.02
500
1.34p.0s
1.190.04
0.99g.0s
1.10a.08
1.06p.or
0.17o.12
0.07o.01
0.07a.0z
50
1.18p.11
/
1.396.01
0.79o.31
0.54a 38
0.54g
200
200
1.08g.06
1.37p.02
0.88
0.10p 02
0. 0.10g 02
50
1.21
1.410.00
0.83o.30
0.740. 36
0.74o a6
500
200
1.1
1.40p.01
0.89o.19
0. 0.188.14
0. 0.180.14
30
1.30p.07
1.31o.00 06
1.21g.09
1.32o.or or
1.22o.07
0.95o.32
1.00p.15
1.026.16
20
200
1.06p
1.17o.0s
1.12g.06
1.18g ee
1.17p.07
0.31o.43
0. 0,18n
0.18g
500
1.00p.04 /
1.08g
1.03p.04
1.06,.00 04
1.08g.05
0.140.31
0.07o 02
0.07,
50
1.250.08
1.330.04
1.240.16
1.07o
1.1
50
200
1.30 1.30p.om
1.336.04
1.22p.04
1.30g.0b
1.18p.04
0.690.43
0.27u.30
0. 0.28u su
500
1.09p.03 /
1:21o.05
1.15p.03
1.18g.06 06
1.15p.0s
0.47o.1e
0.09 06
0.09g o6
50
1.280.07
1.370.02
1.24o.16
1.070
1.130.15
(do-
2)
100
200
1.40p.o3
1.390.02
1.23g.05
1.340.0a
1.270.04
0.890.36
0.35o.88 35
0.37o.36
500
1.266.
1.330
1.22g.03
1.26 04 1
1.15p.02
0.53o.48
0.11o 18
0.1lg.18
50
1.30p.07
1.39g.01
1.220.
1.090.
1.090.15
200
200
1
1.26g.05
1.35p.02
1.250.20
0.340.39
0.34o 39
50
1. 1.32p.07
1.41p.00
1.240.22
1.13u
1.13g 20
500
200
1. 1.28p.05
1.39p.o1
1.220.25
0.28p. 37
0.28n 37
Table 4: The average sin-theta distance between Åo and Ao over 100 repeats for Models 7-9, with p {20, 50, 100, 200, 500} and E (50, 200, 500). We compare two versions of our method with six com- peting approaches. For each setting, we also present 10 times the standard error for cach method as a subscript. The "/" entries indicate that the method is not applicable.
C.2 DIMENSION do UNKNOWN
In this subsection, our primary goal is to evaluate the performance of Algorithm I when the dimen- sion of the final projection, do. is chosen via Algorithm We therefore focus on a single application of our random projection ensemble algorithm. The double dimension reduction technique in Algo- rithm 3 is not considered here, as it requires a prespecified target dimension. For comparison, we also include the results from the competing methods, with their projection dimension selected by the corresponding recommended default approaches-specifiaal)
the marginal dimension test from the dr R package for SIR and pHd: generalised vals-vaaidatiin, which is available as part of the mave package (for MAVE); and the corresponding GitHub links referenced in the previous sub- section (for gKDR and drMARS for these two, we set the maximum projection dimension to be ). The DR method is only applicable with prespecified dg. and it is therefore excluded from this subsection.
[V))).
In the experiments below, reuse the Models 1a-9 described carlier. Our proposal, labelled RPE in the boxplots, follows the recommendations from Section 3 for Algorithm|| with do selected by Algorithm using R = 10000. Competing methods utilize their respective default values for tuning parameters.
25

---

Published as a conference paper at ICLR 2026
Moddd
Model
Model
oche
SMERS
Mive
Model
MASE
SMARS
pd
MiVE
Modal
9001
p76
pton
Figure 13: Boxplots of do from Algorithm 2 alongside competing methods 100 repeats for models 1a-9. present the results for 200 and p = 50. The red horizontal lines indicate true dimension do-
The results in Figure show that, with the exception of SIR and pHd, most methods including ours tend be somewhat liberal choosing dg. in the sense that they typically select do dg. By contrast, SIR and pHd typically select dg = in these examples, indicating no signal is detected as mentioned above, both methods appear to suffer from the curse of dimensionality.
To further assess the performance when do is unknown, we consider two complementary metrics. We introduce the false positives measure
drr(Ag, Ao) [[ - A )
(16)
and we recall the false negatives measure from Appendix A.5
dFN (F,, Ao) = [1 A A
(17)
In both definitions, Åo is a x do projection, where do need not equal do. The metric drp measures false positives, the sense that it quantifies the extent to which Ag contains regions of the ambient p-dimensional space orthogonal to Ag- In particular, if S(A,) S S(pp), then dyr (a,, Ag) = 0. In contrast, deN measures false negatives and quantifies the amount of the space spanned by Ao that is missed by the projection Ag. Indeed, we have d Ao)
when S(Ao) S S(Ao).
Note, however, that small value of dFp or dpN on its own does not guarantee good performance: selecting no directions yields dyp ambient space yields drN but large dpp. In some of our simulations, SIR and pHd effectively behave like the former, returning only a nearly trivial subspace and thus having small dyp but very large drx. Broadly speaking, increasing do typically decreases dex at the expense of increasing dyp. Only when dp = dg can both dep and deN simultaneously (or be close to) zero. Finally, we also have the identity
but large deN, while selecting the full p-dimensional
2 (Ao, An) ( Ao) dex Ao)-
26

---

Published as a conference paper at ICLR 2026
Model 2
Model
Model
5.1
jtd
HAVE
09A15 RIE
5t1
MAVz
056115
54.
ptd
MAVE
GRA15
Model 4
Moded
Model
22
5.
(H4
HOVE
yCOR
0M6R5 RIE
5k
(44
HOVE
09665
KPE
56
(H]
MAVE
APE
Model
Model
23
5h1
MoVe
Rife
5/1
Make
prdn
05615
5h.
More
gr0x
Figure 14: Boxplots of (Ap,Ao, Ao) for the different methods when dy is unknown, for models 1a-9 over 100 simulations, with n = 200. 50. Our RPE method uses Algorithm with do chosen via Algorithm 2 For competing methods, is chosen by the corresponding approaches described at the beginning of this subsection.
27

---

Published as a conference paper at ICLR 2026
Model 1
Model
Model
22
5.1
jfd
HAVE
S5OR,
09A15 RIEE
5/t
MAVE
09611
5h.
MAVE
GRA15
Modil 4
Model
Model
00
05
5.6
(H4
HAVE
yC04
0M6R5 RIE
5k
(41
HOVE
6PF
5k
(M]
MAVE
KPE
Model
Model
8-
5h
MoVE
p508. 086615
RifE
5/
Make
Fife
5h.
More
grbx
Figure 15: Boxplots of drN (Ao, Ao) for the different methods when da is unknown, for models 1a-9 over 100 simulations, with n 200. 50. Our RPE method uses Algorithm with do chosen via Algorithm 2 For competing methods, do is chosen by the corresponding approaches described at the beginning of this subsection.
Figures 14 and 15 present the results of dFp and dEN. respectively, for different methods across the nine models. When do is frequently selected equal dg, namely Models 3, 8 and 9, we observe excellent performance terms of both false positives and false negatives. Models and 4, our algorithm tends to identify the one-dimensional signal direction by selecting the individual coor- dinate axes that contribute to the true projection. these cases, we perform well in terms of the false-negative metric, indicating that nearly all of the signal is captured. However, this comes at the cost of higher false positives, though our approach remains competitive with the other methods considered. For Models and 6, we often take dg larger than do. thereby including some noise directions alongside the signals. Nevertheless, outperforming the other methods considered-eespecially in Model 2, where the competing methods fail to capture the signal directions-although in Model this advantage comes the cost of higher false positives. Finally, in Models 5 and 7, our algorithm tends to select the and e2 directions separately while missing 63, and thus also fails to combine these signals in the optimal way. As result, neither false-positive nor false-negative is close to zero. Nevertheless, our method remains competitive with the other algorithms considered these examples.
approach still excels in terms of false negatives,
Tables 5 and 6 show the full results for the false-positive and false-negative metrics for ad- ditional experiments corresponding to those in Section C.2 We present the results for p {20, 50, 100, 200, 500} with 200. Results for E {50, 500} are broadly similar and are omitted for brevity.
28

---

Published as a conference paper at ICLR 2026
Setting
Competing Methods
Our Methods
Model
p
SIR
pHd
MAVE
gKDR
drMARS
RPE
20
0.06p.23
0.07v.17
0.24g.04
1.480.60
1.150.76
1.000.00
50
0.04
0.00p.00
0.93p.08
2.220.33
1.540.s0
1.00p.04
la (do = 1)
100
p.00p.00
0.00p.00
0.970.04
2.870.26
1.820.84
1.070.21
200
0.970.01
2.840.26
1.750.90
1.000.00
500
0.980.03
2.760.26
1.730.92
1.01o.07
20
0.88p.11 11
0.98 02
1.590.23
1.460.44
1.470.45
0.02o.00
50
0.98o.02
0.99o.01
1.63o.15
2.080.52
1.820.57
0.04o.14
(dg = 1)
100
0.99p.02
1.00p.01
1.62u.14
2.590.11
1.86g.75
0.120.31
200
1.63u.15
2.760.39
2.190.58
0.01g.00
500
1.59p.16
2.700.33
2.070.56
0.01p.00
20
0.466.25
0.04o.1 13
0.92o.30
1.230.52
1.060.63
0.10p.02
50
0.19p.40
0.00p.00
1.340.21
1.940.46
1.500.73
0.190.29
(do 2)
100
0.00p.00
0.00p.00
1.350.20
2.550.41
1.540.84
0.350.47
200
1.400.19
2.620.31
1.860.64
0.120.10
500
1.420.18
2.750.26
1.880.65
0.070.06
20
0.32p.19
0.00u.ov
1.06u.35
1.32p.61
1.090.67
1.420.00
50
0.09
0.00p.00
1.270.21
2. 2.100.49
1.500.81
1.440.09
4 (do 1)
100
0.00p.00
0.00p.00
1.280.20
2.490.45
1.720.83
1.540.18
200
1.280.20
2.540.31
1.770.70
1.400.07
500
1.270.18
2.660.33
1.82o.71
1.360.
20
0.350.21
0.01o. 10
0.980.34
1.450.47
0.970.55
1.0lo.01
50
0.099.28
0.00p.00
1.400.19
2.060.39
1.470.67
1.05p.13
5 (do 2)
100
0.00
0.00p.00
1.41g.18
2.520.40
1.55g.72
1.190.25
200
1.44p.17
2.61o.39
1.88p.65
1.000.04
500
1.46u.18
2.670.33
1.84g.66
0.980.07
20
0.03p.15
0.00p.00
1.260.29
0.850.37
0.590.57
1.01o.16
50
0.00p.00
0.00p.00
1.550.23
1.470.56
1.100.79
1.970.56
(do 3)
100
0.00p.o0
0.00p.00
1.560.22
2.550.50
1.360.73
2.470.62
200
1.6lo.22
2.870.30
1.590.67
0.64g.37
500
1.52o.19
2.91o.31
1.580.
0.650.37
20
0.02p.13
0.00p.00
1.31g.32
1.00p.53
0.870.57
0.990.44
50
0.02p.13
0.00p.00
1.61o.17
1.91
1.230.
1.63p.65
7 (do 2)
100
0.00p.00
0.00p.o0
1.630.16
2.560.41
1.52o.82
2.26p.79
200
1.630.15
2.60g.39
1.640.68
0.640.27
500
1.62o.14
2.75o.33
1.900.60
0.61o.26
20
0.37 0.34
0.11o.17
0.480.24
1.560.40
0.990.63
0.08p.02
50
0.120.35
0.00p.00
0.82o.25
2.230.35
1.42o.71
0.09p.09
(do = 2)
100
0.00
0.00p.00
0.90
2.65g.s5
1. 1.62o.80
0.180.31
200
0.97u.22
2. 7.
1.830.61
0.10p.02
500
1.04g.20
2.73p.23
1.890.65
0.11g.06
20
0.29
0.00p.00
1.450.21
1.63p.10
0.860.61
0.11o.14
50
0.10p.29
0.00p.00
1.530.16
2.230.38
1.260.s0
0.280.36
(do 2)
100
0.00p.00
0.00p.00
1.52g.16
2.670.39
1.740.81
0.41o.51
200
1.53a.17
2.520.45
1.90p.68
0.060.04
500
1.5lo.18
2.51o.39
2.00o.67
0.04o.02
Table 5: The average of (Spn)) S(A)) over 100 repeats of the experiment for Models 1a-9, with E [20, 50, 100, 200, 500) and v. 200. For each setting, we also present 10 times the standard error for each method as subscript. The "f" entries indicate that the method is not applicable.
29

---

Published as a conference paper at ICLR 2026
Setting
Competing Methods
Our Methods
Model
p
SIR
pHd
MAVE
gKDR
drMARS
RPE
20
0.99p.04
0.88p.23
0.24g.04
0.25p.06
0.02o.02
0.05u.01
50
1.00p.00
1.00p.00
0.93p.08
0.69o.15
0.02o.o1
0.04g.01
la (do = 1)
100
1.00p.00
1.00p.00
0.970.04
0.830.07
0.03o.02
0.04p.00
200
0.970.01
0.940.03
0.120.15
0.050.01
500
0.98p.03
0.980.01
0.12o.10
0.06p.01
20
0.88p.11
0.98p.or 02
0.71o.10
0.81o.12
0.070.04
0.02o.00
50
0.98o.02
0.99o.01
0.85o.o6
0.850.06
0.080.04
0.02o.00
(do = 1)
100
0.99p.02
1.00p.01
0.89g.0s
0.91o.04
0.090.01
0.02g.00
200
0.90p.os
0.94p.02
0.14g.0s
0.01g.00
500
0.93p.04
0.98p.o1
0.150.14
0.01p.00
20
1.070.16
1.40p.or
0.69o.12
0.88
0.430.31
0.10p.02
50
1.40p.05
1.41p.00
1.020.09
0.960.08
0.480.28
0.10p.03
(do 2)
100
1.40
1.41p.00
1.05g.09
1.120.06
0.56g.30
0.10p.09
200
1.080.os
1. 1.260.03
0.930.31
0.150.18
500
1.130.07
1.34o.01
0.970.30
0.640.43
20
0.35p.19
1.00p.u0
0.27u.0s
0.35p.u7
0.18p.08
0.1lo.03
50
0.96p.
1.00p.00
0.40p.o5
0.440.05
0.20 0.260.10
0.09p.02
4 (do 1)
100
1.00p.00
1.00p.00
0.43p.06
0.58o.05
0.32p.12
0.09p.02
200
0.460.06
0.750.03
0.71o.22
0.11o.12
500
0.5Lo.06
0.89o.02
0.680.25
0.190.21
20
1.1lo.13
1.41o.00
0.94o.14 14
0.88o.14
0,550.40
0.93o.23
50
1.41o.04
1.41o.00
1.130.06
1.020.06
0.6lo.34
0.97o.16
5 (do 2)
100
1.41p.00
1.41p.00
1.16g.0s
1.15p.04
0.89
0.95g.17
200
1.180.06
1. 1.25p.u2
1.16p.29
1.01g.03
500
1 1.220.05
1.34p.01
1.150.21
1.020.05
20
1.730.02
1.73p.00
1.260.11
1.370.22
1.040.56
0.630.25
50
1.710.00
1.730.00
1.550.07
1.630.11
1.220.18
0.93g.35
(do 3)
100
1.73p.00
1.73p.00
1.57o.07
1.630.05
1.280.50
1. 1.150.28
200
1.58a.0s
1.670.02
1.560.24
1. 1.55g.17
500
1.61o.07
1.70p.01
1.580.20
1.590.15
20
1.41p.00
1.41p.00
0.98u.15
1.01o.11
0.61g.39
0.560.20
50
1.410.01
1.41p.00
1.230.06
1.27o.11
0.61o.40
0.62p.23
7 (do 2)
100
1.41
1.41
1.280.0s
1.32o.05
0.81
0.80n.26
200
1.300.07
1.37o.02
1.080.30
1.070.24
500
1.340.06
1.40o.01
0.980.32
1.120.22
20
1.25p.16
1.31o.16
0.84p.19
0.73
0.18p.18
0.08o.02
50
1.4lo.02
1.41o.00
1.07o.09
0.960.10
0.19
0.08o.o1
(do = 2)
100
1.41p.00
1.41p.00
1.11g.09
1.15g.06
0.25g.10
0.070.01
200
1.11g.10
1.30p.03
0.60p.26
0.10p.02
500
1.140.09
1.370.01
0.63p.27
0.35p.37
20
1. 1.24p.15
1.41o.08
1.08p.07
1.110.10
0.3Lo.
0.260.34
50
1.40p.01
1.41p.00
1.200.05
1.130.01
0.680.44
0.290.33
(do 2)
100
1.41
1.41p.00
1.21o.05
1.220.04
0.850.37
0.34g.36
200
/
1.240.05
1.31o.02
1.140.28
0.450.45
500
1.260.06
1.36p.01
1.080.34
0.440.46
Table 6: The average of (SonSS(A),,SAAA over 100 repeats of the experiment for Models Ia-9, with E [20, 50, 100, 200, 500) and v. 200. For each setting, we also present 10 times the standard error for each method as subscript. The "Y" entries indicate that the method is not applicable.
C.3 COMPUTATIONAL COST
Table 7 reports the wall-clock runtimes in seconds for one run of Model la across different values of n and p. We compare all competing methods considered and our method (label 'RPE') under three compute budgets (single core, 10 cores, 20 cores), all measured on the same 3.20 GHz Intel 19-
30

---

Published as a conference paper at ICLR 2026
14900KF computer. The multi-core columns exhibit that our method achieves roughly 8x speedup at 10 cores and 10-12x at 20 cores.
Setting
Competing Methods
Our Methods
p
16
SIR
plid
MAVE
gKDR
drMARS
RPE
RPE(10 cores)
RPE(20 cores)
50
0.007
0.007
0.051
0.009
0.015
0.138
104.819
14.940
11.677
20
200
0.007
0.005
0.408
0.007
0.078
0.211
135.390
18.826
14,208
500
0.008
0.006
0.961
0.008
0.381
0.276
177.824
24.095
17.888
50
0.047
0.026
0.185
311.972
40.821
29.769
50
200
0.011
0.007
0.827
0.013
0.224
0.178
422.228
54.199
38.174
500
0.014
0.009
2.988
0.020
1.006
0.635
594.113
74.270
50.941
50
0.047
0.109
0.188
702.826
89.105
62.658
100
200
0.032
0.016
0.820
0.018
0.778
0.841
970.345
121.088
82.853
500
0.033
0.020
7.945
0.031
3.327
1.162
1477.162
181.603
116.395
Table Runtime (seconds) for one run of Model la, with p E {200 50, 100) and n E {50, 200, 5003. The "/" entries denote that the method is not applicable.
While the classical SDR methods are faster in wall-clock time, many become inapplicable when p 2 and are typically ineffective in higher-dimensioon regimes. By contrast, random projection- based method delivers substantial improvements across settings and, in several scenarios, is the only approach achieving non-trivial performance (see Tables 2-4, and also the real data analysis in Appendix D). Thus, the additional computation is justified when accuracy matters, and parallel execution keeps runtimes manageable under the recommended defaults.
Modell
Mode12
Mcde13
20100
Nelvvi
0110
03
18 Sussing Time stazond)
B1.ging Time (5222001)
Pedel4
Mldde1
Mice1
01120
0.4
50
Bumning (ssooddl
Bumning Thmm
Buining Time (azzedd
Modkl7
Pidde1
Midel9
123
oefomanne 02 6.4
BAmIiin Timm (5AC0151
Ramning Time taccendo
Rusmning Time 112220051
Figure 16: Performance in terms of sin theta distance versus running time of RPF-based method (20 cores of the same machine) across nine models, with n 500 and p 100. In models 4, and RPE2 is used. In Models and 8, single application RPE with QLS as the hase method used. the other models, a single application of RPE with the MARS default is used.
In cases where computational resources may be limited, we can significantly reduce the computa- tional cost and still remain competitive with (and often superior to) existing approaches. For some
31

---

Published as a conference paper at ICLR 2026
examples, could slightly reduce M and L, use simpler projection distributions, and / or adopt less computationally intensive base method. To demonstrate this, in Figure 16 we present the run- times for some additional experiments with Q = QC (as opposed to Q = 20 + 22 since sparsity in projections is natural in high dimensional settings), M 5p (rather than 10p) and L vary- ing from 10 to 200. We see that, while the performance improves as (and therefore the runtime) increases, we can typically outperform (or come very close to) all existing methods with relatively low computational resources (smaller L); the more computationally intensive default settings used in Table 7 are aimed achieving near-optimal performance. Notably, when using larger values of L in combination with O& in these experiments, RPE-based methods outperform drMARS in Models 5-7, whereas under the default settings in Tables 3 and drMARS was slightly better.
D REAL DATA APPLICATION
In this section, we compare the performance of our methods with competing approaches using real- world data. We make use of the following four datasets from the UC Irvine Machine Learning Repository:
Supereonductivity Hamidiihh 2018): This dataset consists of measurements on 81 phys- ical and atomic covanates for 21 .263 superconductors. The response variable is the critical temperature of the superconductor. We use all p covariates and randomly select a subsample of n = 1000 observations for training, and use the remainder for testing.
2. Communities and Crima RRedmond & Baveja 2002): This dataset contains 1,994 obser- vations 122 predictive variables relating socioeconomic and police data. The response variable of interest is the per capita rate of violent crimes. For simplicity, we remove vari- ables with missing values, leaving p 105 predictive covariates in our experiment, and take a random sample of n = 1329 observations for training (the remainder for testing).
3. Residential Building RRffii & Adeli, 2016); This dataset includes construction cost, sales prices, project variables, and economic variables corresponding real estate single- family residential apartments in Tehran, Iran. There are 372 observations and 105 variables in total. We use all p 105 predictive variables and take random sample of = 248 observations for training (the remainder for testing). There are two options for the response variable-sales price and construction cost-and we include both in our experiments.
4. Geographical Origin of Music Zhou al. 2014 vations on 116 variables relating to audio features of piece of music. We use the latitude of the piece's origin as the response variable. In this case, all p 116 variables are used, and a subsample n = 706 observations taken for training (the remainder for testing).
This dataset consists of 1,059 obser-
In these experiments, since the true projection Ao is unknown, we measure the performance of dif- ferent techniques based on regression prediction accuracy after projecting the data. More precisely, we fixed the projection dimension do to either or 10, and then apply the methods described in Section C.1 using only the training data. The interpretation that each method selects its de most important directions. We then apply MARS the projected training data after each method. The performance is evaluated using the root mean squared error (RMSE) on the corresponding projected test set. The results are presented in Table where our methods exhibit competitive performance across all four datasets considered here.
https: cs uci http::/c.......cccccacaaaaaaaa/a gtt.g/gggrrgapcc coricaltooiinnnn fo
We use cross-validation to determine the degree, with the maximum order of interactions set four.
32

---

Published as a conference paper at ICLR 2026
Dimension reduction technique
Our Methods
Dataset
No proj
do
SIR
pHd
MAVE
DR
gKDR
drMARS
RPE
RPE2
16.5
30.2
19.3
29.9
18.1
18.8
18.5
17.0
17.4
17.3
27.8
18.0
22.2
17.9
20.4
17.0
17.0
3
0.127
0.212
0.140
0.184
0.132
0.138
0.134
0.125
2
0.180
10
0.137
0.187
0.135
0.136
0.132
0.140
0.134
0.126
295
1198
386
158
175
3(sales)
150
10
295
699
295
429
211
125
125
70.4
203
47.6
161
136
80.0
27.6
3ccons.)
32.1
10
68.4
301
42.4
47.1
29.6
24.0
24.0
16.9
18.1
16.7
17.7
17.0
16.8
15.8
25.7
10
17.9
17.5
15.8
17.0
19.4
16.2
15.8
Table The RMSE test four real-data experiments. As a baseline, we also include the result applying MARS directly to the original data (i.e., no projection applied, denoted as "No proj"). For datasets 3 and 4, the DR method yielded very large RMSEs, so we omitted the results (denoted by "P").
DECLARATION OF LLM USAGE
We used large language models solely to aid and polish writing (grammar, clarity, and phrasing) and code presentation (readability and consistency), without altering the logic. LLMs were not used for research ideation, retrieval, or discovery. The authors take full responsibility for the content of this paper.
33


Reviews:
[{"review_id": "c7ZjPEjoyc", "reviewer": "Reviewer_8gwY", "summary": "The paper introduces Random Projection Ensemble Dimension Reduction (RPEDR). The method projects data onto multiple low-dimensional random subspaces, evaluates each projection within its group using a base regression model, selects the best projection in each group, and aggregates the winners via SVD. This procedure adaptively identifies informative directions, yielding stable estimates at modest computational cost. The authors also study design choices, projection distributions and the number of projections, and provide theory showing that the expected estimation error decreases monotonically with the number of projection groups.", "questions": "1. In Algorithm 2 (line 6), how is the threshold $T_{\\ell} > 1/2$ chosen? Could alternative thresholds (other than $1/2$) be used?\n\n2. In Algorithm 3, $L = 200$ is used while many experiments set $n = 200$, which may suggest nontrivial computational cost. Can a smaller $L$ be employed without materially degrading performance? If so, is there guidance or a sensitivity analysis for selecting $L$ as a function of (n,p,d)?\n\n3. Beyond Gaussian and Cauchy projections, are other projection families, e.g., sub-Gaussian or Hadamard-structured transforms, compatible with the framework?\n\n4. I recommend adding a  framework diagram to illustrate the workflow and the relationships among Algorithms 1–3.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* Proposes a novel, flexible, and theoretically principled framework for dimension reduction in high-dimensional regression via random-projection ensembles.\n* Adaptively selects informative projection directions, delivering more stable estimates.\n* Provides theoretical guarantees showing that the expected estimation error decreases monotonically as the number of projection groups increases.", "weaknesses": "I am not deeply familiar with the SDR literature; my expertise lies more in random-projection methods for computational efficiency. Consequently, I cannot fully determine whether the comparisons with existing SDR approaches are sufficiently comprehensive or representative of the current state of the art. The following comments are therefore offered as provisional observations:\n\n* Theorem 1 relies on the assumption $L = \\infty$, which limits its practical value for guiding the selection of $L$ in real applications. Could similar theoretical guarantees be derived for finite $L$?\n\n* The introduction frames the problem as high-dimensional ($n<p$), yet Appendices analyzes parameter choices primarily for $n>p$. This inconsistency weakens empirical support for the stated high-dimensional motivation. Additional experiments or discussion for the $n<p$ regime would improve credibility.\n\n* Algorithm 3 is described as an enhanced version of Algorithm 1, but it is not clear whether Theorem 1, or a corresponding variant, extends to this setting. Clarifying the theoretical validity of the two-stage procedure would strengthen the methodological contribution.\n\n* Algorithm 3 fixes $L = 200$ without justification. Is this choice empirically optimal, or merely heuristic? How sensitive is performance to $L$ across different $n, p, d$? \n\n* While the paper repeatedly emphasizes the flexibility and efficiency of the proposed method, it lacks an explicit characterization of computational complexity, \nThe authors are encouraged to clearly specify the time and memory complexity under typical parameter settings and provide a comparison with established methods like MAVE and drMARS to better illustrate the computational advantages and trade-offs of the proposed approach.\n\n* Although Appendix D mentions real-data experiments, the main text presents only simulation results. It is recommended to include at least one illustrative real-data application n the main paper to enhance the work’s practical relevance and overall impact.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "dXwmbhH6KD", "reviewer": "Reviewer_MBR7", "summary": "The main goal of this paper is to accurately estimate the central mean subspace $S(A_0).$ To achieve that, paper develops random-projection ensemble dimension reduction (RPEDR) framework for high-dimensional regression problems. RPEDR first generates independent random projections and selects the optimal projection based on mean-squared error in each group. After averaging outer products of selected projection matrices to get $\\hat{\\Pi}$, final estimated dimension reduction directions and relative importance of each direction are obtained from SVD of $\\hat{\\Pi}$. The paper then derives the expected estimation error\nbound of proposed method to show the effects of multiple parameters. Finally, extensive empirical results present the effectiveness of new method.", "questions": "As the main assumption in this work, the existence and uniqueness of the CMS, a more detailed description about it is needed, like the exact conditions about this assumption.", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The proposed method is parallelizable, making it applicable to distributed systems.  \n2. The paper not just shows empirical results but also provides theoretical bound for estimation error.", "weaknesses": "1. The superiority of proposed method is not so promising. According to experiment results in Table 2-7, drMARS sometimes shows better performances and needs less runtime. On the other hand, proposed RPE2 and RPE are not always optimal and consume much more runtime. This weakens the contribution of this work.\n2. This framework has multiple parameters, including L, M and d, and each needs special tuning, which makes this framework less practical.\n3. The computation cost of proposed method may bring concerns because of the long runtime.\n4. Considering the main goal of this work is dimension reduction, the reviews about this topic is not clear enough to show existing methods and their relationships with proposed method.\n5. Lack of comparison of theoretical bounds. Since paper derives the error bound, it's necessary to compare with existing bounds.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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": "GLxxAGl17D", "reviewer": "Reviewer_FpME", "summary": "The paper proposes Random-Projection Ensemble Dimension Reduction (RPE), a flexible sufficient dimension reduction (SDR) framework. It samples many low-dimensional random projections, selects the best in each group by held-out regression performance, and aggregates selected projections via the average of outer products followed by SVD; the singular values also drive a data-dependent estimate of the target dimension (Algorithm 2) and a “double pass” refinement (Algorithm 3) when further compression is desired. The authors give a finite-sample bound that decays as L^(-1/2)in the number of groups (Theorem 1), extensive simulations (nine models), and real-data results claiming competitive or superior performance versus SIR, pHd, MAVE, DR, gKDR, and drMARS.", "questions": "1. Related work is missing several recent SDR directions that fit your setting. Such as BC (Huang, H. et.al. (2024)), KSDR (Liu, B. et.al. (2024)) and neural networks based SDR (Xu, S. et.al. (2025)).\n2. The introduction foregrounds the high-dimensional p>n regime, but the main text only highlights a p<n case (e.g., n=200,p=50). Please surface at least one representative p>n setting (e.g., n=50,p=100) in the main paper—moving the results from the appendix if needed. Otherwise, the empirical narrative does not match the paper’s stated scope.\n3. In the proof of Theorem 1, you state that {P_(l,*) }_(l=1)^Lare independent, but marginally they are not: each P_(l,*) is a data-dependent function of the same sample. If the intended claim is conditional independence given the observed dataset D, please state and use it explicitly. Moreover, all concentration steps should be written with P(⋅|D)and E[⋅|D]. \n4. The paper does not report the projection Frobenius distance d_F in unknown-d_0 settings; please include this results. \n5. Please add a sensitivity study varying n_1/n∈{0.5,0.67,0.8}, since n_1 trades off selection variance versus validation reliability. \n\nReference:\n1. Huang, H. H., Yu, F., & Zhang, T. (2024). Robust sufficient dimension reduction via α-distance covariance. Journal of Nonparametric Statistics, 1-16.\n2. Liu, B., & Xue, L. (2024). Sparse kernel sufficient dimension reduction. Journal of Nonparametric Statistics, 1-24.\n3. Xu, S., & Yu, Z. (2025, April). Neural Networks Perform Sufficient Dimension Reduction. In Proceedings of the AAAI Conference on Artificial Intelligence (Vol. 39, No. 20, pp. 21806-21814).", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The method cleanly combines random projections, held-out selection, and SVD-based aggregation, making it plug-and-play with different projection families and base regressors while yielding a stable subspace estimator (with an optional double-pass refinement and a built-in dimension-selection procedure). \n2. The finite-sample analysis isolates an “infinite-simulation” term and a sampling error that shrinks as L^(-1/2), clarifying how performance improves with the number of projection groups.\n3. Extensive simulations spanning multiple p,n settings and diverse baselines, plus real-data evaluations, show competitive or superior accuracy; the algorithms and defaults are presented clearly enough to reproduce results.", "weaknesses": "1. The Contributions section is overly long and diffuses the main claims. The current text mixes method mechanics into the contributions. The Contributions should state only the new ideas.\n2. The paper discusses how to choose L,M, and d, but provides no guidance for n_1 (the within-group training split used for projection selection).\n3. Some arguments in the proof are confusing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "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:

Decision: Accept (Poster)