# Is there a more efficient way to compute the log-density of a Multivariate normal with a Symmetric Woodbury Covariance?
**URL:** https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336
**Category:** General Usage
**Created:** [May 4, 2023, 5:09pm UTC](https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336 "2023-05-04T17:09:45Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@schwob](https://discourse.julialang.org/u/schwob)
#### Post date: [May 4, 2023, 5:09pm UTC](https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336/1 "2023-05-04T17:09:45Z")
Suppose that I have \mathbf{y}\sim N(\boldsymbol{\mu}, \mathbf{A} + \mathbf{UCU'}), where \mathbf{A} and \mathbf{C} are diagonal matrices. Fortunately, the covariance term follows the pattern of a sparse symmetric Woodbury matrix.
I coded up the following function to evaluate the log-density of a multivariate normal:
```julia
function lpdfMVN_Woodbury(y, mu, A, U, C)
k = length(mu)
Const = -k/2*log(2*pi)
invA = inv(A)
invAU = invA*U
inner_term = inv(C) + U'*invAU
log_determinant = logdet(inner_term) + logdet(C) + logdet(A)
wood_inv = invA - invAU*inv(inner_term)*(invAU)'
ymu = y - mu
return[Const - log_determinant/2 - ymu'*wood_inv*ymu/2]
end
```
While this is much quicker than `logpdf(MvNormal(mu, A + U*C*U'), y)`, it’s still a bit slower than I’d like it to be. Does anyone have any suggestions to speed up this evaluation?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [May 4, 2023, 7:19pm UTC](https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336/3 "2023-05-04T19:19:17Z")
Can you profile the code to see which lines take longer on a typical example?
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [May 4, 2023, 7:49pm UTC](https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336/4 "2023-05-04T19:49:31Z")
If they aren’t already, it may help if `A` and `C` are of the `Diagonal` type. Their inverses will take less memory and operations with them will be faster (no need to check for diagonality for `inv` and `logdet`, no need to add off-diagonals for `+`).
While `inv` and `logdet` on diagonal matrices are very cheap, they have a real cost on non-diagonal matrices like `inner_term`. I would replace a few lines with the following:
```julia
inner_term = cholesky(Hermitian(inv(C) + U'*invAU))
log_determinant = logdet(inner_term) + logdet(C) + logdet(A)
wood_inv = invA - invAU * (inner_term \ invAU') # maybe wrap in Hermitian(...)
```
`logdet` and `inv` are both expensive operations that become very easy once `inner_term` has been factorized, so calling `cholesky` ensures that this only has to be done once rather than twice. You may need to apply a pivoting scheme if you run into definiteness issues with `cholesky`, or consider a `bunchkaufman` factorization instead.
Stylistically, I might replace some of the other `inv(X)*Y` statements with `X \ Y`, but for diagonal `X` (and since you use the inverses in additions later) I don’t foresee this making a big difference. There are further games you could play with in-place functions for your matrix operations, but they make things a bit less readable and may not be that impactful. You can look into `LinearAlgebra.mul!` to that end.
---
### Author:  [@schwob](https://discourse.julialang.org/u/schwob)
#### Post date: [May 4, 2023, 9:12pm UTC](https://discourse.julialang.org/t/is-there-a-more-efficient-way-to-compute-the-log-density-of-a-multivariate-normal-with-a-symmetric-woodbury-covariance/98336/5 "2023-05-04T21:12:13Z")
Great idea. The two lines that take noticeably longer are the computations for `log_determinant` and `wood_inv`. However, these appear to be sped up using @mikmoore 's suggestion!