15 ms·

LU Factorization and Linear Systems for Programers

131 points by disaster01 9y ago

deleted 9y ago
[deleted]
hprotagonist 9y ago
Golub's "Matrix Computations" remains a must-read reference text here: http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%20Loan%20C.F.-%20Matrix%20Computations.pdf http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%2...
- krona 9y ago
  Thanks. Theres a 4th edition from 2013, which I've just bought.
- stephencanon 9y ago
  Do note that Golub & Van Loan is very much a reference text, however; it is not a great choice if you're just learning the subject (it covers everything, but without much depth and without much exposition).
  - FabHK 9y ago
    True. J. W. Demmel, Applied Numerical Linear Algebra, is a gentler introduction. However, if you're masochist enough to actually implement any of the algorithms, Golub & van Loan is a great reference. (Though, you really shouldn't implement it yourself except for didactic purposes - just use LAPACK/BLAS, which has been debugged for decades, and deals with all the special cases you're ignoring (underflow/overflow/nans/zeros/infs/...)) EDIT: Oh, and there's the excellent Numerical Linear Algebra by Trefethen and Bau, as kxyvr mentions. EDIT EDIT: Funny, on amazon the top reviews for both books mentioned above are identical. Seems like I'm not the only one having trouble keeping them apart... :-) (Demmel is more of an introduction, FWIW)
    - jacobolus 9y ago
      If anyone wants an even more introductory book, this one (not yet published, but PDF drafts available) looks nice: https://web.stanford.edu/~boyd/vmls/ https://web.stanford.edu/~boyd/vmls/
    - hcho3 9y ago
      I'd also recommend Fundamentals of Matrix Computations by Watkins for a great introduction. It helped immensely when I had to teach myself matrix algorithms as part of my undergraduate research. https://www.amazon.com/Fundamentals-Matrix-Computations-David-Watkins/dp/0470528338 https://www.amazon.com/Fundamentals-Matrix-Computations-Davi...
    - stephencanon 9y ago
      Demmel is also, well, Applied. He goes into (some of the) the nitty-gritty implementation details, where Trefethen and Bau stay at a higher level of algorithmic abstraction.
dagss 9y ago
For programmers, the really interesting part of dense linear algebra is how to achieve high performance, as blocking techniques have to be used to amortize loads from memory to cache. Google for Goto's "Anatomy of high-performance matrix multiplication", one of my favorite programming texts. Also the papers underlying the development of the "Elemental" library for distributed dense linear algebra is worth a look.
- gmiller123456 9y ago
  To save people a few clicks: [1] https://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/gotoPaper.pdf https://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/...
- FabHK 9y ago
  > how to achieve high performance You really just want to use LAPACK/BLAS, no? (That's what the Neanderthal library mentioned in the article does, btw, and basically linear algebra libraries for other languages, too. If not, you probably shouldn't use it...) http://neanderthal.uncomplicate.org http://neanderthal.uncomplicate.org
  - dagss 9y ago
    Of course you should but the blocking and amortization techniques have much wider usecases beyond linear algebra. Matrix multiplicaton is a useful "simple" example to study; of course you don't write your own if what you want to do is covered by BLAS/LAPACK. For instance I have written high performance spherical harmonic transforms, and these techniques apply then.
- FabHK 9y ago
  BTW, I think one important take-away for programmers is: Don't invert matrices. For example, the solution for the least squares (=linear regression) problem is beta^hat = (X'X)^-1 X'y and you might think that the best way to compute it is to compute X'X, invert it, do the multiplication, yada yada yada. But it's really better to just ask your library (Matlab, LAPACK, ...) to solve the problem for you, and it'll probably do a QR decomposition. For small problems, it doesn't really matter, but for big problems you'll be both faster and more accurate.
leeoniya 9y ago
interactive demo of LU-like and QR-like decomposition of affine matrices: http://frederic-wang.fr/decomposition-of-2d-transform-matrices.html http://frederic-wang.fr/decomposition-of-2d-transform-matric... i mostly copy-pasted the js code from this page as a contribution to this lib back in the day: https://github.com/epistemex/transformation-matrix-js https://github.com/epistemex/transformation-matrix-js
dragandj 9y ago
The software used in the tutorial: https://github.com/uncomplicate/neanderthal https://github.com/uncomplicate/neanderthal http://neanderthal.uncomplicate.org http://neanderthal.uncomplicate.org
kxyvr 9y ago
By the way, if anyone is interested in good open source opportunities, computational linear algebra is nowhere near a solved problem and there's good opportunity for impactful contribution. The computational challenges of the algebra versus factorizations is one angle. Dense versus sparse is another. Shared memory parallelization vs distributed memory vs GPUs is another. Even on the GPU, there are different strategies depending on whether or not the entire matrix fits on a single GPU or if we have to use multiple GPUs. Incomplete or multilevel direct methods used as effective preconditioners for iterative methods are also important. Hell, even efficient direct techniques embedded in indirect solvers is important. Part of the way to get started would be to look at something like a general numerical linear algebra book like Numerical Linear Algebra from Trefethen and Bau. There are better computational algorithms than what they present, but they do a good job at introducing important factorizations and why we care about them. Then, have a look at Tim Davis' book Direct Methods for Sparse Linear Systems. The codes in that book are online. Then, try to reimplement these algorithms in other languages, parallelize them, or make them better. These are good algorithms, but there are better and Tim's more recent codes are actively used by both MATLAB and Octave. Then, look for missing routines in open source libraries. For example, I just did a quick look and MAGMA currently lists missing routines between them and LAPACK. Anyway, it's not a field for everyone, but it's one that good architecture and parallelization knowledge can have a positive impact. Nearly all engineering codes depend on good solvers, so the impact is wide.
- FabHK 9y ago
  Though, the classic dense methods have been implemented and optimised and ported to death with LAPACK/BLAS, haven't they. What you're talking about is parallelisation, moving to GPU, and modern (combinatorial) methods for sparse systems, and that's fairly cutting edge, and not trivial to implement/port. You'll need to have a pretty good understanding of the language and its paradigmatic use, and of linear algebra, and of modern computer architecture. However, I assume that you're right in that there might still be some low-hanging fruit. BTW, Julia is an awesome language with excellent LA support, and it's nice in that most algorithms are coded in Julia itself (unlike the two-language situation in Python, Clojure (AFAIK), etc.)
SFjulie1 9y ago
For once I am happy to had an education because this was clearer when exposed by my teacher 25 years ago. Yep, in order to engineer microchip we had to learn how to code SPICE and were never aware of Djikstra BS or math but had to deal with real life capacitors and resistors and how to model them.
flor1s 9y ago
If you want to learn about this stuff but need a more gentler introduction, check out Robert van de Geijn's MOOC at http://www.ulaff.net http://www.ulaff.net (if there is currently no session, just download the lecture notes and use those, they have lecture videos embedded).
compumike 9y ago
For a gentler introduction to LU factorization and how useful the decomposition is for efficiently re-solving the same system multiple times, I'd offer up the "Systems of Equations" section of the "Ultimate Electronics" book I've been working on: https://www.circuitlab.com/textbook/systems-of-equations/ https://www.circuitlab.com/textbook/systems-of-equations/ I tried to go back and forth between the 5x+2y=3 style equations that most people are familiar with and the matrix forms.

131 points by disaster01 9y ago

deleted 9y ago

[deleted]

Golub's "Matrix Computations" remains a must-read reference text here: http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%20Loan%20C.F.-%20Matrix%20Computations.pdf http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%2...

krona 9y ago
Thanks. Theres a 4th edition from 2013, which I've just bought.
stephencanon 9y ago
Do note that Golub & Van Loan is very much a reference text, however; it is not a great choice if you're just learning the subject (it covers everything, but without much depth and without much exposition).
- FabHK 9y ago
  True. J. W. Demmel, Applied Numerical Linear Algebra, is a gentler introduction. However, if you're masochist enough to actually implement any of the algorithms, Golub & van Loan is a great reference. (Though, you really shouldn't implement it yourself except for didactic purposes - just use LAPACK/BLAS, which has been debugged for decades, and deals with all the special cases you're ignoring (underflow/overflow/nans/zeros/infs/...)) EDIT: Oh, and there's the excellent Numerical Linear Algebra by Trefethen and Bau, as kxyvr mentions. EDIT EDIT: Funny, on amazon the top reviews for both books mentioned above are identical. Seems like I'm not the only one having trouble keeping them apart... :-) (Demmel is more of an introduction, FWIW)
  - jacobolus 9y ago
    If anyone wants an even more introductory book, this one (not yet published, but PDF drafts available) looks nice: https://web.stanford.edu/~boyd/vmls/ https://web.stanford.edu/~boyd/vmls/
  - hcho3 9y ago
    I'd also recommend Fundamentals of Matrix Computations by Watkins for a great introduction. It helped immensely when I had to teach myself matrix algorithms as part of my undergraduate research. https://www.amazon.com/Fundamentals-Matrix-Computations-David-Watkins/dp/0470528338 https://www.amazon.com/Fundamentals-Matrix-Computations-Davi...
  - stephencanon 9y ago
    Demmel is also, well, Applied. He goes into (some of the) the nitty-gritty implementation details, where Trefethen and Bau stay at a higher level of algorithmic abstraction.

dagss 9y ago

For programmers, the really interesting part of dense linear algebra is how to achieve high performance, as blocking techniques have to be used to amortize loads from memory to cache. Google for Goto's "Anatomy of high-performance matrix multiplication", one of my favorite programming texts. Also the papers underlying the development of the "Elemental" library for distributed dense linear algebra is worth a look.

gmiller123456 9y ago
To save people a few clicks: [1] https://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/gotoPaper.pdf https://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/...
FabHK 9y ago
> how to achieve high performance You really just want to use LAPACK/BLAS, no? (That's what the Neanderthal library mentioned in the article does, btw, and basically linear algebra libraries for other languages, too. If not, you probably shouldn't use it...) http://neanderthal.uncomplicate.org http://neanderthal.uncomplicate.org
- dagss 9y ago
  Of course you should but the blocking and amortization techniques have much wider usecases beyond linear algebra. Matrix multiplicaton is a useful "simple" example to study; of course you don't write your own if what you want to do is covered by BLAS/LAPACK. For instance I have written high performance spherical harmonic transforms, and these techniques apply then.
FabHK 9y ago
BTW, I think one important take-away for programmers is: Don't invert matrices. For example, the solution for the least squares (=linear regression) problem is beta^hat = (X'X)^-1 X'y and you might think that the best way to compute it is to compute X'X, invert it, do the multiplication, yada yada yada. But it's really better to just ask your library (Matlab, LAPACK, ...) to solve the problem for you, and it'll probably do a QR decomposition. For small problems, it doesn't really matter, but for big problems you'll be both faster and more accurate.

leeoniya 9y ago

interactive demo of LU-like and QR-like decomposition of affine matrices: http://frederic-wang.fr/decomposition-of-2d-transform-matrices.html http://frederic-wang.fr/decomposition-of-2d-transform-matric... i mostly copy-pasted the js code from this page as a contribution to this lib back in the day: https://github.com/epistemex/transformation-matrix-js https://github.com/epistemex/transformation-matrix-js

dragandj 9y ago

The software used in the tutorial: https://github.com/uncomplicate/neanderthal https://github.com/uncomplicate/neanderthal http://neanderthal.uncomplicate.org http://neanderthal.uncomplicate.org

kxyvr 9y ago

By the way, if anyone is interested in good open source opportunities, computational linear algebra is nowhere near a solved problem and there's good opportunity for impactful contribution. The computational challenges of the algebra versus factorizations is one angle. Dense versus sparse is another. Shared memory parallelization vs distributed memory vs GPUs is another. Even on the GPU, there are different strategies depending on whether or not the entire matrix fits on a single GPU or if we have to use multiple GPUs. Incomplete or multilevel direct methods used as effective preconditioners for iterative methods are also important. Hell, even efficient direct techniques embedded in indirect solvers is important. Part of the way to get started would be to look at something like a general numerical linear algebra book like Numerical Linear Algebra from Trefethen and Bau. There are better computational algorithms than what they present, but they do a good job at introducing important factorizations and why we care about them. Then, have a look at Tim Davis' book Direct Methods for Sparse Linear Systems. The codes in that book are online. Then, try to reimplement these algorithms in other languages, parallelize them, or make them better. These are good algorithms, but there are better and Tim's more recent codes are actively used by both MATLAB and Octave. Then, look for missing routines in open source libraries. For example, I just did a quick look and MAGMA currently lists missing routines between them and LAPACK. Anyway, it's not a field for everyone, but it's one that good architecture and parallelization knowledge can have a positive impact. Nearly all engineering codes depend on good solvers, so the impact is wide.

FabHK 9y ago
Though, the classic dense methods have been implemented and optimised and ported to death with LAPACK/BLAS, haven't they. What you're talking about is parallelisation, moving to GPU, and modern (combinatorial) methods for sparse systems, and that's fairly cutting edge, and not trivial to implement/port. You'll need to have a pretty good understanding of the language and its paradigmatic use, and of linear algebra, and of modern computer architecture. However, I assume that you're right in that there might still be some low-hanging fruit. BTW, Julia is an awesome language with excellent LA support, and it's nice in that most algorithms are coded in Julia itself (unlike the two-language situation in Python, Clojure (AFAIK), etc.)

SFjulie1 9y ago

For once I am happy to had an education because this was clearer when exposed by my teacher 25 years ago. Yep, in order to engineer microchip we had to learn how to code SPICE and were never aware of Djikstra BS or math but had to deal with real life capacitors and resistors and how to model them.

flor1s 9y ago

If you want to learn about this stuff but need a more gentler introduction, check out Robert van de Geijn's MOOC at http://www.ulaff.net http://www.ulaff.net (if there is currently no session, just download the lecture notes and use those, they have lecture videos embedded).

compumike 9y ago

For a gentler introduction to LU factorization and how useful the decomposition is for efficiently re-solving the same system multiple times, I'd offer up the "Systems of Equations" section of the "Ultimate Electronics" book I've been working on: https://www.circuitlab.com/textbook/systems-of-equations/ https://www.circuitlab.com/textbook/systems-of-equations/ I tried to go back and forth between the 5x+2y=3 style equations that most people are familiar with and the matrix forms.