GOPHERSPACE.DE - P H O X Y
gophering on hngopher.com
HN Gopher Feed (2017-09-27) - page 1 of 10
 
___________________________________________________________________
LU Factorization and Linear Systems for Programers
78 points by disaster01
http://dragan.rocks/articles/17/Clojure-Numerics-2-General-Linea...
___________________________________________________________________
 
[deleted]
 
[deleted]
 
dagss - 4 hours 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.
 
  FabHK - 1 hours 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.
 
    R_haterade - 2 minutes ago
    Excellent John D. Cook piece on the same subject:
    https://www.johndcook.com/blog/2010/01/19/dont-invert-that-m...
 
  gmiller123456 - 4 hours ago
  To save people a few clicks:[1]
  https://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/...
 
  FabHK - 1 hours ago
  > how to achieve high performanceYou 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
 
dragandj - 4 hours ago
The software used in the tutorial: https://github.com/uncomplicate/
neanderthalhttp://neanderthal.uncomplicate.org
 
hprotagonist - 5 hours ago
Golub's "Matrix Computations" remains a must-read reference text
here:
http://web.mit.edu/ehliu/Public/sclark/Golub%20G.H.,%20Van%2...
 
  krona - 4 hours ago
  Thanks. Theres a 4th edition from 2013, which I've just bought.
 
  stephencanon - 2 hours 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 - 1 hours 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)
 
leeoniya - 4 hours ago
interactive demo of LU-like and QR-like decomposition of affine
matrices: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
 
kxyvr - 1 hours 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 - 1 hours 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.)
 
    dragandj - 1 hours ago
    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.
 
      FabHK - 1 hours ago
      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.
 
        one-more-minute - 55 minutes ago
        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.
 
          FabHK - 35 minutes ago
          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.
 
          dragandj - 44 minutes ago
          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?
 
    kxyvr - 6 minutes ago
    Yes, I believe that dense methods on single core architectures
    have been well ported.  Of course, even there, it wasn't until
    I believe Eigen that we were able to moderately easily run
    automatic differentiation over the algebra and factorizations
    to get their derivatives.  Though, technically, yes, we could
    have run ADIFOR over LAPACK and hated ourselves in the
    process.It's been some time since we've seen significant speed
    increases for single core computers, so as new chips moved to
    multicore or GPU processing, the field followed.  That said,
    even moderately basic operations such as matrix multiplication
    weren't fully implemented until maybe three years.  I'm sure
    there were other groups, but that's when NVIDIA released
    cuBLAS-XT, which allowed multiplication of matrices that didn't
    fit on the GPU.So, yes, I agree that this is nontrivial and
    difficult.  At the academic level, numerical linear algebra
    tends to be in the domain of applied math departments, so I
    think that sometimes well suited people in CS don't get
    introduced to the field.  That's partly why I enjoy seeing
    articles on HN that discuss it.
 
  FabHK - 30 minutes ago
  > so the impact is wideYeah - in Shrek (2001), there's a scene
  where the eponymous ogre takes a mud bath in the swamp. When I
  studied this linear algebra stuff, one of my (more advanced)
  fellow students told me that their research group had "done" the
  mud... :-)