Stable and fast computation of the squared euclidean distance matrix

Question

Let's say I want to compute the matrix $M$ of the squared euclidean distances between each pair of vectors $(x, y)$ belonging to two sets $X$ and $Y$ respectively. The sets of vectors $X$ and $Y$ have size $m$ and $n$. Each vector has $k$ floating point coordinates.

As far as I can tell, there are mainly two ways to implement the computation of an euclidean distance matrix. The direct formula, which is slow but accurate, and the expanded formula which can be fast but numerically inaccurate. I'm wondering if there would be a way to compute the squared euclidean distance matrix in a stable and fast way.

The direct formula compute for each pair of vector $x \in X$ and $y \in Y$: $\sum_i (x_i - y_i)^2$. It is pretty accurate and could even be combined a with compensated summation algorithm. It is however relatively slow as this is an $\mathcal{O}(mnk)$ algorithm and can only take advantage of level 1 BLAS routines.

On the other hand, the expanded formula compute $\sum_i{x_i^2} + \sum_i{y_i^2} - 2\sum_i{x_i y_i}$. It is pretty fast as it only makes 2/3 of the operations of the direct formula and can use the BLAS matrix multiplication for the last term. Unfortunately, it is numerically very unstable and often produce negative result.

Intuitively, the computation of an euclidean distance matrix seems very redundant when $k \ll min(m, n)$. Once we know the distance from a vector $x \in X$ to $k+1$ vectors from $Y$, then the distance from $x$ to all the other vectors $y \in Y$ are pretty much determined. But I can't find a useful way to exploit this redundancy to make less computations.

Side note: the expanded formula also gives us an upper bound on the minimal reachable complexity of the computation of the euclidean matrix. Since its complexity is basically that of the matrix multiplication, it is theoretically possible to compute the euclidean distance matrix in $\mathcal{O}(n^{2.3737})$. But maybe there's a structure to the squared euclidean distance matrix that would make this bound practical. Who knows?