Fast and Locally Adaptive Bayesian Quantile Smoothing using Calibrated Variational Approximations

Let $y_{i}=\theta_{i}+\varepsilon_{i}\quad(i=1,\dots,n)$ be a sequence model, where $y_{i}$ is an observation, $\theta_{i}$ is a true location and $\varepsilon_{i}$ is a noise. The estimate of $\ell_{1}$ trend filtering (Kim et al., 2009) is given by solving the optimization problem

$\displaystyle\hat{\theta}=\arg\min_{\theta}\|y-\theta\|_{2}^{2}+\lambda\|D^{(k+1)}\theta\|_{1},$

(1)

where $y=(y_{1},\dots,y_{n})^{\top}$, $\theta=(\theta_{1},\dots,\theta_{n})^{\top}$, $D^{(k+1)}$ is a $(n-k-1)\times n$ difference operator matrix of order $k+1$, and $\lambda>0$ is a tuning constant. Depending on the different order $k$, we can express various smoothing such as piecewise constant, linear, quadratic, and so forth (see e.g. Tibshirani, 2014). A fast and efficient optimization algorithm for the problem (1) was proposed by Ramdas and Tibshirani (2016).

Recently, Brantley et al. (2020) proposed quantile trend filtering, defined as the optimization problem

$\displaystyle\hat{\theta}_{p}=\arg\min_{\theta}\rho_{p}(y-\theta)+\lambda\|D^{(k+1)}\theta\|_{1},$

(2)

where $\lambda>0$ is a tuning constant and $\rho_{p}(\cdot)$ is a check loss function given by

$\displaystyle\rho_{p}(r)=\sum_{i=1}^{n}r_{i}\{p-1(r_{i}<0)\},\quad 0<p<1.$

(3)

To solve the problem (2), Brantley et al. (2020) proposed a parallelizable alternating direction method of multipliers (ADMM) algorithm, and also proposed the selection of smoothing parameters $\lambda$ using a modified criterion based on the extended Bayesian information criterion.

To conduct Bayesian inference of the quantile trend, we often use the following model:

$\displaystyle y_{i}=\theta_{pi}+\varepsilon_{i},\quad\varepsilon_{i}\sim\mathrm{AL}(p,\sigma^{2}),\quad i=1,\ldots,n,$

(4)

where $\theta_{pi}$ and $\sigma^{2}$ are unknown parameters, $p$ is a fixed quantile level, and ${\rm AL}(p,\sigma^{2})$ denotes the asymmetric Laplace distribution with the density function

$\displaystyle f_{\mathrm{AL}(p)}(x)=\frac{p(1-p)}{\sigma^{2}}\exp\left\{-\rho_{p}\left(\frac{x}{\sigma^{2}}\right)\right\},$

where $p$ is a fixed constant which characterizes the quantile level, $\sigma^{2}$ (not $\sigma$) is a scale parameter, and $\rho_{p}(\cdot)$ is a check loss function defined by (3). However, we often handle multiple observations per grid point in practice. Hereafter, we considered the following model which accounts for multiple observations per grid point:

$\displaystyle y_{ij}=\theta_{p}(x_{i})+\varepsilon_{ij},\quad\varepsilon_{ij}\sim\mathrm{AL}(p,\sigma^{2}),\ \ \ i=1,\ldots,n,\quad j=1,\dots,N_{i},$

(5)

where $\theta_{p}(x)$ is a $p$-th quantile in the location $x$, $n$ is the number of locations data are observed, and $N_{i}$ is the amount of data for each location $x_{i}$. It is a natural generalization of the sequence model (4) (see also Heng et al., 2022). Note that the model (5) is a nonparametric quantile regression with a single covariate, and the proposed approach can easily be generalized for additive regression with multiple covariates. For simplicity in notation, we dropped the subscript $p$ from $\theta$ for the remainder of the paper.

We then introduce shrinkage priors on differences. We define the $(k+1)$th order difference operator $D$ as

$\displaystyle D=\left(\begin{matrix}I_{k+1}&O\\ \lx@intercol\hfil D_{n}^{(k+1)}\hfil\lx@intercol\end{matrix}\right),$

where $I_{k+1}$ is $(k+1)\times(k+1)$ identity matrix, $O$ is zero matrix and $D_{n}^{(k+1)}$ is $(n-k-1)\times n$ standard difference matrix that is defined by

$\displaystyle D_{n}^{(1)}=\left(\begin{matrix}1&-1&&\\ &\ddots&\ddots&\\ &&1&-1\end{matrix}\right)\in\mathbb{R}^{(n-1)\times n},\quad D_{n}^{(k+1)}=D_{n-k}^{(1)}D_{n}^{(k)}.$

We consider flexible shrinkage priors on $D\theta$, and the priors are represented by

$\displaystyle D\theta\mid\tau^{2},\sigma^{2},w\sim N_{n}(0,\sigma^{2}W)\ \ \text{with}\ \ W={\rm diag}(w_{1}^{2},\ldots,w_{k+1}^{2},\tau^{2}w_{k+2}^{2},\dots,\tau^{2}w_{n}^{2}),$

where $w=(w_{1},\ldots,w_{n})$ represents local shrinkage parameters for each element in $D\theta$ and $\tau^{2}$ is a global shrinkage parameter. Since the $D$ is non-singular matrix, the prior of $\theta$ can be rewritten as

$\displaystyle\theta\mid\tau^{2},\sigma^{2},w\sim N_{n}(0,\sigma^{2}(D^{\top}W^{-1}D)^{-1}).$

(6)

Note that since $(D\theta)_{i}=\theta_{i}\sim N(0,w_{i}^{2})$ for $i=1,\dots,k+1$, $w_{i}$ ($i=1,\dots,k+1$) is not related to shrinkage of difference. For this reason, we assumed the conjugate inverse gamma distribution $\mathrm{IG}(a_{w_{i}},b_{w_{i}})$ for $w_{i}^{2}$. For $i=k+2,\dots,n$, we considered two types of distribution; $w_{i}\sim{\rm Exp}(1/2)$ and $w_{i}\sim C^{+}(0,1)$. These priors were motivated from Laplace or Bayesian lasso prior (Park and Casella, 2008) and horseshoe prior (Carvalho et al., 2010), respectively. Regarding the other parameters, we assigned $\sigma^{2}\sim{\rm IG}(a_{\sigma},b_{\sigma})$ and $\tau\sim C^{+}(0,C_{\tau})$. The default choice of hyperparameters is $a_{\sigma}=b_{\sigma}=0.1$ and $C_{\tau}=1$.

We first derived a Gibbs sampler by using the stochastic representation of the asymmetric Laplace distribution (Kozumi and Kobayashi, 2011). For $\varepsilon_{ij}\sim\mathrm{AL}(p,\sigma^{2})$, we have the following stochastic expression:

$$\varepsilon_{ij}=\psi z_{ij}+\sqrt{\sigma^{2}z_{ij}t^{2}}u_{ij},\quad\psi=\frac{1-2p}{p(1-p)},\quad t^{2}=\frac{2}{p(1-p)},$$

where $u_{ij}\sim N(0,1)$ and $z_{ij}\mid\sigma^{2}\sim\mathrm{Exp}(1/\sigma^{2})$ for $i=1,\ldots,n$ and $j=1,\dots,N_{i}$. From the above expression, the conditional likelihood function of $y_{ij}$ is given by

$\displaystyle p(y_{ij}\mid\theta_{i},z_{ij},\sigma^{2})=(2\pi t^{2}\sigma^{2})^{-1/2}z_{ij}^{-1/2}\exp\left\{-\frac{(y_{ij}-\theta_{i}-\psi z_{ij})^{2}}{2t^{2}\sigma^{2}z_{ij}}\right\}.$

(7)

Under the prior (6), the full conditional distributions of $\theta$ and $z_{i}$ are given by

	$\displaystyle\theta\mid y,z,\sigma^{2},\gamma^{2}\sim N_{n}\left(A^{-1}B,\sigma^{2}A^{-1}\right),$
	$\displaystyle z_{ij}\mid y_{ij},\theta_{i},\sigma^{2}\sim\mathrm{GIG}\left(\frac{1}{2},\frac{(y_{ij}-\theta_{i})^{2}}{t^{2}\sigma^{2}},\frac{\psi^{2}}{t^{2}\sigma^{2}}+\frac{2}{\sigma^{2}}\right),\ \ i=1,\ldots,n,\ j=1,\dots,N_{i},$

respectively, where

	$\displaystyle A$	$\displaystyle=D^{\top}W^{-1}D+\frac{1}{t^{2}}\mathrm{diag}\left(\sum_{j=1}^{N_{1}}z_{1j}^{-1},\ldots,\sum_{j=1}^{N_{n}}z_{nj}^{-1}\right),$
	$\displaystyle B$	$\displaystyle=\left(\frac{1}{t^{2}}\sum_{j=1}^{N_{1}}\left(\frac{y_{1j}}{z_{1j}}-\psi\right),\dots,\frac{1}{t^{2}}\sum_{j=1}^{N_{n}}\left(\frac{y_{nj}}{z_{nj}}-\psi\right)\right)^{\top}$

and $\mathrm{GIG}(a,b,c)$ denotes the generalized inverse Gaussian distribution. The full conditional distribution of the scale parameter $\sigma^{2}$ is given by

$\displaystyle\sigma^{2}\mid y,\theta,z,w,\tau^{2}\sim\mathrm{IG}\left(\frac{n+3N}{2}+a_{\sigma},\alpha_{\sigma^{2}}\right),$

where $N=\sum_{i=1}^{n}N_{i}$ is the number of observed values and

$\displaystyle\alpha_{\sigma^{2}}=\sum_{i=1}^{n}\sum_{j=1}^{N_{i}}\frac{(y_{ij}-\theta_{i}-\psi z_{ij})^{2}}{2t^{2}z_{ij}}+\frac{1}{2}\theta^{\top}D^{\top}W^{-1}D\theta+\sum_{i=1}^{n}\sum_{j=1}^{N_{i}}z_{ij}+b_{\sigma}.$

By using the augmentation of the half-Cauchy distribution (Makalic and Schmidt, 2015), the full conditional distribution of the global shrinkage parameter $\tau^{2}$ is given by

$\displaystyle\tau^{2}\mid\theta,w,\sigma^{2},\xi\sim\mathrm{IG}\left(\frac{n-k}{2},\frac{1}{2\sigma^{2}}\sum_{i=k+1}^{n}\frac{\eta_{i}^{2}}{w_{i}^{2}}+\frac{1}{\xi}\right),\quad\xi\mid\tau^{2}\sim\mathrm{IG}\left(\frac{1}{2},\frac{1}{\tau^{2}}+1\right),$

where $\xi$ is an augmented parameter for $\tau^{2}$. For $i=1,\dots,k+1$, we assumed the prior $\mathrm{IG}(a_{w_{i}},b_{w_{i}})$ for $w_{i}$, and then the full conditional distribution of $w_{i}$ is given by

$\displaystyle w_{i}^{2}\mid\theta,\sigma^{2}\sim\mathrm{IG}\left(\frac{1}{2}+a_{w_{i}},\frac{\eta_{i}^{2}}{2\sigma^{2}}+b_{w_{i}}\right),\ \ i=1,\dots,k+1.$

The full conditional distribution of local shrinkage parameter $w_{i}$ ($i=k+2,\dots,n$) depends on the choice of the prior, either Laplace or horseshoe prior.

• -

  (Laplace-type prior)   The full conditional distributions of $\theta$, $z_{i}$, and $\sigma^{2}$ were already derived. For Laplace-type prior, we set $\tau^{2}=1$ and assume $w_{i}\mid\gamma^{2}\sim\mathrm{Exp}(\gamma^{2}/2)$ for $i=k+2,\dots,n$. Then we have $(D\theta)_{i}\sim\mathrm{Lap}(\gamma)$. Noting that $\gamma\sim C^{+}(0,1)$, sampling from the standard half-Cauchy prior is equivalent to $\gamma^{2}\mid\nu\sim\mathrm{IG}(1/2,1/\nu)$ and $\nu\sim\mathrm{IG}(1/2,1)$ using the augmentation technique (Makalic and Schmidt, 2015). Hence, the full conditional distributions of $w_{i}$, $\gamma^{2}$ and $\nu$ are, respectively, given by

  	$\displaystyle w_{i}^{2}\mid\theta,\sigma^{2},\gamma^{2}\sim\mathrm{GIG}\left(\frac{1}{2},\frac{\eta_{i}^{2}}{\sigma^{2}},\gamma^{2}\right),$
  	$\displaystyle\gamma^{2}\mid w,\nu\sim\mathrm{GIG}\left(n-k-\frac{3}{2},\frac{2}{\nu},\sum_{i=k+2}^{n}w_{i}^{2}\right),\quad\nu\mid\gamma^{2}\sim\mathrm{IG}\left(\frac{1}{2},\frac{1}{\gamma^{2}}+1\right).$

• -

  (Horseshoe-type prior)   The full conditional distributions of $\theta$, $z_{i}$, $\sigma^{2}$ and $\tau^{2}$ were already derived. For Horseshoe-type prior, we assume $w_{i}\sim C^{+}(0,1)$ for $i=k+2,\dots,n$. By using the representation $w_{i}^{2}\mid\nu_{i}\sim\mathrm{IG}(1/2,1/\nu_{i})$ and $\nu_{i}\sim\mathrm{IG}(1/2,1)$, the full conditional distributions of $w_{i}$ and $\nu_{i}$ are, respectively, given by

  $\displaystyle w_{i}^{2}\mid\theta,\sigma^{2},\tau^{2},\nu_{i}\sim\mathrm{IG}\left(1,\frac{1}{\nu_{i}}+\frac{\eta_{i}^{2}}{2\sigma^{2}\tau^{2}}\right),\quad\nu_{i}\mid w_{i}\sim\mathrm{IG}\left(\frac{1}{2},\frac{1}{w_{i}^{2}}+1\right).$

The MCMC algorithm presented in Section 2.3 can be computationally intensive when the sample size is large. For the quick computation of the posterior distribution, we derived the variational Bayes approximation (e.g. Blei et al., 2017; Tran et al., 2021) of the joint posterior. The idea of the variational Bayes method is to approximate an intractable posterior distribution by using a simpler probability distribution. Note that the variational Bayes method does not require sampling from the posterior distribution like MCMC, and it searches for an optimal variational posterior by using the optimization method. In particular, we employed the mean-field variational Bayes (MFVB) approximation algorithms that require the forms of full conditional distributions as given in Section 2.3.

The variational distribution $q^{*}(\theta)\in\mathcal{Q}$ is defined by the minimizer of the Kullback-Leibler divergence from $q(\theta)$ to the true posterior distribution $p(\theta|y)$

$\displaystyle q^{*}=\arg\min_{q\in\mathcal{Q}}\mathrm{KL}(q\|p(\cdot\mid y))=\arg\min_{q\in\mathcal{Q}}\int q(\theta)\log\frac{q(\theta)}{p(\theta\mid y)}d\theta.$

(8)

If $\theta$ is decomposed as $\theta=(\theta_{1},\ldots,\theta_{K})$ and parameters $\theta_{1},\ldots,\theta_{K}$ are mutually independent, each variational posterior can be updated as

$\displaystyle q(\theta_{k})\propto\exp(\mathrm{E}_{\theta_{-k}}[\log p(y,\theta)])=\exp\left(\int q(\theta_{k})\log p(y,\theta)d\theta_{-k}\right),\ \ \ \ k=1,\ldots,K,$

where $\theta_{-k}$ denotes the parameters other than $\theta_{k}$ and $\mathrm{E}_{\theta_{-k}}[\cdot]$ denotes the expectation under the probability density given parameters except for $\theta_{k}$. Such a form of approximation is known as the MFVB approximation. If the full conditional distribution of $\theta_{k}$ has a familiar form, the above formula is easy to compute. According to the Gibbs sampling algorithm in Section 2.3, we used the following form of variational posteriors:

$\displaystyle q(\theta,z,\sigma^{2},\tau^{2},\xi)=q(\theta)q(z)q(\sigma^{2})q(\tau^{2})q(\xi),$

where

$\displaystyle\begin{split}&q(\theta)\sim N_{n}(A^{-1}B,(E_{1/\sigma^{2}}A)^{-1}),\ \ \ \ q(z_{ij})\sim\mathrm{GIG}\left(\frac{1}{2},\alpha_{z_{ij}},\beta_{z_{ij}}\right),\\ &q(\sigma^{2})\sim\mathrm{IG}\left(\frac{n+3N}{2}+a_{\sigma},\alpha_{\sigma^{2}}\right),\ \ \ \ q(\tau^{2})\sim\mathrm{IG}\left(\frac{n-k}{2},\alpha_{\tau^{2}}\right),\\ &q(\xi)\sim\mathrm{IG}\left(\frac{1}{2},E_{1/\tau^{2}}+1\right).\end{split}$

(9)

For $i=1,\dots,k+1$, we assume the prior $\mathrm{IG}(a_{w_{i}},b_{w_{i}})$ for $w_{i}$, and then the variational distribution of $w_{i}$ is given by

$\displaystyle q(w_{i}^{2})\sim\mathrm{IG}\left(\frac{1}{2}+a_{w_{i}},\frac{1}{2}E_{\eta_{i}^{2}}E_{1/\sigma^{2}}+b_{w_{i}}\right).$

The variational distributions of the other parameters depended on the specific choice of the distributional form of $\pi(w_{i})$ $(i=k+2,\dots,n)$, which are provided as follows:

• -

  (Laplace-type prior)   The variational distributions for $w_{i}^{2}$ ($i=k+2,\dots,n$), $\gamma^{2}$ and $\nu$ are given by

  	$\displaystyle q(w_{i}^{2})\sim\mathrm{GIG}\left(\frac{1}{2},\alpha_{w_{i}^{2}},E_{\gamma^{2}}\right),$
  	$\displaystyle q(\gamma^{2})\sim\mathrm{GIG}\left(n-k-\frac{3}{2},2E_{1/\nu},\sum_{i=k+2}^{n}E_{w_{i}^{2}}\right),\quad q(\nu)\sim\mathrm{IG}\left(\frac{1}{2},E_{1/\gamma^{2}}+1\right),$

• -

  (Horseshoe-type prior)   The variational distributions for $w_{i}^{2}$ and $\nu_{i}$ ($i=k+2,\dots,n$) are given by

  $\displaystyle q(w_{i}^{2})\sim\mathrm{IG}(1,\alpha_{w_{i}^{2}}),\quad q(\nu_{i})\sim\mathrm{IG}\left(\frac{1}{2},E_{1/w_{i}^{2}}+1\right),$

To obtain the variational parameters in each distribution, we update the parameters by using the coordinate ascent algorithm (e.g. Blei et al., 2017). The two proposed variational algorithms based on the above variational distributions are given in Algorithm 1 and 2. Note that we set $\epsilon=10^{-4}$ as the convergence criterion in the simulation study, $e_{i}$ is a unit vector that the $i$th component is 1, $d_{i}^{\top}$ is the $i$th row of difference matrix $D$, and $K_{c}(\cdot)$ is the modified Bessel function of the second kind with order $c$ in Algorithms 1 and 2.