Numerically stable way to compute conditional covariance matrix using linalg.solve
Please don't inv—it's not as bad as most people think, but there's easier ways: you mentioned how np.linalg.solve(A, b) equals A^{-1} . b, but there's no requirement on what b is. You can use solve to solve your question, A - np.dot(B, np.linalg.solve(D, C)).
(Note, if you're doing blockwise matrix inversion, C is likely B.transpose(), right?)
I know that the recommendation is not to use linalg.inv and use linalg.solve when inverting matrices. This makes sense when I have situation like Ax = b and I want to get x, but is there a way to compute something like: A - B * D^{-1} * C without using linalg.inv? Or what is the most numerically stable way to deal with the inverse in the expression?
Thanks!