1
$\begingroup$

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

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

2 Answers 2

6
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ @johnhenry Now with this essential step form Wolfgang you can do Cholesky factorisation what is advised in rchilton1980 answer. $\endgroup$
    – likask
    Nov 8, 2017 at 11:35
3
$\begingroup$

You might be better served by either the LDL' decomposition or the Cholesky decomposition (in the event that C's are positive definite in addition to symmetric, they probably are). Though all of the algorithms have cubic complexity, the constants are much better for these "linear" solvers than the SVD (which is more of a "spectral" solver). They are also easier to parallelize (Cholesky is particularly easy). These decompositions are readily available in LAPACK or MKL. The respective factorization routines are DSYTRF (for LDL') and DPOTRF (for Cholesky). You probably don't want to form the inverse explicity, you can instead use the decomposition to apply the "action" of the inverses to your forcing data.

$\endgroup$
3
  • $\begingroup$ there is a problem with Cholesky decomposition. I cannot use it, because my matrix happens to be singular (that's why I am using the SVD in the first place). I am editing my post to take this into account, thanks for your input $\endgroup$
    – johnhenry
    Nov 7, 2017 at 18:53
  • $\begingroup$ So perhaps we are really talking about a pseudoinverse, then? I believe you could still use either LDL' or Cholesky to find the inverse of that "island" of nonzero values. Then you would zero-pad the triangular factor/inverse back into the full pseudoinverse, in the same fashion as the original C was padded. Unless that inner island is itself singular, too. $\endgroup$ Nov 7, 2017 at 19:32
  • $\begingroup$ There are also schemes to reorder/pivot the cholesky decomposition to make it suitable/rank-revealing for semidefinite systems. $\endgroup$ Nov 7, 2017 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.