Note issue with formula

[maresb]

The following formula doesn't hold in general. For that you'd need matrices to be simultaneously diagonalizable.

[aseyboldt]

You are right... I guess I'm lucky I never used that one for anything. The square root expression is right though:

def geomean(A, B):
    A_sqrt = linalg.sqrtm(A)
    A_isqrt = linalg.inv(A_sqrt)
    return A_sqrt @ linalg.sqrtm(A_isqrt @ linalg.inv(B) @ A_isqrt) @ A_sqrt

np.allclose(geomean(covx, covalpha), C)

This is where I compute it in nutpie: https://github.com/pymc-devs/nuts-rs/blob/main/src/low_rank_mass_matrix.rs#L409 Please don't find a mistake in there ;-)

[maresb]

It turns out you can solve this.

If $A=\mathrm{Cov}[\alpha_i]$ and $B=\mathrm{Cov}[x_i]$, then

$$\Sigma = A^{-1/2} \left( A^{1/2} B A^{1/2} \right)^{1/2} A^{-1/2}$$

[aseyboldt]

That was the "or some expression with lots of matrix square roots (todo)." in the paper :D Cutting corners always gets punished

[maresb]

Where is geomean from?

[aseyboldt]

I think it was this: https://www.lix.polytechnique.fr/~nielsen/MIGBOOKWEB/9783642302312-c2.pdf

[maresb]

I mean the code for geomean. Is that copy-pasted from somewhere?

[aseyboldt]

No, I just wrote the code. Same solution you just posted.