How can I use matrix operation to represent the following result? I need to derive a pairwise offset matrix from the two known data array $X_{n \times p}$ and $Y_{m \times p}$, in which the offset are derived from values at corresponding locations $x_{ij}$ and $y_{ij}$ in the two array. The following figure shows an example:

How can I derive the offset matrix $D_i$ using matrix operations (vectorized code)? I wll also appreciate that you can give me an example using a code platform, such as Python, R or Matlab. Thank you.

You can get all the columns of $D_1$ in a list by doing
cols <- lapply(1:m, function(i) X[1, ] - Y[i, 1])

Then you bind them by doing do.call(cbind, cols).
Similar for $D_2$.
