As I understand it the function (of two symmetric positive definite matrices) is

```
JL(X,Y) = log( det( (X+Y)/2)) - log( det( X*Y))/2
= log( det( (X+Y)/2)) - (log(det(X)) + log(det(Y)))/2
```

I reckon the way to do this is to compute the cholesky factorisations of X, Y and (X+Y)/2, that is find lower triangular matrices L,M,N so that

```
X = L*L'
Y = M*M'
(X+Y)/2 = N*N'
```

(I’m sorry to say I don’t know matlab, but there is bound to be a function to compute the cholesky factorisation).

and then, using the properties of determinants and logs we have

```
JL(X,Y) = 2*log( det(N)) - log( det(L)) - log(det(M)))
```

For a lower triangular matrix K we have

```
log( det( K)) = log ( Prod( K[i,i], i=1..n))
= Sum ( log( K[i,i]), i=1..n))
```

I think it might be better to use the second form, even though it is slower, as there could be a risk of underflow in the product if the matrices are large, and perhaps the best thing to do would be to compute

```
JL(X,Y) = Sum ( 2*log( N[i,i]) - log( L[i,i]) - log( M[i,i]), i=1..n)
= Sum( log( N[i,i]*N[i,i] / (L[i,i]*M[i,i])), i=1..n)
```

Note that all of this assumes that the cholesky factorisations can be computed, i.e. that the matrices X and Y are strictly positive definite, but then the function JL is only defined in that case.

solved Jensen-Bregman LogDet divergence [closed]