Fast top-n from samples of many different normal distributions Let's say I have 1 million normal distributions, each with a different mean and stdev. I want to sample from each distribution and take the top 10 samples (for a Thompson sampling application.)
Is there a faster approximate way to sample than generating all million samples and taking the top n?
This process gets repeated many times, so I am happy to precompute something if it helps with the amortized cost.

For identically distributed values there is an extremely efficient method.  For differently distributed values, though, you're out of luck unless there is perhaps some special structure to the means and variances.  There aren't even analytical solutions to the problem once $N$ exceeds $2.$  (See my answer at https://stats.stackexchange.com/a/579392/919 for the case $N=2.$)
For the record, and because it might help inspire hybrid solutions, here is the method for iid values.

Consider any continuous distribution function $F,$ such as a Normal CDF, and let $f$ be its density. Given $N\ge 1$ iid values from $F$ and $d$ between $1$ and $N,$ let $X_{[1,N]}\ge X_{[2,N]} \ge \cdots \ge X_{[d,N]}$ be its top $d$ values in descending order.
The joint distribution of $X = X_{[d+1,N]}$ and $Y = X_{[d,N]}$ is continuous with density
$$g_d(x,y) = C(N,d) F(x)^{N-d-1} f(x) f(y) (1-F(y))^{d-1}\,\mathcal{I}(x\le y)$$
where $C$ normalizes the distribution (to integrate to unity) and $\mathcal{I}$ is the indicator function.  It follows (from the definition of conditional distributions) that the conditional density of $X$ given $Y=y$ is proportional to
$$g_{d;X\mid Y}(x\mid y) \propto F(x)^{N-d-1} f(x)\mathcal{I}(x\le y) \propto\frac{\mathrm d}{\mathrm d x} F(x)^{N-d}\mathcal{I}(x\le y).$$
This makes the sequence $X_{[1,N]}, X_{[2,N]}, \ldots, X_{[d,N]}$ a Markov Chain with easily computable transitions.  They are easy to compute because, to generate a random variable from $g_{d;X\mid Y}$ given $Y=y,$ all we need to is generate a random uniform variable $U$ in the range $(0, 1)$ and find the solution $x$ to $F(x) = U^{1/(N-d)}F(y).$  For many distributions, such as the Normal ones, this solution is readily computed using the percentage point function $F^{-1}.$  That is,

Given that $y$ is a realization of the $d^\text{th}$ largest value of $N$ iid values, a realization of the $d+1^{\text{st}}$ largest value is given by $$x\mid y = F^{-1}\left(U^{1/(N-d)} F(y)\right).$$

Starting the chain with $y=\infty,$ where $F(y)=1,$ and proceeding for $d$ steps gives a realization of the top $d$ values.

Here, as a concrete example, is an implementation of this algorithm in R.  It returns an array where each row gives a realization of $(X_{[1,N]}, X_{[2,N]}, \ldots, X_{[d,N]}).$
rhighest <- function(n, N = 1, d = 1, pf = pnorm, qf = qnorm) {
  q <- rep(1, n)
  sapply(seq(N, N-d+1), function(k) qf(q <<- (runif(n)^(1/k) * q)))
}

To generate, say, 100,000 independent realizations of the top $d=10$ of $N=10^6$ Normal variables takes a small fraction of a second:
X <- rhighest(1e5, 1e6, 10)

A histogram of, say, its tenth column shows the distribution of the tenth highest:
hist(X[, 10], freq=FALSE)


A scatterplot of its first and tenth columns show the extent to which those variables are correlated.  (The diagonal line marks the locus of equal values; the points must lie beneath it.)

To check the correctness of this algorithm, I generated 20,000 sets of $N=10$ iid Normal variates, sorted each set, and picked out the $d=1,2,4,$ and $8$ highest in each set (a "brute force" simulation).  I also used rhighest for comparison (the "efficient" simulation).  Here are the QQ plots of their quantiles:
N <- 10
n.sim <- 2e4
d <- 8
set.seed(17)
X <- t(apply(matrix(rnorm(N*n.sim), N), 2, sort, decreasing = TRUE))
Y <- rhighest(n.sim, N, d)

qs <- seq(0, 1, length.out=1e3) # Quantiles for QQ plots
for (d in c(1, 2, 4, 8)) {
  xlim <- range(c(X[, d], Y[, d]))
  plot(quantile(X[, d], qs), quantile(Y[, d], qs), xlim=xlim, ylim=xlim, asp=1,
       cex=0.75, xlab="Brute Force", ylab="Efficient",
       main=bquote(paste("QQ Plot for ", d==.(d))))
  abline(0:1, col=hsv(0.01, 0.6, 0.8))
}


They agree closely.
