A Principal Square Response Forward Regression Method for Dimension Reduction
Zheng Li Yunhao Wang Wei Gao Hon Keung Tony Ng





Multifractaility, topology and anomalous Hall conductivity on a 30 degrees twisted bilayer honeycomb lattice
    Grigory Bednik
    
============================================================================================================



Multifractaility, topology and anomalous Hall conductivity on a 30 degrees twisted bilayer honeycomb lattice
    Grigory Bednik
    
============================================================================================================



§ ABSTRACT 
Dimension reduction techniques, such as Sufficient Dimension Reduction (SDR), are indispensable for analyzing high-dimensional datasets. This paper introduces a novel SDR method named Principal Square Response Forward Regression (PSRFR) for estimating the central subspace of the response variable Y, given the vector of predictor variables X. We provide a computational algorithm for implementing PSRFR and establish its consistency and asymptotic properties. Monte Carlo simulations are conducted to assess the performance, efficiency, and robustness of the proposed method. Notably, PSRFR exhibits commendable performance in scenarios where the variance of each component becomes increasingly dissimilar, particularly when the predictor variables follow an elliptical distribution. Furthermore, we illustrate and validate the effectiveness of PSRFR using a real-world dataset concerning wine quality. Our findings underscore the utility and reliability of the PSRFR method in practical applications of dimension reduction for high-dimensional data analysis.
  


Keywords: central  subspace; principal square response forward regression;  regression model; sufficient dimension reduction.




§ INTRODUCTION


The convergence of technological advancements, the digitalization of society, the big data paradigm, computational capabilities, and the rise of machine learning and artificial intelligence has fueled the rapid growth of high-dimensional data across numerous domains such as personalized medicine <cit.>, computer vision <cit.>, econometrics <cit.>, and causal inference <cit.>. Sufficient dimension reduction (SDR), a statistical method to extract essential information from high-dimensional data while reducing the dimensionality of the data, stands out as a pivotal tool in the analysis of high-dimensional data. The primary goal of SDR is to identify linear combinations of the independent variables that capture all the relevant information about the conditional distribution of the response variable Y given the vector of predictor variables X, i.e., Y |X. By reducing the dimensionality of the data while retaining the essential information, SDR enables more efficient and effective analysis of high-dimensional datasets. In many practical applications, the underlying parametric model is often unknown. In such cases, <cit.> proposed a general model that assumes an ideal scenario where the high-dimensional vector of predictor variables, X, can be reconstructed from low-dimensional projections for the purpose of regressing Y on X:

    Y=g(β_1^⊤X,β_2^⊤X,…,β_k^⊤X,ε),

where Y is the one-dimensional response variable, X=(X_1,…,X_p)^⊤ is  𝑝-dimensional vector of predictors, g is a ℝ^k+1→ℝ unknown link function, β_i's are unknown non-random column vectors, and ε is the error term which is independent of X. 

Let B=( β_1,…, β_k) ∈ℝ^p× k (k≤ p) be a p× k matrix with columns β_i, i=1,…,k. Then, Y depends on X only through B^⊤X, and the purpose of SDR is to find a matrix B such that 

    Y⊥⊥X|B^⊤X,

where ⊥⊥ represents independence. The space Span(B), Spanned by these linear combinations, is often called the effective dimension reduction (EDR) space. Note that B always exists because B degenerates into an identity matrix when k = p, and it is not unique because if Eq. (<ref>) is true, then Y ⊥⊥X|PB^⊤X for any nonsingular matrix P. Thus, the identifiable parameter here is the subspace Span(B) rather than B itself. <cit.> introduced the central dimension reduction subspace (the smallest) to address the uniqueness problem of EDR space. If a space is an EDR space and this space is contained in any EDR space, then this space is called the central dimension reduction subspace, denoted as S_Y|X. 
Considering the mean function of regression E(Y|X), the purpose of SDR is to find a p × k matrix B such that 

    Y⊥⊥E(Y|X) |B^⊤X,

where Span(B) is called the mean dimension reduction space <cit.>, which is equivalent to 

    E(Y|X)=E(Y|B^⊤X)=h(β_1^⊤X,…,β_k^⊤X).

Similarly, if a subspace is a mean dimension reduction space and it is contained in any mean dimension reduction space, then this subspace is called the central mean dimension reduction subspace (the smallest), denoted as S_E(Y|X). As expressed in <cit.>, S_E(Y|X)⊆S_Y|X, that gives Eq. (<ref>) from Eq. (<ref>). 


In general, the methods to obtain the SDR estimator S_Y|X can be classified into two categories: inverse regression and forward regression <cit.>. 

For inverse regression methods, the sliced inverse regression (SIR) was first proposed by <cit.>, where “inverse regression" refers to the conditional expectation E(X| Y) with Var{E(X| Y)} contained in S_Y|X. This is achieved by assuming a linear conditional mean (LCM) for the basis matrix B, which serves as a fundamental assumption for numerous dimension reduction techniques.
Given that S_Y|X remains invariant when B is multiplied by any k× k full-rank matrix, this property holds for all possible unknown β_i in practice. This equivalence extends to X following an elliptically symmetric distribution.

Inspired by the SIR, other “inverse regression” methods designed for estimating  S_Y|X have been studied. These methods include the sliced average variance estimate (SAVE) based on second-order conditional moments <cit.>,

parametric inverse regression  <cit.>, canonical correlation estimator  <cit.>, contour regression  <cit.>, inverse regression estimator (IRE) <cit.>, principal fitted components <cit.>, likelihood acquired directions <cit.>, directional reduction <cit.>, elliptically contour inverse predictors <cit.>, elliptical sliced inverse regression  <cit.>,  generalized
kernel-based inverse regression <cit.>, and functional SDR estimators <cit.>.

For forward regression methods that focus on the conditional distribution of Y given X, i.e., Y|X. <cit.> first introduced the ordinary least squares (OLS) as a dimension reduction method,  known for its intuitive nature and straightforward algorithm. However, the primary drawback of the OLS method lies in its capability to identify only one vector, and its performance suffers notably when the dimension of S_E(Y|X) exceeds one.  <cit.> proposed principal Hessian directions (PHD) by finding the Hessian matrix for E(Y|X) with the application of Stein's lemma.  
Based on utilizing an iterative approach to estimate S_E(Y|X) with the PHD method, 
<cit.> introduced the 
iterative Hessian transformation (IHT) method. This method was further studied by  <cit.>. Additionally, <cit.> proposed the generalized PHD (GPHD), and <cit.> proposed the adjusted PHD (APHD) methods for mixture multivariate skew elliptical distributions and non-Gaussian predictors, respectively. The IHT, GPHD, and APHD methods can be applied to a wider range of scenarios due to their less restrictive conditions compared to the PHD method. Other forward regression dimension reduction methods, such as the minimum average variance estimator (MAVE) <cit.>,  sliced regression  <cit.>,  ensemble of minimum average variance estimators <cit.>, semiparametric dimension reduction method <cit.>, outer-product-gradient method (OPG) <cit.> and optimal SDR <cit.>, have been developed in the literature.






This paper introduces a principal square response forward regression (PSRFR) estimator for dimension reduction. The proposed approach leverages the OLS estimator from a fresh angle, thus resolving the challenge of OLS's limited capability to recover only one dimension. 
In contrast to the PHD method, we demonstrate that the proposed PSRFR approach may be applicable under the assumption of elliptical distributions. Moreover, the PSRFR method offers greater simplicity and intuitiveness compared to the IHT method. Additionally, it can identify more central subspace directions in certain scenarios than the PHD and IHT methods.

The rest of this paper is organized as follows. In Section <ref>, we propose the PSRFR estimator and derive its consistent and asymptotic normal distribution. A simulation study is conducted in Section <ref>, demonstrating that the PSRFR estimator outperforms existing methods when the variance of each component of X is significantly different under the assumption of an elliptical distribution, and it also showcases its robustness. In Section <ref>, we investigate the Wine Quality dataset through a real data analysis. Section <ref> provides some concluding remarks for this paper.




§ PROPOSED METHODOLOGY


In this section, we introduce the proposed PSRFR estimator for Span(B), which is the smallest dimension reduction subspace we are interested in, and examine its associated theoretical properties.
Here, we assume that the structural dimensions (the rank of matrix B) are known, we do not attempt to determine the dimension of the central subspace in this paper. And hence, estimating Span(B) is equivalent to estimating the direction of the basis of  Span(B). For convenience, we consider that B satisfies B^⊤B=_k, where _k is a k-dimensional identity matrix.




 §.§ PSRFR Estimator


To introduce the proposed PSRFR estimator, we start with the OLS estimator under the elliptical distributions with the following assumption. 

 

The distribution of X is an elliptical distribution with mean E(X)=0 and variance-covariance matrix Var(X)=Σ_X.

The following Lemma <ref> shows that E(YX) fall in the central mean dimension reduction subspace S_E(Y|X). 


For Y and X that satisfy Eq. (<ref>), Assumption <ref> implies

    E(YX)=Σ_XΛ,

where Λ=(λ_1, …,λ_k)^⊤ is a constant vector. 
 
The proof of Lemma <ref> is provided in the Appendix . <cit.>, <cit.> <cit.>,  <cit.>, and <cit.> similarly provide results analogous to those in Lemma <ref>.  


Note that the OLS estimator is a vector representation of a basis in S_E(Y|X), which is only a precise estimator when the structure dimension is one. According to Eq. (<ref>), E(YX) is a linear combination of {β_i}_i=1^k,  
which indicates that  E(YX) lies in the hyperplane Spanned by {β_i}_i=1^k. To obtain a complete set of the basis of  S_E(Y|X), it is required to determine a set of standard orthogonal basis for S_E(Y|X). Without loss of generality, we assume that X has mean 0 and an identity matrix as the variance-covariance matrix, and {(Y_j,X_j)}_j=1^n is a set of independent and identically distributed samples from Eq. (<ref>). 

The crux of the issue lies in estimating the basis of S_E(Y|X) using {Z_j}_j=1^n, where Z_j=Y_jX_j. This can be intuitively solved by minimizing the sum of distances from all sample points to the hyperplane spanned by the basis matrix B. Given that B^⊤B=_k, the distance from Z_j to the hyperplane Span(B) can be expressed as follows:


    d_j =Z_j-BB^⊤Z_j_2^2 
    = Z_j^⊤Z_j-Z_j^⊤BB^⊤Z_j,

where ||·||_2 represents L_2 norm. Then, the minimization problem to obtain the estimator for Span(B) can be expressed as

    min_B∑_j=1^n d_j 
       =
    min_B∑_j=1^n(Z_j^⊤Z_j-Z_j^⊤ BB^⊤Z_j)
    
       =min_{β_i}_i=1^k∑_j=1^n (Z_j^⊤Z_j- ∑_i=1^kZ_j^⊤β_iβ_i^⊤Z_j ),

which is equivalent to the maximization problem

    max_{β_i}_i=1^k∑_j=1^n∑_i=1^kZ_j^⊤β_iβ_i^⊤Z_j   =max_{β_i}_i=1^k∑_j=1^n∑_i=1^kβ_i^⊤Z_jZ_j^⊤β_i
       =max_{β_i}_i=1^k∑_i=1^kβ_i^⊤ ( ∑_j=1^nZ_jZ_j^⊤ )β_i.

It can be shown that the solution of the above optimization problem is the first k eigenvectors corresponding to the first k eigenvalues of ∑_j=1^nZ_jZ_j^⊤. 

For the practical situation that the mean E(X) and variance-covariance matrix of X, Σ_X, are unknown, 
the sample mean X and the sample variance S_n based on the observed sample can be used to approximate E(X) and Σ_X, respectively. Then, the data can be transformed as Z_j=S_n^-1Y_j(X_j-X̅). 

The aforementioned results describe how the PSRFR estimator of Span(B) can be obtained. Hence, a set of standard orthogonal bases in Span(B) can also be estimated similarly. The algorithm to obtain the PSRFR estimator, namely the PSRFR algorithm, is presented in Algorithm 1.



At the population level, the PSRFR estimator is based on the p × p positive-defined matrix

    Σ_X^-1E[Y^2{X-E(X)}{X-E(X)}^⊤]Σ_X^-1,

which is associated with the principal components analysis (PCA) and the PHD method, as described below.

For PCA, the principal components of X are defined as the set of linear combinations of X that have the largest variances. The first principal component of X can be obtained by solving the maximization problem

    max_||γ ||_2=1γ^⊤E[{X-E(X)}{X-E(X)}^⊤]γ,

in which the solution is the first eigenvector of E[{X-E(X)}{X-E(X)}^⊤]. Similarly, the first k principal components of X are the first k eigenvectors of Σ_X. After obtaining the principal components, a regression model such as the one presented in Eq. (<ref>) can be constructed. 
However, projecting the p-dimensional predictor variables onto lower dimensions first might inadvertently result in an inaccurate relationship between Y and the original X. In contrast, the proposed PSRFR method systematically establishes the connections between X and Y at the beginning of the dimension reduction process.



For the PHD method, it is based on the Hessian matrix of twice differentiable regression function E(Y|X), denoted as H_m(X). By the chain rule, we have  

    H_m(X) =∂^2E(Y|X)/∂X∂X^⊤=B∂^2E(Y|B^⊤X)∂ (B^⊤X)∂ (X^⊤B)B^⊤.

Let H_m(X)=E{H_m(X)}. If X follows a normal distribution,  H_m(X) can be transformed into an easily solvable form by Stein's lemma, Σ_X^-1Σ_YXXΣ_X^-1,
where 

    Σ_YXX=E[{Y-E(Y)}{X-E(X)}{X-E(X)}^⊤],

which is known as the response-based (y-based) version of the PHD method. In contrast to the PHD approach, besides extending the normal distribution assumption to the elliptical distribution assumption, the proposed PSRFR replaces Y with Y^2 in Eq. (<ref>), potentially leading to notable effects on the estimators.




    The major difference between the proposed PSRFR approach and the  PHD approach is the term  E(Y^2|X) and E(Y|X) involved according to the law of iterated expectation. Consider the following model in <cit.>: 
    Y=f(β_1^⊤X)+g(β_2^⊤X)ε,

    where ε is dependent on X with zero mean,  β_1,β_2∈ℝ^p, and f and g are unknown link functions. 
    Since E(Y|X)=f(β_1^⊤X) and  E(Y^2|X)=f^2(β_1^⊤X)+g^2(β_2^⊤X)Var(ε), the PSRFR method can identify more central subspace directions than the PHD method. Moreover, 
    
    E(Y^2|X)=Var(Y|X)+E(Y|X)E(Y|X)

    according to the definition of conditional variance. Here, E(Y^2|X) contains Var(Y|X), and not just E(Y|X).



 Remark <ref> shows the PSRFR might identify more central subspace directions in the above example under the model in Eq. (<ref>); hence, the EDR space estimated by the PSRFR is no longer limited to the central mean subspace. Furthermore, we are interested in  

    E(Y^2|X)=E(Y^2|B^⊤X)=H(β_1^⊤X, …,β_k^⊤X).



Due to the non-uniqueness of B and for ease of expression in the remainder of the paper, Span(B) is used to represent the dimension reduction subspace in Eq. (<ref>), where Span(B) ⊆S_Y|X, which makes it easy to derive Eq. (<ref>) from Eq. (<ref>).   




 §.§ Asymptotic Properties


In this subsection, we will prove that {β̂_i}_i=1^k converge to a set of standard orthogonal basis for Span(B) under some mild conditions. 


Under the model in Eq. (<ref>), the inequality 

    E{Y^2(β^⊤Σ_X^-1X)^2}>E{Y^2(α^⊤Σ_X^-1X)^2}

holds, where β∈Span(B), α∈Span(B)^⊥, β_2=α_2=1, and the symbol ⊥ represents the orthogonal complement.
 
In Assumption <ref>, the term on the left-hand side of the inequality in Eq. (<ref>) 

    E{Y^2(β^⊤Σ_X^-1X)^2}=(β^⊤BΛ)^2+Var(Yβ^⊤Σ_X^-1X),
 
which represents the sum of a fixed term (β^⊤BΛ)^2  
and the variance of some random variables, 
and the term of the right-hand side of the inequality in Eq. (<ref>) 

    E{Y^2(α^⊤Σ_X^-1X)^2}=Var(Yα^⊤Σ_X^-1X)

represents the variance of some random variables.  



For Model in Eq. (<ref>), under Assumptions <ref> and <ref>, the first k eigenvectors corresponding to the first k eigenvalues of  E(ZZ^⊤) are the basis of Span(), where Z=Σ_X^-1YX, and E(ZZ^⊤)-Σ_X^-1E() is contained in Span(), where E()=E( Y^2 w_k+1^2), w_k+1 is the (k+1)-th element of  
  W=X=(w_1,…,w_k,w_k+1,…,w_p)^⊤, and =(β_1,…,β_k,α_k+1,…,α_p)^⊤≡( ,)^⊤ is an orthogonal matrix.

The proof of Theorem <ref> is provided in Appendix . 

Although E(ZZ^⊤) does not fall completely in Span(B), a set of basis can still be found through the eigendecomposition of E(ZZ^⊤), namely, 

    E(ZZ^⊤)=[ _1 _0 ][ Ψ_1   0;   0 Ψ_0 ][ _1^⊤; _0^⊤ ]=_1Ψ_1_1^⊤+_0Ψ_0_0^⊤,

where = ( _1,_0 ) is a p-dimensional orthogonal matrix, Ψ_1 and Ψ_0 are diagonal matrices with dimensions k and p-k, respectively. The diagonal elements in Ψ_1 and Ψ_0 are ordered from large to small.  
In the proof of Theorem <ref>, we show that 

    E(ZZ^⊤)=BΓ _1B^⊤+AΓ _2A^⊤,

where Γ_1 is a k × k positive-defined matrix with the (i, j)-th elements is E{ H(w_1,…,w_k) w_iw_j}, Γ _2 is a (p-k) × (p-k) diagonal matrix with diagonal elements E(G).
 Consider that Γ_1 is not a diagonal matrix in Eq. (<ref>), we can rewrite BΓ _1B^⊤=BVΦV^⊤B^⊤ by eigendecomposition, where V and Φ are k-dimensional orthonormal and diagonal matrices, respectively. 
Then, it is sufficient to show Span(B)=Span(BV)=Span(Q_1), where the pivotal challenge arises in identifying Span(B) through eigendecomposition of E(ZZ^⊤). 
Under Assumption <ref>, the eigenvalues corresponding to the eigenvectors of Γ_1 exceed those corresponding to the eigenvectors of Γ_2, which guarantees that the first k eigenvectors of E(ZZ^⊤) corresponding to the first k eigenvalues must be the basis of Span(B).
Thus, the basis can be determined by finding the first k eigenvectors corresponding to the first k eigenvalues. 



 For the model in Eq. (<ref>), if X follows a multivariate Gaussian distribution, then  Σ_X^-1E[{Y^2-E(Y^2)}XX^⊤]Σ_X^-1 is contained in Span(B).
In the case of Remark <ref>, the proposed PSRFR method can still identify more central subspace directions. More details are provided in the proof of Theorem <ref> presented in Appendix . 
Moreover, based on the idea that square response is contained in the PSRFR method, if one first transforms the response variable Y into Y^2 and then applies those SDR methods that focus on the E(Y|X), 
more central subspace directions can be identified in the case of Remark <ref>. 


Following Theorem <ref>, the asymptotic properties of the estimator {β̂_i}_i=1^k can be established, and the results are given in Theorem <ref>.



For the model in Eq. (<ref>), under Assumptions <ref> and <ref>, if  E{E(Y|X)^2}
<∞ holds, then

    Span(β̂_1,…,β̂_k)Pr⟶Span(β_1,…,β_k),

and Span() converges to Span() at rate n^1/2. 
In addition, if  Var{vec(ZZ^⊤) } exists, then 

    √(n) [ vec(𝒵̂)- vec{E(ZZ^⊤)} ]L⟶N[ 0, Var{vec(ZZ^⊤) }],

 where vec(·) is the operator that maps a symmetric matrix to a vector by stacking the 
main diagonal and the elements below the main diagonal by columns, i.e., if S is a p × p symmetric matrix 

    S = 
    [ s_11 s_12 s_13    ⋯ s_1p; s_12 s_22 s_23    ⋯ s_2p; s_13 s_23 s_33    ⋯ s_3p;    ⋮    ⋮    ⋮    ⋯    ⋮; s_1p s_2p s_3p    ⋯ s_pp;      ],

then vec(S) = (s_11, …, s_1p, s_22,…, s_2p, s_33, …, s_3p, …, s_pp)^⊤. 



The proof of Theorem <ref> is provided in Appendix . 
Note that 

    E(Y^2)=E{E(Y^2|X)}=E{H(β_1^⊤X,⋯,β_k^⊤X)}≥E{E(Y|X)^2},

hence, in Theorem <ref>, instead of the condition E{E(Y|X)^2}<∞, we can use E(Y^2)<∞.  
Because 𝒵̂, ZZ^⊤, and E(ZZ^⊤) are p × p symmetric matrices, the dimensions of vec(𝒵̂), vec(ZZ^⊤ ), and vec{E(ZZ^⊤)} are p × (p+1)/2. 
Notice that Var{vec(ZZ^⊤) } is a p(p+1)/2 by p(p+1)/2 matrix. We represent the Var{vec(ZZ^⊤) } by the Kronecker product ⊗ in the proof of Theorem <ref>,  the dimension of (ZZ^⊤)⊗(ZZ^⊤) is p^2 by p^2. Although the Kronecker product representation has a degenerate variance (all the elements in symmetric matrices ZZ^⊤ are used), this representation is for notation convenience only as our concern is the convergence of each element of variance of the random matrix ZZ^⊤.


Theorem <ref> shows that  Span(B̂) are √(n)-consistency estimators of Span(B) by the law of large number. Although the first k eigenvectors corresponding to the first k eigenvalue of E(ZZ^⊤) are not a set of original basis in Eq. (<ref>), they represent the same space. Moreover, 𝒵̂=n^-1∑_j=1^n(Z_jZ_j^⊤) is a √(n)-consistency estimator of E(ZZ^⊤). 
Additionally, if Var{vec(ZZ^⊤) } exists, the asymptotic normality property is obtained by the central limit theorem. 



§ MONTE CARLO SIMULATION STUDY


In this section, we evaluate the performance of the proposed PSRFR method by using a Monte Carlo simulation study. To measure the distance between the true subspace Span(B) and the corresponding estimator Span(B̂) for B̂=(β̂_̂1̂, …,β̂_̂k̂), we consider the trace correlation defined as  <cit.> 

    R   = trace (P_BP_B̂ )k,

where P_B=B(B^⊤B)^-1B^⊤ denotes the projection matrix. 
Without loss of generality, we assume B̂ is a column-orthogonal matrix due to the property of B. Otherwise, the Gram-Schmidt ortho-normalization method can be used and will not change the subspace. Then, the trace correlation based on the estimator B̂ in Eq. (<ref>) can be calculated as

    R  =trace (B̂^⊤BB^⊤B̂ )k.

Here, the trace correlation R can be used to evaluate and compare the performance of different estimation methods. The trace correlation is a value between 0 and 1, and a larger value of R indicates a better estimator B̂.  In the following subsections, we consider X follows elliptical distributions in Section <ref> to investigate the validity of the proposed PSRFR method, and X follows non-elliptical distributions in Section <ref> to evaluate the robustness of different methods. 



 §.§ Elliptical distributions

 

In this subsection, two types of elliptical distributions,
normal and non-normal distributions are considered as the distribution of the predictor variables X.



  §.§.§ Normal distribution



First, we compare the proposed PSRFR method to the PHD and IHT methods under the model in <cit.> described in Remark <ref>. The following two models are considered in the simulation study: 
 

  * Model [N1]

    Y=β_1^⊤X+β_2^⊤X·ε.


  * Model [N2]

    Y=sin(β_1^⊤X)+(|β_2^⊤X+1|)^1/2·ε.
 
 
We consider p=dim(X) = 10, X∼N_10(0, Σ_norm),  ε∼N(0,1),  β_1=(1,0,0,…,0) and β_2=(0,1,0,…,0), where Σ_norm is a diagonal matrix with the diagonal elements (1,2,3, …, 10), and the sample sizes are n = 100, 300 and 500.


We use the trace correlation as a comparison criterion and also consider the cosine similarity criteria. 
Specifically, the cosine similarity for β_1 is defined as  

    |cos_1|=max{|β̂_1^⊤β_1|/β̂_1^⊤, |β̂_1^⊤β_2|/β̂_1^⊤},

which is the absolute value of the cosine of the closest true direction to β̂_1^⊤. Similarly, the cosine similarity for β_2 is defined as  

    |cos_2|=max{|β̂_2^⊤β_2|/β̂_2^⊤, |β̂_2^⊤β_1|/β̂_2^⊤}.

Larger values of |cos_1| and |cos_2| indicate a better estimator B̂. 

 The computer program written in R <cit.> for the implementation of the proposed PSRFR method is provided in Appendix . The PHD and IHT methods are implemented in the R packages dr  <cit.> and itdr <cit.>, respectively.
The averages and standard deviations (SDs) of 
R, |cos_1| and |cos_2| for the proposed PSRFR method, the PHD method, and the IHT method based on 1000 simulations are reported in Table <ref>. 

 



 

From the results in Table <ref>, the performance of the methods considered here improved with the increase in sample size. 
We observe that the PSRFR method can identify the whole subspace more accurately and estimate each direction well compared with the PHD and IHT methods, and the IHT is more accurate at recognizing the first direction.

In addition to the PHD and IHT methods, we further consider comparing the PSRFR method to the SIR, SAVE, and IRE methods under three different models with normally distributed predictor variables studied in Example 3 of <cit.>, <cit.>, and Example 3 of <cit.>: 
  

  * Model [N3]

    Y=(4+β_1^⊤X)· (β_2^⊤X+2)+σε.


  * Model [N4]

    Y=β_1^⊤X/{0.5+(β_2^⊤X+3)^2}+σε.


  * Model [N5]

    Y=( β_1^⊤X)^2+ | β_2^⊤X | +σε.


We consider σ=0.5, the aforementioned settings for X, σ, β_1, and β_2, and the number of slices to be H = 10. 
 The SIR, SAVE, and IRE methods are also implemented in the R package dr <cit.>.  Table 2 reports the averages and standard deviations of the trace correlation R based on 1000 simulations.


From Table <ref>, once again, the performance of the methods considered here improved with the increase in sample size. The PSRFR method performs well in almost all the settings considered here, with different variances of each predictor component. Compared to the PHD and IHT methods, the PSRFR underperforms under model [N3], especially for sample size n = 100, since model [N3] satisfies the assumptions for the PHD method is developed based on the normality of the predictor variables, and the IHT method is based on the PHD under the elliptically distributed predictor variables assumption.
The PSRFR method performs better than the other methods considered here under models [N4] and [N5], which indicates that the PSRFR method is an effective method for estimating Span(B) in the multivariate normal case with different variances of the predictor variables.




  §.§.§ Non-normal distributions



In this subsection, we consider the predictor variables following different non-normal multivariate elliptical distributions. Specifically, the multivariate Student's t and multivariate power exponential distributions are considered: 


   Multivariate Student's t distribution <cit.>: 
 
A p-dimensional random vector X is said to be distributed as a multivariate Student's t distribution with degrees of freedom ν, mean vector μ, and positive-definite symmetric matrix Σ if its joint probability density function is given by

    f_t(X)     =    Γ((ν+p)/2)/(πν)^p/2Γ(ν/2)|Σ|^1/2
          ×[1+1/ν(X-μ)^⊤Σ^-1(X-μ)]^-(ν+p)/2,  X∈ℝ^p.

As ν→∞, the limiting form is the multivariate normal distribution. Hence, multivariate Student's t distribution with small degrees of freedom deviates significantly from multivariate normal distribution, especially in the tail areas. In the special case of ν =1, the multivariate Student's t distribution is a multivariate Cauchy distribution. Notice that the variance-covariance matrix of the multivariate Student's t distribution is given by ν/(ν-2)Σ for ν>2. Hence, for multivariate Student's t distribution with degrees of freedom 1 and 2, the variance-covariance matrix does not exist. 

   Multivariate power exponential distribution <cit.>:

A p-dimensional random vector X is said to be distributed as a multivariate power exponential distribution with kurtosis parameters β>0, mean vector μ and positive-definite symmetric matrix Σ if its joint probability density function is given by

    f_PE(X)    =    pΓ(p/2)/Γ(1+p/2β)π^p/22^1+p/2β|Σ|^1/2
          ×[-1/2((X-μ)^⊤Σ^-1(X-μ))^β],  X∈ℝ^p.

In particular, the multivariate Laplace and multivariate normal distributions are special cases of the multivariate power exponential distribution when the kurtosis parameter β=0.5 and β=1, respectively. Therefore, the kurtosis parameter β can be viewed as the disparity between power exponential distribution and normal distribution.
 

In the simulation study, we consider the predictor variables following the multivariate Student's t distributions with degrees of freedom ν = 2 and 3, the multivariate Cauchy distribution (i.e., multivariate Student's t distributions with degree of freedom ν = 1), and the multivariate power exponential distribution with kurtosis parameters β = 0.5 and 5. We utilize the R packages mvtnorm <cit.> and LaplacesDemon <cit.> to simulate random vectors from the multivariate Student's t and the multivariate power exponential distributions, respectively. 








The following four models under non-normal elliptical distributed predictor variables are considered: 
 

  * Model [NN1]: 

    Y=(4+β_1^⊤X)+ (β_2^⊤X+2)·σε^2.

This model is from Example 2 of <cit.>.

  * Model [NN2]: 

    Y=(|4+β_1^⊤X|)^1/2·(|β_2^⊤X+2|)^1/2 +σε.

This model is based on model [NN1] with a slow-growing power function of degree 1/2. 

  * Model [NN3]: 

    Y=(|β_1^⊤X|)^1/2 +(|β_2^⊤X·ε|)^1/2 +σε.

This model is also based on model [NN1] with a slow-growing power function of degree 1/2. 

  * Model [NN4]: 

    Y=0.4·(β_1^⊤X)+3·sin ( β_2^⊤X/4  )+σε.

This model is motivated by Example 1 of <cit.>.
 


Following the settings in Section <ref>, we consider dim(X) = 10, σ=0.5, ε∼N(0,1), β_1=(1,0,…,0) and β_2=(0,1,…,0). We set μ = 0 and Σ = Σ_ellp  is a diagonal matrix with elements (1, 6, 11, 16, 21, 26, 31, 36, 41, 46) in the multivariate Student's t and power exponential distributions. 

The averages and standard deviations of the trace correlation R for different methods based on 1000 simulations are reported in Tables <ref>–<ref> for models [NN1] – [NN4], respectively.


The results in Tables <ref>–<ref> show that the PSRFR method outperforms other methods in most of the models and settings when the predictor variables follow a non-normal elliptical distribution, especially when the distribution has heavier tails compared to the multivariate normal distribution (i.e., the multivariate Student's t  distribution with small degree of freedom ν, and the multivariate power exponential distribution with large kurtosis parameter β). 

The PHD and SAVE methods exhibit comparable performance to the PSRFR method when the predictor variables follow the multivariate power exponential distribution with kurtosis parameter β = 0.5 since the multivariate power exponential distribution with small kurtosis parameter behaves similarly to multivariate normal distribution. 














 §.§ Non-elliptical Distribution for Robust Analysis



In this subsection, we perform a simulation study to investigate whether the PSRFR method can effectively identify Span(B) when the predictor variables do not follow an elliptical distribution. Under the non-elliptical distribution situations, we consider comparing the proposed PSRFR method to the MAVE and the OPG methods proposed by <cit.> for identifying the central mean subspace. The MAVE and the OPG methods require a differentiable link function but have no strict distributional assumptions for the predictor variables. The OPG and MAVE methods are implemented in the R package MAVE <cit.>. 

Before considering the non-elliptical distribution situations, we compare the 
OPG and MAVE methods to the PSRFR methods under model [NN1] and [NN4] when the predictor variables follow an elliptical distribution. Specifically, in Table <ref>, we present the averages and standard deviations of the trace correlations based on 1000 simulations for PSRFR, OPG, and MAVE methods when the predictor variables follow a multivariate normal distribution or multivariate Student's t distribution with degrees of freedom ν = 3 with the setting described in Sections <ref> and <ref>. 
The results in Table <ref> show that the PSRFR method performs better than the OPG and MAVE methods when the predictor variables follow an elliptical distribution.

For the non-elliptical distribution situation, the predictor variables X are generated from a mixture of multivariate normal and multivariate uniform distribution in (-3, 3) with mixture probabilities 0.8 and 0.2, respectively, i.e.,

    X∼ 0.8N_10(0, Σ_norm)+0.2 U_10(-3,3),

where U_10(-3,3) denotes a 10-dimensional multivariate uniform distribution in (-3, 3), each component of which independently follows uniform distribution from -3 to 3. The following three models are considered:
 

  * Model [NE1]:

    Y=β_1^⊤X/{0.5+(β_2^⊤X+1.5)^2}+σε.


  * Model [NE2]:

    Y=β_1^⊤X·(β_2^⊤X+1) + σε.


  * Model [NE3]:

    Y=0.4·(β_1^⊤X)+3·sin ( β_1^⊤XX^⊤β_2/4  )+σε.

 
Following the settings in Sections <ref> and <ref>, we consider σ=0.5, ε∼N(0,1), β_1=(1,0,…,0),  β_2=(0,1,…,0) and 
Σ_norm is a diagonal matrix with elements (1, …, 10). 


The averages and standard deviations of the trace correlations based on 1000 simulations for PSRFR, OPG, and MAVE methods when the predictor variables follow the distribution in Eq. (<ref>) are presented in Table <ref>. From Table <ref>, the PSRFR performs about 10% worse relative to the OPG and MAVE method in terms of the trace correlation when n = 100, and the differences decrease 
to less than 5% when the sample size increases to n = 500. 
Although the OPG and MAVE methods perform better than the PSRFR method under the non-elliptical distribution situations, these simulation results demonstrate that the proposed PSRFR method is robust to the underlying distribution of the predictor variables. 
The proposed PSRFR method 
can effectively identify central subspaces, even the distribution of the predictor variables deviations from the elliptical distribution. 










 §.§ Effect of the dimension of predictor variables



In the simulation studies presented in Sections <ref> and <ref>, the dimension of predictor variables is considered as p = (X) = 10 by following the classical works on SDR. In this subsection, we examine the performance of the proposed PSRFR method when the dimension of predictor variables p is larger than 10. 

In this simulation study, we consider p=30 and 40 with 
β_1=(1,0,0,…,0) and β_2=(0,1,0,…,0) under model [N4] with multivariate normal distributed predictor variables and under model [NN3] with multivariate Cauchy distribution, where the matrix 
Σ_norm is a diagonal matrix with elements  (1, 1, 1, 2, 2, 2, …,
10, 10, 10) for p = 30 and  (1, 1, 1, 1, 2, 2, 2, 2, …, 10, 10, 10, 10) for p = 40, and the matrix and Σ_ellp 
is a diagonal matrix with elements (1, 1, 1, 6, 6, 6, 11, 11, 11, …,
46, 46, 46) for p = 30 and    (1, 1, 1, 1, 6, 6, 6, 6, …, 
 46, 46, 46, 46) for p = 40. 
We compare the performance of the proposed PSRFR method with the SIR and IHT methods when p = 30 and 40. Table <ref> presents the averages and standard deviations of the trace correlation R for PSRFR, SIR, and IHT methods based on 1000 simulations.






From Table <ref>, we observe that the performances of PSRFR, SIR, and IHT methods depreciate when the dimension of X increases and the sample size decreases. This observation is consistent with the intuitive realization. Moreover, we observe that the PSRFR method still performs reasonably well (R close to or greater than 0.7 in most cases) with the dimension of X being 30 and seems less sensitive to the increase in dimensionality when compared to the SIR and IHT methods. 

 

 §.§ Effect of the general basis vectors and different noise levels


In the simulation studies presented in Sections <ref> ,<ref> and <ref>, the number of basis vectors is considered as two and the noise level is 0.5. In this subsection, we examine the performance of the proposed PSRFR method when the basis vector becomes general under different noise levels.

In this simulation study, we consider σ=2 and σ=4 with 
β_1=(1/√(2),1/√(2),0,
…,0), β_2=(1/√(2),-1/√(2),0,…,0), β_3=(0,0,1/√(2),1/√(2), 0,…,0) and β_4=(0,0,1/√(2),-1/√(2),0,…,0) under following model with 10-dimensional multivariate normal distributed predictor variables in Sections <ref>. 

    
  * Model: 
    
    Y=sin(β_1^⊤X+4)+exp(β_2^⊤X )+(β_3^⊤X)^2+|β_4^⊤X| +σε.



We compare the performance of the proposed PSRFR method with the PHD, SIR, SAVE, IHT, OPG and MAVE methods. Table <ref> presents the averages and standard deviations of the trace correlation R based on 1000 simulations.



From Table <ref>, we observe that the more general basis vectors with different noise levels do not have much effect on the performance of the PSRFR method.




 In summary, as demonstrated by the simulation results in Sections <ref>–<ref>, the proposed PSRFR method exhibits promising performance, particularly evident when the variances of individual components diverge, and the tails of predictor variable distributions become heavier. Notably, PSRFR maintains its robustness and reliability even in scenarios where predictor variable distributions deviate from the elliptical distribution assumption and for large dimensions of the predictor variables. This versatility renders PSRFR applicable to a broader spectrum of real-world data analyses, enhancing its practical utility.



§ REAL DATA ANALYSIS


The Wine Quality dataset presented in <cit.> is a publicly available data set that contains two sub-datasets: the red wine and white wine data sets. 
The response variable is the quality of wine (QL), and the 11 predictor variables are: fixed acidity (FA); volatile acidity (VA); citric acid (CA); residual sugar (RS);  chlorides (CL); free sulfur dioxide (FSD); total sulfur dioxide (TSD);  density (DS); pH value (PH);  sulphates (SP); and alcohol level (AH). <cit.> studied the relative importance of 11 predictor variables for wine quality using the support vector machine approach and pointed out that a regression approach on the two sub-datasets can be used. Here, we consider the regression approach and apply the proposed PSRFR method to obtain the relative importance with the first 1599 and 800 observations in the red and white wine sub-datasets, respectively.  
And  we compare to the result obtained by the support vector machine approach.

First, we assess the normality of the 11 predictor variables using the hypothesis testing approach and graphical approach, namely the Shapiro-Wilk test and normal quantile-quantile (Q-Q) plot, after standardizing the data. 
Table <ref> presents the Shapiro-Wilk statistics for predictor variables in the red wine and white wine data sets. The normal Q-Q plots are depicted in Figure <ref> and <ref> for the red wine and white wine data sets, respectively. From Table <ref> and the normal Q-Q plots in Figure <ref> and <ref>, except for variables TSD and PH in the white wine data set, all the other variables in both data sets are not likely to follow a normal distribution.  

To assess the symmetry of the distributions of these 11 predictor variables in the red wine and white wine data sets, we provide the comparative boxplots in Figure <ref> and <ref>. From Figure <ref> and <ref>, we observe that 
the variables in the red wine and white wine data sets are asymmetric and have heavy-tailed characteristics.

Considering that the PSRFR method performs well when the variances of the predictors are significantly different from each other, we diagonalize the data by performing an eigen-decomposition on the sample covariance matrix and multiply the centered data by the corresponding eigenvectors. Then, we apply the PSRFR method to address the regression problem at hand. Once the PSRFR estimator is obtained, the next step involves determining the dimension of the central subspace. We consider evaluating the proportion that an eigenvalue accounts for all the eigenvalues to determine the dimension of the central space, which is similar to the way of determining the number of principal components in PCA. The results show that the proportions of the first eigenvalues in the red wine and white wine data sets account for all the eigenvalues 
are 0.9995 and 0.9997, respectively, which indicate that the dimension of central space can be considered as 1. 

Next, our objective is to estimate the basis β̂_1 of the central space. To account for relative importance, we take the absolute value of each component in β̂_1 and arrange them in descending order, from largest to smallest. In the red wine data set, the relative importance in descending order is AH, PH, SP, DS, FSD, TSD, CL, RS, CA, VA, and FA. In the white wine data set, the relative importance in descending order is AH, SP, TSD, DS, CL, FSD, RS, PH, CA, FA, VA. These results align with the findings in <cit.>, in which PH, SP, and AH are also identified as important variables in the red wine data set, while RS and CA were relatively unimportant. In the white wine data set, SP and AH are relatively important, while FA and PH are relatively unimportant. 

We observe that both SP and AH are highly important variables in both red wine and white wine data sets. <cit.> provided a physiological explanation about this and suggested that an increase in sulfates may be associated with fermenting nutrition, which plays a crucial role in enhancing the wine aroma. The significance of AH in wine is evident. Furthermore, the importance of pH value (PH) in red wine surpasses that in white wine. Although SP and AH are consistently identified as important variables, it is noteworthy that SP holds the highest importance in <cit.>, and AH emerges as the most important variable in our findings. <cit.> suggests that an increase in alcohol content tends to result in higher-quality wine. 
Furthermore, considering that AH holds the highest significance in our results and DS is influenced by the proportion of AH and sugar content, it can be inferred that DS may have a greater importance than initially suggested in the study by <cit.>. 





§ CONCLUDING REMARKS


In this paper, we propose a principal square response forward regression (PSRFR) method, a novel approach for dimension reduction tailored for high-dimensional, elliptically distributed data. Drawing inspiration from the Ordinary Least Squares (OLS) method, PSRFR is devised to handle the complexities of such datasets effectively.

The core principle of PSRFR lies in leveraging the amalgamated information from both predictor and response variables, which typically concentrates around a central subspace. Unlike the OLS method, which tends to recover only a single dimension, PSRFR aims at capturing a comprehensive estimate of this central subspace. The PSRFR method achieves this by minimizing the distance between data points and the central subspace, thus identifying multiple central subspace directions, surpassing the capabilities of methods like PHD and IHT methods when the predictor variables follow an elliptical distribution. Moreover, this paper presents a fundamental theorem affirming the efficacy of PSRFR in achieving substantial dimension reduction. Additionally, we provide a simple algorithm for implementing PSRFR. We delve into the asymptotic behavior of the PSRFR estimator, elucidating its convergence rate in high-dimensional scenarios.

Our simulation results underscore the superiority of PSRFR in enhancing estimation accuracy for elliptically distributed data with varying component variances. This improvement is facilitated through a simple data transformation process. Overall, the proposed PSRFR method furnishes invaluable tools for dimension reduction and analysis of high-dimensional, elliptically distributed data, exhibiting resilience even in the face of deviations from the elliptical distribution. 

Note that determining the dimension of the central subspace is a critical issue, which we leave for further study.



 §.§ Competing interests

The authors declare there are no conflict of interests.



 §.§ Funding

The research is supported by the National Natural Science Foundation of China (No. 12371263). 




§ APPENDIX



 §.§ Proof of Lemma <ref>




For  B=(β_1,…,β_k), without loss of generality, we consider E(X) = 0 (otherwise, we can consider transforming  X by X - E(X)). 


Case 1:  Σ_X=_p. In this case, we have 

    E(YX)   =E{E(YX|X)}
       =E{X·E(Y|X)}
    =E{h(β_1^⊤X,⋯,β_k^⊤X)·X}.

Let C be the p× p orthogonal matrix with the first k rows be β_i^⊤, i∈{1,…,k}. 
Then, define 

    C =(β_1,…,β_k,α_k+1,…,α_p)^⊤≡(B, A)^⊤,

and let

    W   =CX=
    (β_1^⊤X,⋯,β_k^⊤X, α_k+1^⊤X,⋯,α_p^⊤X)^⊤
       =(w_1,⋯,w_k,w_k+1,⋯,w_p)^⊤.

Hence,

    E(YX)   =E[h(β_1^⊤X,⋯,β_k^⊤X)·X]
       =C^⊤ E[h(β_1^⊤X,…,β_k^⊤X)· CX]
       =C^⊤ E[h(w_1,…,w_k)·W].

Note that W = CX also follows elliptical distributions with with mean E(W) = 0 and variance-covariance matrix C_pC^⊤=_p, where _p is a p-dimensional identity matrix. By Theorem 7 in <cit.>, we can obtain 

    E(w_j| w_1, …,w_k)=0,  j∈{k+1,…,p},

hence,

    E   {h(w_1,…,w_k) … w_j}=E{E ( h(w_1,…,w_k) · w_j| w_1,…, w_k ) }
       =E{ h(w_1,…,w_k) ·E(w_j| w_1, …,w_k) }=0,  j∈{k+1,…,p}.

Then,

    E   (YX)=C^⊤E{h(w_1,…,w_k)·W}
       =C^⊤[E{h(w_1,…,w_k)· w_1},…,E{h(w_1,…,w_k)· w_k},0,…,0]^⊤
    
    
    
    
    
    
    
    
       =∑_i=1^kβ_iE{h(w_1,…,w_k)· w_i}
    =∑_i=1^k_pβ_iE{h(w_1,…,w_k)· w_i}=Σ_XBΛ,

where Λ=(λ_1,…,λ_k)^⊤, 
λ_i=E{h(w_1,…,w_k)· w_i}.


Case 2: Σ_X=Σ. Let X^*=Σ^-1/2X, then X^* follows  elliptical distribution with mean 0 and variance-covariance matrix _p. According to Eq. (<ref>), we have

    E(Y|X)=E(Y|Σ^1/2X^*)=h(X^*^⊤Σ^1/2β_1,⋯,X^*^⊤Σ^1/2β_k).

Hence,

    E(YX)=Σ^1/2E(YX^*)
    =Σ^1/2Σ_X^*Σ^1/2BΛ= Σ_XBΛ.





 §.§ Proof of Theorem <ref>


Without loss of generality, we consider  E(X)=0.
 

Case 1:  Σ_X=_p. In this case, 

    E(Y^2|X)=H(β_1^⊤X,…,β_k^⊤X),

then,

    E(Y^2XX^⊤)   =E{E(Y^2XX^⊤|X) }
       =E{XX^⊤E(Y^2|X) }=E{XX^⊤ H(β_1^⊤X,…,β_k^⊤X) }.

Similarly, let C be the orthogonal matrix with the first k rows as β_i^⊤, i∈{1,…,k}, then we define 

    C=(β_1,…,β_k,α_k+1,…,α_p)^⊤≡(B, A)^⊤,

and hence,

    E(Y^2XX^⊤)   =E{H(β_1^⊤X, …,β_k^⊤X)·XX^⊤}
       =C^⊤E{H(β_1^⊤X,…,β_k^⊤X)·CXX^⊤C^⊤}C
       =C^⊤E{H(w_1,…,w_k)·WW^⊤}C.

Note that W=CX also follows elliptical distributions with with mean E(W)=0 and variance-covariance matrix CI_pC^⊤=_p. Let

    W=( w_1,…,w_k,w_k+1,…,w_p)^⊤=(W_(1)^⊤,W_(2)^⊤)^⊤.
 
By Corollary 8 in <cit.>, we can obtain

    E ( W_(2)W_(2)^⊤|W_(1) )
        =diag{E ( w_k+1^2|W_(1) ),…,E ( w_p^2|W_(1) )}I_p-k
       =E ( w_k+1^2|W_(1) )I_p-k

because {w_k+1,…,w_p} have the same status. Therefore,

    E(Y^2XX^⊤)   =C^⊤ E[H(w_1,…,w_k)·WW^⊤]C
       =C^⊤ E[ E [ H(w_1,…,w_k)·WW^⊤|W_(1)] ]C
       =C^⊤ E [ H(w_1,…,w_k)· E(WW^⊤|W_(1))  ]C
       =C^⊤[ Γ_1   0;   0 Γ_2 ]C
       =BΓ _1B^⊤ +AΓ _2A^⊤ ,

where 

    Γ_1=[ E{   H(w_1,…,w_k )w_1^2}                        ⋯ E{  H(w_1,…,w_k) w_1w_k};                        ⋮                        ⋱                        ⋮; E{  H(w_1,…,w_k) w_kw_1}                        ⋯  E{  H(w_1,…,w_k) w_k^2} ]_k× k,



    Γ_2   =[ E{E (H· w_k+1^2|W_(1) ) }                         ⋯                         0;                         ⋮                         ⋱                         ⋮;                         0                         ⋯   E{E (H· w_p^2|W_(1) ) } ]_(p-k)× (p-k)
       =E{E (H(w_1,…,w_k)· w_k+1^2|W_(1) )}I_p-k
       =E{ H(w_1,…,w_k)· w_k+1^2}I_p-k
       =E{E(Y^2|X)· w_k+1^2}I_p-k
       =E(Y^2· w_k+1^2)I_p-k=E(G)I_p-k.

Next, note that
E{Y^2(β^⊤X)^2}=β^⊤E(Y^2XX^⊤)β=β^⊤BΓ_1B^⊤β is the eigenvalue of Γ_1,
and E{Y^2(α^⊤X)^2} is the eigenvalue of Γ_2. Assumption <ref> assures that the eigenvalue of Γ_1 is larger, hence, the first k eigenvectors corresponding to the first k eigenvalues of E(Y^2XX^⊤) are the basis of S_E(Y|X).

Furthermore, E(Y^2XX^⊤) can be rewritten as

    E(Y^2XX^⊤)   =BΓ _1B^⊤+AΓ _2A^⊤
       =BΓ _1B^⊤ +E(G)AA^⊤  +E(G)BB^⊤ 
    -E(G)BB^⊤
       =B{Γ _1-E(G)_k}B^⊤ +E(G)_p,

and we can obtain

    E(Y^2XX^⊤)-E(G)_p=B{Γ _1-E(G)_k}B^⊤≡BM.

Because Γ _1-E(G)_k is a positive definite matrix by Assumption <ref>, 
rank(M)=k and rank(BM)=rank(B) according to the results in  Section A4.4 of <cit.>, i.e., E(Y^2XX^⊤)-E(G)_p is contained in the linear subspace Spanned by the basis matrix B.


Case 2: Σ_X=Σ. Let X^*=Σ^-1/2X, where X^* follows elliptical distribution with mean 0 and variance-covariance matrix _p, then

    E(Y^2|X)=E(Y^2|Σ^1/2X^*)=H(X^*^⊤Σ^1/2β_1,⋯,X^*^⊤Σ^1/2β_k),

and 

    E(Y^2XX^⊤ )   =Σ^1/2E(Y^2X^*X^*^⊤)Σ^1/2
       =Σ^1/2Σ^1/2BΓ _1B^⊤Σ^1/2Σ^1/2+Σ^1/2Σ^1/2AΓ _1A^⊤Σ^1/2Σ^1/2
       =ΣBΓ _1B^⊤Σ+ΣAΓ _1A^⊤Σ
       =Σ_XB{Γ _1-E(G)_k}B^⊤Σ_X+Σ_XE(G)_p.

We can obtain

    E(Y^2XX^⊤ )-E(G)_pΣ_X=Σ_XB{Γ _1-E(G)_k}B^⊤Σ_X.

Similarily, the eigenvalue of Γ_1 is larger and Γ _1-E(G)_k is a positive definite matrix under Assumption <ref>. Therefore, the first k eigenvectors corresponding to the first k eigenvalues of E(ZZ^⊤ ) are the basis of Span() and
Σ_X^-1E(Y^2XX^⊤ )Σ_X^-1-Σ_X^-1E(G)_p is contained in the linear subspace Spanned by the basis matrix B.

When X∼N_p(μ,Σ), 
without loss of generality, 
we can obtain   E( W_(2)W_(2)^⊤|W_(1)) = I_p-k by transforming X^*=Σ^-1/2X, then 

    E   [E{H(w_1,…,w_k)· w_k+1^2|W_(1)}]
       =E{ H(w_1,…,w_k) E ( w_k+1^2|W_(1) )}=E(Y^2).

Furthermore, let Ỹ^2=Y^2-E(Y^2), then E(Ỹ^2)=0 and

    Σ^-1E(Ỹ^2XX^⊤ )Σ^-1= BΓ _1B^⊤≡B.

This means that the non-zero k eigenvectors corresponding to the non-zero k eigenvalues of Σ^-1E(Ỹ^2XX^⊤ )Σ^-1 are the basis of Span(B). 



 §.§ Proof of Theorem <ref>


By the Law of Large Number, we have

    1n∑_i=1^nY_i(X_i-X̅)=1n∑_i=1^nY_i(X_i-μ)+Y̅(μ-X̅)Pr⟶E(YX)=Σ_XBΛ.

By Lemma <ref> and S_nPr⟶Σ_X, we have

    S_n^-11n∑_i=1^nY_i(X_i-X̅)Pr⟶Σ_X^-1Σ_XBΛ=BΛ.

Let

    Z_i=S_n^-1{Y_i(X_i-X̅)},

then,

    Ẑ=1n∑_i=1^nZ_iPr⟶BΛ

with convergence rate of n^1/2, and 

    𝒵̂=1n∑_i=1^nZ_iZ_i^⊤Pr⟶E(ZZ^⊤)

with convergence rate n^1/2. 

From Theorem <ref>, the first k eigenvectors corresponding to the first k eigenvalues of E(ZZ^⊤ ) are the basis of Span(B). Consequently, the first k eigenvectors corresponding to the first k eigenvalues of 𝒵̂, β̂_1, …,β̂_k, converge to the corresponding rotational basis for Span(B) with convergence rate of n^1/2.

Since the elements of vec(𝒵̂) are moment estimators of the elements of ZZ^⊤, by the central limit theorem,  we have 

    √(n) [ vec(𝒵̂)- vec{E(ZZ^⊤)} ]

converges in distribution to a multivariate normal random vector with mean vector 0 and variance-covariance matrix Var{vec(ZZ^⊤ ) }. Here, we derive the specific form of Var{vec(ZZ^⊤ ) }. Using the techniques in the proof of Theorem <ref> in Appendix A.2, we can obtain the fourth conditional moment of Y given X as

    E(Y^4|X)=U(β_1^⊤X,…,β_k^⊤X).

Given E(Y^2XX^⊤ ), the essence of obtaining the variance-covariance matrix is calculating the fourth moment. For notation convenience, we use the formation of the Kronecker product, which contains the fourth moments. Following the definition of W and the result in the proof of Theorem <ref>, we have 

    Var{vec(ZZ^⊤) }=Var{ Y^2vec(XX^⊤) }
       =E{ Y^2vec(XX^⊤)vec(XX^⊤)^⊤  Y^2}
                -E{ Y^2vec(XX^⊤)}×E{ Y^2vec(XX^⊤)}^⊤
       =E[ E{ Y^4vec(XX^⊤)vec(XX^⊤)^⊤|X}]
               -E[ E{ Y^2vec(XX^⊤)|X}]×E[ E{ Y^2vec(XX^⊤)|X}]^⊤
       ∝E{(XX^⊤)⊗(XX^⊤) U(β_1^⊤X,…,β_k^⊤X)  }
               -vec{E(Y^2XX^⊤)}×vec{E(Y^2XX^⊤)}^⊤ ,

where ∝ denotes the dimensional inequality and ⊗ denotes the Kronecker product. As the Kronecker product can incorporate all the elements of the variance-covariance matrix of the random matrix ZZ^⊤,  it does not affect the convergence of elements even though the dimensions are different. 

We can obtain E(Y^2XX^⊤) from the proof of Theorem <ref> in Appendix . Then,  we require the following: 

    E   { Y^4(XX^⊤) ⊗(XX^⊤)}
       =E{(XX^⊤) ⊗(XX^⊤) U(β_1^⊤X,…,β_k^⊤X)  }
       =E{(C^⊤WW^⊤C)⊗(C^⊤WW^⊤C) U(β_1^⊤X,…,β_k^⊤X)  }
       =E[ U(w_1,…,w_k)E{(C^⊤WW^⊤C)⊗(C^⊤WW^⊤C)|W_(1)}]
       =E{ U(w_1,…,w_k)(C^⊤⊗C^⊤)E(WW^⊤⊗WW^⊤|W_(1))(C⊗C) }
       =(C^⊤⊗C^⊤)E{ U(w_1,…,w_k)E(WW^⊤⊗WW^⊤|W_(1)) }(C⊗C).

As a result, we have E( w_iw_jw_sw_t|W_(1)), where i,j,s,t∈{1,…,p}. Considering Σ_X=_p, the form of the elements in E(WW^⊤⊗WW^⊤|W_(1)) can be expressed as 

    E ( w_iw_jw_sw_t|W_(1) )={
     w_iw_jw_sw_t,     i,j,s,t ∈{1,…,k},
     
    w_iw_jE ( w_sw_t|W_(1) ),     i,j ∈{1,…,k}, s=t ∈{k+1,…,p},
    
    w_iE ( w_jw_sw_t|W_(1) ),     i ∈{1,…,k},
    j=s=t ∈{k+1,…,p},
    E ( w_iw_jw_sw_t|W_(1) ),     i=j=s=t ∈{k+1,…,p},
     
    0,    otherwise.
    .

Hence, E{ Y^4(XX^⊤) ⊗(XX^⊤)}, and the same technique can be applied when Σ_X=Σ.



 §.§ R code for the proposed PSRFR




PSRFR <- function(X, y, r) {

  n <- nrow(X)          ## Number of observations
  
  dim <- ncol(X)        ## Dimensionality of X
  
  ## Create data matrix combining X and y
  
  data <- cbind(X, t(y))
  
  ## Mean calculations
  
  my <- mean(y)
  
  mx <- colMeans(X)     ## Vector of means for each column in X
  
  ## Centering matrix X
  
  mxx <- matrix(mx, nrow = n, ncol = dim, byrow = TRUE)
  
  z <- (X - mxx) * data[, dim + 1]  ## Computing z
  
  ## Covariance of X and related calculations
  
  sigmax <- Cov(X)
  
  invSigmax <- solve(sigmax)
  
  K <- invSigmax %*% (t(z) %*% z / n) %*% invSigmax
  
 ## Eigenvalue decomposition
  
  Kv <- eigen(K)
  
  ## Return a list of eigenvalues and eigenvectors
  
  return(list(values = Kv$values, vectors = Kv$vectors[,1:r]))
  
}














apalike

