I have a spline-based model where I have a set of control points in 3D space. I want to compute pairwise distances between every point on a regular 3D grid and these control points. Is there an efficient algorithm to do this in terms of space (i.e. not having to calculate the full distance matrix)?
One thing is that the number of control points is much less than the grid points, (e.g. by 3 orders of magnitude). Approximations can be considered. Bonus points if the algorithm is differentiable.
$ \def\bbR#1{{\mathbb R}^{#1}} \def\o{{\tt1}} \def\L{\left}\def\R{\right} \def\LR#1{\L(#1\R)} \def\BR#1{\Big(#1\Big)} \def\diag#1{\operatorname{diag}\LR{#1}} \def\Sym#1{\operatorname{Sym}\!\BR{#1}} \def\qiq{\quad\implies\quad} \def\ddt{\frac{d}{dt}} \def\dd#1{\frac{d #1}{dt}} $ The matrices $$\eqalign{ &X = \left[\matrix{x_1&x_2&\ldots&x_m}\right]&\in \bbR{\ell\times m} \\ &Y = \left[\matrix{y_1&y_2&\ldots&y_n}\right]&\in \bbR{\ell\times n} \\ }$$ have a different number of columns, but the same number of rows $\big($in your case $\ell=3\big)$.
The distance (squared) between the columns of these matrices in component form is $$\eqalign{ D_{ik}^2 &= \|x_i - y_k\|^2 \;=\; (x_i-y_k)^T(x_i-y_k) \\ }$$ or in matrix form $$\eqalign{ D\odot D \;=\; (X\odot X)^TJ_Y + J_X^T(Y\odot Y) - 2X^TY \;\in\; \bbR{m\times n} \\ }$$ where $\{J_X,J_Y\}$ are all-ones matrices with the same dimensions as $\{X,Y\}$ and $\{\odot\}$ denotes the elementwise/Hadamard product.
This can also be written using diag operators and all-ones vectors. $$\eqalign{ g &= \diag{X^TX},\quad w&=\o_m &\in \bbR{m} \\ h &= \diag{Y^TY},\quad u&=\o_n &\in \bbR{n} \\ }$$ $$\eqalign{ D\odot D &= gu^T + wh^T - 2X^TY \\\\ }$$
The above expressions are easy to differentiate, e.g. $$\ddt\big(D\odot D\big) = 2D\odot\dd D$$
