# Geometric ergodicity of Gibbs samplers for Bayesian error-in-variable
regression

Austin<EMAIL_ADDRESS>[ Department of Statistics,
University of Warwick, Coventry, UK

(4 2023)

###### Abstract

Multivariate Bayesian error-in-variable (EIV) linear regression is considered
to account for additional additive Gaussian error in the features and
response. A 3-variable deterministic scan Gibbs samplers is constructed for
multivariate EIV regression models using classical and Berkson errors with
independent normal and inverse-Wishart priors. These Gibbs samplers are proven
to always be geometrically ergodic which ensures a central limit theorem for
many time averages from the Markov chains. We demonstrate the strengths and
limitations of the Gibbs sampler with simulated data for large data problems,
robustness to misspecification and also analyze a real-data example in
astrophysics.

60J22,

62F15,

Bayesian statistics,

Convergence analysis,

Error-in-variable models,

Gibbs sampling,

Markov chain Monte Carlo,

Measurement error models,

###### keywords:

[class=MSC]

###### keywords:

###### Contents

  1. 1 Introduction
  2. 2 General Gibbs Sampler for EIV regression
    1. 2.1 Bayesian EIV regression with errors in the features
    2. 2.2 Bayesian EIV regression with errors in the response and features
  3. 3 Simulations
    1. 3.1 Limitations of the Gibbs sampler in large problem sizes
    2. 3.2 Robustness to model misspecification
  4. 4 Real-data example: measurement error in astrophysics
  5. 5 Conclusion and Future Directions

##  1 Introduction

Many problems in astrophysics (Feigelson and Babu, 1992; Hilbe, de Souza and
Ishida, 2017; Kelly, 2012; Stefanski, 2000) and epidemiology (Achic et al.,
2018; Buonaccorsi, 2010; Carroll et al., 2006; Clayton et al., 1992) among
other areas of science (Groß, 2016; Pollice et al., 2019; Tang, Li and Tang,
2017) involve error in variables (EIV) which classical linear regression does
not take into account. EIV can occur in many situations such as measurement
error in data collection (Hilbe, de Souza and Ishida, 2017; Kelly, 2012),
discrepancies between the data distribution and the model (Carroll et al.,
2006; Buonaccorsi, 2010), or purposeful adversarial attacks against the data
(Goodfellow, Shlens and Szegedy, 2015; Szegedy et al., 2014). Not
surprisingly, multiple critical issues arise in parameter estimation and
statistical inference when ignoring additional errors in the data such as poor
predictive performance (Goodfellow, Shlens and Szegedy, 2015), statistical
bias (Damgaard, 2020; Kröger, Hoffmann and Pakpahan, 2016; Vidal and Iglesias,
2008), and estimators fail to be consistent (Michalek and Tripathi, 1980).

Bayesian approaches develop a strategy for additional error in the variables
by constructing a new model incorporating additional error. We consider
multivariate Bayesian EIV linear regression (Charisse Farr et al., 2020;
Dellaportas and Stephens, 1995; Fang et al., 2017; Huang, 2010; Mallick and
Gelfand, 1996; Muff et al., 2015; Richardson S, 1993; Rodrigues and Bolfarine,
2007; Torabi et al., 2021; Vidal and Arellano-Valle, 2010) accounting for
additive Gaussian error in the features (covariates) and response. We assume
the variability of the additive Gaussian error is known beforehand which
arises often in astrophysics in the presence of known instrumentation error
(Hilbe, de Souza and Ishida, 2017; Kelly, 2012). Alternative approaches to EIV
models attempt to correct existing parameter estimation methods such as least
squares or method of moments with weighting and other techniques (Fuller,
1987; Stefanski and Carroll, 1985). Several other strategies for EIV modeling
are discussed in more comprehensive treatments on the topic (Buonaccorsi,
2010; Carroll et al., 2006; Fuller, 1987).

We write $x\sim N_{d}(m,C)$ to mean the $d$-dimensional normal distribution
with mean $m$ and symmetric, positive-definite (SPD) covariance matrix $C$. We
also write $x\sim\mathcal{W}_{d}^{-1}(\nu,V)$ to be the inverse-Wishart
distribution with positive integer degrees of freedom $\nu\geq d$ and scale
SPD matrix $V\in\mathbb{R}^{d\times d}$. Let $\text{vec}(A)$ denote the
vectorization of a matrix $A$ by stacking the columns. Let
$(Y_{i},X_{i},Z_{i})_{i=1}^{n}$ be independent and identically distributed
(i.i.d.) where the response $Y_{i}$ takes values in $\mathbb{R}^{m}$ along
with features $X_{i}$ taking values in $\mathbb{R}^{p}$ and fixed, known
features $Z_{i}\in\mathbb{R}^{q}$ where $m,n,p,q$ are positive integers. Let
$\theta=\text{vec}(\Theta)\in\mathbb{R}^{qm}$,
$\beta=\text{vec}(\mathcal{B})\in\mathbb{R}^{pm}$, and SPD matrix
$\Sigma\in\mathbb{R}^{m\times m}$ be unknown regression and covariance
parameters respectively. We introduce new parameters
$\mathcal{A}=(\mathcal{A}_{1},\ldots,\mathcal{A}_{n})^{T}$ with
$\mathcal{A}_{i}\in\mathbb{R}^{p}$ to model additional error in $X_{i}$ using
classical or Berkson errors (Berkson, 1950). The classical error model
specifies $X_{i}|\mathcal{A}_{i}$ and the Berkson error model (Berkson, 1950)
assumes instead a data-dependent prior on $\mathcal{A}_{i}|X_{i}$. When there
is additional error in $X_{i}$, the EIV linear regression model for $i\in
1,\ldots,n$ is i.i.d. with

$\displaystyle Y_{i}|\mathcal{A}_{i},\theta,\beta,\Sigma\sim
N_{m}(\Theta^{T}Z_{i}+\mathcal{B}^{T}\mathcal{A}_{i},\Sigma)$ (1a)
$\displaystyle X_{i}|\mathcal{A}_{i}\sim N_{p}(\mathcal{A}_{i},V_{i})\text{
(Classical)}\text{\hskip 11.38092pt or \hskip
11.38092pt}\mathcal{A}_{i}|X_{i}\sim N_{p}(X_{i},V_{i})\text{ (Berkson)}$ (1b)

where the SPD matrices $V_{i}\in\mathbb{R}^{p\times p}$ are known. When there
is also additional error in the responses $Y_{i}$, we assume an i.i.d.
hierarchical regression model with

$\displaystyle Y_{i}|\mathcal{V}_{i}\sim N_{m}(\mathcal{V}_{i},U_{i})$ (2a)
$\displaystyle\mathcal{V}_{i}|\mathcal{A}_{i},\theta,\beta,\Sigma\sim
N_{m}(\Theta^{T}Z_{i}+\mathcal{B}^{T}\mathcal{A}_{i},\Sigma)$ (2b)
$\displaystyle X_{i}|\mathcal{A}_{i}\sim N_{p}(\mathcal{A}_{i},V_{i})\text{
(Classical)}\text{\hskip 11.38092pt or \hskip
11.38092pt}\mathcal{A}_{i}|X_{i}\sim N_{p}(X_{i},V_{i})\text{ (Berkson)}$ (2c)

where $U_{i}\in\mathbb{R}^{m\times m}$ are known SPD matrices.

We will be interested in the posterior for both models (1) and (2) using
independent normal and inverse-Wishart priors on the parameters
$(\mathcal{A},\theta,\beta,\Sigma)$. The independent prior choice is a popular
choice in Bayesian regression models with and without measurement error
(Carroll et al., 2006; Dellaportas and Stephens, 1995; Ekvall and Jones, 2021;
Rajaratnam and Sparks, 2015). For the EIV regression models (1) and (2), the
independent priors are chosen

$\displaystyle\Sigma\sim\mathcal{W}^{-1}(a_{0},B_{0})$ (3a)
$\displaystyle(\theta,\beta)^{T}\sim N_{m(q+p)}(j_{0},J_{0})$ (3b)

where $a_{0}\geq m$ is a positive integer, $B_{0}\in\mathbb{R}^{m\times m}$ is
a SPD matrix, $j_{0}\in\mathbb{R}^{m(q+p)}$ and SPD matrix
$J_{0}\in\mathbb{R}^{m(q+p)\times m(q+p)}$. The classical and Berkson error
models assume either

$\displaystyle\mathcal{A}_{i}\sim N_{p}(k_{i},K_{i})\text{
(Classical)}\text{\hskip 11.38092pt or \hskip 11.38092pt}\mathcal{A}_{i}\text{
flat prior}\text{ (Berkson)}$ (4)

where $k_{i}\in\mathbb{R}^{p}$ and $K_{i}\in\mathbb{R}^{p\times p}$ are SPD
matrices. For example, an exposure model (Gustafson, 2003) utilized often in
epidemiology would assume classical errors with a data-dependent prior on each
$\mathcal{A}_{i}$ depending on $Z_{i}$. In the Berkson error model, each
$\mathcal{A}_{i}|X_{i}$ is already specified and it is natural to assume an
improper flat prior on each $\mathcal{A}_{i}$.

Previous work has proposed Gibbs sampling (Geman and Geman, 1984) to draw
samples from the posterior, denoted by $\Pi_{n}$, in Bayesian EIV regression
models (Bhadra and Carroll, 2016; Carroll et al., 2006; Dellaportas and
Stephens, 1995; Richardson S, 1993). However, trustworthy estimation from a
Gibbs sampler requires the Markov chain to converge to the posterior
distribution at a sufficiently fast rate. Consider a vector-valued function
$f$ with $\int\left\lVert f\right\rVert^{2+\delta}d\Pi_{n}<\infty$ for some
$\delta\in(0,\infty)$ and denote $\bar{f}_{m}$ as the time average of $m$
samples from the Gibbs sampler. In order to be confident in the estimator
$\bar{f}_{m}$ in applications, a standard error and confidence interval are
essential. A Gibbs sampler is geometrically ergodic if initialized at points,
its marginal distribution is converging to $\Pi_{n}$ at an exponential rate in
total variation. Geometrically ergodic Gibbs samplers provide rich theoretical
guarantees which are of practical relevance in applications. These Gibbs
samplers satisfy a central limit theorem (Chan and Geyer, 1994; Jones, 2004),
that is,

$\sqrt{m}\left(\bar{f}_{m}-\int fd\Pi_{n}\right)$

is asymptotically normally distributed and under suitable assumptions, the
covariance in this normal distribution can be consistently estimated (Flegal
and Jones, 2010). Further pertinent tools to ensuring reliable estimation such
as estimates of the effective sample size, consistent confidence ellipsoids,
and consistent confidence intervals for quantile estimation are also available
(Doss et al., 2014; Vats, Flegal and Jones, 2019a).

To the best of our knowledge, the rate of convergence for Gibbs sampling in
EIV regression models has not been previously investigated. Related approaches
have instead proposed variational Bayesian methods (Bresson et al., 2021;
Pham, Ormerod and Wand, 2013) and the integrated nested Laplace approximation
(INLA) (Håvard Rue, 2009; Muff et al., 2015). We construct a general density
which in special cases, is the posterior for the $4$ Bayesian EIV regression
models (1) and (2) using the independent normal and inverse-Wishart prior
choice on the parameters (3) and (4). Our main contribution constructs a
3-variable deterministic scan Gibbs sampler for this general density, and we
show it is always geometrically ergodic using a drift and minorization
condition (Hairer and Mattingly, 2011; Meyn and Tweedie, 2009). The 3-variable
Gibbs sampler we construct can be simulated efficiently on a computer without
the need for complex Metropolis-Hastings or rejection sampling steps at each
iteration.

The organization of this paper is as follows. In Section 2, we construct a
general EIV regression density and construct a 3-variable Gibbs sampler for
this density. We show the Gibbs sampler is always geometrically ergodic and
apply this to the $4$ multivariate Bayesian EIV regression models presented in
this introduction. Section 3 studies the algorithm empirically where we
demonstrate limitations of the Gibbs sampler with simulated data for large
data problems and also the behavior of the Gibbs sampler under model
misspecification. Section 4 studies a real-data example in astrophysics to
study supermassive black hole mass (Harris, Poole and Harris, 2014; Hilbe, de
Souza and Ishida, 2017). Finally in Section 5, we discuss our results and
future research directions.

##  2 General Gibbs Sampler for EIV regression

For positive integers $p$, define $p$-norms by
$\left\lVert\cdot\right\rVert_{p}$ and the Frobenius norm by
$\left\lVert\cdot\right\rVert_{F}$. Let $\otimes$ denote the Kronecker
product. The posteriors for the Bayesian EIV regression models (1) and (2)
using independent prior choices (3) and (4) for both classical and Berkson
errors share a common general form which we study in this section. The
posterior densities for these Bayesian EIV regression models are special cases
of the density (5) but will differ depending on the EIV modeling choice
illustrated in the subsequent sections. For $i\in 1,\ldots,n$, define hyper-
parameters $a_{0}\in(0,\infty)$, SPD matrix $B_{0}\in\mathbb{R}^{m\times m}$,
$c_{0}\in\mathbb{R}^{m(p+q)}$, SPD matrices $C_{0}\in\mathbb{R}^{m(p+q)\times
m(p+q)}$, $D_{i}\in\mathbb{R}^{p\times p}$, and $d_{i}\in\mathbb{R}^{p}$,
$R=(R_{1},\ldots,R_{n})^{T}\in\mathbb{R}^{n\times m}$, and
$M\in\mathbb{R}^{n\times q}$. For $\mathcal{A}\in\mathbb{R}^{n\times p}$,
$\theta=\text{vec}(\Theta)\in\mathbb{R}^{mq}$,
$\beta=\text{vec}(\mathcal{B})\in\mathbb{R}^{mp}$, SPD matrix
$\Sigma\in\mathbb{R}^{m\times m}$, define the density

$\displaystyle\pi_{n}(\mathcal{A},\theta,\beta,\Sigma)$ (5a)
$\displaystyle\propto\left(\frac{1}{\det(\Sigma)}\right)^{(n+a_{0}+1+m)/2}\exp\left[-\frac{1}{2}\text{tr}(\Sigma^{-1}B_{0})\right]$
(5b) $\displaystyle\hskip
11.38092pt\times\exp\left[-\frac{1}{2}\sum_{i=1}^{n}(R_{i}-\Theta^{T}M_{i}-\mathcal{B}^{T}\mathcal{A}_{i})^{T}\Sigma^{-1}(R_{i}-\Theta^{T}M_{i}-\mathcal{B}^{T}\mathcal{A}_{i})\right]$
(5c) $\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-d_{i})^{T}D_{i}^{-1}(\mathcal{A}_{i}-d_{i})\right)$
(5d) $\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}((\theta,\beta)^{T}-c_{0})^{T}C_{0}^{-1}((\theta,\beta)^{T}-c_{0})\right).$
(5e)

We will construct a 3-variable deterministic scan Gibbs sampler using a
specific update order for the density (5). We also derive the conditional
densities for the Gibbs sampler which can be sampled directly. Initialize
$(\theta_{0},\beta_{0},\Sigma_{0})$ and
$\mathcal{A}_{0}=(\mathcal{A}_{1,0},\ldots,\mathcal{A}_{n,0})$ from an initial
distribution. For $t\in 1,\ldots$, first generate

$\displaystyle\Sigma_{t}|\mathcal{A}_{t-1},\theta_{t-1},\beta_{t-1}$
$\displaystyle\sim\mathcal{W}^{-1}\left(n+a_{0},(R-M\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})^{T}(R-M\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})+B_{0}\right)$

Next, generate $(\theta_{t},\beta_{t})^{T}|\mathcal{A}_{t-1},\Sigma_{t}\sim
N_{m(p+q)}(c_{n}(\mathcal{A}_{t-1},\Sigma_{t}),C_{n}(\mathcal{A}_{t-1},\Sigma_{t}))$
where

$\displaystyle
C_{n}(\mathcal{A}_{t-1},\Sigma_{t})=\left(\Sigma^{-1}_{t}\otimes\begin{pmatrix}M&\mathcal{A}_{t-1}\end{pmatrix}^{T}\begin{pmatrix}M&\mathcal{A}_{t-1}\end{pmatrix}+C_{0}^{-1}\right)^{-1}$
$\displaystyle
c_{n}(\mathcal{A}_{t-1},\Sigma_{t})=C_{n}(\mathcal{A}_{t-1},\Sigma_{t})\left(\left[\Sigma^{-1}_{t}\otimes\begin{pmatrix}M&\mathcal{A}_{t-1}\end{pmatrix}^{T}\right]\text{vec}(R)+C_{0}^{-1}c_{0}\right).$

Finally, generate independently

$\mathcal{A}_{i,t}|\theta_{t},\beta_{t},\Sigma_{t}\sim
N_{p}(d_{n,i}(\theta_{t},\beta_{t},\Sigma_{t}),D_{n,i}(\theta_{t},\beta_{t},\Sigma_{t}))$

where

$\displaystyle
D_{n,i}(\beta_{t},\Sigma_{t})=\left(\mathcal{B}_{t}\Sigma^{-1}_{t}\mathcal{B}_{t}^{T}+D_{i}^{-1}\right)^{-1}$
$\displaystyle
d_{n,i}(\theta_{t},\beta_{t},\Sigma_{t})=D_{n,i}(\beta_{t},\Sigma_{t})\left[D_{i}^{-1}d_{i}+\mathcal{B}_{t}\Sigma_{t}^{-1}\left(R_{i}-\Theta_{t}^{T}M_{i}\right)\right]$

to obtain $\mathcal{A}_{t}=(\mathcal{A}_{1,t},\ldots,\mathcal{A}_{n,t})^{T}$.

For points $(\mathcal{A},\theta,\beta,\Sigma)$ and
$(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})$, the
Gibbs sampler has Markov transition density

$\displaystyle
p\left((\mathcal{A},\theta,\beta,\Sigma),(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})\right)$
$\displaystyle=\pi_{n}(\mathcal{A}^{\prime}|\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})\pi_{n}(\theta^{\prime},\beta^{\prime}|\mathcal{A},{\Sigma}^{\prime})\pi_{n}({\Sigma}^{\prime}|\mathcal{A},\theta,\beta)$

and Markov transition kernel defined for suitable sets $B$ by

$\displaystyle
P\left((\mathcal{A},\theta,\beta,\Sigma),B\right)=\int\int\int_{B}p\left((\mathcal{A},\theta,\beta,\Sigma),(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})\right)d{\mathcal{A}}^{\prime}d{\theta^{\prime}}d{\beta}^{\prime}d{\Sigma}^{\prime}.$

The Markov kernel at larger iteration times $t\geq 2$ is defined recursively
with $P^{1}\equiv P$ by

$P^{t}\left((\mathcal{A},\theta,\beta,\Sigma),B\right)=\int
P^{t-1}\left(\cdot,B\right)dP\left((\mathcal{A},\theta,\beta,\Sigma),\cdot\right).$

We will use the following drift function defined by

$V(\mathcal{A},\theta,\beta)=\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-d_{i})^{T}D_{i}^{-1}(\mathcal{A}_{i}-d_{i})+\frac{1}{2}(\theta,\beta)C_{0}^{-1}(\theta,\beta)^{T}$

combined with a minorization condition to show there is a $\rho\in(0,1)$ and
$M_{0}\in(0,\infty)$ so that for any initialization
$\mathcal{A},\theta,\beta,\Sigma$,

$\displaystyle\sup_{|\varphi|\leq 1+M_{0}V}\left|\int\varphi
dP^{t}\left((\mathcal{A},\theta,\beta,\Sigma),\cdot\right)-\int\varphi
d\Pi_{n}\right|\leq M(\mathcal{A},\theta,\beta)\rho^{t}$ (6)

where
$M(\mathcal{A},\theta,\beta)=2+M_{0}V(\mathcal{A},\theta,\beta)+M_{0}\int
Vd\Pi_{n}$ (Hairer and Mattingly, 2011). The condition (6) implies the Gibbs
sampler is geometrically ergodic. We now state our main result.

######  Theorem 1.

The 3-variable deterministic scan Gibbs sampler
$(\mathcal{A}_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for the
general density (5) is geometrically ergodic.

###### Proof.

Using a special property of the Gibbs sampler, it will be sufficient to
develop a drift and minorization condition based only on the marginal chain
$(\mathcal{A}_{t},\theta_{t},\beta_{t})_{t}$ (Roberts and Rosenthal, 2001,
Example 3.6). In particular, we will use the special property of this Gibbs
Markov kernel $P$ that for suitable sets $B$, $P(\cdot,B)$ is a function of
only the parameters $(\mathcal{A},\theta,\beta)$ and does not depend on
$\Sigma$. We first show a minorization condition. For $\ell\in(0,\infty)$,
define the function $g_{r}$ by

$g_{r}(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})=\inf_{V(\mathcal{A},\theta,\beta)\leq\ell}\pi_{n}(\mathcal{A}^{\prime}|\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})\pi_{n}(\theta^{\prime},\beta^{\prime}|\mathcal{A},{\Sigma}^{\prime})\pi_{n}({\Sigma}^{\prime}|\mathcal{A},\theta,\beta)$

and the constant $Z_{g_{\ell}}=\int
g_{\ell}(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})d\mathcal{A}^{\prime}d\theta^{\prime}d\beta^{\prime}d{\Sigma}^{\prime}.$
The drift function $V$ is a continuous, strongly convex function on a closed,
convex domain so its sublevel sets are closed and bounded (Nesterov, 2018,
Corollary 3.2.2). For fixed
$\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime}$, the function

$(\mathcal{A},\theta,\beta)\mapsto\pi_{n}(\theta^{\prime},\beta^{\prime}|\mathcal{A},{\Sigma}^{\prime})\pi_{n}({\Sigma}^{\prime}|\mathcal{A},\theta,\beta)$

is continuous and achieves its minimum over compact sets. Thus, $Z_{g_{\ell}}$
is not $0$ and we can define the probability measure

$\nu_{\ell}(\cdot)=Z_{g_{\ell}}^{-1}\int_{\cdot}g_{\ell}(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime}))d\mathcal{A}^{\prime}d\theta^{\prime}d\beta^{\prime}d{\Sigma}^{\prime}.$

For any $\ell\in(0,\infty)$ and any suitable set $B$,

$\displaystyle\inf_{\begin{subarray}{c}\Sigma\in\mathbb{R}^{m\times m},\\\
V(\mathcal{A},\theta,\beta)\leq\ell\end{subarray}}P\left((\mathcal{A},\theta,\beta,\Sigma),B\right)$
$\displaystyle=\inf_{V(\mathcal{A},\theta,\beta)\leq\ell}\int_{B}\pi_{n}(\mathcal{A}^{\prime}|\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})\pi_{n}(\theta^{\prime},\beta^{\prime}|\mathcal{A},{\Sigma}^{\prime})\pi_{n}({\Sigma}^{\prime}|\mathcal{A},\theta,\beta)d\mathcal{A}^{\prime}d\theta^{\prime}d\beta^{\prime}d{\Sigma}^{\prime}$
$\displaystyle\geq\int_{B}g_{\ell}(\mathcal{A}^{\prime},\theta^{\prime},\beta^{\prime},{\Sigma}^{\prime})d\mathcal{A}^{\prime}d\theta^{\prime}d\beta^{\prime}d{\Sigma}^{\prime}$
$\displaystyle=Z_{g_{\ell}}\nu_{\ell}(B).$

It remains to show a drift condition. Fix
$\mathcal{A}_{0},\theta_{0},\beta_{0}$, and fix $i\in 1,\ldots,n$. Since
$D_{i}$ is SPD, let $D_{i}=D_{i}^{1/2}D_{i}^{1/2}$,
$D_{i}^{-1}=D_{i}^{-1/2}D_{i}^{-1/2}$ where $D_{i}^{1/2},D_{i}^{-1/2}$ are
SPD. Using the identity

$d_{n,i}(\theta,\beta,\Sigma)=d_{i}+\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}\mathcal{B}\Sigma^{-1}\left(R_{i}-\Theta^{T}M_{i}-\mathcal{B}^{T}d_{i}\right)$

and taking the expectation with respect to
$\mathcal{A}_{i}|\mathcal{A}_{0},\theta_{0},\beta_{0},\theta,\beta,\Sigma$

$\displaystyle\mathbb{E}\left[\frac{1}{2}\left\lVert
D_{i}^{-1/2}(\mathcal{A}_{i}-d_{i})\right\rVert_{2}^{2}\big{|}\mathcal{A}_{0},\theta_{0},\beta_{0},\theta,\beta,\Sigma\right]$
(7a) $\displaystyle=\frac{1}{2}\left\lVert
D_{i}^{-1/2}\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}\mathcal{B}\Sigma^{-1}\left(R_{i}-\Theta^{T}M_{i}-\mathcal{B}^{T}d_{i}\right)\right\rVert_{2}^{2}$
(7b) $\displaystyle\hskip
11.38092pt+\frac{1}{2}\text{tr}\left[D_{i}^{-1/2}\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}D_{i}^{-1/2}\right].$
(7c)

Using singular value decomposition from (Horn and Johnson, 2012, Theorem
2.6.3), choose matrices $U_{i}\in\mathbb{R}^{p\times p}$,
$V_{i}\in\mathbb{R}^{m\times m}$ with $U_{i}^{T}U_{i}=U_{i}U_{i}^{T}=I_{p}$
and $V_{i}^{T}V_{i}$ $=V_{i}V_{i}^{T}=I_{m}$ and a rectangular diagonal matrix
$\Sigma_{i}\in\mathbb{R}^{p\times m}$ with diagonal nonnegative singular
values $(\sigma_{i,k})_{k}$ so that
$D_{i}^{1/2}\mathcal{B}\Sigma^{-1/2}=U_{i}\Sigma_{i}V_{i}^{T}$. Then

$\displaystyle\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}$
$\displaystyle=D_{i}^{1/2}(D_{i}^{1/2}\mathcal{B}\Sigma^{-1/2})(D_{i}^{1/2}\mathcal{B}\Sigma^{-1/2})^{T}+I_{p})^{-1}D_{i}^{1/2}$
$\displaystyle=D_{i}^{1/2}U_{i}(\Sigma_{i}\Sigma_{i}^{T}+I_{p})^{-1}U_{i}^{T}D_{i}^{1/2}.$
(8)

Using (8) and properties of the trace

$\displaystyle\frac{1}{2}\text{tr}\left[D_{i}^{-1/2}\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}D_{i}^{-1/2}\right]$
$\displaystyle=\frac{1}{2}\text{tr}\left[\left(\Sigma_{i}\Sigma_{i}^{T}+I_{p}\right)^{-1}U_{i}^{T}U_{i}\right]$
$\displaystyle\leq\frac{p}{2}.$

For $x\in[0,\infty)$ and $a\in(0,\infty)$, we have the inequality

$\displaystyle\frac{x}{x^{2}+a}\leq\frac{1}{2\sqrt{a}}.$ (9)

Using inequalities (8) and (9), the matrix norm is sub-multiplicative, and
$\left\lVert U_{i}\right\rVert_{2}=1$, we have

$\displaystyle\left\lVert
D_{i}^{-1/2}\left(\mathcal{B}\Sigma^{-1}\mathcal{B}^{T}+D_{i}^{-1}\right)^{-1}\mathcal{B}\Sigma^{-1/2}\Sigma^{-1/2}\right\rVert_{2}^{2}$
$\displaystyle=\left\lVert
U_{i}(\Sigma_{i}\Sigma_{i}^{T}+I_{p})^{-1}U_{i}^{T}D_{i}^{1/2}\mathcal{B}\Sigma^{-1/2}\Sigma^{-1/2}\right\rVert_{2}^{2}$
$\displaystyle\leq\left\lVert(\Sigma_{i}\Sigma_{i}^{T}+I_{p})^{-1}\Sigma_{i}\Sigma^{-1/2}\right\rVert_{2}^{2}$
$\displaystyle\leq\left[\frac{\sigma_{i}}{\sigma_{i}^{2}+1}\right]^{2}\left\lVert\Sigma^{-1/2}\right\rVert_{2}^{2}$
$\displaystyle\leq\frac{\left\lVert\Sigma^{-1}\right\rVert_{2}}{4}.$

Define the matrix $\tilde{X}=(d_{1},\ldots,d_{n})^{T}$. Applying these upper
bounds to (7) and combining for each $i\in 1,\ldots,n$,

$\displaystyle\mathbb{E}\left[\frac{1}{2}\sum_{i=1}^{n}\left\lVert
D_{i}^{-1/2}(\mathcal{A}_{i}-d_{i})\right\rVert_{2}^{2}\big{|}\mathcal{A}_{0},\theta_{0},\beta_{0},\theta,\beta,\Sigma\right]$
$\displaystyle\leq\frac{\left\lVert\Sigma^{-1}\right\rVert_{2}}{8}\left\lVert
R-M\Theta-\tilde{X}\mathcal{B}\right\rVert_{F}^{2}+\frac{pn}{2}.$

By convexity, for every $x,y$, $\left\lVert x-y\right\rVert_{2}^{2}\leq
2\left\lVert x\right\rVert_{2}^{2}+2\left\lVert y\right\rVert_{2}^{2}.$ Since
$C_{0}$ is SPD, let $C_{0}=C_{0}^{1/2}C_{0}^{1/2}$,
$C_{0}^{-1}=C_{0}^{-1/2}C_{0}^{-1/2}$ where $C_{0}^{1/2},C_{0}^{-1/2}$ are
SPD. Using convexity, and the matrix norm is sub-multiplicative, we have

$\displaystyle\frac{1}{2}\left\lVert
R-M\Theta-\tilde{X}\mathcal{B}\right\rVert_{F}^{2}$
$\displaystyle\leq\left\lVert
R\right\rVert_{F}^{2}+\left\lVert\begin{pmatrix}M&\tilde{X}\end{pmatrix}C_{0}^{1/2}\right\rVert_{2}^{2}\left\lVert
C_{0}^{-1/2}(\theta,\beta)^{T}\right\rVert_{2}^{2}.$

Therefore,

$\displaystyle\mathbb{E}\left[\frac{1}{2}\sum_{i=1}^{n}\left\lVert
D_{i}^{-1/2}(\mathcal{A}_{i}-d_{i})\right\rVert_{2}^{2}\big{|}\mathcal{A}_{0},\theta_{0},\beta_{0},\theta,\beta,\Sigma\right]$
(10a)
$\displaystyle\leq\frac{\left\lVert\Sigma^{-1}\right\rVert_{2}}{8}\left\lVert
R-M\Theta-\tilde{X}\mathcal{B}\right\rVert_{F}^{2}+\frac{pn}{2}$ (10b)
$\displaystyle\leq\left\lVert\Sigma^{-1}\right\rVert_{2}\frac{\left\lVert
R\right\rVert_{F}^{2}}{4}+\left\lVert\Sigma^{-1}\right\rVert_{2}\frac{\left\lVert\begin{pmatrix}M&\tilde{X}\end{pmatrix}C_{0}^{1/2}\right\rVert_{2}^{2}}{4}\left\lVert
C_{0}^{-1/2}(\theta,\beta)^{T}\right\rVert_{2}^{2}+\frac{pn}{2}.$ (10c)

Now taking the expectation with respect to
$\theta,\beta|\mathcal{A}_{0},\theta_{0},\beta_{0},\Sigma$

$\displaystyle\mathbb{E}\left[\frac{1}{2}\left\lVert
C_{0}^{-1/2}(\theta,\beta)^{T}\right\rVert_{2}^{2}\big{|}\Sigma,\mathcal{A}_{0},\theta_{0},\beta_{0}\right]$
$\displaystyle=\frac{1}{2}\left\lVert
C_{0}^{-1/2}c_{n}(\mathcal{A}_{0},\Sigma)\right\rVert_{2}^{2}+\frac{1}{2}\text{tr}(C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)C_{0}^{-1/2}).$

Using singular value decomposition (Horn and Johnson, 2012, Theorem 2.6.3),
choose matrices $U\in\mathbb{R}^{mn\times mn}$, $V\in\mathbb{R}^{m(p+q)\times
m(p+q)}$ with $U^{T}U=UU^{T}=I_{mn}$ and $V^{T}V=VV^{T}=I_{m(p+q)}$ and a
rectangular diagonal matrix $\Sigma_{\mathcal{A}_{0}}\in\mathbb{R}^{mn\times
m(p+q)}$ with diagonal nonnegative singular values
$(\sigma_{\mathcal{A}_{0},k})_{k}$ so that
$\Sigma^{-1/2}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}C_{0}^{1/2}=U\Sigma_{\mathcal{A}_{0}}V^{T}$.
We then have

$\displaystyle C_{n}(\mathcal{A}_{0},\Sigma)$
$\displaystyle=\left(\Sigma^{-1}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}^{T}\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}+C_{0}^{-1}\right)^{-1}$
$\displaystyle=C_{0}^{1/2}V\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}V^{T}C_{0}^{1/2}.$
(11)

Using (11) and properties of the trace

$\displaystyle\frac{1}{2}\text{tr}(C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)C_{0}^{-1/2})$
$\displaystyle\leq\frac{1}{2}\max_{k}\left[\left(\sigma_{\mathcal{A}_{0},k}^{2}+1\right)^{-1}\right]\text{tr}(V^{T}V)$
$\displaystyle\leq\frac{m(p+q)}{2}.$

Using convexity,

$\displaystyle\frac{1}{2}\left\lVert
C_{0}^{-1/2}c_{n}(\mathcal{A}_{0},\Sigma)\right\rVert_{2}^{2}$
$\displaystyle=\frac{1}{2}\left\lVert
C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)\left[\left[\Sigma^{-1}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}^{T}\right]\text{vec}(R)+C_{0}^{-1}c_{0}\right]\right\rVert_{2}^{2}$
$\displaystyle\leq\left\lVert
C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)\left[\Sigma^{-1/2}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}^{T}\right]\left[\Sigma^{-1/2}\otimes
I_{mn}\right]\right\rVert_{2}^{2}\left\lVert R\right\rVert_{F}^{2}$
$\displaystyle\hskip 14.22636pt+\left\lVert
C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)C_{0}^{-1}c_{0}\right\rVert_{2}^{2}.$

Using the inequality (9) and the identity (11),

$\displaystyle\left\lVert
C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)\left[\Sigma^{-1/2}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}\right]^{T}\left[\Sigma^{-1/2}\otimes
I_{mn}\right]\right\rVert_{2}^{2}$ $\displaystyle=\left\lVert
V\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}V^{T}\left[\Sigma^{-1/2}\otimes\begin{pmatrix}M&\mathcal{A}_{0}\end{pmatrix}C_{0}^{1/2}\right]^{T}\left[\Sigma^{-1/2}\otimes
I_{mn}\right]\right\rVert_{2}^{2}$ $\displaystyle=\left\lVert
V\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}\Sigma_{\mathcal{A}_{0}}^{T}U^{T}\right\rVert_{2}^{2}\left\lVert\Sigma^{-1}\right\rVert_{2}$
$\displaystyle\leq\left\lVert
V\right\rVert_{2}^{2}\left\lVert\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}\Sigma_{\mathcal{A}_{0}}^{T}\right\rVert_{2}^{2}\left\lVert
U^{T}\right\rVert_{2}^{2}\left\lVert\Sigma^{-1}\right\rVert_{2}$
$\displaystyle\leq\max_{k}\left(\frac{\sigma_{\mathcal{A}_{0},k}}{\sigma_{\mathcal{A}_{0},k}^{2}+1}\right)^{2}\left\lVert\Sigma^{-1}\right\rVert_{2}$
$\displaystyle\leq\frac{\left\lVert\Sigma^{-1}\right\rVert_{2}}{4}.$

Using (11),

$\displaystyle\left\lVert
C_{0}^{-1/2}C_{n}(\mathcal{A}_{0},\Sigma)C_{0}^{-1}c_{0}\right\rVert_{2}^{2}$
$\displaystyle=\left\lVert
V\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}V^{T}C_{0}^{-1/2}c_{0}\right\rVert_{2}^{2}$
$\displaystyle\leq\left\lVert
V\right\rVert_{2}^{2}\left\lVert\left(\Sigma_{\mathcal{A}_{0}}^{T}\Sigma_{\mathcal{A}_{0}}+I_{m(p+q)}\right)^{-1}\right\rVert_{2}^{2}\left\lVert
V^{T}\right\rVert_{2}^{2}c_{0}^{T}C_{0}^{-1}c_{0}$ $\displaystyle\leq
c_{0}^{T}C_{0}^{-1}c_{0}.$

Combining the upper bounds

$\displaystyle\mathbb{E}\left[\frac{1}{2}\left\lVert
C_{0}^{-1/2}(\theta,\beta)^{T}\right\rVert_{2}^{2}\big{|}\Sigma,\mathcal{A}_{0},\theta_{0},\beta_{0}\right]$
$\displaystyle\leq\frac{\left\lVert
R\right\rVert_{F}^{2}}{4}\left\lVert\Sigma^{-1}\right\rVert_{2}+c_{0}^{T}C_{0}^{-1}c_{0}+\frac{m(p+q)}{2}.$
(12)

Now using (10) and (12) and taking the iterated expectation with respect to
$\theta,\beta|\mathcal{A}_{0},\theta_{0},\beta_{0},\Sigma$,

$\displaystyle\mathbb{E}\left[\frac{1}{2}\sum_{i=1}^{n}\left\lVert
D_{i}^{-1/2}(\mathcal{A}_{i}-d_{i})\right\rVert_{2}^{2}\big{|}\Sigma,\mathcal{A}_{0},\theta_{0},\beta_{0}\right]$
(13a) $\displaystyle\leq\left\lVert\Sigma^{-1}\right\rVert_{2}$ (13b)
$\displaystyle\hskip 11.38092pt\times\left[\frac{\left\lVert
R\right\rVert_{F}^{2}}{4}+\frac{\left\lVert\begin{pmatrix}M&\tilde{X}\end{pmatrix}C_{0}^{1/2}\right\rVert_{2}^{2}m(p+q)}{4}+\frac{\left\lVert\begin{pmatrix}M&\tilde{X}\end{pmatrix}C_{0}^{1/2}\right\rVert_{2}^{2}}{2}c_{0}^{T}C_{0}^{-1}c_{0}\right]$
(13c) $\displaystyle\hskip
11.38092pt+\left\lVert\Sigma^{-1}\right\rVert_{2}^{2}\frac{\left\lVert\begin{pmatrix}M&\tilde{X}\end{pmatrix}C_{0}^{1/2}\right\rVert_{2}^{2}\left\lVert
R\right\rVert_{F}^{2}}{8}+\frac{pn}{2}.$ (13d)

Since $\Sigma^{-1}$ has a Wishart distribution and using properties of the
trace,

$\displaystyle\mathbb{E}\left\lVert\Sigma^{-1}\right\rVert_{2}$
$\displaystyle\leq\text{tr}\left[\mathbb{E}\left(\Sigma^{-1}\right)\right]$
$\displaystyle=(n+a_{0})\text{tr}\left[\left((R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})^{T}(R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})+B_{0}\right)^{-1}\right]$
$\displaystyle\leq(n+a_{0})tr[B_{0}^{-1}].$

Similarly, we use the second moment formula of the Wishart (Letac and Massam,
2004) to get the upper bound,

$\displaystyle\mathbb{E}\left\lVert\Sigma^{-1}\right\rVert_{2}^{2}$
$\displaystyle\leq\text{tr}\left[\mathbb{E}\left(\Sigma^{-2}\right)\right]$
$\displaystyle=(n+a_{0})\text{tr}\left[\left((R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})(R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})^{T}+B_{0}\right)^{-1}\right]^{2}$
$\displaystyle\hskip 14.22636pt+(n+a_{0})(n+a_{0}+1)$ $\displaystyle\hskip
14.22636pt\times\text{tr}\left[\left((R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})(R-M\Theta_{0}-\mathcal{A}_{0}\mathcal{B}_{0})^{T}+B_{0}\right)^{-2}\right]$
$\displaystyle\leq(n+a_{0})(n+a_{0}+1)tr[B_{0}^{-1}]^{2}.$

Taking the iterated expectation with respect to
$\Sigma|\mathcal{A}_{0},\theta_{0},\beta_{0}$ in (13), there is a constant
$L\in(0,\infty)$ so that the drift condition is satisfied with

$\mathbb{E}\left[V(\mathcal{A},\theta,\beta)|\mathcal{A}_{0},\theta_{0},\beta_{0}\right]\leq
L.$

∎

###  2.1 Bayesian EIV regression with errors in the features

Using Theorem 1, we develop geometrically ergodic Gibbs samplers for Bayesian
EIV regression with additive Gaussian error in the features. For the
remainder, we write the observed data as
$Y=(y_{1},\ldots,y_{n})^{T}\in\mathbb{R}^{n\times m}$,
$X=(x_{1},\ldots,x_{n})^{T}\in\mathbb{R}^{n\times p}$, and
$Z=(Z_{1},\ldots,Z_{n})^{T}\in\mathbb{R}^{n\times q}$. Consider the Bayesian
EIV regression (1) with Berkson errors and priors (3) and (4). We will write
the posterior density $\pi_{n}$ for this Bayesian model as

$\displaystyle\pi_{n}(\mathcal{A},\theta,\beta,\Sigma)$
$\displaystyle\propto\left(\frac{1}{\det(\Sigma)}\right)^{(n+a_{0}+m+1)/2}\exp\left[-\frac{1}{2}\text{tr}[\Sigma^{-1}B_{0}]\right]$
$\displaystyle\hskip
14.22636pt\times\exp\left[-\frac{1}{2}\sum_{i=1}^{n}(y_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})^{T}\Sigma^{-1}(y_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})\right]$
$\displaystyle\hskip
14.22636pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-x_{i})^{T}V_{i}^{-1}(\mathcal{A}_{i}-x_{i})\right)$
$\displaystyle\hskip
14.22636pt\times\exp\left(-\frac{1}{2}((\theta,\beta)^{T}-j_{0})^{T}J_{0}^{-1}((\theta,\beta)^{T}-j_{0})\right).$

This posterior density is a special case of the general density (5) choosing
$M\equiv Z$, $R\equiv Y$, $c_{0},C_{0}\equiv j_{0},J_{0}$, and
$d_{i},D_{i}\equiv x_{i},V_{i}$.

We can define a 3-variable deterministic scan Gibbs sampler which generates a
Markov chain
$(\mathcal{A}_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for this
posterior density as a special case of the Gibbs sampler constructed in
Section 2. Initialize $(\mathcal{A}_{0},\theta_{0},\beta_{0},\Sigma_{0})$ and
for $t\in 1,\ldots$,

  1. 1.

Generate
$\Sigma_{t}|\mathcal{A}_{t-1},\theta_{t-1},\beta_{t-1}\sim\mathcal{W}^{-1}\left(n+a_{0},B_{n,t}\right)$
where

$B_{n,t}=(Y-Z\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})^{T}(Y-Z\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})+B_{0}$

  2. 2.

Generate $(\theta_{t},\beta_{t})^{T}|\mathcal{A}_{t-1},\Sigma_{t}\sim
N_{m(p+q)}(j_{n,t},J_{n,t})$ where

$\displaystyle
J_{n,t}=\left(\Sigma^{-1}_{t}\otimes\begin{pmatrix}Z&\mathcal{A}_{t-1}\end{pmatrix}^{T}\begin{pmatrix}Z&\mathcal{A}_{t-1}\end{pmatrix}+J_{0}^{-1}\right)^{-1}$
$\displaystyle
j_{n,t}=J_{n,t}\left[\left[\Sigma^{-1}_{t}\otimes\begin{pmatrix}Z&\mathcal{A}_{t-1}\end{pmatrix}^{T}\right]\text{vec}(Y)+J_{0}^{-1}j_{0}\right]$

  3. 3.

Generate $\mathcal{A}_{i,t}|\theta_{t},\beta_{t},\Sigma_{t}\sim
N_{p}(k_{n,i,t},K_{n,i,t}),i\in 1,\ldots,n$ where

$\displaystyle
K_{n,i,t}=\left(\mathcal{B}_{t}\Sigma^{-1}_{t}\mathcal{B}_{t}^{T}+V_{i}^{-1}\right)^{-1}$
$\displaystyle
k_{n,i,t}=K_{n,i,t}\left[V_{i}^{-1}x_{i}+\mathcal{B}_{t}\Sigma^{-1}_{t}\left(y_{i}-\Theta_{t}^{T}Z_{i}\right)\right].$

Applying Theorem 1, we have the following result.

######  Corollary 1.

The 3-variable Gibbs sampler
$(\mathcal{A}_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for the
posterior in Bayesian EIV regression (1) with Berkson errors and priors (3)
and (4) is geometrically ergodic.

Now consider Bayesian EIV regression (1) with additive Gaussian error in
$X_{i}$ using classical errors and priors (3) and (4). The posterior has
density

$\displaystyle\pi_{n}(\mathcal{A},\theta,\beta,\Sigma)$
$\displaystyle\propto\left(\frac{1}{\det(\Sigma)}\right)^{(n+a_{0}+m+1)/2}\exp\left[-\frac{1}{2}\text{tr}[\Sigma^{-1}B_{0}]\right]$
$\displaystyle\hskip
14.22636pt\times\exp\left[-\frac{1}{2}\sum_{i=1}^{n}(y_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})^{T}\Sigma^{-1}(y_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})\right]$
$\displaystyle\hskip
14.22636pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-k_{i}^{\prime})^{T}(V_{i}^{-1}+K_{i}^{-1})(\mathcal{A}_{i}-k_{i}^{\prime})\right)$
$\displaystyle\hskip
14.22636pt\times\exp\left(-\frac{1}{2}((\theta,\beta)^{T}-j_{0})^{T}J_{0}^{-1}((\theta,\beta)^{T}-j_{0})\right)$

where
$k_{i}^{\prime}=(V_{i}^{-1}+K_{i}^{-1})^{-1}\left[V_{i}^{-1}x_{i}+K_{i}^{-1}k_{i}\right].$
The posterior density is also a special case of the general density (5) when
$Z\equiv M$, $R\equiv Y$, and $c_{0},C_{0}\equiv j_{0},J_{0}$, and
$d_{i},D_{i}\equiv k_{i}^{\prime},(V_{i}^{-1}+K_{i}^{-1})^{-1}$.

We define a 3-variable deterministic scan Gibbs sampler similarly. Initialize
$(\mathcal{A}_{0},\theta_{0},\beta_{0},\Sigma_{0})$ and for $t\in 1,\ldots$,

  1. 1.

Generate
$\Sigma_{t}|\mathcal{A}_{t-1},\theta_{t-1},\beta_{t-1}\sim\mathcal{W}^{-1}\left(n+a_{0},B_{n,t}\right)$

  2. 2.

Generate $(\theta_{t},\beta_{t})^{T}|\mathcal{A}_{t-1},\Sigma_{t}\sim
N_{m(p+q)}(j_{n,t},J_{n,t})$

  3. 3.

Generate $\mathcal{A}_{i,t}|\theta_{t},\beta_{t},\Sigma_{t}\sim
N_{p}(k^{\prime}_{n,i,t},K^{\prime}_{n,i,t}),i\in 1,\ldots,n$ where

$\displaystyle
K^{\prime}_{n,i,t}=\left(\mathcal{B}_{t}\Sigma^{-1}_{t}\mathcal{B}_{t}^{T}+V_{i}^{-1}+K_{i}^{-1}\right)^{-1}$
$\displaystyle
k^{\prime}_{n,i,t}=K^{\prime}_{n,i,t}\left[V_{i}^{-1}x_{i}+K_{i}^{-1}k_{i}+\mathcal{B}_{t}\Sigma^{-1}_{t}\left(y_{i}-\Theta_{t}^{T}Z_{i}\right)\right].$

We also have the following as a direct result of Theorem 1.

######  Corollary 2.

The 3-variable Gibbs sampler
$(\mathcal{A}_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for the
posterior in Bayesian EIV regression (1) with classical errors and priors (3)
and (4) is geometrically ergodic.

###  2.2 Bayesian EIV regression with errors in the response and features

Similar to the previous section, we develop geometrically ergodic Gibbs
samplers for Bayesian EIV regression with additional additive Gaussian error
in the features and response. Consider the Bayesian EIV regression (2) with
Berkson errors in $X_{i}$ and additional error in $Y_{i}$ along with priors
(3) and (4). Let
$\mathcal{V}=(\mathcal{V}_{1},\ldots,\mathcal{V}_{n})^{T}\in\mathbb{R}^{n\times
m}$ and $\nu=\text{vec}(\mathcal{V})$, and let
$U_{0}=\text{blockdiag}(U_{i})\in\mathbb{R}^{mn\times mn}$. The Bayesian
posterior $\Pi_{n}$ has density

$\displaystyle\pi_{n}(\mathcal{A},\nu,\theta,\beta,\Sigma)$
$\displaystyle\propto\left(\frac{1}{\det(\Sigma)}\right)^{(n+a_{0}+m+1)/2}\exp\left[-\frac{1}{2}\text{tr}[\Sigma^{-1}B_{0}]\right]$
$\displaystyle\hskip
11.38092pt\times\exp\left[-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{V}_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})^{T}\Sigma^{-1}(\mathcal{V}_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})\right]$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{V}_{i}-y_{i})^{T}U_{i}^{-1}(\mathcal{V}_{i}-y_{i})\right)$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-x_{i})^{T}V_{i}^{-1}(\mathcal{A}_{i}-x_{i})\right)$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}((\theta,\beta)^{T}-j_{0})^{T}J_{0}^{-1}((\theta,\beta)^{T}-j_{0})\right).$

This posterior density is a special case of the density (5) when redefining
$\tilde{\theta}\equiv(\nu,\theta)^{T}$,
$M\equiv\begin{pmatrix}-I&Z\end{pmatrix}$, $r\equiv 0$, $c_{0}=(Y,j_{0})^{T}$,

$C_{0}\equiv\begin{pmatrix}U_{0}&0\\\ 0&J_{0}\end{pmatrix},$

and $d_{i},D_{i}\equiv x_{i},V_{i}$.

We define a 3-variable deterministic scan Gibbs sampler which generates a
Markov chain
$(\mathcal{A}_{t},\nu_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for
this posterior density. Initialize $\mathcal{A}_{0}$, $\nu_{0}$, $\theta_{0}$,
$\beta_{0}$, $\Sigma_{0}$ and for $t\in 1,\ldots$,

  1. 1.

Generate
$\Sigma_{t}|\mathcal{A}_{t-1},\nu_{t-1},\theta_{t-1},\beta_{t-1}\sim\mathcal{W}^{-1}\left(n+a_{0},B^{\prime}_{n,t}\right)$
where

$B_{n,t}^{\prime}=(\mathcal{V}_{t-1}-Z\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})^{T}(\mathcal{V}_{t-1}-Z\Theta_{t-1}-\mathcal{A}_{t-1}\mathcal{B}_{t-1})+B_{0}$

  2. 2.

Generate $(\nu_{t},\theta_{t},\beta_{t})^{T}|\mathcal{A}_{t-1},\Sigma_{t}\sim
N_{p+q}(j^{\prime}_{n,t},J^{\prime}_{n,t})$ where

$\displaystyle
J^{\prime}_{n,t}=\left(\Sigma^{-1}_{t}\otimes\begin{pmatrix}-I&Z&\mathcal{A}_{t-1}\end{pmatrix}^{T}\begin{pmatrix}-I&Z&\mathcal{A}_{t-1}\end{pmatrix}+\begin{pmatrix}U_{0}^{-1}&0\\\
0&J_{0}^{-1}\end{pmatrix}\right)^{-1}$ $\displaystyle
j^{\prime}_{n,t}=J_{n,t}^{\prime}\begin{pmatrix}U_{0}^{-1}&0\\\
0&J_{0}^{-1}\end{pmatrix}(\text{vec}(Y),j_{0})^{T}$

  3. 3.

Generate $\mathcal{A}_{i,t}|\nu_{t},\theta_{t},\beta_{t},\Sigma_{t}\sim
N_{p}(k^{\prime\prime}_{n,i,t},K^{\prime\prime}_{n,i,t}),i\in 1,\ldots,n$
where

$\displaystyle
K^{\prime\prime}_{n,i,t}=\left(\mathcal{B}_{t}\Sigma^{-1}_{t}\mathcal{B}_{t}^{T}+V_{i}^{-1}\right)^{-1}$
$\displaystyle
k^{\prime\prime}_{n,i,t}=K^{\prime\prime}_{n,i,t}\left[V_{i}^{-1}x_{i}+\mathcal{B}_{t}\Sigma^{-1}_{t}\left(\mathcal{V}_{i,t}-\Theta_{t}^{T}Z_{i}\right)\right].$

Using Theorem 1, we have the following result.

######  Corollary 3.

The 3-variable Gibbs sampler
$(\mathcal{A}_{t},\nu_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for
Bayesian EIV regression (2) with Berkson errors and priors (3) and (4) is
geometrically ergodic.

Now consider the Bayesian EIV regression (2) with classical errors in $X_{i}$
and additional error in $Y_{i}$ with priors (3) and (4). The posterior
$\Pi_{n}$ for this Bayesian model has density

$\displaystyle\pi_{n}(\mathcal{A},\nu,\theta,\beta,\Sigma)$
$\displaystyle\propto\left(\frac{1}{\det(\Sigma)}\right)^{(n+a_{0}+m+1)/2}\exp\left[-\frac{1}{2}\text{tr}[\Sigma^{-1}B_{0}]\right]$
$\displaystyle\hskip
11.38092pt\times\exp\left[-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{V}_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})^{T}\Sigma^{-1}(\mathcal{V}_{i}-\Theta^{T}Z_{i}-\mathcal{B}^{T}\mathcal{A}_{i})\right]$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{V}_{i}-y_{i})^{T}U_{i}^{-1}(\mathcal{V}_{i}-y_{i})\right)$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}\sum_{i=1}^{n}(\mathcal{A}_{i}-k_{i}^{\prime})^{T}(V_{i}^{-1}+K_{i}^{-1})(\mathcal{A}_{i}-k_{i}^{\prime})\right)$
$\displaystyle\hskip
11.38092pt\times\exp\left(-\frac{1}{2}((\theta,\beta)^{T}-j_{0})^{T}J_{0}^{-1}((\theta,\beta)^{T}-j_{0})\right).$

This posterior density is also a special case of the density (5) when
redefining $\tilde{\theta}\equiv(\nu,\theta)^{T}$,
$M\equiv\begin{pmatrix}-I&Z\end{pmatrix}$, $R\equiv 0$, $c_{0}=(y,j_{0})^{T}$,

$C_{0}\equiv\begin{pmatrix}U_{0}&0\\\ 0&J_{0}\end{pmatrix},$

and $d_{i},D_{i}\equiv k_{i}^{\prime},(V_{i}^{-1}+K_{i}^{-1})^{-1}$.

We define a 3-variable deterministic scan Gibbs sampler similarly. Initialize
$(\mathcal{A}_{0},\nu_{0},\theta_{0},\beta_{0},\Sigma_{0})$ and for $t\in
1,\ldots$,

  1. 1.

Generate
$\Sigma_{t}|\mathcal{A}_{t-1},\nu_{t-1},\theta_{t-1},\beta_{t-1}\sim\mathcal{W}^{-1}\left(n+a_{0},B^{\prime}_{n,t}\right)$

  2. 2.

Generate $(\nu_{t},\theta_{t},\beta_{t})^{T}|\mathcal{A}_{t-1},\Sigma_{t}\sim
N_{p+q}(j^{\prime}_{n,t},J^{\prime}_{n,t})$

  3. 3.

Generate $\mathcal{A}_{i,t}|\nu_{t},\theta_{t},\beta_{t},\Sigma_{t}\sim
N_{p}(k^{\prime\prime\prime}_{n,i,t},K^{\prime\prime\prime}_{n,i,t}),i\in
1,\ldots,n$ where

$\displaystyle
K^{\prime\prime\prime}_{n,i,t}=\left(\mathcal{B}_{t}\Sigma_{t}\mathcal{B}_{t}^{T}+V_{i}^{-1}+K_{i}^{-1}\right)^{-1}$
$\displaystyle
k^{\prime\prime\prime}_{n,i,t}=K^{\prime\prime\prime}_{n,i,t}\left[V_{i}^{-1}x_{i}+K_{i}^{-1}k_{i}+\mathcal{B}_{t}\Sigma^{-1}_{t}\left(\mathcal{V}_{i,t}-\Theta_{t}^{T}Z_{i}\right)\right].$

Using Theorem 1, we have the following result.

######  Corollary 4.

The 3-variable Gibbs sampler
$(\mathcal{A}_{t},\nu_{t},\theta_{t},\beta_{t},\Sigma_{t})_{t=0}^{\infty}$ for
Bayesian EIV regression (2) with classical errors and priors (3) and (4) is
geometrically ergodic.

##  3 Simulations

###  3.1 Limitations of the Gibbs sampler in large problem sizes

Theoretically, we developed a qualitative convergence result for the Gibbs
sampler in Bayesian multivariate EIV regression. It is important in practice
to understand the relationship between scaling of the problem size and the
estimation reliability from the Gibbs sampler. We look at artificially
generated data to empirically demonstrate the dependence of the Gibbs sampler
when the the dimension of the response $m$ and the dimension of the features
$p$ are increasing in configurations $(m,p)=(1,1),(2,7),(3,7)$ in the Bayesian
posterior. Artificial data is generated according to the multivariate EIV
Berkson linear regression model for $i=1,\ldots,50$ with

$\displaystyle\Sigma\sim\mathcal{W}^{-1}(m,10^{-3}I)$
$\displaystyle(\theta,\beta)^{T}\sim N_{m+mp}(0,10^{3}I)$ $\displaystyle
X_{i}|\mathcal{A}_{i}\sim N_{p}(\mathcal{A}_{i},.2I)$ $\displaystyle
Y_{i}|\mathcal{A}_{i},\Theta,\mathcal{B},\Sigma\sim
N_{m}(\Theta^{T}1+\mathcal{B}^{T}\mathcal{A}_{i},\Sigma)$

We simulate $T=10^{5}$ MCMC realizations from the Gibbs sampler in each
configuration using $10^{4}$ realizations for burn-in and analyze diagnostics
for $\beta_{t}$ taking values in $\mathbb{R}^{mp}$. We independently replicate
the simulation $5$ times to reduce variability. Geometric ergodicity
guarantees the properly scaled and summed samples from the Gibbs sampler

$\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\left[\beta_{t}-\int\beta\Pi_{n}(d\beta)\right]\to
N(0,\Sigma_{*})$

as $T\to\infty$ in distribution where $\Sigma_{*}$ is a SPD covariance matrix.
The central limit theorem can be seen to hold in this case due to the
established drift condition with the drift function. Figure 1(a) and Figure
1(b) plot the largest and smallest eigenvalues of a batch means estimate to
the multivariate standard error matrix ${\Sigma_{*}}^{1/2}$ in the central
limit theorem. This simulation shows an increase in the largest eigenvalue in
larger size problems which can lead to suggesting more iterations are needed
for appropriate estimation in practice. Figure 1(c) plots the multivariate
effective sample size (Vats, Flegal and Jones, 2019b). We can also see a
relatively sharp decrease in the estimation of the effective sample size as
the problem size increases suggesting the algorithm should be run for many
iterations even in moderately sized problems. This simulation demonstrates
that even though the algorithm is always geometrically ergodic and can scale
reasonably well to larger problem sizes, the Gibbs sampler generally requires
many more iterations for reliable estimation even in moderately sized
problems. Simulation code is available using Python (Van Rossum and Drake Jr,
1995) and (Harris et al., 2020) for matrix calculations at
https://github.com/austindavidbrown/BayesEIV.

(a)

(b)

(c)

Figure 1: (a), (b) Largest and smallest eigenvalues of the MCMC standard error
matrix targeting the average of $\beta$ for iterations of the Gibbs sampler
and (c) the multivariate effective sample size for iterations of the Gibbs
sampler

###  3.2 Robustness to model misspecification

Although the multivariate Bayesian model for EIV accounts for additional error
in the features, this error can be misspecified. In particular, the error
$X_{i}|\mathcal{A}_{i}$ from the model in Section 3.1 may be a multivariate
$t$ distribution with heavier tails in practical problems. We are interested
to empirically study the robustness of the convergence of the Gibbs sampler to
misspecification in this modeling error. Denote $t_{d}(v,m,V)$ as a
multivariate $t$ distribution in dimension $d$ with $v$ degrees of freedom,
location vector $m$, and scale matrix $V$. With $df$ denoting the degrees of
freedom, artificial data is generated according to the misspecified
multivariate EIV Berkson linear regression model for $i=1,\ldots,50$ with

$\displaystyle\Sigma\sim\mathcal{W}^{-1}(3,10^{-3}I)$
$\displaystyle(\theta,\beta)^{T}\sim N_{3+9}(0,10^{3}I)$ $\displaystyle
X_{i}|\mathcal{A}_{i}\sim\text{t}_{3}(\text{df},\mathcal{A}_{i},.2I)$
$\displaystyle Y_{i}|\mathcal{A}_{i},\Theta,\mathcal{B},\Sigma\sim
N_{3}(\Theta^{T}1+\mathcal{B}^{T}\mathcal{A}_{i},\Sigma).$

We look at compare more dispersed tail behavior with $df=2$ and less dispersed
tail behavior with $df=10$. We replicate a simulation $5$ times where we
simulate $T=10^{5}$ MCMC realizations from the Gibbs sampler in each
configuration using $10^{4}$ realizations for burn-in and analyze diagnostics
for $\beta_{t}$. Similar to the previous simulation in Section 3.1, Figure
2(a) and Figure 2(b) plot the largest and smallest eigenvalues of a batch
means estimate to the multivariate standard error matrix and Figure 2(c) plots
the multivariate effective sample size (Vats, Flegal and Jones, 2019b). We can
see similar behavior in the estimation from the Gibbs sampler based on the
both smaller and larger degrees of freedom. The simulation results suggest the
Gibbs sampler is reasonably robust to misspecification of the tails in the
error distribution of the features for $X_{i}|\mathcal{A}_{i}$.

(a)

(b)

(c)

Figure 2: (a), (b) Largest and smallest eigenvalues of the MCMC standard error
matrix (c) the multivariate effective sample size for iterations of the Gibbs
sampler

##  4 Real-data example: measurement error in astrophysics

We look at Bayesian EIV linear regression proposed and analyzed in (Harris,
Poole and Harris, 2014; Hilbe, de Souza and Ishida, 2017). The dataset
consists of the central galaxy supermassive black hole mass and the stellar
bulge velocity dispersion from $n=46$ different galaxies (Harris, Poole and
Harris, 2014). The response $Y_{i}$ is the logarithm of the observed central
black hole mass and the predictor variable $X_{i}$ is the logarithm of the
observed velocity dispersion. The measurement errors are known beforehand and
denoted by $\epsilon_{Y_{i}}$ and $\epsilon_{X_{i}}$ for both the response and
predictor variables. The EIV linear regression model studied in (Hilbe, de
Souza and Ishida, 2017) folllows

$\displaystyle\sigma^{2}\sim\text{Inverse-gamma}(10^{-3},10^{-3})$
$\displaystyle\alpha\sim N_{1}(0,10^{3}),\beta\sim N_{1}(0,10^{3})$
$\displaystyle\mathcal{A}_{i}\sim N_{1}(0,10^{3}),X_{i}|\mathcal{A}_{i}\sim
N_{1}(\mathcal{A}_{i},\epsilon_{X_{i}}^{2})$
$\displaystyle\mathcal{V}_{i}|\mathcal{A}_{i},\theta,\beta,\sigma^{2}\sim
N_{1}(\theta+\mathcal{A}_{i}\beta,\sigma^{2}),Y_{i}|\mathcal{V}_{i}\sim
N_{1}(\mathcal{V}_{i},\epsilon_{y_{i}}^{2}).$

We generate $10^{5}$ MCMC realizations from the Gibbs sampler. Figure 3 plots
the autocorrelation, estimates to the standard errors in the central limit
theorem, and effective sample sizes from these realizations. The
autocorrelations are computed up to lag $20$. Overall, we see the Gibbs
sampler performs well. However, the standard error and effective sample size
plots suggest that even though the Gibbs sampler is geometrically ergodic,
many iterations are still recommended even in low dimensions. These figures
suggest empirical diagnostics for the regression parameter $\beta$ as a
reasonable choice as opposed to the other parameters $\alpha,\sigma^{2}$ to
determine the reliability of the algorithm in practice.

(a)

(b)

(c)

Figure 3: (a) Autocorrelation for each regression parameter (b) MCMC standard
errors for the regression parameters (c) MCMC effective sample size plots for
each regression parameter

##  5 Conclusion and Future Directions

We showed using a $3$-variable deterministic scan Gibbs sampler to sample the
posterior in $4$ different multivariate Bayesian EIV regression models with
additive Gaussian errors and independent priors is always geometrically
ergodic. This is of pragmatic importance to practitioners as trustworthy
estimation from a Gibbs sampler is dependent on the speed of convergence of
the Markov chain. More specifically, time averages from the Markov chains have
many practically relevant theoretical guarantees such as a central limit
theorem. Secondly, these Gibbs samplers can be simulated efficiently without
the need for complex, intermediate Metropolis-Hastings or rejection sampling
steps. One drawback, however, is our convergence analysis is qualitative as we
do not construct an explicit convergence rate.

There are many future research directions in studying the convergence of Gibbs
samplers in EIV models. It appears reasonable that some Gibbs samplers for
generalized linear models such as the Pólya-Gamma sampler will also be
geometrically ergodic (Choi and Hobert, 2013; Polson, Scott and Windle, 2013;
Wang and Roy, 2018). It seems also interesting to look at alternative errors
in the variables such as non-Gaussian or non-additive errors.

[Acknowledgments] Thanks to Galin L. Jones for helpful guidance and insightful
comments for developing this article. Thanks to the anonymous referees for
recommendations to improve and extend some results from a previous version of
this manuscript.

## References

  * Achic et al. (2018) [author] Achic, Betsabé G. BlasB. G. B., Wang, TianyingT., Su, YaY., Kipnis, VictorV., Dodd, KevinK. and Carroll, Raymond J.R. J. (2018). Categorizing a continuous predictor subject to measurement error. Electronic Journal of Statistics 12 4032 – 4056.
  * Berkson (1950) [author] Berkson, JosephJ. (1950). Are There Two Regressions? Journal of the American Statistical Association 45 164–180.
  * Bhadra and Carroll (2016) [author] Bhadra, AnindyaA. and Carroll, Raymond J.R. J. (2016). Exact sampling of the unobserved covariates in Bayesian spline models for measurement error problems. Statistics and Computing 26 827–840.
  * Bresson et al. (2021) [author] Bresson, GeorgesG., Chaturvedi, AnoopA., Rahman, Mohammad ArshadM. A. and Shalabh (2021). Seemingly unrelated regression with measurement error: estimation via Markov Chain Monte Carlo and mean field variational Bayes approximation. The International Journal of Biostatistics 17 75–97.
  * Buonaccorsi (2010) [author] Buonaccorsi, J. P.J. P. (2010). Measurement Error: Models, Methods, and Applications, 1 ed. Chapman and Hall/CRC.
  * Carroll et al. (2006) [author] Carroll, R. J.R. J., Ruppert, D.D., Stefanski, L. A.L. A. and Crainiceanu, C. M.C. M. (2006). Measurement Error in Nonlinear Models: A Modern Perspective, 2 ed. Chapman and Hall/CRC.
  * Chan and Geyer (1994) [author] Chan, Kung SikK. S. and Geyer, Charles J.C. J. (1994). Discussion: Markov Chains for Exploring Posterior Distributions. The Annals of Statistics 22 1747–1758.
  * Charisse Farr et al. (2020) [author] Charisse Farr, A.A., Mengersen, KerrieK., Ruggeri, FabrizioF., Simpson, DanielD., Wu, PaulP. and Yarlagadda, PrasadP. (2020). Combining Opinions for Use in Bayesian Networks: A Measurement Error Approach. International Statistical Review 88 335 – 353.
  * Choi and Hobert (2013) [author] Choi, Hee MinH. M. and Hobert, James P.J. P. (2013). The Pólya-Gamma Gibbs sampler for Bayesian logistic regression is uniformly ergodic. Electronic Journal of Statistics 7 2054 – 2064.
  * Clayton et al. (1992) [author] Clayton, DGD. et al. (1992). Models for the analysis of cohort and case-control studies with inaccurately measured exposures. Statistical models for longitudinal studies of health 301–331.
  * Damgaard (2020) [author] Damgaard, ChristianC. (2020). Measurement Uncertainty in Ecological and Environmental Models. Trends in Ecology and Evolution 35 871–873.
  * Dellaportas and Stephens (1995) [author] Dellaportas, PetrosP. and Stephens, David A.D. A. (1995). Bayesian Analysis of Errors-in-Variables Regression Models. Biometrics 51 1085–1095.
  * Doss et al. (2014) [author] Doss, Charles RC. R., Flegal, James MJ. M., Jones, Galin LG. L. and Neath, Ronald CR. C. (2014). Markov chain Monte Carlo estimation of quantiles. Electronic Journal of Statistics 8 2448–2478.
  * Ekvall and Jones (2021) [author] Ekvall, Karl OskarK. O. and Jones, Galin L.G. L. (2021). Convergence analysis of a collapsed Gibbs sampler for Bayesian vector autoregressions. Electronic Journal of Statistics 15 691–721.
  * Fang et al. (2017) [author] Fang, XingX., Li, BofengB., Alkhatib, HamzaH., Zeng, WenxianW. and Yao, YibinY. (2017). Bayesian inference for the Errors-In-Variables model. Studia Geophysica et Geodaetica 61 1573–1626.
  * Feigelson and Babu (1992) [author] Feigelson, Eric DE. D. and Babu, Gutti JG. J. (1992). Linear regression in astronomy. II. The Astrophysical Journal 397 55–67.
  * Flegal and Jones (2010) [author] Flegal, James M.J. M. and Jones, Galin L.G. L. (2010). Batch means and spectral variance estimators in Markov chain Monte Carlo. The Annals of Statistics 38 1034 – 1070.
  * Fuller (1987) [author] Fuller, Wayne A.W. A. (1987). Measurement Error Models. John Wiley.
  * Geman and Geman (1984) [author] Geman, StuartS. and Geman, DonaldD. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence 6 721–741.
  * Goodfellow, Shlens and Szegedy (2015) [author] Goodfellow, IanI., Shlens, JonathonJ. and Szegedy, ChristianC. (2015). Explaining and Harnessing Adversarial Examples. International Conference on Learning Representations.
  * Groß (2016) [author] Groß, MarcusM. (2016). Modeling body height in prehistory using a spatio-temporal Bayesian errors-in variables model. AStA Advances in Statistical Analysis 100 289–311.
  * Gustafson (2003) [author] Gustafson, PaulP. (2003). Measurement error and misclassification in statistics and epidemiology: impacts and Bayesian adjustments. CRC Press.
  * Hairer and Mattingly (2011) [author] Hairer, MartinM. and Mattingly, Jonathan C.J. C. (2011). Yet Another Look at Harris’ Ergodic Theorem for Markov Chains. Seminar on Stochastic Analysis, Random Fields and Applications VI 63\.
  * Harris, Poole and Harris (2014) [author] Harris, G.G., Poole, G.G. and Harris, W.W. (2014). Globular clusters and supermassive black holes in galaxies: further analysis and a larger sample. Monthly Notices of the Royal Astronomical Society 438 2117–2130.
  * Harris et al. (2020) [author] Harris, Charles R.C. R., Millman, K. JarrodK. J., van der Walt, Stéfan JS. J., Gommers, RalfR., Virtanen, PauliP., Cournapeau, DavidD., Wieser, EricE., Taylor, JulianJ., Berg, SebastianS., Smith, Nathaniel J.N. J., Kern, RobertR., Picus, MattiM., Hoyer, StephanS., van Kerkwijk, Marten H.M. H., Brett, MatthewM., Haldane, AllanA., Fernández del Río, JaimeJ., Wiebe, MarkM., Peterson, PearuP., Gérard-Marchant, PierreP., Sheppard, KevinK., Reddy, TylerT., Weckesser, WarrenW., Abbasi, HameerH., Gohlke, ChristophC. and Oliphant, Travis E.T. E. (2020). Array programming with NumPy. Nature 585 357–362. 10.1038/s41586-020-2649-2
  * Hilbe, de Souza and Ishida (2017) [author] Hilbe, Joseph M.J. M., de Souza, Rafael S.R. S. and Ishida, Emille E. O.E. E. O. (2017). Bayesian Models for Astrophysical Data: Using R, JAGS, Python, and Stan. Cambridge University Press.
  * Horn and Johnson (2012) [author] Horn, Roger A.R. A. and Johnson, Charles R.C. R. (2012). Matrix Analysis. Cambridge University Press.
  * Huang (2010) [author] Huang, Hung-JenH.-J. (2010). Bayesian Analysis of Errors-in-Variables Growth Curve Models. PhD Dissertation.
  * Håvard Rue (2009) [author] Håvard Rue, Nicolas ChopinN. C. Sara Martino (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 71 319–392.
  * Jones (2004) [author] Jones, Galin LG. L. (2004). On the Markov chain central limit theorem. Probability surveys 1 299–320.
  * Kelly (2012) [author] Kelly, Brandon CB. C. (2012). Measurement error models in astronomy. In Statistical challenges in modern astronomy V 147–162. Springer.
  * Kröger, Hoffmann and Pakpahan (2016) [author] Kröger, HannesH., Hoffmann, RasmusR. and Pakpahan, EduwinE. (2016). Consequences of measurement error for inference in cross-lagged panel design—the example of the reciprocal causal relationship between subjective health and socio-economic status. Journal of the Royal Statistical Society. Series A (Statistics in Society) 179 607–628.
  * Letac and Massam (2004) [author] Letac, GérardG. and Massam, HélèneH. (2004). All Invariant Moments of the Wishart Distribution. Scandinavian Journal of Statistics 31 295–318.
  * Mallick and Gelfand (1996) [author] Mallick, Bani K.B. K. and Gelfand, Alan E.A. E. (1996). Semiparametric errors-in-variables models A Bayesian approach. Journal of Statistical Planning and Inference 52 307–321.
  * Meyn and Tweedie (2009) [author] Meyn, Sean P.S. P. and Tweedie, Richard L.R. L. (2009). Markov Chains and Stochastic Stability, 2 ed. Cambridge University Press, USA.
  * Michalek and Tripathi (1980) [author] Michalek, Joel E.J. E. and Tripathi, Ram C.R. C. (1980). The Effect of Errors in Diagnosis and Measurement on the Estimation of the Probability of an Event. Journal of the American Statistical Association 75 713–721.
  * Muff et al. (2015) [author] Muff, StefanieS., Riebler, AndreaA., Held, LeonhardL., Rue, HåvardH. and Saner, PhilippeP. (2015). Bayesian analysis of measurement error models using integrated nested Laplace approximations. Journal of the Royal Statistical Society. Series C (Applied Statistics) 64 231–252.
  * Nesterov (2018) [author] Nesterov, YuriiY. (2018). Lectures on Convex Optimization, 2 ed. Springer International Publishing.
  * Pham, Ormerod and Wand (2013) [author] Pham, Tung HT. H., Ormerod, John TJ. T. and Wand, Matthew PM. P. (2013). Mean field variational Bayesian inference for nonparametric regression with measurement error. Computational Statistics and Data Analysis 68 375–387.
  * Pollice et al. (2019) [author] Pollice, AlessioA., Jona Lasinio, GiovannaG., Rossi, RobertaR., Amato, MarianaM., Kneib, ThomasT. and Lang, StefanS. (2019). Bayesian measurement error correction in structured additive distributional regression with an application to the analysis of sensor data on soil–plant variability. Stochastic Environmental Research and Risk Assessment 33 747–763.
  * Polson, Scott and Windle (2013) [author] Polson, Nicholas G.N. G., Scott, James G.J. G. and Windle, JesseJ. (2013). Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables. Journal of the American Statistical Association 108 1339–1349.
  * Rajaratnam and Sparks (2015) [author] Rajaratnam, BalaB. and Sparks, DougD. (2015). MCMC-Based Inference in the Era of Big Data: A Fundamental Analysis of the Convergence Complexity of High-Dimensional Chains. preprint arXiv:1508.00947.
  * Richardson S (1993) [author] Richardson S, Gilks WRG. W. (1993). A Bayesian approach to measurement error problems in epidemiology using conditional independence models. American journal of epidemiology 138 430–42.
  * Roberts and Rosenthal (2001) [author] Roberts, Gareth OG. O. and Rosenthal, Jeffrey SJ. S. (2001). Markov chains and de-initializing processes. Scandinavian Journal of Statistics 28 489–504.
  * Rodrigues and Bolfarine (2007) [author] Rodrigues, JosemarJ. and Bolfarine, HelenoH. (2007). Bayesian inference for an extended simple regression measurement error model using skewed priors. Bayesian Analysis 61 349–364.
  * Stefanski (2000) [author] Stefanski, Leonard AL. A. (2000). Measurement error models. Journal of the American Statistical Association 95 1353–1358.
  * Stefanski and Carroll (1985) [author] Stefanski, Leonard A.L. A. and Carroll, Raymond J.R. J. (1985). Covariate Measurement Error in Logistic Regression. The Annals of Statistics 13 1335 – 1351.
  * Szegedy et al. (2014) [author] Szegedy, ChristianC., Zaremba, WojciechW., Sutskever, IlyaI., Bruna, JoanJ., Erhan, DumitruD., Goodfellow, IanI. and Fergus, RobR. (2014). Intriguing properties of neural networks. International Conference on Learning Representations.
  * Tang, Li and Tang (2017) [author] Tang, Nian-ShengN.-S., Li, De-WangD.-W. and Tang, An-MinA.-M. (2017). Semiparametric Bayesian inference on generalized linear measurement error models. Statistical Papers 58 1091–1113.
  * Torabi et al. (2021) [author] Torabi, MahmoudM., Ghosh, MalayM., Myung, JiyounJ. and Steel, MarkM. (2021). Measurement error in linear regression models with fat tails and skewed errors. Communications in Statistics - Theory and Methods 0 1–20.
  * Van Rossum and Drake Jr (1995) [author] Van Rossum, GuidoG. and Drake Jr, Fred LF. L. (1995). Python reference manual. Centrum voor Wiskunde en Informatica Amsterdam.
  * Vats, Flegal and Jones (2019a) [author] Vats, DootikaD., Flegal, James M.J. M. and Jones, Galin L.G. L. (2019a). Multivariate output analysis for Markov chain Monte Carlo. Biometrika 106 321–337.
  * Vats, Flegal and Jones (2019b) [author] Vats, DootikaD., Flegal, James MJ. M. and Jones, Galin LG. L. (2019b). Multivariate output analysis for Markov chain Monte Carlo. Biometrika 106 321-337. 10.1093/biomet/asz002
  * Vidal and Arellano-Valle (2010) [author] Vidal, IgnacioI. and Arellano-Valle, Reinaldo B.R. B. (2010). Bayesian inference for dependent elliptical measurement error models. Journal of Multivariate Analysis 101 2587–2597.
  * Vidal and Iglesias (2008) [author] Vidal, IgnacioI. and Iglesias, PilarP. (2008). Comparison between a measurement error model and a linear model without measurement error. Computational Statistics and Data Analysis 53 92–102.
  * Wang and Roy (2018) [author] Wang, XinX. and Roy, VivekanandaV. (2018). Geometric ergodicity of Pólya-Gamma Gibbs sampler for Bayesian logistic regression with a flat prior. Electronic Journal of Statistics 12 3295–3311.


# Robust Beamforming Design for Rate Splitting Multiple Access-Aided MISO
Visible Light Communications

Shuai Ma†, Guanjie Zhang†, Zhi Zhang†, and Rongyan Gu†
† Sch. of Inform. and Cont. Eng., China University of Mining and Technology,
Xuzhou 221116

###### Abstract

In this paper, we focus on the optimal beamformer design for rate splitting
multiple access (RSMA)-aided multiple-input single-output (MISO) visible light
communication (VLC) networks. First, we derive the closed-form lower bounds of
the achievable rate of each user, which are the first theoretical bound of
achievable rate for RSMA-aided VLC networks. Second, we investigate the
optimal beamformer design for RSMA-aided VLC networks to maximize the sum rate
under the optical and electrical power constraints. In addition, we show that
the proposed RSMA-aided networks can achieve superior performance compared
with space-division multiple access (SDMA) and non-orthogonal multiple access
(NOMA).

###### Index Terms:

Visible light communication, beamformer design.

##  I Introduction

By utilizing the deployed light emitting diode (LED) as the transmitter, VLC
can provide illumination and high-speed communication simultaneously, without
introducing interference to radio-frequency (RF) communications[2, 1]. Rate-
splitting multiple access (RSMA) was proposed [3], which performs linearly
precoded rate splitting at the transmitter and successive interference
cancellation (SIC) at the receivers. RSMA is a promising solution for VLC
networks. Two existing works have initiated the study of the RSMA design in
VLC networks [4, 5]. However, the achievable rate of RSMA-aided VLC networks
is still unknown due to the distinct characteristics of VLC. To be specific,
both peak and average optical power in VLC networks should be limited for eye
safety and practical illumination considerations [6, 7].

Motivated by the limitations of existing works on RSMA-aided VLC networks[4,
5, 3, 6, 7], we focus on the two key fundamental issues of RSMA-aided VLC
networks: identifying achievable rates and developing both optimal and robust
beamformer design. Specifically, by considering practical power constraints,
i.e., peak optical power constraint, average optical power constraint, and
average electrical power constraint, we derive the lower bound of achievable
rates of users in RSMA-aided VLC networks. Based on the derived lower bound of
achievable rates, we investigate the optimal beamformer design for RSMA-aided
VLC networks to maximize the sum rate under the optical and electrical power
constraints of LEDs.Simulation results show that the proposed algorithms of
RSMA-aided VLC networks can achieve superior performance compared with several
baseline schemes.

##  II System Model of Downlink RSMA-aided VLC Network

For the considered downlink RSMA-aided VLC network, the VLC base station (BS)
equipped with $N$ LEDs simultaneously serves $K$ single-access point users by
adopting RSMA scheme. Specifically, at BS, the message ${\text{M}_{k}}$
intended to user-$k$ is split into a common message $\text{M}_{k,0}$ and a
private message $\text{M}_{k}^{p}$, $\forall k\in{\cal{K}}$. Then, all the
common messages of $K$ users $\left\\{{{{\rm{M}}_{k,0}}}\right\\}_{k=1}^{K}$
are combined into a super common message ${\text{M}_{0}}$, i.e.,
${{\text{M}}_{0}}{\rm{}}\buildrel\Delta\over{=}\left\\{{{{\rm{M}}_{1,0}},\ldots,{{\rm{M}}_{K,0}}}\right\\}$.
In addition, these $K+1$ signals $\left\\{{{s_{i}}}\right\\}_{i=0}^{K}$
satisfy $\left|{{s_{i}}}\right|\leq{A_{i}}$,
$\mathbb{E}\left\\{{{s_{i}}}\right\\}=0$ and
$\mathbb{E}\left\\{{s_{i}^{2}}\right\\}={\varepsilon_{i}}$, $\forall
i\in{\cal{I}}$, where $A_{i}>0$ denotes the signal amplitude. Therefore, the
transmitted signal vector
${\bf{x}}={\left[{{x_{1}},\ldots,{x_{N}}}\right]^{T}}$ of VLC BS is given by

$\displaystyle{\bf{x}}=\sum\limits_{i=0}^{K}{{{\bf{w}}_{i}}{s_{i}}+}{\bf{b}}.$
(1)

The average electrical power and optical power of the transmitted signal
${\bf{x}}$ are respectively given by

$\displaystyle{P_{\rm{e}}}={\mathbb{E}}\left\\{{{{\left\|{\bf{x}}\right\|}^{2}}}\right\\}=\sum\limits_{i=0}^{K}{{\varepsilon_{i}}{{\left\|{{{\bf{w}}_{i}}}\right\|}^{2}}+N{b^{2}}},$
(2a) $\displaystyle
P_{\rm{o}}^{{\rm{ave}}}={\mathbb{E}}\left\\{{\bf{x}}\right\\}=Nb.$ (2b)

Let $I_{\rm{H}}$ denotes the maximum permissible current of LEDs, the
beamformer ${{\bf{w}}_{i}}$ should satisfy

$\displaystyle\sum\limits_{i=0}^{K}{{A_{i}}{\bf{w}}_{i}^{T}{{\bf{e}}_{n}}+b\leq{I_{\rm{H}}},\forall
n\in{\cal{N}}},$ (3)

where ${\bf{e}}_{n}$ is an unit vector with the $n$th element equal to 1.

The optical channel between LED and user is dominated by line-of-sight (LOS)
link, while diffuse links can be neglected [8, 9]. Specifically, let
${{\bf{g}}_{k}}\buildrel\Delta\over{=}{\left[{{g_{k,1}},\ldots,{g_{k,N}}}\right]^{T}}$
denote the channel gain vector of user-$k$, where ${g_{k,n}}$ is the channel
gain between the $n$th LED and the user-$k$. Therefore, at user-$k$, the
received signal $y_{k}$ is given as

$\displaystyle{y_{k}}=\underbrace{{\bf{g}}_{k}^{T}{{\bf{w}}_{0}}{s_{0}}}_{{\rm{common}}\;{\rm{signal}}}$
$\displaystyle+\underbrace{{\bf{g}}_{k}^{T}{{\bf{w}}_{k}}{s_{k}}}_{{\rm{private}}\;{\rm{signal}}}$
$\displaystyle+\underbrace{\sum\limits_{i=1,i\neq
k}^{K}{{\bf{g}}_{k}^{T}{{\bf{w}}_{i}}{s_{i}}}}_{{\rm{interference}}}+{\bf{g}}_{k}^{T}{\bf{b}}+{z_{k}},$
(4)

where ${z_{k}}$ is the received noise which follows the Gaussian distribution
with mean zero and covariance $\sigma_{k}^{2}$. The term
${\bf{g}}_{k}^{T}{\bf{b}}$ denotes the DC component, which can be removed by
the capacitor.

Thus, the residual received signal of user-$k$ after SIC process can be
represented

$\displaystyle
y_{k}^{{\rm{SIC}}}=\underbrace{{\bf{g}}_{k}^{T}{{\bf{w}}_{k}}{s_{k}}}_{{\rm{desired}}\;{\rm{private}}\;{\rm{signal}}}+\underbrace{\sum\limits_{i=1,i\neq
k}^{K}{{\bf{g}}_{k}^{T}{{\bf{w}}_{i}}{s_{i}}}}_{{\rm{interference}}}+{\bf{g}}_{k}^{T}{\bf{b}}+{z_{k}}.$
(5)

##  III DOWNLINK RSMA-aided VLC NETWORK

Let ${R_{k,c}}$ denotes the achievable rate of decoding common signal
${{s_{0}}}$, and its lower bound is given by

$\displaystyle{R_{k,{\rm{c}}}}\geq$
$\displaystyle\frac{1}{2}{\log_{2}}\left({\frac{{2\pi\sigma_{k}^{2}+\sum\limits_{i=0}^{K}{{{\left|{{\bf{g}}_{k}^{T}{{\bf{w}}_{i}}}\right|}^{2}}{e^{1+2\left({{\alpha_{i}}+{\gamma_{i}}{\varepsilon_{i}}}\right)}}}}}{{2\pi\sigma_{k}^{2}+2\pi\sum\limits_{j=1}^{K}{{{\left|{{\bf{g}}_{k}^{T}{{\bf{w}}_{j}}}\right|}^{2}}{\varepsilon_{j}}}}}}\right).$
(6)

Let ${R_{k,{\rm{p}}}}$ denote the achievable rate of decoding private signal
$s_{k}$, and its lower bound is given by

$\displaystyle{R_{k,{\rm{p}}}}\geq$
$\displaystyle\frac{1}{2}{\log_{2}}\left({\frac{{2\pi\sigma_{k}^{2}+\sum\limits_{i=1}^{K}{{{\left|{{\bf{g}}_{k}^{T}{{\bf{w}}_{i}}}\right|}^{2}}{e^{1+2\left({{\alpha_{i}}+{\gamma_{i}}{\varepsilon_{i}}}\right)}}}}}{{2\pi\sigma_{k}^{2}+2\pi\sum\limits_{j=1,j\neq
k}^{K}{{{\left|{{\bf{g}}_{k}^{T}{{\bf{w}}_{j}}}\right|}^{2}}{\varepsilon_{j}}}}}}\right).$
(7)

To ensure that ${{s_{0}}}$ is successfully decoded by all users, the
transmission rate of common message should not exceed
${R_{\rm{c}}}=\min\left\\{{{R_{1,{\rm{c}}}},\ldots,{R_{K,{\rm{c}}}}}\right\\}$.
Let ${c_{k}}\geq 0$ denotes one portion of achievable rate of common message
of user-$k$, where $\sum\nolimits_{k=1}^{K}{{c_{k}}={R_{\rm{c}}}}$ [3].
Therefore, the achievable rate of user-$k$, denoted as $R_{k}$, can be
expressed by

$\displaystyle R_{k}=c_{k}+{R_{k,{\rm{p}}}},~{}\forall k\in\cal{K}.$ (8)

###  III-A Optimal Beamformer Design

Mathematically, the sum rate maximization problem of RSMA-aided VLC networks
can be formulated as the following optimization problem

$\displaystyle\mathop{\max}\limits_{\scriptstyle\left\\{{{{\bf{w}}_{i}}}\right\\}_{i=0}^{K}\hfill\atop\scriptstyle\left\\{{{c_{k}}}\right\\}_{k=1}^{K}\hfill}$
$\displaystyle\sum\limits_{k=1}^{K}{{R_{k}}}$ (9a)
$\displaystyle{\rm{s}}{\rm{.t.}}$
$\displaystyle\sum\limits_{k=1}^{K}{{c_{k}}\leq{R_{k,{\rm{c}}}}},~{}{c_{k}}\geq
0,\forall k\in\cal{K},$ (9b)
$\displaystyle\sum\limits_{i=0}^{K}{{{\left\|{{{\bf{w}}_{i}}}\right\|}^{2}}}{\varepsilon_{i}}\leq{P_{\rm{t}}},$
(9c)
$\displaystyle\sum\limits_{i=0}^{K}{{A_{i}}}{\bf{w}}_{i}^{T}{{\bf{e}}_{n}}\leq\min\left\\{{b,{I_{\rm{H}}}-b}\right\\},\forall
n\in{\cal{N}}.$ (9d)

Since the objective function is non-convex, problem (9) is computationally
intractable. To transform problem (9) into tractable and equivalent problem,
we first introduce a number of auxiliary variables as follows

$\displaystyle{\widehat{\bf{w}}}\buildrel\Delta\over{=}{\left[{{\bf{w}}_{{}_{0}}^{\rm{T}},\cdots,{\bf{w}}_{{}_{K}}^{\rm{T}}}\right]^{\rm{T}}},$
(10a)
$\displaystyle\widehat{\bf{d}}\buildrel\Delta\over{=}{\left[{\varepsilon_{{}_{0}}^{1/2},\cdots,\varepsilon_{{}_{K}}^{1/2}}\right]^{\rm{T}}}\otimes{{\bf{1}}_{N}},$
(10b)
$\displaystyle{{\bf{a}}_{n}}\buildrel\Delta\over{=}{\left[{{A_{0}},\ldots{A_{K}}}\right]^{T}}\otimes{{\bf{e}}_{n}},$
(10c)
$\displaystyle{{\bf{G}}_{{\rm{c}},k}}\buildrel\Delta\over{=}{\rm{diag}}\left\\{{{\tau_{0}},\cdots,{\tau_{K}}}\right\\}\otimes\left({{{\bf{g}}_{k}}{\bf{g}}_{k}^{T}}\right),$
(10d)
$\displaystyle{\overline{\bf{G}}_{{\rm{c}},k}}{\rm{}}\buildrel\Delta\over{=}2\pi{\rm{diag}}\left\\{{0,{\varepsilon_{1}},\cdots,{\varepsilon_{K}}}\right\\}\otimes\left({{{\bf{g}}_{k}}{\bf{g}}_{k}^{T}}\right),$
(10e)
$\displaystyle{\widehat{\bf{G}}_{{\rm{c}},k}}\buildrel\Delta\over{=}{\rm{diag}}\left\\{{0,{\tau_{1}},\cdots,{\tau_{K}}}\right\\}\otimes\left({{{\bf{g}}_{k}}{\bf{g}}_{k}^{T}}\right),$
(10f)
$\displaystyle{\widehat{\bf{G}}_{{\rm{p}},k}}\buildrel\Delta\over{=}2\pi{\rm{diag}}\left\\{{0,{\varepsilon_{1}},\cdots{\varepsilon_{i-1}},0,{\varepsilon_{i+1}},\cdots,{\varepsilon_{K}}}\right\\}\otimes\left({{{\bf{g}}_{k}}{\bf{g}}_{k}^{T}}\right),$
(10g)

where ${{\bf{1}}_{N}}$ is a $N\times 1$ vector with all the elements equal
$1$, and
${\tau_{i}}={e^{1+2\left({{\alpha_{i}}+{\gamma_{i}}{\varepsilon_{i}}}\right)}},\forall
i\in{\cal{I}}$.

Based on the definitions (10), problem (9) can be equivalently reformulated as
follows

$\displaystyle\mathop{\max}\limits_{\scriptstyle\left\\{{{t_{k}},{c_{k}}}\right\\}_{k=1}^{K}\hfill\atop\scriptstyle\;\;\;\;{\bf{\hat{w}}}\hfill}$
$\displaystyle\sum\limits_{k=1}^{K}{{t_{k}}}$ (11a)
$\displaystyle{\rm{s}}.{\rm{t}}.$
$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{{\widehat{\bf{w}}}^{T}}{{\widehat{\bf{G}}}_{{\rm{c}},k}}\widehat{\bf{w}}}\right)$
$\displaystyle-\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{{\widehat{\bf{w}}}^{T}}{{\widehat{\bf{G}}}_{{\rm{p}},k}}\widehat{\bf{w}}}\right)\geq{t_{k}}-{c_{k}},$
(11b)
$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{{\widehat{\bf{w}}}^{T}}{{\bf{G}}_{{\rm{c}},k}}\widehat{\bf{w}}}\right)$
$\displaystyle-\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{{\widehat{\bf{w}}}^{T}}{{\overline{\bf{G}}}_{{\rm{c}},k}}\widehat{\bf{w}}}\right)\geq\sum\limits_{k=1}^{K}{{c_{k}}},$
(11c) $\displaystyle{c_{k}}\geq 0,\forall k\in{\cal K},$ (11d)
$\displaystyle{\left\|{{\bf{\hat{w}}}\odot\left({\widehat{\bf{d}}{{\widehat{\bf{d}}}^{T}}}\right)}\right\|^{2}}\leq{P_{\rm{t}}},$
(11e)
$\displaystyle{\widehat{\bf{w}}^{T}}{{\bf{a}}_{n}}{\bf{a}}_{n}^{T}\widehat{\bf{w}}\leq\min\left\\{{b^{2},{{\left({{I_{\rm{H}}}-b}\right)}^{2}}}\right\\},\forall
n\in{\cal{N}},$ (11f)

where $t_{k}$ is an auxiliary variable, $\forall k\in{\cal K}$.

Then, we employ the SDR technique to relax the constraints of problem (11).
Constraints (11b) and (11c) can be respectively rewritten as follows

$\displaystyle\underbrace{\frac{1}{2}{{\log}_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\widehat{\bf{G}}}_{{\rm{c}},k}}{\rm{)}}}\right)}_{{\rm{concave}}}$
$\displaystyle-\underbrace{\frac{1}{2}{{\log}_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\widehat{\bf{G}}}_{{\rm{p}},k}}{\rm{)}}}\right)}_{{\rm{concave}}}\geq{t_{k}}-{c_{k}},$
(12a)
$\displaystyle\underbrace{\frac{1}{2}{{\log}_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\bf{G}}_{{\rm{c}},k}}{\rm{)}}}\right)}_{{\rm{concave}}}$
$\displaystyle-\underbrace{\frac{1}{2}{{\log}_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\overline{\bf{G}}}_{{\rm{c}},k}}{\rm{)}}}\right)}_{{\rm{concave}}}\geq\sum\nolimits_{k=1}^{K}{{c_{k}}}.$
(12b)

Then, by using the first-order Taylor series expansion, we have

$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\widehat{\bf{G}}}_{{\rm{c}},k}}{\rm{)}}}\right)-{L_{{\rm{p}},k}}\left({\widehat{\bf{W}}}\right)\geq{t_{k}}-{c_{k}},$
(13a)
$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}\widehat{\bf{W}}{{\bf{G}}_{{\rm{c}},k}}{\rm{)}}}\right)-{L_{{\rm{c}},k}}\left({\widehat{\bf{W}}}\right)\geq\sum\limits_{k=1}^{K}{{c_{k}}},$
(13b)

where ${L_{{\rm{p}},k}}$ and ${L_{{\rm{c}},k}}$ are linear functions of
variable ${\widehat{\bf{W}}}$, which are given by

$\displaystyle{L_{{\rm{p}},k}}\left({\widehat{\bf{W}}}\right)=$
$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}{{\widehat{\bf{W}}}^{\left[m\right]}}{{\widehat{\bf{G}}}_{{\rm{p}},k}}{\rm{)}}}\right)$
$\displaystyle+\frac{{{\rm{Tr}}\left({{{\widehat{\bf{G}}}_{{\rm{p}},k}}(\widehat{\bf{W}}-{{\widehat{\bf{W}}}^{\left[m\right]}})}\right)}}{{\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}{{\widehat{\bf{W}}}^{\left[m\right]}}{{\widehat{\bf{G}}}_{{\rm{p}},k}}{\rm{)}}}\right)2\ln
2}},$ (14a) $\displaystyle{L_{{\rm{c}},k}}\left({\widehat{\bf{W}}}\right)=$
$\displaystyle\frac{1}{2}{\log_{2}}\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}{{\widehat{\bf{W}}}^{\left[m\right]}}{{\overline{\bf{G}}}_{{\rm{c}},k}}{\rm{)}}}\right)$
$\displaystyle+\frac{{{\rm{Tr}}\left({{{\overline{\bf{G}}}_{{\rm{c}},k}}(\widehat{\bf{W}}-{{\widehat{\bf{W}}}^{\left[m\right]}})}\right)}}{{\left({2\pi\sigma_{k}^{2}+{\rm{Tr(}}{{\widehat{\bf{W}}}^{\left[m\right]}}{{\overline{\bf{G}}}_{{\rm{c}},k}}{\rm{)}}}\right)2\ln
2}},$ (14b)

where ${\widehat{\bf{W}}^{\left[m\right]}}$ is a feasible point obtained from
the $m$th iteration.

Therefore, at the $\left({m+1}\right)$th iteration, the convex approximation
form of problem (11) is given as

$\displaystyle\mathop{\max}\limits_{\scriptstyle\left\\{{{t_{k}},{c_{k}}}\right\\}_{k=1}^{K}\hfill\atop\scriptstyle\;\;\;\widehat{\bf{W}}\hfill}$
$\displaystyle\sum\limits_{k=1}^{K}{{t_{k}}}$ (15a) $\displaystyle\rm{s.t.}$
$\displaystyle\eqref{DC_form1_a},~{}\eqref{DC_form1_b}$
$\displaystyle{\rm{Tr}}\left({\widehat{\bf{W}}\odot(\widehat{\bf{d}}{{\widehat{\bf{d}}}^{T}})}\right)\leq{P_{\rm{t}}},$
(15b)
$\displaystyle{\rm{Tr}}\left({\widehat{\bf{W}}{{\bf{a}}_{n}}{\bf{a}}_{n}^{T}}\right)\leq\min\left\\{{b^{2},{{\left({{I_{\rm{H}}}-b}\right)}^{2}}}\right\\},\forall
n\in{\cal{N}},$ (15c)
$\displaystyle\widehat{\bf{W}}\succeq{\bf{0}},{{c}_{k}}\geq 0,\forall
k\in{\cal{K}}.$ (15d)

##  IV SIMULATION RESULTS AND DISCUSSION

Consider a RSMA-aided multi-user VLC network deployed in a room with a size of
$7\times 7\times 5{{\rm{m}}^{3}}$, where the unit of distances is meter.

Fig. 1 (a) shows the sum rate (bits/sec/Hz) of $U_{1}$ and $U_{2}$ versus the
transmitted power budget $P_{\rm{t}}$ (dB) of the RSMA, NOMA and SDMA schemes.
It can be observed that the sum rate of the three schemes increases as
$P_{\rm{t}}$ increases, and the sum rate of RSMA is higher than that of NOMA
and SDMA schemes. Furthermore, Fig. 1 (b) depicts the sum rate of two users
$U_{3}$ and $U_{4}$ for different transmitted power budget $P_{\rm{t}}$ of the
RSMA, NOMA and SDMA schemes. From Fig. 1 (b), we observe that the RSMA scheme
always outperforms both the NOMA and SDMA schemes, which is similar to that in
Fig. 1 (a).

(a)

(b)

Figure 1: (a) The sum rate (bits/sec/Hz) versus ${P_{\rm{t}}}({\rm{dB}})$ with
perfect CSIT , where $N=4$, $K=2$ uses are located at $U_{1}$ and $U_{2}$; (b)
The sum rate $\left({{\rm{bit/sec/Hz}}}\right)$ versus ${P_{\rm{t}}}$
$\left({{\rm{dB}}}\right)$ with perfect CSIT, where $N=4$, $K=2$ uses are
located at $U_{3}$ and $U_{4}$;

##  V Conclusion

In this work, we addressed optimal and robust beamformer design allocation for
RSMA-aided VLC networks. Specifically, we derived the lower bound of
achievable rates of RSMA-aided VLC networks. Moreover, based on the derived
bound of achievable rates, we investigated the optimal beamformer design to
maximize the sum rate under the optical and electrical power constraints of
LEDs. Simulation results illustrated the superiority of the proposed
beamformer design and provide useful insights on the design of RSMA-aided VLC
networks.

## References

  * [1] P. H. Pathak, X. Feng, P. Hu, and P. Mohapatra,  “Visible light communication, networking, and sensing: A survey, potential and challenges,”  IEEE Commun. Surveys Tuts., vol. 17, no. 4, pp. 2047–2077, Sep. 2015.
  * [2] A. Jovicic, J. Li, and T. Richardson,  “Visible light communication: opportunities, challenges and the path to market,”  IEEE Commun. Mag., vol. 51, no. 12, pp. 26–32, Dec. 2013.
  * [3] Y. Mao, B. Clerckx, and V. Li,  “Rate-splitting multiple access for downlink communication systems: Bridging, generalizing, and outperforming SDMA and NOMA,”  EURASIP J. Wireless Commun. and Network., pp. 1–54, May. 2018.
  * [4] S. Tao, H. Yu, Q. Li, Y. Tang, and D. Zhang,  “One-layer rate-splitting multiple access with benefits over power-domain NOMA in indoor multi-cell visible light communication networks,”  in Proc. IEEE/CIC Int. Conf. Commun. China, pp. 1–7, Jun. 2020\.
  * [5] S. Naser, P. C. Sofotasios, L. Bariah, W. Jaafar, S. Muhaidat, M. Al-Qutayri, and O. A. Dobre,  “Rate-splitting multiple access: Unifying NOMA and SDMA in MISO VLC channels,”  IEEE Open J. Veh. Tech., vol. 1, pp. 393–413, Oct. 2020.
  * [6] A. Jovicic, J. Li, and T. Richardson,  “Visible light communication:Opportunities, challenges and the path to market,”  IEEE Commun. Mag., vol. 51, no. 12, pp. 26–32, Dec. 2013.
  * [7] P. H. Pathak, X. Feng, P. Hu, and P. Mohapatra,  “Visible light communication, networking, and sensing: A survey, potential and challenges,”  IEEE Commun. Surveys Tuts., vol. 17, no. 4, pp. 2047–2077, Sep. 2015.
  * [8] T. Fath and H. Haas,  “Performance comparison of MIMO techniques for optical wireless communications in indoor environments,”  IEEE Trans.Commun, vol. 61, no. 2, pp. 733–742, Feb. 2013.
  * [9] T. Q. Wang, Y. A. Sekercioglu, and J. Armstrong,  “Analysis of an optical wireless receiver using a hemispherical lens with application in MIMO visible light communications,”  J. Lightw. Technol., vol. 31, no. 11, pp. 1744–1754, Jun. 2013\.