spacedome

A small writeup on the Fast Fourier Transform and implementing it in functional languages.

Thanks for writing this up, Chris!

I took a break from Julia a year or two ago because of some of these issues, one of the big ones being I didn't want to write and maintain a set of non-allocating LAPACK wrappers for iterative solvers, but the memory churn was killing my performance. So, so glad FastLapackInterface and LinerarSolve are a thing now, and the MKL situation is much easier with this trampoline development, makes me want to start working on Julia solvers again.

It does feel difficult to write performant Julia if you don't put a lot of effort to stay "in the know" as a lot of this knowledge is very dispersed, but I guess it makes sense as the language is still changing quite rapidly.

I love Julia and choose to work in it almost exclusively, but I agree with the points in the article. I've run into a lot of issues just writing numerical linear algebra type algorithms.

Even core, and not quite core but maintained by core dev, libraries like Distributed.jl and IterativeSolvers.jl can feel pretty rough. For example IterativeSolvers has had strange type issues and not allowed multiple right hand sides for linear solves, for years, afaik due to some aspects of the type system and some indecision in the linalg interface. DistributedArrays still is very poorly documented and looks like it hasn't been touched in 3 years.

I've run into problems when I need more explicit memory management, for example none of the BLAS/LAPACK routines have interfaces for the work arrays, so you either get reallocation or have to rewrite the ccall wrapper yourself. It can also be hard to tell where the memory allocation is happening.

My most recent problem had been with Distributed and DistributedArrays, where everything is fine if you just want a basic parallel mapreduce, but has been a huge pain past that. It's not even clear to me if Distributed/DistributedArrays has been more or less abandoned in favor of MPI.jl, which for me removes most of the benefit of writing in julia, since you then have to run it through MPI. There is an MPI sort of interface for DistributedArrays but that part is not well documented and looks like more of an afterthought.

My use case isn't even that complex, I just want to persistantly store some matrices across the nodes, run some linear algebra routines on them every iteration and send an update across the nodes, then collect at the end. If anyone has any idea how to do this correctly in Distributed or DistributedArrays or can point me to some examples that would be amazing because it has been taking me forever to piece it together.

Not going to stop using Julia but there are many basic things even just in a scientific computing workflow that still feel like they were rushed and they can really take the wind out of your sails.

Why not just use fixed point? That's specifying how much precision you need in a fixed way and has the benefit of being quite simple

You can do operator overloading in fortran to add a new number type, don't remember the details but I wrote a quaternion library once.

I meant doing day trades of long term puts, but yes I agree any money here is a huge gamble.

It looks good, I'll definitely give it a real try if I start trading options more, the order flow data alone looks like it would make it worthwhile. Do you incorporate L2 data? Couldn't find that anywhere, only thing that seems like its missing.

You can make big money trading the volitility, if you get really lucky. The price of some puts I looked at went up 500%+ during the drop today. Would not try this personally lol

Looking right now some of the options a few months out have implied volitility of 1000%+, your best bet might be selling them off on big crashes instead of waiting it out? Some of the puts I looked at went up 500%+ after this dip

What kind of premiums did you pay on those puts though? Last I looked the prices were pretty insane.

What do you think of tradytics? Looks really interesting, though I mainly trade OTC, which I sort of doubt it has much support for.

You absolutely can use regular jupyter notebooks for julia! Pluto has some advantages, like being stored as a normal julia file. The julia startup time issues affect both.

You are making a lot of ontological and epistemological assumptions that are contentious in the philosophy of math. Not saying you are wrong in thinking this, metaphysical questions don't necessarily have answers, but many would not agree with you.

I agree this is a common sentiment among mathematicians, but this is a very modern perspective. If you look back 100 years ago to Hilbert, there was less distinction between physicists and mathematicians, much less the pure/applied rift that now exists. Arnol'd (who is referenced above) was one of the mathematicians who tried to keep this unity alive.

As an 'academic' who has done plenty of physical labor, I find this argument reductive and offensive. You can disagree with the author without painting this negative picture of them.

On demand printing has also been horrible for textbooks/monographs, Springer being one of the worst. Very few copies are printed, but they serve as an important means of preserving this knowledge, and books that do not survive a single reading do not instill confidence in the longevity of cheap print on demand. They aren't any less expensive now either. I also buy most books used, often to avoid these terrible newer printings.

It is entirely possible, but is not trivial, especially for the user who then needs to know "arcane knowledge" of BLAS/LAPACK work array sizes and flags. There was some discussion about this on github, but it sort of trailed off without a real resolution. I think it is considered too complicated/niche to be in base, and was recommended to be an external library, but nobody (myself included) seems particularly interested in what amounts to maintaining a fork of the entire BLAS package. The base devs would have more insight, but this is my view from the outside at least.

Yes, many of these are "in-place" but will still allocate. I typically use the geev!/ggev!/geevx! routines, if you look at the source code you will see that the work arrays are still allocated inside the call. The in-place here (unfortunately) means only that the input is overwritten, not that there is no allocation.

Higher level BLAS operations, such as solving a linear system, or computing svd/eigen, cannot be done in-place the way matrix multiplication can, and require additional memory of a predetermined, fixed size, called the work array. This cannot be pre-allocated in Julia, as there is no interface to do so in LinearAlgebra, so these BLAS calls will always allocate memory.

As much as I like Julia, I think "trivial to write allocation free code" is a bit of an overstatement. Depending on what you are doing, it can be difficult, for example iteratively calling any of the LinearAlgebra methods, since there is no interface for preallocating work arrays (doing the ccall on BLAS yourself is not a fun work-around). It is also not always clear why something is allocating, even with the debugging tools.

This is good to know, thanks Chris. I mostly solve sparse PDEs, being able to always use LU makes everything much simpler, especially with quaternions.

I found this approach to work well for linear systems, here is some rough code I used for the representation map (note the jmag/kmag functions are part of my implementation, not sure what the equivalent is for Quaternions.jl).

  function cmatrix(Q::AbstractVecOrMat{Quaternion{T}}) where {T}
      [complex.( real.(Q),  imag.(Q)) complex.( jmag.(Q), kmag.(Q));
       complex.(-jmag.(Q),  kmag.(Q)) complex.( real.(Q), -imag.(Q))]
  end

  function qmatrix(C::AbstractMatrix{Complex{T}}) where {T}
      n, m = size(C)
      quat.(C[1:n÷2, 1:m÷2], C[1:n÷2, m÷2+1:m])
  end

My experience with various linear solvers (which a diffeq solver typicaly relies on) is that if they assume commutivity anywhere at all, which they often do, they almost definitely will not work for quaternions. Even if a derivation of the algorithm can be done with non-commutivity in mind, the implementation typically is not. Many NLA solvers are based on orthogonality transforms, which do not translate directly to quaternions, and even solvers which use only inner products and mat-vec multipication like BiCG do not really work as-is for quaternion matrices. Luckily dense linear solves with LU are fine! One can usually use the complex matrix expansion for quaternions and solve that instead without much difficulty though, and for computing eigenvalues this is perhaps even preferable, as it gives a canonical representation. I have not tried to use DifferentialEquations.jl for quaternion problems, so maybe they have figured some of it out, but it is non-trivial.

I think the functionality of Julia's sparse arrays is mostly on par with scipy.sparse, it's just a bit rough around some edges still, and spread out into non-base packages (e.g. into IterativeSolvers.jl, ARPACK.jl, MUMPS(*).jl etc) of various maturity. Having the sparse arrays built-in might (hopefully) lead to a large ecosystem of interoperable sparse packages, but as with much of the Julia ecosystem this is still a work in progress.

Personally, working on sparse eigensolvers, my only complaint is that IterativeSolvers.jl is still somewhat lacking (though the ARPACK.jl bindings are more comparable to what scipy does), but otherwise cannot imagine going back to python for this kind of low-level numerical research. Having code run comparably to fortran speed (barring some memory annoyances) is huge for working on numerical algorithms.

I have only seen this happen with packages that depend on non-julia libraries, such as ARPACK, but the move to providing binaries with BinaryBuilder should fix this.

The person I was responding to seemed to have read my comment and thought I had an issue with raw performance. My point above is that writing iterative code that makes many LAPACK calls becomes difficult to write in Julia because there is no way to manage the memory in this situation other than ccall-ing everything yourself, at which point I would rather write FORTRAN. I work on eigenvalue solvers, so it is all more or less just wrapping a bunch of calls to BLAS/LAPACK. As you say, the main bottleneck is in the BLAS calls anyways, but the excess allocation in Julia can really slow you down. ARPACK is maybe a bad example because its all just mat-vec, when you need to do something like compute an SVD every iteration, thats where you run into issues.

This is great, a lot of the performance in BLAS comes from memory management. This does not solve my problem though. Sometimes you want to manually control memory, and Julia does not make it easy to do so. In particular when interfacing with BLAS/LAPACK though Julia, memory management is quite ugly, mainly due to the need for allocating work arrays, and a pure Julia implementation of BLAS is unlikely to fix this. I don't think writing eg ARPACK performantly in Julia is imposible, just painful to the point that writing it in FORTRAN starts to make sense (to me).

They are solving the linear system Ax=b, written in Julia/Matlab usually as A\b, specifically for a fixed size 5x5 matrix. The Grassmann algebras can give a way of representing the matrix and vectors, and a comparison is being made to the representation (and solver) given by StaticArrays.jl, which just uses typical matrices and arrays but optimized for small sizes. No real motivation to use Grassmann algebras if your only concern is generally solving a linear system, other than this method seeming to be faster, which mostly (imo) points to the StaticArrays.jl method needing to be optimized.

Edit: This post also does not mention numerical stability of the Grassmann.jl method, which is a real concern in practice, even for such small systems.

This has been my general experience writing numerical linear algebra methods in Julia. Optimizing a specific algorithm, eliminating unnecessary allocations is the first thing I do, and can give large performance gains, especially for iterative methods.

It can be a bit annoying, and can result in much less readable code (eg having to explicitly write things like mul!(C, A, B)). It ends up looking a lot more like the FORTRAN it is meant to replace, and if you aren't careful you lose the ability to use generic types, though multiple dispatch is a great solution. The worst case I have found is iteratively calling linear algebra solvers (Eig, SVD, LU) which do not have any option for preallocating and reusing work arrays. The only way to do this now is to call the BLAS/LAPACK routines directly with ccall, which is a huge pain. There was an attempt to fix this with PowerLAPACK.jl, but it seems abandoned, so for the moment writing optimized methods in FORTRAN/C and calling from Julia is sometimes still preferable.

Julia does quite well though, so this only seems necessary for things like core NLA libraries, for example I would not try rewriting ARPACK in Julia for a performance gain. The gain in flexibility for writing these in Julia is absolutely huge though, and I definitely recommend it. The Grassmann.jl library has many examples of what a language like Julia makes possible, or DifferentialEquations.jl.

Worth looking at also is Approximation Theory and Approximation Practice, which is centered around Chebfun. Almost everything by Trefethen is unusually well written, Spectra and Psuedospectra is one of my favorites.