Mahalanobis distance in a matrix form The Mahalanobis distance between two vectors $x_i$ and $y_j$ is given by:
$$d_{ij}(x_i, y_j)^2 = (x_i-y_j)^TQ^{-1}(x_i-y_j)$$
Is there a vectorized way to represent the entries $d_{ij}$ in a matrix form $D$?
Here is my try:
\begin{equation}
\begin{split}
d_{ij}(x_i, y_j)^2 &= (x_i-y_j)^TQ^{-1}(x_i-y_j) \\
 & = \langle x_i-y_j, Q^{-1}(x_i-y_j) \rangle  \\
 & = \langle x_i, Q^{-1}(x_i-y_j) \rangle - \langle y_j, Q^{-1}(x_i-y_j) \rangle \\
 & = x_i^TQ^{-1}x_i + y_j^TQ^{-1}y_j - 2x_i^TQ^{-1}y_j
\end{split}
\end{equation}
but still have elements entries, not a matrix.
 A: Given two sets of vectors $\{x_i,\,y_j\},\,$ construct two matrices having these vectors as columns
$$\eqalign{
X &= \big[\matrix{x_1&x_2&\ldots&x_m}\big]
 &\in{\mathbb R}^{\ell\times m} \\
Y &= \big[\matrix{y_1&y_2&\ldots&y_n}\big]
 &\in{\mathbb R}^{\ell\times n} \\
}$$
Then the $\,m\times n\,$ matrix of the (squared) Euclidean distances between the vectors can be expressed in terms of either the vectors or the matrices
$$\eqalign{
&E_{ik} = \|x_i-y_k\|^2 &= (x_i-y_k)^T(x_i-y_k) \\
&E = \left(X\odot X\right)^TJ_Y &+ J_X^T\left(Y\odot Y\right) - 2X^TY \\
}$$
where $\odot$ denotes the Hadamard product and $(J_X,J_Y)$
are all-ones matrices the same size as $(X,Y)$, respectively.
The Cholesky factorization $\,Q^{-1}=LL^T\;$ can be used to modify the vectors and matrices
$$\eqalign{
x'_i &= L^Tx_i \qquad X'=L^TX \\
y'_k &= L^Ty_k \qquad Y'=L^TY \\
}$$
from which the matrix of (squared) Mahalanobis distances can be calculated.
To calculate the matrix specified in the question, apply the element-wise square root
$$E=D\odot D \quad\implies\quad D=E^{\odot 1/2}$$
A: Let $Z$ be a matrix whose $i$th column is $x_i - y_i$. Then $Z^\top Q^{-1} Z$ is a matrix whose $(i,j)$ entry is $d(x_i, y_j)^2$.
A: The expression given by @greg in his answer allows us to write the following expression $$E = (LX \odot LX)^TJ_Y + J_X^T(LY \odot LY) - 2X^TL^TLY$$
Here is my trial to calculate the derivative of E w.r.t to X, Y, and L.

Let $E_1 = (LX \odot LX)^TJ_Y$, $E_2 = J_X^T(LY \odot LY)$, and $E_3= 2X^TL^TLY$
For the $1^{st}$ term we have:
\begin{equation}
\begin{split}
d(E_1) &= d(LX \odot LX)^TJ_Y + (LX \odot LX)^TdJ_Y\\
& = (2(LX)^T \odot d((LX)^T))J_Y  \\
vec(d(E_1)) & = 2(J_Y \otimes I)(Kvec(LX)\odot ((I\otimes L)Kdx+ (I\otimes X^T)Kdl))\\
&= 2(J_Y \otimes I)(KDiag(vec(LX)))((I \otimes L)Kdx + (X^T \otimes I)Kdl)
\end{split} 
\end{equation}
Then:
\begin{equation}
\begin{split}
& \frac{d(vec(E_1))}{d(vec(X))}= 2(J_Y \otimes I)(KDiag(vec(LX)))((I \otimes L)K)\\
& \frac{d(vec(E_1))}{d(vec(Y))}= 0, \\
& \frac{d(vec(E_1))}{d(vec(L))}= 2(J_Y \otimes I)(KDiag(vec(LX)))(X^T \otimes I)K).
\end{split} 
\end{equation}
For the $2^{nd}$ term we have:
\begin{equation}
\begin{split}
d(E_2) &= d(J_X^T)(LY \odot LY) + J_X^Td(LY \odot LY)\\
& = J_X^T(2LY \odot (dLY + LdY)) \\
& = J_X^T(2LY \odot dLY + 2LY \odot LdY) \\
vec(d(E_2)) &= 2(I \otimes J_X^T)(vec(LY) \odot ((I \otimes L)d(vec(Y)) + (Y^T \otimes I)d(vec(L))) \\
&= 2(I \otimes J_X^T)(Diag(vec(LY))((I \otimes L)dy + (Y^T \otimes I)dl))
\end{split} 
\end{equation}
So
\begin{equation}
\begin{split}
&\frac{d(vec(E_2))}{d(vec(X))}=0, \\
&\frac{d(vec(E_2))}{d(vec(Y))} = 2(I \otimes J_X)^T(Diag(vec(LY))(I \otimes L), \\
&\frac{d(vec(E_2))}{d(vec(L))} = 2(I \otimes J_X)^T(Diag(vec(LY))(Y \otimes I)^T.
\end{split} 
\end{equation}
For the $3^{rd}$ term we have:
\begin{equation}
\begin{split}
d(E_3) & = 2(d(X^T)L^TLY + X^Td(L^TL)Y + X^TL^TLdY)\\
vec(d(E_3)) & = 2((Y^TL^TL \otimes I)d(vec(X^T) +  \\ & \quad (Y^T \otimes X^T)d(vec(L^TL) + (I \otimes X^TL^TL)d(vec(Y)))  \\
\end{split} 
\end{equation}
where, $ d(vec(L^TL)) = ((L^T \otimes I)K + (I \otimes L^T))d(vec(L))$
Thus:
\begin{equation}
\begin{split}
&\frac{d(vec(E_3))}{d(vec(X))}= 2(Y^TL^TL \otimes I)K = 2(LY \otimes I)^T(L \otimes I)K, \\
&\frac{d(vec(E_3))}{d(vec(Y))} = 2(I \otimes X^TL^TL) = 2(I \otimes LX)^T(I \otimes L), \\
&\frac{d(vec(E_3))}{d(vec(L))} = 2(Y^T \otimes X^T)((L^T \otimes I)K + (I \otimes L^T)).
\end{split} 
\end{equation}
Now, by putting the three terms together we obtain
$dE = dE_1 + dE_2 + dE_3$
Holding Y and L constant (i.e. $dY=0$, $dL=0$) yields:
\begin{equation}
\begin{split}
\frac{d(vec(E))}{d(vec(X))} & = 2(J_Y \otimes I)(KDiag(vec(LX))((I \otimes L)K) \\
& - 2(LY \otimes I)^T(L \otimes I)K
\end{split}
\end{equation}
Holding X and L constant (i.e. $dX=0$, $dL=0$) yields:
\begin{equation}
\begin{split}
\frac{d(vec(E))}{d(vec(Y))} & = 2(I \otimes J_X)^T(Diag(vec(LY))(I \otimes L) -  2(I \otimes LX)^T(I \otimes L) \\
& = 2((I \otimes J_X)^T(Diag(vec(LY)) - (I \otimes LX)^T)(I \otimes L)
\end{split}
\end{equation}
Holding X and Y constant (i.e. $dX=0$, $dY=0$) yields:
\begin{equation}
\begin{split}
\frac{d(vec(E))}{d(vec(L))} & = 2(J_Y \otimes I)(KDiag(vec(LX))(X^T \otimes I)K) \\
& + 2(I \otimes J_X)^T(Diag(vec(LY))(Y \otimes I)^T \\
& - 2(Y^T \otimes X^T)((L^T \otimes I)K + (I \otimes L^T)). 
\end{split} 
\end{equation}
