21–30 of 48 posts
I also strongly recommend the course Convex Optimization. Might help some people avoid neural nets when they don't need 'em ;)
On page 2 of http://www.ti3.tuhh.de/paper/rump/Ru08a.pdf there is an innocent-looking matrix for which the numerically computed inverse is off from the true inverse by 47 orders of magnitude.
However, one of the reasons he's given is not correct: Druinsky and Toledo have shown (http://arxiv.org/abs/1201.6035) that -- despite the very widespread belief to the contrary -- solving a linear system be calculating the inverse can be as accurate (though not nearly as efficient) as solving it directly.
They even speculate as to how this myth was able to sustain it for so long. Namely, there is a certain error bound for the solution as calculated by employing the inverse which is quite simple to deduce and which predicts a pretty large error. However, this bound is too conservative. There is a better bound, which ensures that the error is of the same kind as for the direct method. However, this bound is much harder to derive.
I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…
A x = b
A => QR # QR decomposition
Q R x = b
R x = Q' b # Q is inverted by transposing
x = backsubstitute(R, Q' b)
I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…
I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…
Incidentally, while solving `Ax = b` in general is definitely equivalent to inverting `A` (simply by taking `b` to be basis vectors, so that the solutions `x` are the columns of the inverse matrix in that basis), it is certainly possible to solve it in particular cases without inverting `A`; for example, the solution to `Ax = b`, where `b` is the first column of `A` in a particular basis, is the first vector `x` of that basis—no inversion required (nor even computation of `b`)!
I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…
No.
>How else would you solve Ax = b for x?
Call dgetrs instead of dgetri.
You use the factorisation directly (not the inverse) to get x from b. This is done by calling dgetrf to get the factorisation, then calling dgetrs to get the solution. You don't need to call dgetri.
I'd do an example, but there are other replies to your comment that go through the mathematics.
Well it turns out there is something like division for matrices and it's called factorization (followed by back-substitution) and there are actually many ways to do it, not just one (each having different numerical properties and complexity).
So whenever I see A^-1 * b in an equation I view it as something like b / A and I know there's the choice of "division methods" to use for this.
Is there more context around what exactly he's getting at? In games we use matrix inverses all the time(and since they're usually a specialized affine transform there's cheap shortcuts to calculate them) for stiff like hit detection or animation.
In my research, specifically, that spatial mesh can have as many as 10^7 grid points. For linear systems of this size, any exact solution at all is impractical - we exclusively employ iterative approximations built from terms of the form (A^n) b since matrix-vector products are about the only things we can compute in a reasonable amount of time.
I get that inverting a matrix can be computationally intensive, but what's the alternative? "What if you have to solve Ax = b for a lot of different b‘s? Surely then it’s worthwhile to find A-1. No. The first time you solve Ax = b, you factor A and save that factorization. Then when you solve for the next b, the answer comes much faster." I don't get it, although I'm not a numerics expert so maybe I'm missing somethi…
These solvers construct a solution by adding up terms of the form c_n * A^n * y and truncating the series when the residuals are small enough. This is especially efficient in PDEs, where problem matrices are sparse so matrix-vector products are efficient (i.e. O(n)), and in practice only a few terms of the series are necessary.
I am not a mathematician, but simply based on the fact that beautiful/simpler answers to questions seem to get closer to the ultimate truth, the current state of large-dimension matrix inversion is far from the ideal possible solution, because it is ugly as hell.
Arguably LU / LDLT factorizations before solving the direct system are the more elegant / beautiful / simpler answers if you take the algorithm at face value. However, large-dimension matrix inversion has problems outside of just beauty. Part of the problem is that we don't have fast, perfect representations of real or complex numbers. At best we have double precision floating points, which means that when you take t…
We'll probably never have perfect representation of real numbers, but Gustafson's unums seem to be the next best thing that can actually make a difference. His examples of how unums can solve the trickiest numerical challenges are impressive.
http://insidehpc.com/2015/03/slidecast-john-gustafson-explai...