PaperlayerView exact revision on arXiv

Source-linked AI summary

Sufficient Dimesion Reduction via Generalized Stein's Lemma

Ye Tian

arXiv:2608.15121v1stat.MLcs.LGstat.ME

TL;DR

For multivariate responses, sufficient dimension reduction must capture both predictor relationships with individual outcomes and their joint dependence. This paper uses the generalized Stein’s lemma to estimate the central subspace, and reports consistent outperformance of competitors in simulations and a real-data application.

  • Problem

    Multi-response sufficient dimension reduction must capture both predictor relationships with individual outcomes and joint response dependence, which conditional-mean-based methods may miss.

    Source: §1 · p. 3

  • Method

    The method forms a cross-moment matrix from the multivariate response and predictor marginal score, then uses its singular value decomposition and rank selection to recover the central subspace.

    Source: §1 · p. 4

  • Results

    The method consistently outperforms existing competitors in extensive simulations and a real-data application.

    Source: §1 · p. 4

  • Takeaways & Limitations

    The central subspace is a well-defined target even when the model’s dimension-reduction matrix is not identifiable.

    Source: §1 · p. 3

    SVD-based recovery requires the population matrix M1 to have full row rank, which may fail when too few response components contribute independent information.

    Source: §6 · pp. 18–19

Abstract

from arXiv · show

Sufficient dimension reduction (SDR) seeks the minimal subspace of the predictors that captures the full conditional distribution of the response, which is known as the central subspace (CS). When the response is multivariate, the problem becomes considerably more challenging, particularly when the sample size is limited. Existing methods face different limitations:inverse regression approaches rely on strong distributional assumptions and matrix inversion, and their multi-response extensions suffer from severe slice sparsity; forward regression methods depend on computationally intensive iterative smoothing whose cost grows with the response dimension; and deep learning-based approaches demand large amounts of labeled data. To circumvent these shortcomings, we propose an SDR framework based on the generalized Stein's lemma. Our method constructs a cross-moment matrix between the multivariate response and the marginal score function of the predictors, and recovers the CS via its singular value decomposition. The proposed method does not rely on the linearity condition, avoids matrix inversion as well as iterative smoothing, and can leverage unlabeled data when available. We establish convergence guarantees for the proposed estimator under standard regularity conditions. Moreover, we propose a practical rank-selection algorithm to estimate the dimension of the CS. Extensive simulation studies and a real data application demonstrate that the proposed methods consistently outperform existing approaches across a variety of settings, particularly in moderate-dimensional, label-scarce scenarios with high noise levels.

1. Introduction

The central subspace is the minimal predictor subspace capturing the full conditional distribution; for multivariate responses, recovery must account for joint response dependence, not just individual components.

Introduction

  • Under a multi-response multi-index model with a smooth link and independent zero-mean additive error, the CS equals the CMS, enabling CS recovery through CMS estimation.The CMS captures only the conditional mean and can miss joint response dependence such as conditional covariances or higher-order interactions. Although B is not identifiable, its subspace is invariant and is the target estimand.Source: §1 · p. 3
  • SIR requires linearity and covariance-matrix inversion, while MSIR also suffers severe slice sparsity when sample size is limited relative to response dimension.The linearity condition can fail for skewed or multimodal predictors, and covariance inversion can be unstable or infeasible when n is not substantially larger than p.Source: §1 · p. 3
  • Forward methods such as MAVE and OPG avoid linearity assumptions but require response-wise gradient smoothing; csOPG theory assumes univariate responses and does not account for cross-response covariance.The computational cost of smoothing grows with response dimension q (Xia, 2007).Source: §1 · pp. 3–4
  • DNN-based SDR can capture complex nonlinear dependencies but requires large labeled samples, motivating focus on mildly nonlinear, label-scarce settings.The introduction identifies mutual information-based sufficient representation learning (Zheng et al., 2022) and belted and ensembled neural networks as examples.Source: §1 · p. 4
  • The proposed method forms a cross-moment matrix of the multivariate response and marginal predictor score s(x) = −∇x ln p(x), then uses its SVD and rank selection to recover the CS.The construction avoids matrix inversion and iterative smoothing while relying on weaker assumptions than existing SDR approaches.Source: §1 · p. 4

2. Methodology

The method estimates the central subspace from a cross-moment of the multivariate response and predictor score, using generalized Stein’s lemma to avoid specifying the link function. It estimates the score from data and selects the subspace dimension with an adaptive procedure.

2.2. the Model

  • Multi-response diversity lets a first-order Stein identity recover the central subspace when the link-derivative matrix has full row rank, avoiding second-order score functions.Under the model, the central subspace and central mean subspace coincide with the column space of B; the orthogonal parameterization of B preserves that space.Source: §2.2 · p. 5

2.3.1. Plug-in Estimator

  • The plug-in estimator forms the cross-moment using an estimated predictor score, then estimates B from the top-r left singular vectors of that matrix.This allows score-based estimation when the predictor density—and hence its true score—is unknown.Source: §2.3.1 · p. 5

2.3.2. Score Estimators

  • Kernel score estimators offer convex optimization and small-sample reliability but scale poorly, whereas DNN estimators capture complex nonlinearities with abundant data but lack pointwise error guarantees.The paper uses a Tikhonov-regularized estimator with a curl-free IMQ kernel (TRE) and a DNN estimator.Source: §2.3.2 · p. 6

2.3.3. Score Estimation with Unlabeled Data

  • Because score estimation is unsupervised, unlabeled covariates can augment score training, but training and evaluating on the same labeled covariates creates statistical coupling.The paper considers splitting, supervised-only, and full-data configurations; only Ssplit for DNNs has theoretical guarantees, while coupled DNN configurations are included for empirical comparison.Source: §2.3.3 · pp. 6–7

2.3.4. Rank Selection

  • Since the oracle singular-value threshold is unknown and the noise distribution is nonstandard, adaptive rank selection distinguishes heavy- from light-tailed covariates.For heavy tails it uses cumulative explained variance; for light tails it uses response permutations to test singular values sequentially, retaining at least one direction.Source: §2.3.4 · p. 7

3. Theoretical Analysis

The theory gives finite-sample guarantees for data splitting and asymptotic guarantees for the coupled TRE, then uses cross-moment accuracy to establish subspace estimation and rank-recovery results. It also explains when the guarantees apply and how practical rank selection addresses the unavailable oracle threshold.

Data splitting

  • Data splitting has a finite-sample error bound for arbitrary independent score estimators under sub-Gaussianity and score-error convergence, in both fixed- and growing-dimensional settings.The bound propagates the score estimator’s convergence rate and adds finite-labeled-sample error; a parametric score rate with similarly sized labeled and unlabeled samples yields a parametric cross-moment rate up to logarithmic factors.Source: §3 · pp. 8–9

Coupled TRE

  • For coupled TRE, consistency is established asymptotically under an RKHS source condition in fixed dimensions when the labeled sample grows faster than the square root of the total sample; it fails when β ≤ 1/2.The faster rate occurs only when labeled and total sample sizes are of the same order (β = 1); slower labeled-sample growth makes additional variance terms dominant.Source: §3 · pp. 8, 10

Subspace estimation

  • Cross-moment error bounds yield subspace estimation guarantees when the r-th singular value of M is bounded below, ensuring identifiability of its column space.The subspace discrepancy is measured by orthogonal Procrustes distance, equivalent up to a constant factor to the sine of the principal angles.Source: §3 · pp. 8, 10, 11

Finite-sample behavior

  • Coupled TRE can remain competitive in small samples because RKHS regularization, smoothness adaptation, and the curl-free kernel’s gradient-field prior constrain score estimation and pass to the cross-moment estimator.The cross-moment estimator is a weighted average of score evaluations at labeled points, so these structural benefits matter especially when the labeled sample is small.Source: §3 · p. 11

Rank selection

  • With a suitable oracle threshold and a lower bound on σr(M), hard thresholding recovers the true rank with high probability for both configurations when 1/2 < β ≤ 1.The practical procedure instead uses an energy proportion or permutation-test significance level, selected automatically according to covariate skewness and excess kurtosis, because estimation errors are correlated across entries.Source: §3 · pp. 11–12

4. Simulation Study

The simulation study evaluates central-subspace recovery and rank selection. The proposed score-based methods outperform traditional baselines overall, with performance depending on covariate distribution, score estimator, and use of unlabeled data.

4. Simulation Study

  • The experiments fix covariate and response dimensions at 100 and 10, respectively, and use low signal-to-noise settings; recovery experiments vary labeled sample size while fixing unlabeled sample sizes.The fixed unlabeled sample sizes are 10,000 for NN and 3,500 for TRE.Source: §4 · pp. 12–13, §4.1 · p. 13

4.1. Performance of the Proposed Method in Recovering the CS

  • The proposed methods achieve substantially lower central-subspace estimation distances than SIR, csOPG, KSIR, and MLP across configurations, and score-based methods outperform the baselines overall.SIR is generally the strongest of the traditional baselines. The results are summarized in Table 1.Source: §4.1 · pp. 13–14, Table 1 · p. 15
  • TRE Ssup generally performs best under Gaussian covariate distributions, whereas Sfull TRE can modestly improve recovery under heavy tails when labeled samples are extremely scarce.The Sfull TRE advantage diminishes as the labeled sample size grows; with larger labeled samples, Ssup TRE is comparable to or slightly better than Sfull TRE.Source: §4.1 · p. 13
  • TRE Ssplit consistently underperforms Sfull TRE, and this gap persists as labeled sample size grows because splitting discards labeled data during score estimation.Source: §4.1 · pp. 13–14
  • Score-based methods exploit nonlinear information in the predictor score, available from unlabeled data or TRE structural priors, whereas classical baselines learn response-link nonlinearities from labeled data.This difference helps explain why score-based methods perform well when labeled samples are limited.Source: §4.1 · p. 14
  • For limited-label problems, Ssup TRE is the recommended default except for complex covariate distributions with extremely scarce labels, when Sfull TRE offers a modest improvement; Ssplit NN is the recommended DNN default.Ssplit NN performs comparably to Sfull NN while retaining theoretical guarantees.Source: §4.1 · p. 14

4.2. Performance of the Rank-Selection Procedure

  • For rank selection, proposed methods remain stable across light-tailed settings and robust overall under heavy tails, while traditional methods degrade, including SIR selecting the maximum rank and MLP nearly always failing.csOPG and KSIR are competitive at ranks 1 and 3 under light tails but degrade markedly at rank 5.Source: §4.2 · p. 16
  • Within the proposed rank selectors, TRE generally has lower errors and smaller standard deviations than DNN, except under heavy tails at rank 3, where DNN configurations achieve lower errors.Split and full-data configurations achieve the lowest errors and greatest stability; Ssup NN has failure rates exceeding 80% in certain Gaussian settings.Source: §4.2 · p. 16

5. Real Date Analysis

On the M1 Patch-seq data, TRE variants achieve the strongest predictive performance and rank estimates near 8–9, while unlabeled data add no evident benefit and labeled sample size matters more.

Real Date Analysis

  • The analysis uses 1,213 adult-mouse motor-cortex neurons, 16 electrophysiological responses, and 100 cross-validated genes with normalized, log-transformed, and standardized values.The neurons span excitatory and inhibitory types; gene selection follows the original study’s procedure.Source: §5 · pp. 16–17
  • The evaluation compares supervised, split unlabeled-score, and combined-data configurations across training proportions and labeled fractions, assessing projected features by neural-network test MSE over 100 repetitions.The configurations use labeled data for score estimation and cross-moment computation, unlabeled data for score estimation, or combined labeled and unlabeled data for both.Source: §5 · pp. 17–18
  • SIR and MLP select ranks often around 11–12, while csOPG and KSIR select ranks that are too small; none provides a reliable rank estimate.The study reports that csOPG and KSIR have higher MSE at their selected ranks than at the reference rank.Source: §5 · p. 18
  • TRE configurations select ranks around 8–9 with MSE close to the reference; NN selects only 3–7 directions, with markedly higher MSE than using the reference rank, confirming underestimation.Source: §5 · p. 18
  • Ssplit-NN performs comparably to Sfull-NN and is recommended among DNN variants, while Ssup-NN underperforms with insufficient labeled data; increasing the labeled fraction improves most methods.Increasing the total training proportion from 0.6 to 0.8 yields only small test-MSE differences across methods.Source: §5 · p. 18
  • TRE methods achieve the lowest test MSE overall, with Ssup-TRE lowest and Sfull-TRE close behind; the simulated benefit of unlabeled data does not carry over to this application.The authors suggest that Patch-seq structure or limited sample size may explain the lack of added value from unlabeled samples.Source: §5 · p. 18

6. Limitations and Future Work

The method’s recovery can fail when the population matrix M1 is rank-deficient, and its current framework covers only continuous covariates. Future work could address rank deficiency and extend the method to mixed continuous-discrete predictors.

6. Limitations and Future Work

  • Recovery requires M1 to have full row rank r; too few informative response components or strong correlations can reduce its rank and bias dimension or subspace estimates.Symmetry can make columns of M1 vanish; robust procedures for rank-deficient M1 are a future direction.Source: §6 · pp. 18–19
  • At training ratios 0.6 and 0.8, csOPG and KSIR success rates were 0.74 and 0.84, respectively; csOPG occasionally produced NaN subspace estimates at γ = 0.1.Reported results exclude unsuccessful trials; rank selection remained feasible at γ = 0.1.Source: §6 · p. 18
  • The framework assumes continuous covariates because Stein’s identity requires differentiability of ln p(x), excluding discrete or categorical predictors.Discrete Stein’s methods such as Concrete Score Matching could extend the pipeline to mixed continuous-discrete predictors.Source: §6 · p. 19

Appendix A. Proofs of Main Results

The appendix supplies technical probability and RKHS tools, then uses moment identities and concentration steps in proofs of the paper’s main results.

Appendix A.1. Technical Lemmas

  • Centered products of sub-Gaussian variables are mean-zero sub-exponential, with ψ1 norm bounded by a constant times the product of their ψ2 norms.The stated bound is ∥z∥ψ1 ≤ C∥x∥ψ2∥y∥ψ2 for z=xy−E(xy).Source: Appendix A · p. 22

Appendix A.2. Proof of Lemma 2.1

  • The proof of Lemma 2.1 removes the error term from E[y_j s_i(x)] and stacks the moments to obtain M=BM1.The displayed derivation shows E[y_j s_i(x)] = E[f_j(B^T x)s_i(x)] before stacking across predictor coordinates.Source: Appendix A · p. 22

Appendix A.3. Proof of Theorem 3.3

  • Theorem 3.3’s proof bounds two terms separately, using sub-exponential Bernstein concentration, Jensen’s and Cauchy–Schwarz inequalities, and a union bound.The displayed intermediate probability bounds use 1−δ/2 events before the union bound yields the desired result.Source: Appendix A · p. 23

Appendix A.4. Proof of Theorem 3.7

  • Theorem 3.7’s proof introduces Hermitian positive-definite operator-valued RKHS properties and the kernel ridge regression estimator setup.The appendix also mentions a dataset formed by replacing one sample with an independent copy.Source: Appendix A · pp. 23–24

Appendix A.4.2. Properties of TRE with Curl-Free IMQ Kernels

This appendix analyzes TRE with curl-free IMQ kernels by defining the kernel and establishing regularity and convergence properties, including a high-probability guarantee for sufficiently large samples. It also identifies inconsistency of the semi-supervised estimator when the labeled-sample fraction exponent β is at most 1/2.

Appendix A.4.2. Properties of TRE with Curl-Free IMQ Kernels

  • The curl-free IMQ kernel uses ϕ(r) = (1 + r^2/c^2)^−1/2, whose corresponding scalar kernel is the IMQ kernel kI.The appendix also states that all third-order mixed partial derivatives of kI are globally bounded.Source: Appendix A.4.2 · p. 25

Appendix A.4.3. Proof of Theorem 3.7

  • For β ≤ 1/2, the semi-supervised estimator is not consistent: its error diverges as n^(1/2−β) for β < 1/2 and stays O(1) at β = 1/2.The truncation bias and McDiarmid variance prevent convergence to zero.Source: Appendix A.4.2 · p. 33
  • For 1/2 < β < 1, the truncation bias is absorbed by the McDiarmid variance term, and the dominant variance and bias terms are balanced to obtain a convergence rate.The supplied passage identifies the resulting rate but does not provide its expression.Source: Appendix A.4.2 · p. 33
  • For β = 1, the truncation bias is absorbed by V1, and balancing V1 with B1 gives a rate matching the classical nonparametric score-estimation rate.The comparison is to Theorem 4.2 of Zhou et al. (2020).Source: Appendix A.4.2 · p. 34

Appendix A.4.2. Properties of TRE with Curl-Free IMQ Kernels

  • Choosing ρ=1/4 and applying the bias bound from Zhou et al. (2020) with failure probability δ/2 yields a guarantee with probability at least 1−δ for sufficiently large n.The resulting constant Cδ absorbs constants that depend on δ through concentration-inequality thresholds.Source: Appendix A.4.2 · p. 34

Proof of Lemma Appendix A.2 … Proof of Lemma Appendix A.9

The appendix proofs establish probability-tail bounds and RKHS stability inequalities, then verify boundedness properties of the IMQ kernel and its derivatives used in the analysis.

Proof of Lemma Appendix A.2

Proof of Lemma Appendix A.4

  • The RKHS argument uses the reproducing property and Cauchy–Schwarz to derive the stated function bound from the kernel’s diagonal quadratic form.It rewrites the squared RKHS norm as ⟨f(x), K(x,x)f(x)⟩.Source: Proof of Lemma Appendix A.4 · p. 38

Proof of Lemma Appendix A.5

  • The proof bounds the effect of replacing one training observation on regularized-risk minimizers using the kernel bound and convexity.It defines risks on the full and replaced datasets and combines the minimizers’ optimality inequalities.Source: Proof of Lemma Appendix A.5 · pp. 38–39

Proof of Lemma Appendix A.6

  • A bound on the kernel’s divergence, together with Lemma A.5, bounds the first term, and combining terms yields the stated result.The proof assumes sup_x,y ∥div_x K(x,y)^⊤∥2 ≤ Mζ.Source: Proof of Lemma Appendix A.6 · p. 41

Proof of Lemma Appendix A.7

  • The proof shows the IMQ kernel’s third-order mixed partial derivatives are uniformly bounded by checking their behavior as r approaches zero and infinity.At infinity the derivatives decay, while near zero potentially singular terms cancel, leaving bounded terms.Source: Proof of Lemma Appendix A.7 · p. 42

Proof of Lemma Appendix A.8

Proof of Lemma Appendix A.9

  • The IMQ curl-free kernel’s diagonal limit is c^-2 I_p, independent of x, so its supremum over x has the same value.The proof cites the radial curl-free kernel representation of Zhou et al. (2020).Source: Proof of Lemma Appendix A.9 · p. 44

Proof of Lemma Appendix A.10 … Appendix B.2. TRE (for Ssup

The appendix bounds the score estimator and derives kernel regularity bounds, then records fixed implementation settings and architectures for the neural and kernel-based estimators. It also describes how the MLP baseline produces candidate directions.

Proof of Lemma Appendix A.10

  • The proof bounds the score estimator by comparing its objective with the feasible zero function, then converts an RKHS norm bound to a pointwise Euclidean bound.The pointwise step uses the reproducing property and Cauchy–Schwarz.Source: Proof of Lemma Appendix A.10 · p. 45
  • The proof derives Lipschitz bounds for the kernel and its divergence from bounds on derivatives of the underlying kernel.The kernel bound uses the mean value theorem; the divergence components are linear combinations of third-order derivatives.Source: Proof of Lemma Appendix A.10 · pp. 46–47

Appendix B. Implementation Details of Neural and Kernel-Based Estimators

  • All experiments fix random seeds at 42 plus the repetition index and keep hyperparameters constant across experiments.The appendix summarizes training hyperparameters for the DNN score estimator, TRE, and MLP baseline.Source: Appendix B · p. 48

Appendix B.1. DNN Score Estimator (for Ssup

  • The denoising score-matching network uses a p-to-256 input projection, two residual blocks, and a 256-to-p output projection.The input projection is followed by ReLU and LayerNorm; the residual blocks use LayerNorm, skip connections, and ReLU.Source: Appendix B.1 · p. 48

Appendix B.2. TRE (for Ssup

  • Table B.5 lists training parameters for TRE, but the supplied passage does not provide their values.Source: Table B.5 · p. 48
  • The MLP baseline has two hidden layers of widths 256 and 128, and its rescaled first-layer weights yield candidate directions via right singular vectors.Its output layer maps 128 units to q without an activation.Source: Appendix B.2 · p. 49

Appendix B.3.2. Training Parameters · Appendix B.4. Evaluation MLP (for assessing qualities of estimated CSs) · Appendix B.4.1. Architecture

The appendix reports training hyperparameters for the MLP baseline and describes a single-hidden-layer MLP used only to assess the quality of estimated central subspaces. The evaluator is trained on projected data, early-stopped on validation data, and maps r inputs through 32 hidden units to q outputs.

Appendix B.3.2. Training Parameters

Appendix B.4. Evaluation MLP (for assessing qualities of estimated CSs)

Appendix B.4.1. Architecture

Appendix B.4.2. Training Parameters

Appendix B.4.2 presents a table of training hyperparameters for the evaluation MLP.

Appendix B.4.2. Training Parameters

  • Table B.7 lists the evaluation MLP’s training hyperparameters.The supplied passage identifies the table’s scope but does not provide any hyperparameter values.Source: Table B.7 · p. 50

Appendix C. Remaining details of the Simulation Study and Real Data Analysis

Appendix C describes three rank-selection procedures that adapt to covariate-tail behavior: ARS uses cross-moment singular values, while SIR and MLP use method-specific quantities and tests.

Algorithm 2 Adaptive Rank Selection (ARS)

  • ARS detects tail behavior from skewness or kurtosis, then selects rank by cumulative singular-value energy for heavy tails or response-permutation quantiles for light tails.For light tails, it compares cross-moment singular values with permutation quantiles and retains at least rank one.Source: Appendix C · p. 51

Algorithm 3 Rank Selection for SIR

  • SIR uses cumulative energy of slice-averaged singular values for heavy tails, but for light tails compares them with permutation quantiles after random-direction projection and slicing.Tail behavior is detected using feature skewness or excess kurtosis.Source: Appendix C · p. 52

Algorithm 4 Rank Selection for MLP

  • MLP rank selection uses singular values of trained first-layer weights; heavy tails use an energy threshold, while light tails use conservative energy and a first-direction-dominance rule.The light-tail procedure also handles zero total variance by selecting rank one.Source: Appendix C · p. 53
Loading 2608.15121v1…