There is a simple efficient solution. It uses the ideas common to all rank-sum tests, such as the Wilcoxon tests. This answer derives the solution and provides an R implementation.
The code in the question simulates data that have a vanishingly small chance of exhibiting any ties at all between a sample of group $a$ and a sample of group $b$, so let's assume there exist no ties.
Let there be $m$ elements in $a$ (and they can have ties among each other) and $n$ elements in $b$ (which also may have ties). Let $A$ be the random variable modeled by drawing one element randomly and uniformly from group $A$ and similarly let $B$ be the random variable for one draw from group $B$. The desired value (as I interpret the question) is the chance that $A$ exceeds $B$.
Notice that the test of whether a realization of $A$ exceeds one of $B$ is a simple comparison. Thus, the problem is unchanged if we replace all elements of $a$ and $b$ by their indexes when the two sets are sorted in increasing order. These indexes are their ranks, provided that ties are resolved in some arbitrary manner (that is, do not average the ranks of any groups of ties).
For example, let $a = (0,2)$ and $b = (1,1,3)$. Sorting the combined two (multi)sets gives the sequence $(0,1,1,2,3)$. The indexes of the values coming from $a$ are $1$ and $4$ while the indexes of the values from $b$ are $2, 3,$ and $5$.
Compute the chance $\Pr(A\gt B)$ by summing over the possible values of $A$, each of which has the probability $1/m$. Let the ranks of these values be $r_1 \lt r_2 \lt \cdots \lt r_m$. Suppose $r_i$ is picked. Then the chance, conditional on this selection, that $B$ has a smaller value equals the number of smaller values in $b$ divided by $n$. The number of smaller values altogether in both $a$ and $b$ is, by definition, $r_i-1$, but we know exactly $i-1$ of them (namely, $r_1, r_2, \ldots, r_{i-1}$) are in $a$. Thus
$$\Pr(A \gt B | A = r_i) = \frac{1}{n}\left(r_i-1 - (i-1)\right) = \frac{1}{n}\left(r_i-i\right)$$
entailing
$$\Pr(A\gt B) = \sum_{i=1}^m \Pr(A\gt B | A=r_i)\Pr(A=r_i) = \frac{1}{mn}\sum_i \left(r_i-i\right).$$
In the example, $\Pr(A\gt B) = \frac{1}{2\times 3}((1-1) + (4-2)) = \frac{2}{6}$ and (reversing the roles of $a$ and $b$ as a quick check) $\Pr(B\gt A) = \frac{1}{3\times 2}((2-1) + (3-2) + (5-3)) = \frac{4}{6} = 1 - \frac{2}{6}$ as one would expect.
This calculation (when implemented as a general-purpose algorithm) requires sorting all $m+n$ values to find their ranks and then summing either $m$ or $n$ values (for efficiency, one would pick whichever is smaller). Therefore the computational burden is $O((m+n)\log(m+n)),$ and can be reduced to $O(\min(m\log(n), n\log(m)))$ when the larger of $a$ and $b$ is already sorted. That's pretty efficient.
When there are ties between elements of $a$ and $b$, the idea to reduce the question to rank sums and then compute a sum over conditional probabilities still works, but the calculations of the conditional probabilities get more complicated.
R Code
R will have trouble with calculations that overflow its integer data type. The following solution handles that possibility.
prob <- function(a, b, ...) {
# Returns chance that a random sample of `a` will exceed one of `b`
# (assuming no ties between elements of `a` and `b`)
# Optional args are passed to `rank` to control handling of NAs and
# how to resolve any ties.
m <- length(a); n <- length(b)
if (m < n) {
r <- rank(c(a,b), ...)[1:m] - 1:m
} else {
r <- rank(c(a,b), ...)[(m+1):(m+n)] - 1:n
}
s <- ifelse ((n+m)^2 > 2^31, sum(as.double(r)), sum(r)) / (as.double(m)*n)
return (ifelse(m < n, s, 1-s))
}
To emulate the data in the question, let's simulate sets of normally distributed values. When the elements of $a$ come from a normal distribution with mean $\mu$ and standard deviation $\sigma$ and those of $b$ come from a normal distribution with mean $\nu$ and SD $\sigma$, we may analytically compute that $\Pr(A\gt B)$ (prior to simulating the elements) equals $\Phi(\frac{\mu-\nu}{\sigma\sqrt{2}})$ where $\Phi$ is the cumulative standard normal distribution function. This enables us to test prob, as in the following:
set.seed(17)
m <- 10^6; n <- 10^4
mu.a <- 0; mu.b <- -2
a <- rnorm(m, mu.a)
b <- rnorm(n, mu.b)
system.time(print(prob(a,b), digits=5))
system.time(print(1 - prob(b,a), digits=5))
print(pnorm((mu.a - mu.b)/sqrt(2)), digits=5)
The output is
> system.time(print(prob(a,b), digits=5))
[1] 0.92124
user system elapsed
0.51 0.00 0.51
> system.time(print(1 - prob(b,a), digits=5))
[1] 0.92124
user system elapsed
0.51 0.00 0.52
> print(pnorm((mu.a - mu.b)/sqrt(2)), digits=5)
[1] 0.92135
It shows that the computation time does not depend on the order in which a and b are provided to prob. The computation time of $1/2$ second is reasonably quick (for over one million numbers). The close agreement of $0.92124$ and $0.92135$ is evidence in favor of the correctness of this solution.
This solution can easily be iterated over groups using the usual R idioms for looping.
Rcode! -- What isdata.table? There are no variablesaorbto make sense of the second through fourth lines, either. $\endgroup$