In this post I found a very similar probem to the one I have, but not a satisfactory answer for my purposes.
I have a set of matrices $C_\ell$. They are exactly symmetric by construction. Unfortunately, they are also singular. Their structure is similar to a block matrix.
Here a "picture" of how they look like:
0000000000
0000000000
00------00
00------00
00------00
00------00
0000000000
0000000000
(and they are square, although ti doesn't seem like that from the picture). The dashed lines represent values different from 0.
What I need to do is to calculate the Fisher matrix from them. The formula for the Fisher matrix looks like this $$F_{\alpha \beta} = \sum_\ell \frac{1}{2} Tr [ C^{-1}_\ell C_{\ell,\alpha} C_\ell^{-1} C_{\ell,\beta}]$$ where $C_{\ell,\alpha}$ denotes derivative of the $C_\ell$ w.r.t the parameter $\alpha$.
So far I have been using a Singuar Value Decomposition to operate the inversion of $C_\ell$'s. But recently I discovered that, as the dimension of my $C_\ell$'s increases (now it has reached 900x900
) this method gives very imprecise results.
I am implementing everything in C++
. Any idea on how to perform better these inversions? Speed is not an issue (as long as the operation remains feasible in human times), I really care about precision