Linear systems of equations are pervasive in computational statistics. One special system I have encountered (e.g., in factor analysis) is the system
where $$A=D+ B \Omega B^T$$ Here $D$ is a $n\times n$ diagonal matrix with a strictly positive diagonal, $\Omega$ is an $m\times m$ (with $m\ll n$) symmetric positive semi-definite matrix, and $B$ is an arbitrary $n\times m$ matrix. We are asked to solve a diagonal linear system (easy) that has been perturbed by a low-rank matrix. The naive way to solve the problem above is to invert $A$ using Woodbury's formula. However, that doesn't feel right, since Cholesky and QR factorizations can usually speed up solution of linear systems (and normal equations) dramatically. I recently came up on the following paper, that seems to take the Cholesky approach, and mentions numerical instability of Woodbury's inversion. However, the paper seems in draft form, and I could not find numerical experiments or supporting research. What is the state of the art to solve the problem I described?