Update Magnitude Function to include Bastian's solver

[lizliz]

Pasting in the original message from @Pseudomanifold Bastien Rieck in here. The idea is to modify the magnitude function computation to do this version as an option.

---

Sorry, don't know how to properly add to this request...but here's some code that calculates magnitude using a Cholesky decomposition.

from scipy.linalg import cho_factor
from scipy.linalg import solve_triangular

[...]

M = np.exp(-D)
c, lower = cho_factor(M)
x = solve_triangular(c, np.ones(M.shape[0]), trans=1)

return x.T @ x

The idea is that the similarity matrix $\zeta$ is positive definite (under nice conditions), so it affords a Cholesky decomposition, which factorises as $L L^T$ for $L$ a lower-triangular matrix. Magnitude can be written as $\textrm{Mag}(X) = \mathbb{1}^T \zeta^{-1} \mathbb{1}$, i.e. a quadratic form (since we are taking row sums). Thus, calculating magnitude is equivalent to calculating $x^T x$ with $x$ defined by $Lx = \mathbb{1}$.

(We have described this a little bit in Appendix B of a recent preprint)

Hope that helps!

Originally posted by @Pseudomanifold in #89 (comment)