Hacker News Reader

Comment by mwfunk | Hacker News Reader

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 something. He says, "the first time you solve Ax = b, you factor A and save that factorization."

Great. But don't you then use the factorization to determine the inverse of A (LAPACK example: call dgetrf to get the factorization, then call dgetri to compute the inverse using the factorization)? Maybe what he's talking about is going over my head. How else would you solve Ax = b for x, if not by computing the inverse of A (assuming A is a 3x3 or 4x4)?

View on HN

"How else would you solve Ax = b for x, if not by computing the inverse of A (assuming A is a 3x3 or 4x4)?"

For instance, there a fast iterative solvers. They start with an initial guess for the solution (which can be an arbitrary vector) and successively generate more accurate approximations.

These iterative solvers don't need access to the inverse of A. Furthermore, many of those solvers don't even need access to A: They only need access to a subroutine which, given a vector u, returns Au.

Many matrices arising in applications (for instance one often reduces the problem of solving partial differential equations to the problem of solving a (very huge) linear system of equations) are sparse (have only few non-zero entries). For these it can make sense to write specialized code which specifically calculates Au.

If you want to read more about this, then a good starting point is https://en.wikipedia.org/wiki/Kaczmarz_method.

You calculate the factorization (QR, LU, or similar), but you never calculate the inverse of the factorization:

    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)

Q was the transpose of it's inverse, and so you don't get any rounding error when you transpose it. R didn't need to be inverted because it was in a form where back substitution lets you apply the inverse of R to be without needing the inverse of R.

MATLAB's linsolve uses your strategy! (http://in.mathworks.com/help/matlab/ref/linsolve.html)

Linear models in R are also solved by QR factorization.

For example, you factorize A into LU, so that you have LUx = b.

Now you can let Ux = y, solve the system Ly = b for y, then solve Ux = y for x.

Solving each triangular system is fast: O(n^2).

If you need to do this with a different b, you can do it again, and you already have A factorized, which was the expensive part, i.e. O(n^3).

(Depending on the type of problem, you might want to factorize the original matrix in a different way, e.g. QR or Cholesky or whatever.)

IIRC, the HP-15C used LU decomposition in its matrix math routines.

Sure, assuming that A is full-rank and small, then its perfectly reasonable to invert A with adjoint(A)/discriminant(A). But this is very expensive, and sensitive to roundoff for larger A.

Lets say you are using QR factorization. Then R\b is easy to compute by back-substitution, because R is triangular. Q\b is also easy, because Q is orthogonal, and Q's inverse is its transpose, so you just calculate transpose(Q)*b.

Same thing for LDL^T factorization, LU factorization, etc. Its easy because back-substitution with a triangular matrix is easy. In the case of LDL^T factorization, D is diagonal, and again D\b is easy (its a bunch of rank-1 division operations).

>Great. But don't you then use the factorization to determine the inverse of A (LAPACK example: call dgetrf to get the factorization, then call dgetri to compute the inverse using the factorization)?

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.

In PDEs, we typically use something called a Krylov subspace iterative solver.

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 think that the idea is not that you don't take steps that, taken together, amount to inverting `A`, but rather that you don't just try to invert `A`. Thus, I think that the idea is to do something like "Compute `A = LU`, and then invert `L` and `U`"—which of course can be used to invert `A`, but which isn't exactly the same as "Compute `A^{-1}`" (using Cramer's rule, or some other absurdity).

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`)!

In LAPACK you would use dgetrs to solve for x directly with the LU factorization.

http://www.amazon.com/Applied-Numerical-Linear-Algebra-Demme...

Numerical linear algebra has a long and vast history. You really do need specialized information that your undergrad linear algebra classes would not have imparted.

For instance, the condition number is important: https://en.wikipedia.org/wiki/Condition_number

You would solve e.g.

min_x* ( ||A x* - b|| (+\lamba regularization(x*) )

Use a Pseudoinverse or a decomposition instead. A great book to read is Matrix Computations by Golub/van Loan.

Reply on news.ycombinator.com