full conditional distribution for $\eta_i$ is: $$\begin{aligned}
& -\frac{1}{2} \Big[ \eta_i^T(\frac{\omega \omega^T}{\sigma^2_y}+\Lambda^T \Psi^{-1}\Lambda + I - 2 \frac{ \Omega Y_i }{\sigma_y^2})\eta_i  - 2 \eta_i ^T (\Lambda^T \Psi^{-1}X_i + \frac{\omega Y_i}{\sigma_y^2})\Big] - \\ 
& -\frac{1}{2} \Big[ \frac{2 \eta_i^T \omega \eta_i^T \Omega \eta_i}{\sigma_y^2} + \frac{(\eta_i^T \Omega \eta_i)^2}{\sigma_y^2}  \Big] + C,\end{aligned}$$ where $C$ is a normalizing constant. We update the factors $\eta_i$ using the Metropolis-Adjusted Langevin Algorithm (MALA) [@grenander1994representations], [@roberts1996exponential]. Sampling the factors is the main computational bottleneck of our approach since we have to update $n$ vectors, each of dimension $k$. The overall MCMC algorithm and the MALA step are summarized in *Algorithm 1*.

[ *Step 1* Sample $\eta_i$, $i = 1, \cdots,n$ via Metropolis-Hastings using as a proposal distribution a $N\big(\eta_i + \hspace*{0.9cm}\frac{1}{2}\nabla_{\eta_i } log(\pi(\eta_i | \hbox{ --- }) \big), \epsilon I_k)$. *Step 2* Sample the main effects coefficients $\omega$ from a multivariate normal distribution: $$\begin{aligned}
\pi(\omega |  \hbox{ --- }  ) \sim N \bigg( (\frac{\eta^T \eta}{\sigma^2} + I_n/100)^{-1} \eta \big(y - diag(\eta \Omega \eta)\big)/\sigma^2,(\frac{\eta^T \eta}{\sigma^2} + I_n/100)^{-1} \bigg)\end{aligned}$$\
where $\eta$ is the matrix with rows equal to $\eta_i$. *Step 3* Sample upper triangular part of $\Omega$, namely $\Omega^{\text{U}}$, from a multivariate normal distribution: $$\begin{aligned}
\pi(\Omega^{\text{U}} |  \hbox{ --- }  ) 