I'm a bit rusty on my linear algebra. Without inverting the matrix, what's the algorithm for solving Ax=b? Factorization is mentioned but I'm not sure what it means in this case.
Gaussian elimination. Remember where you subtract multiples of one row from another, in order to get zeros below the diagonal? That way you can back-solve a system of equations by reading off solutions starting with the bottom-right corner.
The result of carrying out the elimination steps is a triangular matrix with all zeros below the main diagonal. Call it U, for "upper".
The elimination steps themselves can be encoded as a matrix, rather than describing them in words. For instance, "subtract 2 times row 1 from row 2" is the matrix:
[ 1 0 0]
[-2 1 0]
[ 0 0 1]
So, to solve the system, you can perform elimination on A x = b, which gives you two matrices (U, and the elimination steps L). In some sense, you've factored A in to L and U.
Example:
Using Strang's terminology, let's call the matrix that "eliminates (i.e. sets to 0) the (i,j) entry" Eij
So, to solve a 3x3 system A x = b, you want to do elimination as: (E32 E31 E21) A = U
I.e.
1) Take A and subtract some multiple of row 2 from row 1, to put 0 in (2,1)
2) Take the result and subtract some multiple of row 3 from row 1 to put 0 in (3,1)
3) Take the result and subtract some multiple of row 3 from row 2 to put 0 in (3,2)
Now the goal of elimination is to get an uppper-triangular matrix U. I.e. if you were solving some system of equations in x, y, z, then because U is triangular, you could just read off the answer to z. Then plug that into the equation for row 2 to solve for y, then plug that into the equation for row 1 to solve for x. That's back-substitution.
BUT, let's return to our operations. We have:
E32 E31 E21 A = U
(E32 E31 E21) A = U
inverse(E32 E31 E21) (E32 E31 E21) A = inverse(E32 E31 E21) U
I A = inverse(E32 E31 E21) U
Now, it just so happens that the inverse of (E32 E31 E21) is lower-triangular. I.e. all entries above the diagonal are zero, so let's call it L.
I A = L U
A = L U
LU-decomposition. Or "factorization", since we've factored A into two matrices: L, and U.
Yes, but it's a little more complicated than that. In most cases you cannot trust the solution of straight Gaussian elimination. You need to take additional measures to avoid excessive error growth, such as using partial pivoting.
LU decomposition is the most common choice, which is a similar operation to matrix inversion in that you're iterating row (or column) operations on both sides of an equation. Numerical Recipes has a great section on this, IIRC.
But that kind of obscures the fact that, while algorithms like LU decomposition are faster than a pure matrix inversion, they're faster by a constant factor of about 2 or so. It's not a fundamentally different algorithm, just a different set of tricks to solve a bunch of simultaneous equations. If you aren't hugely performance constrained and already have a matrix inversion routine, you'd never be bothered to re-write your equation solver just for this speedup.
No. They're both O(N^3) in the dimension of the matrix. And it's more stable only in the sense that it's doing fewer operations. There's actually a great trick (again, in Numerical Recipes) when doing inversion to spend a few more (constant factor) cycles to iterate a correction to the limit of the precision of your number representation.
Really, check Numerical Recipes. I like their section a lot, though I don't have my copy handy to give you a cite.
The asymptotics are completely different for sparse and banded matrices. And NR really isn't a good reference unless you are looking for a mildly flawed 30 year old analysis.
Numerical Recipes has both sound theory exposition and working code (yes, with some bugs here and there). This blog post has neither. Nor, to my knowledge, does any other single reference out there. It's easily the best "how to think about numerics work" text I know; if you have a better one please point me there.
It's easy for a wonk to flame about something as pedestrian as a how-to text for working scientists, but it's really not helpful.
In addition to those suggestions, here are a few more relevant to the present discussion. For dense matrices and how to think about linear algebra algorithms, I highly recommend Trefethen and Bau "Numerical Linear Algebra". For sparse direct solvers, Tim Davis' book (http://www.ec-securehost.com/SIAM/FA02.html) is good, and you can transition from the "learning"-level implementation to Umfpack which is his production-quality solver (it's what is behind Matlab's backslash).
For iterative solvers, Saad is pretty standard. For multiphysics solvers, I don't think any book can match the Knoll and Keyes' 2004 review on Jacobian-free Newton-Krylov methods (it's very accessible).
The GSL has better implementations for many general-purpose algorithms in NR. SciPy is great if you work in Python. For scalable linear and nonlinear solvers, look at PETSc.
I guess I'll have to go back to my lecture notes. There was a reason why we did LU-decomposition (and similar things) instead of inverting matrices directly.
As others have said, Gaussian elimination, or if the matrix is not invertible (square and singular or not square) then the Moore-Penrose pseudoinverse can be used to find the least squares solution.
Except that the normal equations are rarely a good way to compute the action of the pseudoinverse (though explicitly computing it would be even worse). QR is the workhorse for under- and over-determined dense systems, with SVD as a fallback for especially nasty ones.
Comments
I'm a bit rusty on my linear algebra. Without inverting the matrix, what's the algorithm for solving Ax=b? Factorization is mentioned but I'm not sure what it means in this case.
Gaussian elimination. Remember where you subtract multiples of one row from another, in order to get zeros below the diagonal? That way you can back-solve a system of equations by reading off solutions starting with the bottom-right corner.
The result of carrying out the elimination steps is a triangular matrix with all zeros below the main diagonal. Call it U, for "upper".
The elimination steps themselves can be encoded as a matrix, rather than describing them in words. For instance, "subtract 2 times row 1 from row 2" is the matrix:
[ 1 0 0]
[-2 1 0]
[ 0 0 1]
So, to solve the system, you can perform elimination on A x = b, which gives you two matrices (U, and the elimination steps L). In some sense, you've factored A in to L and U.
Example:
Using Strang's terminology, let's call the matrix that "eliminates (i.e. sets to 0) the (i,j) entry" Eij
So, to solve a 3x3 system A x = b, you want to do elimination as: (E32 E31 E21) A = U
I.e.
1) Take A and subtract some multiple of row 2 from row 1, to put 0 in (2,1)
2) Take the result and subtract some multiple of row 3 from row 1 to put 0 in (3,1)
3) Take the result and subtract some multiple of row 3 from row 2 to put 0 in (3,2)
Now the goal of elimination is to get an uppper-triangular matrix U. I.e. if you were solving some system of equations in x, y, z, then because U is triangular, you could just read off the answer to z. Then plug that into the equation for row 2 to solve for y, then plug that into the equation for row 1 to solve for x. That's back-substitution.
BUT, let's return to our operations. We have:
E32 E31 E21 A = U
(E32 E31 E21) A = U
inverse(E32 E31 E21) (E32 E31 E21) A = inverse(E32 E31 E21) U
I A = inverse(E32 E31 E21) U
Now, it just so happens that the inverse of (E32 E31 E21) is lower-triangular. I.e. all entries above the diagonal are zero, so let's call it L.
I A = L U
A = L U
LU-decomposition. Or "factorization", since we've factored A into two matrices: L, and U.
Yes, but it's a little more complicated than that. In most cases you cannot trust the solution of straight Gaussian elimination. You need to take additional measures to avoid excessive error growth, such as using partial pivoting.
see http://en.wikipedia.org/wiki/Pivoting and http://en.wikipedia.org/wiki/Gaussian_elimination
LU decomposition is the most common choice, which is a similar operation to matrix inversion in that you're iterating row (or column) operations on both sides of an equation. Numerical Recipes has a great section on this, IIRC.
But that kind of obscures the fact that, while algorithms like LU decomposition are faster than a pure matrix inversion, they're faster by a constant factor of about 2 or so. It's not a fundamentally different algorithm, just a different set of tricks to solve a bunch of simultaneous equations. If you aren't hugely performance constrained and already have a matrix inversion routine, you'd never be bothered to re-write your equation solver just for this speedup.
As far as I know LU decomposition is more than a constant factor faster than inversion. And it is numerically more stable.
No. They're both O(N^3) in the dimension of the matrix. And it's more stable only in the sense that it's doing fewer operations. There's actually a great trick (again, in Numerical Recipes) when doing inversion to spend a few more (constant factor) cycles to iterate a correction to the limit of the precision of your number representation.
Really, check Numerical Recipes. I like their section a lot, though I don't have my copy handy to give you a cite.
The asymptotics are completely different for sparse and banded matrices. And NR really isn't a good reference unless you are looking for a mildly flawed 30 year old analysis.
Numerical Recipes has both sound theory exposition and working code (yes, with some bugs here and there). This blog post has neither. Nor, to my knowledge, does any other single reference out there. It's easily the best "how to think about numerics work" text I know; if you have a better one please point me there.
It's easy for a wonk to flame about something as pedestrian as a how-to text for working scientists, but it's really not helpful.
You are probably familiar with this page
http://www.fceia.unr.edu.ar/~fisicomp/apuntes/biblios/wnotnr...
and the suggested alternatives
http://www.fceia.unr.edu.ar/~fisicomp/apuntes/biblios/altnr....
In addition to those suggestions, here are a few more relevant to the present discussion. For dense matrices and how to think about linear algebra algorithms, I highly recommend Trefethen and Bau "Numerical Linear Algebra". For sparse direct solvers, Tim Davis' book (http://www.ec-securehost.com/SIAM/FA02.html) is good, and you can transition from the "learning"-level implementation to Umfpack which is his production-quality solver (it's what is behind Matlab's backslash). For iterative solvers, Saad is pretty standard. For multiphysics solvers, I don't think any book can match the Knoll and Keyes' 2004 review on Jacobian-free Newton-Krylov methods (it's very accessible).
The GSL has better implementations for many general-purpose algorithms in NR. SciPy is great if you work in Python. For scalable linear and nonlinear solvers, look at PETSc.
> And NR really isn't a good reference unless you are looking for a mildly flawed 30 year old analysis.
Could you recommend a better book?
They have old editions of Numerical Recipes online: http://www.nrbook.com/a/bookcpdf.php
Sadly they're using the horrid 'File open' encryption on the PDFs.
I guess I'll have to go back to my lecture notes. There was a reason why we did LU-decomposition (and similar things) instead of inverting matrices directly.
As others have said, Gaussian elimination, or if the matrix is not invertible (square and singular or not square) then the Moore-Penrose pseudoinverse can be used to find the least squares solution.
Except that the normal equations are rarely a good way to compute the action of the pseudoinverse (though explicitly computing it would be even worse). QR is the workhorse for under- and over-determined dense systems, with SVD as a fallback for especially nasty ones.