4 ms·

Clojure Numerics: Orthogonalization and Least Squares

118 points by dragandj 9y ago

dragandj 9y ago
Source code of the library is at: https://github.com/uncomplicate/neanderthal https://github.com/uncomplicate/neanderthal
Fede_V 9y ago
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 9y ago
  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 9y ago
    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 9y ago
      It is not used on objects on the java heap. No copying occurs.
      - tmyklebu 9y ago
        Aha, makes sense. (Now I see that at http://neanderthal.uncomplicate.org/articles/tutorial_native.html. http://neanderthal.uncomplicate.org/articles/tutorial_native...) 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 9y ago
        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 9y ago
        > 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 9y ago
        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.
  - ilammy 9y ago
    > Neanderthal works with barest primitives I see what you did there with naming.
_sdegutis 9y ago
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 9y ago
  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).
  - mac01021 9y ago
    Does neanderthal claim feature parity with numpy and ND4J?
    - dragandj 9y ago
      Do numpy and ND4J have feature parity among themselves?
      - mac01021 9y ago
        Fine... I should have said "or" instead of "and". I think it's a valid question in the context of your code size comparison.
        dragandj 9y ago
        The point is that Neanderthal is not a port of either of those libraries. Sure, it does not have all the features of numpy or ND4J, but, on the other hand, it has many features that neither numpy nor ND4J have.
geokon 9y ago
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 :)

118 points by dragandj 9y ago

Source code of the library is at: https://github.com/uncomplicate/neanderthal https://github.com/uncomplicate/neanderthal

Fede_V 9y ago

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 9y ago
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 9y ago
  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 9y ago
    It is not used on objects on the java heap. No copying occurs.
    - tmyklebu 9y ago
      Aha, makes sense. (Now I see that at http://neanderthal.uncomplicate.org/articles/tutorial_native.html. http://neanderthal.uncomplicate.org/articles/tutorial_native...) 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 9y ago
        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 9y ago
        > 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 9y ago
        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.
- ilammy 9y ago
  > Neanderthal works with barest primitives I see what you did there with naming.

_sdegutis 9y ago

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 9y ago
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).
- mac01021 9y ago
  Does neanderthal claim feature parity with numpy and ND4J?
  - dragandj 9y ago
    Do numpy and ND4J have feature parity among themselves?
    - mac01021 9y ago
      Fine... I should have said "or" instead of "and". I think it's a valid question in the context of your code size comparison.
      - dragandj 9y ago
        The point is that Neanderthal is not a port of either of those libraries. Sure, it does not have all the features of numpy or ND4J, but, on the other hand, it has many features that neither numpy nor ND4J have.

geokon 9y ago

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