I have thousands of subjects, for each of whom I have a fitted gamma distribution, with the parameters estimated from each subject's data. It is easy to look at the distribution for one subject (say qq-plots, etc.) to get an idea of how good is the fit. But how can I do this on a large scale, for all subjects?
I'll suggest using representative plots. Pull 16 or 20 subjects, and show their QQ-plots in 4x4 or 4x5 chart. Sometimes, you can plot several subjects in the same plot. This doesn't substitute other ways of representing the fits, but on the other hand I don't think you can avoid this step either. It's used a lot in panel (longitudinal) data analysis. You really need to see the representative plots.
See, Fig.12-1.3 in this book. It's not the distributions, but the same idea: show the sample plots for subjects.
You can get fancy and draw 3d plots, of course, or contour plots, where x-axis is subject, but these are sometimes hard to analyze visually. They may reveal important patterns though.
UPDATE You can also show the histogram of Kolmogorov-Smirnov statistics. It's true that the critical values are expensive to compute, but the statistics itself is easy to compute. So, you can obtain KS-stat for each subject, and show the histogram of obtained values. This will give you a great visual cue as to how the gamma distribution fits in general. It's almost like bootstrapping.
I hope that understood your situation and the question correctly. Considering your data set's number of distributions, visual exploratory approaches (such as QQ plots, which you mentioned) are not feasible in this case. Therefore, you have to resort to analytical approaches, such as goodness-of-fit (GoF) tests, as some have already mentioned in the comments above.
Since you have informed that distribution parameters are estimated from data, I assume that you have used or plan to use one of distribution fitting approaches. One of the most popular fitting approaches (along with least squares, to a lesser degree) is maximum likelihood estimation (MLE), which is generally easy to perform, for example, using function
MASS. However, depending on your particular data, fitting via
fitdistr() might not be so trivial. Some people prefer
fitdistrplus, as they consider it more advanced or useful.
After this straightforward step, you need to validate the estimation results, using one or more of the following GoF tests for continuous data (considering their pros and cons): chi-square (via binning), Kolmogorov-Smirnov (via corrected tables for critical values or Monte Carlo simulation, which I'm listing here just for completeness, as you are trying to avoid this), Anderson-Darling, Lilliefors, Cramér–von Mises and Watson. In terms of performance, the problem gets reduced to performing a relatively large number of non-parametric GoF tests, which IMHO is achievable either via doing it on a more powerful hardware (i.e., renting Amazon EC2 virtual instance), or via parallelizing code.
Returning to the essence of your question, my idea of possible approaches is to aggregate results either via bootstrapping (similarly to the one presented in this excellent answer), or some kind of averaging approach, similar to ensemble methods (for example, take a look at this research paper).