Improving the local scoring algorithm using gradient sampling

TL;DR

This paper introduces an enhanced local scoring algorithm that employs gradient sampling to handle non-smooth optimization problems, improving stability and applicability to real-world data analysis tasks.

Contribution

The authors adapt the gradient sampling technique to the local scoring algorithm, enabling it to effectively optimize non-differentiable objective functions.

Findings

Improved numerical stability in local scoring algorithms.

Successful application to quantile regression and peaks-over-threshold methods.

Demonstrated effectiveness on retail and temperature datasets.

Abstract

We adapt the gradient sampling algorithm to the local scoring algorithm to solve complex estimation problems based on an optimization of an objective function. This overcomes non-differentiability and non-smoothness of the objective function. The new algorithm estimates the Clarke generalized subgradient used in the local scoring, thus reducing numerical instabilities. The method is applied to quantile regression and to the peaks-over-threshold method, as two examples. Real applications are provided for a retail store and temperature data analysis.

Equations

$$\underset{\mathbf{x}\in D}{\min }f(\mathbf{x})$$

$$\mathbf{x}\leftarrow \mathbf{x}+ts(\mathbf{d}),\text{mbox}or\mathbf{x}\leftarrow s(\mathbf{x}+ts(\mathbf{d})),$$

$$\bar{\partial }f(\mathbf{x})=\text{mbox}conv\left\{\underset{j}{\lim }\nabla f({\mathbf{z}}^j):{\mathbf{z}}^j\to \mathbf{x},{\mathbf{z}}^j\in D\right\},$$

$${\bar{\partial }}_{\epsilon }f(\mathbf{x})=\text{mbox}conv\left[\bar{\partial }f\{B(\mathbf{x},\epsilon )\}\right],$$

$$G_{\epsilon }(\mathbf{x})=\text{mbox}clconv\left[\nabla f\{B(\mathbf{x},\epsilon )\cap D\}\right],$$

$${\hat{G}}_{\epsilon }(\mathbf{x})=\text{mbox}conv\{\nabla f(\mathbf{x}),\nabla f({\mathbf{x}}_1),\dots ,\nabla f({\mathbf{x}}_m)\},$$

$${\hat{g}}_{\epsilon }(\mathbf{x})=\arg \min \left[\Vert g\Vert :g\in \text{mbox}conv\left\{{\hat{G}}_{\epsilon }(\mathbf{x})\right\}\right].$$

$$\sum _{i=1}^n{\rho }_{\alpha }\{q_{\alpha }({\mathbf{w}}_i);y_i\},$$

$${\rho }_{\alpha }(q;y)=(1-\alpha )(y-q{)}^{-}+\alpha (y-q{)}^{+},$$

$$\displaystyle\nabla\rho_{\alpha}(q;y)=\left\{\begin{array}[]{rl}1-\alpha,&\mbox{ if }y-q<0,\\ -\alpha,&\mbox{ if }y-q>0.\end{array}\right.$$

$$\ell (\sigma ,\kappa )=\sum _{i=1}^n-\log {\sigma }_i-\left(1+1/{\kappa }_i\right)\log \left(1+{\kappa }_i{y}_i/{\sigma }_i\right).$$

$${\sigma }_i=h\left\{{\beta }_0+\sum _{j=1}^k{f}_j(w_{ij})\right\},i=1,\dots ,n,$$

$${\theta }_i=\theta ({\mathbf{w}}_i)={\beta }_0^{\theta }+\sum _{j=1}^p{f}_j^{\theta }(w_{ij}),{\zeta }_i=\zeta ({\mathbf{w}}_i)={\beta }_0^{\zeta }+\sum _{j=1}^p{f}_j^{\zeta }(w_{ij}),$$

$$\theta_{i}=\left\{\begin{array}[]{lr}(c_{\alpha}^{-\kappa_{i}}-1)\sigma_{i}/\kappa_{i},&\kappa_{i}\neq 0,\\ -\sigma_{i}\log(c_{\alpha}),&\kappa_{i}=0.\end{array}\right.,\quad\zeta_{i}=\left\{\begin{array}[]{lr}(\theta_{i}+\sigma_{i})/(1-\kappa_{i}),&\kappa_{i}\neq 0,\\ \theta_{i}+\sigma_{i},&\kappa_{i}=0.\end{array}\right.$$

$${\nabla }_{\Theta }\ell =\frac{\partial \ell }{\partial \Theta }=\frac{\partial {\Lambda }^T}{\partial \Theta }\frac{\partial \ell }{\partial \Lambda }={\left(\frac{\partial {\Theta }^T}{\partial \Lambda }\right)}^{-1}\frac{\partial \ell }{\partial \Lambda }={\left({\nabla }_{\Lambda }\Theta \right)}^{-1}{\nabla }_{\Lambda }\ell .$$

$${\Lambda }^{\ast }=\Lambda ({\Theta }^{\ast })=\Lambda (\Theta +t\mathbf{d})\approx \Lambda +t{\left({\nabla }_{\Lambda }\Theta \right)}^{-1}\mathbf{d}.$$

$$q_{\alpha }(d_i,h_j)={\alpha }_0+{\gamma }_{ij}{d}_i{h}_j,$$

$$\theta_{i}^{\alpha_{1}}=\left\{\begin{array}[]{lr}(c_{\alpha_{1}}^{-\kappa_{i}}-1)\sigma_{i}/\kappa_{i},&\kappa_{i}\neq 0,\\ -\sigma_{i}\log(c_{\alpha_{1}}),&\kappa_{i}=0.\end{array}\right.,\quad\theta_{i}^{\alpha_{2}}=\left\{\begin{array}[]{lr}(c_{\alpha_{2}}^{-\kappa_{i}}-1)\sigma_{i}/\kappa_{i},&\kappa_{i}\neq 0,\\ -\sigma_{i}\log(c_{\alpha_{2}}),&\kappa_{i}=0.\end{array}\right.$$

$${\hat{g}}_{\epsilon }(\mathbf{x})=\arg \min \left[\Vert g\Vert :g\in \text{mbox}conv\left\{{\hat{G}}_{\epsilon }(\mathbf{x})\right\}\right].$$

$$f(\mathbf{x}+t\mathbf{d})<f(\mathbf{x})-\beta t{\mathbf{d}}^T\nabla f(\mathbf{x}).$$

$$\mathbf{z}=r_1{\mathbf{z}}_1+\cdots +r_m{\mathbf{z}}_l=Z\mathbf{r},\mathbf{r}⪰0,\sum _{j=1}^m{r}_j=1,$$

$$\underset{\mathbf{r}\in {\mathbb{R}}^m}{\min }$$

$$s.t.$$

$${\nabla }_{\Theta }\ell ({\Theta }^{\ast })=\frac{\partial \ell }{\partial \Theta }({\Theta }^{\ast })=\frac{\partial {\Lambda }^T}{\partial \Theta }({\Theta }^{\ast })\frac{\partial \ell }{\partial \Lambda }\left({\Lambda }^{\ast }\right)$$

Full text

Keywords: Gradient sampling; local scoring; peaks-over-threshold; quantile regression

1 Introduction

Among the statistical applications based on the optimization of an objective function (e.g., negative log-likelihood, risk function, etc.), more and more are based on a function that can be non-differentiable and/or non-convex. In addition, the parameter estimates are often constrained to have some predefined properties (e.g., being linear, additive in covariates, smooth, etc.). This paper introduces the use of the gradient sampling algorithm in such cases.

In essence, the gradient sampling algorithm, introduced by Burke et al. (2002), is a descent algorithm where the gradient direction is replaced by a stochastic approximation of the Clarke subdifferential. As a result, the descent algorithm is stabilized and can be used for non-convex and/or non-differentiable objective functions. To illustrate, in this paper, we propose the use of the gradient sampling approximation for the local scoring algorithm of Hastie & Tibshirani (1986) or for an extension of it, the generalized additive models for location scale and shape (GAMLSS) of Rigby & Stasinopoulos (2005).

In general, a descent algorithm aims at solving the optimization problem:

[TABLE]

with steps of the form $\mathbf{x}\leftarrow\mathbf{x}+t\mathbf{d}$, moving from the current point $\mathbf{x}$ to the next one. The direction $\mathbf{d}$ and the step size $t>0$ should both be selected such that $f$ is decreasing enough at each step. A common choice for $\mathbf{d}$ is the Newton step $-\{\nabla^{2}f(\mathbf{x})\}^{-1}\nabla f(\mathbf{x})$ and $t=1$. Although probably the best choice for convex second-differentiable functions, it is impossible to use in more complex situations, and we do not consider it here. The most common alternative choice is the gradient descent with $\mathbf{d}=-\nabla f(\mathbf{x})$, and $t$ being selected using a line search algorithm. A gradient descent can be applied to more general functions than the Newton descent. However it is much slower to converge.

The local scoring algorithm is a descent algorithm where the final solution should have some predefined properties, often smoothness and additivity in covariates. To achieve this, the descent direction $\mathbf{d}$ is projected onto a set with these properties at each step of the algorithm. The updating rule of the current solution $\mathbf{x}$ is of the form:

[TABLE]

where $s$ refers to the projection. Replacing the gradient descent direction $\mathbf{d}$ with the gradient sampling approximation looks natural here.

We present the introduction of the use of the gradient sampling algorithm in two cases. The first is an additive quantile regression. It provides an alternative algorithm to that of Portnoy & Koenker (1997). The second is a smooth peaks-over-threshold (POT) when the smoothness is imposed on the return levels at different levels instead of the canonical parameters of the generalized Pareto distribution underlying the method. It provides an alternative to GAMLSS, bringing an appreciable stabilization to the maximum likelihood algorithm in that complex setting.

The remainder of the paper is organized as follows. In Section 2, we introduce the gradient sampling descent algorithm in general, and combined with the local scoring algorithm for quantile additive models and for POT methods. Two real data cases one in retail store management and the other in the environment context are presented in Section 3. Section 4 consists of a discussion. Algorithms are described in details in the Appendix.

2 The gradient sampling descent algorithm (GSDA)

2.1 The algorithm

The gradient sampling algorithm was first introduced by Burke et al. (2002) and has been further developed by Burke et al. (2005), Kiwiel (2007), and Kiwiel (2010), for example. An up-to-date presentation can be found in Bagirov et al. (2014), chap. 13. The idea is to build a gradient descent algorithm where the descent direction is the Clarke subdifferiential. However, the latter cannot be computed exactly because one has to determine the subgradient set, which is feasible only in specific cases. Burke et al. (2002) have therefore proposed the gradient sampling algorithm to approximate it.

For the sake of completeness, below we quote some definitions following Kiwiel (2007), adapted to our statistical context. Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a locally Lipschitz continuous function, continuously differentiable111In the Lebesgue sense, that is the set of point where $f$ is non-continuously differentiable is of null Lebesgue measure. on an open dense set $D\in\mathbb{R}^{n}$. The Clarke subdifferential of $f$ at $\mathbf{x}\in\mathbb{R}^{n}$ is the set:

[TABLE]

where conv$(A)$ is the convex hull of $A$. A point $\mathbf{x}$ is called stationary for $f$ if $\mathbf{0}\in\bar{\partial}f(\mathbf{x})$. In particular, when $f$ is $C^{1}$, then $\bar{\partial}f(\mathbf{x})$ reduces to $\{\nabla f(\mathbf{x})\}$, and a stationary point satisfies $\nabla f(\mathbf{x})~{}=~{}0$. The Clarke $\varepsilon$-subdifferential is defined by:

[TABLE]

where $B(\mathbf{x},\varepsilon)$ is a ball of radius $\varepsilon\geq 0$ centered at $\mathbf{x}$. A point $\mathbf{x}$ is called $\varepsilon$-stationary for $f$ if $\mathbf{0}\in\bar{\partial}_{\varepsilon}f(\mathbf{x})$. For a closed convex set $G$, let $\mbox{Proj}(\mathbf{0}|G)$ be the minimum norm element of $G$, or, equivalently, the projection of $\mathbf{0}$ onto $G$. The Clarke generalized subgradient of $f$ at $\mathbf{x}$ is defined by $g_{\varepsilon}(\mathbf{x})=\mbox{Proj}\{\mathbf{0}|\bar{\partial}_{\varepsilon}f(\mathbf{x})\}$.

Because $g_{\varepsilon}(\mathbf{x})$ cannot be explicitly computed in general, the gradient sampling algorithm approximates it. Denoted by cl conv$(A)$, the closed convex hull of $A$, the Clarke $\varepsilon$-subdifferential $\bar{\partial}_{\varepsilon}f(\mathbf{x})$ is approximated by:

[TABLE]

since $G_{\varepsilon}(\mathbf{x})\subset\bar{\partial}_{\varepsilon}f(\mathbf{x})$, and $\bar{\partial}_{\varepsilon_{1}}f(\mathbf{x})\subset G_{\varepsilon_{2}}(\mathbf{x})$ for any $0\leq\varepsilon_{1}<\varepsilon_{2}$. To put that in action, the gradient sampling simulates an independent sample $\{\mathbf{x}_{i}\}_{i=1}^{m}$ in ${\cal B}(\mathbf{x},\varepsilon)$ to approximate $G_{\varepsilon}(\mathbf{x})$ by

[TABLE]

where $m\geq n+1$. Next, $g_{\varepsilon}(\mathbf{x})$ is approximated by $\hat{g}_{\varepsilon}(\mathbf{x})=\mbox{Proj}\{\mathbf{0}|\hat{G}_{\varepsilon}(\mathbf{x})\}$. Finally, the gradient descent algorithm uses the opposite of the Clarke generalized subgradient, $-\hat{g}_{\varepsilon}(\mathbf{x})$, as the descent direction.

Henceforth, we will refer to the Gradient Sampling Descent Algorithm using the initialism GSDA. The full algorithm is reported in Appendix A. The main step is the computation of $\hat{g}_{\varepsilon}(\mathbf{x})$, as follows

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $n$-vector of 0. 2. 2.

Set $\hat{G}_{\varepsilon}(\mathbf{x})=\{\nabla f(\mathbf{x}),\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{1}),\ldots,\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{m})\}$. 3. 3.

Find

[TABLE]

We can see that this step is fairly easy to implement in practice because it only requires to simulate in a ball of small radius around the current $\mathbf{x}$ and to solve a problem that happens to be quadratic and easy to solve.

Before detailing the GSDA further, let us illustrate how it works at a single step in Figure 1. Borrowed from Overton (2015), we consider the function $f(\mathbf{x})=10(x_{2}-x_{1}^{2})+(1-x_{1})^{2}$. The minimum (diamond) is at $(1,1)^{T}$. The bold gray line indicates the non-differentiability set $\Omega:=\{x_{2}=x_{1}^{2}\}$. The gradient descent $-\nabla f(\mathbf{x})$ is the dashed arrow. The short arrows are the sampled gradient descents $-\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{1}),\ldots,-\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{m})$. The gradient sampling descent $-\hat{g}_{\varepsilon}(\mathbf{x})$ is the solid bold arrow. The lengths of all the arrows have been modified to facilitate readability. The four plots show cases where the current solution $\mathbf{x}$ (big dot) is either above or below $\Omega$, and either far from it or close to it. In the two top plots, where $\mathbf{x}$ is far from $\Omega$, $\nabla f(\mathbf{x})$ and $\hat{g}_{\varepsilon}(\mathbf{x})$ are very close: the gradient sampling step is like a descent step when differentiability is guaranteed. In the two bottom plots, where $\mathbf{x}$ is close to $\Omega$, the gradient descent alternates up (right plot) and down (left plot), while the approximation descent is more robust to the current position of $\mathbf{x}$.

Therefore, in practice, using the gradient sampling at differentiable points close to non-differentiable area on one hand avoids the typical zigzag behavior of descent algorithms. On the other hand, using the gradient sampling comes at the cost of several gradient computations, making it inappropriate in large dimensions or complex cases without further adaptations.

To go into further details, the algorithm requires solving a sub-problem (2). This can be written as a quadratic problem under linear constraints that can be efficiently solved using classical quadratic programming e.g., function solveQP of R package quadprog (Turlach & Weingessel 2013). We have noted that, at certain iterations of the algorithm, the sub-problem (2) may be numerically unstable or difficult to compute. In such case, $\hat{g}_{\varepsilon}(\mathbf{x})$ can be conveniently replaced with the average of $\{\nabla f(\mathbf{x}),\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{1}),\ldots,,\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{m})\}$. Although no formal proof exists, this is intuitive because any stable vector pointing toward $\hat{G}_{\varepsilon}(\mathbf{x})$ could be used as an approximation of $g_{\varepsilon}(\mathbf{x})$. The details for solving (2) are reported in Appendix A.

In the two next sections, we show how to apply the association between GSDA and the local scoring algorithm to quantile additive models and to POT methods. In both cases, it consists in using a step as in (1), where $\mathbf{d}$ is $-\hat{g}_{\varepsilon}(\mathbf{x})$.

2.2 GSDA for quantile additive models

Quantile regression is now a well-known technology that has been introduced in detail, for example in Koenker (2005). An implementation in R is available from package quantreg with the function rqss in the context of additive models for quantile regression (Koenker 2016). Of note, it implements a Frisch–Newton interior point method (Portnoy & Koenker 1997). We now show how using the GSDA provides a simple alternative algorithm.

Let $y_{1},\ldots,y_{n}$ be $n$ observations, independent conditional on $k$-dimensional covariates $\mathbf{w}_{1},\ldots,\mathbf{w}_{n}$, where $\mathbf{w}_{i}=(w_{i1},\ldots,w_{ik})^{T}$. In its sample version, the quantile additive model aims at finding $q_{\alpha}(\mathbf{w}_{i})=\alpha_{0}+\sum_{j=1}^{k}f_{j}(w_{ij})$, a quantile function at level $0<\alpha<1$, minimizing:

[TABLE]

where the risk function is:

[TABLE]

$z^{+}=\max\{z,0\}$, and $z^{-}=-\min\{z,0\}$. As for any additive model, the functions $f_{j}$ need to satisfy some identifiability constraints $\sum_{i}f_{j}(w_{ij})=0$. At $y\neq q$, the derivative of the risk function is:

[TABLE]

The non-differentiability set of $\rho_{\alpha}$ makes the classical gradient descent unstable to use in a local scoring. The GSDA offers a conceptually very simple alternative. From the current $\mathbf{q}_{\alpha}=\{q_{\alpha}(\mathbf{w}_{1}),\ldots,q_{\alpha}(\mathbf{w}_{n})\}^{T}$ , the algorithm moves along $-\hat{g}_{\varepsilon}(\theta)$ after a smoothing operation. The main steps are as follows:

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $n$-vector of 0. 2. 2.

Set $\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})=\left\{\nabla\rho_{\alpha}(\mathbf{q}_{\alpha})+\sum_{k=1}^{m}\nabla\rho_{\alpha}(\mathbf{q}_{\alpha}+\varepsilon\mathbf{u}_{k})\right\}/(m+1)$.

Then

Set $\mathbf{d}^{*}=-{\tt gam}(\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\sim{\tt lo}(w_{1})+\ldots+{\tt lo}(w_{k}))$ and $\mathbf{d}=\mathbf{d}^{*}/\|\mathbf{d}^{*}\|$. 2. 4.

Update $\mathbf{q}_{\alpha}\leftarrow\mathbf{q}_{\alpha}+t\mathbf{d}$, where $t$ is appropriately selected.

See Appendix A for details, particularly the selection of $t$. The calculation of $\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})$ can be replaced with the full program (2). The initialization of $\mathbf{q}_{\alpha}$ can be a constant quantile or a more sophisticated estimate. For the smoothing part $-{\tt gam}(\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\sim{\tt lo}(w_{1})+\ldots+{\tt lo}(w_{k}))$, we mimic standard R coding and refer the interested reader to Wood (2006) for more details. Note that the smoothing part could be replaced by ${\tt lm}(\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\sim w_{1}+\cdots+w_{n})$, for example. This would provide a quantile regression model.

2.3 GSDA for peaks-over-threshold method

The POT method, first proposed by Davison & Smith (1990), is now a standard approach used in extreme value analysis. It has been further developed and refined in many directions. For a review of the methodologies proposed in the non-stationary cases, including additive models, see, for example, Chavez-Demoulin & Davison (2005).

The POT method assumes that observed excesses, $y_{1},\ldots,y_{n}$, above a sufficiently high threshold $u$ are independent conditionally on a set of $k$-covariates, $\mathbf{w}_{1},\ldots,\mathbf{w}_{n}$, and are distributed according to a generalized Pareto distribution (GPD) with scale $\sigma_{i}=\sigma(\mathbf{w}_{i})$ and shape $\kappa_{i}=\kappa(\mathbf{w}_{i})$. The log-likelihood function is:

[TABLE]

The domain of definition is $1+\kappa_{i}y_{i}/\sigma_{i}>0$, and $\sigma_{i}>0$, for $i=1,\ldots,n$. We may assume an additive predictor linked to the parameter such as e.g. Chavez-Demoulin & Davison (2005),

[TABLE]

where $h(\cdot)$ is the link function, $\exp(\cdot)$ in this case, and $f_{j}$ are smooth functions satisfying identifiability constraints, $\sum_{i}f_{j}(w_{ij})=0$. More generally, Rigby & Stasinopoulos (2005) consider any kind of reasonable link functions. In many applications, inference for both the so-called value-at-risk and the expected shortfall is of main importance. For a level $\alpha$, consider:

[TABLE]

where the value-at-risk $\theta_{i}\equiv\theta(x_{i})$ and the expected shortfall $\zeta_{i}\equiv\zeta(x_{i})$, both at level $1-\alpha$, are respectively:

[TABLE]

Typically, $1-\alpha=99\%$. The scale factor is $c_{\alpha}=\alpha/\mbox{Pr}(Y>u)$, and $\zeta_{i}$ is only defined for $\kappa_{i}<1$. In the environment context, the value-at-risk is called return level. The interest is in modelling return levels at two different levels, as we illustrate in Section 3, for the U.S. maximal temperatures.

The GSDA can be used to fit the POT model with value-at-risk additive in the covariates, because it will bring some appreciable numerical stability. For this, the log-likelihood and its derivatives are computed at each step. Denoting $\Theta=(\theta^{T},\zeta^{T})^{T}$ and $\Lambda=(\eta^{T},\kappa^{T})^{T}$, the derivatives can be conveniently computed with:

[TABLE]

This provides an explicit formula (i.e., computationally tractable) for the gradient of $\ell$ in $\Theta$. It should be noted at that point that computing $\ell$ and its derivatives given $\Lambda$ is easy, while it is complex given $\Theta$ since one needs to inverse $\Theta(\Lambda)$.

At each step, $\Theta$ moves along $\mathbf{d}=-s\{\hat{g}_{\varepsilon}(\Theta)\}$, where $s$ is a suitable transformation for matching the constraint. However, moving in the $\Theta$-space is painful because computing likelihoods and their derivatives is computationally difficult. Instead, we can move in the $\Lambda$-space, using a first-order approximation,

[TABLE]

Overall, the main steps in the algorithm are as follows:

Compute $M=\nabla_{\Lambda}\Theta$. 2. 2.

Compute the approximation $\hat{g}_{\varepsilon}(\Theta)=\left(\hat{g}_{\varepsilon,\theta}^{T}(\Theta),\hat{g}_{\varepsilon,\zeta}^{T}(\Theta)\right)^{T}$.

1. (a)

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $2n$-vector of 0. 2. (b)

Set $\hat{g}_{\varepsilon}(\Theta)=M^{-1}\left\{\nabla_{\Lambda}\ell(\Lambda)+\sum_{k=1}^{m}\nabla_{\Lambda}\ell(\Lambda+\varepsilon\mathbf{u}_{k})\right\}/(m+1)$ .

Then,

Set $\mathbf{d}^{*}_{\theta}=-{\tt gam}(\hat{g}_{\varepsilon,\theta}(\Theta)\sim{\tt lo}(w_{1})+\ldots+{\tt lo}(w_{k})).$ 2. 4.

Set $\mathbf{d}^{*}_{\zeta}=-{\tt gam}(\hat{g}_{\varepsilon,\zeta}(\Theta)\sim{\tt lo}(w_{1})+\ldots+{\tt lo}(w_{k})).$ 3. 5.

Set $\mathbf{d}^{*}=((\mathbf{d}^{*}_{\theta})^{T},(\mathbf{d}^{*}_{\zeta})^{T})^{T}$, $\mathbf{d}=\mathbf{d}^{*}/\|\mathbf{d}^{*}\|$. 4. 6.

Update $\Lambda\leftarrow\Lambda+tM^{-1}\mathbf{d}$, where $t$ is appropriately selected.

The complete algorithm is reported in Appendix A, in particular some justification for step 2(b), which is an approximation.

3 Real data

In this section, we illustrate the use of the gradient sampling algorithm in practice. The first example illustrates the use of the gradient sampling algorithm to estimate a semi-parametric quantile for alimentary products of an European retail store. The second shows the estimation of semi-parametric POT return levels at two different levels in the context of temperature data.

3.1 European retail store

Shelf replenishment is a large area of operations research (Khouja & Goyal 2008). Few works study intradaily sales, although in highly frequented stores some shelves of alimentary products need to be replenished several times a day. This is the case for the European retail store we consider, which is located in the railway station of a big city in Europe. The store is open every day of the week from 6 am to 11 pm. The data consist of hourly sales of three different products from November 1, 2012 to November 23, 2014, that is 743 days times 17 hours. The products are “Butter croissants”, an “Energy drink”, and “Milk” (one liter). The respective daily sales over the period are shown in Figures 2, 3, and 4. The patterns are different from one product to the other, and Sunday has, in all cases, a clear specific pattern. The croissants are sold mainly early in the morning, with a decreasing trend throughout the day, although a small sales peak occurs again around 5 pm. The energy drinks are basically sold constantly throughout day (except on Sunday), with some peaks appearing at breakfast and lunch time. The sales are much more important on Sunday, especially in the afternoon. The milk sales show an increasing trend until early in the evening for all days. If the store wants to guarantee, at a level of 90%, that the product shelf is never empty at any time of the day, a quantity of interest is the 90% quantile of sales. For each product, we use the gradient sampling algorithm given in Section 2.2 to estimate a 90%-quantile additive model of the form:

[TABLE]

where $\alpha=90\%$, $d_{i}$ represents the $i^{th}$-day of the week and $h_{j}$ the $j^{th}$-hour of the day, and $j=6,\ldots,23$. The algorithm quickly converges and the results are the black points in Figure 2 for the butter croissants, Figure 3 for the energy drink, and Figure 4 for the one-liter containers of milk. The interaction between day of the week and hour of the day allows the 90%–quantile estimate to capture the intradaily and intraweek patterns.

3.2 U.S. maximal temperature

We investigate the global change of hot days over the period 1950–2004 at three different stations in the U.S. The data were extracted from the U.S Cooperative Observer Program (COOP) and consist of daily temperature maxima measured over the period 1950–2004 at different stations in Alabama, California, and Colorado. The dataset can be freely downloaded from the Climate and Global Dynamics Division (2010). Our aim is to detect the presence of non-stationarity in the return level at two different levels $1-\alpha_{1}=95\%$, and $1-\alpha_{2}=99\%$, and its dependence on the station. We detrend the observations of each station by applying a smoothing spline through the 20089 data points and removing the smoothed mean to the observations. Figure 5 shows the daily maxima (points) observed at the three sites in 1963 (in grey) and 2000 (black). The line represents the corresponding estimated smoothing trend. The smoothed mean corresponding to year 2000 seems slightly below the one of year 1963. This is later confirmed by the return levels estimates.

The POT method is the standard quantitative risk methodology when very high quantile calculations are required. Following the POT approach, we fix a constant threshold for each site that corresponds to the 90%-quantile ot the detrended data and pool together the resulting three site series each of size 2009 over threshold. We fit the GPD return level models as in Section 2.3 with:

[TABLE]

where $\theta_{i}^{\alpha_{l}}=\theta_{i}^{\alpha_{l}}(s_{i},t_{i})$ $l=1,2$ correspond respectively to the $1-\alpha_{1}=0.95$ and $1-\alpha_{2}=0.99$ return levels depending on covariates $s_{i}$ specifying the site and $t_{i}$ the year of the $i$th value. For simplicity and because of the choice of the threshold, we set $c_{\alpha_{l}}=\alpha_{l}/0.01$, although the probability of exceeding the threshold in the denominator could have been estimated by maximum likelihood estimation, which would have led to a value close to 0.01. The estimated return levels are shown in Figure 6. The estimated model uses a GAM with 10 degrees of freedom for the time covariate and an indicator for the site. It shows a global decreasing trend in the sizes of the detrended data above threshold, with higher return levels at 95% and 99% for Colorado. The peak of return levels early in the 1960s is in accordance with the findings of Meehl et al. (2009). The fact that Colorado has a higher return level of size over threshold can be explained by the greater variability of its exceedances over threshold, probably due to altitude. The ability of the methodology to allow a simultaneous estimation of two flexible return levels within the algorithm offers an important clarification in terms of interpretation. First, the two returns levels are likely to behave similarly at different levels. Second, their time-varying structure is directly guided by the data itself and does not resulting from the time-varying forms of the GPD parameters estimated, as would occur in other standard non-stationary models.

4 Discussion

In this paper, the gradient sampling algorithm of Burke et al. (2002) is adapted with the local scoring to solve quite complex estimation problems based on the optimization of an objective function (a log-likelihood and a risk function) under some constraints of additivity and smoothness of the parameters. The aim is to overcome non-differentiability and/or non-smoothness of the objective function. In addition, the algorithm can be used for non-convex objective functions at the cost of local convergence.

The algorithm is very flexible and easy to implement. It can be used where a gradient descent could be used if the objective function were differentiable. It can also be used for differentiable functions to reduce numerical instability. The cost is computational because it requires a very large number of calculations of the gradient. This probably precludes its use in several high-dimensional cases. At the time of writing, we are not aware of any improvement to that aspect, although we can be confident that hybrid algorithms will soon bring an efficient solution.

The algorithm was adapted for a quantile additive model, where it offers an alternative algorithm to the Frisch–Newton interior point method of Portnoy & Koenker (1997), and for the POT model with parameters being smooth and additive in the covariates. Because the final results are maximum likelihood estimators, we have not included any confidence interval calculations, as these are not any different from any that would be obtained with another optimization method.

Even with its drawbacks, the gradient sampling represents a convenient alternative to existing algorithms in many situations where optimization is challenging, and even the only solution in many cases. Therefore, it should be part of the statistician’s toolbox.

Acknowledgements

This work was supported by the Swiss National Science Foundation.

Appendix A Appendix

A.1 Gradient sampling descent algorithm

The gradient sampling algorithm of Kiwiel (2007) is reported below, borrowed from Overton (2015).

Fix the sampling size $m\geq n+1$, a line search parameter $0<\beta<1$, and the reduction factors $\mu\in(0,1)$ and $\lambda\in(0,1)$. 2. 2.

Initialize the solution $\mathbf{x}$, the radius $\varepsilon>0$ and the tolerance $\tau>0$. 3. 3.

Compute the approximation $\hat{g}_{\varepsilon}(\mathbf{x})$.

1. (a)

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $n$-vector of 0. 2. (b)

Set $\hat{G}_{\varepsilon}(\mathbf{x})=\{\nabla f(\mathbf{x}),\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{1}),\ldots,\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{m})\}.$ 3. (c)

Find

[TABLE] 4. 4.

If $\|\hat{g}_{\varepsilon}(\mathbf{x})\|\leq\tau$, update $\varepsilon\leftarrow\mu\varepsilon$ and $\tau\leftarrow\lambda\tau$, go to 3. 5. 5.

If $\|\hat{g}_{\varepsilon}(\mathbf{x})\|>\tau$, do a backtracking line search along $\mathbf{d}=-g_{\varepsilon}(\mathbf{x})/\|g_{\varepsilon}(\mathbf{x})\|$, diminishing $t\in\{1,1/2,1/4,\ldots\}$ until the Armijo’s condition is satisfied,

[TABLE] 6. 6.

Update $\mathbf{x}\leftarrow\mathbf{x}+t\mathbf{d}$ and go to 3.

The overall algorithm is stopped when one or several convergence criteria are satisfied, typically, when $\tau$ and $\varepsilon$ reach a predefined threshold.

To solve sub-problem (7), recall that any element of $\mbox{conv}(G)$ can be written as a unique convex combination of vectors at the edges of $\mbox{conv}(G)$, say $G_{e}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{m}\}$. In more details, let $\mathbf{z}\in\mbox{conv}(G)$, and then there exists an unique $\mathbf{r}$ such that:

[TABLE]

where $Z$ is the matrix with columns $\{\mathbf{z}_{1},\ldots,\mathbf{z}_{l}\}$. Eddy (1977) shows how to derive the columns of $Z$. The method is efficiently implemented in chull of R. Because there exits $\mathbf{r}^{*}$ such that $g=Z\mathbf{r}^{*}$ for each $g\in\mbox{conv}\left\{\hat{G}_{\varepsilon}(\mathbf{x})\right\}$, the sub-problem (7) can be written as the quadratic problem under linear constraints:

[TABLE]

This minimization problem can be efficiently solved using classical quadratic programming (e.g., function solveQP of R package quadprog; Turlach & Weingessel (2013)). To further simplify this step, we have noted that, at certain iterations of the algorithm, the above problem may be numerically unstable or difficult to compute. In such case, $\hat{g}_{\varepsilon}(\mathbf{x})$ can be conveniently replaced by the average $(m+1)^{-1}\{\nabla f(\mathbf{x})~{}+~{}\sum_{j=1}^{m}\nabla f(\mathbf{x}+\varepsilon\mathbf{u}_{j})\}$. Although we have not found a formal proof, this is intuitive since any stable vector pointing toward $\hat{G}_{\varepsilon}(\mathbf{x})$ could be used as an approximation of $g_{\varepsilon}(\mathbf{x})$.

Finally, the Armijo’s conditions $(\ref{linesearch})$ can be replaced by $f(\mathbf{x}+t\mathbf{d})<f(\mathbf{x})-\beta t\|\hat{g}_{\varepsilon}(\mathbf{x})\|$. Indeed, as an approximation of $\nabla f(\mathbf{x})$, $\hat{g}_{\varepsilon}(\mathbf{x})$ can be used in the right-hand side of $(\ref{linesearch})$, giving $\mathbf{d}^{T}\hat{g}_{\varepsilon}(\mathbf{x})=\|\hat{g}_{\varepsilon}(\mathbf{x})\|$, see, for instance, Overton (2015).

A.2 GSDA for quantile additive models

The complete algorithm when applied to quantile additive models is given below.

Fix the sampling size $m\geq n+1$, a line search parameter $0<\beta<1$, and the reduction factors $\mu\in(0,1)$ and $\lambda\in(0,1)$. 2. 2.

Initialize the radius $\varepsilon>0$ and the tolerance $\tau>0$. 3. 3.

Initialize the solution $\mathbf{q}_{\alpha}$. 4. 4.

Compute the approximation $\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})$.

1. (a)

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $n$-vector of 0. 2. (b)

Set $\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})=(m+1)^{-1}\left\{\nabla\rho_{\alpha}(\mathbf{q}_{\alpha})+\sum_{k=1}^{m}\nabla\rho_{\alpha}(\mathbf{q}_{\alpha}+\varepsilon\mathbf{u}_{k})\right\}.$ 5. 5.

If $\|\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\|\leq\tau$, update $\varepsilon\leftarrow\mu\varepsilon$ and $\tau\leftarrow\lambda\tau$, go to 4. 6. 6.

If $\|\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\|>\tau$.

1. (a)

Set $\mathbf{d}^{*}=-{\tt gam}(\hat{g}_{\varepsilon}(\mathbf{q}_{\alpha})\sim{\tt lo}(w_{1})+\ldots+{\tt lo}(w_{k})).$ and $\mathbf{d}=\mathbf{d}^{*}/\|\mathbf{d}^{*}\|$. 2. (b)

Diminish $t\in\{1,1/2,1/4,\ldots\}$ until $\rho_{\alpha}(\mathbf{q}_{\alpha}+t\mathbf{d})<\rho_{\alpha}(\mathbf{q}_{\alpha})-\beta t\mathbf{d}^{T}\nabla\rho_{\alpha}(\mathbf{q}_{\alpha})$. 7. 7.

Update $\mathbf{q}_{\alpha}\leftarrow\mathbf{q}_{\alpha}+t\mathbf{d}$ and go to 4.

As in Section A.1, the algorithm can be adapted, for example, replacing the average in step 4(b) with the solution to sub-problem (7).

A.3 GSDA for POT

The complete algorithm when applied to POT is given below.

Fix the sampling size $m\geq 2n+1$, a line search parameter $0<\beta<1$, and the reduction factors $\mu\in(0,1)$ and $\lambda\in(0,1)$. 2. 2.

Initialize the radius $\varepsilon>0$ and the tolerance $\tau>0$. 3. 3.

Initialize the solution $\Lambda=(\eta^{T},\kappa^{T})^{T}$, $\Theta=\Theta(\Lambda)=(\theta^{T},\zeta^{T})^{T}$. 4. 4.

Compute $M=\nabla_{\Lambda}\Theta$. 5. 5.

Compute the approximation $\hat{g}_{\varepsilon}(\Theta)=\left(\hat{g}_{\varepsilon,\theta}^{T}(\Theta),\hat{g}_{\varepsilon,\zeta}^{T}(\Theta)\right)^{T}$.

1. (a)

Sample $\mathbf{u}_{1},\ldots,\mathbf{u}_{m}$ on the unit ball ${\cal B}(\mathbf{0},1)$, where $\mathbf{0}$ is the $2n$-vector of 0. 2. (b)

Set $\hat{g}_{\varepsilon}(\Theta)=M^{-1}\left\{\nabla_{\Lambda}\ell(\Lambda)+\sum_{k=1}^{m}\nabla_{\Lambda}\ell(\Lambda+\varepsilon\mathbf{u}_{k})\right\}/(m+1).$ 6. 6.

If $\|\hat{g}_{\varepsilon}(\Theta)\|\leq\tau$, update $\varepsilon\leftarrow\mu\varepsilon$ and $\tau\leftarrow\lambda\tau$, go to 4. 7. 7.

If $\|\hat{g}_{\varepsilon}(\Theta)\|>\tau$.

1. (a)

Set $\mathbf{d}^{*}_{\theta}=-{\tt gam}(\hat{g}_{\varepsilon,\theta}(\Theta)\sim{\tt lo}(x_{1})+\ldots+{\tt lo}(x_{k})).$ 2. (b)

Set $\mathbf{d}^{*}_{\zeta}=-{\tt gam}(\hat{g}_{\varepsilon,\zeta}(\Theta)\sim{\tt lo}(x_{1})+\ldots+{\tt lo}(x_{k})).$ 3. (c)

Set $\mathbf{d}^{*}=((\mathbf{d}^{*}_{\theta})^{T},(\mathbf{d}^{*}_{\zeta})^{T})^{T}$, $\mathbf{d}=\mathbf{d}^{*}/\|\mathbf{d}^{*}\|$. 4. (d)

Diminish $t\in\{1,1/2,1/4,\ldots\}$ until $\ell(\Lambda+tM^{-1}\mathbf{d})<\ell(\Lambda)-\beta tM^{-1}\mathbf{d}^{T}\nabla_{\Lambda}\ell(\Lambda)$. 8. 8.

Update $\Lambda\leftarrow\Lambda+tM^{-1}\mathbf{d}$, $\Theta=\Theta(\Lambda)$, and go to 4.

Focusing on step 5(b), it should be noted that $\hat{g}_{\varepsilon}(\Theta)$ should be the average of vectors like $\nabla_{\Theta}\ell(\Theta^{*})$ where $\Theta^{*}=\Theta+\varepsilon\mathbf{u}$, $\mathbf{u}$. We use the following approximation

[TABLE]

It means that we admit a common matrix $M$ although it could be computed at each sampled points. In addition, given that the distribution of $\mathbf{u}$ on the unit ball ${\cal B}(\mathbf{0},1)$ is not crucial to the algorithm, a further simplification can be obtained by replacing $\left\{\frac{\partial\Theta^{T}}{\partial\Lambda}(\Lambda)\right\}^{-1}\mathbf{u}$ by $\mathbf{u}$. Finally, like for the general GSDA, the average in step 5(b) could be replaced by the solution to sub-problem (7).