All the Maths

@mikolalysenko

From a few conversations with @mikolalysenko I learned that languages like R and libraries like scipy which are known for amazing math are just binding to the same old Fortran libraries (and in a few cases C). That said, I don't think a binding layer will work for JavaScript, especially if we want to use them in the browser, but I do think we could use emscripten to transpile them in to JS, hopefully asm.js.

I got as far as writing most of a script to pull down calgo, which is laid out incredibly horribly, extract them all and generate package.json files which we could eventually checkin to git.

My attempts at getting emscripten to work ended at getting dragonegg to find the right llvm :(

Anyway, my instinct is that the best strategy is to:

Pull down all these ancient libraries in to a repo and check them in to git (they aren't in any source control as is).
Write scripts which can automate creating a package.json and individual npm package for each algorithm.
Write benchmarks to get the VM engineers competing for performance, hopefully ending in superior performance to the original Fortran libraries.

Thoughts?

I know that @groundwater is interested in this as I imagine @substack and @maxogden are as well. Additionally, I saw @rwaldron throw out some links on twitter which I think were related to crazy math in JS.

Oh, forgot to link to the other two big algorithm collections: lapack and arpack.

I'm unclear about the goal. Is it to provide BLAS and/or LAPACK/LINPACK bindings for people to use in their applications? Encourage the V8 people to work on numerical performance? Both?

There's a strawman for native SIMD support and I believe Intel is committed to bringing that to V8.

(Spidermonkey nightly already implements the strawman but I don't think it has made its way into a release yet. Or if it did, it's probably behind a flag.)

There's also River Trail (strawman), another Intel effort that aims to bring parallel numerical computing to JS.

Mozilla calls it PJS and I believe they have declared an intent to implement it but I don't know what the status is. This meta-bug is the tracking bug, if you're interested.

@bnoordhuis

@bnoordhuis I think all the enhancements that Mozilla and Intel have been working on will be great and I'm looking forward to them but R and scipy don't have them and are seeing adoption in one area that JS isn't at the moment. The only they have but we don't is access to these ancient libraries. If we were able to transpile them we might see some of that adoption, if they are slow the VM engineers will likely make them faster or work a little faster getting these other features in and improving emscripten or these libraries directly to use them.

The obstacles to getting high performance numerical code working in JS are relatively clear:

Lack of SIMD data types, which prevents using high performance CPU features
Lack of value types or pass-by-value semantics, which severely restricts optimizations based on cache and memory bandwidth considerations
Limited multithreading and shared memory support, which restricts the effective memory bandwidth available in JS
No viable pathway for either porting or linking to legacy code bases

The first three of these issues are really language design problems, and not things that we are going to solve at the node level. It is up to the ES committee to get these real performance impacting issues fixed, and the best we can do is put pressure on them to make it a priority over syntactic twiddling.

Still even if the ES committee never sorts these problems out, I would like to point out that at least in regards to the first 3 points JS is certainly no worse than Python, which suffers from exactly the same draw backs. And yet, people use Python extensively for numerical computing, even though an equally strong case could be made that the semantics of the language do not make it suitable for the task.

Accepting the first 3 items as a given, we are left with what I believe is the real problem here:

The main blocking issue to getting numerical computing working in JavaScript is interoperability with C and Fortran code bases.

Interoperability is really a broader problem that affects many applications (for example crypto could benefit a lot here as well), but it is felt in numerical computing especially acutely due to the size, sophistication and fundamental importance of legacy numerical code.

In my dream of how things could work, I would love to see a more modular/npm-centric solution to this problem. Something where we could just npm install lapack and then require it and call things like DGESVD directly -- without having to install or set up a bunch of external native code, using apt or whatever. The resulting bundles ought to work both in node.js and in the browser, and while we are at it on native environments there could be a fallback to the original Fortran where possible.

I think ruling out making native bindings to existing libraries sells Node short of possibilities. Especially in terms of performance. TBH I don't care about making this browser compatible, because even with all the proposed optimizations (which will take a while to come out) it'll still be much slower than using native methods.

V8 has dome some great work to make calling into C++ from JS really cheap. Here's an example that just sums an array of uint32_t.

// run.js
var runMe = require('./build/Release/addon').runMe;
var smalloc = require('smalloc');
var stypes = smalloc.Types;
var ITER = 1e5;
var SIZE = 1024;
var b = smalloc.alloc(SIZE, {}, stypes.Uint32);

for (var i = 0; i < SIZE; i++)
  b[i] = i;

function simpleSum(arr, len) {
  var ret = 0;
  for (var i = 0; i < len; i++)
    ret += arr[i];
  return ret;
}

var t = process.hrtime();
for (var i = 0; i < ITER; i++) {
  //runMe(b);
  //simpleSum(b, SIZE);
}
t = process.hrtime(t);
console.log(((t[0] * 1e9 + t[1]) / ITER).toFixed(1) + ' ns/op');

#include <v8.h>
#include <node.h>

using v8::FunctionCallbackInfo;
using v8::Handle;
using v8::Local;
using v8::Object;
using v8::Value;

void RunMe(const FunctionCallbackInfo<Value>& args) {
  Local<Object> obj = args[0].As<Object>();
  uint32_t* arr =
      static_cast<uint32_t*>(obj->GetIndexedPropertiesExternalArrayData());
  size_t len = obj->GetIndexedPropertiesExternalArrayDataLength();
  uint32_t ret = 0;

  for (size_t i = 0; i < len; i++)
    ret += arr[i];

  args.GetReturnValue().Set(ret);
}

void init(Handle<Object> exports) {
  NODE_SET_METHOD(exports, "runMe", RunMe);
}

NODE_MODULE(addon, init)

So, simpleSum() runs in 851 ns while runMe() runs in 90 ns.

Now, the math libraries themselves should be able to be done without much build pain and, TBH, if computationally complex operations are being done on my server I'll care way more about taking a minute to build a library and save the CPU cycles.

As a data scientist in learning (who use Python for that), with a huge love for JavaScript, I can name three reasons why I doubt I will ever get to use JavaScript for actual numerical computing.

It is not possible to create a Uint8Array or Buffer with more than 1024^3 elements. If you think that is a huge number, just imagine that you have 1024 picture with a 1024 x 1024 resolutions. It is really not that many pictures and not that high a resolution.
Even if that is solved, there is still the 64bit integer problem. After 4 * 1024^3 elements, accessing items beyond that index becomes unreliable. I have seen some suggestions on how to implement 64bit integers, but they all seams too obscure to use for simple indexing.
There is no operator overloading, for people to be interested in numerical computing there has to be someone to use it. Those are typically people with a more mathematical background and expect A * B to involve a matrix multiplication. You can implement it with objects and function, but it quickly becomes quite ugly to look at. mul(mul(inv(mul(t(X), X)), t(X)), y) (not numerically stable solution to OLS). As you see the distance between the arguments and the operator quickly becomes to long as opposed to np.inv(X.T * X) * X.T * y. I have seen some suggest methods, but that is not extendable and the I think the mathematical syntax is lost, making operation ordering difficult to keep track of. X.T.mul(X).inv().mul(X.T).mul(y).

@AndreasMadsen

@AndreasMadsen

It is not possible to create a Uint8Array or Buffer with more than 1024^3 elements.

After 4 * 1024^3 elements, accessing items beyond that index becomes unreliable.

For technical clarity, it's actually 1024^3 - 1 bytes. That aside, it is easily possible to attach more than that to an object but not have it accessible by array index. Slightly more complicated, it is possible to access up to 2^52 bytes by numeric index but at a performance hit.

There is no operator overloading

True, and this could be a challenge. Technically someone could create their own parser and instead of operating in native JS operations they instead pass in a string. Ugly, yes, but it could get the job done.

One reason why I implemented the smalloc module for v0.12 was for scientific computing. It's easy to integrate into a little C++ and call existing C/C++ math libraries.

I like the idea of extending the Buffer concept to implement a BigNum type, which could be further extended to implement things like Vector and Matrix types. If I recall correctly, it might even be possible to override operators at the C++ level to have the prettier syntax for working with those more complex types.

@Qard

@Qard

it might even be possible to override operators at the C++ level

Do you mean JS operators?

Yeah, I thought I remembered seeing functions to define operator behaviour
at C++ level. I could be wrong. I've never tried it myself.

@Qard

@Qard Not that I remember. @bnoordhuis You know of some such black magic?

I don't think there's any such magic. About the only thing you can do in C++ that you can't do from JS is make objects callable as functions (like how RegExp used to be.)

Indeed. You appear to be right. After some quick grepping of the source, I'm finding no such convenience. I must've been remembering some other VM. :(

It is not possible to create a Uint8Array or Buffer with more than 1024^3 elements. If you think that is a huge number, just imagine that you have 1024 picture with a 1024 x 1024 resolutions. It is really not that many pictures and not that high a resolution.

In ES6 it is possible to create a typed array with up to 2^53 - 1 elements. If this is not implemented in V8 (or perhaps, not implemented in the version of V8 that ships with Node) than it is worth filing a bug to make sure it gets fixed.

Even if that is solved, there is still the 64bit integer problem. After 4 * 1024^3 elements, accessing items beyond that index becomes unreliable.

You should be able to access up to 2^53 - 1. Not quite 2^64 - 1, but certainly a lot more than 4 * 1024^3.

There is no operator overloading, for people to be interested in numerical computing there has to be someone to use it.

This is where things get tougher, and require tricky ESfuture-track work. The current thought on all things in this design space is found in https://github.com/nikomatsakis/typed-objects-explainer/; operator overloading is the least-developed of these proposals, unfortunately. But the good news, at least, is that everyone on TC39 agrees this must be done. It's just a matter of finding the right people and time to design an operator overloading extension to the JavaScript language, and making sure it integrates with all the related topics like value types and typed objects and user-defined literals.

Lack of SIMD data types, which prevents using high performance CPU features

This is an accepted ESfuture-track proposal with good detail already, shipping in SpiderMonkey and with a patch open for V8 by Intel.

Lack of value types or pass-by-value semantics, which severely restricts optimizations based on cache and memory bandwidth considerations

(See above about value types)

Limited multithreading and shared memory support, which restricts the effective memory bandwidth available in JS

This gets tricky. We'd really like to avoid shared-memory multithreading and that whole concurrency model, with data races etc. But there is definitely room for improving the existing web worker and transferable infrastructure into something with more room for performance. I don't remember if this is exactly right, but I think at one point I was convinced you could allow a typed array (or similar) to be shared among multiple workers as long as only one of them had write access.

No viable pathway for either porting or linking to legacy code bases

Emscripten seems perfect for this?

I think the language needs infix notation, not operator overloading.

a + b // the + symbol means addition... or maybe something completly different?
a `add` b // equals to add(a, b)

@mikeal

@mikeal said:

but R and scipy don't have them and are seeing adoption in one area that JS isn't at the moment

R does generate SSE/SIMD instructions depending on the implementation of specific libraries or routines. But, primarily, R achieves parallel vectorization by means of CPU multiprocessing, and does so inherently any time R's native array processing syntax is used. Python also has a threading model that nicely supports non-I/O bound numerical calculations.

JS will need a decent threading model to support parallel vectorization at high- and low-levels.

Maybe some WebWorker + messaging-for-shared-state high-level abstractions could be provided for acheiving the heavy lifting of easy parallelism that is necessary for numerical computation.

@derekm

@derekm

R does generate SSE/SIMD instructions depending on the implementation of specific libraries or routines. But, primarily, R achieves parallel vectorization by means of CPU multiprocessing, and does so inherently any time R's native array processing syntax is used.

I am not sure what exactly you mean here, but to me all of R's array processing seems very sequential. See e.g. double binary operators here: https://github.com/wch/r-source/blob/170a2a4a030e84a2437aee1a8d15c310a3b22527/src/main/arithmetic.c#L904 and the macros here: https://github.com/wch/r-source/blob/fbf5cdf29d923395b537a9893f46af1aa75e38f3/src/include/R_ext/Itermacros.h#L32

Am I missing something? How is this inherently parallel?

@eush77

@eush77

Operator overloading is rarely useful for regular use cases, and bringing more complexity to a single language because of some edge cases doesn't make any sense.

Overloading is very useful for numerical work. It's common to want to define a custom type that behaves like a number syntactically but with special behavior - arbitrary-precision arithmetic, polynomials, dual numbers for automatic differentiation, probability distributions, placeholder variables for modeling optimization problems, etc. But it's not very pleasant to do when overloading is tacked-on to a language that is otherwise not designed for it. Using JavaScript as a compilation target for a better-suited language seems like a better approach when it comes to numerical computing in the browser. It would likely save some work to pick an LLVM-based language to start from - Julia, Rust, or maybe even Swift would be reasonable choices. I've heard GHCJS works very well which would be another option, somewhat lacking in the array libraries or ease-of-use categories though.

@derekm

Python also has a threading model that nicely supports non-I/O bound numerical calculations.

Huh? Python's threading model for numerical calculations is "do threading (and all of the actual work) in C."

@gaborcsardi

all of R's array processing seems very sequential

Any time R calls out to BLAS, you get multithreading (and SIMD) by using OpenBLAS or MKL. You can run into issues with excessive fork-join, but I understand packages like plyr are designed to decompose problems at a higher level to mitigate this.

@tkelman

@tkelman:

Any time R calls out to BLAS, you get multithreading (and SIMD) by using OpenBLAS or MKL.

Right. But R does not use BLAS for array operations. AFAIK, these are the only operations it uses BLAS/LAPACK for: https://github.com/wch/r-source/blob/d75f39d532819ccc8251f93b8ab10d5b83aac89a/src/main/names.c#L926 (the do_lapack stuff). This is SVD, a linear solver, plus some decompositions, and nothing else.

but I understand packages like plyr are designed to decompose problems at a higher level to mitigate this.

Yes, some packages work around sequential R, but you need to do extra work to use them. It's not that you just load the package and then array operations are parallel.

IMO there is nothing inherently parallel in R or its packages. (One could actually override the default array operations to make them SIMD parallel, but I don't know of any package that does this. It is too much work, I guess.)

But R does not use BLAS for array operations. AFAIK, these are the only operations it uses BLAS/LAPACK for

Ah, good to know, you know much more about the internals of R than I do. I'm a little surprised to not see gemv or gemm (matrix-vector or matrix-matrix multiplication, for those not familiar with the Fortran acronyms http://www.netlib.org/blas/blasqr.pdf) there.

One could actually override the default array operations to make them SIMD parallel, but I don't know of any package that does this.

Some packages use Rcpp to write extensions in C++ where you can do whatever you want in your own code, right? But maybe it's considered bad form to have your R package substantially change default behaviors.

Ah, good to know, you know much more about the internals of R than I do. I'm a little surprised to not see gemv or gemm (matrix-vector or matrix-matrix multiplication,

You are actually right, the above ones are only the ones called directly from R code. I missed a couple more, they are here:

~/works/r-source/src/main (trunk)$ grep F77_CALL *
array.c:        F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c:    F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
array.c:    F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c:    F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);
array.c:    F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c:    F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c:    F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
main.c:    F77_CALL(intpr)("dummy", &i, &i, &i);

So there is indeed gemm as least (if your matrices do not contain NA). Btw. as most of BLAS/LAPACK is included in R, you can call BLAS/LAPACK from packages as well, and some packages do that, e.g. there are calls from the base stats package, but this again just a handful of functions.

So the situation is slightly better than what I said above, but not much, and array operations are certainly sequential.

Some packages use Rcpp to write extensions in C++ where you can do whatever you want in your own code, right?

Right. With some limits. (Btw. you don't need Rcpp for C/C++ code in R package, base R has an API for it. Rcpp just makes it easier.)

Okay... I only recently began prototyping stuff in R, it is hard not
to notice the differences between vector operations and for loops, either
in the literature or in your own code.

The literature is a little misleading claiming things in R having
vectorizations or being vector operations or being a vectorized array
processing language. When they say vectorization they must mean in syntax
only and not in operation or processing, their claims having nothing to do
with vector processors.

Performance differences must be solely interpreted R on for-loop syntax vs.
compiled C loops on so-called "vector" syntax.

On Thursday, December 11, 2014, Gabor Csardi notifications@github.com
wrote:

Ah, good to know, you know much more about the internals of R than I do.
I'm a little surprised to not see gemv or gemm (matrix-vector or
matrix-matrix multiplication,

You are actually right, the above ones are only the ones called directly
from R code. I missed a couple more, they are here:

~/works/r-source/src/main (trunk)$ grep F77_CALL *
array.c: F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c: F77_CALL(zgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
array.c: F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc);
array.c: F77_CALL(dgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c: F77_CALL(zgemm)(transa, transb, &ncx, &ncy, &nrx, &one,
array.c: F77_CALL(dsyrk)(uplo, trans, &nr, &nc, &one, x, &nr, &zero, z, &nr);
array.c: F77_CALL(dgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c: F77_CALL(zgemm)(transa, transb, &nrx, &nry, &ncx, &one,
array.c: F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
main.c: F77_CALL(intpr)("dummy", &i, &i, &i);

So there is indeed gemm as least (if your matrices do not contain NA).
Btw. as most of BLAS/LAPACK is included in R, you can call BLAS/LAPACK from
packages as well, and some packages do that, e.g. there are calls from the
base stats package, but this again just a handful of functions.

So the situation is slightly better than what I said above, but not much,
and array operations are certainly sequential.

Some packages use Rcpp to write extensions in C++ where you can do
whatever you want in your own code, right?

Right. With some limits. (Btw. you don't need Rcpp for C/C++ code in R
package, base R has an API for it. Rcpp just makes it easier.)

—
Reply to this email directly or view it on GitHub
#1 (comment)
.

Performance differences must be solely interpreted R on for-loop syntax vs.
compiled C loops on so-called "vector" syntax.

I believe so. Plus there is the potential, that they will be parallel some time in the future.

Great to see people working to make JavaScript better for numerical computing!

I do not see any compelling argument why R, a language which I am quite familiar with as an aspiring statistician, would be better suited for numerical work than JavaScript. However, the biggest disadvantage of JS at the moment is the lack of appropriate, well maintained libraries.

For example, it would be great to have a well-implemented linear-algebra package. I do not know of any such package currently. There are a couple node modules binding to C or C++ libraries (such as the mentioned eigenjs or armadillo, which I have been working on), but they are all unfinished. And the pure JS implementations often lack a lot of functionality, e.g. matrix decompositions.

Concerning operator-overloading, one could use sweet.js macros to implement it, at the cost of having to compile one's code.

Besides for node.js, it would be nice to have access to such libraries in the brower, too. For example, I recently wanted to create a d3 visualization plotting some simulation of logistic regression results. The ability to easily carry out the model fitting inside the browser would have been of great help in this case.

So porting some of the ancient math libraries to JS via emscripten sounds like a good idea to me.
So far, I played around a bit with it and tried to port parts of the GNU Scientific Library (GSL) to JS via emscripten. While the process is trivial for functions accepting or returning standard types such as floats or doubles, the documentation is lacking in how to properly convert functions accepting or returning custom C structs (the Embind and WebIDL facilities are only tailored at C++ code).

For example, consider the following artificial example:

#include <math.h>

struct Vector3d{
  double x;
  double y;
  double z;
};

double length(struct Vector3d vec){
  double vec_sum = vec.x*vec.x + vec.y*vec.y + vec.z*vec.z;
  return sqrt(vec_sum);
}

struct Vector3d createVec(double x, double y, double z){
  struct Vector3d ret;
  ret.x = x;
  ret.y = y;
  ret.z = z;
  return ret;
}

I am very new to emscripten, so I have no idea how I could port such code as to make it callable from JS, in the sense that I want to end up with a constructor function than can create 3d vectors and a length function to calculate their length.

Has anyone figured out how to do this?

From the very small amount of time I've spent playing around with emscripten, it seems that it will compile a complete program + runtime environment into a javascript "program" that can be run in a browser. It does not have a facility for separately compiling one C/C++ function idiomatically into javascript.

One of my goals with femscripten (see above) is to make it possible (in the special case of reasonably pure functions) to produce javascript functions from C/C++/FORTRAN source, without having to link in the entire runtime+memory model that emscripten uses in order to emulate the complete program environment.

The nice thing is that emscripten does provide facilities to expose functions directly to JS. I created a git repo with the gsl code and my shell script which exports a large number of functions for use in JavaScript. Here is the link:
https://github.com/Planeshifter/gsl-js
It is also immediately possible to require the created *.js file in node, as in var gsl = require("./output.js"). This will give you access to all the exported gsl functions listed in https://github.com/Planeshifter/gsl-js/blob/master/exported_functions

I have not investigated the performance hit incurred because of the overhead of having the emscripten runtime in the back, but the claim on their website is that it should be reasonably fast.

Oh, awesome. I hadn't found that yet. I will have to try that with the FORTRAN frontend, it may save me a lot of work. Thanks for pointing it out.

Quick update: I asked on the emscripten repo about how to deal with structs, and meanwhile figured out how to export the C example I posted above via Embind: emscripten-core/emscripten#3083. So I do not think that there are any fundamental barriers preventing us from exposing entire C libraries to JS via emscripten.

I played around a bit and ended up transpiling all of the complex number functionality of the GSL library into a node package. The GitHub repo with some documentation is here: https://github.com/Planeshifter/gsl-complex-js. Installation via npm: npm install gsl-complex.
All functions are there, and using emscripten to export them ended up being straightforward (see the bottom of the complex.cpp file for the Embind bindings I had to create). I only had to create getter and setter functions for the real and imaginary part of a complex number in the original C code.

Navigation Menu

Repository navigation

All the Maths #1

The main blocking issue to getting numerical computing working in JavaScript is interoperability with C and Fortran code bases.

33 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

All the Maths #1

Description

Activity

mikeal commented on Aug 8, 2014

bnoordhuis commented on Aug 9, 2014

mikeal commented on Aug 11, 2014

mikolalysenko commented on Aug 11, 2014

The main blocking issue to getting numerical computing working in JavaScript is interoperability with C and Fortran code bases.

trevnorris commented on Aug 13, 2014

AndreasMadsen commented on Oct 7, 2014

trevnorris commented on Oct 7, 2014

Qard commented on Oct 9, 2014

trevnorris commented on Oct 9, 2014

Qard commented on Oct 9, 2014

trevnorris commented on Oct 9, 2014

bnoordhuis commented on Oct 9, 2014

Qard commented on Oct 9, 2014

domenic commented on Oct 10, 2014

domenic commented on Oct 10, 2014

33 remaining items

madbence commented on Dec 9, 2014

derekm commented on Dec 11, 2014

gaborcsardi commented on Dec 11, 2014

tkelman commented on Dec 12, 2014

gaborcsardi commented on Dec 12, 2014

tkelman commented on Dec 12, 2014

gaborcsardi commented on Dec 12, 2014

derekm commented on Dec 12, 2014

gaborcsardi commented on Dec 12, 2014

Planeshifter commented on Dec 15, 2014

smikes commented on Dec 15, 2014

Planeshifter commented on Dec 15, 2014

smikes commented on Dec 15, 2014

Planeshifter commented on Dec 17, 2014

Planeshifter commented on Dec 17, 2014

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions