Globally-optimized nearest neighbors lookup
Here's one approach -- not necessarily the fastest.  This can be viewed as an instance of the assignment problem, or equivalently, finding a minimum weight perfect matching of a bipartite graph.  Imagine a bipartite graph with $\mathcal{A}$ being the set of vertices on the left, $\mathcal{B}$ being the set of vertices on the right, and an edge between every pair of vertices  whose weight is equal to the distance between those points.  Now you are looking for a perfect matching of minimum weight.  Once the graph is built, there are standard algorithms for this, such as the Hungarian algorithm or linear programming.

The first step in such an approach is to compute all-pairs distances, i.e., the distance between every point in $\mathcal{A}$ and every point in $\mathcal{B}$.  This takes $O(nma)$ time.  Based on your original question, it appears that $n \approx 20$, $m \approx 50,000$, and $a \approx 800$, so that part should be very fast.

Unfortunately, the Hungarian algorithm might be too slow.  The running time in this setting will be $O((n+m)^3)$, which might be too slow for your needs.  I think a tighter analysis shows that its running time is actually $O(n^3m)$, which might be fast enough for your needs, but I'm not certain about that.

If the Hungarian algorithm is too slow, here is a heuristic that you could try.  In each stage, we'll have a non-optimal solution and iteratively try to improve it.  (You can start with the greedy solution, to get the process started.)  Suppose our current solution assigns $f(b_j)$ to the point $b_j$.  Let $t$ be a threshold.  Filter the graph to include only edges $(a_i,b_j)$ where $\Delta(a_i,b_j) \le \Delta(f(b_j),b_j) + t$, and delete all isolated left vertices.  Now find a minimum weight perfect matching in that graph.  That will be a new solution that is at least as good as the previous solution.  In each stage, increase $t$ according to some schedule.  If you keep $t$ small, this won't necessarily give you the optimal solution, but it will make the Hungarian algorithm more efficient because it works with a smaller graph.

Theoretically, you can find a value for $t$ that enables you to prove that this will find the optimal solution (namely, a lower bound for the sum can be obtained by matching each element $b_j$ to its nearest neighbor in $\mathcal{A}$, ignoring the constraint that the elements of $\mathcal{C}$ be unique, and summing those distances; now if you let $t$ be at least the difference between the cost of the current solution and that lower bound, and apply this scheme, if it converges -- yields the same solution you started with -- then you know you've found the optimal solution).  However I don't know whether that will be useful in practice.

Anyway, this is one approach.  It might not be the best or fastest, because it doesn't take into account the geometrical structure.  But it's something you could try, to see if it meets your needs.