# Matrix inverse and multivariate normal distribution

**URL:** <https://discourse.julialang.org/t/matrix-inverse-and-multivariate-normal-distribution/118724>\
**Category:** Statistics\
**Tags:** statistics, linearalgebra, optimization, math\
**Created:** [August 28, 2024, 4:47pm UTC](https://discourse.julialang.org/t/matrix-inverse-and-multivariate-normal-distribution/118724 "2024-08-28T16:47:01Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![fmor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fmor/32/211337_2.png) [@fmor](https://discourse.julialang.org/u/fmor)\
**Post date:** [August 28, 2024, 4:47pm UTC](https://discourse.julialang.org/t/matrix-inverse-and-multivariate-normal-distribution/118724/1 "2024-08-28T16:47:01Z")

</div>

Hello, I am fitting a spatio-temporal-and-informed-with-covariates bayesian model in Julia. I have a regression vector with prior distribution \vec{\beta}\_t \sim \mathcal{N}\_p(\vec{b\_0},s^2I) (where p is the number of covariates, I the identity matrix). At every fitting iteration it has (as the other parameters of the model have) to be updated. The updating rule derive from the full conditional derivation, which corresponds to the following computations:  
f(\vec{\beta}\_t|-) \sim \exp \left\{ \vec{\beta}\_t^T \underbrace{\left( \frac{1}{s^2}I + \sum\_i \frac{\vec{x}\_{it}\vec{x}\_{it}^T}{\sigma^2\_{it}}\right)}\_{A\_\star}\vec{\beta}\_t -2 \underbrace{\left( \frac{\vec{b\_0}}{s^2} + \sum\_i \frac{\text{stuff}\_{it} \vec{x}\_{it}}{\sigma^2\_{it}}\right)}\_{\vec{b\_\star}}\cdot \vec{\beta}\_t \right\}  
which is proportional to

\propto \exp \left\{ (\vec{\beta}\_t -\vec{b}\_\star)^T A\_\star (\vec{\beta}\_t -\vec{b}\_\star)\right\}  
which is the kernel of another multivariate gaussian with covariance matrix _the inverse_ (sadly) of A\_\star (you can see why the inverse from the picture).

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/a/eab786048d88da64ad350f1277632d9e5c9a1253.png)

That math translates into the following code:

```julia
sum_Y = zeros(p)
A_star = I(p)/s2_beta
for j in 1:n
	X_jt = @view Xlk_covariates[j,:,t]
	sum_Y += (Y[j,t] - muh_iter[Si_iter[j,t],t] - eta1_iter[j]*Y[j,t-1]) * X_jt / sig2h_iter[Si_iter[j,t],t]
	A_star += (X_jt * X_jt') / sig2h_iter[Si_iter[j,t],t]
end
b_star = beta0/s2_beta + sum_Y
Am1_star = inv(Symmetric(A_star))
# Symmetric needed to solve an inversion and hermitian problem
# https://discourse.julialang.org/t/isposdef-and-eigvals-do-not-agree/118191/10
beta_iter[t] = rand(MvNormal(Am1_star*b_star, Am1_star))

```

which would work fine but for the inversion of the matrix required, which appears to be a major problem. Here is for example of an A\_star matrix obtained during the fit and her inverse:

```julia
A_star = [89.53344398111983 18.907579544068415 -22.32917044622947 -26.313826364434874; 18.907579544068415 34.821884016362716 -18.103513531335015 -30.76131476315723; -22.32917044622947 -18.103513531335015 199.9111114786475 21.85563450347075; -26.313826364434874 -30.76131476315723 21.85563450347075 75.35950798333245]
# 89.5334 18.9076 -22.3292 -26.3138
# 18.9076 34.8219 -18.1035 -30.7613
# -22.3292 -18.1035 199.911 21.8556
# -26.3138 -30.7613 21.8556 75.3595
cond(A_star)
# 12.015947669881697, which doesn't seem too bad to me
A*inv(A_star) 
# 0.130423 -0.0446254 0.00778601 0.0250668
# -0.0446254 0.47336 0.0190656 0.172111
# 0.00778601 0.0190656 0.0531558 -0.00491498
# 0.0250668 0.172111 -0.00491498 0.21313

```

So I was looking for some paths to avoid that inversion.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [August 28, 2024, 5:07pm UTC](https://discourse.julialang.org/t/matrix-inverse-and-multivariate-normal-distribution/118724/2 "2024-08-28T17:07:00Z")

</div>

I don’t actually know much, but does this help:

> **[Precision (statistics)](https://en.wikipedia.org/wiki/Precision_(statistics))**
>
> In statistics, the precision matrix or concentration matrix is the matrix inverse of the covariance matrix or dispersion matrix, 
>   
>     
>       
> P
> =
>         
> Σ
>           
> −
> 1
>           
>         
>       
>     
> {\\displaystyle P=\\Sigma ^{-1}}
>   
> .
> For univariate distributions, the precision matrix degenerates into a scalar precision, defined as the reciprocal of the variance, 
>   
>     
>       
> p
> =
>         
>           
> Other su...

> One particular use of the precision matrix is in the context of [Bayesian analysis](https://en.wikipedia.org/wiki/Bayesian_analysis) of the [multivariate normal distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution): for example, Bernardo & Smith prefer to parameterise the multivariate normal distribution in terms of the precision matrix, rather than the covariance matrix, because of certain simplifications that then arise.[[10]](https://en.wikipedia.org/wiki/Precision_(statistics)#cite_note-10) For instance, if both the [prior](https://en.wikipedia.org/wiki/Prior_probability) and the [likelihood](https://en.wikipedia.org/wiki/Likelihood_function) have [Gaussian](https://en.wikipedia.org/wiki/Gaussian_function) form, and the precision matrix of both of these exist (because their covariance matrix is full rank and thus invertible), then the precision matrix of the [posterior](https://en.wikipedia.org/wiki/Posterior_probability) will simply be the sum of the precision matrices of the prior and the likelihood.

---

<div class="post-metadata">

**Author:** ![jkbest2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jkbest2/32/7350_2.png) [@jkbest2](https://discourse.julialang.org/u/jkbest2)\
**Post date:** [August 28, 2024, 6:07pm UTC](https://discourse.julialang.org/t/matrix-inverse-and-multivariate-normal-distribution/118724/3 "2024-08-28T18:07:09Z")

</div>

The [`MvNormalCanon` distribution](https://juliastats.org/Distributions.jl/stable/multivariate/#Distributions.MvNormalCanon) is parameterized in terms of the precision matrix (inverse covariance).