A gradient search maximization algorithm for the asymmetric Laplace likelihood

Abstract

The asymmetric Laplace likelihood naturally arises in the estimation of conditional quantiles of a response variable given covariates. The estimation of its parameters entails unconstrained maximization of a concave and non-differentiable function over the real space. In this note, we describe a maximization algorithm based on the gradient of the log-likelihood that generates a finite sequence of parameter values along which the likelihood increases. The algorithm can be applied to the estimation of mixed-effects quantile regression, Laplace regression with censored data, and other models based on Laplace likelihood. In a simulation study and in a number of real-data applications, the proposed algorithm has shown notable computational speed.

Keywords:

• quantile regression
• Laplace regression
• mixed-effects quantile regression
• asymmetric Laplace distribution
• direct search optimization algorithms

1. Introduction

An increasing number of recently proposed methods for the estimation of conditional quantiles of a continuous outcome variable given covariates make use of the asymmetric Laplace distribution.[1–6] We describe a maximization algorithm for the asymmetric Laplace likelihood and compare its large-sample performance with that of the Frisch–Newton (FN) interior-point method for the estimation of quantile regression. Simple extensions of the presented algorithm have been applied to the estimation of Laplace regression with censored data [7] and mixed-effects quantile regression as an alternative to the Gibbs sampler.[8,9]

Let y_i, i=1, … , n, be a sample of continuous variables and x_i be k-dimensional vectors of corresponding covariates. We consider the regression model:

$y_i=x_i^{\prime }{\beta }_p+u_i,$

where the residual u_i follows an unspecified distribution with $P(u_i\le 0)=p$ for a given probability p∈(0, 1). If model (1) is correct, then x_i′β_p is the p-quantile of the conditional distribution of y_i given x_i. The unknown regression coefficient ${\beta }_p\in {\mathbb{R}}^v$ depends on p and can be estimated by minimizing a weighted sum of absolute residuals with computationally efficient methods.[10–12]

It is known [13] that this minimization problem can be seen as the maximization of a likelihood function by assuming that the regression residual u_i follows the asymmetric Laplace density function $f(u_i)=p(1-p)\exp \{u_i(I_{u_i\le 0}-p)\}$. Here and throughout, I_A denotes the indicator function of the set A.

The log-likelihood function of model (1)

$l({\beta }_p)=n\log \{p(1-p)\}+\sum _{i=1}^n(y_i-x_i^{\prime }{\beta }_p)(I_{y_i\le x_i^{\prime }{\beta }_p}-p)$

is continuous at all ${\beta }_p\in {\mathbb{R}}^v$. Its first derivative, $\mathrm{\partial }l({\beta }_p)/\mathrm{\partial }{\beta }_p$, is continuous everywhere except at the points β_p such that $y_i=x_i^{\prime }{\beta }_p$ for at least one observation $i\in \{1,2,\dots ,n\}$. We define the gradient

$s({\beta }_p)=-\sum _{i=1}^n{x}_i(I_{y_i\le x_i^{\prime }{\beta }_p}-p).$

The function, s(β_p), ${\mathbb{R}}^v\mapsto {\mathbb{R}}^v$, equals the first derivative $\mathrm{\partial }l({\beta }_p)/\mathrm{\partial }{\beta }_p$ at all points of the parameter space where the latter is defined.

2. Gradient-search algorithm

In this section, we describe an algorithm for the unconstrained maximization of the log-likelihood function, l(β_p), with respect to the parameter, ${\beta }_p\in {\mathbb{R}}^v$. Briefly, from a current parameter value, the algorithm searches the positive semi-line in the direction of the gradient s(β_p) for a new parameter value at which the likelihood is larger. The algorithm stops when the change in the log-likelihood is less than a specified tolerance. Convergence is guaranteed by the continuity and concavity of the log-likelihood as shown in Section 3.

The algorithm can be summarized by the following steps:

1. (1) Set k=0 and the initial values ${\beta }_p^0\in {\mathbb{R}}^v$ and δ⁰>0

2. (2) If $l({\beta }_p^k)>l\{{\beta }_p^k+{\delta }^{k}s({\beta }_p^k)\}$

   • (a) then set ${\delta }^{k+1}=a{\delta }^k$

   • (b) else if $\{l({\beta }_p^{k+1})-l({\beta }_p^k)\}>ϵ$

     • (i) then set ${\beta }_p^{k+1}={\beta }_p^k+{\delta }^{k}s({\beta }_p^k)$; ${\delta }^{k+1}=b{\delta }^k$

     • (ii) else return ${\beta }_p^{k+1}$; stop

3. (3) Set k=k+1; go to step 2.

The algorithm requires setting the following initial values: the starting value of the parameter, ${\beta }_p^0\in {\mathbb{R}}^v$; the initial step-length, δ⁰>0; the factor for shortening the step-length, a∈(0, 1); the factor for expanding the step-length, b≥1; and the tolerance for the change in the log-likelihood, ε>0. In addition, the constant p∈(0, 1) and the data y_i and x_i are required to evaluate the functions l(β_p) and s(β_p).

3. Convergence of the algorithm

This section discusses the convergence of the gradient search (GS) algorithm defined in Section 2 and states some properties of the log-likelihood function l(β_p) defined in Equation (2).

Definition 1

Let $\mathrm{\Omega }=\{{\beta }_p^{\ast }\in {\mathbb{R}}^v:l({\beta }_p^{\ast })\ge l({\beta }_p),\mathrm{\forall }{\beta }_p\in {\mathbb{R}}^v\}$ be the set of maximizers of the log-likelihood function.

The log-likelihood function l(β_p) is a sum of continuous, concave functions and therefore itself continuous and concave. Its continuity and concavity guarantee that a unique maximum exists and that the set of maximizers, $\mathrm{\Omega }\subset {\mathbb{R}}^v$, is compact. The well-known results stated in Lemma 1 is instrumental in establishing convergence of the GS algorithm in Theorem 1.

lemma 1

For any sequence of parameters $\{{\beta }_p^k\}$ such that $l({\beta }_p^k)<l({\beta }_p^{k+1})\mathrm{\forall }k\in \mathbb{N}$ $\underset{k\to \mathrm{\infty }}{lim}{\beta }_p^k={\beta }_p^{\ast }\in \mathrm{\Omega }.$

Proof The statement in Lemma 1 follows by the continuity and concavity of the log-likelihood function l(β_p) defined in Equation (2).

The following Theorem 1 shows that the sequence of parameter values ${\beta }_p^{k+1}={\beta }_p^k+{\delta }^{k}s({\beta }_p^k)$ iteratively produced by Step 2.b.i of the GS algorithm satisfies the condition in Lemma 1.

Theorem 1

For any given point ${\beta }_p\in {\mathbb{R}}^v,$ such that $y_i\ne x_i^{\prime }{\beta }_p$ for all $i\in \{1,\dots ,n\},$ there exists a value δ₀>0 such that $l({\beta }_p)\le l({\beta }_p+\delta s({\beta }_p))$ for every $0<\delta <{\delta }_0$.

Proof The directional derivative of l(β_p) in the direction v≠0 evaluated at β_p is $\begin{array}{rl}{\left.\frac{\mathrm{\partial }}{\mathrm{\partial }\delta }l({\beta }_p+\delta v)\right|}_{\delta =0}& ={\left.\frac{\mathrm{\partial }}{\mathrm{\partial }\delta }\sum _{i=1}^n(y_i-x_i^{\prime }{\beta }_p-\delta x_i^{\prime }v)(I_{y_i\le x_i^{\prime }{\beta }_p+\delta x_i^{\prime }v}-p)\right|}_{\delta =0}\\ & =-\sum _{i=1}^n{x}_i^{\prime }v(I_{y_i\le x_i^{\prime }{\beta }_p}-p{)}^{I_{y_i\ne x_i^{\prime }{\beta }_p}}(I_{x_i^{\prime }v\ge 0}-p{)}^{I_{y_i=x_i^{\prime }{\beta }_p}}.\end{array}$ Assuming that $y_i\ne x_i^{\prime }{\beta }_p$ for all $i\in \{1,\dots ,n\}$, the directional derivative in the direction of the gradient $v=s({\beta }_p)$ evaluated at β_p ${\left[-\sum _{i=1}^n{x}_i(I_{y_i\le x_i^{\prime }{\beta }_p}-p)\right]}^{\prime }\left[-\sum _{i=1}^n{x}_i(I_{y_i\le x_i^{\prime }{\beta }_p}-p)\right]=\parallel s({\beta }_p){\parallel }^2$ is non-negative and greater than that in any other direction $v\in {\mathbb{R}}^v$ such that $\parallel v\parallel =\parallel s({\beta }_p)\parallel$.

If $\parallel s({\beta }_p^k){\parallel }^2>0$, then the directional derivative at ${\beta }_p^k$ in the direction of the gradient s(β_p) is positive and there exists a δ₀>0 such that the log-likelihood function increases in a δ₀-neighbourhood of ${\beta }_p^k$ in the direction of s(β_p); that is, $\mathrm{\exists }{\delta }_0>0:l({\beta }_p^k)<l({\beta }_p^k+{\delta }^{k}s({\beta }_p^k))=l({\beta }_p^{k+1})\mathrm{\forall }\delta \in (0,{\delta }_0)$. Step 2a of the GS algorithm shortens the step-length δ^k by a factor a∈(0, 1) until ${\delta }^k<{\delta }_0$. If $\parallel s({\beta }_p^k){\parallel }^2=0$, then $\parallel s({\beta }_p^k)\parallel =0$. The first-order condition $\parallel s({\beta }_p^k)\parallel =0$, along with the continuity and concavity of l(β_p) over the entire space ℝ^v, implies that ${\beta }_p^k$ is a maximizer of the log-likelihood, i.e. ${\beta }_p^k\in \mathrm{\Omega }$.

Theorem 1 assumes that $y_i\ne x_i^{\prime }{\beta }_p^k$ for all observations $i\in \{1,\dots ,n\}$. With the GS algorithm defined in Section 2, the event that y_i is numerically equal to $x_i^{\prime }{\beta }_p^k$ for an observation i occurs with probability zero. If it did occur, the GS would move off of it in the direction s(β_p) defined in Equation (3).

While meeting the condition of Lemma 1 guarantees that a parameter sequence $\{{\beta }_p^k\}$ convergences to a maximizer, it does not guarantee that it converges efficiently. Indeed, in an infinite number of iterations, any algorithm that generates a sequence of increasing likelihood values would converge to the maximum. For example, a trivial compass search along the orthogonal directions of ℝ^v would be a convergent, albeit very inefficient, alternative to the GS. Efficiency of the GS algorithm ensues from its moving from the current parameter value in the direction of maximum local increase.

4. A simulation study

We performed a simulation study and compared the proposed GS method with the FN interior point method in large samples. In this case, the latter method has been shown to outperform the Barrodale–Roberts method.[11]

With the statistical program R, we generated 7200 samples, which consisted of 100 samples under each of 72 scenarios obtained from the combination of 3 quantiles, 2 sample sizes, 4 regression models, and 3 distributions for the regression residual. The three quantiles were $p\in \{0.50,0.75,0.90\}$. The two sample sizes were $n\in \{{10}^4,{10}^5\}$. The four regression models were (1) y_i=i, with $i\in \{1,\dots ,n\}$, (2) $y_i=1+x_i^{(1)}+x_i^{(2)}+u_i$, (3) $y_i=1+x_i^{(3)}+\cdots +x_i^{(10)}+u_i$, and (4) $y_i=-30+5000x_i^{(11)}+2x_i^{(12)}+0.5x_i^{(13)}+(0.5+x_i^{(11)}-0.5x_i^{(12)}+0.5x_i^{(13)})u_i$; where $x_i^{(j)}\sim$ Normal(0,1), $j\in \{1,\dots ,10\}$, $x_i^{(11)}\sim$ Uniform(0,1), $x_i^{(12)}\sim$ Bernoulli(0.5), and $x_i^{(13)}\sim$ T(3). The three distributions of the regression residual were u_i ∼ Normal(0,1), u_i ∼ T(3), and u_i ∼ Lognormal(0,1).

The GS algorithm was implemented as an R function with the following starting values: ${\beta }_p^0$ equal to the least-squares estimate, δ⁰ equal to the sample standard deviation of y, a=0.5, b=1.25, and $ϵ={10}^{-10}$. For the FN method, we used the R function rq.fit.fnb, which called compiled Fortran code (quantreg package, version 5.05).

In each generated sample, we estimated the regression coefficients with the GS and the FN algorithm. We compared the methods on the following measures:

1. (1) ‘FN-to-GS time ratio’: the ratio between mean CPU time of FN and mean CPU time of GS. In all models, the timing of GS included least-squares estimation for the initial values ${\beta }_p^0$. For models (2)–(4), it also included the standardization of all covariates and subsequent rescaling of the parameter estimates.

2. (2) ‘Condition on residuals met’: the proportion of samples in which the number of negative (R₋) and the number of positive (R₊) residuals satisfied the condition $R_{-}\le np\le n-R_{+}$.[14, Theorem 3.4]

3. (3) ‘Max ${\log }_{10}|l({\beta }_p^{\mathrm{GS}})-l({\beta }_p^{\mathrm{FN}})|$’: the maximum of the logarithm to base 10 of the absolute difference between the maximum log-likelihood values achieved in each sample.

4. (4) ‘Max $|{\beta }_p^{\mathrm{GS}}/{\beta }_p^{\mathrm{FN}}-1|$’: the maximum of the absolute relative difference of the parameter estimates in each sample.

5. (5) ‘Observed bias’: the difference between the estimated and the true parameter values averaged over the simulated samples.

The results of the simulation are given in Table 1. The observed bias was virtually zero at all quantiles in all the considered scenarios and it was not reported. The entries in Table 1 are averages over 300 samples comprising 100 replicates for each of the 3 distributions of the regression residual. No notable differences were observed across distributions. The performance of the two methods was comparable with respect to maximum likelihood achieved and proportion of sample in which the condition on residuals was met. The GS algorithm was faster than FN, although its advantage shrank as the number of covariates increased or the sample size decreased. Optimizing the R programming code and compiling it into machine language might further improve its computational speed.

Table 1.  Performance of the GS and FN algorithm for quantile regression estimation.

			Condition on residuals met
Sample		FN-to-GS
size	Quantile	time ratio	GS	FN	Max ${\log }_{10}|l({\beta }_p^{\mathrm{GS}})-l({\beta }_p^{\mathrm{FN}})|$	Max $|{\beta }_p^{\mathrm{GS}}/{\beta }_p^{\mathrm{FN}}-1|$
Model 1: yi=i, with $i\in \{1,\dots ,n\}$
n=10^4	p=0.50	9.59	1.00	1.00	−Inf	0.00
	p=0.75	9.96	1.00	1.00	−Inf	0.00
	p=0.90	10.09	1.00	1.00	−Inf	0.00
n=10^5	p=0.50	14.01	1.00	1.00	−Inf	0.00
	p=0.75	13.70	1.00	1.00	−11.44	0.00
	p=0.90	15.11	1.00	1.00	−Inf	0.00
Model 2: $y_i=1+x_i^{(1)}+x_i^{(2)}+u_i$
n=10^4	p=0.50	2.02	0.54	0.52	−7.16	0.00
	p=0.75	2.59	0.54	0.61	−6.89	0.00
	p=0.90	2.41	0.52	0.54	−6.20	0.01
n=10^5	p=0.50	3.98	0.39	0.58	−7.98	0.00
	p=0.75	4.59	0.43	0.52	−8.13	0.00
	p=0.90	4.13	0.46	0.57	−8.21	0.00
Model 3: $y_i=1+x_i^{(3)}+x_i^{(4)}+x_i^{(5)}+x_i^{(6)}+x_i^{(7)}+x_i^{(8)}+x_i^{(9)}+x_i^{(10)}+u_i$
n=10^4	p=0.50	0.96	0.35	0.40	−6.64	0.00
	p=0.75	1.06	0.38	0.41	−6.14	0.00
	p=0.90	1.04	0.37	0.35	−5.85	0.01
n=10^5	p=0.50	1.82	0.28	0.40	−8.25	0.00
	p=0.75	2.22	0.30	0.39	−7.78	0.00
	p=0.90	2.28	0.29	0.43	−7.55	0.00
Model 4: $y_i=-30+5000x_i^{(11)}+2x_i^{(12)}+0.5x_i^{(13)}+(0.5+x_i^{(11)}-0.5x_i^{(12)}+0.5x_i^{(13)})u_i$
n=10^4	p=0.50	1.01	0.42	0.53	−6.49	0.00
	p=0.75	1.96	0.50	0.47	−6.76	0.01
	p=0.90	2.28	0.52	0.48	−6.14	0.12
n=10^5	p=0.50	1.16	0.17	0.45	−7.30	0.00
	p=0.75	3.10	0.38	0.37	−8.27	0.00
	p=0.90	4.15	0.39	0.34	−6.91	0.01

Entries are averages over 300 simulated samples, 100 from each of 3 different distributions for the regression residual, u_i ∼ Normal(0,1), u_i ∼ T(3), and u_i ∼ Lognormal(0,1). The covariates are $x_i^{(j)}\sim$ Normal(0,1), $j\in \{1,\dots ,10\}$, $x_i^{(11)}\sim$ Uniform(0,1), $x_i^{(12)}\sim$ Bernoulli(0.5), and $x_i^{(13)}\sim$ T(3).

Following an insightful comment from an anonymous reviewer, we checked sensitivity of the algorithm with respect to its initial parameter values ${\beta }_p^0$. We re-ran the entire simulation study twice. The first time we set all elements of ${\beta }_p^0$ equal to 0; the second time we set them equal to 100. While GS took longer to converge with the unreasonable starting values than with the least-squares estimates, it was still generally faster than FN. The final estimates for the parameter β_p and achieved maximum likelihood seemed largely unaffected by the choice of the initial values and remained comparable to those reported in Table 1.

5. Remarks

The GS algorithm presented in this note is similar to the Newton–Raphson algorithm, in that it moves from the current parameter value in the direction of the gradient, and to direct search methods, in that it adjusts iteratively the step-length multiplier. Unlike methods based on vertex-search in a simplex, as Barrodale and Roberts’ method,[10] the algorithm proposed in this paper searches the parameter space without constraints. As Newton–Raphson and other optimization algorithms, GS only approximates the optimal solution within a given tolerance.

Perhaps one of the most useful features of the GS algorithm is its wide applicability. GS may be regarded as a general algorithm for unconstrained optimization of any objective function that is continuous, concave, and first-order differentiable everywhere except at a set of point with measure zero. Unlike Newton–Raphson, GS does not require a second derivative, which is generally unavailable in Laplace-based likelihoods. Thanks to its flexibility, the algorithm can be applied to the estimation of Laplace regression with censored data, mixed-effects quantile regression, and other models based on Laplace likelihood or the optimization of objective function with similar features.

We recommend setting the initial parameter, ${\beta }_p^0$, at reasonable values. While GS appears to be rather insensitive to the choice of the initial parameter value, reasonable starting values generally improve both convergence speed and accuracy. Least-squares estimates have been satisfactory in our experience. The choice of the convergence criteria and tolerance levels may substantially impact convergence. Stricter criteria may improve maximization accuracy but increase computing time. Generally, we recommend setting ε as small, a as close to 1, and b as large, as computational time and user's patience allow. The values used in the simulation study ($ϵ={10}^{-10}$, a=0.5, and b=1.25) gave a comparable performance with the FN method. These values have also proved adequate in a number of recent applications of the GS algorithm to the estimation of the coefficients of Laplace regression with censored data.[15–17among others]

The steps of the GS algorithm listed in Section 2 can be modified in several ways. The following are some examples. The convergence criteria in step 2.b can include checking change from the previous iteration in the parameter estimates or in the signs of the residuals. In step 2.b.i, the step-length can be set to the initial value, ${\delta }^{k+1}={\delta }^0$. The step-length multipliers, a and b, may depend on the sample size. When searching for a new parameter value along the direction of the gradient, the algorithm could find the root of the directional derivative by bisection or by considering the sign of the residuals. This would reduce the number of iterations and increase the computational burden at each iteration.

In the simulation reported in Table 1 and in a number of real-data applications, the GS algorithm has shown a remarkable computational speed, which may be a desirable feature when analysing large samples. In small samples (e.g. n=100), the GS, FN, and Barrodale–Roberts algorithms converge so quickly as to make any possible slight differences in computational time inconsequential from most practical standpoints.