Nearest-center assignment without an N x K x D distance tensor

Given N points and K centers in D dimensions, assign each point to its nearest center, as in the assignment step of k-means. The obvious broadcast builds an N x K x D array of differences, which runs out of memory at scale. Rewrite it so it never builds that array (for example by expanding the squared distance into dot products, or by working in chunks), and discuss memory, speed and numerical accuracy.