I'm interesting in solving a linear IV-GMM (see page 5 and 6 for background). The solution takes the form $$ \hat{\beta} = (X'ZWZ'X)^{-1} X'ZWZ'y $$
where $W$ is a positive definite weighting matrix and $Z'X$ has full rank. I'm interested in calculating this multiple times for different values of $y$, so I'm happy to compute an expensive factorisation of $X$, $Z$ and $W$ up-front and then reuse that.
As $W$ is positive definite then $X'ZWZ'X$ is also positive definite so this hints at the Cholesky decomposition but I'm aware that this is discouraged due to numerical stability issues.
One idea I had was to use the Cholesky decomposition $W=LL' = U'U$ and then use the QR decomposition of $UZ'X$, which would simplify to $$ R\hat{\beta} = Q'UZ'y. $$
Any hints or good references as to how to solve this would be appreciated.