Adaptive Capped Least Squares
This package includes two methods applied to minimize the adaptive capped least squares loss: randomized gradient descent method and gradient descent method with initials obtained from CPLEX.
Suppose we observe data vectors (xi,yi) that follow a linear model yi=xiTβ+εi, i=1,...n, where yi is a univariate response, xi is a d-dimensional predictor, β denotes the vector of regression coefficients, and εi is a random error. We propose the adpative capped least squares loss, ℓ(x)=x2/2 if |x| ≤ τ; τ2/2, if |x| > τ, where τ=τ(n) > 0 is referred to as the adaptive capped least squares parameter. The proposed methods are applied to find β that minimizes L(β)= n-1 ∑ ℓ(yi-xiT β ).
Install ACLS from GitHub:
install.packages("devtools")
library(devtools)
devtools::install_github("rruimao/ACLS/ACLS")
library(ACLS)
Install Rcplex if the user wants to use CPLEX to obtain the intial or the solution to the problem:
See the INSTALL file from https://cran.r-project.org/web/packages/Rcplex/index.html.
There are four functions in this package:
- GD: Gradient descent method with user-defined initial.
- RGD: Randomized gradient descent method.
- picksamples: Pick subsamples to obtain the initial.
- cplexcoef: Get the coefficients for the problem solved by CPLEX.
We present two examples: random generated data with y-outliers and random generated data with x-outliers and y-outliers.
we generate contaminated random errors εi from a mixture of normal distribution 0.9N(0,1)+0.1N(10,1) and xi's are independently and identically distributed (i.i.d.) from N(0,Σ) where Σ=0.5|j-k|. We set β* =(0,3,4,1,2,0)T to generate yi. We provide one example of this type, "ex_1.Rdata", and it can be downloaded from example file.
load("ex_1.Rdata")
Y<-ex_1$Y
X<-ex_1[,!(names(ex_1) %in% "Y")]
We first apply our proposed methods, randomized gradient descent (denoted by ACLS) and gradient descent with intials obtained from CPLEX (denoted by ACLS-h), to the data.
We randomly generate 10 initials β* ~ Unif(B2(τ)), where Unif(B2(τ)) is a uniform distribution on the l2-ball B2(τ)={x: ||x||2 ≤ τ }. This method finds the initial that provides the smallest adaptive capped least squares loss.
tau<-sqrt(n)/log(log(n))
iter<-10
beta.rgd<-RGD(X,Y,tau,iter)
We first sample 30% data and use CPLEX to find the optimal solution as the initial, then we apply gradient descent method with this initial.
library("Rcplex")
n_ratio<-0.3
n_sp<-round(n*n_ratio)
Sample<-picksamples(X,Y,n_ratio)
Coef_sp<-cplexcoef(Sample$X_sp,Sample$Y_sp,tau)
ans_sp<-Rcplex(Coef_sp$f,Coef_sp$A,Coef_sp$b,Coef_sp$Q,Coef_sp$lb,Coef_sp$ub,vtype=Coef_sp$vtype)
beta_0<-ans_sp$xopt[(3*n_sp+1):(3*n_sp+d+1)]
beta.gdc<-GD(beta_0,tau,X,Y)
We then compare mean square errors (MSEs), defined as || βˆ-β*||2, of these two methods with MSEs of ordinary least squares method (OLS), Huber method with adaptive resistant parameter (denoted by AHR) and least trimmed squares method (LTS). Results obtained from the CPLEX on the whole dataset (denoted by ACLS-C) are treated as benchmark.
install.packages("robustbase")
library("robustbase")
install.packages("robustreg")
library("robustreg")
Z.OLS<-lm(Y~X1+X2+X3+X4+X5,data=ex_1)
beta.OLS<-Z.OLS$coefficients
Z.Huber<-robustRegH(Y~X1+X2+X3+X4+X5,data=ex_1,tune=tau)
beta.Huber<-Z.Huber$coefficients
Z.LTS<-ltsReg(X[,-1],Y,intercept=TRUE,adjust=TRUE)
beta.LTS<-Z.LTS$coefficients
Coef<-cplexcoef(X,Y,tau)
ans<-Rcplex(Coef$f,Coef$A,Coef$b,Coef$Q,Coef$lb,Coef$ub,vtype=Coef$vtype)
beta.cplex<-ans$xopt[(3*n+1):(3*n+d+1)]
We summarize the MSEs of all methods in the following table.
OLS | AHR | LTS | ACLS | ACLS-h | ACLS-C | |
---|---|---|---|---|---|---|
MSE | 4.3016 | 4.3016 | 0.5026 | 0.0626 | 0.0625 | 0.0636 |
we generate contaminated random errors εi from a mixture of normal distribution 0.9N(0,1)+0.1N(10,1) and xi's are independently and identically distributed (i.i.d.) from N(0,Σ) where Σ=0.5|j-k|. We then add a random perturbation vector zi ~ N(10 × 1d-1,Id-1) to each covariate xi in the contaminated samples. We also use β* =(0,3,4,1,2,0)T and use uncontaminated xi to generate yi. We provide one example of this type, "ex_2.Rdata", and it can be downloaded from example file.
load("ex_2.Rdata")
Y<-ex_2$Y
X<-ex_2[,!(names(ex_1) %in% "Y")]
We use the same code in the first example with the new data to get estimators for all methods. We also collect the MSEs in the following table.
OLS | AHR | LTS | ACLS | ACLS-h | ACLS-C | |
---|---|---|---|---|---|---|
MSE | 18.6446 | 18.6446 | 0.7033 | 0.3558 | 0.3555 | 0.3569 |
Sun, Q., Mao, R. and Zhou, W.-X. Adaptive capped least squares. Paper
Huber, P. J. (1973). Robust regression: asymptotics, conjectures and Monte Carlo. The Annals of Statistics 1 799-821. Paper
Rousseeuw, P. J. and Driessen, K. V. (1999). A fast algorithm for the minimum covariance determinant estimator. Technometrics 41 212-223. Paper
Rousseeuw, P. J. (1984). Least median of squares regression. Journal of the American Statistical Association 79 871-880. Paper
Sun, Q., Zhou, W.-X. and Fan, J. (2020). Adaptive Huber regression. Journal of the American Statistical Association 115 254-265. Paper