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

Multivariate sufficient dimension reduction must capture joint response dependence, which makes central-subspace estimation challenging. This paper uses generalized Stein’s lemma and reports consistent outperformance of existing competitors in simulations and a real-data application.

  • Problem

    Multivariate responses require sufficient dimension reduction methods to capture joint dependence beyond conditional means, complicating recovery of the central subspace.

  • Method

    The method forms a cross-moment matrix between the multivariate response and predictor score function, then estimates the central subspace using SVD and rank selection.

  • Results

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

  • Takeaways & Limitations

    The framework provides a practical approach to estimating the central subspace without matrix inversion or iterative smoothing.

  • Takeaways & Limitations

    The method can underestimate dimension or bias subspace estimates when the population cross-moment matrix is rank-deficient.

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 introduction frames multi-response SDR as a challenge because the central subspace captures dependence beyond the conditional mean, while existing inverse, forward, and deep-learning methods face restrictive assumptions, computational costs, or labeling demands. It proposes a generalized-Stein framework that estimates the central subspace from a response–score cross-moment matrix and SVD, with rank selection and weaker assumptions.

  • The central subspace is the most parsimonious predictor subspace capturing the full conditional relationship between response and predictors.
  • Multi-response SDR is difficult because the central mean subspace captures only E(y | x), omitting conditional covariances and higher-order response interactions.
  • Existing methods are limited by distributional assumptions, computational costs that grow with response dimension, weak multi-response theory, or large labeled-data requirements.Inverse regression uses response slicing, forward regression relies on nonparametric smoothing, and DNN approaches require substantial labeled data.
  • The proposed method constructs a cross-moment matrix between the multivariate response and marginal predictor score function, then recovers the central subspace through singular value decomposition.Its dimension is determined by a rank-selection procedure.
  • The framework avoids matrix inversion and iterative smoothing, uses weaker assumptions, and consistently outperforms existing competitors in simulations and a real-data application.

2. Methodology

The methodology models the multivariate response through a minimal central subspace and estimates it using a generalized Stein cross-moment with estimated predictor scores. It also addresses score estimation with labeled and unlabeled data and adaptively selects the central-subspace rank.

  • Model and central subspace: The multi-response multi-index model defines the central subspace as the column space of a minimal-rank matrix B satisfying y ⊥⊥ x | Bᵀx.The model uses predictor x ∈ R^p, response y ∈ R^q, link function F, and independent zero-mean error.
  • Generalized Stein estimation: Generalized Stein’s lemma separates B from link-function derivatives, enabling central-subspace estimation without specifying F.The method orthogonally parameterizes B without changing its column space and relies on stated existence and boundary conditions.
  • Generalized Stein estimation: The multivariate response permits first-order Stein recovery of the central subspace, avoiding the slower second-order approaches used in single-response multi-index settings.Second-order methods require reliable estimation of higher-order scores and converge substantially more slowly.
  • Score estimation and SVD: The estimated cross-moment uses a data-driven score estimator, and the top-r left singular vectors of its estimate provide a plug-in estimator of B.The score need not be known because the predictor density is rarely known in applications.
  • Score estimation and SVD: Kernel and DNN score estimators offer complementary tradeoffs: kernels suit small samples, whereas DNNs capture complex nonlinearities but depend heavily on large samples.The implementation uses TRE, a Tikhonov estimator with a curl-free IMQ kernel, and a DNN estimator.
  • Rank selection: An adaptive rank-selection strategy distinguishes light- and heavy-tailed covariates because the oracle threshold and noise distribution are unknown.Tail behavior is assessed using sample skewness and excess kurtosis for each feature.

3. Theoretical Analysis

The theoretical analysis establishes finite-sample and asymptotic guarantees for cross-moment estimation under data splitting and coupled TRE configurations, and proves high-probability recovery of the central-subspace rank by oracle hard thresholding.

  • Data splitting configuration: Independent score estimation enables a finite-sample cross-moment error bound under mild score-convergence and sub-Gaussianity assumptions.The bound applies to arbitrary score estimators and any sample size.
  • Coupled TRE: For coupled TRE, an explicit RKHS analysis yields asymptotic consistency rather than an equally general finite-sample bound because score estimation and evaluation data are dependent.Theorem 3.7’s faster rate occurs only when nℓ is of the same order as n; otherwise additional variance terms dominate.
  • Rank selection: Oracle hard thresholding of the cross-moment singular values recovers the true CS rank r with high probability under suitable regularity conditions.The practical adaptive procedure is evaluated empirically, while correlated estimation errors invalidate naive entry-wise thresholding.
  • Data splitting configuration: The data-splitting bound propagates the score estimator’s convergence rate, attaining the parametric rate up to logarithmic factors when both sample sizes scale parametrically.Slower nonparametric score convergence is inherited, with an additional nℓ^-1/2 term.

4. Simulation Study

The simulation study evaluates CS recovery and rank selection across competing methods, covariate distributions, and labeled-data regimes. Proposed score-based methods generally outperform traditional baselines, with performance depending on the estimator, data configuration, and tail behavior.

  • Experimental design: The study compares TRE and DNN score estimators, each under Ssplit, Ssup, and Sfull configurations, against multi-response SIR, csOPG, KSIR, and an MLP baseline.Experiments use p = 100 predictors, q = 10 responses, and unlabeled samples for the NN and TRE estimators.
  • CS recovery: TRE generally outperforms DNN, while Ssup_TRE is strongest for simple covariate distributions and Sfull_TRE helps under complex heavy-tailed distributions with extremely scarce labels.Ssplit_TRE consistently degrades because data splitting discards labeled information during score estimation.
  • CS recovery: For DNN estimation, unlabeled data are essential: Ssplit_NN and Sfull_NN substantially outperform labeled-only configurations, while Ssplit_NN is generally comparable to Sfull_NN.The split configuration is emphasized because the coupled DNN configuration lacks theoretical guarantees.
  • Rank selection: In rank selection, proposed methods remain robust overall, whereas traditional methods deteriorate at r = 5 and perform especially poorly under heavy-tailed covariate distributions.SIR consistently selects the maximum allowable rank and MLP nearly always fails in heavy-tailed settings; NN failure rates exceed 80% in certain Gaussian settings.

5. Real Date Analysis

The real-data analysis applies the proposed framework to Patch-seq data with standardized gene-expression and electrophysiological features under supervised and semi-supervised configurations. TRE configurations select parsimonious ranks and achieve the strongest predictive performance, while unlabeled data provide no clear advantage in this application.

  • Data and preprocessing: The Patch-seq dataset contains 1,213 mouse motor-cortex neurons, with 16 electrophysiological responses and 1,000 selected predictor genes.Cross-validation reduces the predictors to 100 genes; counts undergo sequencing-depth normalization, ln2(x + 1) transformation, and standardization.
  • Experimental design: The evaluation varies training proportions and labeled fractions across supervised, split, and full configurations using independent fitting and hold-out test sets.The labeled fraction γ ranges from 0.1 to 0.5, while the fitting and test proportions are fixed at 0.15 and 0.05.
  • Rank selection: SIR and MLP select ranks around 11 or 12, whereas csOPG and KSIR select ranks that are too small for the nonlinear setting.The results interpret these rank patterns as overfitting or failure to identify a parsimonious representation.
  • Rank selection: TRE configurations select ranks around 8 to 9 with MSE close to the reference, while Ssup-NN selects only 3 to 7 directions and has markedly higher MSE.Ssplit-NN performs substantially better than Ssup-NN and is comparable to Sfull-NN.
  • Predictive performance: TRE methods achieve the lowest test MSE, with Ssup-TRE slightly ahead of Sfull-TRE; unlabeled data do not improve performance in this real-data application.The study attributes the lack of benefit possibly to Patch-seq structure or limited sample size, and reports small differences between training proportions 0.6 and 0.8.

6. Limitations and Future Work

The method’s main limitations are its dependence on a full-row-rank population matrix and continuous covariates. Future work should address rank deficiency and extend the framework to mixed continuous-discrete predictors, while implementation issues remain for some competing methods.

  • Rank and identifiability: SVD-based CS recovery requires the population matrix M1 ∈Rr×q to have full row rank r.At least r response components must provide linearly independent gradient directions.
  • Computational robustness: At γ = 0.1, csOPG and KSIR sometimes produced NaN subspace estimates; success rates were 0.74 at training ratio 0.6 and 0.84 at 0.8.The reported results include only successful trials.
  • Rank and identifiability: Symmetric response components or strongly correlated response functions can reduce M1’s effective rank, causing underestimated CS dimension or biased subspace estimates.If response components are symmetric about zero and B⊤x has a symmetric distribution, corresponding columns of M1 vanish.
  • Covariate types: The framework assumes continuous covariates because Stein’s identity requires differentiability of ln p(x), excluding discrete and categorical predictors.Discrete Stein methods such as Concrete Score Matching could extend the pipeline to mixed continuous-discrete predictors.

Appendix A. Proofs of Main Results

Appendix A establishes technical lemmas and proves the paper’s main theoretical results. The arguments use sub-exponential concentration, generalized Stein identities, union bounds, and reproducing-kernel preliminaries.

  • Appendix A.1. Technical Lemmas: Products of sub-Gaussian variables, after centering, are mean-zero sub-exponential with controlled ψ1 norm.Specifically, ∥z∥ψ1 ≤ C∥x∥ψ2∥y∥ψ2 for an absolute constant C.
  • Appendix A.3. Proof of Theorem 3.3: Theorem 3.3 is proved by conditioning, applying the triangle inequality, and controlling separate terms with sub-exponential Bernstein bounds and Jensen’s inequality.The argument obtains probability-at-least-1 − δ/2 bounds for intermediate terms and combines them using a union bound.
  • Appendix A.4. Proof of Theorem 3.7: The proof of Theorem 3.7 begins by defining Hermitian positive-definite reproducing operator-valued kernels and introducing vector-valued RKHS lemmas.The appendix also defines the kernel ridge regression estimator under a bounded-response condition and considers replacing one sample with an independent copy.

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

This appendix defines the curl-free IMQ kernel and establishes boundedness properties for its associated TRE estimator. It also derives convergence behavior showing inconsistency when the labeled-sample fraction is at most one-half and rate recovery in denser-label regimes.

  • Kernel properties: All third-order mixed partial derivatives of the scalar IMQ kernel are globally bounded.This property is stated in Lemma Appendix A.7.
  • Kernel construction: The curl-free IMQ kernel is obtained from the radial function ϕ(r) = (1 + r2/c2)−1/2 and is denoted KCI.Its scalar counterpart is the IMQ kernel kI(x, y).
  • Estimator properties: The operator norm of KCI(x, x) is uniformly bounded, and the curl-free IMQ TRE estimator is pointwise bounded for λ > 0.The pointwise bound applies to any test point when the estimator is trained on a fixed dataset.
  • Convergence regimes: For β ≤ 1/2, the semi-supervised estimator is not consistent because truncation bias and McDiarmid variance prevent the error from converging to zero.When β < 1/2, the error diverges as n1/2−β; when β = 1/2, it remains O(1).
  • Convergence regimes: For 1/2 < β < 1 and β = 1, balancing the dominant variance and bias terms yields convergence rates, with the β = 1 rate matching the classical nonparametric score-estimation rate.The β = 1 case is explicitly identified as matching Theorem 4.2 of Zhou et al. (2020).

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

The appendix proves auxiliary probability, RKHS, stability, kernel-derivative, and matrix bounds needed for the proposed method. These arguments rely on union bounds, reproducing properties, convexity, asymptotic expansions, and continuity.

  • Proof of Lemma Appendix A.2: Union bounds control predictor norms and complete the probability argument for Lemma Appendix A.2.The proof applies coordinatewise tail control, then extends it over observations and analogously handles the response bound.
  • Proof of Lemma Appendix A.4: The reproducing property and Cauchy–Schwarz inequality establish the RKHS pointwise bound in Lemma Appendix A.4.The argument uses the identity ⟨K(·, x)f(x), K(·, x)f(x)⟩HK = ⟨f(x), K(x, x)f(x)⟩Rp.
  • Proof of Lemma Appendix A.5: Convexity, regularized-risk optimality, and the kernel bound control the difference between minimizers on full and replaced datasets.The proof compares empirical and regularized risks, separates the replaced observation, and uses nonnegative terms for k ≠ j.
  • Proof of Lemma Appendix A.6: The bound in Lemma Appendix A.6 follows by bounding the first term with the derivative supremum and combining it with the second term.The derivative condition is supx,y∈X ∥divxK(x, y)⊤∥2 ≤Mζ.
  • Proof of Lemma Appendix A.7: For the IMQ radial function, third-order partial derivatives decay at infinity, remain finite near zero, and are therefore uniformly bounded.At infinity, the derivative orders are O(r−2), O(r−3), and O(r−4); near zero, potentially singular terms cancel.
  • Proof of Lemma Appendix A.8: For p ≥3, the IMQ-kernel expression is non-negative, while its auxiliary function is bounded and decays as O(u−4).These upper and lower bounds establish supu≥0 Fp(u) ≍p.
  • Proof of Lemma Appendix A.9: The curl-free IMQ kernel has an explicit matrix form whose limit as y →x is (1/c2)Ip, independent of x.The rank-one term vanishes because ∥rr⊤∥op = r2 →0, while the second term converges to (1/c2)Ip.

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

The appendix proves Lemma A.10 by bounding the kernel estimator through RKHS arguments and specifies reproducible implementation details for the DNN, TRE, and MLP estimators. It also describes how TRE extracts candidate directions from rescaled first-layer weights.

  • Proof of Lemma Appendix A.10: The proof decomposes the score estimator and bounds its RKHS norm using feasibility of the zero function and the resulting objective comparison.The argument then uses the reproducing property, Cauchy–Schwarz, Lemma Appendix A.9, and the triangle inequality to obtain pointwise control.
  • Proof of Lemma Appendix A.10: The proof represents the estimator through normal equations involving a block kernel matrix K and vector h.The coefficient vector c satisfies the normal equations, while K is assembled from vector-valued kernel components.
  • Proof of Lemma Appendix A.10: The proof derives Lipschitz bounds for the kernel and its divergence from bounds on third-order derivatives of kCI and the mean value theorem.For the kernel gradient, each component is bounded by M3, and the resulting bound includes 2M3MCIζ/λ^2.
  • Appendix B. Implementation Details of Neural and Kernel-Based Estimators: Training hyperparameters are summarized for the DNN score estimator, TRE, and MLP baseline, with seeds fixed at 42 + repetition index and settings constant across experiments.These choices are stated to support reproducibility.
  • Appendix B.1. DNN Score Estimator (for Ssup: The denoising score-matching network projects p inputs to 256 units, uses two residual blocks, and projects back from 256 to p without output activation.The input projection applies ReLU and LayerNorm; each residual block includes LayerNorm, a skip connection, and ReLU.
  • Appendix B.2. TRE (for Ssup: TRE training parameters are provided in Table B.5.The supplied passage identifies the table as the parameter summary for TRE.
  • Appendix B.2. TRE (for Ssup: The MLP baseline has two hidden layers with widths 256 and 128, followed by an output layer mapping 128 to q without activation.Its layers are fc1: p →256, fc2: 256 →128, and fc3: 128 →q.
  • Appendix B.2. TRE (for Ssup: After training, TRE rescales first-layer weights W1 ∈R256×p to the original input space and uses their right singular vectors as candidate directions.This procedure converts the learned weights into directions for subsequent use.

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

The appendices specify the MLP baseline’s training setup and describe an evaluation MLP used to assess reduced data, with a single-hidden-layer architecture mapping r-dimensional inputs to q outputs.

  • Appendix B.3.2. Training Parameters: Table B.6 reports the training hyperparameters for the MLP baseline.The supplied passage identifies the table but does not list individual hyperparameter values.
  • Appendix B.4. Evaluation MLP (for assessing qualities of estimated CSs): The evaluation MLP is trained on the projected training set.Its training data are the projected observations produced by the estimated reduction.
  • Appendix B.4. Evaluation MLP (for assessing qualities of estimated CSs): The evaluation MLP is early-stopped using the validation set.Validation-based early stopping is part of the model-training procedure.
  • Appendix B.4. Evaluation MLP (for assessing qualities of estimated CSs): The evaluation MLP assesses reduced-data quality rather than estimating directions.The passage explicitly limits its role to evaluating the quality of reduced data.
  • Appendix B.4.1. Architecture: The evaluation MLP uses a single hidden layer.The architecture description identifies one hidden layer before the output layer.
  • Appendix B.4.1. Architecture: Its first layer maps r inputs to 32 units and applies ReLU.The specified layer is fc1: Linear r →32, followed by ReLU.
  • Appendix B.4.1. Architecture: Its second layer maps 32 hidden units to q outputs without an activation function.The specified output layer is fc2: Linear 32 →q (no activation).

Appendix B.4.2. Training Parameters

This appendix section reports the training hyperparameters used for the evaluation MLP.

  • Table B.7 lists the training hyperparameters for the evaluation MLP.

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

The appendix details adaptive rank-selection procedures that use tail diagnostics to choose between energy-based and permutation-based criteria. ARS operates on cross-moment singular values, while SIR and MLP use method-specific matrices or weights.

  • Adaptive Rank Selection: ARS first computes feature-wise skewness and excess kurtosis, then detects heavy tails when smax > s0 or kmax > k0.The algorithm sets smax = maxj |sj| and kmax = maxj kj.
  • Adaptive Rank Selection: Under heavy tails, ARS selects the smallest k whose cumulative singular-value variance reaches η, defaulting to rmax otherwise.The singular values are taken from the sample cross-moment matrix and padded with zeros when necessary.
  • Adaptive Rank Selection: Under light tails, ARS permutes response rows m times, estimates α-quantile thresholds for each singular-value position, and retains consecutive values exceeding them.The returned rank is constrained to at least 1.
  • Rank Selection for SIR: SIR applies the same heavy-tail diagnostic, using slice-averaged matrix singular values for energy selection and a SIR-specific permutation test otherwise.The light-tail branch permutes responses, projects them onto a random direction, slices the whitened covariates, and compares singular values with permutation thresholds.
  • Rank Selection for MLP: MLP rank selection evaluates singular values of the trained first-layer weight matrix, using heavy-tail energy selection and a conservative light-tail energy rule.The light-tail rule handles zero total variance, first-direction dominance via ρ, and cumulative variance before bounding the estimated rank.
Loading 2608.15121v1…