Robust and Sparse Regression via $\gamma$-divergence

In high-dimensional data, many sparse regression methods have been proposed. However, they may not be robust against outliers. Recently, the use of density power weight has been studied for robust parameter estimation and the corresponding divergences have been discussed. One of such divergences is the $\gamma$-divergence and the robust estimator using the $\gamma$-divergence is known for having a strong robustness. In this paper, we consider the robust and sparse regression based on $\gamma$-divergence. We extend the $\gamma$-divergence to the regression problem and show that it has a strong robustness under heavy contamination even when outliers are heterogeneous. The loss function is constructed by an empirical estimate of the $\gamma$-divergence with sparse regularization and the parameter estimate is defined as the minimizer of the loss function. To obtain the robust and sparse estimate, we propose an efficient update algorithm which has a monotone decreasing property of the loss function. Particularly, we discuss a linear regression problem with $L_1$ regularization in detail. In numerical experiments and real data analyses, we see that the proposed method outperforms past robust and sparse methods.


Introduction
In high-dimensional data, sparse regression methods have been intensively studied. The Lasso (Tibshirani 1996) is a typical sparse linear regression method with L 1 regularization, but is not robust against outliers. Recently, robust and sparse linear regression methods have been proposed. The robust least angle regression (RLARS) (Khan et al. 2007) is a robust version of LARS (Efron et al. 2004), which replaces the sample correlation by a robust estimate of correlation in the update algorithm. The sparse least trimmed squares (sLTS) (Alfons et al. 2013) is a sparse version of well-known robust linear regression method LTS (Rousseeuw 1984) based on the trimmed loss function with L 1 regularization.
Recently, the robust parameter estimation using density power weight has been discussed by Basu et al. (1998), Basu et al. (2010), Jones et al. (2001), Windham (1995), Fujisawa and Eguchi (2008), Kanamori and Fujisawa (2015), and so on. The density power weight gives a small weight to the terms related to outliers and then the parameter estimation becomes robust against outliers. Among them, the γ-divergence proposed by Fujisawa and Eguchi (2008) is known for having a strong robustness, which implies that the latent bias can be sufficiently small even under heavy contamination. In addition, to obtain the robust estimate, an efficient update algorithm was proposed with a monotone decreasing property of the loss function.
In this paper, we propose the robust and sparse regression problem based on γ-divergence. First, we extend the γ-divergence to regression problem and show a strong robustness under heavy contamination even when outliers are heterogeneous. Next, we consider a loss function based on the γ-divergence with sparse regularization and propose an update algorithm to obtain the robust and sparse estimate. Fujisawa and Eguchi (2008) used a Pythagorian relation on the γ-divergence, but it is not compatible with sparse regularization. Instead of this relation, we use the Majorization-Minimization algorithm (Hunter and Lange 2004). This idea is deeply considered in a linear regression problem with L 1 regularization. Finally, in numerical experiments and real data analyses, we show that our method outperforms other robust and sparse methods.

Regression based on γ-divergence
The γ-divergence was defined for two probability density functions and its properties were investigated by Fujisawa and Eguchi (2008). In this section, the γ-divergence is extended to the regression problem, in other words, de-fined for two conditional probability density functions and the corresponding robust properties are presented.

γ-divergence for regression
We suppose g(x, y), g(y|x) and g(x) are the underlying probability density functions of (x, y), y given x and x, respectively. Let f (y|x) be another parametric conditional probability density function of y given x. Let us define the γ-cross entropy for regression by (2.1) The γ-divergence for regression is defined by ).
( 2.2) The γ-divergence for regression was first proposed by Fujisawa and Eguchi (2008) and many properties were already shown. However, we adopt the definition (2.2), which is slightly different from the past one, because (2.2) satisfies the Pythagorian relation approximately (see Sec 2.3).
Theorem 2.1. We can show that The proof is in the Appendix. In what follows, we refer to the regression based on γ-divergence as the γ-regression.

Estimation for γ-regression
Let f (y|x; θ) be the conditional probability density function of y given x with parameter θ. The target parameter can be considered by When g(y|x) = f (y|x; θ * ), we have θ * γ = θ * . Let (x 1 , y 1 ), . . . , (x n , y n ) be the observations randomly drawn from the underlying distribution g(x, y). Using the formula (2.1), the γ-cross entropy for regression, d γ (g(y|x), f (y|x; θ); g(x)), can be empirically estimated bȳ By virtue of (2.3), we define the γ-estimator bŷ In a similar way to in Fujisawa and Eguchi (2008), we can show the consistency ofθ γ to θ * γ under some conditions.

Robust properties for γ-regression
In this subsection, the robust properties are presented from two viewpoints of latent bias and Pythagorian relation. The latent bias was discussed in Fujisawa and Eguchi (2008) and Kanamori and Fujisawa (2015), which is described later. Using the results obtained there, the Pythagorian relation is shown in Theorems 2.2 and 2.3. Let f * (y|x) = f θ * (y|x) = f (y|x; θ * ) and δ(y|x) be the target conditional probability density function and the contamination conditional probability density function related to outliers, respectively. Let ǫ and ǫ(x) denote the outlier ratios which are independent and dependent of x, respectively. Under homogeneous and heterogeneous contaminations, we suppose that the underlying conditional probability density function can be expressed as and let Here we assume that which implies that ν f θ * ,γ (x) ≈ 0 for any x (a.e.) and illustrates that the contamination conditional probability density function δ(y|x) lies on the tail of the target conditional probability density function f (y|x; θ * ). For example, if δ(y|x) is the dirac function at the outlier y † (x) given x, then we have , which should be sufficiently small because y † (x) is an outlier. In this subsection, we show that θ * γ − θ * is expected to be small even if ǫ or ǫ(x) is not small. To make the discussion easier, we prepare the monotone transformation of the γ-cross entropy for regression bỹ

Homogeneous contamination
Here we provide the following proposition, which was given in Kanamori and Fujisawa (2015).
Theorem 2.2. Let ν = max {ν f θ ,γ , ν f θ * ,γ }. Then, the Pythagorean relation among g(y|x), f (y|x; θ * ), f (y|x; θ) approximately holds: The proof is in the Appendix. The Pythagorean relation implies that the minimization of the divergence from f (y|x; θ) to the underlying conditional probability density function g(y|x) is approximately the same as that to the target conditional probability density function f (y|x; θ * ). Therefore, under homogeneous contamination, we can see why our proposed method works well in terms of minimization of γ-divergence.

Heterogeneous contamination
Under heterogeneous contamination, we assume that the parametric conditional probability density function f (y|x; θ) is a location-scale family given by where s(y) is a probability density function, σ is a scale parameter and q(x; ξ) is a location function with a regression parameter ξ, e.g., q(x; ξ) = ξ T x. Then, we can obtain That does not depend on the explanatory variable x. Here, we provide the following proposition, which was given in Kanamori and Fujisawa (2015).
The second term can be approximated to be 0 from the condition ν f θ ,γ ≈ 0 and ǫ(x) < 1 as follows: We see from Proposition 2.2 and (2.4) that Therefore, under heterogeneous contamination in a location-scale family, it can be expected that the latent bias θ * γ − θ * is small even if ǫ(x) is not small. Moreover, we can show the following theorem, using Proposition 2.2.
Then, the following relation among g(y|x), f (y|x; θ * ), f (y|x; θ) approximately holds: The proof is in the Appendix. The above is slightly different from a conventional Pythagorean relation, because the base measure changes from g(x) to (1 − ǫ(x))g(x) in part. But, it also implies that the minimization of the divergence from f (y|x; θ) to the underlying conditional probability density function g(y|x) is approximately the same as that to the target conditional probability density function f (y|x; θ * ). Therefore, under heterogeneous contamination in a location-scale family, we can see why our proposed method works well in terms of minimization of γ-divergence.

Estimation for sparse γ-regression
The empirical estimation of the γ-cross entropy with a penalty term can be given by where P (θ) is a penalty for parameter θ and λ is a tuning parameter for the penalty term. As an example of penalty term, we consider L 1 (Lasso, Tibshirani 1996), elasticnet (Zou and Hastie 2005), and so on. The sparse γ-estimator can be proposed bŷ (3.1) To obtain the minimizer, we propose the iterative algorithm by Majorization-Minimization algorithm (MM algorithm) (Hunter and Lange 2004).

MM algorithm for sparse γ-regression
The MM algorithm is constructed as follows. Let h(η) be the objective function. Let us prepare the majorization function h M M satisfying where η (m) is the parameter of the m-th iterative step for m = 0, 1, 2, . . .. Let us consider the iterative algorithm by Then, we can show that the objective function h(η) monotonically decreases at each step, because We construct the majorization function for the sparse γ-regression by the following inequality: are positive. The inequality (3.2) holds from Jensen's inequality. Here, we take . We can propose the majorization function as follows: f (y l |x l ; θ (m) ) γ and const is a term which does not depend on the parameter θ.
The first term on the original target function h(θ) is a mixture type of densities, which is not easy to be optimized, while the first term on h M M (θ|θ (m) ) is a weighted log-likelihood, which is often easy to be optimized.

Sparse γ-linear regression
Let f (y|x; θ) be the conditional density with θ = (β 0 , β, σ 2 ), given by where φ(y; µ, σ 2 ) is the normal density with mean parameter µ and variance parameter σ 2 . Suppose that P (θ) is the L 1 regularization ||β|| 1 . Then, after a simple calculation, we have This function is easy to be optimized by an update algorithm. For a fixed value of σ 2 , the function h M M is almost the same as Lasso except for the weight, so that it can be updated using the coordinate decent algorithm with a decreasing property of the loss function. For a fixed value of (β 0 , β T ) T , the function h M M is easy to be minimized. Consequently, the update formula can be obtained by ( 3.7) It should be noted that h M M is convex with respect to parameter β 0 , β and has the global minimum with respect to parameter σ 2 , but not convex with respect to some of them, so that the starting value of update formula is important. This issue is often important in robust statistics and discussed in Sec 4.4. Moreover, we explain how to implement coordinate decent algorithm. In practice, we also use active set strategy (Friedman et al. 2007). Specifically, for a given β (m) , we only update the non-zero coordinates of β (m) in the active set, until they are converged. Then, the non-active set parameter estimates are updated once. When they remain zero, the coordinate descent algorithm stops.

Robust cross-validation
In sparse regression, a regularization parameter is often selected via some criterion. Cross-validation is often used for selecting the regularization parameter. However, the ordinal cross-validation will fail due to outliers. Therefore, we propose the robust cross-validation based on the γ-cross entropy. Let θ γ be the robust estimate based on the γ-cross entropy. The ordinal cross validation is based on the squared error and it can also be constructed using the KL-cross entropy with normal density. The cross-validation based on the γ-cross entropy can be given by is the γ-estimator deleting the i-th observation and γ 0 is an appropriate tuning parameter. We can also adopt the K-fold cross validation to reduce the computational task (Hastie et al. 2010).
Here, we give a small modification on the above. We often focus only on the mean structure for prediction, not on the variance parameter. Therefore, in this paper,θ

Redescending Property
First, we review a redescending property on M-estimation (see, e.g., Maronna et al. 2006), which is often used in robust statistics. Suppose that the estimating equation is given by n i=1 ζ(z i ; θ) = 0. Letθ be a solution of the estimating equation. The bias caused by outlier z o is expressed asθ n=∞ − θ * , wherê θ n=∞ is the limiting value ofθ and θ * is the true parameter. We hope the bias is small even if the outlier z o exists. Under some conditions, the bias can be approximated to ǫIF(z o ; θ * ), where ǫ is a small outlier ratio and IF(z; θ * ) is the influence function. The bias is expected to be small when the influence function is small. The influence function can be expressed as IF(z; θ * ) = Aζ(z; θ * ), where A is a matrix independent of z, so that the bias is also expected to be small when ζ(z o ; θ * ) is small. In particular, the estimating equation is said to have a redescending property if ζ(z; θ * ) goes to zero as ||z|| goes to infinity. This property is favorable in robust statistics, because the bias is expected to be sufficiently small when z o is very large.
Recall that the estimate of the sparse γ-linear regression is the minimizer of the loss function Then, the estimating equation is given by This can be expressed by the M-estimation formula given by We can easily show that as ||y|| goes to infinity, φ(y; β 0 + x T β, σ 2 ) goes to zero and φ(y; β 0 + x T β, σ 2 )s(y|x; θ) also goes to zero. Therefore, the function ψ(y|x; θ) goes to zero as ||y|| goes to infinity, so that the estimating equation has a redescending property.

Numerical experiment
In this section, we compare our method (Sparse γ-linear regression) with the representative sparse linear regression method, Least absolute shrinkage and selection operator (Lasso) (Tibshirani 1996), and the robust and sparse regression methods, sparse least trimmed squares (sLTS) (Alfons et al. 2013) and robust least angle regression (RLARS) (Khan et al. 2007).

Regression models for simulation
We used the simulation model given by The sample size and the number of explanatory variables were set to be n = 100 and p = 100, 200, respectively. The true coefficients were given by β 1 = 1, β 2 = 2, β 4 = 4, β 7 = 7, β 11 = 11, β j = 0 for j ∈ {0, . . . , p}\{1, 2, 4, 7, 11}. (4.1) We arranged a broad range of regression coefficients to observe sparsity for various degrees of regression coefficients. The explanatory variables were generated from a normal distribution N(0, Σ), where the covariance matrix Σ = (σ ij ) 1≤i,j≤p was given by σ ij = ρ |i−j| . We generated 100 random samples. Outliers were incorporated into simulations. We investigated two cases for outlier ratio (ǫ = 0.1 and 0.3) and two cases for outlier pattern: (a) The outliers were generated around the middle part of the explanatory variable, where the explanatory variables were generated from N(0, 0.5 2 ) and the error terms were generated from N(20, 0.5 2 ). (b) The outliers were generated around the edge part of the explanatory variable, where the explanatory variables were generated from N(−1.5, 0.5 2 ) and the error terms were generated from N(20, 0.5 2 ).

Performance measure
The root mean squared prediction error (RMSPE) and mean squared error (MSE) were examined to verify predictive performance and fitness of regres-sion coefficient: where y * i and x * i (i = 1, . . . , n) is the test sample generated from the simulation model without outliers and β * j 's are the true coefficients. The true positive rate (TPR) and true negative rate (TNR) were also reported to verify the sparsity:

Comparative methods
In this subsection, we explain three comparative methods; Lasso, RLARS, sLTS. The Lasso is performed by R-package "glmnet". The regularization parameter λ lasso is selected by grid search via cross-validation in "glmnet". We used "glmnet" by default.
The RLARS is performed by R-package "robustHD". This is a robust version of LARS (Efron et al. 2004). The optimal model is selected via BIC by default.
The sLTS is performed by R-package "robustHD". The sLTS has the regularization parameter λ sLT S and the fraction parameter α of squared residuals used for trimmed squares. First, the regularization parameter λ sLT S is selected by grid search via BIC. The number of grids is 40 by default. However, we considered that this would be small under heavy contamination. Therefore, we used 80 grids under heavy contamination. Next, the fraction parameter α is 0.75 by default. In the case of α = 0.75, the ratio of outlier is less than 25%. We considered this would be small under heavy contamination and large under low contamination in terms of statistical efficiency. Therefore, we used 0.65, 0.75, 0.85 as α under low contamination and 0.50, 0.65, 0.75 under heavy contamination.

Details of our method 4.4.1 Initial points
In our method, we need an initial point to obtain the estimate, because we use the iterative algorithm proposed in Sec 3.3. The estimate of other conventional robust and sparse regression methods would give a good initial point. For another choice, the estimate of the RANSAC (random sample consensus) algorithm would also give a good initial point. In this experiment, we used the estimate of the sLTS as an initial point.

How to choose tuning parameters
In our method, we have to choose some tuning parameters. The parameter γ in γ-divergence was set to 0.1 or 0.5. The parameter γ 0 in robust crossvalidation was set to 0.5. In our experience, the result via Rocv is not sensitive to the selection of γ 0 when γ 0 is enough large, e.g. γ 0 = 0.5, 1. The parameter λ of L 1 regularization is often selected via grid search. We used 50 grids on range [0.05λ 0 , λ 0 ] with log scale, where λ 0 is an estimate of λ which would shrink regression coefficients to zero. Table 1 is the low contamination case with outlier pattern (a). For the RM-SPE, our method outperformed other comparative methods (the oracle value of the RMSPE is 0.5). For the TPR and TNR, the sLTS showed a similar performance to our method. The Lasso presented the worst performance, because it is sensitive to outliers. Table 2 is the heavy contamination case with outlier pattern (a). For the RMSPE, our method outperformed other comparative methods except in the case (p, ǫ, ρ) = (100, 0.3, 0.2) for sLTS with α = 0.5. The Lasso also presented a worse performance and furthermore the sLTS with α = 0.75 showed the worst performance due to a lack of truncation. For the TPR and TNR, our method showed the best performance.    Table 3 is the low contamination case with outlier pattern (b). For the RMSPE, our method outperformed other comparative methods (the oracle value of the RMSPE is 0.5). For the TPR and TNR, the sLTS showed a similar performance to our method. The Lasso presented the worst performance, because it is sensitive to outliers. Table 4 is the heavy contamination case with outlier pattern (b). For the RMSPE, our method outperformed other comparative methods. The sLTS with α = 0.5 showed the worst performance. For the TPR and TNR, it seems that our method showed the best performance.

Real data analyses
In this section, we use two real datasets to compare our method with comparative methods in real data analysis. We show the best result of comparative methods among some parameter situations (e.g., Sec 4.3), unless otherwise noted.

NCI-60 cancer cell panel
We applied our method and comparative methods to regress protein expression on gene expression data at the cancer cell panel of the National Cancer Institute. Experimental conditions were set in the same way as in Alfons et al. (2013) as follows. The gene expression data were obtained with an Affymetrix HG-U133A chip and normalized GCRMA method, resulting in a set of p = 22, 283 explanatory variables. The protein expressions based on 162 antibodies were acquired via reverse-phase protein lysate arrays and log 2 transformed. One observation had to be removed since all values were missing in the gene expression data, reducing the number of observations to n = 59. Then, the KRT18 antibody was selected as the response variable because it had the largest MAD among 162 antibodies, i.e., KRT18 may include a large number of outliers. Both the protein expressions and the gene expression data can be downloaded via the web application CellMiner (http://discover.nci.nih.gov/cellminer/). As a measure of prediction performance, the root trimmed mean squared prediction error (RTMSPE) was computed via leave-one-out cross-validation given by where (e) 2 [1:n] ≤ · · · ≤ (e) 2 [n:n] are the order statistics of e 2 and h = ⌊(n + 1)0.75⌋. The choice of h is important because it is preferable for estimating prediction performance that trimmed squares does not include outliers. We set h in the same way as in Alfons et al. (2013), because the sLTS detected 13 outliers in Alfons et al. (2013). In this experiment, we used the estimate of the RANSAC algorithm as an initial point instead of the sLTS because the sLTS required high computational cost in such a high dimensional data.  Table 5 shows that our method outperformed other comparative methods for the RTMSPE in high dimensional data. Our method presented the smallest RTMSPE with the second smallest number of explanatory variables. The RLARS presented the smallest number of explanatory variables, but a much larger RTMSPE than our method.

Protein homology dataset
We applied our method and comparative methods to protein sequences dataset used for KDD-Cup 2004. Experimental conditions were set in the same way as in Khan et al. (2007) as follows. The whole dataset consists of n = 145, 751 protein sequences which has 153 blocks corresponding to native protein. Each data point in particular block is a candidate homologous protein. There were 75 variables in the data set: the block number (categorical) and 74 measurements of protein features. The first protein feature was used as the response variable. Then, 5 blocks with a total of n = 4, 141 protein sequences were selected because they contained the highest proportions of homologous proteins (and hence the highest proportions of potential outliers). The data of each block was split into two almost equal parts to get a training sample of size n tra = 2, 072 and a test sample of size n test = 2, 069. The number of explanatory variables was p = 77, consisting of 4 block indicators (variables 1-4) and 73 features. The whole protein, train and test dataset can be downloaded from http://users.ugent.be/~svaelst/software/RLARS.html. As a measure of prediction performance, the root trimmed mean squared prediction error (RTMSPE) was computed for the test sample given by where (e) 2 [1:ntest] ≤ · · · ≤ (e) 2 [ntest:ntest] are the order statistics of e 2 and h = ⌊(n test + 1)0.99⌋ , ⌊(n test + 1)0.95⌋ or ⌊(n test + 1)0.9⌋. In this experiment, we used the estimate of the sLTS as an initial point.  Table 6 shows that our method outperformed other comparative methods for the RTMSPE. Our method presented the smallest RTMSPE with the largest number of explanatory variables. It might seem that other methods gave the smaller number of explanatory variables than necessary.