Efficient Computation of Matrix Determinants and Inversion

[Jammy2211]

We have an efficient data structure for our F matrix, where each value is stored in a dictionary. We can do the same for the H matrix trivially. Unfortunately, we currently need to convert this to a matrix to compute the determinant and solve the set of linear equations of these matrices.

Thus, we now want to use our data structure to perform the above tasks. We need the determinant of two matrices (H and F+H) and need to solve for one matrix (S, from D and F+H). For both tasks, we need to compute the Cholesky decomposition of our matrices, e.g. F+H = LL and H = LL, which represents them as the multiplication of a diagonal matrix (i.e. all zeros in the upper or lower diagonal) multiplied by its conjugate transpose. There are standard algorithms for doing this.

The determinant is now easy to compute, it is simply the sum of the natural log of all diagonal terms of L, multiplied by 2:

Det_F = 2.0* sum_i log (L(i,i))

The inversion requires more care, using forward and back substitution . I haven't yet wrapped my head around how exactly that works, but the c# code supplied in the link is only 30 lines long so it can't be too tricky (right?...).

[Jammy2211]

No doubt, the question might arise as to whether we want to compute F and H at all, or directly just compute their Cholesky decompositions.

[Jammy2211]

Well, numpy linear algrebra is a lot faster than I could of ever expected, so no need to do this.