There are two real matrices: $X, Y$ with $X$ being of dimension $n_1$ x $p$, $Y$ of dimension $n_2$ x $p$.
The goal is to form the matrix $D$ of dimension $n_1$ x $n_2$ where each element $d_{ij}$ is computed as the (L2-) norm of the difference of row i of X and row j of Y:
$$d_{ij} = \lVert x_i - y_j \rVert_2^2$$
My question is how to derive a matrix formula for $D$. I.e
$$D = g(X,Y)$$
Here is my attempt at starting: I can expand $\lVert v\rVert_2^2$ as $v^Tv$ where $v = x_i - y_j$ to obtain:
$$d_{ij} = x_i^Tx_i + y_j^Ty_j - 2x_i^Ty_j$$
But I'm not sure how to proceed.
You just need to "matricize" your result a bit.
Let $X\in {\mathbb R}^{p\times n_1}$, $Y\in {\mathbb R}^{p\times n_2}$, $D\in {\mathbb R}^{n_1\times n_2}$, $v_1\in {\mathbb R}^{n_1}$, and $v_2\in {\mathbb R}^{n_2}$; where the elements of the $v$-vectors are all ones.
Then $$ D = {\rm diag}(X^TX)\,v_{2}^T + v_{1}\,({\rm diag}(Y^TY))^T - 2 X^TY $$ Note that I defined my $X,Y$ as the transpose of yours.
(My brain works better in terms of column-vectors.)