0

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!

StatsNoob
  • 360
  • 1
  • 5
  • 15
  • Unless there is some special underlying structure, inverting is probably your best bet. What are those variables, especially, what are their dimensions? – Nico Schertler Dec 24 '19 at 22:10
  • The specific application I have in mind is getting the conditional variance of multivariate Gaussian. A,B,C,D corresponds to the block matrices of the original covariance matrix of unconditional Gaussian. So A and D would be symmetric positive definite for example. The dimension could vary depending on the uses each time I suppose. – StatsNoob Dec 24 '19 at 22:17

1 Answers1

0

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?)

Ahmed Fasih
  • 6,458
  • 7
  • 54
  • 95
  • Can you back up your claim (that using `solve` is better than `inv`) by some reference or experiments? I'm not saying that it's not - I'm just curious. – Nico Schertler Dec 25 '19 at 17:42
  • It's conventional wisdom isn't it? [1](https://stackoverflow.com/q/1419580/500207) [2](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/) [3](https://www.r-bloggers.com/dont-invert-that-matrix-why-and-how/) but see also [4](https://arxiv.org/abs/1201.6035) whose abstract opens with "Several widely-used textbooks lead the reader to believe that solving a linear system of equations Ax = b by multiplying the vector b by a computed inverse inv(A) is inaccurate. Virtually all other textbooks on numerical analysis and numerical linear algebra advise against using computed inverses". – Ahmed Fasih Dec 25 '19 at 18:06
  • These all apply to solving linear systems, where `b` is just a vector. In fact, that last reference even says that calculating the inverse is not less accurate than solving the linear system. Here in this case, we would, however, need to solve several linear systems to get the result. That's why it's not exactly the same thing. And I am not sure how much of the analysis translates. – Nico Schertler Dec 25 '19 at 18:12
  • Whether you're solving for vector `b` or for matrix `C` as in op's question, it doesn't matter: `solve` will calculate `D`'s Cholesky decomposition (since `D` comes from a covariance matrix, it will almost assuredly be hermitian) and then use that factorization to back-solve `b` or `C` or any right-hand-side, just like `inv`. Now. It's true that Druinsky & Toledo show that `inv` isn't as bad as people make it out, and it's unlikely that op's matrixes will have the pathological issues that cause `inv` to fail, but humans are herd animals and most people will unthinkingly shout about `inv`. – Ahmed Fasih Dec 25 '19 at 18:19
  • 1
    If you absolutely insist I give you a technical reason to use `solve` instead of `inv` even when your data is non-pathological, here is one: `invD * C` is a full dense matrix-multiplication, O(N^3) assuming `D` and `C` are both N by N. Meanwhile `scipy.linalg.solve_triangular(cholD, C)` back-substitutes a **triangular** matrix, which is *half* the work, O(0.5 N^3). (One could still argue that GEMM (GEneral Matrix-Matrix Multiplication) implementations are insanely hyper-optimized, and so the first might wind up being faster than the second on a real computer… ) – Ahmed Fasih Dec 25 '19 at 18:36