# Inverting big symmetric and singular matrices

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

• In what way was that other Q/A unsatisfactory? Upon review I think perhaps my answer only told you things you'd already read Nov 7, 2017 at 17:01

Think of your matrix as block-diagonal with blocks $$C_\ell = \begin{pmatrix} 0 & & \\ & D_\ell & \\ & & 0\end{pmatrix}.$$ Then it is clear that $$C_\ell^{-1} = \begin{pmatrix} 0^{-1} & & \\ & D_\ell^{-1} & \\ & & 0^{-1}\end{pmatrix}$$ where of course $0^{-1}$ is not well defined. But that doesn't matter because when you compute a term such as $$C_\ell^{-1}C_{\ell,\alpha}$$ you will get $$C_\ell^{-1}C_{\ell,\alpha} = \begin{pmatrix} 0^{-1}0 & & \\ & D_\ell^{-1}D_{\ell,\alpha} & \\ & & 0^{-1}0\end{pmatrix}$$ which shows that you only ever multiply the ill-defined terms with another zero matrix. So $$C_\ell^{-1}C_{\ell,\alpha} = \begin{pmatrix} 0 & & \\ & D_\ell^{-1}D_{\ell,\alpha} & \\ & & 0\end{pmatrix}$$ and similarly for the second product: $$C_\ell^{-1}C_{\ell,\alpha}C_\ell^{-1}C_{\ell,\beta} = \begin{pmatrix} 0 & & \\ & D_\ell^{-1}D_{\ell,\alpha}D_\ell^{-1}D_{\ell,\beta} & \\ & & 0\end{pmatrix}$$ and consequently $$\text{trace}\left(C_\ell^{-1}C_{\ell,\alpha}C_\ell^{-1}C_{\ell,\beta}\right) = \text{trace}\left(D_\ell^{-1}D_{\ell,\alpha}D_\ell^{-1}D_{\ell,\beta}\right).$$ This is perfectly well defined as long as $D_\ell$ is invertible.