Pairwise distance matrix Suppose $x_1,\dots,x_n \in \mathbb{R}^d$. Is there a vectorized way of representing the square distance matrix $D_{ij} = \lVert x_i - x_j \rVert^2$?
 A: We could make use of $\|x_i-x_j\|^2=\|x_i\|^2+\|x_j\|^2-2(x_i \cdot x_j)$, the Gramian matrix $G=[\langle x_i, x_j\rangle]$, its diagonal elements vector $g=\operatorname{diag}(G)$, and the vector of all ones $\mathbf{1}$ to write
$$
D = g\mathbf{1}^T + \mathbf{1}g^T - 2G.
$$
A: The first answer gives a nice formula for the distance matrix using the Gramian matrix.
Here I will develop the answer a little more.
We can write the distance between two vectors $x_i$ and $x_j$ in terms of their inner product as follow:
$$ d(x_i,x_j)^2 = \langle x_i-x_j, x_i-x_j\rangle = \langle x_i,x_i\rangle + \langle x_j, x_j\rangle - 2\langle x_i, x_j\rangle \\ 
= \|x_i\|^2 + \|x_j\|^2 -2\langle x_i, x_j\rangle \quad\qquad\qquad$$
In matrix form, we have $$D = [d(x_i, x_j)^2] $$
$$= \begin{pmatrix}
0 & \langle x_1,x_1\rangle + \langle x_2, x_2\rangle - 2\langle x_1, x_2\rangle & \cdots & \langle x_1,x_1\rangle + \langle x_n, x_n\rangle - 2\langle x_1, x_n\rangle \\
\langle x_2,x_2\rangle + \langle x_1, x_1\rangle - 2\langle x_2, x_1\rangle & 0 & \cdots & \langle x_2,x_2\rangle + \langle x_n, x_n\rangle - 2\langle x_2, x_n\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_n,x_n\rangle + \langle x_1, x_1\rangle - 2\langle x_n, x_1\rangle & \langle x_n,x_n\rangle + \langle x_2, x_2\rangle - 2\langle x_n, x_2\rangle & \cdots & 0 
\end{pmatrix} $$
Which can be written as
$$D =
\begin{pmatrix}
\langle x_1,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_n, x_n\rangle \\
\langle x_1,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_n, x_n\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_1,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_n, x_n\rangle 
\end{pmatrix} 
+
\begin{pmatrix}
\langle x_1,x_1\rangle &  \langle x_1, x_1\rangle & \cdots & \langle x_1, x_1\rangle \\
\langle x_2,x_2\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_2, x_2\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_n,x_n\rangle &  \langle x_n, x_n\rangle & \cdots & \langle x_n, x_n\rangle 
\end{pmatrix}
- 2
\begin{pmatrix}
\langle x_1,x_1\rangle &  \langle x_1, x_2\rangle & \cdots & \langle x_1, x_n\rangle \\
\langle x_2,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_2, x_n\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_n,x_1\rangle &  \langle x_n, x_2\rangle & \cdots & \langle x_n, x_n\rangle 
\end{pmatrix}
$$
Which has the form $$D = diag(G)1_n^T + 1_ndiag(G)^T - 2G,$$
where $$ G= 
\begin{pmatrix}
\langle x_1,x_1\rangle &  \langle x_1, x_2\rangle & \cdots & \langle x_1, x_n\rangle \\
\langle x_2,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_2, x_n\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_n,x_1\rangle &  \langle x_n, x_2\rangle & \cdots & \langle x_n, x_n\rangle 
\end{pmatrix}
, diag(G) = 
\begin{pmatrix}
\langle x_1,x_1\rangle \\
\langle x_2,x_2\rangle \\
\vdots  \\
\langle x_n,x_n\rangle
\end{pmatrix}
, 1_n=\begin{pmatrix}
1 \\
1 \\
\vdots  \\
1
\end{pmatrix}$$
D can be given also as $$D = 1_n^Tdiag(G) + diag(G)^T1_n - 2G,$$
where $$ G= 
\begin{pmatrix}
\langle x_1,x_1\rangle &  \langle x_1, x_2\rangle & \cdots & \langle x_1, x_n\rangle \\
\langle x_2,x_1\rangle &  \langle x_2, x_2\rangle & \cdots & \langle x_2, x_n\rangle \\
\vdots  & \vdots  & \ddots & \vdots  \\
\langle x_n,x_1\rangle &  \langle x_n, x_2\rangle & \cdots & \langle x_n, x_n\rangle 
\end{pmatrix}
, \\
diag(G) = 
\begin{pmatrix}
\langle x_1,x_1\rangle & \langle x_2,x_2\rangle &\cdots & \langle x_n,x_n\rangle
\end{pmatrix}
, \\ 1_n=\begin{pmatrix}
1 & 1 \cdots  1
\end{pmatrix}$$
The last point, if the vectors $x_i, x_j$ are normalized, that is
$$\langle x_i,x_i\rangle=\langle x_j,x_j\rangle=1 $$
We get
$$D = 2(1_{nxn}- G)$$
