orthoDr: semiparametric dimension reduction via orthogonality constrained optimization

orthoDr is a package in R that solves dimension reduction problems using a orthogonality constrained optimization approach. The package serves as a unified framework for many regression and survival analysis dimension reduction models that utilize semiparametric estimating equations. The main computational machinery of orthoDr is a first-order algorithm developed by Wen&Yin (2013) for optimization within the Stiefel manifold. We implement the algorithm through Rcpp and OpenMP for fast computation. In addition, we developed a general-purpose solver for such constrained problems with user-specified objective functions, which works as a drop-in version of optim(). The package also serves as a platform for future methodology developments along this line of work.


Introduction
Dimension reduction is a long-standing problem in statistics and data science. While the traditional principal component analysis (Jolliffe, 1986) and related works provide a way of reducing the dimension of the covariates, the term "sufficient dimension reduction" is more commonly referring to a series of regression works originated from the seminal paper on sliced inverse regression (Li, 1991). In such problems, we observe an outcome Y ∈ R, along with a set of covariates X = (X 1 , . . . , X p ) T ∈ R p . Dimension reduction models are interested in modeling the conditional distribution of Y given X, while their relationship satisfies, for some p × d matrix B = (β 1 , . . . , β p ), where represents any error terms and h, with a slight abuse of notation, represents the link function using X or B T X. One can easily notice that when d, the number of columns in B, is less than p, a dimension reduction is achieved, in the sense that only a d dimensional covariate information is necessary for fully describing the relationship (Cook, 2009). Alternatively, this relationship can be represented as (Zeng and Zhu, 2010) which again describes the sufficiency of B T X. Following the work of Li (1991), a variety of methods have been proposed. An incomplete list of literature includes Cook and Weisberg (1991); Cook and Lee (1999); Yin and Cook (2002); Chiaromonte et al. (2002); Zhu et al. (2006); Li and Wang (2007); Zhu et al. (2010b,a); Cook et al. (2010); Lee et al. (2013); Cook and Zhang (2014); Li and Zhang (2017). For a more comprehensive review of the literature, we refer the readers to Ma and Zhu (2013b). One advantage of many early developments in dimension reduction models is that only a singular value decomposition is required to obtain the reduced space parameters B through inverse sliced averaging. However, this comes at a price of assuming the linearity assumption (Li, 1991), which is almost the same as assuming that the covariates follow an elliptical distribution (Li and Dong, 2009;Dong and Li, 2010). Moreover, some methods require more restrictive assumptions on the covariance structure (Cook and Weisberg, 1991). Many methods attempt to avoid these assumptions by resorting to nonparametric estimations. The most successful ones include Xia et al. (2002) and Xia (2007). However, recently a new line of work started by Zhu (2012b,a, 2013a) shows that by formulating the problem into semiparametric estimating equations, not only we can avoid many distributional assumptions on the covariates, the obtained estimator of B also enjoys efficiency. Extending this idea, Sun et al. (2017) developed a framework for dimension reduction in survival analysis using a counting process based estimating equations. The method performs significantly better than existing dimension reduction methods for censored data such as Li et al. (1999); Xia et al. (2010) and Lu and Li (2011). Another recent development that also utilizes this semiparametric formulation is Zhao et al. (2017), in which an efficient estimator is derived.
Although there are celebrated theoretical and methodological advances, estimating B through the semiparametric estimating equations is still not a trivial task. Two challenges remain: first, by a careful look at the model definition 1, we quickly noticed that the parameters are not identifiable unless certain constraints are placed. In fact, if we let A be any d × d full rank matrix, then (BA) T X preserves the same column space information of B T X, hence, we can define h * ((BA) T X, ) accordingly to retain exactly the same model as (1). While traditional methods can utilize singular value decompositions (SVD) of the estimation matrix to identify the column space of B instead of recovering each parameter (Cook and Lee, 1999), it appears to be a difficult task in the semiparametric estimating equation framework. One challenge is that if we let B change freely, the rank of the B matrix cannot be guaranteed, which makes the formulation meaningless. Hence, for both computational and theoretical concerns, Ma and Zhu (2012b) resorts to an approach that fixes the upper d × d block of B as an identity matrix, i.e., B = (I d×d , B * T ) T , where B * is a (p − d) × d matrix that sits in the lower block of B. Hence, in this formulation, only B * needs to be solved. While the solution is guaranteed to be rank d in this formulation, as pointed out by Sun et al. (2017), this approach still requires correctly identifying and reordering of the covariate vector x such that the first d entries are indeed important, which creates another daunting task. Another challenge is that solving semiparametric estimating equations requires the estimation of nonparametric components. These components need to be computed through kernel estimations, usually the Nadaraya-Watson type, which significantly increases the computational intensity of the method considering that these components need to be recalculated at each iteration of the optimization. Up to date, these drawbacks remain as the strongest criticism of the semiparametric approaches. Hence, although enjoying superior statistical asymptotic properties, are not as attractive as a traditional sliced inverse type of approaches such as Li (1991) and Cook and Weisberg (1991).
The goal of our orthoDr package is to develop a computationally efficient optimization platform for solving the semiparametric estimating equation approaches proposed in Ma and Zhu (2013a), Sun et al. (2017) and possibly any future work along this line. Revisiting the rank preserving problem of B mentioned above, we can essentially set a constraint that where I is a d × d identity matrix. A solution of the estimating equations that satisfies the constraint will correctly identify the dimensionality-reduced subspace. This is known as optimizing on the Stiefel manifold, which is a class of well-studied problems (Edelman et al., 1998). A recent R development (Martin et al., 2016) utilizes quasi-Newton methods such as the well known BFGS method on the Riemannian manifold (Huang et al., 2018). However, second order optimization methods always require forming and storing large hessian matrices. In addition, they may not be easily adapted to penalized optimization problems, which often appear in high dimensional statistical problems Zhu et al. (2006); Li and Yin (2008). On the other hand, first-order optimization methods are faster in each iteration, and may also incorporate penalization in a more convenient way Wen et al. (2010). By utilizing the techniques developed by Wen and Yin (2012), we can effectively search for the solution in the Stiefel manifold, and this becomes the main machinery of our package. Further incorporating the popular Rcpp (Eddelbuettel and François, 2011) and RcppArmadillo (Eddelbuettel and Sanderson, 2014) toolboxes and the OpenMP parallel commuting, the computational time for our package is comparable to state-of-the-art existing implementations (such as ManifoldOpthm), making the semiparametric dimension reduction models more accessible in practice.
The purpose of this article is to provide a general overview of the orthoDr package (version 0.6.2) and provide some concrete examples to demonstrate its advantages. orthoDr is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=orthoDr and GitHub at https://github.com/teazrq/orthoDr. We begin by explaining the underlying formulation of the estimating equation problem and the parameter updating scheme that preserves orthogonality. Next, the software is introduced in detail using simulated data and real data as examples. We further demonstrate an example that utilizes the package as a general purpose solver. We also investigate the computational time of the package compared with existing solvers. Future plans for extending the package to other dimension reduction problems are also discussed.

Model description Counting process based dimension reduction
To give a concrete example of the estimating equations, we use the semiparametric inverse regression approach defined in Sun et al. (2017) to demonstrate the calculation. Following the common notations in the survival analysis literature, let X i be the observed p dimensional covariate values of subject i, Y i = min(T i , C i ) is the observed survival time, with failure time T i and censoring time C i , and is observed. We are interested in a situation that the conditional distribution of failure time T i |X i depends only on the reduced space B T X i . Hence, to estimate B, the estimating equation is given by where the operator vec(·) is the vectorization of matrix. Several components are estimated nonparasitically: the function ϕ(u) is estimated by sliced averaging, where u is chosen such that there are hn number of observations lie between u and u + u. The conditional mean function E X|Y ≥ u, B T X = z is estimated through the Nadaraya-Watson kernel estimator In addition, the the conditional hazard function at any time point u can be estimated by However, this substantially increase the computational burden since the double kernel estimator requires O(n 2 ) flops to calculate the hazard at any given u and z. Instead, an alternative version using Dabrowska (1989) can greatly reduce the computational cost without compromising the performance. Hence, we estimate the conditional hazard function by which requires only O(n) flops. In the above equations (5), (6) and (8), h is a pre-specified kernel bandwidth and K h (·) = K(·/h)/h, where K(·) is the Gaussian kernel function. By utilizing the method of moments estimators (Hansen, 1982) and noticing our constraint for identifying the column space of B, solving for the solution of the estimating equations (4) is equivalent to Essentially all other semiparametric dimension reduction models described in Ma and Zhu (2013a), and more recently Ma and Zhang (2015) (2017) and many others can be estimated in the samimilar fashion as the above optimization problem. However, due to the difficult in the constrains and the purpose of identifiability, all of these methods resort to either fixing the upper block of the B matrix as an identity matrix or adding a penalty of B T B − I F to preserve the orthogonality constraint. There appears to be no existing method that solves (9) directly. Here, we utilize Wen and Yin (2012)'s approach which can effectively tackle this problem.

Orthogonality preserving updating scheme
The algorithm works in the same fashion as a regular gradient decent, except that we need to preserve the orthogonality at each iteration of the update. As described in Wen and Yin (2012), given any feasible point B 0 , i.e., B 0 T B 0 = I, which can always be generated randomly, we update B 0 as follows. Let the p × d gradient matrix be Then, utilizing the Cayley transformation, we have with the orthogonality preserving property B T new B new = I. Here, A = GB 0 T − B 0 G T is a skewsymmetric matrix. It can be shown that {B new (τ)} τ≥0 is a descent path. Similar to line search algorithms, we can then find a proper step size τ through a curvilinear search. Recursively updating the current value of B, the algorithm stops when the tolerance level is reached. An initial value is also important for the performance of nonconvex optimization problems. A convenient initial value for our framework is the computational efficient approach developed in Sun et al. (2017), which only requires a SVD of the estimation matrix.

The R package orthoDr
There are several main functions in the orthoDr package: orthoDr_surv, ortho_reg and ortho_optim. They are corresponding to the survival model described perviously (Sun et al., 2017), the regression model in Ma and Zhu (2012b), and a general constrained optimization function, respectively. In this section, we demonstrate the details of using these main functions, illustrate them with examples.

Semiparametric dimension reduction models for survival data
The orthoDr_surv function implements the optimization problem defined in Equation (9), where the kernel estimations and various quantities are implemented and calculated within C++. Note that in addition, the method defined previously, some simplified versions are also implemented such as the counting process inverse regression models and the forward regression models, which are all described in Sun et al. (2017). These specifications can be made using the method parameter. A routine call of the function orthoDr_surv proceed as orthoDr_surv(x, y, censor, method, ndr, B.initial, bw, keep.data, control, maxitr, verbose, ncore) • x: A matrix or data.frame for features (numerical only).
• y: A vector of observed survival times.
• censor: A vector of censoring indicators.
• method: The estimating equation method used.
-"forward": forward regression model with one structural dimensional.
• ndr: The number of structural dimensional. For method = "dn" or "dm", the default is 2. For method = "forward" only one structural dimension is allowed, hence the parameter is suppressed.
• B.initial: Initial B values. Unless specifically interested, this should be left as default, which uses the computational efficient approach (with the CPSIR() function) in Sun et al. (2017) as the initial. If specified, must be a matrix with ncol(x) rows and ndr columns. The matrix will be processed by Gram-Schmidt if it does not satisfy the orthogonality constrain.
• bw: A kernel bandwidth, assuming each variables have unit variance. By default we use the Silverman rule-of-thumb formula Silverman (1986) to determine the bandwidth This bandwidth can be computed using the silverman(n,d) function in our package.
• keep.data: Should the original data be kept for prediction? Default is FALSE.
• control: A list of tuning variables for optimization, including the convergence criteria. In particular, epsilon is the size for numerically approximating the gradient, ftol, gtol, and btol are tolerance levels for the objective function, gradients, and the parameter estimations, respectively, for judging the convergence. The default values are selected based on Wen and Yin (2012) .
• verbose: Should information be displayed? Default is FALSE.
• ncore: Number of cores for parallel computing when approximating the gradients numerically. The default is the maximum number of threads.
We demonstrate the usage of orthoDr_surv function by solving a problem with generated survival data.
• s2: A matrix for the second column space (e.g., B).
• method: -"dist": the Frobenius norm distance between the projection matrices of the two given matrices, where for any given matrix B, the projection matrix P = B(B T B) −1 B T .
-"trace": the trace correlation between two projection matrices tr(P P)/d, where d is the number of columns of the given matrix.
-"canonical": the canonical correlation between B T X and B T X.
• x: The design matrix X (default = NULL), required only if method = "canonical" is used.
We compare the accuracy of the estimations obtained by the method ="dm" and "dn". Note that the "dm" method enjoys double robustness property of the estimating equations, hence the result is usually better.

Semiparametric dimension reduction models for regression
The orthoDr_reg function implements the semiparametric dimension reduction methods proposed in Ma and Zhu (2012b). A routine call of the function orthoDr_reg proceed as orthoDr_reg (x, y, method, ndr, B.initial, bw, keep.data, control, maxitr, verbose, ncore) • x: A matrix or data.frame for features (numerical only).
• y: A vector of observed continuous outcome.
• method: We currently implemented two methods: the semiparametric sliced inverse regression method ("sir"), and the semiparametric principal Hessian directions method ("phd").
-"sir": semiparametric sliced inverse regression method solves the sample version of the estimating equation -"phd": semiparametric principal Hessian directions method that estimates B by solving the sample version of • ndr: The number of structural dimensional (default is 2).
• B.initial: Initial B values. For each method, the initial values are taken from the corresponding traditional inverse regression approach using the dr package. The obtained matrix will be processed by Gram-Schmidt for orthogonality.
• bw, keep.data, control, maxitr, verbose and ncore are exactly the same as those in the orthoDr_surv function.
To demonstrate the usage of orthoDr_reg, we consider the problem of dimension reduction by fitting a semi-PHD model proposed by Ma and Zhu (2012b).

Parallelled gradient approximation through OpenMP
The estimation equations of the dimension reduction problem in the survival and regression settings usually have a complicated form. Especially, multiple kernel estimations are involved, which results in difficulties in taking derivatives analytically. As an alternative, numerically approximated gradients are implemented using OpenMP. A comparison between a single core and multiple cores (4 cores) is given in the following example. Results from 20 independent simulation runes are summarized in Table 1. The data generating procedure used in this example is the same as the survival data used in Section 2.3.1. All simulations are performed on an i7-4770K CPU. R> t0 = Sys.time() R> dn.fit = orthoDr_surv(dataX, Y, Censor, method = "dn", ndr = ndr, + ncore = 4, control = list(ftol = 1e-6)) R> Sys.time() -t0 General solver for orthogonality constrained optimization ortho_optim is a general purpose optimization function that can incorporate any user defined objective function f (and gradient function if supplied). The usage of ortho_optim is similar to the widely used optim() function. A routine call of the function proceed as ortho_optim(B, fn, grad, ..., maximize, control, maxitr, verbose) The R Journal Vol. 11/2, December 2019 ISSN 2073-4859 • B: Initial B values. Must be a matrix, and the columns are subject to the orthogonality constrains. It will be processed by Gram-Schmidt if not orthogonal.
• fn: A function that calculates the objective function value. The first argument should be B. Returns a single value.
• grad: A function that calculate the gradient. The first argument should be B. Returns a matrix with the same dimension as B. If not specified, a numerical approximation is used.
• maximize: By default, the solver will try to minimize the objective function unless maximize = TRUE.
• The parameters maxitr, verbose and ncore works in the same way as introduced in the previous sections.
To demonstrate the simple usage of ortho_optim as a drop-in function of optim(), we consider the problem of searching for the first principle component for a data matrix.
We found that "LRBFGS" and our orthoDr package usually achieve the best performance, with functional value decreases the steepest in the log scale. In terms of computing time, "LRBFGS" and orthoDr performers similarly. Although "LRTRSR1" has similar computational time, its functional value falls behind. This is mainly because the theoretical complexity of second-order algorithms is similar to first order algorithms, both are of order O(p 3 ). However, it should be noted that for a semiparametric dimension reduction method, the major computational cost is not due to the parameter updates, rather, it is calculating the gradient since complicated kernel estimations are involved. Hence, we believe there is no significant advantage using either "LRBFGS" or our orthoDr package regarding the efficiency of the algorithm. However, first order algorithms may have an advantage when developing methods for penalized high-dimensional models.

Examples
We use the Concrete Compressive Strength (Yeh, 1998) dataset as an example to further demonstrate the orthoDr_reg function and to visualize the results. The dataset is obtained from the UCI Machine Learning Repository.
Concrete is the most important material in civil engineering. The concrete compressive strength is a highly nonlinear function of age and ingredients. These ingredients include cement, blast furnace slag, fly ash, water, superplasticizer, coarse aggregate, and fine aggregate. In this dataset, we have n = 1030 observation, 8 quantitative input variables, and 1 quantitative output variable. We present the estimated two directions for structural dimension and further plot the observed data in these two directions. A non-parametric kernel estimation surface is further included to approximate the mean concrete strength.

Discussion
Using the algorithm proposed by Wen and Yin (2012) for optimization on the Stiefel manifold, we developed the orthoDr package that serves specifically for semi-parametric dimension reductions problems. A variety of dimension reduction models are implemented for censored survival outcome and regression problems. In addition, we implemented parallel computing for numerically appropriate the gradient function. This is particularly useful for semi-parametric estimating equation methods because the objective function usually involves kernel estimations and the gradients are difficult to calculate. Our package can also be used as a general purpose solver and is comparable with existing manifold optimization approaches. However, since the performances of different optimization approaches could be problem dependent, hence, it could be interesting to investigate other choices such as the "LRBFGS" approach in the ManifoldOptim package.
Our package also serves as a platform for future methodology developments along this line of work. For example, we are currently developing a personalized dose-finding model with dimension reduction structure (Zhou and Zhu, 2018). Also, when the number of covariates p is large, the model can be over-parameterized. Hence, applying a L 1 penalty can force sparsity and allow the model to handle high-dimensional data. To this end, first-order optimization approaches can have advantages over second-order approaches. However, persevering the orthogonality during the Cayley transformation while also preserve the sparsity can be a challenging task and requires new methodologies. Furthermore, tuning parameters can be selected through a cross-validation approach, which can be implemented in the future.