
Clojure Numerics: Orthogonalization and Least Squares - dragandj
http://dragan.rocks/articles/17/Clojure-Numerics-5-Orthogonalization-and-Least-Squares
======
dragandj
Source code of the library is at:
[https://github.com/uncomplicate/neanderthal](https://github.com/uncomplicate/neanderthal)

------
_sdegutis
What’s the benefit of writing this in Clojure, instead of in Java and
designing it such that calling it from Clojure feels natural? Does it use any
Clojure specific features or have strict Clojure dependencies?

~~~
dragandj
In Java, it would require at least an order of magnitude more code.
Additionally, you'd have to know exactly how to do it from the beginning,
since it is much more difficult to experiment. High performance computing is
not something that is fabulously well documented, so you _have_ to poke around
a lot to discover how things work. To see what I'm talking about, you should
compare the codebase of neanderhtal (20.000 lines of code) with roughly
similar libraries in Java (ND4J is something like 400.000 lines of code), or
even Python (numpy is roughly 250.000).

~~~
tejinderss
Bit offtopic. Would you recommend clojure as a web development programming
language? I am from a python background and want to learn a new language.
Previously thought of Rust and Haskell, rust seems to be too low level and
carries mental burden of manual memory management and haskell seems too
complicated for me. Thanks in advance

~~~
dragandj
I would. Web dev is (for better or worse :) one of the main applications of
Clojure. See further at
[https://clojurescript.org/](https://clojurescript.org/) and
[https://pragprog.com/book/dswdcloj2/web-development-with-
clo...](https://pragprog.com/book/dswdcloj2/web-development-with-clojure-
second-edition)

------
Fede_V
I'm not very familiar with Clojure, but are the numerics done with unboxed
numbers? In the end, what every language that wants to do scientific
programming right is expose some sort of c-array like data structure that you
can pass onto c/fortran code.

~~~
dragandj
Do not worry. Neanderthal works with barest primitives, and uses Intel MKL,
cuBLAS, CLBlast, and custom kernels at the low level. This is practicaly as
fast as you can get in a general library.

~~~
tmyklebu
How do you use MKL and friends on objects in the Java heap? It seems like
you'd have to copy to the native heap, do your linear algebra operation, and
then copy the result back from the native heap.

~~~
dragandj
It is not used on objects on the java heap. No copying occurs.

~~~
tmyklebu
Aha, makes sense. (Now I see that at
[http://neanderthal.uncomplicate.org/articles/tutorial_native...](http://neanderthal.uncomplicate.org/articles/tutorial_native.html.))
Doesn't this make it very painful to do unusual operations on matrices? For
example, solving linear systems Ax=b where A is of the form

    
    
      [ 1                                             ]
      [ a_2 b_1  1                                    ]
      [ a_3 b_1  a_3 b_2  1                           ]
      [ ...      ...      ...      ...                ]
      [ a_n b_1  a_n b_2  a_n b_3  ...  a_n b_{n-1} 1 ].

~~~
dragandj
I do not understand what the form of that matrix has to do with java or native
heap. This seems to me as a completely orthogonal issue. As for the triangular
form, this is supported in Neanderthal, as well as the rest of special
structural sparse shapes.

~~~
tmyklebu
> I do not understand what the form of that matrix has to do with java or
> native heap. This seems to me as a completely orthogonal issue. As for the
> triangular form, this is supported in Neanderthal, as well as the rest of
> special structural sparse shapes.

You support triangular matrices. However, you can solve a linear system of the
form I gave in linear time, while it takes quadratic space and time to form
the corresponding triangular matrix and do a triangular solve.

Not all special matrices are supported by BLAS/LAPACK. Other common examples
might be block Toeplitz/Hankel matrices for which fast multiplication and fast
solvers are available. In order to support special (not in BLAS/LAPACK) matrix
operations naturally, you'd want natural, no-extra-copying access to the
vectors within Java or Clojure so that you can write the good algorithm
manually, as you'd do in C or Fortran.

~~~
dragandj
Sure! If you need to efficiently implement that yourself, you can use OpenCL
to implement the kernels on the CPU, or OpenCL & CUDA on the GPU (If that
makes sense). Check out ClojureCUDA and ClojureCL.

------
geokon
Thanks for writing this series. I'm really interested in this stuff and I'm
sorta writing my own version where I'm implementing the algorithms in lisp:
[http://geokon-gh.github.io/linearsystems/](http://geokon-
gh.github.io/linearsystems/) (I'm sure it's quite amateurish compared to what
you're doing)

Nothing quite forces you to understand the issues deeply like implementing it
in code!

I haven't gotten to LS/QR yet (just getting the LU working correctly was quite
tricky), but when I do I'll make sure to reference your series :)

