Accelerating Python Code

Speaker 1 [00:03]

So this talk is about accelerating Python code. And it's more for those people of you who have a scientific background, not really a software engineering background where you're used to use different programming languages all the time anyway. So those maybe coming from mathematics or physics and generating models to tackle a complex problem, love that language, Python, because it's so simple, so easy to use, and see there are problems in getting the ideas to production code. Myself, I'm being a physicist with an applied physics background. I'm happily with Rosen for more than 20 years now, working in different roles there in development and also management. I'm a Python user since Python 1.4, which is a long time ago, and currently I'm working on embedded electronics and have a software stack running there which is Python-based. I'm a Vim user, just to get me in the right category, and I do like the print function a lot. So what's the motivation that you have here? about Python as a language that we all love for prototyping, and as I said, when it comes to getting your idea into production, there might be concerns from other people which are on the way between you and your success in doing so. Being security concerns that you just didn't do the authentication process. But all there is no matching infrastructure for the workflows, like, well, will it run on Windows then? And attitudes about the language itself, when people say you really want us to use an interpreted untyped languages that has GIL issues, well, the good news is there's low hanging fruits to help us. And we see that. So what we're covering today, what we're not covering today, there will be, there's not a machine learning talk. There will be no chat GPT or GPT whatever, no agents and no transformers, sorry for that. It's just plain Python things, numerical issues, that's all we do today. What we will be addressing is we give an overview of the possible solutions to speed up simple numerical calculations to overcome Python's main limitations right there. Prepare a strategy for scaling an execution model. And we were looking for a simple numerical formula to do that with a notable workload that we can create with that, which is well suited to try to scale on the path from very small systems to very, very large systems. For the impatient and maybe for your management, the right thing is what you can report when you get home. Everything will be just better after that. So what about theory that we use? We're using a thing called the Mandelbrot set to calculate. I'm happy that at least at the talks that I attended, nobody used that before, at least this time. And it belongs to the area of chaos and fractals. So what is a chaos? In the scientific sense, mainly stands for unpredictable behavior in systems, in complex systems. Think of the butterfly effect, probably most of you know that, that the wing of a butterfly somewhere in the world might affect the weather on the very other side of the world. That might be a very theoretical example, but it's gotten quite famous in that area. And there is fractals, and they show very complex behaviour, very complex shapes, can be found in natural structures like we've seen here in the photos, and people worked on that for a very long time to find simple mathematical formulas that help to describe that. So it's a very, very complex structure that we see right in the middle, the Romanesco, the vegetable, and if you were to describe how that looks, make a model that describes it very well, it is very, very complicated to do. And we find these structures in nature all over the place and probably there is no complex code in nature that resembles that. So there must be a much more simple solution and there is. There is one man that got famous behind it which is Benoit Mandelbrot. He is a Polish-born French-American mathematician. He studied mathematics both in France and the United States and he had a broad interest in various scientific topics if you look up his work it's lots of areas where you contributed he is a very lucky guy being employed at IBM and sometime he became a fellow at IBM which is a job every one of us probably envy is a lot because that is you get a lot of money that you can spend on things that you think are important for the company and for the future. And IBM just supports you in doing so. So that is really a thing. Because of that, he was able to get to significant computing resources quite early in his work. Where no one else had any computers at hand, he was able to do so. And he was working on the idea or to study roughness in nature. There was something really describing the think of the clouds that you see in the sky where it's just not a circle or a distorted circle, it's something very much more complex than that, or coastlines. And he came up with the idea to describe this thing as roughness and to find something easy describing that roughness. He was the one that made the phrase fractal to show these complex structures where self-similarity is a thing where you look at a large scale, look in more detail and find the same structure that you've seen in the large scale and the small scale again. But he was not the first one to do so. He was extending a work that is more than 100 years older from the French mathematics Fatou and Julien, but he studied a specific subset of what the work before was done, and that led to a set of numbers in the complex plane that was named after him, the Mandelbrot set. So what it looks like formally, we can do a formal definition of that, and the Mandelbrot Z is starting with a number zero, Z zero being zero, the next number you can calculate like squaring that number and add a number C to it. And he was using that new interpretation based on the work of Fatou and Julia, then using the coordinates and the complex plane for that specific number to add to that formula. So, all these numbers that fulfilling a special rule like those in iteration do not blow up. They do converge. All those numbers, they belong to the Mandelbrot set. So, let's have a look at some examples for that. And there It's all gone. I don't know why. Let's see if we can find another. Always have a backup. Unpredictable outcome. It's literally named like that. So back again. Well, if you want to iterate that, maybe in your head, that formula, then it says take Take a number that comes to your mind and update that number by squaring it and adding that original number to it. Repeat it a few times and see where everything goes. So let's do that. Write a little Python function that supports us, call it calculate, be a bit responsive, and all it gets are two parameters. One is that number that we are thinking of, and we will limit things a bit when we're starting it, depending on what might be going on. We start off with using that very number as x0, iterate the formula by squaring it, adding the original value, and print the outcome. And we're looking at that for some numbers that immediately come to our mind, like all those numbers, 0, 1.5, these things. And see where things really go when iterating that function. And there's things missing. Let me see if we can have that again. That's not looking too well. I don't know what happens here really. Issues in loading the data. Locally. Sorry for that. You know, if you try these things days before the talk, over and over again, hours before that things just work and then you're standing here and things no longer work for whatever reason that I don't know. Okay, that's looking better. We're coming there. There we are. So here I have some tabs that show outcomes from these numbers that come to our mind like starting with zero, what will happen then? Well, everything stays at zero. You're squaring zero, outcome will be zero, adding zero, nothing fancy happens. Then you're going for number one, that you can do in your head, you're using one, squaring it, it's one, adding one is two, next iteration will be squaring two, it's going to four, adding one, it's five, and things will explode up quite quickly, really. Next thing to show is, well, let's half it down, when we see one is larger than zero, getting any larger will probably also blow up. So let's get a little smaller than that, let's say use a half, and then yet still we see things tend to blow up again. We are at 10 to the 32 in the last number. That's a lot. Let's half it another time. 1 over 4. That looks as though we can get some conversion here. Not quite sure. Let's have a few more numbers to look at. Well, yeah, that pretty much looks as though there might be an end when repeating this. Let's get some further. Halving it again. Looking at a few more numbers. Probably it is getting somewhere around probably half in the end. Not sure. We need to repeat a lot more. Let's halve it a bit more and see whether at least probably that is going to converge as well. What about on the other side of the, of our zero point, minus dot 1.5, and that probably also is going to converge quite good, next wing, okay, will also converge, minus a half, probably also converge, we probably get the idea, minus 1, well, that's not expected, right? But if we use square in our head, things Things get clearer, like minus 1 squared, that's 1, removing 1 gets us to 0, squaring it will be 0, removing 1 is minus 1 again, and we start from scratch. That will lead us literally nowhere. And if we move further down the road, then things get even more complicated. We see some minus values here, a positive one, minus values, a positive one. That pretty much seems chaotic. And you can go on. Minus two, that will literally have an attractive point. Minus two square will be four, subtracting the two will be two again, and that will literally leave you where you are. And one step further down the road and things will explode again. So that is by playing in your head, but you still do not get the right idea how everything looks like. Let's do it graphically. The idea being to now iterate on the complete real numbers plane and see how the outcome looks like in that case. And there we see all those areas that are quite interesting, actually. We have stable points. We have a bifurcation, how it's called, that things split up, being stable, into even more branches. And it's getting chaotic afterwards, but it's getting stable again, and chaotic afterwards. That's quite interesting. But the idea that Mandelbrot had was to extend this to the complex plane and see where that leads. So let's, all we need to do is change the code so that it's no longer working on real numbers but on complex numbers. Get the imaginary part in and try your luck with that. So what about the imaginary part being zero? That gets our bifurcation diagram that we just saw. Adding an imaginary part changes things quite significantly. That's even more chaotic. It's the same things that we see going on, we have stable points, we have splitting up in branches, we have chaotic regions, quite a thing going on. It can go on on the imaginary axis, it's getting chaotic earlier, and it's getting more complex and more complex, and even we see things are totally blowing off if imaginary numbers get larger and larger. On the other side, having negative imaginary numbers, that gets us some symmetry, that's something that we've seen before, so there might be some symmetry involved. Yet we still do not have the good idea of everything looks like. So let's have a look when we do an animation on how these things develop, when we vary the imaginary numbers. Here we have colorized all the real numbers like a table of contents, it's the bifurcation diagram for the imaginary part being zero. And what really happens if you do the complex multiplications is that you have your rotation in the complex plane. So multiplication of complex numbers actually are a rotational process. And to visualize this, I've added another dimension. So like when the real axis on the x-axis, we have the real part. On the Y axis we have the outcome for the imaginary part and if you iterate all these things then we use the Z axis to show where the complex value, the imaginary value will end up. If we play that animation then we might get an idea how these things develop when changing the imaginary numbers that are involved. These are all the fancy things that you can do with I pi volume, by the way, which is quite a nice library to be using for these things and quite fast. But still when it comes to just seeing the set of numbers that are involved, We don't have the right idea. So that's we probably need another approach. We had an idea now how complex numbers involve in this iteration, how things change, but we still have not the right idea how everything looks like. So let's adjust our implementation just a bit. This is the very naive, the most naive implementation you can think of in Python. If you start doing Python, then probably that is what you would write down. At least that was almost what I've written down years and years back. So what we're doing here is we create, based on some interesting parameters, our view for the x and y axis, standing for the real and the imaginary part. This here gives us a list of numbers, pretty much like numpy linspace would generate them. And then all you need to be doing is iterating the complex plane with the imaginary parts and the real parts, building the complex number, using the C0 as the starting point and the number that we're adding in our iteration process, squaring it, adding the number, and have some conversion criteria being checked where we know things will explode. Then we know the number of iterations it took, which we can add to our row of complications ending up in our complete set of iterations, our 2D array, and then we can have a look at how this Mandelbrot set really looks like when colorizing it. For some reason I can't show the whole thing, but you get the idea. So all the pixels painted in black here, they do belong to the Mandelbrot set, and all the other colored areas just show the number of iterations it took until deciding that things will blow up when iterating. Quite nice. Complex features, fractal indeed, and a very, very simple formula behind it. on some machine that I worked on that took 4.7 seconds to compute. Well, the next thing when it comes to speeding up things, we probably want to prepare a bit and wrapping things, you know, that refactoring thing. And you see there is different purposes in the code. Let's refactor them out, give them a meaning. And one thing that you can do is immediately see that the centre part, really the iterative process can be directly wrapped into a function. And then apply the rest. And that took 3.6 seconds. That's interesting. It's already speed up. Anyone knows why? That's right. Exactly that's the thing. Beforehand everything was in the global context. And if the Python interpreter starts to interpret our script, then he puts all the variables in the global context and needs to find them there. And if there's a function, he starts to byte compile it, putting all the variables that are needed or addressed there into the byte compile, and he knows that they are in that scope, no longer needs to search for them in the global scope, which takes some time. So that's already a first step that everyone can do. When you're used to do scripting, use functions. And you will have a significant, a really significant increase of speed. Well, the next thing, we have a look at some implementation, so we prepare some benchmarking scenario where we wrap all the code that is necessary into one single function that provides us an API. So, we know whenever we need to calculate the model broad set, we have some bounding box described by this X min, X max, Y min, Y max. There is a parameter for max iterations. We say that's where we stop. And there will be a resolution where we say how many pixels we like to calculate. that we will see more often throughout this talk. The next thing that we can imagine is and there was the trick in the old days that everyone did, is can we benefit from symmetry? We've seen that mirroring reflection, the top part was exactly like the bottom part. And that's probably something that we can address in Python quite easily and quite elegant actually by just using our func tools cache decorator. And all we need to do here to benefit from this is in our main function that is calculating everything, wrap it by using the cache decorator and afterwards using the absolute value of our imaginary part. By doing that the effect is that when this set of parameters x, y and max iterations occurs the first time time, then really the function is executed and we have the result, but it's also stored in the cache. And as soon as it is retrieved the second time, then no longer is the function calculated, but just the result that's already known is retrieved from the cache. So that should give us a significant speed up by a factor of probably two, right? So let's see we are. It's a lot quicker, but it's not really what we expected. It's not half that. So any idea why that is the case? There's some sort of cache overhead. It's not the main problem here. There's an issue. There's a hint down below. And Miro talked about it in the talk before. It's floating point precision issues. Let's use some few numbers to illustrate that. The way that I create the axis, like numpy linspace does it, is like so. And if we get some numbers, they look very symmetric, minus 1.25, plus 1.25, but already the next number just matches at the beginning. At the very end, they no longer match. So it's not a parameter set that you will find in the cache, thus it will be recalculated. It's exactly what's happening. We have a lot of cache misses with that. Luckily, we can correct it quite easily. And the idea is to just reverse the list of numbers that we have, invert it, and add them all together, averaging them. So if there's an error in floating point numbers involved, it will be symmetrical everywhere. And if we do that, then we are correctly benefiting from our symmetry. You see that I put the trick in here. And then we are at 2.1 seconds. That's almost half what we spent before. So that's the best thing we could do for naive Python code. That's specialized. There are ways to generalize this. A bit more complex. it will be working. So I've made a benchmarking engine for this, so that we can try all our applications with different views, because maybe the results are just depending on the workload, so that's the idea. And there is a GitHub repository showing all these things, we'll post about that later on, so you can see all the implementations in detail. What What I've done is running it on a machine, all these different views, being very simple or more complex, using all these engines as a plug-in. And here are the results for the Naive implementations and the wraps. And because of the nature of the benchmarking engine using a function approach, we no longer have this naïve global based implementation effect anymore, so that's why the naïve inner func and outer func thing are literally the same numbers. But we see that naïve caching and optimised caching get us down, and the base view is the main Mandelbrot set view that we've seen before, and there are more complex scenarios like a medium view or a a slide view where things get different. We will see the results for other implementations later on. So for optimised implementations, next thing that you probably want to do is, as we have a numerical problem, use NumPy. Things change drastically when doing so. So NumPy is a classical library of choice to perform the array-based calculations, it gets you to the idea to no longer work on a single number, one by one, but using the complete array, put it in our iteration process and see where things go. It uses optimised libraries behind the scenes in the older days and on some machines architecture still, they are Fortran-based libraries like Open Blast, for example, and on Intel machines you can use the math kernel library which might be a bit faster than that. So how do things change? So now we have the numpy linspace that you're probably all familiar with to generate our axis and we create the complex array right from the beginning by using our real numbers and adding and switching the direction of our imaginary numbers, multiplying by our complex number, and then we have a field of complex numbers that covers the complete interesting complex plane. And then we have the inner function, which pretty much looks like the Naive implementation before, but is working on a complete array. And now, as we've seen that there are numbers that explode quite quickly and will give you overflow errors, we need a strategy just to work on the necessary subset. And that's what we're doing by a binary mask, we see where the convergence criteria is fulfilled, that makes our to-do list, and then we are just working on the subset of the array, iterating our formula. And provide back a complete iterations array afterwards for colorization. Will that work? Probably it will. That gives us the complete iteration process again. So, the benchmark results, they change, and it roughly gives us the speedup of a factor of two with not really much change in it. You will have more benefit if things get larger, more complex than that. Next thing that comes to mind is, and that is a very modern approach, is jitting. So, Can we JIT it? Yes, we can. Probably can we JIT it. There are a number of things that you can use for that. One that immediately should come to your mind is using PyPy. And that is a different implementation. It is an implementation of the Python language in a different way. Like CPython is the implementation that we use to use for Python in C or Jython is a variant where everything has been done in Java, PyPy is a variant implementing it, the Python language and using a JIT approach. And the JIT approach is the idea behind this is to, when the code is executed, get an idea what's going on and using optimized strategies instead. Another approach is using Numba. And Numba can just be used to decorate functions that are there. It has a tracing JIT so it will also have a look at what code is executed and what is it aimed for. It will see where things might be parallelised and use optimised machine code for that. It will use the code blocks that a C compiler would be using to compile C code and instead use that for running your Python code. And things do not change a lot. So this is the NAIF implementation that we had already. With the three nested follows. This ugly thing. And all you need to do to accelerate it is to decorate it using number JIT and that will give you a significant speed up as we see in a moment. Another idea is to start to think in kernels like switching the problem around, no longer having the idea to iterate a whole array of things but to have something that can be applied to each and every element of the array. So you write a kernel, we have to be a bit more specific here telling Numba about the data types that are involved, and then it can compile for you a very optimised version that you can use. Another new kid in town is Jaxx, pretty much aiming for the same thing, and it belongs to the tool set of deep learning and machine learning probably, but as we know for all of this, we do need numerics, so we probably want to have a look at that as well. I personally wasn't able to get it running on the machine, so that's why sad that things will be missing in the benchmark, but things are simple because all the libraries that are evolving today will use the API that most of us are familiar with already based on NumPy. So that's already the best idea to look for those libraries that implement it because the change will be very, very subtle as we see. It's nothing really different here but just one line and we could have done that in NumPy as well because JAX doesn't allow in place modifications. So we have to have a copy but that is very, very optimised. And the benchmark results, they are stunning. The use of PyPy and NumPy and Numba and Numba vectorised, they bring you significantly down to our reference for the naive implementation around three seconds. So, we are now at half a second for the same thing. And that does not vary a lot with the workload. So, the points you see in the diagrams, they hardly move but the scales on the axis, they do. And there are all of core implementations and I have to be quick with what I'm saying because we are lacking time a lot. The idea being that you have multiple cores available in your machines these days and Python being stuck with the global interpreter lock hardly can use that. But there's ways to overcome this by applying the divide and conquer principle so that you can split up your workload and our tool, our formula that we use is ideal for that because no pixel is dependent on the outcome on the other from any neighbouring pixel so that's a good thing and you simply split up the work and have a sub process running for you and there are various methods in the Python standard library that allow that, one being the oldest being multiprocessing, so you have a multiprocessing pool that helps you doing that, and a more modern approach which is concurrent futures, where you can choose between a thread pool and a process pool. And if you look at the results in the end, then you see a pretty significant speed speed up by that as well. Again, we have here the naive code being at three seconds, and it can do all sorts of combinations of things that we've seen before to just let it run in parallel on the number of calls that you have, and it will give you a significant feedback speedup. Many libraries allow you to execute that code also on the GPU, which is massively parallel in that case. Same strategy, you need to think of getting to the idea to have small functions working on a subset of the data that you have. QPy is the thing that is probably the best to start for NVIDIA GPUs, it's from the context of NVIDIA, and it resembles the Numpy interface again, so it's a very, very slight change, the whole code just doesn't change the namespace and then you can start off accelerating a lot. Numba can do the same thing for you. You have a target for that that allows you to compile the code for the GPU from numpy import CUDA, CUDA JIT and then you have a kernel that you can apply and will give you a significant speed up. And there's things like PyTorch that again use the same API, it's a pretty good start, nothing really changes and you get there. And the benchmark results, they are stunning. That gives you to near zero seconds. It's just and it happens. Usually the first calculation that you're starting, they are costly, but after that things will get a lot better. So if you have repetitive workloads to do, then you will benefit strongly from that. And if you're used to using the concurrent futures interface, then it's just a very, very small step to use many computers for that, like the distributed approach to do things. And there, the newest kid in town probably is Dask, and that's what I've brought with me today, like a very small Raspberry Pi base cluster, which is doing all these things. You might have a look at it later on. And you can start using Dask on your local machine. It also follows the concurrent futures interface, the idea being to split your work into pieces, map your function performing the work, and putting all things together afterwards. And then you can, if that works, scale out, use many clusters. You have to start your cluster by a helper function. You have to prepare it like your secure shell based login should work, you should use the same architecture on these things, but then fun begins, all you need to change from going local to remotely is to use that scheduler address, and then you're good to go. And then you can enjoy the wonderful DAST dashboard and follow what work has been done for you all day. And that is really literally what you can stay all day. It's very, very satisfying to do so. So that brings me sadly to the end. Time is left. So what is to say? What is the idea? I can think of three scenarios, I think, where you can look at one being you have your local machine and it should accelerate on that. Then go for concurrent futures. Learn how you can divide and conquer and have some fun. thing is the GPU, and then you can benefit that most libraries today resemble the same interface and make it easy for you to switch. Once you're used to using NumPy, it will be easy to use the rest. And then there's a thing where you have many nodes around, you can persuade your management to buy in, and then probably Dask is the first thing to go because you came from that strategy and going to concurrent futures and then you can scale out. And that's That's all I've got to say today, so thank you very much for listening, being here. And I'm open for questions and probably and hopefully have answers.

Speaker 2 [38:56]

Couple of questions actually as a starter what tool

Speaker 1 [39:02]

I was hoping someone was asking for that. That is Jupyter book. So the idea being to use a missed markdown. When you're reading docs, like literally read the docs, then mostly you see a structure like that. With all the table of contents in there and the pages there. All I need to do is have a little CSS to make the fonts larger and make it look like slides. hopefully people will have an eye on that and support this more natively in the future. That would be cool.

Speaker 2 [39:37]

So, the next question would be what are some examples that you could use the idea of fractals to solve some real world problems.

Speaker 1 [39:46]

No. Not a single one. And the idea is that this is a very simple way to generate arbitrary workloads easily. And that is, for me, it's an ideal candidate for exploring the properties that a system like the Raspberry Pi cluster has. When you have an idea how the formula works and how the workload on a single machine is, things change drastically if you go for distributed approaches like this because you have the network bandwidth involved, all the communication and management over it involved, and you need to play with these things to get an idea what the best setup would be. If you do this with production code and strategies, it might be two or three steps too far already. This is a very manageable and easy setup, So that's where I started.

Speaker 2 [40:42]

Okay, thank you. The next one would be how do you figure out that you're being naive and miss something? And how do you balance that with over complex over? complexification

Speaker 1 [40:53]

Say it again, and do I figure out how being naive?

Speaker 2 [40:57]

How do you know when you're being naive about your implementation, and how do you find the right balance of accelerating your code and not over-complexify it?

Speaker 1 [41:07]

So, first of all, I love naive code, because it is so easy to follow. You can literally read English text in Python to see what's going on. And for many people, it's like code smell, if they see that, because they know better ways how to do things. But if you start off with the most efficient implementations first, you get a hard time in debugging, probably. So I tend to start with many for loops to see what's going on. There's a very simple way to accelerate if that's necessary. Actually, number has the best options on naive code. It will be the easiest. It's hard for number to accelerate efficient solutions like using NumPy, for example. But NumPy is the next step that I I use because everything else resembles the API. So if it comes from for numerical calculations, that's my way.

Speaker 2 [42:11]

Okay. Thank you. One question would be whether the notebooks are accessible after your

Speaker 1 [42:17]

They will. I will put everything in a GitHub repository. I'll find a way to post where to find everything. Because there's lots more here that we didn't have the time to cover.

Speaker 2 [42:28]

Okay, one more. Are non-NVIDIA GPUs also supported by some libraries?

Speaker 1 [42:34]

Actually, yes. There is things going on for Apple, for example, and the AMD chips and Apple Silicon with the metal backend. It's supported by some libraries. But for reasons that I don't really understand, most things really focus on NVIDIA. That's rather sad, actually. Yeah, I agree. For many reasons.

Speaker 2 [43:02]

And the last one very inspiring. Do you know polars and have you tried it?

Speaker 1 [43:09]

Yeah, sure. It is a complete different beast, I think. Polars, for me at least, I see it in the context of pandas and replacing pandas, so if I organize data in a table format, that's pretty much the thing to go forward. Not really in that context here.

Speaker 2 [43:31]

Okay, thank you so much again for the talk.