
LU Factorization and Linear Systems for Programers - disaster01
http://dragan.rocks/articles/17/Clojure-Numerics-2-General-Linear-Systems-and-LU-Factorization
======
kxyvr
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
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.)

~~~
radarsat1
> 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
that I just don't know about, but it seems to me there can't be _that_ many
optimisation shortcuts for matrix multiplications and the like that require
intimate, secret knowledge of the hardware.

~~~
egocodedinsol
The hard part isn’t getting it basically right. The hard part is knowing where
subtle bugs are.

I remember my numerical computing textbook had one major takeaway: almost
always write your own matrix algorithms; almost never use them.

An extraordinary amount of work has gone into linear algebra computation and
much of it for mission critical stuff. It might not be “rocket science” but
it’s close.

~~~
kxyvr
This is an important point, but I don't think it's insurmountable for a
broader group of people.

Though technically not required, I think it helps to know a good deal of
mathematical theory as well in order to check properties that should hold,
which helps in tracking down subtle bugs. In my experience, tracking these
bugs is a royal pain in the ass and writing good tests for these properties is
tedious. Even then, there's a lot of funny stuff that goes on. One of my
favorites is that it's possible to take the Choleski factorization of an
indefinite matrix in some situations. That means it's not safe to use a
Choleski to check for whether or a not a matrix is positive definite for
mission critical functions:

[https://scicomp.stackexchange.com/a/26223](https://scicomp.stackexchange.com/a/26223)

But, anyway, mostly I contend that writing these algorithms is extremely
tedious, so there's not that many people with the patience. I also still
maintain that there's room for contribution here. Although most often these
routines are written in C and Fortran, the parallel tooling for these
languages is less sophisticated than that of modern languages. I've been
wondering whether a language like Rust would make it easier to code and
parallelize some of these algorithms.

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

~~~
FabHK
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.

~~~
R_haterade
Excellent John D. Cook piece on the same subject:
[https://www.johndcook.com/blog/2010/01/19/dont-invert-
that-m...](https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/)

~~~
conistonwater
That's a slightly different issue: the issue with solving normal equations
directly is that the condition number of X^\top X is the square of the
condition number of X, so it's much more significant than lu-vs-inv. It can
matter a great deal for ill-conditioned linear regression problems where some
features resemble each other.

------
dragandj
The software used in the tutorial:
[https://github.com/uncomplicate/neanderthal](https://github.com/uncomplicate/neanderthal)

[http://neanderthal.uncomplicate.org](http://neanderthal.uncomplicate.org)

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

~~~
stephencanon
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
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)

~~~
hcho3
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-
Davi...](https://www.amazon.com/Fundamentals-Matrix-Computations-David-
Watkins/dp/0470528338)

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

------
leeoniya
interactive demo of LU-like and QR-like decomposition of affine matrices:

[http://frederic-wang.fr/decomposition-of-2d-transform-
matric...](http://frederic-wang.fr/decomposition-of-2d-transform-
matrices.html)

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)

