0
$\begingroup$

I was wondering what would be the most efficient algorithm to solve quantile regression with a single predictor and no intercept? I tried a Brent line search, but unfortunately that's no faster than the interior point method implemented in quantreg (see below). So I was wondering if there would be potentially faster ways to do a quantreg model with a single predictor? Ideally I would also like to be able to use warm starts, which would exclude interior point methods...

EDIT: In this paper, "cqrReg: An R Package for Quantile and Composite Quantile Regression and Variable Selection", on page 14 I see a coordinate descent algorithm mentioned to solve quantile regression (implemented in R in https://github.com/cran/cqrReg/blob/master/R/QR.cd.R and https://github.com/cran/cqrReg/blob/master/src/QRCD.cpp; qrfit(x, y, tau, beta0, method="cd")). This might fit the bill, except that the current implementation always fits an intercept, which I don't want. Would anyone know if that algorithm could be adapted to also allow models to be fit without intercept, and what changes in the Rcpp code this would entail?

set.seed(1)
n=500
p=1
a=rnorm(n*p, mean = 10, sd = 1)
x=matrix(a,n,p)
beta=rnorm(p,1,1)
beta=matrix(beta,p,1)
y=x%*%beta-matrix(rnorm(n,0.1,1),n,1)
tau=0.01 # example tau, here to get underapproximation

library(microbenchmark)

# quantreg solution
library(quantreg)
microbenchmark(microbenchmark(rq.fit.br(x, y, tau = tau)$coefficients)) # quantreg, 149 µs, beta=0.84

# using Brent line search, using [0, least square estimate] search interval
# function giving log(quantile loss)
logquantloss = function(beta, x, y, tau) { r=y-beta*x   
                                           log(sum(r*(r>0) - (1-tau)*r)) } 
microbenchmark(optimize(logquantloss, 
                       interval=c(0, mean(x*y)/mean(x*x)), # least square beta = maximum
                       x=x, y=y, tau=tau)$minimum) # 184 µs, beta=0.84
$\endgroup$
  • 1
    $\begingroup$ When all regressors are positive, clearly the solution is given by a suitable quantile of the $y_i/x_i,$ which is found with $O(n)$ effort. The generalization to data that straddle zero shouldn't be much harder. $\endgroup$ – whuber Feb 21 at 21:53
  • $\begingroup$ Ha yes of course - stupid that I hadn't thought of that... microbenchmark(quantile(y/x,tau)) runs in 139 µs and indeed gives 0.844 in my example... Next question then would be if you would happen to know of the fastest Rcpp quantile() implementation in R that runs in O(n) time? $\endgroup$ – Tom Wenseleers Feb 21 at 22:55
  • $\begingroup$ I don't--but it's easy to code yourself. Emulate a Quicksort algorithm or use a heap. $\endgroup$ – whuber Feb 22 at 12:24
  • $\begingroup$ Thanks - this Rcpp quantile function here, stackoverflow.com/questions/26786078/…, is already 4x faster than quantile() but will keep on looking. Quicksort is only O(n log n) efficient though isn't it, not O(n)? $\endgroup$ – Tom Wenseleers Feb 22 at 14:35
  • 1
    $\begingroup$ Thanks for the pointers - found what I need here: yiqun-dai.blogspot.com/2017/03/… and that runs in 6 µs as opposed to 139 µs for my example :-) Great! $\endgroup$ – Tom Wenseleers Feb 22 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.