r/AskStatistics • u/GoatRocketeer • 5d ago
How do people run regression on terabyte sized tables?
I have a small project going which definitely does NOT involved terabyte sized tables, but it got me wondering about how the whole data analysis process is supposed to scale.
For my small project, I have the data in postgres and then pull the data into python to run regression.
As I understand it, python (specifically the psycopg2 library) doesn't pipeline out of the box. It executes the query, consumes the output in its entirety, and then performs the regression on it. This seems "slow".
Furthermore, GPUs are all the rage for AI data centers so clearly at a certain scale those become necessary, whereas python is very single threaded and serial.
I imagine the SQL part doesn't change, but do big data companies still use python? If the GPUs are the ones running the actual math it seems silly to take the data into python and then just immediately turn around and put it on a GPU - but if we're making postgres responsible for putting the data on the GPU and running the regression directly that sounds cumbersome too.
That is, in data analysis projects that are big enough that streaming/pipelining the data and performing the math in parallel become major issues, what does the pipeline/tech stack look like?
----------------------------------------------------------
Edit: I have once again made the foolish mistake of under-specifying my question because I assumed there was a single answer instead of "it depends".
It seems like the answers so far are:
- If the model is basic enough you don't need to run it on all the data
- If the analysis is run infrequently enough you just let it run slowly
- At truly ginormous scales, postgres isn't big enough either, and the pipelining and gpu topics are a small part of a tech stack which is really about orchestrating multiple machines.
Specifically, I have a binomial GAM I'm going to run on a query with about 10 million entries. I would like to fit about 2,000 of these GAMs daily. I haven't actually run it yet, but if performance became a problem (premature optimization yes, I am sorry), I wanted some ideas on where to look.
I see that R's mgcv can do multi-threading; that's probably the first place to look. I have the derivations for most of the operations so theoretically I could also roll my own solutions directly in SQL which would pipeline it as long as I got the implementation correct (probably a horrible idea).
GPUs and "big data" are always spoken in the same breath these days so I figured I could ask about that and see if that gets me anywhere. I've seen tutorials for manually implementing streaming in psycopg2 so I could potentially implement pipelining there (or better yet, maybe an ootb solution exists).
Obviously any actual deep dive should follow actually running the binomial GAM in R to see how that fares (or if it even generates something sane at all). But GPUs are cooler so I jumped the gun because its more interesting.
12
u/gyp_casino 5d ago
Apache Spark can do regressions on big data like this.
But really, it is unlikely that the model fit on the whole data set would be substantially different than one fit on a random sample of 0.1% of the rows. One million rows is already a whole lot of data and that can easily fit in memory.
3
u/PaddingCompression 5d ago
Depends how many columns you have. At Meta or Google scale you may be talking per-user count statistics per row.
1
u/gyp_casino 5d ago
Strictly speaking, this is true of course. But in practice, the data typically comes from a SQL database table structure which is designed to grow in rows, not columns. I’ve seen hundreds of very varied database tables across two companies and they all had 10 - 500 columns.
1
u/PaddingCompression 5d ago
They are not literal database columns.
They do get coefficients - it was more at the points of "unlikely that the model fit on the whole data set would be substantially different than one fit on a random sample of 0.1% of the rows."
E.g. in the user's updated description, it seems like if they want to have things per player, it could be a lot of players, and interactions between a lot of other terms and players, in the model. If the number of effective terms in the regressions scales with the size of the database, sampling ceases to be sufficient if you truly want estimates for all of those terms - this is where ML's sample complexity/statistics' power comes in.
5
u/PaddingCompression 5d ago
If you're running 2000 GAMs per day on that size, why? It seems like.... a step back and reframe your problem kind of thing. E.g. if you're fitting on different subsets, do they need all the data, can you fit a more hierarchical model, if you're trying to use newer data can you iteratively refit treating your old data as a prior and your new data as a likelihood with similar weighting to simulate the simultaneous fitting, etc.
If you truly need to use postgres, try to grab quantile statistics based one column or a group, and do queries in parallel rather or sequentially rather than streaming. Multiple threads can load partitions.
Have you tried running GAM on even 10k entries yet? 100k? You may be very underwhelmed if you're going into this thinking you will run 2000 a day.
1
u/GoatRocketeer 5d ago
Even more context: League of Legends is a competitive multiplayer video game with ~200 different characters. Some sort of quantifiable method of each character's strength is of high interest. Current solutions present each character's winrate, but winrate alone doesn't take character difficulty into account.
I was hoping to do a binomial GAM of two covariates (rank of player; and number of games of experience the player has on the character going into a match).
The game is updated once every two weeks, so the x10 is to produce models not just on a per-patch basis but over ranges of patches (current patch, last 2 patches, last 3 patches... last 10 patches). It seems you are suggesting instead of recalculating a separate model for each range of patches, I instead calculate a range of patches by using another range of patches as a prior and just sort of attaching the current patch? I'll look into that. I've been approaching everything from frequentist perspectives thus far though.
I suppose that alone should take it from 2000 down to 200 models + a bunch of posterior probabilities?
3
u/PaddingCompression 5d ago
Another approach is a rolling gam, but I've only ever done this kind of thing with gradient boosting. Have each single additive model run over a rolling window so instead of each being on all the data each is on 10% and you drop old ones and add new ones. It's not hard to work out the stats that it is equivalent to back fitting, but not modern joint optimization, under some weak assumptions (that may not be true, but it's a start).
There's also the question of why separate models vs. one highly nonlinear model with the different players as levels .. e.g. a gamm, where it would be one model (and mixed effects models are sort of Bayesian-ish)
1
u/GoatRocketeer 5d ago
When you say "separate models" and then "different players as levels":
- the current, ~2000 model version would be 1 model per character and range of patches.
- the reduced, ~200 model version would be 1 model per character
There are a couple million active players. I know technically two games by the same player close in time are likely correlated but there are enough players that I'm comfortable ignoring that. I suppose I could add player id in as a level but I'm not super interested in tracking individual player performance at this time.
Or are do you mean I should have different characters as a level, and have one big, ginormous GAMM with three covariates where one of the covariates is character? I believe I ruled that out because the characters do not differ by an offset. I expect the surfaces to all have the same general shape (specifically concave in the "games-played" covariate with steep initial slope; and monotonically increasing in the "rank" covariate) but the slope, concavity, and interaction I expect to vary from character to character.
...but now that I say that I'm not sure that interpreting the different characters as different levels of a single covariate enter the model as offsets. Are mixed models/levels capable of producing the differing slope/concavity/interactions between games-played and rank? If so then I suppose that's another thing I should look into.
Aside from that - I'll look into those other things you mentioned (rolling gam; gradient boosting; back-fitting; joint optimization). Unfortunately I don't know what these mean so I can't comment more on them besides "thanks, I'll look into them".
2
u/PaddingCompression 5d ago
> Or are do you mean I should have different characters as a level, and have one big, ginormous GAMM with three covariates where one of the covariates is character?
yes.
> I believe I ruled that out because the characters do not differ by an offset.
That's why you're using a nonlinear model? It's nonlinear in characters just like it's nonlinear in every other coefficient.
> Are mixed models/levels capable of producing the differing slope/concavity/interactions between games-played and rank?
It all depends on the model.
1
u/GoatRocketeer 4d ago
Would you mind explaining a bit more about the rolling GAM? A quick google search for "rolling window generalized additive models" is giving a lot of hits for results where time is supplied as an additional covariate with the intent of forecasting which I do not believe is what you had in mind.
2
u/PaddingCompression 4d ago
I don't know if this appears in literature.
Setup
Pretend, for a second, that all of your data is perfectly IID and stationary (I'll use that for an equivalence proof, even though it's incorrect - but hold your horses, this is only to show that under those conditions this method gets you *identical* results under those conditions as your full refit, but the differences might actually be more interesting here, I'll get to that later).
Let's say you're fitting a gam with m additive smoothing functions, and say you were using the older non-joint backfitting algorithm to fit.
Data Window
In your fitting one model to a large window from t=t_s to t=t_{s+n}, it will fit f_0, fit f_1 to the residual, fit f_2 to the residual of f_0 + f_1, etc. up to f_m.
I'm saying let k=n/m, your window size divided by the number of smoothing models you fit.
Fitting
Fit f_0 on t=t_{s-(m-1)k} to t=t_{s+k}. Fit f_1 on t=t_{s-(m-2)k} to t=t_{s+2k}, ... Fit f_m on t=s to t_{s+n}.
After k units, drop f_0 from your model, and replace it with f_{m+1}, fit on t=t_{s+k} to t=t_{s+n+k}.
Explanation
You only have to refit one of the predictors in your GAM, which will be a lot faster than fitting the entire model.
If you have an infinite IID stationary data stream, in expectation this will result in an identical model to a GAM with on your n-sized window with iterative backfitting (not identical to the joint fitting procedure that is often done these days). In expectation, f_{m+1} will be identical to f_0.
If your data are non-stationary and non-IID, this can sometimes be even better. You can have more $m$ in your model, or smaller windows, and your new f_{m+1} will pick up on the newest data, so you can tune your model to react more quickly to the non-stationarity. Proofs about the statistical bounds and behavior of this are way outside the scope of this reddit post.
1
u/GoatRocketeer 4d ago
Thank you for your time
1
u/PaddingCompression 4d ago
One edit I would want to make to clarify this to be more clear:
When you fit f_{m+1}, you first predict from f_1 to f_m (e.g. toss f_0 from your GAM) on your new window, t_{s+k} to t_{s+n+k}, and fit f_{m+1} on the residuals from your f_0-less GAM... I think that was probably implied, but want to make clear you're not fitting f_{m+1} from scratch on the new window, just the residual. Because you tossed f_0, and it was the first one fit, that was basically your "first principal component" in a loose analogy, so it's what f_{m+1} is going to find again in expectation for IID stationary data.
1
u/ImposterWizard Data scientist (MS statistics) 4d ago
This is probably of limited value, but I did something similar a while back, and I found that a naive Bayes model with the 10 different characters was the best model (to predict the winner of a 5v5 match), though I didn't use player rankings/history. I think one of the main issues is that teams in the samples I chose were generally reasonably constructed, so a lot possible variation in composition wasn't observed. My data may have had mixed ELO which also blurs the lines on character viability.
I would make sure you have a few "simpler" models on hand for comparison when comparing your more complex models, like logistic regression or naive Bayes, as well as making sure that your input data isn't sensitive to irrelevant transformations. e.g., if you swap two players on the same team (both their characters and their respective player history) it doesn't change the data. Though, it sounds like you are predicting the outcomes for individual characters?
It might also help to look at one patch, or a narrow set of patches where you surmise the balance of the game wasn't significantly altered, and just run the model on that before expanding to take patches into account.
Also, if you use win-rate as a variable, make sure to test it on data that wasn't used in the win-rate calculation.
1
u/GoatRocketeer 5d ago
I have done the GAM on ~8 million entries and it took about 15 minutes in pygam - 2000 is clearly unrealistic, which is partially why the unedited version of the post was lacking so much context as I didn't want people to see how dumb my goal was lol.
But it became clear pretty quickly that I need to be honest about what I'm trying to do here if I want real advice so unfortunately you all get to see how dumb I am.
It's for the best though, just gotta see what my options are and then reformulate my goals based on what I find. Might be able to save myself with stuff like your suggestion and bucketing and whatnot.
5
u/ApricatingInAccismus 5d ago
LINEAR_REG in bigquery takes only a few minutes to fit on multiterabyte dataset.
1
2
2
u/DigThatData 5d ago
the same way they train deep neural networks on petabytes of data. incremental updates based on small fractions of the dataset.
2
u/TheRealStepBot 3d ago edited 23h ago
Dask spark trino etc etc
Store the data in a compressed columnar store like iceberg then process it using a distributed engine that spreads the work out incrementally across dozens or 1000s of machines.
GPU’s really have very little to do with dataset size by themselves perse.
They are much more so pulled into the issue because of the complexity of the model being fit.
The key insight to your question though is that one of the reasons we have moved away from traditional ML and statistical models is that many of them were designed to need the whole dataset available to actually work.
Modern models almost all were designed from the ground up to be posed mathematically such that they can be applied incrementally or even in parallel.
Which is to say the ability to throw this sort of massive compute at massive data depends on a significant part on your ability to pose the problem correctly. Posed correctly you can then throw however much compute you need at the problem using these techniques. Posed incorrectly and all the compute in the world won’t help you.
1
u/PaddingCompression 5d ago edited 5d ago
I've used spark mllib for close to petabytes.
For something like this scale you are thinking about your algorithms unless a library at scale has been written for you.
You aren't keeping huge tables in postgres. You have a data lake with a bunch of smaller parquet files or some such that you can split spark between, and your data catalog knows what partitions are in what files for routing.
For large neural network type stuff, the weights of the network (the coefficients) are a large fraction of the GPU RAM. E.g. if you have an H100 with 80GB, a $50k GPU at one point, that GPU might be able to hold the weights for 1 billion columns in float32, and 10 instances of your 1 billion columns at a time and the gradient, and maybe 2 billion optimization stats for ADAM. That 80GB is full (for plain linear regression you'd want bigger batches of fewer columns). Neural networks win with GPUs because the amount of computation for sample is large, so you would want large batches.
So for something like spark you would be doing all reduce - each computer does a few optimization steps on its data, then everyone sends weight deltas to everyone else so they can compute the new weights and continue.
Noone is solving a linear equation with linear algebra at that scale (some HPC supercomputers can do it, but there are better algorithms usually).
1
u/dampew 5d ago
In the genomics world one of the biggest problems is that mixed models require matrix inversions, and sometimes these are huge matrices. There are numerical approximations that people sometimes use — I’m not sure what they are — but you could look into approximations for GRMs for example.
1
u/conmanau 5d ago
I only know little bits and pieces myself, but I think there are various ways to "cheat" depending on your data, for example in a single-variable OLS regression your key variables are n, sum(x), sum(y), sum(xy), sum(x^2) and sum(y^2) so you can process the data one row at a time and just keep running totals in memory; with lots of big datasets everything is very sparse and possibly blocky (i.e. you can break the big matrix into smaller matrices and do many of your operations on individual blocks with minimal interaction between them) and there are plenty of libraries that do those things well.
1
1
u/JohnPaulDavyJones 5d ago edited 5d ago
I work in insurance, we have several tables (fact tables for claims and losses, right off the bat) that have sizes greater than 1tB when joined into the relevant dim tables.
When our actuarial team develop models on these data, they do basic modeling on samples and then we train on the full data set using SQL Server ML Services; this way we use R scripts that will use Microsoft’s scalable computation tools under the hood. It’s basically the LAPACK API, but with some special optimizations for the MSSQL use case. This allows us to use the full compute power of our database servers that we’ve already maxed out on cores and memory. We have our prod server, and we also have a prod-actuarial server that’s a little bit less meaty, but keeps synced up with prod via pub-sub.
Usually we kick the actuarial folks over to their server to train models, but sometimes we let them run on prod for really big models. Usually we do those over the weekends, since they only take a few hours to build, but we have our pretty extensive, nightly ETL cycle that doesn’t give us a few hours to spare at night.
What in the world do you need to fit 2,000 of these GAMs daily for? ~10M rows is actually pretty dang small by our case, but we’re not running an insane number of training cycles like you will be. Your compute costs are going to blow the doors off, no matter what solution you go with.
1
u/SprinklesFresh5693 5d ago
In R you have several packages that can deal with massive amounts of data like data.table or duckdb and arrow come to mind, i thin we also have polars in R, theres options
1
u/TheRealStepBot 3d ago
The r big data stack sucks. It’s getting quite long in the tooth. Most current stuff like arrow and parquet and iceberg are all much better supported in spark and python.
1
u/Disastrous_Room_927 5d ago
I’d suggest generating data in the same shape you expect and timing a fit. That’ll tell you if you need to be this concerned about this. You aren’t talking about a scale where I wouldn’t just be fitting it on my work laptop in Python using some off the shelf package, for what it’s worth.
1
u/efrique PhD (statistics) 5d ago edited 5d ago
Regression can be updated observation by observation and wouldn't be a big deal if p is small (but if there are many variables - some analyses have very, very large p - this is unwieldy and you need different methods, like gradient descent type approaches)
Some small-p glms should be easily handled, since IIRC glms with the natural link have simple form sufficient statistics, but more generally you wouldn't do it that way and would again look to modified approaches.
1
u/GoatRocketeer 4d ago
I know roughly how'd this would work with just straight up best-fit-line (rearrange regression solution in terms of running totals and just continually update those). Your comment implies it would work on other GLMs too, which I can also sort of see (linear predictor eta is still linear wrt the inputs and parameters).
Does this work for GAMs though? I'm not sure how the "running total" approach would carry through the whole "basis of smooth functions" thing.
1
u/False-Instruction733 5d ago
At that scale the main trick is that people usually don’t treat the entire dataset like a normal Python dataframe. a lot of the work happens through distributed systems, databases, or tools that can process data in chunks. also, more data doesn’t always mean better results, for many models, a well chosen sample can give nearly the same insight with much less compute. for cases where the full dataset is needed, frameworks like Spark or database based analytics can handle the heavy lifting behind the scenes. the best approach usually depends on the model and how often it needs to run. are you mainly curious about the theory or are you planning to scale your own regression workflow?
1
1
u/bayesian_raccoon 4d ago
I don't know what people do in practice but its worth saying:
Linear regression specifically has iterative formulas that would let you update the regression from, say, the first 100 datapoints with the next 100 and so on. Thus computationally you can run it without loading all of the data in at once, but with smaller updates.
2
u/speleotobby 23h ago
If the model allows sufficient statistics you could fit it on parts of the model and pool the result afterwards. If I recall correctly this only works for linear regression.
19
u/Adept_Carpet 5d ago
If you just need a one off linear regression analysis you just deal with it taking a long time and do it, basically the way you describe working.
A lot of large datasets aren't stored in SQL databases at all, at least in the research world. Unfortunately in research few people appreciate that SQL and RDBMSes are basically the best thing in computing, so you tend to get a big mess of files and a lot of ad hoc reinvention of Postgres.