# A Privacy Preserving Method with a Random Orthogonal Matrix for ConvMixer
Models.

## Abstract

In this paper, a privacy preserving image classification method is proposed
under the use of ConvMixer models. To protect the visual information of test
images, a test image is divided into blocks, and then every block is encrypted
by using a random orthogonal matrix. Moreover, a ConvMixer model trained with
plain images is transformed by the random orthogonal matrix used for
encrypting test images, on the basis of the embedding structure of ConvMixer.
The proposed method allows us not only to use the same classification accuracy
as that of ConvMixer models without considering privacy protection but to also
enhance robustness against various attacks compared to conventional privacy-
preserving learning.

##  1\. Introduction

Deep learning has been deployed in many applications including security-
critical ones. Generally, data contains sensitive information such as personal
informational, so privacy-preserving methods for deep learning have become an
urgent problem [1]. To achieve privacy-preserving learning, various methods
have been proposed. One of them is Federated Learning (FL) [2], which is a
type of distributed learning. FL allows us to train a model over multiple
participants without directly sharing their raw data. However, FL have not
considered the protection of test data in cloud environments so far. In this
paper, we propose a novel method for protecting visual information on test
images.

To protect visual information on plain images in untrusted cloud environments,
many learnable encryption methods have been studied so far [3]-[13]. Learnable
encryption has to satisfy three requirements in general: (a) having a high
accuracy that is almost the same as that of plain models, (b) being robust
enough against various attacks, and (c) easily updating a secret key. However,
most of existing methods [3]-[11] degrade the accuracy of models due to the
use of encrypted images, and moreover, need to retrain models to update the
key. In contrast, the similarity between block-wise encryption and the
architecture of isotropic networks has been pointed out to enable us to
perfectly stratify the two requirements that the existing methods cannot
[12][13]. Information on embeddings in isotropic networks such as the vision
transformer [14] and ConvMixer [15] is encrypted by random matrixes generated
with secret keys for privacy-preserving learning. However, in the conventional
methods [12][13], simple permutation matrixes are used for image and model
encryption, so encrypted images are not robust enough against various attacks.
Accordingly, we propose the use of a novel random matrix, which is called a
random orthogonal one generated by using the Gram-Schmidt orthonormalization.
The proposed method allows us to enhance the visual protection of images,
while maintaining the same as that of plain models and the easy update of a
secret key.

##  2\. ConvMixer

Before discussing the proposed method, we summarize ConvMixer and its
properties briefly. ConvMixer is mainly used for image classification tasks
and is known for its high classification performance[15]. The structure of
ConvMixer is inspired by the Vision Transformer (ViT)[14]. ViT consists of two
Embedding processes (Patch Embedding and Position Embedding) and a Transformer
structure. On the other hand, ConvMixer consists of a Patch Embedding and a
CNN structure. Figure 2 shows the structure of ConvMixer, which consists of
two main structures: Patch Embedding and ConvMixer Layer. In this paper, we
focus on Patch Embedding. In Patch Embedding, an input image
$x\in\mathbb{R}^{H\times W\times C}$ of height $H$, width $W$, and number of
channels $C$ is divided into patches of size $p\times p$. Each patch is then
transformed into a vector $x_{p}^{i}\in\mathbb{R}^{p^{2}C}$, multiplied by the
learnable filter $E$ and linearly transform it into a vector of $d$-dimensions
by taking the product of $x_{p}^{i}\in\mathbb{R}^{p^{2}C}$ as

$\displaystyle z=[x_{p}^{1}E,...,x_{p}^{i}E,...,x_{p}^{N}E]$ (1)
$\displaystyle z\in\mathbb{R}^{N\times d},E\in\mathbb{R}^{(p^{2}c)\times d}$

In previous studies[12][13], it is known that it is possible to protect the
privacy of test images by transforming the filter $E$ with a secret key. In
this paper, we propose a method to achieve stronger privacy preserving of test
images by using random orthogonal matrices.

Figure 1: Architecture of ConvMixer

##  3\. Proposed Method

###  3.1 Overview

Figure 2: Framework of proposed method

Figure 2 illustrates the framework of the proposed method. The proposed method
aims to protect visual information on test images. To achieve this aim, we
encrypt test images and a transform model by using an random orthogonal
matrix. The framework is summarized as below.

>   * •
>
> A third party (trusted) generates random numbers with a secret key (seed),
> and prepares a random orthogonal matrix $A$ from the random numbers and an
> inverse random orthogonal matrix $A^{-1}$.
>
>   * •
>
> The third party trains a ConvMixer model $\varphi$ with plane images. The
> trained model $\varphi$ is transformed into an encrypted model $\varphi_{T}$
> by using $A^{-1}$.
>
>   * •
>
> The third party provides the random orthogonal matrix $A$ to a client
> (trusted) and model $\varphi_{T}$ to a provider (untrusted).
>
>   * •
>
> The client transforms a test image $x$ into an encrypted image $\hat{x}$ by
> using $A$. After that, the client sends $\hat{x}$ to the provider.
>
>   * •
>
> The provider inputs $\hat{x}$ into model $\varphi_{T}$, and sends back a
> prediction result to the client.
>
>

Even if the provider is not trusted, the client does not give visual
information of test images and matrix $A$ used for image encryption to the
provider. Thus, the client can receive prediction results while maintaining
the privacy preserving of test images.

###  3.2 Test Image Encryption

A test image $x\in\mathbb{R}^{H\times W\times C}$ is transformed into an
encrypted image $\hat{x}\in\mathbb{R}^{H\times W\times C}$ as below.

  1. 1.

Divide $x$ into $N$ blocks with a size of $p\times p$ such that
$B=\\{B_{1},...,B_{N}\\}$, where $p\times p$ is the same size as the patch
size used in a ConvMixer model, and $N$ is $(H\times W)/p^{2}$.

  2. 2.

Flatten each block $B_{i}\in\mathbb{R}^{p\times p\times C}$ into a vector
$x_{p}^{i}\in\mathbb{R}^{p^{2}C}$ as

$\displaystyle x_{p}^{i}=[x_{p}^{i}(1),...,x_{p}^{i}(p^{2}C)].$ (2)

  3. 3.

Generate a encrypted vector $\hat{x}_{p}^{i}\in\mathbb{R}^{p^{2}C}$ by
multiplying vector $x_{p}^{i}$ by matrix
$A\in\mathbb{R}^{(p^{2}C)\times{(p^{2}C)}}$ as

$\displaystyle\hat{x}_{p}^{i}=x_{p}^{i}A.$ (3)

  4. 4.

Rebuild vector $\hat{x}_{p}^{i}$ into block $\hat{B}_{i}$ in the reverse order
of step 2.

  5. 5.

Concatenate $\hat{B}=\\{\hat{B}_{1},...,\hat{B}_{N}\\}$ into an encrypted test
image $\hat{x}$.

###  3.3 Model Encryption

To avoid the performance degradation caused by encryption of test images, $E$
in Eq.(1) is transformed by using $A^{-1}$ as

$\displaystyle E^{\prime}=A^{-1}E.$ (4)

When replacing $E$ and $x_{p}^{i}$ with $E^{\prime}$ and $\hat{x}_{p}^{i}$ ,
respectively, vector z in Eq.(1) is reduced to as

$\displaystyle
z^{\prime}=[\hat{x}_{p}^{1}E^{\prime},...,\hat{x}_{p}^{i}E^{\prime},...,\hat{x}_{p}^{N}E^{\prime}].$
(5)

Thus, by substituting Eqs. (3) and (4) to Eq.(5), we obtain:

$\displaystyle
z^{\prime}=[x_{p}^{1}AA^{-1}E,...,x_{p}^{i}AA^{-1}E,...,x_{p}^{N}AA^{-1}E]$
$\displaystyle=[x_{p}^{1}E,...,x_{p}^{i}E,...,x_{p}^{N}E]=z.$ (6)

From Eq.(5), encrypted model $\varphi_{T}$ allows us to have the same
performance as that of the model trained with plane images, under the use of
encrypted images.

###  3.4 Generation of Random Orthogonal Matrices

A random orthogonal matrix $A$ can be generated by using the Gram–Schmidt
orthonormalization. The procedure for generating $A$ with a size of $n\times
n$ is given as follows.

  1. 1.

Generate an real matrix $R$ with a size of $n\times n$ by using a random
number generator with a seed.

  2. 2.

Calculate $det(R)$, and proceed to 3 if $det(R)\neq 0$. Otherwise, return to
1.

  3. 3.

Compute a random orthogonal matrix $A$ from $R$ by using the Gram-Schmidt
orthogonalization.

In this framework, any regular matrix can be used as A for image encryption.
Several conventional methods for privacy-preserving image classification use
permutation matrices of pixel values, in which many elements have zero values
in matrices as

$\displaystyle A=\begin{bmatrix}0&1&0\\\ 0&0&1\\\ 1&0&0\\\ \end{bmatrix}.$ (7)

In contrast, the proposed random orthogonal matrices include no the zero
values as elements. The use of such matrices allows us not only to more
strongly protect visual information on plain images but to also enhance
robustness against various attacks, while maintaining the same performance as
that of models trained with plain images. In addition, $A^{-1}$ can easily be
calculated as the transposed matrix of A.

##  4\. Experiment Results

To verify the effectiveness of the proposed method, we ran a number of
experiments on the CIFAR-10 dataset.

###  4.1 Setup

We used the CIFAR-10 dataset, which consists of 60,000 color images (dimension
of 32 × 32 × 3) with 10 classes (6000 images for each class) where 50,000
images are for training and 10,000 for testing. ConvMixer was trained and
tested on the CIFAR-10 dataset. In the model setting, we set the patch size to
4, the number of channels after patch embedding to 256, the kernel size of
depth-wise convolution to 7 and the number of ConvMixer layers to 8. models
were trained for 200 epochs with the Adam optimizer, where the learning rate
was 0.001. We also used a random orthogonal matrix with a size of 48 × 48 for
the encryption of test images and models.

###  4.2 Visual Protection Performance

Figure 3 shows an example of images encrypted with a conventional encryption
method[12][13], in which pixel shuffling and negative-positive transformation
are carried out for image encryption, and an example of images encrypted with
the proposed method, where the images had $H\times W\times C=512\times
512\times 3$ as an image size, and the block sizes used for encryption were
$p=8$ and $p=16$. When using an orthogonal matrix for encryption, transformed
pixel values are real values, so (b) in Fig. 3 were displayed after
normalizing the pixel values to the range of [0.1]. From the figures, the
selection of a larger the block size gave smaller visual information. The use
of random orthogonal matrices was also demonstrated to have a stronger visual
protection performance than that of the conventional method. In addition to
visual protection, encrypted images have to be robust enough against various
attacks, which aim to restore visual information from encrypted images. We
already confirmed that images encrypted with the proposed method are more
robust against attacks including jigsaw puzzle solver attacks [16]. In
particular, unlike ViT, ConvMixer models do have position embedding, so the
position of patches cannot be changed. Therefore, privacy-preserving ConvMixer
needs a stronger encryption method than ViT.

|  (a) p = 8 |  (b) p = 8
---|---|---
|  (c) p = 16 |  (d) p = 16
Plane image | (a) | (b)

Figure 3: Example of encrypted images with (a) conventional method[12][13] and
(b) proposed method

###  4.3 Classification Performance

We evaluated the classification performance of the proposed method as shown in
Table 1, where plain and encrypted indicate plane test images and an encrypted
test images, respectively, and plain model and encrypted model are models
trained with plain images, and models trained with encrypted images. Table 1
shows the classification results for each combination. From the table, the
proposed method (the combined use of Encrypted models and encrypted images)
had the same classification accuracy as that of the baseline without privacy
protection (plain models and plain images). Accordingly, the proposed method
can not only protect the visual information of test image, but also classify
the encrypted image without any degradation of classification accuracy.

Table 1: Classification accuracy (%) model\test image | plane | encypted |  |
---|---|---|---|---
plane model | 90.38 | 13.2 |  |
encypted model | 10.27 | 90.38 |  |

##  5\. Conclusion

In this paper, we proposed a novel method for protecting visual informational
on test images under the use of ConvMixer models. The proposed method allows
us to use a random orthogonal matrix for image encryption, and it was
demonstrated not only to enhance the visual protection of images but to also
maintain the same accuracy as that of models trained with plain images.

## Ackowledgment

This study was partially supported by JSPS KAKENHI (Grant Number JP21H01327).

## References

  * [1] H. Kiya, MM. AprilPyone, Y. Kinoshita, S. Imaizumi, and S. Shiota, “An Overview of Compressible and Learnable Image Transformation with Secret Key and its Applications,” APSIPA Transactions on Signal and Information Processing, vol.11, no.1, e11, 2022.
  * [2] H.B. McMahan, E. Moore, D. Ramage, S. Hampson, and B.A.y. Arcas, “Communication-Efficient Learning of Deep Networks from Decentralized Data,” Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, PMLR, vol.54, pp.1273-1282, 2017.
  * [3] I. Nakamura, Y. Tonomura, H. Kiya, “Unitary Transform-Based Template Protection and Its Application to l2-norm Minimization Problems,” IEICE Transactions on Information and Systems, vol.E99-D no.1, pp.60-68, 2016.
  * [4] M. Tanaka, ”Learnable Image Encryption,” 2018 IEEE International Conference on Consumer Electronics-Taiwan (ICCE-TW), pp. 1-2, 2018.
  * [5] W. Sirichotedumrong, T. Chuman, S. Imaizumi and H. Kiya, ”Grayscale-Based Block Scrambling Image Encryption for Social Networking Services,” 2018 IEEE International Conference on Multimedia and Expo (ICME), San Diego, CA, USA, 2018, pp. 1-6.
  * [6] W. Sirichotedumrong and H. Kiya, ”A GAN-Based Image Transformation Scheme for Privacy-Preserving Deep Neural Networks,” 2020 28th European Signal Processing Conference (EUSIPCO), Amsterdam, Netherlands, 2021, pp. 745-749.
  * [7] H. Ito, Y. Kinoshita, MM. Aprilpyone and H. Kiya, “Image to perturbation: An image transformation network for generating visually protected images for privacy-preserving deep neural networks”, IEEE Access, vol.9, pp.64629-64638, 2021.
  * [8] MM. AprilPyone, H. Kiya, “Piracy-resistant DNN watermarking by block-wise image transformation with secret key,” Proceedings of the 2021 ACM Workshop on Information Hiding and Multimedia Security, pp.159-164, 2021.
  * [9] MM. AprilPyone and H. Kiya, ”Privacy-Preserving Image Classification Using an Isotropic Network,” in IEEE MultiMedia, vol. 29, no. 2, pp. 23-33, 1 April-June 2022.
  * [10] Z. Qi, A. MaungMaung, Y. Kinoshita and H. Kiya, ”Privacy-Preserving Image Classification Using Vision Transformer,” 2022 30th European Signal Processing Conference (EUSIPCO), 2022, pp. 543-547.
  * [11] N. Senanayake, R. Podschwadt, D. Takabi, V.D. Calhoun and S. M. Plis“Neurocrypt: Machine learning over encrypted distributed neuroimaging data”, Neuroinform 20, pp.91–108, 2022.
  * [12] H. Kiya, T. Nagamori, S. Imaizumi and S. Shiota “Privacy-Preserving Semantic Segmentation Using Vision Transformer”, Journal of Imaging, vol.8, no.9, 2022.
  * [13] H. Kiya, R. Iijima, AP. Maungmaung, Y. Kinoshita, “Image and model transformation with secret key for vision transformer”, IEICE Transactions on Information and Systems, vol.E106.D, Issue.1, pp.2-11, 2023.
  * [14] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby, “An image is worth 16x16 words: Transformers for image recognition at scale,” arXiv:2010.11929, 2020.
  * [15] A. Trockman and J.Z. Kolter, “Patches are all you need?,” arXiv:2201.09792, 2022.
  * [16] T. Chuman and H. Kiya, ”Security Evaluation of Block-based Image Encryption for Vision Transformer against Jigsaw Puzzle Solver Attack,” 2022 IEEE 4th Global Conference on Life Sciences and Technologies (LifeTech), Osaka, Japan, 2022, pp. 448-451.


2022

[1]Carlos Tadeu Pagani Zanini

[1]Departamento de Métodos Estatísticos, Universidade Federal do Rio de
Janeiro, Rio de Janeiro, RJ, Brazil

2]Instituto de Matemática, Estatística e Ciência da Computação, Universidade
Estadual de Campinas, Campinas, SP, Brazil

# Variational Inference for Bayesian Bridge Regression

<EMAIL_ADDRESS>Helio S. Migon<EMAIL_ADDRESS>Ronaldo Dias
<EMAIL_ADDRESS>*  [

###### Abstract

The bridge approach for regularization of coefficients in regression models
uses $\ell_{\alpha}$ norm, with $\alpha\in(0,+\infty)$, to define a
penalization on large values of the regression coefficients. Particular cases
include the lasso and ridge penalizations. In Bayesian models, the
penalization is enforced by a prior distribution on the coefficients. Although
MCMC approaches are available for Bayesian bridge regression, they can be very
slow for large datasets, specially in high dimensions. This paper develops an
implementation of Automatic Differentiation Variational Inference for Bayesian
inference on semi-parametric regression models with bridge penalization. The
non-parametric effects of covariates are modeled by B-splines. The proposed
inference procedure allows the use of small batches of data at each iteration
(due to stochastic gradient based updates), therefore drastically reducing
computational time in comparison with MCMC. Full Bayesian inference is
preserved so joint uncertainty estimates for all model parameters are
available. A simulation study shows the main properties of the proposed method
and an application to a large real dataset is presented.

###### keywords:

Variational Inference, Bridge Penalization, Bayesian Inference, Splines.

##  1 Introduction and related works

It is not uncommon for an experimenter to be interested in understanding how
covariates might explain a response variable. For this, one can assume a
general non-parametric regression model:

$y_{i}=g(x_{i1},\ldots,x_{iD})+\epsilon_{i},$

where $\epsilon_{i}\sim N(0,\phi^{-1})$ for $i=1,\ldots,n$.
A usual approach to estimate the surface $g$ is to consider the well-known
(General) Additive Model, which briefly means that $g$ can be decomposed as

$g(\mathbf{x}_{1},\ldots,\mathbf{x}_{D})=g_{0}+\sum_{j=1}^{D}g_{j}(\mathbf{x}_{j}),$

where $\mathbf{x}_{j}=(x_{1j},\ldots,x_{nj})^{\top}$, for $j=1,\ldots,D$ and
$i=1,\ldots,n$. Each univariate function $g_{j}$ can be written as a linear
combination of basis functions $B_{j,1},\ldots,B_{jK_{j}}$, such as B-splines,
wavelets, Fourier basis, etc. That is,
$g_{j}(\mathbf{x}_{j})=\sum_{k=1}^{K_{j}}B_{jk}(\mathbf{x}_{j})\beta_{jk}$.

Throughout this work, the basis functions chosen to represent the univariate
function $g_{j}$ are the well-known B-splines. Thus, the surface regression
model is:

$y_{i}=g_{0}+\sum^{D}_{j=1}\sum_{k=1}^{K_{j}}B_{jk}(x_{ij})\beta_{jk}+\epsilon_{i}$
(1)

Regularization plays a very important role in Statistics: penalizing
overcomplex models often reduces the risk of overfitting and produces better
generalization to hold out data. For this, under the frequentist point of
view, one have to solve the following optimization problem:

$\operatorname*{arg\,min}_{\bm{\beta}}\
\lVert\mathbf{y}-g_{0}-\sum_{k=1}^{K_{j}}B_{jk}(\mathbf{x}_{j})\beta_{jk}\rVert^{2}+\sum_{j=1}^{D}\lambda_{j}\;{\mathcal{P}}(\bm{\beta}_{j},\alpha_{j}),$

where $\lambda_{j}$ is the smoothing parameter and
${\mathcal{P}}(\bm{\beta}_{j},\alpha_{j})$ is the roughness penalty term.

Different types of penalties, ${\mathcal{P}}(\bm{\beta}_{j},\alpha)$ can be
applied. For instance,
${\mathcal{P}}(\bm{\beta}_{j},\alpha=1)=\lVert\bm{\beta}_{j}\rVert_{1}$,
$\ell_{1}$-norm,
${\mathcal{P}}(\bm{\beta}_{j},\alpha=2)=\lVert\bm{\beta}_{j}\rVert^{2}_{2}$,
the $\ell_{2}$-norm squared, the smoothing splines penalty,
${\mathcal{P}}(\bm{\beta}_{j},\alpha=2)=\int(\sum_{j}\bm{\beta}_{j}\ddot{B}_{j})^{2},$
and P-spline,
${\mathcal{P}}(\bm{\beta}_{j},\alpha=2)=\sum_{j=k+1}^{D}(\Delta^{k}\bm{\beta}_{j})^{2}$,
where $\Delta^{k}$ is the $k$-th order difference operator.

In particular, after writing (1) in matrix form, the semi-parametric
regression model can be written as:

$\mathbf{y}=\mathbf{X}_{0}\bm{\beta}_{0}+\sum_{j=1}^{D}\mathbf{X}_{j}\bm{\beta}_{j}+\bm{\epsilon},$

where $\bm{\epsilon}\sim N(\mathbf{0},\phi^{-1}\mathbf{I}_{n})$, with
$\sum_{k=1}^{K_{j}}B_{jk}(\mathbf{x}_{j})\beta_{jk}=\mathbf{X}_{j}\bm{\beta}_{j}$,
where $\mathbf{X}_{j}$ is the $n\times K_{j}$ matrix
$\mathbf{X}_{j}=[B_{j1}(\mathbf{x}_{j}),\ldots,B_{jK_{j}}(\mathbf{x}_{j})]$
and $\bm{\beta}_{j}=(\beta_{j1},\ldots,\beta_{jK_{j}})^{\top}$. Note that the
matrix $\mathbf{X}_{0}$ comprises possible parametric covariates and hence can
also be viewed as Partial Splines (Silverman Green, 1994). In particular, if
$\mathbf{X}_{0}=\mathbf{1}$ then $\mathbf{X}_{0}\bm{\beta}_{0}$ plays the role
of $g_{0}$ from equation (1) (Additive Model).

Penalties defined in terms of $\ell_{\alpha}$ norm in the context of Bayesian
regression models include the ridge (Hoerl Kennard 19701, Hoerl Kennard 19702)
($\alpha=2$), the lasso (Park Casella, 2008) ($\alpha=1$), elastic net (Q. Li
Lin, 2010) (convex combination of both $\ell_{1}$ and $\ell_{2}$ norms) and
bridge (Polson ., 2014; Mallick Yi, 2018) ($\alpha>0$), first presented by
Frank Friedman (1993). A possible bridge model for semi-parametric regression
is a regularization problem defined as,

$\operatorname*{arg\,min}_{\bm{\beta}_{0},\
\bm{\beta}}\frac{1}{2}\lVert(\mathbf{y}-\mathbf{X}_{0}\bm{\beta}_{0}-\sum_{j=1}^{D}\mathbf{X}_{j}\bm{\beta}_{j})\|^{2}+\sum^{D}_{j=1}\lambda_{j}\lvert\lvert\bm{\beta}_{j}\lvert\lvert^{\alpha_{j}}_{\alpha_{j}}$
(2)

where $\lambda_{j}>0$ is the penalty term that controls the strength of
shrinkage over $\bm{\beta}_{j}$ and $\alpha_{j}>0$ is the concavity parameter
of the penalty function.

The different choices of $\ell_{\alpha}$ norm for the penalization imply
different forms of shrinking the regression coefficients towards zero. While
ridge regression does not zero out coefficients, the lasso penalization is
capable of producing sparse solutions to the corresponding maximum a
posteriori objective, therefore working as a variable selection procedure
(Tibshirani, 1996, 1997). From a Bayesian perspective, Hastie . (2015),
Casella . (2010) and Leng . (2014), highlight that one can interpret the term
$\lVert\bm{\beta}_{j}\lVert^{\alpha_{j}}_{\alpha_{j}}$ as proportional to the
negative log-prior density of $\bm{\beta}_{j}$ with the contours illustrated
in Figure 1 representing the contours of the prior distribution. The case
$\alpha_{j}<1$ implies a non-convex prior that concentrates more mass along
the coordinates’ axis, producing solutions with fewer nonzero coefficients and
less shrinkage.

Figure 1: Constraint regions $\sum^{p_{x}}_{j=1}\mid\beta_{j}\mid^{\alpha}\leq
1$ for different values of $\alpha$. For $\alpha<1$, the constraint region is
nonconvex. Figure extracted from Hastie . (2015) (adapted).

Other forms of penalization are also relevant. The elastic net provides
sparsity of representation similarly to the lasso while also encouraging a
grouping effect where strongly correlated predictors tend to be in or out of
the model together (Zou Hastie, 2005). In Bayesian regression models, the
horseshoe prior implies a form of variable selection that acts strongly on
coefficients of small magnitude while having little influence on coefficients
of large magnitude (Carvalho ., 2009).

In regression models that include non-parametric effects, it is crucial to
assign some form of penalization in order to calibrate the level of complexity
of the relationships between the covariates and the response to avoid
overfitting. Generalized additive models (GAM) from Hastie Tibshirani (1986)
models the non-parametric effects of covariates on the mean of a response
variable belonging to the exponential family, with one spline being used for
each covariate. For common choices of spline basis functions such as B-spline,
see Eilers Marx (1996). The level of smoothness can be controlled by a
penalization acting on the splines coefficients, see Currie Durban (2002).
Some Bayesian approaches in this context include Lang Brezger (2004) and
Hastie Tibshirani (2000). In the context of bridge penalized linear regression
with Gaussian response, Polson . (2014) and Mallick Yi (2018) develop full
Bayesian inference through fixed dimension MCMC based on variable augmentation
schemes that guarantee conjugacy for the full conditional distributions.
Denison . (1998) and Dias Gamerman (2002) describe reversible jumps MCMC
algorithms for selecting the number of knots in B-splines and P-splines. These
MCMC approaches, however, are computationally very expensive when applied to
large datasets.

This work proposes an alternative to MCMC for approximate inference on
Bayesian bridge semi-parametric regression with B-splines based on the
Automatic Differentiation Variational Inference (ADVI) by Kucukelbir . (2015)
and Kucukelbir . (2017). The proposed inference algorithm jointly estimates
all parameters in the model (including $\alpha_{j}$ and $\lambda_{j}$) and
drastically increases computational speed for large datasets in comparison
with MCMC implementations since stochastic gradient updates involved in the
underlying optimization process small batches of data at each iteration,
rather than the entire dataset as required by MCMC schemes. A key point in
variational inference concerns the choice of the variational family of
distributions, which needs to be tractable while at the same time flexible
enough to approximate the posterior distribution well. Armagan (2009) describe
a mean-field variational inference method for Bayesian bridge regression model
with approximate inference for the bridge parameter ($\alpha$). Alves . (2021)
also proposes a variational approach based on mean-field assumption for
Bayesian inference in regression models with splines, however it is restricted
to Lasso penalization. In contrast, the ADVI does not require the
oversimplifying mean-field assumption for the variational family, as required
by coordinate ascending variational inference (Blei ., 2017), producing very
close approximations to the full posterior distribution, as demonstrated in
the simulation study. Previous gradient based variational procedures, such as
Ranganath . (2014) or Kingma Welling (2013) could be used, but they imply more
restricted forms of dependence in the variational family and produce noisier
estimates of the gradients when compared to ADVI. More recent approaches, such
as Yin Zhou (2018), allow for more flexible variational families, however we
found the simpler ADVI approach to produce already very good results in terms
of computational speed and posterior approximation under the semi-parametric
Bayesian bridge model.

To summarize, the main contribution of the paper is the development of a full
Bayesian inference procedure based on variational inference for semi-
parametric regression with bridge penalization. As strengths of the proposed
methodology, we can list the following:

  1. 1.

The proposed inference method is flexible enough to capture the dependence
structure in the target posterior distribution since the variational family
does not require the mean-field assumption. More specifically, for the mean-
field approach, the joint variational distribution of the parameters is the
product of its marginals. On the other hand, ADVI admits dependence structure
given by a transformation of a multivariate normal distribution with full
covariance matrix.

  2. 2.

It enables full Bayesian inference on semi-parametric regression for large
datasets at drastically lower computing times when compared with more
traditional MCMC implementations due to batch processing of data by the
stochastic gradient updates specially for large data.

The remainder of the paper is structured as follows. The Bayesian bridge semi-
parametric model is described in section 2. Section 3 covers the MCMC scheme
by Mallick Yi (2018), used for comparison with the proposed ADVI for Bayesian
bridge semi-parametric regression defined in section 4. Section 5 concerns a
simulation study with focus on large data. Section 6 contains an application
to real data on energy charges in Brazil and, finally, section 7 presents the
conclusions and future works.

##  2 The Bayesian bridge model for semi-parametric regression

A Bayesian bridge model for multiple covariates, can be written in the form of
a hierarchical model as follows:

$(\mathbf{y}\mid\bm{\beta}_{0},\bm{\beta}_{1},\ldots,\bm{\beta}_{D},\phi)\sim\\\
\sim
N\left(\mathbf{X}_{0}\bm{\beta}_{0}+\sum^{D}_{j=1}\mathbf{X}_{j}\bm{\beta}_{j},\
\phi^{-1}\mathbf{I}_{n}\right),\\\
(\beta_{jk}\mid\lambda_{j},\phi,\alpha_{j})\overset{\mathrm{iid}}{\sim}GG\left(0,\
\lambda^{-\frac{1}{\alpha_{j}}}_{j}\phi^{-\frac{1}{2}},\ \alpha_{j}\right),\\\
\ j=1,\ldots,D,\ k=1,\ldots,K_{j}.$ (3)

where $\mathbf{y}\in\mathbb{R}^{n},$ $N(\bm{\mu},\bm{\Sigma})$ denotes the
multivariate Gaussian distribution with mean $\bm{\mu}$ and covariance matrix
$\bm{\Sigma},$ $GG(\mu,\sigma,\alpha)$ denotes the generalized Gaussian
distribution (Subbotin, 1923) with mean $\mu$, scale $\sigma$ and shape
$\alpha$. The generalized Gaussian density evaluated at $x\in\mathbb{R}$ is
denoted by

$GG(x\mid\mu,\sigma,\alpha)=\frac{\alpha}{2\sigma\Gamma(\alpha^{-1})}\exp\left\\{-\left(\frac{\mid
x-\mu\mid}{\sigma}\right)^{\alpha}\right\\}.$

Equation (3) therefore implies the prior probability density function

$p(\bm{\beta}_{jk}\mid\lambda_{j},\phi,\alpha_{j})=\\\
=\frac{\alpha_{j}\lambda^{\frac{1}{\alpha_{j}}}_{j}\phi^{\frac{1}{2}}}{2\Gamma(\alpha^{-1}_{j})}\exp\left\\{-\lambda_{j}\left(\phi^{\frac{1}{2}}\mid\beta_{j}\mid\right)^{\alpha_{j}}\right\\},$
(4)

so that maximizing the posterior density
$p(\bm{\beta}\mid\mathbf{y},\bm{\lambda},\phi,\bm{\alpha},\bm{\beta}_{0})$
over $\bm{\beta}$ and $\bm{\beta}_{0}$ for fixed
$(\bm{\lambda},\phi,\bm{\alpha})$ is equivalent to solving (2), where
$\bm{\beta}=(\bm{\beta}_{1},\ldots,\bm{\beta}_{D}),\
\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{D}),\
\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{D})$. The particular form of dependency
of $\beta_{jk}$ on $\lambda_{j}$ and $\phi$ expressed by equation (4) avoids
undesired multimodality of the posterior distribution as pointed out by Park
Casella (2008).

The proposed semi-parametric model splits the regressors into matrices
$\mathbf{X}_{j},\ j=1,\ldots,D$ of covariates with their own penalized effects
and a matrix $\mathbf{X}_{0}$ of covariates with unpenalized effects. In the
proposed formulation, $\mathbf{X}_{j}$ represents a spline basis function that
implies non-parametric effects to the corresponding covariates. B-splines were
chosen as the set of spline basis functions to build $\mathbf{X}_{j}$ with the
bridge penalization acting on the B-spline coefficients. Naturally, other
types of spline basis functions could also be used.

Within the Bayesian framework, it is possible to provide joint (approximate)
posterior estimates for all parameters, including $\bm{\alpha}$ and
$\bm{\lambda}$ by specifying a prior distribution on them, as will be shown
further in sections 3 and 4. Therefore, joint posterior inference can be done
in a single-step procedure.

##  3 MCMC on Bayesian bridge regression

Both Polson . (2014) and Mallick Yi (2018) propose MCMC schemes for Bayesian
bridge regression for fixed $\alpha_{j}$ based on different data augmentation
schemes. The algorithm by Mallick Yi (2018) was chosen for comparison with the
proposed ADVI for Bayesian bridge, mostly because of the simplicity of its
gamma-uniform variable augmentation in comparison with the one proposed by
Polson . (2014). Accordingly, one can recover (3) by specifying

$\displaystyle(u_{jk}\mid\alpha_{j},\lambda_{j})\sim
Ga\left(\frac{1}{\alpha_{j}}+1,\ \lambda_{j}\right),\ \ \ k=1,\ldots,K_{j},$
$\displaystyle(\beta_{jk}\mid u_{jk},\phi,\alpha_{j})\sim
Unif\left(-u^{\frac{1}{\alpha_{j}}}_{jk}\phi^{-\frac{1}{2}},\
u^{\frac{1}{\alpha_{j}}}_{jk}\phi^{-\frac{1}{2}}\right),$ $\displaystyle\ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ j=1,\ldots,D.$ (5)

In fact,

$\displaystyle p(\beta_{jk}$ $\displaystyle\mid\lambda_{j},\alpha_{j},\phi)=$
$\displaystyle=\int^{\infty}_{0}p(\beta_{jk}\mid
u_{jk},\phi,\alpha_{j})p(u_{jk}\mid\alpha_{j},\lambda_{j})du_{jk}$
$\displaystyle\propto\exp\\{-\lambda_{j}\phi^{\frac{\alpha_{j}}{2}}\mid\beta_{jk}\mid^{\alpha_{j}}\\}.$

The model specification is completed by assigning independent priors $\phi\sim
Ga(a_{\phi},b_{\phi}),$ $\bm{\beta}_{0}\sim
N(\bm{\mu}_{0},\Sigma_{0}),\lambda_{j}\sim Ga(a_{\lambda},b_{\lambda}),$ where
$Ga(a,b)$ denotes the Gamma distribution with mean $a/b$ and variance
$a/b^{2}$. The prior distribution on $\alpha_{j}$ was defined by taking
$\alpha_{j}=2.5\eta_{j}$ where $\eta_{j}\sim Beta(a_{\eta},b_{\eta})$. The
upper limit for $\alpha_{j}$ in this work is fixed arbitrarily at 2.5 so that
the lasso and the ridge penalizations can be approximately represented as
special cases of bridge regression.

Gibbs sampler can be easily implemented for posterior inference on the
Bayesian bridge regression model. Due to conjugacy results, the full
conditional distributions of all parameters, except $\alpha_{j}$, are
analytically available. Details concerning the MCMC implementation are
described in the appendix A. In section A.1, the Gibbs sampler approach
considering $\alpha_{j}$ fixed is described. To include $\alpha_{j}$ in the
sampling scheme, one can carry out a Metropolis-Hastings step with a
transformed random walk proposal as described in section A.2.1 (when the
auxiliary vector $\mathbf{u}$ is marginalized) or in section A.2.2 when
$\mathbf{u}$ is not marginalized). It was found that the marginalized random
walk proposal produce well mixing Markov chains without much effort to tune
the proposal variance. Naturally, other Metropolis-Hastings proposals could be
used.

##  4 Variational Inference

The variational approach searches among a predefined variational family
$\mathcal{Q}=\\{q_{\bm{\psi}}(\bm{\theta}):\ \bm{\psi}\in\mathcal{V}\\}$
(where its members are densities indexed by the variational parameter
$\bm{\psi}$) for the density $q_{\bm{\psi}^{*}}(\bm{\theta})$ that best
approximates the posterior $p(\bm{\theta}\mid\mathbf{y})$ in terms of
Kullback-Leibler divergence. For the proposed Bayesian semi-parametric bridge
regression,
$\bm{\theta}=(\bm{\beta},\bm{\gamma},\phi,\lambda,\alpha)^{\top}\in\mathbb{R}^{p_{X}}\times\mathbb{R}^{p_{Z}}\times(0,\
1)^{2}\times(0,\ 2.5).$

Therefore, the variational objective is

$\operatorname*{arg\,min}_{\bm{\psi}\in\mathcal{V}}KL\left(q_{\bm{\psi}}(\bm{\theta})\mid\mid
p(\bm{\theta}\mid\mathbf{y})\right),$

which is equivalent to maximize the evidence lower bound (ELBO), i.e.
$\operatorname*{arg\,max}_{\bm{\psi}\in\mathcal{V}}ELBO(\mathbf{y},\bm{\psi}),$
where

$\displaystyle E$ $\displaystyle
LBO(\mathbf{y},\bm{\psi}):=\mathbb{E}_{q_{\bm{\psi}}}[\log
p(\mathbf{y},\bm{\theta})-\log q_{\bm{\psi}}(\bm{\theta})]$
$\displaystyle=\int\left[\log p(\mathbf{y},\bm{\theta})-\log
q_{\bm{\psi}}(\bm{\theta})\right]q_{\bm{\psi}}(\bm{\theta})\ d\bm{\theta}.$

One way to maximize the ELBO is to calculate its gradient
$\nabla_{\bm{\psi}}ELBO(\mathbf{y},\bm{\psi})$ and use it in a stochastic
gradient ascend based algorithm, such as Adam (Kingma Ba (2014)), Adagrad
(Duchi ., 2011), AdaDelta (Zeiler, 2012) and others. The issue with such
gradient ascending methods is evaluating the expectation in the ELBO, which is
often intractable. In this case, a common solution involves rewriting the ELBO
in a way that the gradient operator can switch order with the expectation so
the resulting $\nabla_{\bm{\psi}}ELBO(\mathbf{y},\bm{\psi})$ is an expected
value with respect to the variational distribution. The resulting expectation
can then be estimated via Monte Carlo methods.

This work implements ADVI (section 4.3) to the semi-parametric Bayesian bridge
regression described in section 2. The variational approaches by Ranganath .
(2014) and by Kingma Welling (2013) are briefly presented for comparison and
contextualization. A reader interested in variational inference for semi-
parametric regression models may see (Luts Wand, 2015), (Menictas Wand, 2015),
(Ong ., 2017) and (Wand, 2017) to mention a few.

###  4.1 Score method

Ranganath . (2014) introduced the black box variational inference (BBVI) in
which the score method (also known as ”log trick”) allows the gradient of the
ELBO to be written as an expectation with repect to the variational
distribution. The resulting expectation can be estimated by Monte Carlo as
follows

$\nabla_{\bm{\psi}}ELBO(\mathbf{y},\bm{\psi})=\\\
=\int\nabla_{\bm{\psi}}\left\\{[\log p(\mathbf{y}\mid\bm{\theta})+\log
p(\bm{\theta})+\right.\\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.-\log
q_{\bm{\psi}}(\bm{\theta})]q_{\bm{\psi}}(\bm{\theta})\right\\}d\bm{\theta}\\\
=\mathbb{E}_{q_{\bm{\psi}}}[(\log p(\mathbf{y}\mid\bm{\theta})+\log
p(\bm{\theta})-\log q_{\bm{\psi}}(\bm{\theta}))\times\\\ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \times\nabla_{\bm{\psi}}\log q_{\bm{\psi}}(\bm{\theta})]\\\
\approx\frac{1}{M}\sum^{M}_{m=1}\left[\log
p(\mathbf{y}\mid\bm{\theta}^{(m)})+\log p(\bm{\theta}^{(m)})+\right.\\\
\left.-\log q_{\bm{\psi}}(\bm{\theta}^{(m)})\right]\nabla_{\bm{\psi}}\log
q_{\bm{\psi}}(\bm{\theta}^{(m)}),$ (6)

where
$\bm{\theta}^{(m)}\overset{\mathrm{iid}}{\sim}q_{\bm{\psi}}(\bm{\theta}),\
m=1\ldots,M.$

Evaluating (6) requires (i) the prior, likelihood and variational densities to
be analytically available and (ii) the ability to draw samples from the
variational distribution $q_{\bm{\psi}}(\bm{\theta})$. However, in many
applications, the Monte Carlo estimate from (6) has high variance (Paisley .,
2012), even when control variates and Rao-Blackwellization are used to reduce
variance. It is simple to mitigate the mean field hypothesis in the BBVI
method to some extent by grouping blocks of components within $\bm{\theta}$
and specifying a dependent multivariate variational distribution within each
block. However specifying a single dependent distribution for the entire
$\bm{\theta}$ is typically hard, specially when its components lie in
different subspaces of $\mathbb{R}$ (so a multivariate Gaussian or Student-t
would be inappropriate choices for $q_{\bm{\psi}}$). We found it crucial to
account for full dependence structure in the variational family for the case
of semi-parametric Bayesian bridge model. In Appendix B, we provide the
specifics on how one can apply the BBVI algorithm to the Bayesian bridge
model.

###  4.2 Reparameterization gradient

An alternative way to derive Monte Carlo estimates for the gradient of the
ELBO is described in Kingma Welling (2013). The authors present the
reparameterization trick, which assumes that the parameter vector
$\bm{\theta}$ can be analytically written as
$\bm{\theta}=T_{\bm{\psi}}(\bm{\epsilon})$ where $T_{\bm{\psi}}$ is a
differentiable deterministic transformation involving the variational
parameters $\bm{\psi}$ and $\bm{\epsilon}$ is a random noise required to have
a closed form distribution $q^{*}(\bm{\epsilon})$ that does not depend on
$\bm{\psi}$ and is easy to sample from. For example, if $\bm{\theta}\sim
N(\mathbf{m},\mathbf{S}),$ then we can write
$\bm{\beta}=T_{\mathbf{m},\mathbf{S}}(\bm{\epsilon})=\mathbf{m}+\mathbf{L}\bm{\epsilon},$
with $\mathbf{L}$ being the Choleskey decomposition of the covariance matrix
$\mathbf{S}$ (ie, $\mathbf{L}\mathbf{L}^{\top}=\mathbf{S}$) and
$\bm{\epsilon}\sim N(\bm{0},\mathbf{I}).$ Another example of
reparameterization is the log-Normal distribution: if $\theta\sim\log
N(m,s^{2}),$ one can write
$\theta=T_{m,s}(\epsilon)=\exp\\{m+s\times\epsilon\\}$ with $\epsilon\sim
N(0,1).$

The reparameterization of $\bm{\theta}$ allows one to switch the order of the
gradient with the expectation when deriving the updating equations for
optimizing the ELBO. The gradient of the ELBO under reparameterization becomes

$\nabla_{\bm{\psi}}ELBO(\bm{\psi})=\\\ =\mathbb{E}_{\bm{\epsilon}\sim
q^{*}}\left[\nabla_{\bm{\psi}}\log
p(\mathbf{y},T_{\bm{\psi}}(\bm{\epsilon}))\right]+\\\
-\mathbb{E}_{\bm{\theta}\sim
q_{\bm{\psi}}(\bm{\theta})}[\nabla_{\bm{\psi}}\log
q_{\bm{\psi}}(\bm{\theta})].$

If a random sample $\bm{\epsilon}^{(1)},\ldots,\bm{\epsilon}^{(M)}$ is drawn
from $q^{*}(\bm{\epsilon}),$ the Monte Carlo estimate for the gradient is
obtained as

$\widetilde{\nabla}_{\bm{\psi}}ELBO(\bm{\psi})=\\\
\frac{1}{M}\sum^{M}_{m=1}\nabla_{\bm{\psi}}\log
p(\mathbf{y},T_{\bm{\psi}}(\bm{\epsilon}^{(m)}))\\\
-\frac{1}{M}\sum^{M}_{m=1}\nabla_{\bm{\psi}}\log
q_{\bm{\psi}}(T_{\bm{\psi}}(\bm{\epsilon}^{(m)})).$ (7)

The variance of the Monte Carlo estimates for the reparameterization gradients
tend to exhibit lower variance than the BBVI Monte Carlo estimates, but the
difficulty of assuming a full dependence structure on the variational family
for different entries of $\bm{\theta}$ persists. In Appendix B an
implementation of the reparameterization method for variational inference on
the proposed semi-parametric Bayesian bridge regression model is briefly
described.

###  4.3 Automatic differentiation variational inference

The implementation of ADVI method from Kucukelbir . (2015) and Kucukelbir .
(2017) to the proposed semi-parametric Bayesian bridge model is briefly
described in this section. The ADVI shares similarities with the
reparameterization method of Kingma Welling (2013) since it also makes use of
reparameterization in order to write the gradient of the ELBO as an
expectation to be approximated by Monte Carlo. The key distinction is that the
method works in a transformed parameter space that is (ideally) suitable to be
modeled as a multivariate normal distribution.

Suppose the original parameters $\bm{\theta}$ vary on a subset
$\Theta\subset\mathbb{R}^{d}$ with $\Theta\neq\mathbb{R}^{d}.$ This happens
for example when one or more entries of $\bm{\theta}$ lie in constrained
subsets of $\mathbb{R},$ say $\mathbb{R}^{+}$ or the interval $(0,1)$ for
instance. We consider the transformed parameter vector
$\bm{\xi}=T(\bm{\theta})$ where $T:\Theta\rightarrow\mathbb{R}^{d}$ is a is
diffeomorphism map (differentiable and invertible transformation) such that
$\bm{\xi}$ lies in $\mathbb{R}^{d},$ with no restrictions. For the mean-field
approach, the joint variational distribution of the parameters is the product
of its marginals:
$q_{\bm{\psi}}(\bm{\theta})=\prod^{K}_{j=1}q_{\psi_{j}}(\theta_{j})$ for
$\bm{\theta}=(\theta_{1},\ldots,\theta_{K})$. On the other hand, ADVI admits
dependence structure given by $\bm{\theta}=T^{-1}(\bm{\xi})$ where
$\bm{\xi}\sim N(\mathbf{m},\bm{\Sigma})$ with full covariance matrix
$\bm{\Sigma}$.

In the case of the Bayesian bridge model, we have the original parameters
$\bm{\theta}=(\bm{\beta}_{0},\bm{\beta}_{1},\ldots,\bm{\beta}_{D},\phi,\lambda_{1},\ldots,\lambda_{D},\alpha_{1},\ldots,\alpha_{D})^{\top}\in\mathbb{R}^{K_{0}+K_{1}+\ldots+K_{D}}\times(0,\
+\infty)^{D+1}\times(0,\ 2.5)^{D}$. A possible choice for $T$ is

$T(\bm{\theta})=(\bm{\beta}_{0},\bm{\beta}_{1},\ldots,\bm{\beta}_{D},\log\phi,\log\lambda_{1},\ldots,\log\lambda_{D},\\\
\log\frac{\alpha_{1}}{2.5-\alpha_{1}},\ldots,\log\frac{\alpha_{D}}{2.5-\alpha_{D}})^{\top}.$

In the remainder of this section, we will denote
$\bm{\xi}=(\bm{\xi}_{\bm{\beta}_{0}},\bm{\xi}_{\bm{\beta}_{1}},\ldots,\bm{\xi}_{\bm{\beta}_{d}},\xi_{\phi},\xi_{\lambda_{1}},\ldots,\xi_{\lambda_{D}},\xi_{\alpha_{1}},\ldots,\xi_{\alpha_{D}})^{\top},$
where $\bm{\xi}_{\bm{\beta}_{0}}=\bm{\beta}_{0},\
\bm{\xi}_{\bm{\beta}_{j}}=\bm{\beta}_{j},\ \xi_{\phi}=\log{\phi},\
\xi_{\lambda_{j}}=\log{\lambda_{j}},\
\xi_{\alpha_{j}}=\log\frac{\alpha_{j}}{2.5-\alpha_{j}},$ for $j=1,\ldots,D$.
The joint distribution $p(\mathbf{y},\bm{\theta})$ is defined as in the right
hand side of equation (10) with likelihood
$(\mathbf{y}|\bm{\beta}_{0},\bm{\beta}_{1},\ldots,\bm{\beta}_{D},\phi)\sim
N(\sum^{D}_{j=0}\mathbf{X}_{j}\bm{\beta}_{j},\phi^{-1}\mathbf{I}_{n})$ and
prior
$p(\bm{\theta})=p(\bm{\beta}_{0})p(\phi)\prod^{D}_{j=1}p(\bm{\beta}_{j}|\lambda_{j},\phi,\alpha_{j})p(\lambda_{j})p(\alpha_{j}),$
$(\beta_{j}|\lambda_{j},\phi,\alpha_{j})\overset{\mathrm{iid}}{\sim}GG(0,\lambda_{j}^{-\frac{1}{\alpha_{j}}}\phi^{-\frac{1}{2}},\alpha_{j})$,
$\phi\sim Ga(a_{\phi},b_{\phi}),$ $\bm{\beta}_{0}\sim
N(\bm{\mu}_{0},\Sigma_{0}),\lambda_{j}\sim Ga(a_{\lambda},b_{\lambda})$ with
$\alpha_{j}=2.5\eta_{j}$ where $\eta_{j}\sim Beta(a_{\eta},b_{\eta})$ for
$j=1,\ldots,D.$

The ADVI method redefines the joint model density in terms of the joint
distribution of $\mathbf{y}$ and the transformed parameter $\bm{\xi}$, here
denoted as $\tilde{p}(\mathbf{y},\bm{\xi})$ to distinguish it from the
original joint density $p(\mathbf{y},\bm{\theta})$. It follows that

$\displaystyle\tilde{p}(\mathbf{y},\bm{\xi})$
$\displaystyle=p(\mathbf{y},\bm{\theta})\bigg{\rvert}_{\bm{\theta}=T^{-1}(\bm{\xi})}\times\mid
J_{T^{-1}(\bm{\xi})}\mid$
$\displaystyle=p(\mathbf{y},T^{-1}(\bm{\xi}))\times\mid
J_{T^{-1}(\bm{\xi})}\mid,$

where $J_{T^{-1}(\bm{\xi})}$ represents the Jacobian of the inverse
transformation $T^{-1}:\mathbb{R}^{d}\rightarrow\Theta,$

$T^{-1}(\bm{\xi})=\left(\bm{\xi}_{\bm{\beta}_{0}},\bm{\xi}_{\bm{\beta}_{1}},\ldots,\bm{\xi}_{\bm{\beta}_{D}},e^{\xi_{\phi}},e^{\xi_{\lambda_{1}}},\ldots,e^{\xi_{\lambda_{D}}},\right.\\\
\left.\frac{2.5}{1+e^{-\xi{\alpha_{1}}}},\ldots,\frac{2.5}{1+e^{-\xi{\alpha_{D}}}}\right)^{\top}$

which in this case is
$|J_{T^{-1}(\bm{\xi})}|=2.5^{D}e^{\xi_{\phi}}\prod^{D}_{j=1}e^{\xi_{\lambda_{j}}}e^{-\xi_{\alpha_{j}}}/(1+e^{-\xi_{\alpha_{j}}})^{2}$.
A multivariate Gaussian variational distribution $q(\bm{\xi}\mid\mathbf{L},\
\mathbf{m})=N(\bm{\xi}\mid\mathbf{m},\ \mathbf{L}\mathbf{L}^{\top})$ is
specified for $\bm{\xi}$ and the variational parameters are
$\bm{\psi}=(\mathbf{m},\mathbf{L})$. By reparameterizing
$\bm{\xi}=\mathbf{m}+\mathbf{L}\ \bm{\epsilon},$ the ADVI method enables
calculation of the gradient of the ELBO as an expectation, which can be
approximated via Monte Carlo:

$\nabla_{\bm{\psi}}ELBO(\mathbf{y},\bm{\psi})=\\\
=\nabla_{\bm{\psi}}\mathbb{E}_{q_{\bm{\psi}}(\bm{\xi})}\left[\log\tilde{p}(\mathbf{y},\bm{\xi})-\log
q_{\bm{\psi}}(\bm{\xi})\right]\\\
=\nabla_{\bm{\psi}}\mathbb{E}_{\bm{\epsilon}\sim N(\bf
0,\mathbf{I})}\left[\log p(\mathbf{y},T^{-1}(\bm{\xi}))+\right.\\\
\left.+\log\mid J_{T^{-1}(\bm{\xi})}\mid-\log N(\bm{\xi}\mathchar
24635\relax\;\
\mathbf{m},\mathbf{L}\mathbf{L}^{\top})\bigg{\rvert}_{\bm{\xi}=\mathbf{m}+\mathbf{L}\bm{\epsilon}}\right]\\\
=\mathbb{E}_{\bm{\epsilon}\sim N(\bf
0,\mathbf{I})}\left\\{\nabla_{\bm{\psi}}\left[\log
p(\mathbf{y},T^{-1}(\bm{\xi}))+\log\mid
J_{T^{-1}(\bm{\xi})}\mid\right.\right.\\\ \left.\left.-\log
N(\bm{\xi}\mathchar 24635\relax\;\
\mathbf{m},\mathbf{L}\mathbf{L}^{\top})\bigg{\rvert}_{\bm{\xi}=\mathbf{m}+\mathbf{L}\bm{\epsilon}}\right]\right\\}\\\
\approx\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\left[\log
p(\mathbf{y},T^{-1}(\bm{\xi}))+\log\mid J_{T^{-1}(\bm{\xi})}\mid\right.\\\
\left.-\log N(\bm{\xi}\mathchar 24635\relax\;\
\mathbf{m},\mathbf{L}\mathbf{L}^{\top})\bigg{\rvert}_{\bm{\xi}=\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(\ell)}}\right]$
(8)

where $\bm{\epsilon}^{(\ell)}\overset{\mathrm{iid}}{\sim}N(\bf 0,\mathbf{I})$,
$\ell=1,\ldots,M$.

One can also easily compute the stochastic gradient approximation of (8) for a
random minibatch $\tilde{\mathbf{y}}=\left(y_{i_{1}},\ldots,y_{i_{K}}\right)$
where $i_{1},\ldots,i_{K}$ represent a random subset of size $K$ from
$\\{1,\ldots,n\\}$. There are many ways in which one can draw random
minibatches $\tilde{\mathbf{y}}$, as long as the Monte Carlo estimate for the
gradient based on $\tilde{\mathbf{y}}$ is unbiased for the full gradient based
on $\mathbf{y}$. In this work, we randomly permute the indexes of $\mathbf{y}$
at the beginning of each epoch and pick blocks of $M$ consecutive observations
to form the random batches at each iteration. The Monte Carlo estimate for
stochastic gradient of the ELBO becomes

$\widetilde{\nabla}_{\bm{\psi}}ELBO(\tilde{\mathbf{y}},\bm{\psi})\\\
=\frac{n}{KM}\sum^{M}_{\ell=1}\left[\nabla_{\bm{\psi}}\log
p(\tilde{\mathbf{y}}\mid
T^{-1}(\bm{\xi}))\bigg{\rvert}_{\bm{\xi}=\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(\ell)}}\right]\
+\ \\\ +\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\left[\log
p(T^{-1}(\bm{\xi}))+\log\mid J_{T^{-1}(\bm{\xi})}\mid\right.\\\ \left.-\log
N(\bm{\xi}\mathchar 24635\relax\;\
\mathbf{m},\mathbf{L}\mathbf{L}^{\top})\bigg{\rvert}_{\bm{\xi}=\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(\ell)}}\right].$
(9)

In equation (9), the terms $p(\tilde{\mathbf{y}}\mid T^{-1}(\bm{\xi}))$ and
$p(T^{-1}(\bm{\xi}))$ denote the densities
$p(\tilde{\mathbf{y}}\mid\bm{\theta})$ and $p(\bm{\theta})$ evaluated at
$\bm{\theta}=T^{-1}(\bm{\xi}).$ The multiplicative term $\frac{n}{K}$ rescales
the minibatch gradient to make it an unbiased estimate for the full posterior
$p(\bm{\theta}\mid\mathbf{y})$ based on $n$ observations rather than the batch
size posterior, which is based on $M$ data points only. Algorithm 1 summarizes
the steps of the proposed ADVI algorithm for semi parametric Bayesian bridge.
The expression in equation (9) is central to Algorithm 1 so its details are
fully developed in Appendix D. The update equation in line 5 corresponds to
the basic stochastic gradient ascent, although any gradient based method such
as Adam, Adagrad and others could be used.

Input :

Data: $\mathbf{y},\mathbf{X},\mathbf{Z}$,

Learning rate: $\delta\in(0,1)$,

Batch size: $K\in\mathbb{N}$

Number of Monte Carlo samples: $M\in\mathbb{N},$

Number of iterations: $I\in\mathbb{N}$,

Number of samples from $q_{\bm{\psi}}(\bm{\theta})$: $S\in\mathbb{N}$,

Batch size: $B\in\mathbb{N}.$

.

.

1for  _ $\mbox{iter}\in\\{1,\ldots,I\\}$_ do

2 Pick a random batch
$\tilde{\mathbf{y}}=(\tilde{y}_{i_{1}},\dots,\tilde{y}_{i_{B}})$ from
$\mathbf{y}=(y_{1},\ldots,y_{n}).$

3Sample
$\bm{\epsilon}^{(\ell)}\overset{\mathrm{iid}}{\sim}N(\bm{0},\mathbf{I}),\
\ell\in\\{1,\ldots,M\\}.$

4Evaluate $\widetilde{\nabla}_{\bm{\psi}}ELBO(\tilde{\mathbf{y}},\bm{\psi})$
from equation (9) and Appendix D.

5Gradient ascent iteration:
$\bm{\psi}\leftarrow\bm{\psi}+\delta\times\widetilde{\nabla}_{\bm{\psi}}ELBO(\tilde{\mathbf{y}},\bm{\psi}).$

end for

6Sample $\bm{\epsilon}^{(s)}\overset{\mathrm{iid}}{\sim}N(\bm{0},\mathbf{I}),\
\ s=1,\ldots,S.$

7Generate samples $\bm{\theta}^{(s)}\sim q_{\bm{\psi}}(\bm{\theta})$ by taking
$\bm{\theta}^{(s)}=T^{-1}(\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(s)}),\ \
s=1,\ldots,S.$

Output :

$\bm{\psi}=(\mathbf{m},\mathbf{L})$

$\bm{\theta}^{(s)}=T^{-1}(\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(s)}),\ \
s=1,\ldots,S.$

Algorithm 1 ADVI algorithm for semi parametric Bayesian bridge

The ADVI method was found to suit well to the Bayesian bridge semi-parametric
proposed model, both in terms of computational speed and also in terms of
quality of posterior approximation, due to the full dependence structure
incorporated in the variational family. These findings are documented in the
simulation study and in the application to real data presented in sections 5
and 6.

##  5 Numerical experiments

In this section, posterior estimates under the MCMC approach based on Mallick
Yi (2018) as described in section 3 and the proposed ADVI for semi-parametric
Bayesian bridge, as described in section 4.3 are compared under different
scenarios, focusing on large datasets. We also compare results with the
Bayesian formulations for generalized additive models provided by the R
packages: mgcv (S. Wood 2012; S.N. Wood 2017), brms (Bürkner, 2017), and INLA
(Rue ., 2009). Finally, we exemplify the proposed procedure in a synthetic
dataset with multiple covariates.

###  5.1 Scenario 1: small datasets

In this scenario, multiple small datasets are simulated and used to draw
comparisons between the proposed ADVI for semi-parametric regression as
described in Section 4.3 and the MCMC scheme described in 3. One hundred
datasets were simulated, with $n=100$ observations each, according to the
model defined by $y_{i}\stackrel{{\scriptstyle
ind}}{{\sim}}N(f_{\bm{\beta}}(x_{i}),\ \sigma^{2}),\ i=1,\ldots,n,$ where
$f_{\bm{\beta}}(x)$ represents the non-parametric effects of the covariate
$x$, modeled by a B-splines with coefficients $\bm{\beta}.$ The observational
variance parameter $\sigma^{2}$ was fixed at 1 when simulating the data and
the nodes for the B-spline basis were regularly spaced from -0.066 to 1.066
with 0.033 units of consecutive distance, which implied 34 B-spline
coefficients. We also fixed the covariates $x_{i},\ i=1,\ldots,100$ on a
regular grid over the interval (0, 1). The B-spline coefficients $\beta_{k},\
k=1,\ldots,34$ were sampled as follows:
$\beta_{1},\ldots,\beta_{10}\stackrel{{\scriptstyle iid}}{{\sim}}Nt(5,2),$
$\beta_{16},\ldots,\beta_{25}\stackrel{{\scriptstyle iid}}{{\sim}}Nt(10,2),$
$\beta_{31},\ldots,\beta_{34}\stackrel{{\scriptstyle iid}}{{\sim}}Nt(4,0.25)$
and $\beta_{k}=0\ \forall k\in\\{11,\ldots,15\\}\cup\\{26,\ldots,30\\}$ and
the values were truncated to the nearest integers to generate the splines. The
true simulated values for $\beta_{k},\ k\in\\{1,\ldots,34\\}$ can be seen in
Figure 4 and are shared for all 100 simulated datasets, which implies only one
underlying B-splines curve, along which all the replicas for
$\mathbf{y}=(y_{1},\ldots,y_{n})$ are simulated. Since the B-splines do not
include intercept, the parameter $\bm{\gamma}$ does not need to be included in
the model.

The prior hyperparameters specified for Bayesian inference under both MCMC and
Variational Bayes (VB) approach via ADVI were fixed as
$a_{\alpha}=a_{\lambda}=a_{\phi}=b_{\alpha}=b_{\lambda}=b_{\phi}=1,$ therefore
representing vague prior knowledge about $\alpha,\ \lambda$ and $\phi$. Figure
2 shows the posterior inference for $\mu(x)=f_{\beta}(x),\ x\in(0,\ 1)$ based
on the first replica of the simulated data. MCMC and VB produce very similar
posterior estimates for $\mu(x)$, with noticeable differences only at the
begining and end of the series. The same holds for the chosen alternative
inference methods: Integrated Nested Laplace Approximation (INLA) by Rue .
(2009), the brms R package by Bürkner (2017) (which implements MCMC within
Stan) and the smooth.spline R function.

Next, Figure 3 compares the ADVI with bridge penalization on the spline basis
coefficients as in equation (2), to its analogous Bayesian p-splines
formulation by Lang Brezger (2004) (MCMC) and to the off-the-shelf
smooth.spline R function. The objective is to compare the 3 distinct forms of
penalization: 1) directly on the coefficients (proposed ADVI), 2) on the
second order differences on consecutive coefficients (P-splines) and under no
penalization (smooth.spline). The fitted curves are very similar, with the
proposed ADVI for Bayesian bridge penalization being slightly less smooth than
the other two approaches. Furthermore, the observed mean absolute prediction
error when estimating $\mu=f_{\beta}(x)$ was lower when using the ADVI
(0.4263) when compared to Bayesian P-splines (0.5648) and smoothing splines
(0.5296).

Figure 4 shows that both methods produce very similar posterior marginal
distributions for the spline coefficients and that the marginals capture the
true values of the coefficients. Figure 5 shows MCMC and ADVI posterior
marginals for $\lambda$, $\phi$ and $\alpha$. The marginal posterior
distributions for $\phi$ and $\lambda$ under VB and MCMC are similar, while
$\alpha$ posterior estimates via VB are less dispersed in comparison with MCMC
marginal posterior densities. Finally, joint posterior uncertainties are also
very similar under both aproaches, as examplified in Figures 6 and 7, except
for the bridge penalization parameter $\alpha$. Despite the discrepancies
regarding estimation of $\alpha$ by the VB and MCMC methods, there is very
little difference regarding the goodness of fit to the simulated data (2)
which might indicate that the data brings little information about $\alpha$.

Furthermore, similar conclusions can be drawn from the other 99 replicates.
Posterior point estimates produced by the proposed ADVI and by the baseline
MCMC on all 100 replicas are compared in Figure 8 and Figure 15 in the
appendix. Point estimates for $\beta_{k},\ k=1,\ldots,34$ are very similar
under both methods, except for $\beta_{34}.$ However, a closer look reveals
that the 34-th column of the B-spline regression matrix $\mathbf{X}$ has all
entries equal to 0, except the last one, which makes $\beta_{34}$ very hard to
estimate as it has very little impact on the mean function $f_{\beta}(x).$ The
level of agreement between the point estimates produced by MCMC and the
proposed ADVI for $\lambda$ and $\alpha$ are not as high as for the regression
coefficients $\beta_{k},\ k=1,\ldots,34$.

Figure 2: ADVI and MCMC posterior confidence bands for the true B-spline curve
$\mu=f_{\bm{\beta}}(x)$ simulated in the first replica. Only one curve is
shown for the point estimates of $\mu$ for ease of visualization, since both
methods produced virtually indistinguishable curves. Figure 3:  Posterior
estimates for the true B-spline curve $\mu=f_{\bm{\beta}}(x)$ simulated in the
first replica obtained by ADVI and Bayesian P-splines (MCMC). The classical
point estimate produced by the R function smooth.spline is also shown. Figure
4: Marginal posterior distributions for $\beta_{k},\ k=1,\ \ldots,7$ based on
the first replicate of the simulated data. Red dots denote the true values of
the parameters. The marginal densities for the remaining coefficients can be
seen in Figure 14.  Figure 5: Marginal posterior distributions for $\phi$,
$\lambda$ and $\alpha$ based on the first replicate of the simulated data. Red
dashed curves represent VB and black curves represent MCMC marginal posterior
approximations. Figure 6: ADVI and MCMC bivariate joint posteriors for the
first 10 B-spline coefficients. Figure 7: ADVI and MCMC bivariate joint
posteriors for $\beta_{1},\beta_{2},\phi,\lambda,\alpha$. Figure 8: Point
estimates for all parameters under MCMC and VB for each one of the 100
replicas. The black lines represent the identity function. The scatterplots
for the remaining coefficients can be seen in Figure 15 in the appendix.

###  5.2 Scenario 2: varying sample sizes

Table 1: Number of batches ($n_{b}$), number of epochs ($n_{e}$) and total number of iterations ($n_{i}$) of the proposed VB algorithm according to the sample sizes ($n$) of the simulated datasets. $n$ | 1,000 | 10,000 | 50,000 | 100,000 | 500,000 | 1,000,000
---|---|---|---|---|---|---
$n_{b}$ | 1,000 | 1,000 | 1,000 | 10,000 | 10,000 | 10,000
$n_{e}$ | 2,000 | 200 | 100 | 100 | 100 | 100
$n_{i}$ | 2,000 | 2,000 | 5,000 | 1,000 | 5,000 | 10,000
|  |  |  |  |  |

In this section, the effect of sample sizes in computational times of the MCMC
and ADVI algorithms in the context of estimation of the semi-parametric
Bayesian bridge is investigated, focusing on large datasets. The
specifications for the simulated data are the same as described in section
5.1, except that now the sample sizes of the simulated data vary as
$n\in\\{10^{3},\ 10^{4},\ 5\times\ 10^{4},\ 10^{5},\ 5\times\ 10^{5},\
10^{6}\\}.$ For all six cases, the MCMC was ran for 5000 iterations (except
for n=1000, which required 50000 iterations until we could get reasonable
evidence of convergence), initializing $\bm{\beta}$ and $\phi$ at their OLS
estimates under no regularization. For the VB implementation, we always used
learning rate 0.01 for Adam and $M=100$ Monte Carlo samples to approximate
expectations in the calculation of the ELBO. Distinct batch sizes and number
of epochs were for each sample size according to Table 1, which implies 1000
to 10000 iterations in total depending on the value of $n$ to be compared with
the results obtained after 5000 iterations of MCMC.

Regarding computational time, the proposed ADVI algorithm for semi-parametric
Bayesian bridge regression is orders of magnitude faster than MCMC, as shown
in Table 2. ADVI remais faster even when running more iterations than MCMC and
also when VB times are adjusted to account for the 5000 iterations ran under
MCMC.

In summary, the proposed VB aproach is capacble of reaching convergence faster
than MCMC, using less epochs and providing an accurate approximation of the
posterior.

Table 2: Computational times (in seconds) for MCMC and ADVI implementations according to the size of the simulated data ($n$). *: In this case we ran 50000 mcmc iterations instead of 5000 because of lack of evidence for convergence under MCMC. $n$ | MCMC | ADVI | ADVI
---|---|---|---
|  |  | (5000 iterations)
$10^{3}$ | $219s^{*}$ | $6s$ | $15s$
$10^{4}$ | $64s$ | $6s$ | $15s$
$5\times 10^{4}$ | $294s$ | $15s$ | $15s$
$10^{5}$ | $618s$ | $13s$ | $65s$
$5\times 10^{5}$ | $2,744s$ | $76s$ | $76s$
$10^{6}$ | $4,570s$ | $143s$ | $72s$

###  5.3 Scenario 3: multiple covariates with non-parametric effects

In this scenario, a dataset with $n=1000$ observations
$(y_{i},x_{i1},x_{i2}),\ i=1,\ldots,n$ were simulated from the model
$y_{i}=\beta_{0}+f_{1}(x_{i1}\mathchar
24635\relax\;\tau_{1})+f_{2}(x_{i2}\mathchar
24635\relax\;\tau_{2})+\epsilon_{i},\ \epsilon_{i}\sim N(0,\sigma^{2}),$ where
$f_{i}(\cdot\mathchar 24635\relax\;\tau_{i})$ denotes a realization of a
1-dimensional Gaussian process with mean 0 and covariance function
$Cov:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R},\
Cov(z,z^{\prime})=\exp\left[-\tau^{-1}_{i}(z-z^{\prime})^{2}\right]$. In the
simulation, we used $\tau_{1}=1$ and $\tau_{2}=2$. The objective is to
investigate how the proposed semi-parametric bridge regression model performs
when there are more than one covariate affecting the response in a non-
parametric way. The results presented here assume $n_{1}=n_{2}=100$ knots for
each b-spline non-parametric effect.

Figure 9 shows that the underlying simulated surface for the mean of the
response variable $y$ as a function of $x_{1}$ and $x_{2}$ is accurately
estimated by the MCMC and VB approaches. It is worth to mention that the
adopted VB is much faster than MCMC.

When comparing the model formulation under MCMC, VB, and the point estimates
from GAM, Figure 10 shows that they yield virtually indistinguishable
estimates for the non-parametric effects. While GAM penalizes the second-order
differences of the basis functions’ coefficients, the bridge penalization used
in the MCMC and VB leads to the same fit.

Figure 9: Simulated mean of $y$ as a function of covariates $x_{1}$ and
$x_{2}$ according to section 5.3. Posterior point estimate for the average of
$y$ under the proposed MCMC approach. Figure 10: Comparison of MCMC and GAM
when fitting the simulated data from section 5.3.

##  6 Real data application

This section illustrates the use of the proposed ADVI inference procedure for
semi-parametric Bayesian bridge regression on a large real dataset. The data
consists of hourly measured Energy Charges starting from 2014-04-06 (10pm) to
2022-03-31 (11pm) averaged over stations in the Northern region of Brazil (see
Figure 11). In total, there are $n=69717$ obsevations. The dataset is
maintained by Operador Nacional do Sistema Elétrico (ONS) and can be obtained
at https://dados.ons.org.br/dataset/carga-energia .

The data exhibits strong seasonal patterns with multiple frequencies due to
periodicities in energy consumption according with time of the day, day of the
week, season of the year and possibly more. We model the seasonal harmonics as
the cosine and sine Fourier basis representation for weekly periodicity
(period = $24\times 7=168$) following West Harrison (1997). To capture overall
level changes in the series, we included a cubic B-spline with one knot at
every 100 hours for a total of 700 knots. The B-spline coefficients are
subjected to bridge penalization while the coefficients of the Fourier
harmonics basis functions are not penalized. In total, the resulting covariate
matrix has $p=868$ columns.

Figure 11: Full Energy Charge data (measured hourly).

The proposed ADVI inference approach took approximately 831 seconds (less than
14min) to complete 2000 iterations using 100 Monte Carlo samples to estimate
the gradient of the ELBO. On the other hand, MCMC takes 69239 seconds (19.23
hours) to run the same 2000 iterations. Posterior estimates for the beginning
and end of the series are shown in Figures 12 and 13.

Figure 12: VB 95% credibility bands and posterior mean for the average
response. Only the first 1000 observations of Energy Charge data are shown.

##  7 Conclusions and future work

The present work developed a variational inference procedure based on ADVI for
Bayesian inference in semi-parametric bridge regression models. The use of
small batches of data at each iteration of the training algorithm reduces
computational time in comparison with a more traditional MCMC approach. Full
Bayesian inference is preserved so joint uncertainty estimates for all model
parameters are available. It was verified in the simulation study that the
joint posterior is well approximated by the proposed variational family.

Figure 13: VB 95% credibility bands and posterior mean for the average
response. Only the last 1000 observations of Energy Charge data are shown.

Some directions for future research include (i) the extension to non-Gaussian
distribtutions for the response variable; (ii) expanding flexibility of the
variational family with other approaches such as the semi implicit variational
inference from Yin Zhou (2018) and (iii) considering different spline basis
functions to represent non-parametric effects of covariaties including GAM
models with tensor splines avoiding MCMC and RJMCMC. See details in Denison .
(1998), Dias Gamerman (2002), F. Li Villani (2013) and references therein.

##  8 Acknowledgements

This paper was partially supported by Fapesp Grants (R. Dias) 2018/04654-9,
(R. Dias and H. S. Migon) 2019/10800-0 and (R. Dias) 2019/00787-7.

##  Appendix A Details for the MCMC algorithm

###  A.1 Full conditionals

The original posterior distribution (without variable augmentation) is given
by

$p(\bm{\beta},\bm{\beta}_{0},\bm{\lambda},\phi,\bm{\alpha}\mid\mathbf{y})\propto
p(\mathbf{y}\mid\bm{\beta},\bm{\beta}_{0},\phi)p(\bm{\beta}_{0})\times\\\
\times
p(\bm{\beta}\mid\bm{\lambda},\phi,\bm{\alpha})p(\bm{\lambda})p(\phi)p(\bm{\alpha}),$
(10)

where $\bm{\lambda}=(\lambda_{1},\ldots,\lambda_{D})$ and
$\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{D})$.

The joint posterior distribution with the uniform-gamma variable augmentation
is expressed as

$p(\bm{\beta},\bm{\beta}_{0},\bm{\lambda},\phi,\mathbf{u},\bm{\alpha}\mid\mathbf{y})\propto
p(\mathbf{y}\mid\bm{\beta},\bm{\beta}_{0},\phi)p(\bm{\beta}_{0})\times\\\
\times
p(\bm{\beta}\mid\mathbf{u},\phi,\bm{\alpha})p(\mathbf{u}\mid\bm{\lambda},\bm{\alpha})p(\bm{\lambda})p(\phi)p(\bm{\alpha}).$

With $\bm{\alpha}$ fixed, all full conditional distributions are analytically
available, therefore a straightforward Gibbs sampler scheme can be
implemented. The procedure from Damien Walker (2001) is used to get a
numerically stable sampler for the truncated Gaussian and truncated
exponential distributions that appear in the Gibbs sampler scheme that
follows.

Posterior full conditional distribution of $\bm{\beta}_{j}$

$(\bm{\beta}_{j}\mid\mathbf{y},\bm{\beta}_{0},\phi,\mathbf{u},\bm{\alpha})\sim\\\
N\left(\
(\mathbf{X}^{\top}_{j}\mathbf{X}^{\top}_{j})^{-1}\mathbf{X}^{\top}_{j}(\mathbf{y}-\mathbf{X}_{0}\bm{\beta}_{0}-\sum_{k\neq
j}\mathbf{X}_{k}\bm{\beta}_{k}),\right.\\\
\left.(\mathbf{X}^{\top}_{j}\mathbf{X}_{j})^{-1}\phi^{-1}\right)\mathds{1}_{\mathcal{B}_{j}},$

where

$\mathcal{B}_{j}=\left\\{(\beta_{j1},\ldots,\beta_{j,K_{j}})\in\mathbb{R}^{K_{j}}:\right.\\\
\left.\lvert\beta_{jk}\rvert<u^{\frac{1}{\alpha_{j}}}_{jk}\phi^{-\frac{1}{2}},\
\forall k\in\\{1,\ldots,K_{j}\\}\right\\}$

and $N(\bm{\mu},\bm{\Sigma})\mathds{1}_{\mathcal{A}}$ denotes the multivariate
Gaussian distribution with mean $\bm{\mu}$ and covariance matrix
$\bm{\Sigma}$, truncated in the set $\mathcal{A}$.

Posterior full conditional distribution of $\bm{\beta}_{0}$

$(\bm{\beta}_{0}\mid\mathbf{y},\phi,\bm{\beta})\sim
N\left(\left(\mathbf{X}^{\top}_{0}\mathbf{X}_{0}\phi+\Sigma^{-1}_{0}\right)^{-1}\times\right.\\\
\left.\times\left(\phi\mathbf{X}^{\top}_{0}\left(\mathbf{y}-\sum^{D}_{j=1}\mathbf{X}_{j}\bm{\beta}_{j}\right)+\Sigma^{-1}_{0}\bm{\mu}_{0}\right)\right.,\\\
\left.\left(\mathbf{X}^{\top}_{0}\mathbf{X}_{0}\phi+\Sigma^{-1}_{0}\right)^{-1}\right)$

Posterior full conditional distribution of $\phi$

$(\phi\mid\mathbf{y},\bm{\beta},\bm{\beta}_{0},\mathbf{u},\bm{\alpha})\sim\\\
\sim Ga\left(\frac{n}{2}+\frac{1}{2}\sum^{D}_{j=1}K_{j}+a_{\phi},\
\frac{1}{2}RSS(\bm{\beta},\bm{\beta}_{0})+b_{\phi}\right)\\\
\mathds{1}_{\left(0,\
\min_{j,\ell}\left(u^{\frac{2}{\alpha_{j}}}_{j\ell}\mid\beta_{j\ell}\mid^{-2}\right)\right)}$

where
$RSS(\bm{\beta},\bm{\beta}_{0})=(\mathbf{y}-\sum_{j=1}^{D}\mathbf{X}_{j}\bm{\beta}_{j}-\mathbf{X}_{0}\bm{\beta}_{0})^{\top}(\mathbf{y}-\sum_{j=1}^{D}\mathbf{X}_{j}\bm{\beta}_{j}-\mathbf{X}_{0}\bm{\beta}_{0})$
denotes the residual sum of squares and $Ga(a,b)\mathds{1}_{\mathcal{A}}$
denotes the Gamma distribution with mean $a/b$ and variance $a/b^{2}$
truncated in the set $\mathcal{A}$.

Posterior full conditional distribution of $u_{j\ell}$

$(u_{j\ell}\mid\mathbf{y},\lambda_{j},\bm{\beta},\phi,\alpha_{j})\sim
Exp(\lambda_{j})\mathds{1}_{\left(\mid\beta_{j\ell}\mid^{\alpha_{j}}\phi^{\frac{\alpha_{j}}{2}},\
+\infty\right)},$

where $Exp(\lambda)\mathds{1}_{\mathcal{A}}$ denotes the exponential
distribution with rate parameter $\lambda$ truncated on $\mathcal{A}$.

Posterior conditional distribution of $\lambda$

When sampling from $\lambda_{j}$, we can marginalize out $\mathbf{u}_{j}$.
Such marginalization was found to improve mixing of the Markov chains.

$(\lambda_{j}\mid\mathbf{y},\bm{\beta}_{j},\phi,\alpha_{j})\sim\\\ \sim
Ga\left(a_{\lambda}+\frac{K_{j}}{\alpha_{j}},\
b_{\lambda}+\phi^{\frac{\alpha_{j}}{2}}\sum^{K_{j}}_{\ell=1}\mid\beta_{j\ell}\mid^{\alpha_{j}}\right).$

###  A.2 Estimation of $\alpha_{j}$

There is no possible choice of prior for $\alpha_{j}$ that leads to an
analytically available full conditional distribution on $\alpha_{j}$.
Therefore a Metropolis-Hastings (MH) step is proposed for $\alpha_{j}$. The MH
algorithm requires specification of a proposal distribution
$q(\alpha^{*}_{j}|\alpha^{(i)}_{j})$ where $\alpha^{(i)}_{j}$ denotes the
value of $\alpha_{j}$ in the current iteration $i$ and $\alpha^{*}_{j}$
denotes the proposed value for $\alpha_{j}$ according to the proposal density
$q(\cdot\mid\alpha^{(i)}_{j})$.

####  A.2.1 Marginalized MH proposal

When sampling $\alpha_{j}$, we consider the reduced parameter vector
$\bm{\theta}=(\bm{\alpha},\bm{\beta},\bm{\beta}_{0},\phi,\bm{\lambda})$
marginalizing out the auxiliary variables $\mathbf{u}$. Alternatively,
appendix A.2.2 describes a proposal $q(\cdot\mid\alpha^{(i)}_{j})$ based on
the full parameter vector
$(\bm{\alpha},\bm{\beta},\bm{\beta}_{0},\phi,\bm{\lambda},\mathbf{u})$
including the augmented variables $\mathbf{u}$. It was found crucial to
marginalize $\mathbf{u}$ out of the model when sampling $\alpha_{j}$ in the
MCMC.

We denote by
$\bm{\theta}^{*}=(\alpha^{*}_{j},\bm{\alpha}_{-j},\bm{\beta},\bm{\beta}_{0},\phi,\bm{\lambda})$
the non-augmented parameter vector at the proposed
$\alpha_{j}=\alpha^{*}_{j}$, with $\bm{\alpha}_{-j}$ denoting the vector
$\alpha$ excluding its $j$-th entry and by
$\bm{\theta}=(\bm{\alpha},\bm{\beta},\bm{\beta}_{0},\phi,\bm{\lambda})$ the
parameter vector with $\alpha_{j}$ representing the current iteration. Notice
that the auxiliary variables $\mathbf{u}$ are marginalized out. The proposed
$\alpha^{*}_{j}$ is accepted with probability
$\rho(\alpha^{*}_{j}\mid\alpha_{j})=\min\\{1,r(\alpha^{*}_{j}\mid\alpha_{j})\\}$
where

$\displaystyle r(\alpha^{*}_{j}\mid\alpha_{j})$
$\displaystyle=\frac{p(\mathbf{y}\mid\bm{\theta}^{*})p(\bm{\theta}^{*})q(\alpha_{j}\mid\alpha^{*}_{j})}{p(\mathbf{y}\mid\bm{\theta})p(\bm{\theta})q(\alpha^{*}_{j}\mid\alpha_{j})}$
$\displaystyle=\frac{p(\bm{\beta}_{j}\mid\lambda_{j},\phi,\alpha^{*}_{j})p(\alpha^{*}_{j})q(\alpha_{j}\mid\alpha^{*}_{j})}{p(\bm{\beta}_{j}\mid\lambda_{j},\phi,\alpha_{j})p(\alpha_{j})q(\alpha^{*}_{j}\mid\alpha_{j})}.$

We chose a prior on $\alpha_{j}$ given by $\alpha_{j}=2.5\eta_{j}$ where
$\eta_{j}\sim Beta(a_{\eta},b_{\eta})$. The MH proposal
$q(\alpha^{*}_{j}\mid\alpha_{j})$ is derived from a Gaussian random walk on
$v=\log\frac{\alpha_{j}}{2.5-\alpha_{j}}$, i.e., $q(v^{*}\mid
v)=N(v^{*}\mathchar 24635\relax\;v,\ W)$ where $N(x\mathchar 24635\relax\;\
\mu,\sigma^{2})$ denotes the Gaussian density with mean $\mu$ and variance
$\sigma^{2}$ evaluated at $x$. It is easy to show that the Gaussian random
walk $q(v^{*}\mid v)$ on $v$ implies
$q(\alpha^{*}_{j}\mid\alpha_{j})=N\left(\log\frac{\alpha^{*}_{j}}{2.5-\alpha^{*}_{j}}\mathchar
24635\relax\;\
\log\frac{\alpha_{j}}{2.5-\alpha_{j}},W_{j}\right)\times\frac{2.5}{\alpha^{*}_{j}(2.5-\alpha^{*}_{j})}.$

After simplifications,

$\displaystyle r$
$\displaystyle(\alpha^{*}_{j}\mid\alpha_{j})=\frac{{\alpha^{*}_{j}}^{K_{j}}}{\alpha^{K_{j}}_{j}}\times\frac{\Gamma(\alpha^{-1}_{j})^{K_{j}}}{\Gamma({\alpha^{*}_{j}}^{-1})^{K_{j}}}\times\lambda_{j}^{\frac{K_{j}}{\alpha^{*}_{j}}-\frac{K_{j}}{\alpha_{j}}}\times$
$\displaystyle\times\exp\left\\{-\lambda_{j}\sum^{K_{j}}_{\ell=1}\left(\mid\beta_{j\ell}\mid^{\alpha^{*}_{j}}\phi^{\frac{\alpha^{*}_{j}}{2}}-\mid\beta_{j\ell}\mid^{\alpha}\phi^{\frac{\alpha}{2}}\right)\right\\}\times$
$\displaystyle\times\frac{{\alpha^{*}_{j}}^{a_{\eta}}(2.5-\alpha^{*}_{j})^{b_{\eta}}}{\alpha^{a_{\eta}}_{j}(2.5-\alpha_{j})^{b_{\eta}}}.$

####  A.2.2 Non-marginalized MH proposal

Section A.2.1 described a MH scheme for $\alpha_{j}$ taking the advantage of
the marginalization of the auxiliary variables $\mathbf{u}.$ This section
describes a simpler alternative that does not marginalize $\mathbf{u}$.
However, we could not get well mixing Markov chains by usig the simpler
proposal described here.

We chose a prior on $\alpha_{j}$ given by $\alpha_{j}=2.5\eta_{j}$ where
$\eta_{j}\sim Beta(a_{\eta},b_{\eta})$. The MH proposal
$q(\alpha^{*}_{j}|\alpha_{j})$ is derived from a Gaussian random walk on
$v=\log\frac{\alpha_{j}}{2.5-\alpha_{j}}$, i.e., $q(v^{*}\mid
v)=N(v^{*}\mathchar 24635\relax\;v,W)$ where $N(x\mathchar
24635\relax\;\mu,\sigma^{2})$ denotes the Gaussian density with mean $\mu$ and
variance $\sigma^{2}$ evaluated at $x$. It is easy to show that the Gaussian
random walk $q(v^{*}\mid v)$ on $v$ implies
$q(\alpha^{*}_{j}\mid\alpha_{j})=N\left(\log\frac{\alpha^{*}_{j}}{2.5-\alpha^{*}_{j}}\mathchar
24635\relax\;\
\log\frac{\alpha_{j}}{2.5-\alpha_{j}},W\right)\times\frac{2.5}{\alpha^{*}_{j}(2.5-\alpha^{*}_{j})}.$

After simplifications,

$\displaystyle r(\alpha^{*}_{j}\mid\alpha_{j})$
$\displaystyle=\lambda_{j}^{\left(\frac{K_{j}}{\alpha^{*}_{j}}-\frac{K_{j}}{\alpha_{j}}\right)}\times\frac{\Gamma\left(\frac{1}{\alpha^{*}_{j}}+1\right)^{K_{j}}}{\Gamma\left(\frac{1}{\alpha_{j}}+1\right)^{K_{j}}}\times$
$\displaystyle\times\frac{{\alpha^{*}_{j}}^{a_{\eta}}(2.5-\alpha^{*}_{j})^{b_{\eta}}}{{\alpha_{j}}^{a_{\eta}}(2.5-\alpha_{j})^{b_{\eta}}}\mathds{1}(m\leq\alpha^{*}_{j}\leq
M),$

where

$m=\max\left(\\{0\\}\cup\left\\{\frac{\log(u_{jk})}{\log(\mid\beta_{jk}\mid\phi^{1/2})}:\
\ k\in S_{-}\right\\}\right),$
$M=\min\left(\\{2.5\\}\cup\left\\{\frac{\log(u_{jk})}{\log(\mid\beta_{jk}\mid\phi^{1/2})}:\
\ k\in S_{+}\right\\}\right),$ $S_{-}:=\\{j=1,\ldots,K_{j}:\
0<\mid\beta_{jk}\mid<\phi^{-1/2}\\},$ $S_{+}:=\\{j=1,\ldots,K_{j}:\
\mid\beta_{jk}\mid>\phi^{-1/2}\\}.$

It is important to notice that $S_{-}$ does not include values of $k$ such
that $\beta_{jk}=0$.

##  Appendix B Details of BBVI for semi-parametric Bayesian bridge

The details for implementation of BBVI for the proposed semi-parametric
Bayesian bridge are described in this section. Only the case with $D=1$ is
shown, although the calculations could be easily extended to the multivariate
case where $D\in\mathbb{N}$.

Under mean-field, the proposed marginal variational distributions are
$q(\bm{\beta}\mid m^{\prime}_{\beta},S_{\beta})=N(\bm{\beta}\mid
m^{\prime}_{\beta},S_{\beta}),\ q(\gamma\mid
m^{\prime}_{\gamma},S_{\gamma})=N(\bm{\gamma}\mid
m^{\prime}_{\gamma},S_{\gamma}),\ q(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi})=Ga(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi}),\ q(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda})=Ga(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda}),\ q(\alpha\mid
a^{\prime}_{\alpha},b^{\prime}_{\alpha})=2\times Beta(\alpha\mid
a^{\prime}_{\alpha},b^{\prime}_{\alpha}).$

We now describe the analytical expressions for the gradient of $\log
q_{\bm{\psi}}(\bm{\theta})$ with respect to each variational parameter. In the
equations bellow, $\mbox{dig}(\cdot)$ denotes the digamma function, i.e., the
derivative of the log gamma function.

$\displaystyle\nabla_{S^{-1}_{\beta}}\log q(\bm{\beta}\mid
m^{\prime}_{\beta},S_{\beta})$
$\displaystyle=-\frac{1}{2}(\bm{\beta}-m^{\prime}_{\beta})(\bm{\beta}-m^{\prime}_{\beta})^{\top}$
$\displaystyle{\color[rgb]{1,1,1}{.}}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ +\frac{1}{2}S_{\beta}$
$\displaystyle\nabla_{m^{\prime}_{\beta}}\log q(\bm{\beta}\mid
m^{\prime}_{\beta},S_{\beta})$
$\displaystyle=S^{-1}_{\beta}(\bm{\beta}-m^{\prime}_{\beta})$
$\displaystyle\nabla_{S^{-1}_{\gamma}}\log q(\bm{\gamma}\mid
m^{\prime}_{\gamma},S_{\gamma})$
$\displaystyle=\frac{1}{2}S_{\gamma}-\frac{1}{2}(\bm{\gamma}-m^{\prime}_{\gamma})(\bm{\gamma}-m^{\prime}_{\gamma})^{\top}$
$\displaystyle\nabla_{m^{\prime}_{\gamma}}\log q(\bm{\gamma}\mid
m^{\prime}_{\gamma},S_{\gamma})$
$\displaystyle=S^{-1}_{\gamma}(\bm{\gamma}-m^{\prime}_{\gamma})$
$\displaystyle\nabla_{a^{\prime}_{\phi}}\log q(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi})$ $\displaystyle=\log
b^{\prime}_{\phi}-\mbox{dig}(a^{\prime}_{\phi})+\log\phi$
$\displaystyle\nabla_{b^{\prime}_{\phi}}\log q(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi})$
$\displaystyle=\frac{a^{\prime}_{\phi}}{b^{\prime}_{\phi}}-\phi$
$\displaystyle\nabla_{a^{\prime}_{\lambda}}\log q(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda})$ $\displaystyle=\log
b^{\prime}_{\lambda}-\mbox{dig}(a^{\prime}_{\lambda})+\log\lambda$
$\displaystyle\nabla_{b^{\prime}_{\lambda}}\log q(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda})$
$\displaystyle=\frac{a^{\prime}_{\lambda}}{b^{\prime}_{\lambda}}-\lambda$
$\displaystyle\nabla_{a^{\prime}_{\alpha}}\log q(\alpha\mid
a^{\prime}_{\alpha},b^{\prime}_{\alpha})$ $\displaystyle=-\log
2+\mbox{dig}(a^{\prime}_{\alpha}+b^{\prime}_{\alpha})$
$\displaystyle-\mbox{dig}(a^{\prime}_{\alpha})+\log\alpha$
$\displaystyle\nabla_{b^{\prime}_{\alpha}}\log q(\alpha\mid
a^{\prime}_{\alpha},b^{\prime}_{\alpha})$ $\displaystyle=-\log
2+\mbox{dig}(a^{\prime}_{\alpha}+b^{\prime}_{\alpha})$
$\displaystyle-\mbox{dig}(b^{\prime}_{\alpha})+\log(2-\alpha)$

##  Appendix C Variational family in Bayes bridge reparameterization method

The details for implementation of the reparameterization method for the
proposed semi-parametric Bayesian bridge are described in this section. Only
the case with $D=1$ is shown, although the calculations could be easily
extended to the multivariate case where $D\in\mathbb{N}$.

For the Bayesian lasso model, we have
$\bm{\theta}=(\bm{\beta}^{\top},\bm{\gamma}^{\top},\lambda,\phi)^{\top}$ as
the parameter vector. The proposed marginal variational distributions are
$q(\bm{\beta}\mid m^{\prime}_{\beta},S_{\beta})=N(\bm{\beta}\mid
m^{\prime}_{\beta},S_{\beta}),\ q(\gamma\mid
m^{\prime}_{\gamma},S_{\gamma})=N(\bm{\gamma}\mid
m^{\prime}_{\gamma},S_{\gamma}),\ q(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi})=\log N(\phi\mid
a^{\prime}_{\phi},b^{\prime}_{\phi}),\ q(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda})=\log N(\lambda\mid
a^{\prime}_{\lambda},b^{\prime}_{\lambda}).$ The proposed joint variational
distribution follows the mean-field assumption, i.e.
$q(\bm{\theta})=q(\bm{\beta})q(\bm{\gamma})q(\lambda)q(\phi).$

The entropy of multivariate normal distributions and log-normal distributions
are available in closed form. If $q(\bm{\theta})=N_{d}(\bm{\theta}\mathchar
24635\relax\;\ \mathbf{m},\mathbf{S}),$ where $N_{d}(\bm{\theta}\mathchar
24635\relax\;\ \mathbf{m},\mathbf{S})$ denotes the density of a
$d$-dimensional multivariate normal distribution with mean vector $\mathbf{m}$
and covariance matrix $\mathbf{S}$ evaluated at $\bm{\theta},$ the entropy of
$q(\bm{\theta})$ is $H(q(\bm{\theta})\mathchar
24635\relax\;\mathbf{m},\mathbf{S})=\frac{1}{2}\log\mid 2\pi
e\mathbf{S}\mid=\frac{d}{2}[\log(2\pi)+1]+\sum^{d}_{i=1}\ell_{ii},$ where
$\ell_{ii}$ denotes the $i$-th entry of the diagonal of $\mathbf{L},$ the
Choleskey decomposition of $\mathbf{S}.$ For the log-Normal distribution, if
$q(\theta)=\log N(\theta\mathchar 24635\relax\;m,s),$ meaning $\log\theta\sim
N(m,s^{2}),$ then $H(q(\theta)\mathchar
24635\relax\;m,s)=\log_{2}\left(s\sqrt{2\pi}e^{\mu+0.5}\right).$

##  Appendix D Detailed expression for gradient of the ELBO under ADVI

This section derives in details the approximation of the gradient of the ELBO
in equation (9) to facilitate implementation of Algorithm 1.

Let $\epsilon^{(\ell)}\sim N(\bf 0_{d},\mathbf{I}_{d}),$ $\ell=1,\ldots,M$
where $d=K_{0}+K_{1}+\ldots+K_{D}+2D+1$ is the dimension of the parameter
vector $\bm{\theta}$. Then we compute
$\bm{\xi}^{(\ell)}=\mathbf{m}+\mathbf{L}\bm{\epsilon}^{(\ell)}$ and using the
transformation $T$ defined in section 4.3, we compute the implied Monte Carlo
samples in the original parameter vector $\bm{\theta}$ as

$\displaystyle\bm{\beta}_{0}^{(\ell)}$
$\displaystyle=\bm{\xi}_{\bm{\beta}_{0}^{(\ell)}},$
$\displaystyle\bm{\beta}_{d}^{(\ell)}$
$\displaystyle=\bm{\xi}_{\bm{\beta}_{d}^{(\ell)}},\ d=1,\ldots,D$
$\displaystyle\phi^{(\ell)}$ $\displaystyle=e^{\xi_{{\phi}^{(\ell)}}},$
$\displaystyle\lambda_{d}^{(\ell)}$
$\displaystyle=e^{\xi_{\lambda_{d}^{(\ell)}}},\ d=1,\ldots,D$
$\displaystyle\alpha_{d}^{(\ell)}$
$\displaystyle=\frac{2.5}{1+e^{-\xi_{{\alpha_{d}^{(\ell)}}}}},\ d=1,\ldots,D.$

Starting from equation (9), we get

$\widetilde{\nabla}_{\bm{\psi}}ELBO(\tilde{\mathbf{y}},\bm{\psi})=\\\
=\frac{n}{KM}\sum^{M}_{\ell=1}\left[\nabla_{\bm{\psi}}\log
p(\tilde{\mathbf{y}}\mid T^{-1}(\bm{\xi}^{(\ell)}))\right]\ +\ \\\
+\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\left[\log
p(T^{-1}(\bm{\xi}^{(\ell)}))+\log\mid
J_{T^{-1}(\bm{\xi}^{(\ell)})}\mid\right]\\\
+\nabla_{\bm{\psi}}\mathbb{E}_{\bm{\xi}\sim
N(\mathbf{m},\mathbf{L}\mathbf{L}^{T})}\left[-\log N(\bm{\xi}\mathchar
24635\relax\;\ \mathbf{m},\mathbf{L}\mathbf{L}^{\top})\right]\\\ =\\\
\frac{n}{KM}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\log
N\left(\widetilde{\mathbf{y}}\mathchar 24635\relax\;\
\sum^{D}_{j=0}\widetilde{X}_{j}\bm{\beta}^{(\ell)}_{j},{\phi^{(\ell)}}^{-1}\mathbf{I}\right)+\\\
\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\left[\log
p(\bm{\beta}^{(\ell)}_{0})+\log p(\phi^{(\ell)})+\right.\\\
\left.+\sum^{D}_{j=1}\sum^{K_{j}}_{k=1}\log
p(\bm{\beta}^{(\ell)}_{jk}|\lambda^{(\ell)}_{j},\phi^{(\ell)},\alpha^{(\ell)}_{j})+\right.\\\
+\left.\sum^{D}_{j=1}\log p(\lambda^{(\ell)}_{j})+\sum^{D}_{j=1}\log
p(\alpha^{(\ell)}_{j})\right]+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\bm{\xi}_{\phi}^{(\ell)}+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\nabla_{\bm{\psi}}\left[\xi^{(\ell)}_{\lambda_{j}}+\xi^{(\ell)}_{\alpha_{j}}-2\log\left(1+e^{\xi^{(\ell)}_{\alpha_{j}}}\right)\right]+\\\
+\nabla_{\bm{\psi}}\sum^{d}_{i=1}\mathbf{L}_{ii}\\\ =\\\
\frac{n}{KM}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\log
N\left(\widetilde{\mathbf{y}}\mathchar 24635\relax\;\
\sum^{D}_{j=0}\widetilde{X}_{j}\bm{\beta}^{(\ell)}_{j},{\phi^{(\ell)}}^{-1}\mathbf{I}\right)+\\\
\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\log
N(\bm{\beta}^{(\ell)}_{0}\mathchar
24635\relax\;\bm{\mu}_{0},\bm{\Sigma}_{0})+\\\
\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\log Ga(\phi^{(\ell)}\mathchar
24635\relax\;a_{\phi},b_{\phi})+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\sum^{K_{j}}_{k=1}\nabla_{\bm{\psi}}\log
GG\left(\bm{\beta}^{(\ell)}_{jk}\mathchar
24635\relax\;0,\frac{{\lambda_{j}^{(\ell)}}^{-\frac{1}{\alpha^{(\ell)}_{j}}}}{\sqrt{\phi^{(\ell)}}},\alpha^{(\ell)}_{j}\right)+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\nabla_{\bm{\psi}}\log
Ga(\lambda^{(\ell)}_{j}\mathchar 24635\relax\;a_{\lambda},b_{\lambda})+\\\
\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\nabla_{\bm{\psi}}\log{\alpha^{(\ell)}_{j}}^{a_{\eta}-1}(2.5-\alpha^{(\ell)}_{j})^{b_{\eta}-1}+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\nabla_{\bm{\psi}}\log\phi^{(\ell)}+\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\nabla_{\bm{\psi}}\log\lambda^{(\ell)}_{j}+\\\
+\frac{1}{M}\sum^{M}_{\ell=1}\sum^{D}_{j=1}\nabla_{\bm{\psi}}\left(\log\frac{\alpha^{(\ell)}_{j}}{2.5-\alpha^{(\ell)}_{j}}-2\log\frac{2.5}{2.5-\alpha^{(\ell)}_{j}}\right)+\\\
+\nabla_{\bm{\psi}}\sum^{d}_{i=1}\mathbf{L}_{ii},$ (11)

where the term $\nabla_{\bm{\psi}}\sum^{d}_{i=1}\mathbf{L}_{ii}$ comes from
the fact that the entropy of a multivariate Gaussian random vector
$\mathbf{x}\sim N(\bm{\mu},\bm{\Sigma})$ is
$\mathbb{H}[\mathbf{x}]=\mathbb{E}_{\mathbf{x}\sim
N(\bm{\mu},\bm{\Sigma})}[-\log
N(\mathbf{x}\mid\bm{\mu},\bm{\Sigma})]=\log|2\pi e\bm{\Sigma}|$.

##  Appendix E Further results and comparisons

This section presents more results on the simulation described in section 5.1.

Figure 14 contains the posterior marginal distributions for all parameters of
the model described in section 5.1. The marginals for
$\beta_{1},\ldots,\beta_{7}$ are also shown in Figure 4. As in Figure 4,
Figure 14 shows high level of agreement between marginals obtained via MCMC
and ADVI, with brms having higher discrepancy with respect to the MCMC. INLA
also approximates well the marginal distributions except for a small number of
coefficients (e.g., $\beta_{1},\beta_{2},\beta_{34}$)

Figure 15 complements Figure 8 from section 5.1 with the point estimates for
all parameters obtained via ADVI and MCMC for each one of the 100 simulated
datasets.

Figure 14: Marginal posterior distributions for $\beta_{k},\ k=1,\ \ldots,34$
based on the first replicate of the simulated data. Red dots denote the true
values of the parameters. Figure 15: Point estimates for all parameters under
MCMC and VB for each one of the 100 replicas. The black lines represent the
identity function.

To further address the quality of the variational approximations to the true
posterior distribution, Figure 16 shows the p-values obtained when comparing
(via Kolmogorov-Smirnov tests) marginal posterior predictive distributions
estimated by VB and MCMC. The shapes of the histograms are approximately
uniform with a low prevalence of small p-values, as expected under the null
hypothesis that the distributions under MCMC and VB are the same. For example,
the empirical proportions of p-values below 0.05 are close to 0.05, which
corresponds to the expected proportion of false discoveries when $H_{0}$ is
true and a type 1 error of 0.05 is fixed.

Figure 16: p-values for the Kolmogorov-Smirnov test for comparison of each
marginal posterior predictive distributions under ADVI and MCMC. In red, we
highlight P(type 1 error) = 5%. The expression P(Rej. $H_{0}$) represents the
proportion of times we would reject the null hypothesis of equal distributions
under ADVI and MCMC if we were to reject whenever p-value $<$ 0.05.

## References

  *   * Alves . (2021) alves2021Alves, L., Dias, R. Migon, H.S.  2021\.  Variational Full Bayes Lasso: Knots Selection in Regression Splines Variational full bayes lasso: Knots selection in regression splines.  arXiv preprint arXiv:2102.13548.
  * Armagan (2009) armagan2009Armagan, A.  2009\.  Variational bridge regression Variational bridge regression.  Artificial Intelligence and Statistics Artificial intelligence and statistics ( 17–24).
  * Blei . (2017) blei2017Blei, D.M., Kucukelbir, A. McAuliffe, J.D.  2017\.  Variational inference: A review for statisticians Variational inference: A review for statisticians.  Journal of the American statistical Association112518859–877.
  * Bürkner (2017) brmsBürkner, PC.  2017\.  brms: An R Package for Bayesian Multilevel Models Using Stan brms: An R package for Bayesian multilevel models using Stan.  Journal of Statistical Software8011–28.
  * Carvalho . (2009) carvalho2009Carvalho, C.M., Polson, N.G. Scott, J.G.  2009\.  Handling sparsity via the horseshoe Handling sparsity via the horseshoe.  Artificial Intelligence and Statistics Artificial intelligence and statistics ( 73–80).
  * Casella . (2010) Casella_2010Casella, G., Ghosh, M., Gill, J. Kyung, M.  2010\.  Penalized regression, standard errors, and Bayesian lassos Penalized regression, standard errors, and Bayesian lassos.  Bayesian Analysis52369 – 411.  https://doi.org/10.1214/10-BA607  10.1214/10-BA607
  * Currie Durban (2002) currie2002Currie, I.D. Durban, M.  2002\.  Flexible smoothing with P-splines: a unified approach Flexible smoothing with p-splines: a unified approach.  Statistical Modelling24333–349.
  * Damien Walker (2001) Damien_Walker_2012Damien, P. Walker, S.G.  2001\.  Sampling Truncated Normal, Beta, and Gamma Densities Sampling truncated normal, beta, and gamma densities.  Journal of Computational and Graphical Statistics102206-215.
  * Denison . (1998) denison1998Denison, D., Mallick, B. Smith, A.  1998\.  Automatic Bayesian curve fitting Automatic Bayesian curve fitting.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)602333–350.
  * Dias Gamerman (2002) dias2002Dias, R. Gamerman, D.  2002\.  A Bayesian approach to hybrid splines non-parametric regression A Bayesian approach to hybrid splines non-parametric regression.  Journal of Statistical Computation and Simulation724285–297.
  * Duchi . (2011) duchi2011Duchi, J., Hazan, E. Singer, Y.  2011\.  Adaptive subgradient methods for online learning and stochastic optimization. Adaptive subgradient methods for online learning and stochastic optimization.  Journal of machine learning research127.
  * Eilers Marx (1996) eilers1996Eilers, P.H. Marx, B.D.  1996\.  Flexible smoothing with B-splines and penalties Flexible smoothing with b-splines and penalties.  Statistical science11289–121.
  * Frank Friedman (1993) frank1993Frank, L.E. Friedman, J.H.  1993\.  A statistical view of some chemometrics regression tools A statistical view of some chemometrics regression tools.  Technometrics352109–135.
  * Hastie Tibshirani (1986) hastie1986Hastie, T. Tibshirani, R.  1986\.  Generalized Additive Models Generalized additive models.  Statistical Science13297–310.  [2022-08-19]http://www.jstor.org/stable/2245459
  * Hastie Tibshirani (2000) hastie2000Hastie, T. Tibshirani, R.  2000\.  Bayesian Backfitting Bayesian backfitting.  Statistical Science153196–213.  [2022-06-23]http://www.jstor.org/stable/2676659
  * Hastie . (2015) hastie2015Hastie, T., Tibshirani, R. Wainwright, M.  2015\.  Statistical learning with sparsity Statistical learning with sparsity.  Monographs on statistics and applied probability143143.
  * Hoerl Kennard (19701) hoerl1970bHoerl, A.E. Kennard, R.W.  19701\.  Ridge regression: applications to nonorthogonal problems Ridge regression: applications to nonorthogonal problems.  Technometrics12169–82.
  * Hoerl Kennard (19702) hoerl1970aHoerl, A.E. Kennard, R.W.  19702\.  Ridge regression: Biased estimation for nonorthogonal problems Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics12155–67.
  * Kingma Ba (2014) adamKingma, D.P. Ba, J.  2014\.  Adam: A method for stochastic optimization Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kingma Welling (2013) kingma2013Kingma, D.P. Welling, M.  2013\.  Auto-encoding variational Bayes Auto-encoding variational Bayes.  arXiv preprint arXiv:1312.6114.
  * Kucukelbir . (2015) kucukelbir2015Kucukelbir, A., Ranganath, R., Gelman, A. Blei, D.  2015\.  Automatic variational inference in Stan Automatic variational inference in stan.  Advances in neural information processing systems28.
  * Kucukelbir . (2017) kucukelbir2017Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A. Blei, D.M.  2017\.  Automatic differentiation variational inference Automatic differentiation variational inference.  The Journal of Machine Learning Research181430–474.
  * Lang Brezger (2004) lang2004Lang, S. Brezger, A.  2004\.  Bayesian P-splines Bayesian p-splines.  Journal of computational and graphical statistics131183–212.
  * Leng . (2014) Leng_2014Leng, C., Tran, MN. Nott, D.  2014\.  Bayesian adaptive lasso Bayesian adaptive lasso.  Annals of the Institute of Statistical Mathematics66221–244.
  * F. Li Villani (2013) Li2013Li, F. Villani, M.  2013\.  Efficient Bayesian Multivariate Surface Regression Efficient bayesian multivariate surface regression.  Scandinavian Journal of Statistics404706-723.
  * Q. Li Lin (2010) li2010Li, Q. Lin, N.  2010\.  The Bayesian elastic net The Bayesian elastic net.  Bayesian analysis51151–170.
  * Luts Wand (2015) luts2015Luts, J. Wand, M.P.  2015\.  Variational inference for count response semiparametric regression Variational inference for count response semiparametric regression.
  * Mallick Yi (2018) mallick2018Mallick, H. Yi, N.  2018\.  Bayesian bridge regression Bayesian bridge regression.  Journal of Applied Statistics456988–1008.
  * Menictas Wand (2015) Menictas2015Menictas, M. Wand, M.P.  2015\.  Variational Inference for Heteroscedastic Semiparametric Regression Variational inference for heteroscedastic semiparametric regression.  Australian & New Zealand Journal of Statistics571119-138.
  * Ong . (2017) Ong2017Ong, V.M.H., Mensah, D.K., Nott, D.J., Jo, S., Park, B. Choi, T.  2017\.  A variational Bayes approach to a semiparametric regression using Gaussian process priors A variational Bayes approach to a semiparametric regression using Gaussian process priors.  Electronic Journal of Statistics1124258 – 4296.
  * Paisley . (2012) paisley2012Paisley, J.W., Blei, D.M. Jordan, M.I.  2012\.  Variational Bayesian Inference with Stochastic Search Variational Bayesian inference with stochastic search.  Proceedings of the 29th International Conference on Machine Learning (ICML-12) Proceedings of the 29th international conference on machine learning (icml-12) ( 1367–1374).
  * Park Casella (2008) park2008Park, T. Casella, G.  2008\.  The Bayesian lasso The Bayesian lasso.  Journal of the American Statistical Association103482681–686.
  * Polson . (2014) polson2014Polson, N.G., Scott, J.G. Windle, J.  2014\.  The Bayesian bridge The Bayesian bridge.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)764713–733.
  * Ranganath . (2014) ranganath2014Ranganath, R., Gerrish, S. Blei, D.M.  2014\.  Black Box Variational Inference. Black box variational inference.  Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics. Proceedings of the seventeenth international conference on artificial intelligence and statistics.
  * Rue . (2009) rue2009Rue, H., Martino, S. Chopin, N.  2009\.  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  Journal of the royal statistical society: Series b (statistical methodology)712319–392.
  * Silverman Green (1994) silv:gree:1994Silverman, B.W. Green, P.J.  1994\.  Nonparametric Regression and Generalized Linear Models Nonparametric regression and generalized linear models.  Chapman and Hall (London).
  * Subbotin (1923) subbotin1923Subbotin, M.T.  1923\.  On the law of frequency of error On the law of frequency of error.  Sbornik: Mathematics312296–301.
  * Tibshirani (1996) tibshirani1996Tibshirani, R.  1996\.  Regression shrinkage and selection via the lasso Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological)581267–288.
  * Tibshirani (1997) tibshirani1997Tibshirani, R.  1997\.  The lasso method for variable selection in the Cox model The lasso method for variable selection in the cox model.  Statistics in medicine164385–395.
  * Wand (2017) Wand2017Wand, M.P.  2017\.  Fast Approximate Inference for Arbitrarily Large Semiparametric Regression Models via Message Passing Fast approximate inference for arbitrarily large semiparametric regression models via message passing.  Journal of the American Statistical Association112517137-168.
  * West Harrison (1997) west2006West, M. Harrison, J.  1997\.  Bayesian forecasting and dynamic models Bayesian forecasting and dynamic models.  Springer Series in Statistics, New York:Springer-Verlag.
  * S. Wood (2012) wood2012Wood, S.  2012\.  mgcv: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness estimation mgcv: Mixed gam computation vehicle with gcv/aic/reml smoothness estimation.  R package version 1.8-1.
  * S.N. Wood (2017) wood2017Wood, S.N.  2017\.  Generalized additive models: an introduction with R Generalized additive models: an introduction with R (2nd ).  Boca Raton, FL: CRC Press.
  * Yin Zhou (2018) yin2018Yin, M. Zhou, M.  2018\.  Semi-implicit variational inference Semi-implicit variational inference.  International Conference on Machine Learning International conference on machine learning ( 5660–5669).
  * Zeiler (2012) zeiler2012Zeiler, M.D.  2012\.  Adadelta: an adaptive learning rate method Adadelta: an adaptive learning rate method.  arXiv preprint arXiv:1212.5701.
  * Zou Hastie (2005) zou2005Zou, H. Hastie, T.  2005\.  Regularization and variable selection via the elastic net Regularization and variable selection via the elastic net.  Journal of the royal statistical society: series B (statistical methodology)672301–320.