TULIP: A Toolbox for Linear Discriminant Analysis with Penalties

Linear discriminant analysis (LDA) is a powerful tool in building classifiers with easy computation and interpretation. Recent advancements in science technology have led to the popularity of datasets with high dimensions, high orders and complicated structure. Such datasetes motivate the generalization of LDA in various research directions. The R package TULIP integrates several popular high-dimensional LDA-based methods and provides a comprehensive and user-friendly toolbox for linear, semi-parametric and tensor-variate classification. Functions are included for model fitting, cross validation and prediction. In addition, motivated by datasets with diverse sources of predictors, we further include functions for covariate adjustment. Our package is carefully tailored for low storage and high computation efficiency. Moreover, our package is the first R package for many of these methods, providing great convenience to researchers in this area.


Introduction
Linear discriminant analysis (LDA) is one of the most popular classification method and a cornerstone for multivariate statistics (Michie et al., 1994, e.g). Classical LDA builds a linear classifier based on p-dimensional multivariate predictor X ∈ R p to distinguish K classes and to predict the class label Y ∈ {1, . . . , K}. Despite its simplicity, LDA is shown to be very accurate on many benchmark datasets (Lim et al., 2000;Dettling, 2004;Hand, 2006). Moreover, LDA is easily interpretable and is thus often used as a visualizing tool for exploratory data analysis.
In recent decades, the advancements in science and technology have enabled researchers to collect datasets with increasing sizes and complexity. Such datasets pose challenges to LDA. Four challenges that we tackle with this package are as follows. First, in research areas such as biology, genomics and psychology, we often have more predictors than samples. However, LDA is not applicable on these high-dimensional data, because sample covariance matrix becomes not invertible when the number of predictors exceeds the sample size.
Secondly, when we have a large number of predictors, variable selection is often desired such that we can obtain a sparse classifier involving only a small proportion of the variables. On one hand, Fan and Fan (2008); Bickel and Levina (2008) showed in theory that variable selection is critical for accurate classification. On the other hand, sparse classifiers much easier to interpret in practice. However, LDA generally does not perform variable selection.
Thirdly, contemporary datasets often have complicated structure that renders the linear classifier in LDA inadequate. For example, in the presence of thousands of predictors, it may be inappropriate them to model all of them with the normal distribution. Moreover, research in neuroimaging, computational biology and personalized recommendation produces data in the form of matrices (2-way tensor) or tensors. The analysis of tensor datasets requires considerable modification to the vector-based LDA model.
Last but not least, integrative analysis with multiple data sources are drawing researchers' attention recently. Co-existence of diverse data types, such as vector, matrix and tensor calls for more sophisticated models to integrate the information from them to improve classification accuracy. It is critical to model the dependence among different types of data to reduce the noise level in the data and improve prediction accuracy (Pan et al., 2019).
Motivated by these challenges, many methods have been proposed to generalize LDA to datasets with high dimensions, non-normality and/or higher order predictors. In this package, we implement six methods that generalize LDA to contemporary complicated datasets. All of them are developed under models closely related to the LDA model, and penalties are imposed to achieve classification accuracy and variable selection in high dimensions. These methods include: 1. Direct sparse discriminant analysis (DSDA): DSDA generalizes the classical LDA model to high dimensions when there are only two classes (Mai et al., 2012). It formulates high-dimensional LDA into a penalized least squares problem.
2. Regularized optimal affine discriminant (ROAD): under the same model as DSDA, ROAD fits a sparse classifier by minimizing the classification error under the 1 constraint (Fan et al., 2012).
3. Sparse optimal scoring (SOS) for binary problems: SOS is also developed under the LDA model (Clemmensen et al., 2011). It penalizes the optimal scoring problem (Hastie et al., 1994). We focus on its application in binary problems.  Table 1: Comparison of model settings between models. SOS was originally proposed to deal with both binary and multiclass problems, but we focus on binary problems in the package. Model SeLDA stands for Semi-parametric linear discriminant analysis, which is introduced in Section 2.2.4. Model TDA/CATCH represents tensor discriminant analysis and covariate-adjusted tensor in highdimensions, which are illustrated in Section 2.2.5 and 2.3.6.
4. Semiparametric sparse discriminant analysis (SeSDA): SeSDA assumes a semiparametric model where data transformation can be applied to alleviate the non-normality. In practice, SeSDA finds the data-driven transformation and then performs model-fitting on the transformed data .
5. Multiclass sparse discriminant analysis (MSDA): Instead of focusing on binary problems, MSDA considers the multiclass LDA model . It takes note of the fact that the Bayes' rule can be estimated with minimizing a quadratic loss. To account for the multiclass structure, a group lasso penalty (Yuan and Lin, 2006) is applied to achieve variable selection.
6. Covariate-adjusted tensor classification in high-dimensions (CATCH): CATCH (Pan et al., 2019) is developed for tensor predictors. It takes advantage of the tensor structure to significantly reduce the number of parameters and hence alleviate computation complexity.
See Table 1 for a comparison of these methods. Despite their different model assumptions and formulas, all of them have strong theoretical support and excellent empirical performance. We further note that they can be combined with covariate adjustment when multiple data sources are available. Our package TULIP (Pan et al., 2021) integrates diverse discriminant analysis models and supportive functions to make it a convenient and well-equipped toolbox. It has several notable advantages. First, we not only include functions for model fitting, but also cross validation functions for easy control of the sparsity level, and prediction functions for the prediction of future observations. In addition, we provide covariate adjustment functions that efficiently remove heterogeneity in the predictors and combine information from covariates. Second, our package greatly facilitates the application of DSDA, ROAD and SeSDA for R users, as they do not have public R packages on CRAN outside ours. Third, although MSDA and SOS have been implemented in packages msda  and sparseLDA (Clemmensen and Kuhn, 2016), we carefully modify their algorithms in our implementation to lower storage cost and/or speed up computation.
We acknowledge that many other efforts have been spent on topics closely related to that of our paper. On one hand, by now a large number of high-dimensional discriminant analysis methods have been developed. Some excellent examples include Fan and Fan (2008); Tibshirani et al. (2002); Trendafilov and Jolliffe (2007);Fan et al. (2012); Wu et al. (2009);Cai et al. (2011);Shao et al. (2011);Clemmensen et al. (2011);Witten and Tibshirani (2011);Xu et al. (2015); Niu et al. (2015). On the other hand, in the literature, many works study matrix/tensor regression and classification methods. Many of them impose low rank assumption (Zhou et al., 2013;Kolda and Bader, 2009;Chi and Kolda, 2012;Liu et al., 2017;Li and Schonfeld, 2014;Lai et al., 2013;Zhong and Suslick, 2015;Zeng et al., 2015). All these methods have been reported to have great performance, but a comprehensive study of them is apparently out of the scope of our current paper.
The rest of this paper is organized as follows. We start with a brief overview of discriminant analysis models in Section 2.2. Model estimation and implementation details are discussed in Section 2.3. Section 2.4 contains instructions and examples on the usage of the package. A real data example is given in Section 2.5 to confirm the numerical performance of methods in the package.

Bayes rule for classification
Recall that Y ∈ {1, . . . , K} is the categorical response (class indicator), and we use the generic X to denote the predictor and (potential) additional covariate. Specifically, X = X ∈ R p in classical multivariate discriminant analysis; X = X ∈ R p 1 ×···×p M in tensor discriminant analysis; and X = (X, U) in covariate-adjusted classification settings, where U ∈ R q is additional covariates and X can be either vector or tensor. Our goal is to construct the optimal classifier to distinguish and predict Y based on X under various settings. Denote π k = Pr(Y = k) and f k as the conditional distribution of X within Class k (e.g. f k (X ) = f (X, U | Y = k) is the joint distribution of X and U given Y = k, in presence of U). The optimal classifier, often referred to as the Bayes rule, is thus (1) The Bayes rule achieves the lowest classification error possible (Friedman et al., 2001). Therefore, it is our ultimate goal to estimate the Bayes rule. However, additional model assumptions are often needed for f k to ensure statistical and computational efficiency. Consider the classical LDA setting of X ∈ R p and Y ∈ {1, . . . , K}. To gain intuition, we often assume that within each class, the predictor follows a normal distribution with different means and a common covariance matrix. Then the Bayes rule is a linear function of X and can be straightforwardly estimated.
In the rest of this section, we discuss various statistical models that have been widely studied in the literature, along with the Bayes rules under these assumptions. Specifically, we review the classical LDA model, the semiparametric LDA model, and the tensor discriminant analysis model. We also discuss a general framework for covariate adjustment.

The linear discriminant analysis model (LDA)
Given a multivariate predictor X ∈ R p and Y ∈ {1, . . . , K}, the LDA model assumes that X is normally distributed within each class, i.e, where µ k ∈ R p is the mean of X within class k, and Σ ∈ R p×p is the common within class covariance matrix.
Define β k = Σ −1 (µ k − µ 1 ) for k = 1, · · · , K. The Bayes' rule turns out to be a linear function: The LDA model is simple yet elegant. All the parameters in this model have natural interpretations, while the Bayes rule has a nice linear form. An interesting fact about the Bayes rule in (3) is that it does not explicitly involve the p 2 -dimensional parameter Σ −1 . Instead, Σ −1 is only implicitly included in the discriminant directions β k . Moreover, it can be shown that the Bayes rule is equivalent to first reducing data to X T β 2 , . . . X T β K and then fitting the LDA model on the (K − 1)-dimensional space. Therefore, to estimate the Bayes rule in high dimensions, our interest centers on the estimation of β k . We assume that β k 's are sparse with many elements being zero. Enforcement of this sparsity assumption will facilitate our estimation and naturally lead to variable selection.
Although the Bayes rule is derived under the somewhat restrictive normality and equal covariance assumptions, the discriminant directions β k are still meaningful when data are non-normal, thanks to their geometric properties. It can be shown that, if we project X to β k , k = 1, . . . , K, the separation between classes is maximized over all possible sets of K − 1 linear projections. Consequently, the LDA model is reasonably resistant to model misspecification. However, in some of the cases where the LDA model assumptions are severely violated, one can resort to more flexible models. For example, the quadratic discriminant analysis model (Jiang et al., 2015;Fan et al., 2015;Li and Shao, 2015;Sun and Zhao, 2015) relaxes the equal covariance assumption, while severely non-normal data can be modeled by the semiparametric model to be discussed in Section 2.2.4.

Covariates adjustment
In many real-life problems, we have additional covariates along with the predictors. The covariates play two roles in the classification: it has predictive power on it own, and it also accounts part of the variation in the predictors. For example, in genomics studies, we record not only gene expression levels but also age and clinical measurements. In this case, we may view the gene expression levels as the high-dimensional predictor, and the age and clinical measurements as the covariates. We consider an LDA-type model to incorporate the covariates. In addition to the response Y and the predictor X, we denote the covariates as U ∈ R q . We assume that where φ k ∈ R q is the within-class mean, Ψ ∈ R q×q , Ψ > 0 is the common within class covariance matrix of covariates, and α ∈ R p×q is the dependence of X on U. We refer to this model as the covariate-adjusted LDA (CA-LDA) model. The CA-LDA model is conceptually similar to the CATCH model (Pan et al., 2019) for tensor, which is to be introduced in Section 2.2.5, but the CA-LDA model focuses on vector predictor X rather than tensor predictor.
Obviously, the CA-LDA model reduces to the LDA model in the absence of covariates. With the covariates, the CA-LDA model continues to have natural interpretations. Equation (4) indicates that (U, Y) marginally follow the LDA model. Equation (5) implies that the distribution of X not only depends on Y, but also U through mean dependence. Therefore, within each class, X is linked to U through a linear regression model, while, after we adjust for U, (X, Y) follow the LDA model as well. See Figure 1 for a graphical illustration of the relationship among X, U and Y.
Under the CA-LDA model, the Bayes' rule is where γ k = Ψ −1 (φ k − φ 1 ), β k = Σ −1 (µ k − µ 1 ) and a k = log(π k /π 1 ) − 1 2 γ T k (φ k + φ 1 ) − 1 2 β T k (µ k + µ 1 )) is a scalar that does not involve X or U. Throughout this paper, we assume that U is lowdimensional and does not need variable selection, but X is high-dimensional. In the presence of covariates, X needs to be first adjusted to X − αU before entering the Bayes rule. Similar to the LDA model, we assume that the coefficient of X − αU, β k , is sparse.

The semiparametric LDA model
Although LDA is reasonably resistant to model misspecification, we may still need more flexible models when data are heavily non-normal. The semiparametric linear discriminant analysis (SeLDA) model (Lin and Jeon, 2003) is proposed for this purpose. SeLDA assumes that there exists a set of strictly monotone univariate transformations h 1 , . . . , h p such that For identifiability, we further assume that all the diagonal elements in Σ are 1, and all elements in µ 1 are 0. We also use the shorthand notation h(X) = (h 1 (X 1 ), · · · , h p (X p )). The transformation h is assumed to be unknown and needs to be estimated from data. The SeLDA model assumes that the LDA model is true up to an unknown transformation. It has the same spirit as the well-known Box-Cox transformation, with which model assumptions are relaxed by proper data mapping.
It is easy to see that the LDA model is a special case of the SeLDA model, if we restrict h(X) = X. However, in the SeLDA model, we do not impose any parametric assumptions on h, which leads to great flexibility in practice. We further review a formula for h j that will facilitate its estimation. It can be shown that where Φ is the cumulative distribution function (CDF) of the standard normal random variable, and F kj is the CDF of X j within Class k. Equation (8) will be used in Section 2.3.4. The SeLDA model also amounts to assuming that the data follow the Gaussian copula model within each class (Klaassen and Wellner, 1997;Hoff et al., 2014;Liu et al., 2009).
Although the SeLDA model requires much weaker conditions than the LDA model, it preserves many of the desirable properties. One of them is that the Bayes rule continues to be a linear function of the transformed data h(X): Consequently, just as in the LDA model, when the dimension is high, we assume that β k is sparse to allow accurate estimation.

Tensor discriminant analysis (TDA) and covariate adjustment
The tensor discriminant analysis (TDA) model is proposed for classification based on tensor predictors. We first briefly introduce some standard tensor notation (Kolda and Bader, 2009). See Appendenx A for more rigorous definitions. An M-way tensor is denoted by a multidimensional array A ∈ R p 1 ×···×p M where M ≥ 2, p 1 , . . . , p M are all positive integers. We often need to multiply an M-way tensor C by M matrices along each mode G i , i = 1, . . . , M, denoted by C; G 1 , . . . , G M . For example, in Figure 2 we obtain A = C; G 1 , . . . , G 3 by multiplying a 3-way tensor C with matrices G i along each mode. If G i , i = m are identity matrices and G m is a vector, then we write C× m G m = C; I, . . . , G m , . . . , I .
Further, we say a tensor X ∈ R p 1 ×···×p M follows the tensor normal distribution TN(µ, Σ 1 , . . . , Σ M ) if it can be written as has elements all independently standard normal, µ ∈ R p 1 ×···×p M is the mean tensor, and Σ m ∈ R p m ×p m are covariance matrices. See Figure 3 for an illustration. Figure 2: Tucker decomposition of tensor A. Figure 3: Tensor normal distribution. Now we discuss the tensor discriminant analysis (TDA) model. Consider the M-way tensor predictor X ∈ R p 1 ×···×p M where M ≥ 2 and class label Y ∈ {1, . . . , K}. The TDA model assumes that where µ k ∈ R p 1 ×···p M , Σ m ∈ R p m ×p m is the within-class mean, Σ m > 0 is the common within-class covariance matrix along the m-th mode of the tensor, and 0 < π k < 1 is the prior probability for Class k. Compared to the LDA model, TDA utilizes the tensor normal distribution to model X within each class. By taking advantage of the tensor structure, TDA drastically reduces the number of unknown parameters (Pan et al., 2019). It can be seen that the TDA model requires O(∑ M m=1 p 2 m ) parameters to model the dependence among X. However, if we ignore the tensor structure and assume the LDA model on the vectorized version of X, the covariance matrix has O(∏ M m=1 p 2 m ) parameters. Under the TDA model, the Bayes' rule is is a scalar that does not involve X. It can be seen that the Bayes rule is again a linear function in X, with the linear coefficients B k . In high dimensions, we again impose the sparsity assumption by assuming that many elements in B k are zeros.
Similar to the vector case, when additional covariates are provided, the TDA model can be combined with covariate adjustment. Pan et al. (2019) proposed the CATCH model for this purpose.
In addition to (Y, X), we are given the covariates U ∈ R q . The CATCH model assumes that where φ k ∈ R q is the within-class mean of U, Ψ ∈ R q×q is the within-class covariance of U, and α ∈ R p 1 ×···×p M ×q characterizes the dependence of X on U. The parameters in the CATCH model can be interpreted in the same way as the CA-LDA model in Section 2.2.3.
The Bayes' rule under the CATCH model is where is a scalar that does not involve X or U. Similar to the TDA model, we assume that B k is sparse in high dimensions, but impose no further sparsity assumptions on other parameters.

Methods
In this section, we formally introduce the six methods implemented by the package: DSDA, ROAD, SOS, SeSDA, MSDA and CATCH. Throughout the rest of this paper, we denote Σ as the pooled sample covariance, µ k as the within-class sample mean, n as the sample size, and n k as the sample size in class k. All the methods involve a tuning parameter λ > 0 that controls the amount of sparsity. Hence, when we refer to an estimate β, it should be understood as β(λ), although we suppress λ in most estimates for presentation convenience. We will discuss the tuning parameter in detail in Section 2.3.8.

Direct sparse discriminant analysis (DSDA)
The direct sparse discriminant analysis (DSDA) is proposed for binary classification under the LDA model in (2). Recall that our main interest is in estimating the coefficients β k in the Bayes rule (3). Because DSDA assumes that there are only two classes, it suffices to estimate β = Σ −1 (µ 2 − µ 1 ). In high dimensions, we assume that β is sparse. Let y i = − n 1 n if Y i = 1 and y i = n n 2 if Y i = 2. DSDA first solves the penalized least squares problem where λ > 0 is the tuning parameter, ∑ p j=1 |β j | is the LASSO penalty (Tibshirani, 1996), and β DSDA is our estimate for β. Because of the LASSO penalty, β DSDA is typically sparse. To estimate the Bayes rule, we further estimate the LDA model on the reduced data {Y i , X T i β DSDA } n i=1 . Numerical and theoretical studies show that DSDA consistently estimate the Bayes rule under mild conditions. Also, DSDA can be computed very efficiently, as (15) is a heavily-studied 1 penalized least squares problem. Our implementation utilizes glmnet to solve (15).

Regularized optimal affine discriminant (ROAD)
Regularized optimal affine discriminant (ROAD, Fan et al. (2012)) is another binary penalized discriminant analysis method for high-dimensional data. ROAD estimates β by We remark that Wu et al. (2009) independently proposed the 1 -Fisher's discriminant analysis method that closely resembles ROAD, but the developments of ROAD and the 1 -Fisher's discriminant analysis have different emphasis. ROAD clarifies several theoretical aspects of high-dimensional classification, while 1 -Fisher's discriminant analysis is developed for simultaneous testing for gene pathways. For simplicity, we focus on ROAD in what follows.
In its optimization, the constraint of 1 -norm can be recast as a 1 -penalty with parameter λ. ROAD rewrites (16) as The authors of ROAD proposed to solve (18) by replacing the nonconvex constraint with a quadratic penalty. However, we adopt a different approach to solve (18). It is showed in Mai and Zou (2013) that the solution paths of DSDA and ROAD are equivalent. In other words, for any λ > 0, there . Because DSDA produces a solution path much faster than the original proposal of ROAD, we solve ROAD by first finding the solution path of DSDA for a range of λ, and then find each correspondingλ to recover the solution path of ROAD.

Sparse optimal scoring (SOS) in binary problems
We also implement the successful discriminant analysis method, sparse optimal scoring (SOS, Clemmensen et al. (2011)). We focus on binary problems, where we are able to greatly improve the computation speed. For multiclass problems, SOS can be solved by the R package sparseLDA.
In binary problems, SOS creates a dummy variable Y dm ∈ R n×2 as a surrogate for the categorical response Y, where Y dm ik = 1{Y i = k}. Then SOS estimates coefficient by solving where X is the centered X, and θ ∈ R 2 is the score for the two classes. SOS is a popular penalized discriminant analysis method because of its impressive empirical performance. It can be solved by iteratively minimizing the objective function in (19) over θ and β.
However, we take another approach to solve SOS with lower computation cost. Mai and Zou (2013) showed that β SOS is closely related to the DSDA estimator defined in (15). Letπ y = n y n . We have that Therefore, to solve for β SOS (λ), we first find β DSDA ( λ √π 1π2 ) with DSDA, and rescale it to obtain the SOS solution. This approach avoids iteration between θ and β, and is often faster than the original algorithm for SOS.

Semiparametric sparse discriminant analysis (SeSDA)
SeSDA  fits the SeLDA model in (7) for binary problems. It is expected to have better performance than DSDA when data are heavily non-normal. SeSDA has two steps. First, we find an estimate h for the unknown function h. Second, we apply DSDA on the pseudo data ( h(X), Y).
In what follows, we focus on the estimation of h.
Two estimators have been proposed for h based on (8), the naive estimator and the pooled estimator. Without loss of generality, we assume that Class 1 has more observations than Class 2. DenoteF 1j as the empirical CDF of X j within Class 1. To avoid infinity values at tails, we further WinsorizeF 1j tô F 1j , whereF The naive estimator is shown to consistently estimate h, but in practice it is vulnerable to loss of efficiency, as it only utilizes one class of data. Therefore, the pooled estimator is proposed as a more efficient estimator.

whereĥ
(1) . The pooled estimator is usually more accurate than the naive estimator because it utilizes both classes to form an estimate for h j .

Multiclass sparse discriminant analysis (MSDA)
Up to now, we have focused on binary classifiers. In this section, we discuss a multiclass classifier under the LDA model (2). Assume that K ≥ 2. By the Bayes rule (3), we need to estimate the coefficients β k = Σ −1 (µ k − µ 1 ), k = 2, . . . , K. There is no need to estimate β 1 , as it is zero by definition. As in the binary problems, we continue to assume that the classifier is sparse in high dimensions, in the sense that only a few predictors are relevant to classification. However, this sparsity assumption has slightly different implication in multiclass problems. Note that, for any X j , if any one of β 2j , . . . , β Kj is nonzero, X j is important for classification, as it helps with distinguishing between at least one pair of classes. Therefore, in order for an X j to be unimportant, we have to have β 2j = . . . = β Kj = 0. In other words, the coefficients β 2 , . . . , β K has a group sparsity structure.
The multi-class sparse discriminant analysis (MSDA) has been proposed for fitting a sparse classifier under the context of interest. It takes note of the fact that, on the population level, we have (β 2 , · · · , β K ) = arg min Therefore, in high dimensions, MSDA replaces the parameters with the sample estimates and impose the group sparsity structure through group lasso (Yuan and Lin, 2006). More specifically, MSDA estimates β by ( β 2 , · · · , β K ) = arg min The problem in (23) can be solved by a blockwise coordinate descent algorithm  summarized in Algorithm 1. We refer to Algorithm 1 as the original MSDA algorithm. The R package msda implements such an algorithm. However, the original MSDA algorithm can be demanding on storage for high-dimensional data, because it requires the input of Σ ∈ R p×p . When p is very large, the original MSDA algorithm can be practically inapplicable. Moreover, because of the sparse nature of β, many elements in Σ are never used, and the calculation of them leads to unnecessary computation burden.
Therefore, in our implementation we modify the original MSDA algorithm for lower storage and computation cost for high-dimensional data. Note that Σ is only used in updating rule (25). We take advantage of two properties of this updating rule (25). First, given the natural element-wise property of coordinate descent algorithm, only the j-th column of covariance matrix Σ ·j is needed in each iteration. The full covariance matrix is never used during the computation process. Therefore, it is not necessary to store the huge covariance matrix. Secondly, a large number elements of β are exactly 0. Hence among the column Σ ·j , we only need to compute the rows corresponding to the nonzero coefficients. These facts motivate us to develop the modified MSDA algorithm. The modified MSDA algorithm is largely identical to the original algorithm, but with two important distinctions.
On one hand, in Step 1 we only require the input of δ k but not Σ. On the other hand, Step 3(b) in (25) is replaced withβ where T k = {i : y i = k}. By doing so, we avoid the storage and the computation of the full matrix of Σ.
In computing (26), we further use three tricks to speed up the computation. Firstly, we calculate and store all diagonal elements in the covariance matrix as they will be called multiple times. Secondly, we keep the indexes of nonzero elements in T k and update it every time we observe a new nonzero element. Hence we do not need to check all elements to locate the nonzero ones in each iteration. Thirdly, we update equation (26) by only computing elements corresponding to the nonzero indexes in T k . With these three tricks, the modified algorithm reduces the space complexity from O(p 2 ) to O(p), and is also faster than the original algorithm for large p.

Covariate-adjusted tensor classification in high dimensions (CATCH)
When X is a tensor instead of a vector, we need to fit the TDA model or the CATCH model (in presence of covariates) for better efficiency and accuracy. Pan et al. (2019) proposed the CATCH method to fit both models, but in this section we focus on the CATCH method on the TDA model, where there is no covariate. The inclusion of covariates will be discussed in Section 2.3.7.
Recall that, under the TDA model, we aim to estimate the parameters B k = µ k − µ 1 ; Σ −1 1 , . . . , Σ −1 M . We first rewrite B k as solutions to estimating equations: where for two M-way tensors A, C, A, C = ∑ j 1 ···j M a j 1 ···j M c j 1 ···j M is the inner product of two tensors. To estimate B k , we find the within-class sample mean µ k as the estimate for µ k , and moment-based unbiased estimators Σ m for Σ m ; see the formulas in Appendix C. We further add the group LASSO penalty for variable selection. Therefore, CATCH solves the following problem: CATCH can be solved by a coordinate descent algorithm with an explicit updating formula in each iteration.

Covariates adjustment
When we have additional covariates U, the CA-LDA model or the CATCH model should be fitted. Whether X is a vector or a tensor, a key step for the covariate adjustment is the estimation of α, the dependence of X on U. We use the maximum likelihood estimator (MLE). Denote U k as the sample mean of U within class k and X k as the sample mean of X within class k. Define group-wise centered data X For tensor-variate X i ∈ R p 1 ×···×p M , we let α j 1 ···j M ∈ R q be the regression coefficient of univariate X i,j 1 ···j M on multivariate U i ∈ R q . Then the MLE for α j 1 ···j M is α j 1 ···j M = ( U T U) −1 U T X j 1 ···j M , which can be expressed more explicitly as,  Afterwards, the ensemble of all α j 1 ···j M , α, is our estimator for α. The covariate-adjusted predictor is then obtained as X i − α× M+1 U i .

Selection of the tuning parameter
We recommend selecting the tuning parameter in all methods by cross validation, which is implemented in our package as supportive functions for most of the methods. In cross validation, a sequence of potential tuning parameters is supplied. For each candidate tuning parameter λ, the dataset is random split into L folds. Then we fit L classifiers, each of which is fitted on L − 1 folds of the data and validated on the remaining one fold. The average validation error rate of the L classifiers is used as a measurement of the performance of the corresponding λ. The λ with the smallest average validation error is used in our final model fitting.
If desired, our package can automatically generate a sequence of tuning parameters for all the methods. They will first compute the smallest λ that shrinks all coefficients to zero; this value is taken as the upper bound of tuning range. Then the upper bound is multiplied by a small number to generate the lower bound. Finally, a sequence of tuning parameters is uniformly generated between the lower and the upper bound.

Using the R package
The R package TULIP provides user-friendly functions to fit discriminant analysis model and perform predictions on vector and tensor data. The package can be downloaded through https://cran.rproject.org/web/packages/TULIP or install in R through install.packages( TULIP ). In installing package, the pre-required packages MASS (Venables and Ripley, 2002) for LDA model fitting, packages Matrix (Bates and Maechler, 2016) and tensr (Gerard and Hoff, 2016) for matrix and tensor operations, and the package glmnet (Friedman et al., 2010) for LASSO are also automatically installed. Users do not need to install them separately. To guarantee higher computation efficiency of the package, core algorithms of MSDA and CATCH are implemented in Fortran, which have already been compiled and can also be used directly.
Among all the six methods, there is always a tuning parameter λ to control the size of sparsity. On the implementation aspect, MSDA and CATCH also have parameter dfmax to limit the number of selected variables and will only return the solutions with number of non-zero elements less than dfmax. Furthermore, MSDA has a model option to specify version of implementation between multi.original and multi.modified. The methods are summarized in Table 2.
The functions in the package consists of two parts. One part contains core functions which generate solution paths of all the methods, including functions dsda, road, sos, SeSDA, msda and catch. Since binary classification can be regarded as a special case of multi-class problems, we also embedded DSDA into msda function. See Section 2.4.1 for details. The other part includes supportive functions to perform covariate adjustment, prediction, cross validation and handle some special cases.
To illustrate how to use the functions, we first simulate a binary vector data set named dat.vec with dimension p = 500 and sample size n k = 75. In the data set, we have X i | (Y i = k) ∼ N(µ k , Σ), where µ 1 = 0, µ 2 = Σβ, σ ij = 0.3 if i = j and σ ii = 1, β j = 0.5 for 1 ≤ j ≤ 10 and β j = 0 otherwise. We further generate a testing data set with sample size 1000 from the same distribution. Variables in dat.vec is summarized in  Moreover, we include two real data sets, GDS1615 and colorimetric sensor array data set, in the package to demonstrate usage of the functions. Data set GDS1615 (Burczynski et al., 2006) is a vector data set where observations belong to three classes. The original data set contains 127 observations and 22283 variables. Package msda preprocessed the data by computing F-test statistics of each variable , whose definition is in appendix. Hence only 127 variables are kept in the data set. Colorimetric sesor array data (CSA) was used to show the performance of discriminant analysis method (Zhong and Suslick, 2015). It records information of chemical dyes after exposed to volatile chemical toxicants to identify their classes. It contains 147 observations in 21 classes. For each observation, the predictor is a 36 × 3 matrix. We include two conditions in our dataset, but focus on the Immediately Dangerous to Life or Health (IDLH) condition.

Core functions Function dsda
The following code shows an example of utilizing DSDA. Given the data set dat.vec, we fit DSDA on {X, Y} by specifying the tuning range of parameter λ to be a sequence between [0.005, 0.3]. Hence the function will generate a solution path. Next, we apply predict function on the model and obtain the prediction for each λ and error rate. In the example, we report the minimum error rate and corresponding parameter value.
obj <-dsda(dat.vec$x, y=dat.vec$y, lambda=seq(0.005, 0.3, length.out=20)) pred <-predict(obj, dat.vec$testx) err<-apply(pred, 2, function(x){mean(x!=dat.vec$testy)}) print(min(err)) [1] 0.111 print(obj$lambda[which.min(err)]) [1] 0.02052632 If one wishes, dsda can also be used in a more automatic way. On one hand, it can be called without supplying a sequence of value for tuning parameter. The function will automatically generate a sequence based on data. On the other hand, the prediction can be performed along with model fitting if testing data is supplied. The function dsda will produce the prediction error on the testing data corresponding to each tuning parameter. See the following example.

Function SeSDA
Function SeSDA fits a semiparametric sparse discriminant analysis model on the input vector data. The simulated data dat.vec follows normal distribution within each class. We take an exponential transformation on it to violate the normality assumption. The following example shows that SeSDA achieves error rates 11%. However, if we directly apply DSDA on the data set, the minimum error rate is as high as 15.8%. Therefore, the preprocessing of SeSDA can indeed help to improve performance under this scenario. x <-exp(dat.vec$x) testx <-exp(dat.vec$testx) obj.SeSDA <-SeSDA(x, y=dat.vec$y) pred.SeSDA <-predict(obj.SeSDA, testx) err <-apply(pred.SeSDA, 2, function(x){mean(x!=dat.vec$testy)}) min(err) [1] 0.11 Further, Figure 5 shows how the distribution of the first variable changes after transformation. It is clear that both pooled and naïve transformatins result in approximately normal distribution.

Functions ROAD and SOS
Functions ROAD and SOS can generate equivalent solution paths as ROAD (Fan et al., 2012) and SOS (Clemmensen et al., 2011) methods on binary vector data, respectively. Both of the two models are fit by calling dsda function. Compared to the original package for SOS, sparseLDA, our implementation is usually faster, especially when a solution path or parameter tuning is needed. For example, to fit a solution path with 10 possible values of λs on a toy example with p = 40, our implementation reduces the computation time by half compared to sparseLDA. An example of fitting ROAD and SOS model is as follows. The lambdas passed into ROAD and SOS functions will be directly used by dsda function. The lambdas returned by the two functions are their corresponding parameters in ROAD and SOS model, respectively. Figure 6 shows the relationship between the λ's that generate the same solution.

Function msda
The function msda provides an interface to fit MSDA. Similarly to dsda, without specification of possible values of λ, the function will automatically generate a sequence of λs. Function msda can also perform predictions when testing data is supplied and make adjustments on covariates when covariates exist. We apply msda on GDS1615 data set to as a demonstration. We report the minimum training error, its corresponding parameter value and the number of non-zero variables selected by the model. data(GDS1615) x <-GDS1615$x y <-GDS1615$y set.seed(123456) teindex <-c(sample(which(y==1), sum(y==1)/3), sample (which(y==2), sum(y==2)/3), sample(which(y==3), sum(y==3)/3)) obj <-msda(x [-teindex, ]   We also note that msda has an argument model that can be specified by users to use different algorithms in MSDA. The options for model include binary, multi.original and multi.modified.

Variable
Type Dimension  The option binary can only be used in binary problems. If selected, MSDA is solved by DSDA, which gives the same solution with usually less time. However, using this option in multi-class problems will result in an error. The option multi.original indicates that MSDA is solved by the original algorithm, which requires the calculation of the full covariance matrix. When the dimension is low, multi.original is often efficient. The option multi.modified, on the other hand, solves MSDA with the modified algorithm, where only part of the covariance matrix is calculated in each iteration. This option allows MSDA to be applicable in much higher dimensions. Also, when multi.modified is selected, we suggest using relatively larger tuning parameters to account for the high dimensionality. If unspecified, model is set to be binary in binary problems. If the response variable is multi-class, the function will call multi.original implementation for p ≤ 2000 and multi.modified implementation for p > 2000.
Data set dat.vec can be simulated by code dat.ten<-sim.tensor.cov(1000) Function catch fits a CATCH model on the input tensor data. Covariates are optional for the function and the function will fit a TDA model when there is no covariate. Function catch has already integrated the adjustment step and model fitting step, hence it will automatically adjust for covariates when covariates exist. If one prefers to seperate the adjustment step, he/she can call adjten function to make adjustments and then supply the adjusted predictors into catch, which we will discuss in supportive functions. Similar to the two functions above, function catch can generate a solution path on default or user specified potential values of the parameter. It will also perform prediction when testing data is specified. To make predictions on CATCH model, user can directly apply catch function or separating adjustment and model fitting step and then call the predict function. The following example shows how to fit the model and make prediction when covariates exist. As mentioned above, functions dsda and msda shares the same arguments name for covariates. obj <-catch(dat.ten$x, dat.ten$z, dat.ten$y, dat.ten$testx, dat.ten$testz) pred <-obj$pred err <-apply (pred, 2, function(x){mean(x!=dat.ten$testy)}) min ( An example of applying CATCH to fit model and perform prediction on CSA data is as follows. catch function takes list of multi-dimensional array as input. In the dataset, x is a list of length 148, where each element is a matrix of dimension 36 × 3; y is a vector whose value ranges between 1 and 21. We use default parameter sequence of length 100 and the prediction for each value of parameter is generated. data(csa) x <-csa$IDLH y <-csa$y teindex <-seq (1,

Other functions Two special cases
We provide two more functions for two common problems in practice, binary classification and matrix classification, respectively. First, the function catch_matrix fits CATCH model on matrix data (2-way tensor), which is a special case of catch. The usage of catch_matrix is exactly the same as that of catch, with the only exception that the predictor has to be a matrix instead of higher-order tensor.
Second, our package includes the function dsda.all that integrates cross validation, model fitting and prediction. It requires the input of the training set and testing set. Then the optimal tuning parameter is chosen by cross validation on the training set, and the corresponding testing error is reported. See the following example.

Supportive functions
The package provides functions cv.dsda, cv.msda, cv.SeSDA and cv.catch to perform cross validation. For all of these functions, user can give a sequence of potential values to tune parameter. Otherwise, the function will first fit a model on the entire data set and then perform cross validation on the automatically generated λs from the entire data set. Similar as msda, user can specify which model to use in cv.msda or let the function determine by input data.
Users can also specify the number of folds by the argument nfolds. Another argument lambda.opt has two options "min" and "max". When multiple λs lead to same error rate, "min" will return the smallest tuning parameter with the lowest error rate while "max" will return the largest one. We take cv.dsda and cv.catch as two examples.
obj.dsda <-cv.dsda(dat.vec$x, dat.vec$y, nfolds = 10) obj.catch <-cv.catch(dat.ten$x, dat.ten$z, dat.ten$y, lambda.opt="min") Function adjten and adjvec implement the adjustment step for tensor and vector data, respectively. It takes training tensor/vector, covariate and response as input, and outputs the adjusted tensor/vector and adjustment coefficients α. The adjustement step has already been incorporated into the modeling  fitting functions dsda, msda and catch. When user input covariates along with tensor/vector, the model fitting functions will automatically make the adjustment. But if user do not want to use the automatic prediction in model fitting function and prefer to predict via predict, user need to first make the adjustment to obtain adjustment coefficient gamma, and pass it into function predict. Notice that making adjustment and fitting a model on the adjusted tensor and response without covariate is equivalent as fitting a model by inputting the original tensor, covariate and response labels. Examples of two approaches are given as follows.
obj <-catch(dat.ten$x, dat.ten$z, dat.ten$y, dat.ten$testx, dat.ten$testz) obj.adj <-adjten(dat.ten$x, dat.ten$z, dat.ten$y, dat.ten$testx, dat.ten$testz) obj.fit <-catch(dat.ten$x, dat.ten$z, dat.ten$y) pred <-predict(obj.fit, obj.adj$testxres, dat.ten$z, dat.ten$testz, obj.adj$gamma) There are three prediction functions corresponding to dsda, msda and catch, respectively. All of them can be directly called by predict and the function will recognize which function to use based on the input fitted model object. When covariate exists, user needs to pass the adjustment coefficient obtained from function adjten, the fitted model and testing data altogether to make predictions. Therefore, we encourage user to direct use model fitting functions msda and catch to fit model and predict categorical responses.

Real data example
In this section, we will show the performance of the models by a real data set. We considered the attention deficit hyperactivity disorder (ADHD) data set. The dataset is available on NITRC (http://fcon_1000.projects.nitrc.org/indi/adhd200) (Bellec et al., 2017). It contains three parts of information: s-MRI data which is a 3-D tensor, covariate information including age, gender and handedness which is a vector, and response label. Among all 930 individuals, there are four types of categorical labels: Typically Developing Childrem (TDC), ADHD Combined, ADHD Hyperactive and ADHD Inattentive.
We downsize the tensor to dimension 24 × 27 × 24 and consider two classification scenarios. One is to combine ADHD Hyperactive with the ADHD Combined since there are only 13 subjects in class ADHD Hyperactive. This results in a multi-class problem with three classes. The second one is to further combine TDC and ADHD Inattentive since none of these two categories have hyperactivity symptoms. This give us a binary problem. We split the dataset into a training set and testing set by ratio 8 : 2.
Given the tensor structure and existence of covariates, the most suitable approach is to apply CATCH on that. We also vectorize the tensor into vectors and stack covariates along with the long vector to apply vector methods. For binary case, DSDA is applied. For multi-class case, MSDA model with multi.modified is applied since the dimension is too large to employ multi.original. We also compared with SOS (Clemmensen et al., 2011) by its own package sparseLDA and 1 -GLM (Friedman et al., 2010) by package glmnet.
For each replicate, we perform cross validation on training data and record the classification error on testing data. The entire process was repeated for 100 times and we report the mean and standard error of the error rates. The performance is shown in Table 5.

Discussion
Package TULIP provides a toolbox to fit various sparse discriminant analysis models, including parametric models DSDA, ROAD, SOS for binary vector data, semiparametric model SeSDA for binary vector data, MSDA for multiclass vector data, and CATCH for multiclass tensor data. As a comprehensive toolbox, the package provides prediction and cross validation functions as well.
Meanwhile, the package propose an approach to handle cases when both predictor and covariates are supplied. The predictor can be vector and tensor, while the covariates are usually low-dimensional vectors. Covariates may have an effect on both response and the predictor. Therefore making adjustment and excluding the influence of covariates from predictor is important. The package includes functions to make the adjustments and can be called easily by supplying covariates in the model fitting function.

Appendices Tensor notation
On each dimension, which is named mode, a tensor is composed by vectors of length (p k × 1) called mode-k fiber, defined as A i 1 ···i k−1 I k i k+1 ···i M , I k = 1, . . . p k . Stacking the mode-1 fiber by row gives the vectorization of a tensor vec(A), which is a (∏ m p m × 1) column vector. If we unfold the tensor along the k-th mode, we obtain a matrix A (k) ∈ R p k ×∏ l =k p l .
Denote the mode-k product of a tensor A and a matrix α ∈ R d×p k by A × k α ∈ R, which results in a tensor of dimension p 1 × · · · × p k−1 × d × p k+1 × · · · × p M . Each element of the product is the product of a mode-k fiber of A and a row vector of α. In particular, the mode-k vector product of a tensor A and a vector c ∈ R p k is a (M − 1)-way tensor as a special case when d = 1. The Tucker decomposition of a tensor is defined as A = C × 1 G 1 × 2 · · · × M G M , in short of C; G 1 , . . . , G m . In particular, the vectorization of tucker decomposition has the fact that vec( C; G 1 , . . . , G M ) = (G M ⊗ · · · ⊗ G 1 ) vec(C), where ⊗ denotes Kronecker product. If X = µ + Z; Σ 1/2 1 , . . . , Σ 1/2 M , where Z ∈ R p 1 ×···×p M and all elements of Z independently follow the univariate standard normal distribution, we say X follows a tensor normal distribution X ∼ TN(µ, Σ 1 , . . . , Σ M ). The dependence structure on the j-th mode is measured by Σ j > 0. Hence, vec(X) = vec(µ) + Σ 1/2 vec(Z), where Σ = Σ M ⊗ · · · ⊗ Σ 1 .

Simulation code
Data sets dat.vec and dat.ten are used to illustrate usage of the functions. Detailed model settings are described in Section 2.4. Here are the code to simulate the two data sets.