LU Factorization and Linear Systems for Programers • Hacker News Board

Re: LU Factorization and Linear Systems for Programers

#11

Post by FabHK » Wed, Sep 27, 2017, 9:33 PM UTC

stephencanon wrote: post #9

hprotagonist wrote: post #3

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

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

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)

Re: LU Factorization and Linear Systems for Programers

#12

Post by kxyvr » Wed, Sep 27, 2017, 9:34 PM UTC

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.

Re: LU Factorization and Linear Systems for Programers

#13

Post by FabHK » Wed, Sep 27, 2017, 9:44 PM UTC

dagss wrote: post #5

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.

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.

Re: LU Factorization and Linear Systems for Programers

#14

Post by SFjulie1 » Wed, Sep 27, 2017, 9:46 PM UTC

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.

Re: LU Factorization and Linear Systems for Programers

#15

Post by FabHK » Wed, Sep 27, 2017, 9:50 PM UTC

kxyvr wrote: post #12

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 strategie…

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

Re: LU Factorization and Linear Systems for Programers

#16

Post by dragandj » Wed, Sep 27, 2017, 9:52 PM UTC

FabHK wrote: post #15

kxyvr wrote: post #12

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 strategie…

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, an…

I think you are wrong about Julia, at least when it comes to linear algebra. Julia also relies on BLAS/LAPACK, at least those libraries that strive to be fast.

Re: LU Factorization and Linear Systems for Programers

#17

Post by FabHK » Wed, Sep 27, 2017, 9:59 PM UTC

dragandj wrote: post #16

FabHK wrote: post #15

Earlier quoted context omitted.

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, an…

I think you are wrong about Julia, at least when it comes to linear algebra. Julia also relies on BLAS/LAPACK, at least those libraries that strive to be fast.

you're absolutely right, I shouldn't have said "most". It basically calls BLAS/LAPACK (as every other language, and as one should), and does almost everything beyond that in Julia.

Re: LU Factorization and Linear Systems for Programers

#18

Post by one-more-minute » Wed, Sep 27, 2017, 10:04 PM UTC

FabHK wrote: post #17

dragandj wrote: post #16

Earlier quoted context omitted.

I think you are wrong about Julia, at least when it comes to linear algebra. Julia also relies on BLAS/LAPACK, at least those libraries that strive to be fast.

you're absolutely right, I shouldn't have said "most". It basically calls BLAS/LAPACK (as every other language, and as one should), and does almost everything beyond that in Julia.

To be fair, almost all of those linalg kernels have pure-Julia implementations (mainly so they work across any numeric type). I think "most" is a reasonable description, with a few special cases that reach out to BLAS and co.

Re: LU Factorization and Linear Systems for Programers

#19

Post by dragandj » Wed, Sep 27, 2017, 10:15 PM UTC

one-more-minute wrote: post #18

FabHK wrote: post #17

Earlier quoted context omitted.

you're absolutely right, I shouldn't have said "most". It basically calls BLAS/LAPACK (as every other language, and as one should), and does almost everything beyond that in Julia.

To be fair, almost all of those linalg kernels have pure-Julia implementations (mainly so they work across any numeric type). I think "most" is a reasonable description, with a few special cases that reach out to BLAS and co.

Can you give an example in popular Julia library(ies) that do that, and, very important, to speed comparisons? Are those operations close to the speed found in MKL for example?

Re: LU Factorization and Linear Systems for Programers

#20

Post by FabHK » Wed, Sep 27, 2017, 10:24 PM UTC

one-more-minute wrote: post #18

FabHK wrote: post #17

Earlier quoted context omitted.

you're absolutely right, I shouldn't have said "most". It basically calls BLAS/LAPACK (as every other language, and as one should), and does almost everything beyond that in Julia.

To be fair, almost all of those linalg kernels have pure-Julia implementations (mainly so they work across any numeric type). I think "most" is a reasonable description, with a few special cases that reach out to BLAS and co.

Good point - Julia does LinAlg with BigFloats and other "exotic" numbers, if you wish. But it definitely also integrates LAPACK/BLAS, and Tim Davis' SuiteSparse.