Block diagonal matrix operations with Eigen (C++)

Developers
1

I’m working on generalizing the embedded Laplace for cases where the likelihood admits a non-diagonal Hessian, but rather a block-diagonal Hessian (this ends up including the diagonal and fully dense cases).

Does Eigen support operations with block-diagonal matrices? In particular, I’m interested in cholesky decomposition, multiplication, and solve. A quick search didn’t return anything promising.

I have a good sense of how to write these things out, but I figured I’d check if the work had already been done.

2

They have sparse matrix support which should help here, no?

3

There’s this block sparse file in unsupported Eigen-unsupported: BlockSparseMatrix.h Source File

1 Like
4

Depending on the size you may also be able to use a BandMatrix, though I can’t seem to find the docs related to it besides the below. You might get better performance by just leaving 0’s on the non-block parts of the supers for smaller 2x2 like block diagonals. There’s also not a lot of gurantees on storage order in sparse matrices so though you’ll pay more for non-zero spaces you’ll get better SIMD stuff etc.

https://eigen.tuxfamily.org/index.php?title=SpecialMatrix#skyline.2Fband_matrix

5

Cool, thanks for the pointers. I’ll dig into them.

I tried using dense operations and it is way slower than in the diagonal case, so I’ll definitely need to exploit sparsity. (That said, sparsity is pretty organized, since it’s block diagonal, and there’s no need for re-ordering).

2 Likes
6

related Eigen sparse matrix question: is there an Eigen implementation that gets the inverse of a sparse matrix? need this to compute a piece of the puzzle that is the BYM2 model.
INA’s qinv: qinv: Computes (parts of) the inverse of a SPD sparse matrix in inbo/INLA: Full Bayesian Analysis of Latent Gaussian Models using Integrated Nested Laplace Approximations

7

I think you can use one of the sparse solvers from here like this stackoverflow answer

8

Why do you need the actual inverse, and not just something like an LU-decomposition?

9

need to compute the geometric mean of a sparse matrix, following proposal here: https://arxiv.org/pdf/1601.01180.pdf - paragraph near bottom of p.6

if this were available as a function in the Stan math library, that would be awesome for the folks using the BYM2 model for spatial smoothing of areal data. afaik, the only way this function is available is through package R-INLA, which isn’t available on CRAN. R-INLA’s usefulness is surpassed only by its inscrutability.

2 Likes