# A Novel Bayesian Linear Regression Model for the Analysis of Neuroimaging Data
## 3 Methods
### A Dual Bayesian Linear Regression Model with Feature Selection
#### 3.1.3 Variational Inference  
Once the generative model is defined, we should evaluate the posterior distribution of the variables to estimate their optimum values. Although, in this case, the posterior distribution is not tractable, we can use variational inference together with the mean-field technique [34] to find an approximation to this posterior \(q(\mathbf{\Theta})\approx p(\mathbf{\Theta}|\boldsymbol{y},\mathbf{X})\), where \(\mathbf{\Theta}\) contains all model variables. Then, we can define a Lower Bound (LB) using the Kullback  
Figure 2: Plate diagram for the DBL-FS graphic model. Grey circles denote observed variables, white circles unobserved variables. Model hyperparameters do not have a circle. Leibler divergence between the posterior and its approximation; so, maximizing this LB, we can obtain the optimum model parameters. Therefore, using the mean-field approximation to factorize over the posterior, we obtain:  
\[p(\mathbf{\Theta}|\boldsymbol{y},\mathbf{X})\approx q(\mathbf{\Theta})=q( \mathbf{w})q(\mathbf{a})q(\boldsymbol{\alpha})q(\tau), \tag{8}\]  
and we can determine each approximated distribution by calculating:  
\[\ln(q_{j}^{*})=\mathbb{E}_{-q_{j}}[\ln(p(\mathbf{X},\boldsymbol{y},\mathbf{ \Theta}))]+const, \tag{9}\]  
where \(\mathbb{E}_{-q_{j}}\) implies that we calculate the expectation over all random variables except the \(j\)-th variable, and \(p(\mathbf{X},\boldsymbol{y},\mathbf{\Theta})\) is the joint probability. Therefore, we can apply (9) to the joint probability for each random variable to obtain the model update rules. Firstly, the distribution of the dual weights \(\mathbf{a}\) is:  
\[q(\mathbf{a})=\mathcal{N}(\mathbf{a}|\langle\mathbf{a}\rangle,\boldsymbol{ \Sigma}_{a}), \tag{10}\]  
with mean and variance determined by:  
\[\langle\mathbf{a}\rangle=\langle\tau\rangle\boldsymbol{\Sigma}_{a}\mathbf{K}^ {\intercal}y \tag{11}\]  
\[\boldsymbol{\Sigma}_{a}^{-1}=\mathbf{X}\text{diag}(\langle\boldsymbol{\alpha }\rangle)\mathbf{X}^{T}+\langle\tau\rangle\mathbf{K}^{\intercal}\mathbf{K}, \tag{12}\]  
where \(\text{diag}(\langle\boldsymbol{\alpha}\rangle)\) represents an identity matrix with vector \(\boldsymbol{\alpha}\) as the diagonal. The distribution of variable \(\boldsymbol{\alpha}\) is:  
\[q(\boldsymbol{\alpha})=\ \Gamma(\boldsymbol{\alpha}|a_{\alpha},\mathbf{b}_{ \alpha}), \tag{13}\]  
with parameters  
\[a_{\alpha}=a_{0}^{\alpha}+\frac{D}{2} \tag{14}\]  
\[\mathbf{b}_{\alpha}=b_{0}^{\alpha}+\frac{1}{2}\text{diag}(\mathbf{X}^{T} \langle\mathbf{a}\mathbf{a}^{T}\rangle\mathbf{X}), \tag{15}\]  
where \(a_{0}^{\alpha}\) and \(\beta_{0}^{\alpha}\) are hyperparameters, and the operator diagonal returns a column vector formed by the main diagonal of the matrix. Moreover, the distribution of the noise precision \(\tau\) is given by:  
\[q(\tau)=\ \Gamma(\tau|a_{\tau},b_{\tau}), \tag{16}\]  
with parameters  
\[p(\mathbf{\Theta}|\boldsymbol{y},\mathbf{X})\approx q(\mathbf{\Theta})=q( \mathbf{w})q(\mathbf{a})q(\boldsymbol{\alpha})q(\tau), \tag{8}\]0  
\[p(\mathbf{\Theta}|\boldsymbol{y},\mathbf{X})\approx q(\mathbf{\Theta})=q( \mathbf{w})q(\mathbf{a})q(\boldsymbol{\alpha})q(\tau), \tag{8}\]1  
where \(a_{0}^{\tau}\) and \(\beta_{0}^{\tau}\) are hyperparameters, and \(Tr\{\}\) is the trace operator. See Appendix A for the full development of these mean field distribution approximations. Once we have defined the different distributions, the model updates the different random variables in an iterative coordinate-ascent-like optimization where the distribution of each factor is obtained using (10) to (22). This optimization process is guided by the \(LB\) cost function defined as:  
\[p(\mathbf{\Theta}|\boldsymbol{y},\mathbf{X})\approx q(\mathbf{\Theta})=q( \mathbf{w})q(\mathbf{a})q(\boldsymbol{\alpha})q(\tau), \tag{8}\]2  
where we analyze its convergence to stop the distribution parameters update. See Appendix B for the full development of the LB. For an efficient optimization of the model, in practice, we will work in the dual space updating the Equations (10) to (18).