\begin{align}
\label{eqn:loss}
    L(\beta_1,...,\beta_p,\gamma_1,...,\gamma_p)     &= l(\beta_1,...,\beta_p,\gamma_1,...,\gamma_p;\textbf{Y}) + \text{Penalty} \nonumber \\
    &= \frac{1}{2}(\mathbf{Y}-\beta_0-\mathbf{X}\mathbf{\beta})^T\left(\sigma_0^2\mathbf{R} + \left(\sum_{i=1}^p\gamma_i\mathbf{X_iX_i^T}\right)\odot \mathbf{V} + \sum_{i=1}^p\gamma_iq_i\mathbf{X_iX_i^T}\right)^{-1} \nonumber \\&\qquad(\mathbf{Y}-\beta_0-\mathbf{X}\mathbf{\beta})+\frac{1}{2}\log{\det{\left(\sigma_0^2\mathbf{R} + \left(\sum_{i=1}^p\gamma_i\mathbf{X_iX_i^T}\right)\odot \mathbf{V} + \sum_{i=1}^p\gamma_iq_i\mathbf{X_iX_i^T}\right)}} \nonumber \\&\qquad+\lambda_1\|\beta\|_1+\lambda_2\|\gamma\|_1
\end{align}
