Pairwise distance matrix
Solution 1:
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. $$
Solution 2:
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)$$