Efficient Bayesian shape-restricted function estimation with constrained Gaussian process priors

TL;DR

This paper introduces an efficient Bayesian method for shape-restricted function estimation using constrained Gaussian process priors, leveraging basis coefficient formulations and circulant embedding for improved sampling efficiency.

Contribution

It presents a novel sampling strategy that incorporates shape constraints via basis coefficients and circulant embedding, reducing computational complexity in Bayesian shape-restricted inference.

Findings

Method achieves high accuracy in simulated data

Algorithm demonstrates efficiency on real data examples

Effectively mitigates computational costs of hyperparameter updates

Abstract

This article revisits the problem of Bayesian shape-restricted inference in the light of a recently developed approximate Gaussian process that admits an equivalent formulation of the shape constraints in terms of the basis coefficients. We propose a strategy to efficiently sample from the resulting constrained posterior by absorbing a smooth relaxation of the constraint in the likelihood and using circulant embedding techniques to sample from the unconstrained modified prior. We additionally pay careful attention to mitigate the computational complexity arising from updating hyperparameters within the covariance kernel of the Gaussian process. The developed algorithm is shown to be accurate and highly efficient in simulated and real data examples.

Figures

Equations

$$\displaystyle p(\xi)\propto\,\,\exp\bigg{\{}-\frac{1}{2\sigma^{2}}\|Z-\Phi\xi\|^{2}\bigg{\}}\,\exp\bigg{\{}-\frac{1}{2\tau^{2}}\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi\bigg{\}}\,\,\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi),\quad\xi\in\mathbb{R}^{N},$$

$$y_i=f(x_i)+{ϵ}_i,{ϵ}_i\sim \text{mbox}N(0,{\sigma }^2),i=1,\dots ,n$$

$$\displaystyle\mathcal{C}_{\xi}:\,=\mathcal{C}_{\xi}^{N}=\bigg{\{}\xi\in\mathbb{R}^{N}:~{}\xi_{j}\geq 0,\ j=1,\dots,N\bigg{\}}.$$

$$\pi (\xi )\propto \mathcal{N}(\xi ;0,{\tau }^{2}K)\text{mathbbm}{1}_{{\mathcal{C}}_{\xi }}(\xi ),$$

$$\mathcal{N}_{N}\Big{(}(\Phi^{\mathrm{\scriptscriptstyle T}}\Phi/\sigma^{2}+K^{-1}/\tau^{2})^{-1}\Phi^{\mathrm{\scriptscriptstyle T}}Y,(\Phi^{\mathrm{\scriptscriptstyle T}}\Phi/\sigma^{2}+K^{-1}/\tau^{2})^{-1}\Big{)},$$

$$p(\xi )\propto L(\xi )\mathcal{N}(\xi ;0,\Sigma )$$

$${\xi }^{\prime }=\rho {\nu }_e+\sqrt{1-{\rho }^2}\xi ,{\nu }_e\sim \mathcal{N}(0,\Sigma )$$

$${\xi }^{\prime }={\nu }_e\sin \theta +\xi \cos \theta$$

$$\text{mathbbm}{1}_{{\mathcal{C}}_{\xi }}(\xi )\approx {\mathbb{J}}_{\eta }(\xi )=\prod _{j=1}^N\frac{e^{\eta {\xi }_j}}{1+e^{\eta {\xi }_j}}.$$

$$\tilde{p}(\xi \mid -)$$

$$Z=\bigg{(}Z\big{(}0\big{)},Z\bigg{(}\frac{1}{m}\bigg{)},Z\bigg{(}\frac{2}{m}\bigg{)},\cdots,Z\bigg{(}\frac{m-1}{m}\bigg{)}\bigg{)}^{\mathrm{\scriptscriptstyle T}}$$

$$\bigg{\{}0,\frac{1}{m},\frac{2}{m},\cdots,\frac{m-1}{m}\bigg{\}}\hskip 5.69054pt\equiv\hskip 5.69054pt\bigg{\{}u_{j}:u_{j}=\frac{j}{m};0\leq j<1\bigg{\}},\hskip 8.53581ptm\geq 1$$

$$G=\begin{bmatrix}\gamma(0)&\gamma\big{(}\frac{1}{m}\big{)}&\cdots&\gamma\big{(}\frac{m-1}{m}\big{)}\\ \gamma\big{(}\frac{1}{m}\big{)}&\gamma(0)&\cdots&\gamma\big{(}\frac{m-1}{m}\big{)}\\ \vdots&\vdots&\ddots&\vdots\\ \gamma\big{(}\frac{m-1}{m}\big{)}&\gamma\big{(}\frac{m-2}{m}\big{)}&\cdots&\gamma(0)\end{bmatrix}$$

$$Tu=-(r_1,\dots ,r_{M-1}{)}^{\mathrm{T}}$$

$$T_h{u}^{(h)}=-r^{(h)},h=1,\dots ,M$$

$$(R^{-1}{)}_h=\left[\begin{array}{cc}{E}_h& O_{h\times (M-h)}\\ O_{(M-h)\times h}& O_{(M-h)\times (M-h)}\end{array}\right]\left[\begin{array}{c}{u}^{(h)}\\ O_{(M-h)\times h}\end{array}\right]$$

$$\displaystyle h_{j}(x)=h\bigg{(}\frac{x-u_{j}}{\delta_{N}}\bigg{)},\quad\psi_{j}(x)=\int_{0}^{x}h_{j}(t)\,dt,\quad\phi_{j}(x)=\int_{0}^{x}\int_{0}^{t}h_{j}(u)\,du\,dt;\quad x\in[0,1]$$

$$f(x)-f(0)={\int }_0^x{f}^{\prime }(t)dt.$$

$$f(x)={\xi }_0+\sum _{j=0}^N{\xi }_{j+1}{\psi }_j(x).$$

$$f(x)-f(0)-x{f}^{\prime }(0)={\int }_0^x{\int }_0^t{f}^{\prime \prime }(s)dsdt.$$

$$f(x)={\xi }_0+{\xi }^{\ast }x+\sum _{j=0}^N{\xi }_{j+1}{\varphi }_j(x).$$

$$Y={\xi }_0{1}_n+\Psi \xi +\epsilon ,\epsilon \sim {\mathcal{N}}_n(0,{\sigma }^2{\mathrm{I}}_n),\xi \in {\mathcal{C}}_{\xi }^{N+1},$$

$$Y={\xi }_0{1}_n+{\xi }^{\ast }X+\Phi \xi +\epsilon ,\epsilon \sim {\mathcal{N}}_n(0,{\sigma }^2{\mathrm{I}}_n),\xi \in {\mathcal{C}}_{\xi }^{N+1},$$

$$\pi(Y,\xi_{0},\xi,\sigma^{2},\tau^{2})\propto\big{(}\sigma^{2}\big{)}^{-\frac{n}{2}-1}\exp\bigg{\{}-\frac{1}{2\sigma^{2}}\|Y-\xi_{0}\mathrm{1}_{n}-\Psi\xi\|^{2}\bigg{\}}\\ \,\big{(}\tau^{2}\big{)}^{-\frac{N+1}{2}-1}\,\exp\bigg{\{}-\frac{1}{2\tau^{2}}\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi\bigg{\}}\,\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi)$$

$$\pi(Y,\xi_{0},\xi_{*},\xi,\sigma^{2},\tau^{2})\propto\big{(}\sigma^{2}\big{)}^{-\frac{n}{2}-1}\exp\bigg{\{}-\frac{1}{2\sigma^{2}}\|Y-\xi_{0}\mathrm{1}_{n}-\xi_{*}X-\Phi\xi\|^{2}\bigg{\}}\big{(}\tau^{2}\big{)}^{-\frac{N+1}{2}-1}\exp\bigg{\{}-\frac{1}{2\tau^{2}}\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi\bigg{\}}\,\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi)$$

Full text

Keywords: Circulant embedding; Durbin’s recursion; Elliptical Slice Sampling; Smooth relaxation; Toeplitz

1 Introduction

In diverse application areas, it is often of interest to estimate a function nonparametrically subject only to certain constraints on its shape. Typical examples include (but are not limited to) monotone dose-response curves in medicine (Kelly and Rice, 1990), concave utility functions in econometrics (Meyer and Pratt, 1968), increasing growth curves, non-increasing survival function or ‘U’-shaped hazard function (Reboul, 2005) in survival analysis, computed tomography (Prince and Willsky, 1990), target reconstruction (Lele et al., 1992), image analysis (Goldenshluger and Zeevi, 2006), queueing theory (Chen and Yao, 1993), and circuit design (Nicosia et al., 2008).

A Bayesian framework offers a unified probabilistic way of incorporating various shape constraints and accordingly there is a large literature devoted to Bayesian shape constrained estimation. A general approach is to expand the unknown function in a basis and translating the functional constraints to linear constraints in the coefficient space. Some representative examples include piecewise linear models (Neelon and Dunson, 2004; Cai and Dunson, 2007), Bernstein polynomials (Curtis and Ghosh, 2009), regression splines (Meyer et al., 2011), penalized spines (Brezger and Steiner, 2008), cumulative distribution functions (Bornkamp and Ickstadt, 2009), and restricted splines (Shively et al., 2011) used as the basis. Maatouk and Bay (2017) recently exploited a novel basis representation to equivalently represent various shape restrictions such as boundedness, monotonicity, convexity etc as non-negativity constraints on the basis coefficients. Although originally developed in the context of computer model emulation, the approach of Maatouk and Bay (2017) is broadly applicable to general shape constrained problems. Zhou et al. (2018) adapted their approach to handle a combination of shape constraints in a nuclear physics application to model the electric form factor of a proton. The main idea of Maatouk and Bay (2017) is to expand a Gaussian process using a first or second order exact Taylor expansion, with the remainder term approximated using linear combinations of compactly supported triangular basis functions. A key observation is that the resulting approximation has the unique advantage of enforcing linear inequality constraints on the function space through an equivalent linear constraint on the basis coefficients. In terms of model fitting under a standard Gibbs sampling framework, this necessitates sampling from a high-dimensional truncated multivariate normal (tMVN) distribution.

The problem of sampling from a tMVN distribution is notoriously challenging in high dimensions and a number of solutions have been proposed in the literature. Existing Gibbs samplers for a tMVN distribution sample the coordinates one-at-a-time from their respective full conditional truncated univariate normal distributions (Geweke, 1991; Kotecha and Djuric, 1999; Damien and Walker, 2001; Rodriguez-Yam et al., 2004). While the Gibbs sampling procedure is entirely automated, such one-at-a-time updates can lead to slow mixing, especially if the variables are highly correlated. More recently, Pakman and Paninski (2014) proposed a Hamiltonian Monte Carlo (HMC) algorithm which has drastically improved the speed and efficiency of sampling from tMVNs. However, implementing this algorithm within a larger Gibbs sampler can still lead to inefficiencies if the sample size is large. The second and third authors of this article encountered this challenge in Zhou et al. (2018) with a sample size greater than 1000. A related issue which contributes to the complexity is the $\mathcal{O}(N^{3})$ computation and $\mathcal{O}(N^{2})$ storage requirements for inverting and storing a general $N\times N$ covariance matrix.

In this article, we propose a novel algorithm to exploit additional structure present in the tMVN distributions arising in the aforesaid shape-constrained problems using the basis of Maatouk and Bay (2017). Our approach is based on a novel combination of elliptical slice sampling (ESS; Murray et al. (2010)), circulant embedding techniques, and smooth relaxations of hard constraints. We additionally use Durbin’s recursion to efficiently update hyperparameters within the covariance kernel of the parent Gaussian process. We analyze the per-iteration complexity of the proposed algorithm with and without hyperparameter updates and illustrate through simulated and real data examples that the proposed algorithm provides significant computational advantages while retaining the statistical accuracy. R code to implement the proposed algorithm for monotone and convex function estimation is provided at https://github.com/raypallavi/BNP-Computations. We note that our algorithm and code be trivially adapted to other basis functions.

The rest of the paper is organized as follows. In §2, we revisit the Bayesian shape constrained function estimation problem and describe the novel algorithm for inference. The algorithm is specialized to estimating monotone and convex functions §3 with illustrations in §4. We conclude with a real data application in §5. Various algorithmic and implementation details are provided in an Appendix.

2 Algorithm development

Consider the problem of sampling from a distribution having the following form:

[TABLE]

where, $Z\in\mathbb{R}^{n}$, $\Phi\in\mathbb{R}^{n\times N}$ with $N\leq n$, $\mathcal{C}_{\xi}\subset\mathbb{R}^{N}$ is determined by a set of linear inequality constraints on $\xi$, and $K$ is positive definite matrix. While our methodology generally applies to any such $K$, we are specifically interested in situations where $K$ arises from the evaluation of a stationary covariance kernel on a regular grid.

The distribution (1) arises as a conditional posterior of basis coefficients in many Bayesian nonparametric regression problems where linear shape constraints (such as monotonicity, convexity, or a combination of these; Zhou et al. (2018)) on the regression function are present, and a constrained Gaussian process prior is placed on the coefficient vector $\xi$ in an appropriate basis representation. Sampling from the density (1) is then necessitated within a larger Gibbs sampler to fit the said constrained regression model.

Specifically, suppose we observe response-covariate pairs $(y_{i},x_{i})\in\mathbb{R}\otimes\mathbb{R}^{d}$ for $i=1,\ldots,n$, related by the Gaussian regression model

[TABLE]

where the unknown regression function $f$ is constrained to lie in some space $\mathcal{C}_{f}$, a subset of the space of all continuous functions on $[0,1]^{d}$. When $\mathcal{C}_{f}$ corresponds to the space of monotone or convex functions, Maatouk and Bay (2017) identified a novel basis representation for $f$ which allowed equivalent representations of the aforesaid constraints in terms of the basis coefficients $\xi$ restricted to the positive orthant:

[TABLE]

We provide more details on their basis representation in §3.1 and §3.2. They also considered the case where the regression function is globally bounded between two constants. See also Zhou et al. (2018) where a combination of interpolation, monotonicity, and convexity constraints can be equivalently expressed in terms of linear constraints on the coefficients.

Relating the basis coefficients $\xi$ with the function values and (or) its derivatives, Maatouk and Bay (2017) proposed a constrained Gaussian prior on $\xi$. If the function $f$ was unconstrained, then a Gaussian process (GP) prior on $f$ induces a Gaussian prior on $\xi$, aided by the fact that derivatives of GP are again GPs provided the covariance kernel is sufficiently smooth. A natural idea then is to restrict the induced prior on $\xi$ to the constrained region $\mathcal{C}_{\xi}$,

[TABLE]

which is precisely the specification of Maatouk and Bay (2017). The density in equation (1) is then recognized as the conditional posterior of $\xi$.

In what follows, we shall additionally assume that $K_{jj^{\prime}}=k(u_{j}-u_{j^{\prime}})$ for a positive definite function $k$ and a set of uniform grid points $\{u_{j}\}$. For example, if $d=1$, we have $u_{j}=j/N$ for $j=0,1,\ldots,N$. This is a slight departure from Maatouk and Bay (2017) in the monotone and convex case. Since the derivatives of a stationary GP is generally non-stationary, so is their induced prior on $\xi$ from a parent stationary GP on $f$. We instead directly place a stationary GP on an appropriate derivative of $f$, which results in $K$ having a form as above. While there is little difference between the two approaches operationally, there is a large computational benefit for our approach, as we shall see below.

Returning to (1), a simple calculation yields that $p(\xi)$ is a truncated normal distribution, specifically,

[TABLE]

truncated to $\mathcal{C}_{\xi}$. While one can use off-the-shelf samplers for tMVNs (Pakman and Paninski, 2014) to sample from the above, the intrinsic complexity of sampling from tMVNs coupled with the computation and storage of the inverse of the kernel matrix $K$ contributes to the challenges of executing this sampling step for large $N$. In particular, $(\Phi^{\mathrm{\scriptscriptstyle T}}\Phi/\sigma^{2}+K^{-1}/\tau^{2})$ keeps changing over each MCMC iteration with new updates of $\sigma$ and $\tau$, which requires an $N\times N$ matrix inversion at each iteration while applying any of the existing algorithms. The usual Sherman–Morrison–Woodbury matrix inversion trick does not render beneficial in this case. Barring issues with matrix inversions, implementation of such algorithm will be expensive in terms of storage. In addition, if there are unknown hyperparameters in the covariance kernel that get updated at each iteration within a larger MCMC algorithm, either the inversion has to take place at each step, or one has to pre-store a collection of $K^{-1}$ on a fine grid for the hyperparameters.

In this article, we present a different approach to sample from the density $p$ in (1) which entirely avoids matrix inversions. Our approach is based on three basic building blocks: (i) approximating the indicator function in $p$ with a smooth approximant, (ii) a novel use of elliptical slice sampling (Murray et al., 2010) to avoid sampling from truncated non-Gaussian distribution, and (iii) using highly efficient samplers based on the fast Fourier transform for stationary GPs on a regular grid (Wood and Chan, 1994). We describe the details below, starting with a brief review of elliptical slice sampling.

The elliptical slice sampler is a general technique for sampling from posterior distributions of the form,

[TABLE]

proportional to the product of a zero-mean multivariate Gaussian prior with a general likelihood function $L(\cdot)$. In this context, Metropolis–Hastings proposals

[TABLE]

for $\rho\in[-1,1]$ are known to possess good empirical (Neal, 1999) and theoretical (Cotter et al., 2013) properties. The elliptical slice sampler presents an adaptive and automated way to tune the step-size parameter which guarantees acceptance at each step. Specifically, a new location on the randomly generated ellipse determined by the current state $\xi$ and the auxiliary draw $\nu_{e}$ is produced according to

[TABLE]

where the angle $\theta$ is uniformly generated from a $[\theta_{\min},\theta_{\max}]$ interval which is shrunk exponentially fast until an acceptable state is reached. The only requirement is the ability to evaluate $L$ at arbitrary points, rendering the approach broadly applicable.

Turning to (1), note however that the elliptical slice sampler isn’t immediately applicable as we have a truncated normal prior. As a simple fix-up, we approximate the indicator function $\mathbbm{1}_{\mathcal{C}_{\xi}}(\cdot)$ in (1) by a suitable smooth function. Specifically, assuming $\mathcal{C}_{\xi}$ has the same structure as in (3), we use sigmoid-like approximations $\mathbbm{1}_{(0,\infty)}(x)\approx(1+e^{-\eta x})^{-1}$ for large $\eta>0$ to obtain a smooth approximation $\mathbb{J}_{\eta}(\cdot)$ to $\mathbbm{1}_{\mathcal{C}_{\xi}}(\cdot)$ as

[TABLE]

The parameter $\eta$ controls the quality of the approximation; higher the value of $\eta$, better is the approximation. Experimenting across a large number of simulation scenarios, we find that $\eta=50$ already provides a highly accurate approximation for all practical purpose.

With $\mathbb{J}_{\eta}(\xi)$ defined like this, let us define

[TABLE]

Our goal now is to sample from $\widetilde{p}$ which we shall consider as approximate samples from $p$.

Now we can apply ESS to draw samples from $\widetilde{\pi}(\xi\mid-)$, since treating the quantity in the square brackets in (5) as “redefined likelihood”, $\xi$ has an (untruncated) multivariate Gaussian prior, which we call as the “working prior”. So one just needs to draw from the “working prior” distribution and compute the logarithm of the “redefined likelihood” function. In our case, computing the log-likelihood function has computational cost of $\mathcal{O}(nN)$ and we are to sample $\nu_{e}\sim\mathcal{N}(0,\tau^{2}\,K)$, which is usually $\mathcal{O}(N^{3})$. Note that these computational complexities correspond to a single iteration of the MCMC sampler.

Now under the assumption that the covariance matrix $K$ is obtained from a regular grid in $[0,1]$, sampling from the “working prior” is same as simulating realizations of a stationary Gaussian Process on a regular grid in $[0,1]$. Such a covariance matrix is known to have a Toeplitz structure and the simulation can be carried out using the sampling scheme developed by Wood and Chan (1994) which reduces the overall complexity of the algorithm to a commendable extend. The details of this algorithm is discussed in the following section.

2.1 Sampling from the prior distribution of $\xi$

Sampling from the “working prior” distribution requires sampling from a stationary GP on a regular grid in $[0,1]$ with a Toeplitz structure of the covariance matrix. In such settings, the algorithm of Wood and Chan (1994) based on clever embedding techniques can be readily applied. In particular, they exploit the discrete fast Fourier transform twice to offer substantially reduced compared cost. We briefly discuss some of the key ingredients of the algorithm.

The goal is to sample a random vector of the form

[TABLE]

from a mean-zero Gaussian random process on each of the grid points

[TABLE]

with covariance function $\gamma:\mathbb{R}\to\mathbb{R}$. Then $Z\sim\mathcal{N}_{m}(0,G)$, where

[TABLE]

It is to be noted that $G$ is a Toeplitz matrix, which is equivalent to $\tau^{2}K$ in our notation.

There are two basic steps of this method:

• 1

Embedding $G$ in a circulant covariance matrix $C$ of order $d\times d,\text{where}\hskip 8.53581ptd=2^{g}$, for some integer $g$ and $d\geq 2(m-1)$. The circulant matrix is formed such a way that by construction, $C$ is symmetric and the $m\times m$ submatrix in the top left corner of $C$ is equal to $G$. For details on the embedding technique, one can refer to Wood and Chan (1994).

• 2

Using fast Fourier transform twice to generate Z=\big{(}Z_{0},Z_{1},\cdots,Z_{d-1}\big{)}\sim\mathcal{N}_{d}(0,C). Then due to appropriate construction of $C$, \big{(}Z_{0},Z_{1},\cdots,Z_{m-1}\big{)}\sim\mathcal{N}_{m}(0,G).

Thus the problem of sampling from $\mathcal{N}_{m}(0,G)$ changes to sampling from $\mathcal{N}_{d}(0,C)$, which is much more efficient due to double implementation of super efficient FFT algorithm in the procedure. This is done with computational complexity of $\mathcal{O}(d\log d)$.

It is to be noted that, the circulant matrix $C$ is not guaranteed to be positive definite. In practice, it is usually possible to choose $d$ large enough such that $C$ is nonnegative definite (Wood and Chan, 1994) and in such cases the procedure is exact in principle. But if such a $d$ cannot be found or if it is too large to be practicable, then an approximate method is used to make $C$ nonnegative definite, where the matrix $C$ is split into nonnegative definite and nonpositive definite parts; For more details on the approximate method one can consult Wood and Chan (1994). For our algorithm, we prefer the exact method and use an appropriate $d$ that satisfies all the necessary conditions. Applying this algorithm a sample is drawn from $\mathcal{N}_{N}(0,\tau^{2}\,K)$ with $\mathcal{O}(N\log N)$ (for each MCMC iteration) computations, instead of $\mathcal{O}(N^{3})$ computations, and is used as the proposal $\nu_{e}$ in ESS given by (4).

2.2 Algorithm

We implemented our algorithm with $K$ as a stationary Mátern kernel with smoothness parameter $\nu>0$ and length-scale parameter $\ell>0$. Our method takes design points $X$, observations $Y$, Mátern kernel parameters $\nu$ and $\ell$, $\eta$ as in (5), dimension of the random coefficients $N$ and number of posterior samples $n_{0}$ as inputs and gives $n_{0}$ many posterior samples of $\xi$ as output.

The computation cost for drawing a random sample from the prior distribution usually dominates. But that is not the case here. Since $n>N$, computational cost for computing the log-likelihood, using ESS scheme, dominates which leads to computational complexity of \mathcal{O}\big{(}nN\big{)}, for each MCMC iteration.

2.3 Updating hyperparameters

As already discussed, updating the hyperparameters present in the covariance matrix $K$ is computationally challenging in the absence of any structure in $K$. Any likelihood-based method for updating $\nu$ and $\ell$ (e.g. Metropolis–Hastings) requires computing $K^{-1}$ which leads to $\mathcal{O}(N^{3})$ computational steps and $\mathcal{O}(N^{2})$ storage. Hence the computational complexity per MCMC iteration of Algorithm 1 is always bounded above by $\mathcal{O}(N^{3})$.

However, substantial speed-up is possible in our case as $K$ is a symmetric positive-definite Toeplitz matrix. We turn to a non-trivial but effective approach of finding $K^{-1}$ utilizing inverse Cholesky factor of $K$ using Durbin’s recursion algorithm (Golub and van Loan, 1996) which has a computational complexity of $\mathcal{O}(N^{2})$. The columns of the inverse Cholesky factor of $K$ is obtained by solving Yule–Walker systems. Durbin recursion is a procedure of recursively finding the solution to a system of equations involving a Toeplitz matrix, in particular, it is applicable to Yule–Walker systems. Given real numbers $r_{0},r_{1},\ldots,r_{M-1}$ with $r_{0}=1$ such that T=\big{(}r_{|i-j|}\big{)}\in\mathbb{R}^{M\times M} is positive definite then Durbin’s algorithm computes $u\in\mathbb{R}^{M}$ as a solution of the Yule–Walker problem:

[TABLE]

For more details on Durbin’s recursion for solving Yule–Walker equation, refer to Golub and van Loan (1996).

Now suppose, we have the Cholesky factor $R$ such that $R^{\mathrm{\scriptscriptstyle T}}R=T$ where $R$ is an upper-triangular matrix and the inverse Cholesky factor is given by $R^{-1}$. Therefore, $TR^{-1}=R^{\mathrm{\scriptscriptstyle T}}$ and noting that $R^{\mathrm{\scriptscriptstyle T}}$ is lower-triangular, it is enough to solve only the upper-triangular part of $R^{-1}$. The first $h$ elements of the $h^{th}$ column of $R^{-1}$ can be found as a solution of

[TABLE]

where $T_{h}$ is the $h\times h$ principal submatrix of $T$, $u^{(h)}$ is $h$-dimensional vector of solutions and $r^{(h)}=(r_{1},\ldots,r_{h})^{\mathrm{\scriptscriptstyle T}}$ and each of these $M$ equations can be solved using Durbin’s algorithm mentioned above. Note that, the $h^{th}$ column of $R^{-1}$ denoted by $(R^{-1})_{h}$ is then given by:

[TABLE]

where $E_{h}$ is an exchange matrix of order $h$ with anti–diagonal elements all ones and other elements are all zeros. This approach requires $\mathcal{O}(M^{2})$ computations to find $R^{-1}$.

We considered continuous uniform priors on compactly supported intervals on both $\nu$ and $\ell$, independently of each other. Updating $\nu\mid\xi,-$ and $\ell\mid\xi,-$ using Metropolis–Hastings requires to compute acceptance ratio which involves computation of $\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi$ and $|K|^{-1/2}$ for proposal and current combinations of $(\nu,\ell)$. Using Durbin’s algorithm, we can find $S$ such that $({S^{-1}})^{\mathrm{\scriptscriptstyle T}}S^{-1}=K$ and then $\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi=(S^{\mathrm{\scriptscriptstyle T}}\xi)^{\mathrm{\scriptscriptstyle T}}(S^{\mathrm{\scriptscriptstyle T}}\xi)=\sum_{j=1}^{N}v_{j}$ where $v=S^{\mathrm{\scriptscriptstyle T}}\xi$ and $|K|^{-1/2}=\prod_{j=1}^{N}S_{jj}$. Evidently, computation for $S$ dominates and Durbin’s algorithm allows us to update the hyperparameters in $\mathcal{O}(N^{2})$ computations for each iteration within an MCMC algorithm. Thus per-iteration computational complexity for Algorithm 1 combined with this hyperparameter update technique remains $\mathcal{O}(nN)$ as before.

Figure 1 shows the posterior density plots of $\nu$ (left panel) and $\ell$ (right panel) and comparison with the uniform prior densities. Posterior samples were drawn using Metropolis–Hastings and the previously mentioned computation scheme. Posterior densities suggest that it is possible to learn the hyperparameters through this technique. Moreover, based on numerous simulation studies that we had conducted, the Metropolis–Hastings sampler for the hyperparameter update attained at least 15% acceptance probability.

3 Application to shape constrained estimation

We now return to the constrained Gaussian regression setup in (2), and consider applications of our sampling algorithm to situations when $f$ is a smooth monotone or convex function. We first introduce some notation to define the basis functions employed by Maatouk and Bay (2017).

Let $\{u_{j}\in[0,1],\,j=0,1,\dots,N\}$ denote equally spaced knots on $[0,1]$ with spacing $\delta_{N}=1/N$ and $u_{j}=j/N$. Let

[TABLE]

where $h(x)=(1-|x|)\,\mathbbm{1}_{[-1,1]}(x)$. The collection of functions $\{h_{j}\}$ is called the interpolation basis by Maatouk and Bay (2017), since for any continuous function $f:[0,1]\to\mathbb{R}$, the function $\widetilde{f}(\cdot)=\sum_{j=0}^{N}f(u_{j})\,h_{j}(\cdot)$ approximates $f$ by linearly interpolating between the function values at the knots $\{u_{j}\}$.

The integrated basis $\{\psi_{j}\}$ and $\{\phi_{j}\}$ take advantage of higher-order smoothness. For example, if $f$ is continuously differentiable, then by the fundamental theorem of calculus,

[TABLE]

Expanding $f^{\prime}$ in the interpolation basis implies the model

[TABLE]

Similarly, if $f$ is twice continuously differentiable, we have

[TABLE]

Now expanding $f^{\prime}$ and $f^{\prime\prime}$ in the interpolation basis implies the model

[TABLE]

Maatouk and Bay (2017) showed that under (6), $f$ is monotone non-decreasing if and only if $\xi_{i}\geq 0$ for all $i=1,\ldots,N+1$. Similarly, under (7), $f$ is convex non-decreasing if and only if $\xi_{i}\geq 0$ for all $i=1,\ldots,N+1$. This equivalence relationship between the functional constraint and the linear inequality constraints on the basis coefficients is an attractive feature of the interpolation basis and isn’t shared by many commonly used basis functions.

For an unrestricted $f$, a GP prior on $f$ implies a dependent Gaussian prior on the coefficient vector $\xi=(\xi_{1},\ldots,\xi_{N+1})^{{\mathrm{\scriptscriptstyle T}}}$. A natural idea is to restrict this dependent prior subject to the linear restrictions on the coefficients which results in a dependent tMVN prior. Fitting the resulting model using a Gibbs sampler, the full conditional of $\xi$ assumes the form (1), rendering our Algorithm 1 applicable.

We provide more details regarding the model and prior for the monotone and convex cases separately. Let $X=(x_{1},\ldots,x_{n})^{\mathrm{\scriptscriptstyle T}}$ be the vector of $n$ design points, $Y=(y_{1},\ldots,y_{n})^{\mathrm{\scriptscriptstyle T}}$ be the vector of corresponding responses.

3.1 Monotonicity constraint

We can express (6) in vector notation as

[TABLE]

where recall from (3) that $\mathcal{C}_{\xi}^{m}$ denotes the positive orthant in $\mathbb{R}^{m}$ and $\xi=(\xi_{1},\ldots,\xi_{N+1})^{{\mathrm{\scriptscriptstyle T}}}$. Also, $\Psi$ is an $n\times(N+1)$ basis matrix with $i$th row $\Psi_{i}^{\mathrm{\scriptscriptstyle T}}$ where $\Psi_{i}=(\psi_{0}(x_{i}),\ldots,\psi_{N}(x_{i}))^{\mathrm{\scriptscriptstyle T}}$ and $\mathrm{1}_{n}$ denotes an $n$ dimensional vector of all $1$’s.

The parameter $\xi_{0}\in\mathbb{R}$ is unrestricted, and we place a flat prior $\pi(\xi_{0})\propto 1$ on $\xi_{0}$. We place a tMVN prior on $\xi$ independently of $\xi_{0}$ as $p(\xi)\propto\mathcal{N}(\xi;0,\tau^{2}\,K)\,\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi)$, where $K=(K_{jj^{\prime}})$ with $K_{jj^{\prime}}=k(u_{j}-u_{j}^{\prime})$ and $k(\cdot)$ the stationary Mátern kernel with smoothness parameter $\nu>0$ and length-scale parameter $\ell>0$. To complete the prior specification, we place improper priors $\pi(\sigma^{2})\propto 1/\sigma^{2}\,;\,\,\pi(\tau^{2})\propto 1/\tau^{2}$ on $\sigma^{2}$ and $\tau^{2}$, and compactly supported priors $\nu\sim\mathcal{U}(0.5,1)$ and $\ell\sim\mathcal{U}(0.1,1)$ on $\nu$ and $\ell$. A straightforward Gibbs sampler is used to sample from the joint posterior of $(\xi_{0},\xi,\sigma^{2},\tau^{2},\nu,\ell)$ whose details are deferred to the Appendix. The parameters $\sigma^{2}$, $\tau^{2}$ and $\xi_{0}$ have standard conditionally conjugate updates. The key feature of our algorithm is sampling the high-dimensional parameter $\xi$ using Algorithm 1 and updating $\nu$ and $\ell$ via Metropolis-within-Gibbs using Durbin’s recursion as outlined in §2.3.

Before concluding this section, we comment on a minor difference in our prior specification from Maatouk and Bay (2017), which nevertheless has important computational implications. Since the basis coefficients $\xi_{j},j\geq 1$ target the derivatives $f^{\prime}(u_{j})$, Maatouk and Bay (2017) consider a joint prior on $(\xi_{0},\xi)$ obtained by computing the induced prior on \big{(}f(0),f^{\prime}(u_{0}),\ldots,f^{\prime}(u_{N})\big{)} from a GP prior on $f$, and then imposing the non-negativity restrictions. Since the derivative of a sufficiently smooth GP is again a GP, the joint distribution of \big{(}f(0),f^{\prime}(u_{0}),\ldots,f^{\prime}(u_{N})\big{)} can be analytically calculated. However, one downside is that the derivative of a stationary GP is no longer stationary in general, and thus sampling from the joint Gaussian prior of \big{(}f(0),f^{\prime}(u_{0}),\ldots,f^{\prime}(u_{N})\big{)} cannot take advantage of the embedding techniques for a stationary GP. We instead directly place a prior on $\xi$ induced from a stationary GP prior on $f^{\prime}$ and then imposing the necessary restrictions. Since $\xi_{0}$ is only a single real-valued parameter, we break the dependence between $\xi_{0}$ and $\xi$ in the prior and assign a flat prior on $\xi_{0}$ independent of $\xi$. Our simulations suggest against any loss of efficiency in doing so, while there is a substantial computational gain because Algorithm 1 becomes readily applicable to update $\xi$ with our prior.

3.2 Convexity constraint

The development here proceeds in a similar fashion to the monotone case and we only provide a brief sketch. We can write (7) in vector notation as

[TABLE]

where $\Phi$ is an $n\times(N+1)$ basis matrix with $i$th row $\Phi_{i}^{\mathrm{\scriptscriptstyle T}}$ and $\Phi_{i}=(\phi_{0}(x_{i}),\ldots,\phi_{N}(x_{i}))^{\mathrm{\scriptscriptstyle T}}$. The only additional parameter here from the previous case is $\xi^{\ast}$ to which we assign a flat prior on $\xi^{\ast}$ independent of everything else. For all other parameters, the exact same prior specification is employed as in the previous case. The details of the Gibbs sampler are once again deferred to the Appendix. Once again, the salient feature is the update of $\xi$ using Algorithm 1 and the usage of Durbin’s recursion to update $\nu$ and $\ell$.

4 Simulation Results

We now conduct a simulation study to provide empirical evidence that the proposed algorithm provides substantial speed-ups without sacrificing accuracy. We consider two examples corresponding to a monotone and convex truth respectively. For the monotone case, the true function $f(x)=\log(20x+1)$, also considered in Maatouk and Bay (2017), while in the convex case, the true $f(x)=5(x-0.5)^{2}$. In both cases, we uniformly generated the covariates on $[0,1]$ and added Gaussian noise. For all our simulations, we set $\eta=50$ and the number of knots to be half the sample-size, $N=\lceil n/2\rceil$. With $N=\lceil n/2\rceil$, the computational complexity of Algorithm 1 within a single iteration of MCMC sampler is \mathcal{O}\big{(}n^{2}\big{)}. For the hyperparameters $\nu$ and $\ell$ inside the covariance kernel, we considered either of two situations where these were fixed a priori or updated inside the MCMC algorithm. When fixed, we used default choices of $\nu=0.75$ and $\ell$ chosen so that the correlation at a maximum possible separation between the covariates equals $0.05$. When updated, we employed independent priors $\nu\sim\mathcal{U}(0.5,1)$ and $\ell\sim\mathcal{U}(0.1,1)$, as described earlier. Figure 2 provides an illustration of the point estimation along with uncertainty characterization using our Gibbs sampler for the monotone (left panel) and convex (right panel) cases.

4.1 Cost per-iteration

We first illustrate the gains in per-iteration cost due to Algorithm 1 as opposed to directly sampling from a tMVN. We consider an alternative Gibbs sampler which samples $\xi\in\mathcal{C}_{\xi}^{N+1}$ from its tMVN full-conditional using the Hamiltonian Monte Carlo (HMC) sampler of Pakman and Paninski (2014), implemented in the ** R** package “tmg”. We kept $\nu$ and $\ell$ fixed and all other parameters are updated in the exact same way in either sampler. We did not consider the rejection sampler used by Maatouk and Bay (2017) as it becomes quite inefficient with increasing dimension, with the “tmg” sampler substantially more efficient than the rejection sampler in high dimensions. The combination of Maatouk and Bay (2017) with the “tmg” sampler does not exist in the literature to best of our knowledge, and thus we are being entirely fair in constructing the best possible competing method.

Figure 3 plots the run-time per iteration (in seconds) against the sample size $n$ (varied between $50$ to $1000$ ) for the two approaches, both implemented under identical conditions on a quadcore Intel Core i7-2600 computer with 16 GB RAM. Evidently, Algorithm  1 provides more pronounced improvements for larger $N$. While not shown here, one obtains a similar picture when the hyperparameters $\nu$ and $\ell$ are updated. Naturally, both methods incur additional computation time due to the additional calculations to compute the Metropolis–Hastings ratio. The “tmg” sampler incurs additional burden in this case due to the continuous nature of the priors on $\nu$ and $\ell$, which either require storage of matrix inverses on a fine grid, or new matrix multiplications on the fly at each iteration. For fair comparison, while implementing “tmg” along with hyperparameter updates, matrix inversions within each MCMC iterations were carried out using matrix inversion function of “FastGP” that utilizes Trench’s algorithm and requires $\mathcal{O}(N^{2})$.

4.2 Effective sample sizes

We now investigate the mixing behavior of our Gibbs sampler by computing the effective sample size of the estimated function value at $75$ different test points. The effective sample size is a measure of the amount of the autocorrelation in a Markov chain, and essentially amounts to the number of independent samples in the MCMC path. From an algorithmic robustness perspective, it is desirable that the effective sample sizes remain stable across increasing sample size and/or dimension, and this is the aspect we wish to investigate here. We only report results for the monotonicity constraint; similar behavior is seen for the convexity constraint as well.

We consider 20 different values for the sample size $n$ with equal spacing between $50$ and $1000$. Note that the dimension of $\xi$ itself grows between $25$ and $500$ as a result. For each value of $n$, we run the Gibbs sampler for 12,000 iterations with $5$ randomly chosen initializations. For each starting point, we record the effective sample size at each of the $75$ test points after discarding the first 2,000 iterations as burn-in, and average them over the different initializations. Figure 4 shows boxplots of these averaged effective sample sizes across $n$; the left and right panel corresponds to $\nu$ and $\ell$ being fixed or updated within the MCMC, respectively. In either case, the effective sample sizes remain quite stable across growing $n$. When $\nu$ and $\ell$ are updated, there is a small dip in the median effective sample size which is pretty much expected.

5 Application on a real data set

We considered a car price data111available at https://bitbucket.org/remcc/mbart/src/30f684da12ce5ad2d28a9b3c0a5fb166a5e15697/data/?at=master studied by Chipman et al. (2016) on the sales prices of used Mercedes cars. The data consists of $1000$ observations on the car prices and we took mileage to be the explanatory variable. The relation between price and mileage is monotone decreasing, since a car with higher mileage is expected to have a lower price.

We scaled the mileage to lie within $[0,1]$ and considered the log-transformed sales price as the response. As with the simulation study, we compared the performance of our algorithm to the exact algorithm which samples $\xi$ from its tMVN full conditional using the “tmg” package. To compare the predictive performance, we considered 10 different random splits of the data into training and test tests. On each split, we reserved 900 training points for model fitting and the remaining 100 data points were held out for predictive evaluation. For each of the 10 datasets, both samplers were run for $6000$ iterations with the first $1000$ discarded as burn-in in each case. We also took into account the cases of fixing $\nu$ and $\ell$ or updating them within the MCMC. The average MSPE across the 10 splits for our method is $0.16645$ with standard error $0.0333$, while for the exact sampler the average MSPE is $0.16646$ with standard error $0.0332$ corresponding to the case where $\nu$ and $\ell$ were updated within the MCMC. As a visual illustration, Figure 5 shows the performance of both methods with hyperparameter updates, on the last fold of the $10$-fold validation process. While both algorithms resulted in similar accuracy, it took $0.25$ seconds to run one iteration of our algorithm in comparison to $0.39$ seconds for the exact sampler, corresponding to the case where $\nu$ and $\ell$ were updated in each iteration of the MCMC. For fixed $\nu$ and $\ell$, it took $0.014$ seconds to run one iteration of our algorithm in comparison to $0.055$ seconds for the exact sampler. Thus, as anticipated, our algorithm continues to multi-fold speed-ups without sacrificing on accuracy.

Appendix A Appendices

A.1 Full conditionals

Consider model (8) and the prior specified in section 3.1. The joint distribution is given by:

[TABLE]

Then,

$\xi\mid Y,\xi_{0},\sigma^{2},\tau^{2}$ is truncated multivariate Gaussian truncated on $\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi)$.

$\xi_{0}\mid Y,\xi,\sigma^{2},\tau^{2}\sim\mathcal{N}(\bar{Y}^{*},\sigma^{2}/n)$, where, $\bar{Y}^{*}$ is average of components of $Y^{*}=Y-\Psi\xi$.

\sigma^{2}\mid Y,\xi_{0},\xi,\tau^{2}\sim\mathcal{I}\mathcal{G}\big{(}n/2,\|Y-\xi_{0}\mathrm{1}_{n}-\Psi\xi\|^{2}/2\big{)}

\tau^{2}\mid Y,\xi_{0},\xi,\sigma^{2}\sim\mathcal{I}\mathcal{G}\big{(}(N+1)/2,\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi/2\big{)}

Again, consider model (9) and the prior specified in section 3.2. The joint distribution is given by:

[TABLE]

Then,

$\xi\mid Y,\xi_{0},\xi_{*},\sigma^{2},\tau^{2}$ is truncated multivariate Gaussian truncated on $\mathbbm{1}_{\mathcal{C}_{\xi}}(\xi)$.

$\xi_{0}\mid Y,\xi_{*},\xi,\sigma^{2},\tau^{2}\sim\mathcal{N}(\bar{Y}^{*},\sigma^{2}/n)$, $\bar{Y}^{*}$ is average of components of $Y^{*}=Y-\xi_{*}X-\Phi\xi$.

$\xi_{*}\mid Y,\xi_{0},\xi,\sigma^{2},\tau^{2}\sim\mathcal{N}(\sum_{i=1}^{n}x_{i}y_{i}^{**}/\sum_{i=1}^{n}x_{i}^{2},\sigma^{2}/\sum_{i=1}^{n}x_{i}^{2})$, where $Y^{**}=Y-\xi_{0}\mathrm{1}_{n}-\Phi\xi$.

\sigma^{2}\mid Y,\xi_{0},\xi_{*},\xi,\tau^{2}\sim\mathcal{I}\mathcal{G}\big{(}n/2,\|Y-\xi_{0}\mathrm{1}_{n}-\xi_{0}X-\Phi\xi\|^{2}/2\big{)}

\tau^{2}\mid Y,\xi_{0},\xi_{*},\xi,\sigma^{2}\sim\mathcal{I}\mathcal{G}\big{(}(N+1)/2,\xi^{\mathrm{\scriptscriptstyle T}}K^{-1}\xi/2\big{)}

Algorithm 1 was used to draw samples from the full conditional distribution of $\xi$ while sampling from the full conditionals of $\xi_{0}$, $\xi_{*}$, $\sigma^{2}$ and $\tau^{2}$ are routine.

A.2 R code

We used ** R** for the implementation of Algorithm 1 and Durbin’s recursion to find the inverse of the Cholesky factor, with the computation of the inverse Cholesky factor optimized with ** Rcpp**. We provide our code for implementing the monotone and convex function estimation procedures in §4 and §5 in the Github page mentioned in §1. There are six different functions to perform the MCMC sampling for monotone increasing, monotone decreasing, and convex increasing functions with and without hyperparameter updates. Each of these main functions take $x$ and $y$ as inputs along with other available options, and return posterior samples on $\xi_{0}$, $\xi^{*}$, $\xi$, $\sigma$, $\tau$ and $f$ along with posterior mean and 95% credible interval of $f$ on a user-specified grid. A detailed description on the available input and output options for each function can be found within the function files.