Jared Hoberock

What say you, Jared? It's always difficult at first glance to find the parallelism in a problem. I will say that. And it's um it usually comes as a surprise when you figure out the the solution to the puzzle.

Bryce

Yeah. That that is how I often feel. I feel like I keep telling myself that like we keep we've talked regularly on this podcast about how we so often reach for scans or reductions or transforms. Like those are like, I think in my mind, the three basis operation is scan, reduce, and transform.

Conor

Welcome to ADSP the podcast, episode 241, recorded on May 21st, 2025. My name is Connor, and today with my co-host Bryce, we continue part three of our four-part chat with Jared Hobrock. In this episode, we chat about a tricky parallel scam problem, parallel rotate, and more.

Bryce

Alright, so as I mentioned, as I mentioned, NumPy does have some unfortunate warts. Like there's a lot of weird historic quirks in how a bunch of the APIs work. And um, and you know, it's like object model is not great, where a lot of things give you back copies, and then over the years they tried to make some things give you back views. But so you never really know with a numpy API whether you're getting back a view or a copy. But that is not the thing that I'm gonna pick on today. I I'm gonna I'm gonna pick on one very particular thing that's gonna be near and dear to our hearts, and I get to I get to uniquely pick on both NumPy and something that Connor used to work on, which is Rapids, specifically QDF. So one of the first things that I did when I started writing more Python code, uh, which was like a year or so ago, I've I've always used Python a little bit, uh, but in a very crude way, but I started using it in earnest, you know, one to two years ago. One of the first things I tried to do was to write a user-defined reduction. How do you think I went about doing that, Connor?

Conor

You read a Mark Harris blog and you implemented that?

Bryce

No, no, no, no. Just like imagine that I have a numpy array. What do you think I tried to do to write a user-defined reduction?

Conor

Wait, in numpy?

Bryce

Yeah.

Conor

Uh you looked for a generic reducer fold.

Bryce

That's correct. I tried to do np.reduce and then pass it a lambda. And you know what happened? It didn't work. Because there is no such thing there. And the as I learned more about like, you know, sort of the numpy and python mindset, it it uh I discovered that the thing, the reason that this doesn't exist is because in this sort of numpy worldview, as far as my best guess is that in this numpy worldview, uh where you've got this ranked polymorphism, if you needed to write a user-defined operation, you should really write like a UFunc, which is a UFunc is a function that uh takes things that seem like scalars uh but then can be uh applied to uh a you know to to tensors. Um and that like there's a particular way of doing this application uh uh with a ufunk that makes sense for reductions. Um however I I I always find this answer like very unappealing, that like I just want to be able to like call dot reduce in place. And also I think for parallelism, uh uh having to write your reduction as like a uh a U func um is like a UFunc is a much more general form than writing a s a purely scalar reduction operation. And so for like serial code, maybe it doesn't matter so much, but if you're if you've got a a parallel array API, uh you would prefer to have the much more closed form of hey, you just give me purely a scalar thing, a purely a scalar operation, and then I do the application of this operation across this whole array. So, anyways, NumPy does not have a user-defined reduction uh function, does not have a user-defined scan function, and but because of that, uh pandas never had these things either. And then because pandas didn't have them, pandas is this data frame library. Uh for people who don't know what a data frame is, it's just like a it's just like a 2D, it's a fancy 2D matrix. Okay. Uh it's like a spreadsheet, but like it's really just a fancy 2D matrix. Um maybe that's an unfair classification. But, anyways, that pandas is how you make these fancy 2D matrices for dealing with tabular data.

Conor

The primary difference between a numpy array and a pandas data frame is that your columns can have different types. Uh numpy matrix all has the same type, whereas a data frame, the columns can have different types, which makes it much more amenable to I don't know, data analytics.

Bryce

Yeah. So um so because pandas didn't have it, then QDF doesn't have it. And do you know Bradley Dice, Connor?

Conor

Of course. He gave me a very nice jigglypuff uh sticker when I left the Rapids team.

Bryce

Okay, so good. Um I was at PyCon last week, and me and Bradley were talking uh about a particular interesting problem. And the the basics of the problem was hey, I've got um this uh this data, this like uh a time series data, um, and I want to apply like I have it the data in one column into time series, and I want to apply this like recurrence relation uh to compute some uh uh some value that depends upon the values in this this time series. Specifically, I'm gonna share my screen now. I know, I'm sorry. I'll try to try my best to use words here. Specifically the operation that uh here, let me. Alright, now I'm sharing my screen. Specifically the operation that uh that that I wanted to apply. So uh uh if I'm doing something where I'm computing, you know, a sum for each value in this time series database, then it should be a scan, right? And and essentially over the course of three days, Bradley and I were discussing whether or not this scan could be paralyzed. So this is what the operation looks like. So it's a it's a it's a you know a reduction to style operation where it takes two inputs, and like if the two inputs I have written here are O and I. Um and so the the result of this operation is O plus K times I minus zero uh minus uh O. And so this is another way of expressing it as a recurrence relation here where it's like XI equals XI minus one plus uh ki times XI minus x i minus one.

Conor

Alright, so wait, just just for the listener, because the listener got lost there, guaranteed, because I got lost looking at this code based on your explanation. The relation here is between two adjacent elements. Uh yes and no. Um it says Xi minus XI minus one.

Bryce

You have to think about this in the sense of an accumulation here, where um the the left hand side is the already is the the already computed value, right? Like the left hand side is the value of like you're computing this this new column, right? The left-hand side here is the value that's coming from this like this column. And the the the right hand side is the value from the the other column that you're that you're using as an input.

Jared Hoberock

Uh are you following Jared? Because I'm lost. It looks like some set uh it looks like the kind of thing you would write down if you were trying to figure out how to implement a prefix sum. That's what it looks like to me. Is that not is that not what it is?

Bryce

Yeah, yeah, no, it is.

Conor

It does look like something that is the input to a scam, but I just said you're looking at a uh or so I mean it if it's an if it's the input to a scan, I should see an accumulator and a current element. What I see is XI and XI minus one. Um and then on the top line I see own I. Spell this out to us.

Bryce

Let's let's leave that one alone. Let's just focus on the this binary operation here. This is a binary operation, right?

Conor

Okay, the op o i, where o and i are the two inputs.

Bryce

This is an operator that you could presumably pass as an operator to a scan. Yeah, but uh uh a scan a scan just takes a binary operation. An operation that takes two two things. This is a binary operation. Yes?

Conor

I mean it's definitely a binary operation. Whether it is suitable for a parallel scan is unclear.

Bryce

All right. Is this binary operation? So so we understand what this binary operation is. It's it's you know, it's uh uh it takes two inputs, o and i, and it returns you o plus k times i minus o. So is this operation associative?

Jared Hoberock

I don't remember the rules.

Conor

I was gonna say it's it's entirely unit's entirely unclear to me whether it is.

Bryce

Okay, well, it's simple enough to test, right? So uh associated.

Conor

It's not simple enough because he copy and pasted some like 56 character uh duplicate lines of some expansion of this formula. It's simple enough to test. Yeah. Says the guy who couldn't type it out and had to copy and paste this stuff into his notepad.

Bryce

I typed it out live when I was talking to Jake earlier. You want it? Fine, boom, we'll type it out live. All right, we're gonna type it out live.

Conor

I don't want to type it out live. Copy and paste it and walk us through it.

Bryce

Okay, so how do we test for associativity? Well, associativity means that if I do op of a of or sorry, let me let me let me phrase it another way. Associativity means that A plus parens B plus C is the same as parens B plus C. Or sorry, is this I'm gonna start that over. I mean, as parentheses don't matter.

Conor

I'm gonna leave this in to make sure that you are on record as just absolutely destroying the understandability of this topic.

Bryce

Okay, so A plus open parens B plus C close parens has to be equal to open parens A plus B close parens plus C for an associative operator. Was that clear?

Jared Hoberock

Associativity means, yeah, you can you can um move the parentheses around and you get the same result.

Bryce

And the reason why this is useful for parallelism is because when you can do this, when you can arbitrarily like move the parentheses around, this allows you to essentially arbitrarily introduce intermediate values. So what this means is if you you have an operator that's associative, it means that you can have one thread go and apply this operator to four values, producing a result for that thread. You can have another thread go and apply the operator to its own four values, producing a result. And then you can take those two aggregates from each thread and you can apply the operator to combine them. And that's why we need associativity for parallelism. Unfortunately, however, in the case of this particular operation, as we can see, if we expand out those the the the two sides of this equation of the the op of a comma opc and op of op AB, comma C, this expands out to some you know complex thing which is not going to be equal to each other. And and so that means that this is not associative. And so what do we do? Do we just give up on parallelism? Well, I started thinking because everybody else in the room at the time with me, including my dear colleague Bradley, was like, no, can't be done. Can't paralyze this.

Jared Hoberock

Is it just the minus, is the minus sign is the problem?

Bryce

The minus sign does get to be the problem here, um, ultimately. Like if you think about like like what where does the non-associative part come from? It does, I believe, fundamentally come from the minus. Because if you just had pluses here, um then, you know, because addition is is associative, I think you'd be fine in that case. So I I I you like the best way to distract me for a week is to tell me that you think that something can't be uh paralyzed. So I I I started thinking more about this, and then I came up with the idea of okay, what if we could separate this and we could decompose this operator into two operators uh uh into like a into like an associative part and a non-associative part. And that is where we start going down the rabbit hole. Uh let's see. Um so I started talking to Chat GPT about this. Um and uh we we we came up with this idea of using an affine map here. So we we decompose this operation. So first of all, you you can rewrite this operation as being uh like parens one minus k close parens times z plus k times x, uh, which is what it's done here. And and and so it rewrites it in that form, and then the idea is that uh we can separate uh this operator into uh uh into two parts, and that we can we can represent this as an affine map. So this little two by two matrix here, uh, and and we represent these terms to one minus k and to k. And it's not important that us or the the the audience necessarily understand the specifics of the math here. The key thing is that you take your operation and you go and you look at if I apply that law of associativity and it fails, okay, what is the difference? Like what is like when I when I do that A of B of like A op B op C versus the A op B op C. When I do that comparison and they fail, I just subtract those two out, and then I can figure out what's the correction that's needed. And then I can I can split this operator into an associative component and a non-associative component. And then what I can do is going over here to where we have some code, I can do I can do my operation, in this case a scan, and what I'm gonna do is I'm gonna call scan with a transform iterator as the input, and then a special scan operator, and then a transform output iterator for the output. You with me so far? Connor?

Conor

Um I'm hanging on for dear life. I'm not sure if I'm totally with you, but uh I've I've been following for the to the extent that one can say, I think I'm following.

Bryce

Yeah. So what what do these transform iterators, uh, transform iterator and transform output iterator do? The basic idea here is that you you take your your problem and you you figure out some associative space. Specifically, you find some monoid that you can linearly transform your operator to this monoid space, and then later you can transform it back. And in particular, you want ideally you want a form where the the representation of this monoid, of this the function composition of this monoid, has a very small state. Like you one way you could do this is like for for you could scan where for every element you just built up the expression that you would apply, but that would be very inefficient. But in this particular case, where it's just this linear like recurrence relation, uh the only state that you need is two elements of this little two by two matrix. And so our encode function is this little make affine function here, which takes like a uh a single number a float and transforms it uh into uh uh the uh like the two elements of this matrix that we need to uh store, which is one minus k and then k times x. And then we have our compose operation, uh, which is gonna do some math, some associative math in this monoid space that we've mapped into. And then later, our transform output function is going to apply the initial element to uh uh to this matrix to produce the actual result in the space that we care about. And so this allows us by by splitting this this operation into a associated part and a non-associated part, and then doing this encoding and decoding, we can paralyze this operation. Details of it, not really super important. I don't know if it generates the graphs, but anyways, it works, it's very cool. Make some little nice little diagrams here. The particular like like algorithm that that Bradley and I were talking about here was uh one it's called uh like proportional error control. This is what you'd use for like um like a if you had like an industrial machine that you wanted to like keep on a certain target or like a self-driving car, like you want the car to follow the road, you would use a formula like this. So like imagine that you have a data set where your data set is like time series historical data of like you know a car's position like of like road conditions, and you want to see how well your lane keeping algorithm would work on this historical data. Uh well you could use something like this to in parallel simulate how well your algorithm would perform. You can also do things like uh uh differential equations, stuff like this using this technique. Uh going down to the next part of my long conversation with the LLM here. Um uh the interesting thing, I think, is this is where it's talking about this encode to a monoid space, do stuff in the monoid space decode pattern. But the interesting thing is the classes of operations that uh it says that it came up with as being a good uh uh uh fit for this because they have a small state. So it says things like what we just gave, which is like a uh he calls it an exponential moving average, or like a moving average with some some factor to it. Because the problem that I just described, it turns out it's actually just like some fancy form of a moving average. The other one is doing things like polynomial evaluation, uh, which is a tri a trick that we actually have some thrust examples showing this. Uh this is a more common like like scan and reduce uh uh thing. Uh things like integer carry, like propagate uh ads, um, things like a running minimum with a dropout, uh, which comes up in like machine learning uh like uh layers, and things like a max plus convolution uh all can be implemented in terms of this. Um and there's like a couple other classes of algorithms uh or uh classes of oper of non-associative operators that you can paralyze using this technique of uh encoding each input item into like a little state machine or a little state transformer and then composing those state transformers together and then scanning through uh uh uh and then and then decoding uh afterwards. Like the the idea here is that like instead of like applying your operation, like you you just store in your state that you're scanning like two or three terms, and then you uh the composition of those terms is associative and so you build those up and then at the end when you do your transform output iterator you apply you actually like piece everything together using those pop parts that you've accumulated. Oh I'm very excited about this kind of I'm sure that not not much of that.

Conor

No, I have to say that this is just very upsetting from an It's very upsetting from a non-academic point of view because a lot of what you just explained I was going to ask Jared if he was following because I am not. But then when you got to the end and you mentioned Horner polynomial evaluation it's just like um it's just academic speak for something that's very simple. If people think about a polynomial x squared plus five x minus six or whatever you can think about evaluating that polynomial as a scan because each term in that polynomial has an x with a coefficient and you can evaluate that by multiplying the value of x by scanning that kind of sequence where the coefficient is the number in an array. Which whether or not you understood my explanation just now is so much simpler than saying oh Horner slash polynomial evaluation scam blah blah blah blah blah. Anyways the point being is that like a lot of this stuff is a lot simpler than it seems I got lost in Bryce's explanation. I'm not sure maybe Jared I mean he is the author co-author of Thrust. I mean Nathan's still a Google I'm not sure if Nathan would have hung on but I imagine you know Jared was sitting there being like this is basic stuff why is Bryce so excited? I was sitting there being like uh this seems way too off in the weeds and then it wasn't until that end table that ChatGPT showed me where I was like I did I still didn't understand what you were saying but that stuff in that table I do know which is it's not that like Horner polynomial evaluation is three words that are you know a dollar to explain something that's quite simple. That's all I have to say is that uh I'm confused, but I probably do understand it in some sense.

Bryce

So I will say Connor the reason I I sh the reason I went through all that is because I too was quite confused by all the theory. But then once I got to that point of seeing the other of seeing that table and seeing the other examples I started thinking about a lot of the problems that we've talked about here in particular a lot of the problems involving scans and in particular a lot of the problems where we have built up some sort of custom state for our scanner reduction operator. And then the more like I went back and looked at some of them like the chunk buy example that I use in my thinkparel talk and then like some of the interesting reductions that we talked about on the uh the lost reduction episode and I realized that all of them fit into this like high level you know mathy formulation. And that I think there's some uh there's some like general lesson to be taken from here. And I also think that if you had like a high level like a ray programming language um that where you could do some anal analysis on your operator and determine some of the properties of your operator that you could automate this transformation so that you could um uh you could paralyze a non-associative operator. And then that got me thinking more about whether scans actually really require that your operation be associative. And I think I came to the conclusion that that's probably just a convenienced requirement that we've put that I don't actually think that a scan requires this. A scan probably requires that you have some associative component of your operation or that your operation can be composed decomposed into an associative and non-associative part where the associative part is non-trivial. And as long as you can do that where there's some non-trivial work that you can do in the associative component, then you can paralyze your operation. And so I think that there's actually like I keep encountering people who believe that they have a problem that cannot be paralyzed with a scan and a reduction. And I think that a larger class of those problems can be paralyzed than people think.

Jared Hoberock

What say you Jared It's always difficult uh at first glance to to find the parallelism in a problem I will say that uh and it's um it usually comes as a surprise when you figure out the the solution to the puzzle.

Bryce

Yeah that that is how I often uh feel I feel like I keep telling myself that like we we keep we've talked regularly on on this podcast about how like we so often reach for scans or reductions or transforms. Like those are like I think in my mind the three uh uh critical the three basis operation is scan, reduce and transform. Um and like just last week I was uh um I was talking with someone about uh how we could implement a rotate in this this uh you know this new uh tile programming thing that that we talked about earlier um and and somebody's like oh that that that can that rotate could just in theory be a scan on like a matrix and I was like oh my gosh you're right how did I not realize that why didn't why didn't the first my first reaction should have been oh yes I should try to fit this parallel programming problem into one of these three patterns that I know they almost always fall into but I think it's hard it's like it's hard for us to to see those connections and then it's like I always have these aha moments.

Conor

Well I mean uh I mean uh if we haven't mentioned it on this podcast, I mean a a rotate in parallel which is currently not implemented in Thrust is just a permutation iterator plus a gather. Yeah yeah that's that's the other way yeah a gather is just a scan which is what you just said so a scan on a matrix. I mean Jared is his eyes are pointing up at the ceiling.

Jared Hoberock

I'm not sure it's just that uh because you're mutating as you're um I think you need to synchronize or carefully uh I mean you're remind me how rotate works. You want to uh you want to shift move part of array to the beginning.

Conor

Yeah basically it's it's shifting the first n elements to the back.

Bryce

I mean there's technically a rotate left a rotate right but take n elements and shift those to the back and then the other doing it doing it out with doing it without like a scratch scratch buffer uh would take some thought I think if you're doing I think that Jared's right that if you that you can't do it in place with the scan of the permutation iterator you would have to write to a new uh a new location but I actually think that that's the in the particular case I'm I was thinking about it was perfectly fine because you can't do in-place updates anyways.

Conor

In this in this uh this tile uh programming model everything is a value there's there's no in-place updates of anything I was gonna say is uh uh pending that you are doing a uh a gather to a output sequence it definitely works right like if you're doing a permutation iterator from an input sequence and then a gather operation to an output sequence you don't have to think hard about it correct?

Jared Hoberock

Yeah that's right but I think the interesting case is when you want to rotate in place which definitely requires some synchronization.

Bryce

We don't implement that in um in thrust in cover we don't even implement a rotate out of place in thrust so the the the use case for um for a rotate that I was thinking of is imagine that you've got a tensor um and you want to do you want to like do like adjacent difference um uh on this thing uh but you you have to work with a tensor of that particular size. You can't make smaller tensors. And so you can't like you you have a like an n by n thing. You can't make an n minus one by n minus one things. You can only make n by n things. And so how do I do adjacent difference? Well in that world the way you could do adjacent difference is you have one n by n thing and then you create another n by n thing which is just the first one rotated. And then you can take you can subtract the two together and then you can map you can you ignore the last element which is you know invalid but then the first zero to n minus one elements in each dimension have the the the correct uh uh a result.

Jared Hoberock

You see what I'm saying what would adjacent difference on a tensor mean?

Bryce

Uh like a stencil operation.

Conor

Yeah it just means look at look at adjacent elements in the same row and take the difference. But I was going to say is the library that is yet to be released that I've been working on I mean um I literally ran into this in the in the thrust examples is a simple moving average where the way that they do that is they take the they they do a plus scan on your 1D sequence and then they take the difference between the ith element and the ith plus four element. But you need to do a exclusive scan where you prefix it with a zero. I don't technically have a pre-scan or uh exclusive scan implemented in my library but I do have prepend. So I just prepend the zero and then use that to do the the differences. So what you're saying is you're saying do a rotate one a one rotate and then take the differences but if you had a a prepend zero you could do the same kind of thing where you take the difference. So there's a couple different ways to do it. I mean it's all coming back to what Jared said is that uh in this day and age doing this on 1D sequences is just the wrong way to do it.

Bryce

You should be doing this from a multi-dimensional uh perspective um because so many of these things yeah don't come out ergonomically another interesting thing that I I insight or conclusion I came to from looking at this uh uh proportional error control uh uh algorithm is um this is the second time in the last couple months where I've had um like a a sort of scan by key type algorithm where I had a a complex aggregate state type that I was working with and where when the key when I move to the next key um I have to reset like one of the things in the state. Like okay when the key changes like i in chunk by when the key changes you reset the count of how many things are in the chunk. In this particular algorithm when the key changes like you reset the initial value to uh the the current value. And uh whenever these whenever I've run into these uh instead of using scan by key I just end up using scan and just putting you know the key into this aggregate state. And I've also just found myself instead of using the by key algorithms I more and more I'm just like writing the by key part myself because I almost find it like a little bit uh cleaner uh to just expose like all of the state that I care about uh which is like okay like when am I at the end of a group? Uh uh and then like what am I going to do once I get the end of the group? Um and I just I don't know I th I thought that that that it was interesting that when I get to these more complex cases I have to switch from scan by key to just uh scan and do the do the the the grouping myself well I mean the question for Jared is where did the motivation for the reduce by key and scan by key algorithms come from?

Conor

Was it specifically for these kind of matrix operations where you want to do row wise reductions or is it do you recall where it came from?

Jared Hoberock

I think it I think uh some people were publishing results on segmented scan uh at the time uh that people were interested in and it was clear it was clear it could be useful to have a primitive like that to you know do interesting things and the question was how do you uh how do you how I guess the interesting question is what should the interface be because you could do it like thrust does it or you could do it um cub does it differently they use um like head flags um instead of encoding per element which group each element is in you uh encode a bit where each run begins um or each group begins in in in the array so there's you know trade-offs um depending on how you choose to to phrase the problem but um I think at the time we already had sort by key so there was nice um synergy with a like a scan by key to go along with that you could write nice example codes basically if if it made for a nice example program it was probably a good idea so it went in the library.

Bryce

When did the first uh thrust sort algorithm come around? Because I think uh like the the sort algorithm is sort of what thrust is known for and it's one of the most widely like used.

Jared Hoberock

Like there's some projects that uh like they write their own like reductions but they the only place that they use thrust for is for the sort right uh Michael Garland had a uh summer intern who uh published a merge sort and I remember I inherited that code and tried to make it work in uh like for real um in uh in thrust I don't remember if merge sort preceded Radix sort or not um I don't remember my guess is that so there was a a similar library called uh CUDPP um that Mark Harris worked on with some people at Davis and they published some scan results and I think they also had a RATEX sort result but I don't remember if we maybe adapted the RATEX sort code before the merge sort code or what but um at the time lots of people were publishing how to do X on the GPU so we would just um uh steal from those results and um you know put some polish on them and and publish them in the library.

Conor

Be sure to check these show notes either in your podcast app or at adspthepodcast dot com for links to anything we mentioned in today's episode as well as a link to a GitHub discussion where you can leave thoughts, comments and questions.

Bryce

Thanks for listening we hope you enjoyed and have a great day low quality high quality that is the tagline of our podcast.

Conor

It's not the tagline our tagline is chaos with sprinkles of information