The Delassus operator commonly shows up in the equations of multibody dynamics with contact.

In matrix form, we have matrix M (a mass matrix, SPD, invertible) and a matrix J, not square in general.

The Delassus operator is defined as W = J⋅M⁻¹⋅Jᵀ.

When both J and M are sparse matrices, is there an efficient way to compute W with Eigen?

Thank you in advance.

In practice, computing matrix inverses explicitly is very inefficient, so it's worth trying to avoid that.
Because the actual M⁻¹ isn't of interest (but its product with other matrices is), we don't need to generate it, only its action.
Let M be have dimensions (b, b),
and let J have dimensions (a, b),
so that W has dimensions (a, a).

First, let's look at M⁻¹Jᵀ = P which has dimensions (b, a).
For each column, jᵀ of Jᵀ and p of P, we have the system, Mp = jᵀ.
In other words, we only need to solve "a" number of sparse systems of linear equations to get P.
With a sparse symmetric positive definite mass matrix (coming from a mesh), a multigrid-preconditioned conjugate gradient solver should do really well (converging very close to the true p within a handful of iterations).

In this way, you compute the action of M⁻¹ on Jᵀ. Then just multiply JP = W to get the operator.

I don't know if there's a better way, but this method seems more efficient than the naïve approach.