3 Ways to Speed up Your Regression Modeling in Python

Speaker 1 [00:08]

This is somewhat my personal story as well, how I learned to stay same while fitting very big regression models, and basically I can't stay same because I learned about a couple of techniques to fit these models in a reasonable amount of time, and I'll walk you through them basically. But, yeah, we kind of start at the very beginning of data science today, here the guy who invented everything Gauss linear regression, that's how far back we go, so this is not cutting edge technology, but really things from the 19th century. But for me, this is also basically the beginning. After leaving university, I started working at Trivago. And then we work with data sets that kind of look like this. Here you have a result list of multiple items. And then you scroll down, and there's like a lot of items in a given list. And each user sees a different list. And at the end of the day, we have a data set per user with 35 hotels or even more. We have multiple million users. And then if we want to do list-wise evaluations, you can imagine that data sets very quickly get very big. So if we have one million users and each of them see 35 hotels, we have already a data set of 35 million hotel impressions. And then you could think of having more users, them seeing more items. Data sets get super big. All right. And so this is basically kind of the use case. And today, we kind of have this walking example of me starting at Trivago, and I want to analyze interactions of users with the list. We have this full data set of user item impressions, and we know at which position of the list the item has been displayed. So basically, it's four variables and potentially extremely big data set, because for each user, we have multiple items. And then we just want to know, okay, is the user going to click on a given item or not? And in particular, we want to know, okay, will a user click on an item if we put a red badge on it? Yeah? So by how much has it changed? So typical e-commerce example. So we have this data set, kind of looks like this. We have 2 million users, a kind of universe of only 50 items, which is, of course, simplifying assumptions. And each user sees 10 items and lists it in a given position. And the data set then has four variables, item ID. We have 50 of them. If there is a batch assigned on a given item, that could have been randomized, the item position, and if the item is clicked. Okay. And that's basically what we have. Item ID, category of 50, batch 2, position. We have ten positions, and clicked is a dummy variable 1 or 0. And so then we can fit a linear regression where on the left-hand side we have the dependent variable, X is just a numeric design matrix, D is a special set of features, which is categorical variables that you would have to one-hot encode to actually fit the model, and you would be an error term. Now, then you get the ordinary least-squares estimator, which basically is the standard x prime x, and that you have to invert it, and then you matrix multiply it with x prime y. So basically, very simple formula, and so more or less what all the basic statistical software routines implement, sort of. And the main problem here is that this matrix X here, or X prime is of dimension N times K, so you have this very big matrix multiplication K times M times N times K matrix. So the end product is K times K. So if you have 100 features, you end up with a matrix that is 100 times 100. So you have to do the multiplication. You have to invert it. That is very computational heavy. So that's a little bit annoying. And if your data set is just very big, your statistical routine will take a lot of time doing it. Okay. So we'll fit a regression model today that kind of looks like this, click outs on the left-hand side, then we have a dummy that's a badge, and then we have two fixed effects, categorical variables, which is basically a position effect and an item ID, and we get this data set by the formulaic library, which basically handles all the data pre-processing for us. we get a pretty big design matrix with 2 million users and yeah each of them see 10 items we have 20 million observations and 60 features after the one hot encoding so you could say okay this is reasonably big 60 features of course we could be much more extreme right it could be 1,000 features if we really had all items that we have I think in Trivago we have millions of items if all of them are observed this this gets crazy crazy big but yeah to keep everything under control here we have 60 features and then we do basically we fit an OLS model textbook OLS we just implement the formula I just showed you by hand and so what we have to compute is two matrix products X prime X and X prime Y and then we can invert X prime X and multiply it with X prime Y and then we have our OLS estimator and we'll just do a benchmark how long does it take us on our data set of interest there's a lot of code here basically but we write around 10 seconds to get a result so 10 seconds I'm not sure which area of machine learning or data science you work in but in my area it's a lot of time because a lot of work is interactive and waiting 10 seconds and one model fit a little bit annoying okay so can we do faster so this is the challenge today can we make this run faster and on even bigger data sets can we actually make it feasible to fit these kind of regression models on our local machines without having to go to the cloud or to a distributed system. Okay, so these are the two challenges. Now, what I just showed you is only sort of what is implemented in most statistics packages. For example, if you look at the R documentation, you will learn that the lm function calls another function called lm.fit, which calls another function which is called .lm.fit, which under the hood implements something called a QR decomposition. And so this is basically a procedure to solve the same least squares problems that we have, but it's more numerically stable. So in SciPy, you could access this, for example, through NumPy and SciPy. There's this routine SciPy linear least squares. And if you choose this driver, you will call a QR decomposition. So this will be numerically more stable, which basically but has a disadvantage if you actually benchmark it, whoo, it will take forever. So this is by far the slowest kind of benchmark I'll show you today. It's like if you use it, you have higher security, but you wait and wait and wait and wait. All right. So maybe this doesn't generalize to every problem, but to the one at hand here. Okay. So but coefficients match. So there's a sanity check you'll see throughout. So this is good. We get the same result. Okay. So what else could we try? So this is basically what I asked myself, right? I tried to run a regression on a big data set and I waited and waited and waited. And so what else is out there? And so basically I think there's a few strategies. First is to run batch regression. So we basically split our data set into chunks and kind of run the regression on chunks of the data set and at the end combine the subresults into one big result. We can then also make use of the structure of the data, which we kind of learned about earlier in the GLAMP presentation. So for example, if we have a lot of sparsity in our data, we can use sparse solvers. If we have these fixed effects, we can use an econometric method called the Frisch-Vogt Lobel theorem, or we can then even try something else that's really cool. We can compress the data set into kind of a smaller data set and then fit on the compressed data set using weighted least squares. Okay, so these are the things I'll show you. You can also, of course, try to put everything on the GPU and see if that's faster. But I couldn't do this on my local machine. You might see it here. It's, I think, a laptop from 2018. I didn't think of GPUs back then. Okay, so first, data sparsity. What does it mean? Basically, in our design matrix X, we have a lot of zeros. So here in this example, we have 94% of the data just consists of zeros. And basically, there's optimized data structures that they can use, that you can apply there, and that might help you quite a bit. So just to drive the point really back home, here we have more or less a data frame or a matrix. And there's a lot of zeros in there. So we would call it sparse. All right. So there's a couple of sparse solvers implemented in SciPy sparse. There's also sparse matrix algebra operations implemented there. So why don't we just use these? And so this is basically what happens here. We just reimplement the kind of direct or less fit on NumPy arrays on a SciPy sparse matrix, which we convert outside of here. So it exchanges a little bit, but that's it. We kind of sparsify this operation here, which is the heavy operation, really, that we need to get a handle on. All right. And so how does it look if we then time this? OK, we already do quite a bit better. So this is around a 30% improvement by just making use of the fact that, okay, we have a special data structure, so we can choose a more appropriate algorithm to actually tackle this problem. Okay. Do coefficients match? Yes. They match again. Cool. Okay. Very nice. But that was just easy exercises, right? So you would assume that we get equivalence here. Then in SciPy you find one kind of approximate sparse solver here in LSQR. You can call this this one as well, just to show you the benchmark. It is again taking more time than even our baseline estimation. So maybe that is not the fastest solution for us out there. Coefficients match almost and that is again the reason is that it is just an approximate solver and not an exact solver being employed here. Where these sparse solvers are implemented, for example, in scikit-learn, right? You could call scikit-learn regression, linear regression model based on the sparse matrix, and all the operations would run on sparse matrices, or there's another Python library called fastrack that purely runs on sparse matrices and gets pretty good performance out of that. Okay. That was sparsity. Now what if we don't have sparse data structures? What else can we do? We could run batch regressions. What's the idea here? Basically, the outline is we just chunk our data set into different pieces and then fit sub-regressions on each piece of the data, and we bring it all together at the end. And sometimes strategies like this are called split, apply, combine. So we split the data, we combine the results, we apply a function, we combine the results, we're done. Okay. Now imagine we have n observations and we just chunk it in g groups, say 100 groups, and we have basically group 1 to 100, and this is kind of how it looks like. So we just stack the groups on top of each other from group 1. So here you have multiple observations, and group g you also have multiple observations. And what we can then do is rewrite this whole s equation into basically a summation operation of the X prime X operation by group, then we sum it up and do the inverse on the sum of this problem, and we can do the same with the X prime Y operation. And why is this great? Because all of a sudden we don't have to operate on the full N times K matrix, but just a subset of the N observations. So that helps us. Each individual matrix operation gets a little bit smaller, but we have to do a couple of them. So this is basically what we do here. In principle, we can parallelize this, we can save results to disks if we wanted to, we can read subsets of the data from disks and write to disks in each iteration, and so on. So how would this look in Python code? Basically it's quite simple. We have this batch function here, and we have a for loop, and basically in each iteration of the for loop, we create a subset of the input x tilde, call it xg, the subset of y, YG. Again, compute X prime X, accumulate it in each step, and do the same with X prime Y. And then we have this full object X prime X and X prime Y, which are identical to the basic objects that we computed earlier. We can pass it to our regression function, and that's it. Okay. So that's batch regression. And that works reasonably fast for the example that we have at hand here yeah pretty cool of course we have to check do we get the same result yes we get the same result here this is slightly off for the approximate reason that's the same as a sparse over same as least squares so this is a viable strategy to kind of reduce the dimension of your regression model make it feasible to fit something that was very big on smaller chunks all right we can And of course also run this on the sparse matrix, right, performs kind of equally well, so not much improvement. Okay. So where would this be implemented in the Python ecosystem? If you think of distributed frameworks like Spark or in Dask, they would basically run strategies like the one below. Basically what happens there, they chunk the original data set into different pieces, send them to different runners, evaluate kind of these individual pieces, and bring it all together at the end. So this is what happens very abstractly in SPARC and I think in Dask as well. I cannot fully promise you that because I cannot really read scalar code, but abstractly this is how distributed regression would work. Then another thing you can do is apply an econometrics theorem from the 1930s. This is called the Frisch-Vorg-Lovell theorem and if you study econometrics you hear it and if you study statistics in any other science you won't hear about it but it's actually pretty neat it also makes use of the structure of your data and again here recall that we have this model of interest where we have regular numerical features and these so-called fixed effects which are just categoricals and one thing that is annoying about one hot encoding categorical features is if the routine that does it is not very efficient it it might take some time and get this very big matrix. It takes a lot of memories. So the question is, what if we actually didn't have to do this, right? So if we didn't have to do the one health encoding and just skip this step. And so basically the Frisch-Vogel-Lovitz theorem allows us to do it at the cost of actually not getting parameter estimates of this alpha parameter here. And so the key idea is that we residualize the matrix X and the matrix Y by the fixed effects. And there's a little bit of mathematical or linear algebra theory behind it that I won't go into detail. But basically there's, I think, one page proof that you can look up at Wikipedia that kind of will argue that what I'll show you works. Okay, but just to show you some code, we need one function that is actually doing this residualization or demeaning for us. And we can load it from the PyFixes library, and it's kind of in here, pyfixes.estimation.demean. And then, basically, we write a function that we provide with the design matrix X of numeric features, the variable Y, and just the fixed effect F, which, as we see below here, is just basically a tabular part of a data frame converted to a NumPy array with the position and the item ID. So, we don't do any one-hot encoding. this is an array of dimension n times two and Then everything that we need to do basically is pass this to the demeaning function which PyFixas provides for us which does the demeaning and we get kind of this matrix YX demeaned out and Then we can take the dependent variable the demean dependent variable the demeaned covariate X as well And then just do an or less fit on that Okay. And so that's very nice, because in the example at hand, our fixed effects were basically 49, dimension 59. We really just interested in one target variable. So this thing here, x, is of dimension x times y. We save us a lot of work. And that's cool. Okay. And let's look at the benchmark, how it works. Okay. And we're doing even better than the sparse matrix. We're at roughly 2.7 seconds of runtime for this problem. All right. So, are the coefficients the same? Yes. Again, we get exactly the same result. And so, that's always the coefficient, by the way, I didn't say this, on the kind of batch parameter. Yeah. So, the one I'm really interested in, does it drive users to click on an item? Okay. In Python, this is implemented in two packages. One is called linear models that has absorbing least squares, class, and PyFixers, which mimics the rfxest package. Okay. Then, last talking point for today, and the last strategy, and actually something that's very cool. So compression and sufficient statistics. So the key idea is that we make use of the structure of our data set one more time. And so the question basically is, is there redundant information on our data set that we basically don't have to kind of keep in memory at every point in time. And most importantly, we ask, okay, are there copies of specific rows in our data set that are duplicates of each other? So basically, we show one item to one user in one position, and then the behavior is the same. Maybe we have the same exactly row in our data set multiple times. Then we don't actually have to keep that in our data set, but we can compress our data set, and then use weighted least squares to fit a regression without any loss on the compressed data set. All right. And so, basically, the strategy can be summarized as if you have a big regression, at the end of the day, maybe it's a small regression with particular weights. All right. A personal background for me, by the way, is we used to have a Hadoop cluster, and I tried to get very big data sets out of it and I couldn't yeah didn't work I didn't figure it out so I had to kind of shrink the data set to actually be able to work with it so that's when I started to look into this and so that's a range of strategies the first is called compress and weighted least squares and so key inside again is we have four variables in our data set and what if we just compress our data set by counting how many unique strata of these variables we have and the strata is basically a combination of unique values of these different variables. So here in pandas you basically get the idea we do a group by over the list of variables and then basically do a count how often do these occur and this counts is then basically a frequency. Okay and if we do this in our data set of 20 million observations basically what we see is that the first stratum here with item ID 46 shown in position 7 without a batch not being clicked we observed this around 120,000 times yeah and so in the compressed data set we have one column for it in the uncompressed data set we had 120,000 columns for it okay so this is basically already quite a heavy reduction in how much data we have to work with and this basically at the end of the day, compressed data set that we get is not 20 million observations but 2,000 observations. So massive improvement, very easy to handle an object like this in memory and to actually also do numerical algebra on it. Okay, so yeah, recall we started out with 20 million, now we have around 2K observations. And so then we have to use a slightly different estimator, which is the weighted least squares estimator, which is more or less a generalization of the ordinary least squares estimator. You now have this weight matrix in between here, which is basically a diagonal matrix with off-diagonal elements zero. And on the diagonal, you just have the weights that we've computed, and in our case, it would be the frequency of different columns. Okay. And so here in three lines of Python, this is basically how you would implement it. So super simple. You can just multiply Y and X with the square root of the weights and pass it to an O less routine. And that's it basically. All right. And so just to kind of really drive the point home for one last time, here's from a very nice paper that I encourage you to take a look at if this piques your interest. Basically we have this observational pair A1. Basically we have it twice. And so we have kind of an aggregated data set where we have this count variable that It says, aha, combination M equals A, Y equals 1. We have it two times. All other combinations, we only have it once. So this is basically the compression step. Okay, and then we can fit with at least squares on it. And here's an implementation of this in Pandas, where basically we do just the same thing as before. We have a data set. We group by basically all the variables in our data set, create a count. The count will then basically be equal to the weights that we use in the weighted least squares problem, and then we fit weighted least squares and that's About it basically so implemented in pandas oops unfortunately. This is very slow. Yeah, okay? Admittedly I do some things that are yeah, maybe not Super optimized right there's a concatenation step in my benchmark, but still Okay, but main takeaway is here have the compression basically we get the same Point estimates as well and we can also recover standard errors with a minor correction So this is basically lossless procedure. You get point estimates and standard errors. I love you Gresham model. Okay, cool So when does this work? Basically, if your variables are non-continuous if you have a lot of discrete variables with not super high cardinality And so I think this actually in a lot of e-commerce setups I think this actually holds and you can also discretize if you have continuous features you can discretize, you lose a little bit of precision, but you get quite far doing so. Generally, the lower the cardinality of your values, the better. Because the more you can compress, the more likely you will have duplicates. Also the fewer variables, the better. And yeah, that's it. That's one challenge. Maybe the outcome variable is continuous. And then the question is, okay, is there something else that we can do? And yes, there is something else that we can do. It gets even better. slight adjustment to the strategy, and so this is again from the Wong et al. paper. So instead of just doing a group by and count, so counting frequency, by each stratum we now basically report the sum of a variable and the sum of squares and count the weights. And basically here is the strategy in a nutshell. You compute the sum and maybe also the sum of squares of the dependent variable for basically all strata. Now you take out the sum of square as like the dependent variable out of the group by. And you compute the count for each stratum. And then you can basically fit a regression of basically the variables that you have on the mean of the dependent variable. And then you use weighted least squares. And then you're good to go. And again, this is a lossless procedure. If you run A-B tests and analyze them in SQL sometimes, what you will do or will see is basically the exact same strategy. You group by the treatment, compute the sum of your dependent variable, you count how many observations were in the test and control group, and you do the same thing to get the variance estimate for this. You need the sum of squares of your dependent variable. So here this is basically just a generalization of the same procedure for a regression application. You just have some other covariates in between. Okay. So this is basically generalization of A-B testing in SQL. Okay. And here you kind of see how this works in Pandas. Basically we have a group by about aggregation variables, which are now only the covariates, and then we compute the sum of the dependent variable and the sum of squares of the dependent variable as well as the frequency. And then we have an even smaller data set than before, we have around 1,000 observations left. Okay. And based on that, we can then basically losslessly fit a regression model with weighted least squares. And losslessly means if we have the sum and the sum of squares, we can recover the point estimates and the standard errors at the same time. Okay. And again, implemented in Pandas, also not helping so much. All right. But, again, we get the same result as before. You see it here. So, valid strategy to go forward. Surely, this would be faster in Polars. Of course, this is the Polars benchmark, still not heavily optimized, so you can improve a little bit. And then some of you might ask, okay, if we can do it in Polars, can we do it in DuckDB? Can we potentially even do it out of memory in DuckDB? And so what happens if we move all the Polar's code into DuckDB? And so there's a library for that called duckreg. And so let's see what pops out of that. First, we have to do a lot of work to send all the data to the DuckDB database basically created. So this is what we do on this slide. And then we call a duck regression. Basically it has a formula interface. We fit the same model as before. We fit it just for one go. We don't compute standard errors or nothing, just a pure fit. And let's look at the timings. Okay. So now we are down at four seconds. And so that's pretty cool. Yeah. And so what does this mean for the viability of the strategy? Second principle, if you want to fit a very big regression model and your data, so it's the strategy by having relatively few features that have not crazy high dimensionality, you You can do this compression step out of memory, for example, in a library like DuckDB. Or if you work in the cloud, you could basically do this in the cloud SQL dialect that you're using. Just compress there, compute the sum of squares and the sum of the dependent variable, load the compressed library into Pandas. It might become very slow, and then there you can fit a regression in, yeah, regression package of your desire. a small tweak that you have to do to the standard errors. So, if you just plug the compressed data into stats model, standard errors will not be exactly correct. Darkreg would handle this for you. Okay. So, really cool. You can basically fit all these models even out of memory if you wanted to. Okay. I did a couple of more systematic benchmarks. My main takeaways I'm summarizing on the next slide here. If you can load your data into memory, I think Sparse regression is a super viable strategy, just call it through scikit-learn, and in particular relevant if you want to do prediction afterwards. If you're more interested about inference, then applying the Frisch-Vaugh global theory might be more relevant for you. Here it's relatively easy to compute standard errors for sparse matrices, not necessarily basically earlier in the talk on GLUM, the so-called sandwich matrices, where I mentioned And those you would have to compute here using the sparse regression. And basically these matrices would not necessarily be sparse anymore and would then have needed to be inverted, which takes a little bit of time. Okay. If you don't, yeah, if you cannot fit your data into memory, so there's two things that you can do. Either chunk it up and use a batch regression strategy, you can either roll your own estimator or use a distributed framework like Spark, for example, or you can use compression strategies if the structure of your data actually merits it. Okay, and so, of course, you can also combine some of these strategies in a fruitful way. So, for example, one thing that I do quite a bit is first use a compression strategy to get a frequency weight, And then I use Frisch-Vorg-level on the compressed data with weighted least squares and frequency weights. Okay. Well, these things implemented. Just to run you through this quickly. Standard normal equation. You get it in SPRESS models, PyFixS. You would also get it in Scikit-learn. SPARSE estimators you would get in Scikit-learn and in FastDrag. Batch regression typically implemented in SPARC and DASC. I didn't find an implementation in scikit-learn, so if this exists, please tell me. Sometimes these methods are also called online learning, online regression, but I didn't see it. Fridge-forth level, PyFixed and linear models, and the weighted least squares and sufficient statistic strategy would be implemented in PyFixed in memory through POLARS and out of memory through DAGREC. Okay. And yeah, I promised to talk about GPUs. I didn't really talk about GPUs involved in the PyFixers project. I'll do some advertisement at the end. We have played around a little bit with GPUs. They are basically implementing the FrischWalk-level demeaning algorithm in JAX. And we actually see pretty good results for very big problems. So if you feel like FrischWalk-level is appropriate for you and you have access to a GPU, maybe Maybe you can slice some performance out of using Drax as a back-end for PyFixerst. Okay. And so, yeah. More shameless self-promotion. Yeah. So, we're actively developing both PyFixerst and the dark rig libraries. If you're interested in these, please say hi. In particular, I've tried to rewrite the PyFixerst back-end, which is currently in Numba, to Rust. So, if you're a Rust enthusiast and want to help out, please, yeah, please come and talk to me. Okay. So, you can find both PyFixers and DAGRAC on Github in the PyEconometrics project, and here's a couple of links. All right. Thanks so much. Yeah. Thank you, Alexander. Thank you. So, we don't have much questions for this session, only one question. building the graph, you turn the unstructured days into structured. After this step, can the LLM not only identify the objects in graph? This is probably another question for me, right? This has nothing to do with, you know. So with that, there are no more questions, so I would like to thank you again for your presentation, and we are closing the sessions.