Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

LU Factorization and Linear Systems for Programers

dragan.rocks

Post reply on HN
Print view View original

41–46 of 46 posts

Re: LU Factorization and Linear Systems for Programers

#41

Post by JustFinishedBSG »

radarsat1 wrote: post #29
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…

> 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. Honestly, it might be tricky, but implementing matrix operations is not rocket science either. I find it incredible that so many projects rely on NVidia's proprietary libraries for doing this on GPU. Maybe there is some secret juice th…

> implementing matrix operations is not rocket science either.

Naively implementing them sure, but looking at OpenBLAS/LAPACK/etc code it seems pretty close to rocket science

Re: LU Factorization and Linear Systems for Programers

#42

Post by dragandj »

newen wrote: post #28
dragandj wrote: post #19

Earlier quoted context omitted.

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?

Well if you define a new type, and * and + operations for it, you will get matrix multiplication for free in native julia. I think you can do similar things for a few other algorithms in native julia. I don't think it will be MKL speed though because MKL uses cache size information etc. and is very optimized. Of course julia uses BLAS/Lapack etc for the types it can. GenericSVD.jl is an interesting example of how you…

When I hear "you just do X", I become sceptical of that instantly :)

Re: LU Factorization and Linear Systems for Programers

#43

Post by kxyvr »

dragandj wrote: post #42
newen wrote: post #28

Earlier quoted context omitted.

Well if you define a new type, and * and + operations for it, you will get matrix multiplication for free in native julia. I think you can do similar things for a few other algorithms in native julia. I don't think it will be MKL speed though because MKL uses cache size information etc. and is very optimized. Of course julia uses BLAS/Lapack etc for the types it can. GenericSVD.jl is an interesting example of how you…

When I hear "you just do X", I become sceptical of that instantly :)

As a note, this feature isn't exclusive to Julia. For example, Eigen implements something similar and the requirements for doing so are specified here:

https://eigen.tuxfamily.org/dox-devel/TopicCustomizing_Custo...

The reason to do so is to use a more exotic type that gives additional information. In the simplest case, we may want to use quad precision, which BLAS and LAPACK do not support by default (but it is possible to hack in). Alternatively, it's a good way to run automatic differentiation through the solve, or use infinite precision libraries like GMP, or use interval arithmetic.

So, yes, it's right to be skeptical, but it is possible, some libraries support it, and it's a useful trick to have available.

Re: LU Factorization and Linear Systems for Programers

#44

Post by dragandj »

kxyvr wrote: post #43
dragandj wrote: post #42

Earlier quoted context omitted.

When I hear "you just do X", I become sceptical of that instantly :)

As a note, this feature isn't exclusive to Julia. For example, Eigen implements something similar and the requirements for doing so are specified here: https://eigen.tuxfamily.org/dox-devel/TopicCustomizing_Custo... The reason to do so is to use a more exotic type that gives additional information. In the simplest case, we may want to use quad precision, which BLAS and LAPACK do not support by default (but it is poss…

It's a useful trick, but I was talking about comparison with what MKL and cuBLAS and other Nvidia libraries do. For the float/double/half implementations of "standard" operations that 95% use cases fall into, of which there are hundreds if not thousands of, I haven't seen even a close match.

Of course, if you need something special, there are various techniques. My preference for those things is to skip the middleman and code them in CUDA/OpenCL kernels directly...

OTOH, I'm interested in machine learning applications rather than physics/engineering. Double is an overkill here, float is the sweet spot, and some techniques use half or even less precision. Quad-precision may just not be the case that I need, but for people who do, I suppose they'd have their ways of solving that.

Re: LU Factorization and Linear Systems for Programers

#45

Post by stephencanon »

FabHK wrote: post #11
stephencanon wrote: post #9

Earlier quoted context omitted.

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…

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.

Re: LU Factorization and Linear Systems for Programers

#46

Post by dagss »

FabHK wrote: post #10
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.

> 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

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.

Post reply on HN