There are several answers on similar questions on StackExchange (see for example Pairwise Mahalanobis distance in R).

Centering the data

Let us assume your matrices are "data matrices", i.e, they have the dimensions $n \times N$ where $n$ is the number of observations in a sample and $N$ is the number of variables.

Ideally you should center the data matrices by subtracting the column means from the columns before forming your covariance matrix and inverting it. Let us assume that, in your example, two samples are $$ \mathbf{X} = \begin{bmatrix} 0.678177 & 0.989365 & 0.558944\\ 0.652491 & 0.799086 & 0.514262\\ 0.230560 & 0.299587 & 0.023458 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} 0.015513 & 0.932740 & 0.351038 \\ 0.400391 & 0.515025 & 0.179851 \end{bmatrix} $$ Let us center the data by subtracting the sample means which are $$ E[\mathbf{X}] = \begin{bmatrix} 0.52041 & 0.69601 & 0.36555 \end{bmatrix} $$ and $$ E[\mathbf{Y}] = \begin{bmatrix} 0.20795 & 0.72388 & 0.26544\end{bmatrix} $$ The centered data are $$ \mathbf{X} = \begin{bmatrix} 0.15777 & 0.29335 & 0.19339\\ 0.13208 & 0.10307 & 0.14871\\ -0.28985 & -0.39643 & -0.34210 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594\\ 0.192439 & -0.208857 & -0.085594 \end{bmatrix} $$

Covariance

Then the cross-covariance matrix can be computed using $$ \begin{aligned} \mathbf{S}(\mathbf{X},\mathbf{Y}) & = E\left[(\mathbf{X} - E[\mathbf{X}]) \otimes (\mathbf{Y} - E[\mathbf{Y}])\right] \\ & = E\left[\mathbf{X}\otimes\mathbf{Y} - E[\mathbf{X}]\otimes\mathbf{Y} - \mathbf{X}\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\right] \\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] \end{aligned} $$ In terms of matrix components $$ S_{ij} = E[X_i\,Y_j] - E[X_i]\,E[Y_j] $$ If we assume that the expected value can be estimated by the sample mean, then we can use $$ E[\mathbf{X}] = \frac{1}{n_x}\sum_{i=1}^{n_x} \mathbf{X}^{(i)} ~,~~ E[\mathbf{Y}] = \frac{1}{n_y}\sum_{i=1}^{n_y} \mathbf{Y}^{(i)} $$ and, with $N$ being the number of variables (columns of the two matrices), $$ E[\mathbf{X} \otimes \mathbf{Y}] = \frac{1}{N-1}\sum_{i=1}^N (\mathbf{X} \otimes \mathbf{Y})^{(i)} $$ If the matrices have been centered, $E[X_i] = E[Y_j] = 0$ and we have $$ S_{ij} = E[X_i\,Y_j] \quad \implies \mathbf{S}(\mathbf{X},\mathbf{Y}) = E[\mathbf{X} \otimes \mathbf{Y}]\,. $$

Pooled covariance

The correlations among the data in each sample are estimated using the sample covariance and an estimate of the population covariance is computed using the pooled covariance.

The pooled covariance is given by $$ \mathbf{S} = \frac{n_x}{n_x+n_y} \mathbf{S}(\mathbf{X},\mathbf{X}) + \frac{n_y}{n_x+n_y} \mathbf{S}(\mathbf{Y},\mathbf{Y}) $$

Covariance example

In the example, if $$ \mathbf{X} = \begin{bmatrix} \mathbf{X}^{(1)} \\ \mathbf{X}^{(2)} \\ \mathbf{X}^{(3)} \end{bmatrix} = \begin{bmatrix} [X_1,X_2,X_3]^{(1)} \\ [X_1,X_2,X_3]^{(2)} \\ [X_1,X_2,X_3]^{(3)} \end{bmatrix} \quad \text{and} \quad \mathbf{Y} = \begin{bmatrix} \mathbf{Y}^{(1)} \\ \mathbf{Y}^{(2)} \end{bmatrix} = \begin{bmatrix} [Y_1,Y_2,Y_3]^{(1)} \\ [Y_1,Y_2,Y_3]^{(2)} \end{bmatrix} $$ Then, $$ (\mathbf{X} \otimes \mathbf{X})^{(1)} = (\mathbf{X}^{(1)})^T \cdot \mathbf{X}^{(1)} = \begin{bmatrix} 0.024891 & 0.046281 & 0.030511\\ 0.046281 & 0.086056 & 0.056731\\ 0.030511 & 0.056731 & 0.037399\end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(2)} = (\mathbf{X}^{(2)})^T \cdot \mathbf{X}^{(2)} = \begin{bmatrix} 0.017446 & 0.013614 & 0.019642\\ 0.013614 & 0.010624 & 0.015328 \\ 0.019642 & 0.015328 & 0.022114 \end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(3)} = (\mathbf{X}^{(3)})^T \cdot \mathbf{X}^{(3)} = \begin{bmatrix} 0.084013 & 0.114904 & 0.099157\\ 0.114904 & 0.157153 & 0.135616 \\ 0.099157 & 0.135616 & 0.117030 \end{bmatrix} $$ Adding these and dividing by $N - 1 = 2$ gives $$ \mathbf{S}(\mathbf{X},\mathbf{X}) = \begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} $$ Note that $$ \mathbf{X}^T \cdot \mathbf{X} = (\mathbf{X} \otimes \mathbf{X})^{(1)} + (\mathbf{X} \otimes \mathbf{X})^{(2)} + (\mathbf{X} \otimes \mathbf{X})^{(3)} \,. $$ Similarly, $$ \mathbf{S}(\mathbf{Y},\mathbf{Y}) = \begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$

Pooled covariance example

Then the pooled covariance is $$ \mathbf{S} = \frac{3}{5}\begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} + \frac{2}{5}\begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$ Therefore $$ \mathbf{S} = \begin{bmatrix} 0.067531 & 0.020286 & 0.031615 \\ 0.020286 & 0.111047 & 0.076604 \\ 0.031615 & 0.076604 & 0.058824 \end{bmatrix} $$

Inverse of the pooled covariance matrix

Now compute the inverse of $\mathbf{S}$ which clearly has to be a square matrix. This is the most computationally intensive part of the calculation and requires special attention for high dimensional data. $$ \mathbf{S}^{-1} = \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix} $$

Mahalanobis distance

Now to find the distance (squared) between $\mathbf{X}^{(i)}$ and $\mathbf{Y}^{(j)}$ we use the relation $$ d^2_{ij} = (\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})\cdot\mathbf{S}^{-1} \cdot(\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})^T $$ First compute $\mathbf{Z} := \mathbf{X}^{(i)} - \mathbf{Y}^{(j)}$. Both vectors have to be the same length for this operation to work. Define $\mathbf{T} := \mathbf{S}^{-1}$. Then $$ d^2_{ij} = \mathbf{Z} \cdot \mathbf{T} \cdot \mathbf{Z}^T $$ In terms of indices, with $N$ as the length of the vector (number of variables), $$ d^2_{ij} = \sum_{i=1}^N \sum_{j=1}^N Z_i T_{ij} Z_j $$ That should give you the distances you seek.

Mahalanobis distance example

Let us find the distance between the vectors $$ \mathbf{X}^{(2)} = \begin{bmatrix} 0.13208 & 0.10307 & 0.14871 \end{bmatrix} $$ and $$ \mathbf{Y}^{(1)} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594 \end{bmatrix} $$ We have, $$ \mathbf{Z} = \mathbf{X}^{(2)} - \mathbf{Y}^{(1)} = \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} $$ Therefore, $$ \begin{aligned} (d_{21})^2 &= \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix}\begin{bmatrix} 0.324521 \\ -0.105784 \\ 0.063114 \end{bmatrix} \\ &= 3.4722 \,. \end{aligned} $$

Samples with one observations each

Merge the two samples into one and just use the sample covariance instead of the pooled covariance. The result will not make much dense and a simpler distance measure is preferable.

There are several answers on similar questions on StackExchange (see for example Pairwise Mahalanobis distance in R).

Centering the data

Let us assume your matrices are "data matrices", i.e, they have the dimensions $n \times p$ where $n$ is the number of observations in a sample and $p$ is the number of variables.

Ideally you should center the data matrices by subtracting the column means from the columns before forming your covariance matrix and inverting it. Let us assume that, in your example, two samples are $$ \mathbf{X} = \begin{bmatrix} 0.678177 & 0.989365 & 0.558944\\ 0.652491 & 0.799086 & 0.514262\\ 0.230560 & 0.299587 & 0.023458 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} 0.015513 & 0.932740 & 0.351038 \\ 0.400391 & 0.515025 & 0.179851 \end{bmatrix} $$ Let us center the data by subtracting the sample means which are $$ E[\mathbf{X}] = \begin{bmatrix} 0.52041 & 0.69601 & 0.36555 \end{bmatrix} $$ and $$ E[\mathbf{Y}] = \begin{bmatrix} 0.20795 & 0.72388 & 0.26544\end{bmatrix} $$ The centered data are $$ \mathbf{X} = \begin{bmatrix} 0.15777 & 0.29335 & 0.19339\\ 0.13208 & 0.10307 & 0.14871\\ -0.28985 & -0.39643 & -0.34210 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594\\ 0.192439 & -0.208857 & -0.085594 \end{bmatrix} $$

Covariance

Then the cross-covariance matrix can be computed using $$ \begin{aligned} \mathbf{S}(\mathbf{X},\mathbf{Y}) & = E\left[(\mathbf{X} - E[\mathbf{X}]) \otimes (\mathbf{Y} - E[\mathbf{Y}])\right] \\ & = E\left[\mathbf{X}\otimes\mathbf{Y} - E[\mathbf{X}]\otimes\mathbf{Y} - \mathbf{X}\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\right] \\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] \end{aligned} $$ In terms of matrix components $$ S_{ij} = E[X_i\,Y_j] - E[X_i]\,E[Y_j] $$ If we assume that the expected value can be estimated by the sample mean, then we can use $$ E[\mathbf{X}] = \frac{1}{n_x}\sum_{i=1}^{n_x} \mathbf{X}^{(i)} ~,~~ E[\mathbf{Y}] = \frac{1}{n_y}\sum_{i=1}^{n_y} \mathbf{Y}^{(i)} $$ and, with $N = n_x + n_y$ if the $\mathbf{X}$ and $\mathbf{Y}$ are different matrices and $N = n_x$ if $\mathbf{X} = \mathbf{Y}$, $$ E[\mathbf{X} \otimes \mathbf{Y}] = \frac{1}{N-1}\sum_{i=1}^N (\mathbf{X} \otimes \mathbf{Y})^{(i)} $$ If the matrices have been centered, $E[X_i] = E[Y_j] = 0$ and we have $$ S_{ij} = E[X_i\,Y_j] \quad \implies \mathbf{S}(\mathbf{X},\mathbf{Y}) = E[\mathbf{X} \otimes \mathbf{Y}]\,. $$

Pooled covariance

The correlations among the data in each sample are estimated using the sample covariance and an estimate of the population covariance is computed using the pooled covariance.

The pooled covariance is given by $$ \mathbf{S} = \frac{n_x}{n_x+n_y} \mathbf{S}(\mathbf{X},\mathbf{X}) + \frac{n_y}{n_x+n_y} \mathbf{S}(\mathbf{Y},\mathbf{Y}) $$

Covariance example

In the example, if $$ \mathbf{X} = \begin{bmatrix} \mathbf{X}^{(1)} \\ \mathbf{X}^{(2)} \\ \mathbf{X}^{(3)} \end{bmatrix} = \begin{bmatrix} [X_1,X_2,X_3]^{(1)} \\ [X_1,X_2,X_3]^{(2)} \\ [X_1,X_2,X_3]^{(3)} \end{bmatrix} \quad \text{and} \quad \mathbf{Y} = \begin{bmatrix} \mathbf{Y}^{(1)} \\ \mathbf{Y}^{(2)} \end{bmatrix} = \begin{bmatrix} [Y_1,Y_2,Y_3]^{(1)} \\ [Y_1,Y_2,Y_3]^{(2)} \end{bmatrix} $$ Then, $$ (\mathbf{X} \otimes \mathbf{X})^{(1)} = (\mathbf{X}^{(1)})^T \cdot \mathbf{X}^{(1)} = \begin{bmatrix} 0.024891 & 0.046281 & 0.030511\\ 0.046281 & 0.086056 & 0.056731\\ 0.030511 & 0.056731 & 0.037399\end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(2)} = (\mathbf{X}^{(2)})^T \cdot \mathbf{X}^{(2)} = \begin{bmatrix} 0.017446 & 0.013614 & 0.019642\\ 0.013614 & 0.010624 & 0.015328 \\ 0.019642 & 0.015328 & 0.022114 \end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(3)} = (\mathbf{X}^{(3)})^T \cdot \mathbf{X}^{(3)} = \begin{bmatrix} 0.084013 & 0.114904 & 0.099157\\ 0.114904 & 0.157153 & 0.135616 \\ 0.099157 & 0.135616 & 0.117030 \end{bmatrix} $$ Adding these and dividing by $n_x - 1 = 2$ gives $$ \mathbf{S}(\mathbf{X},\mathbf{X}) = \begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} $$ Note that $$ \mathbf{X}^T \cdot \mathbf{X} = (\mathbf{X} \otimes \mathbf{X})^{(1)} + (\mathbf{X} \otimes \mathbf{X})^{(2)} + (\mathbf{X} \otimes \mathbf{X})^{(3)} \,. $$ Similarly, $$ \mathbf{S}(\mathbf{Y},\mathbf{Y}) = \begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$

Pooled covariance example

Then the pooled covariance is $$ \mathbf{S} = \frac{3}{5}\begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} + \frac{2}{5}\begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$ Therefore $$ \mathbf{S} = \begin{bmatrix} 0.067531 & 0.020286 & 0.031615 \\ 0.020286 & 0.111047 & 0.076604 \\ 0.031615 & 0.076604 & 0.058824 \end{bmatrix} $$

Inverse of the pooled covariance matrix

Now compute the inverse of $\mathbf{S}$ which clearly has to be a square matrix. This is the most computationally intensive part of the calculation and requires special attention for high dimensional data. $$ \mathbf{S}^{-1} = \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix} $$

Mahalanobis distance

Now to find the distance (squared) between $\mathbf{X}^{(i)}$ and $\mathbf{Y}^{(j)}$ we use the relation $$ d^2_{ij} = (\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})\cdot\mathbf{S}^{-1} \cdot(\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})^T $$ First compute $\mathbf{Z} := \mathbf{X}^{(i)} - \mathbf{Y}^{(j)}$. Both vectors have to be the same length for this operation to work. Define $\mathbf{T} := \mathbf{S}^{-1}$. Then $$ d^2_{ij} = \mathbf{Z} \cdot \mathbf{T} \cdot \mathbf{Z}^T $$ In terms of indices, with $N$ as the length of the vector (number of variables), $$ d^2_{ij} = \sum_{i=1}^N \sum_{j=1}^N Z_i T_{ij} Z_j $$ That should give you the distances you seek.

Mahalanobis distance example

Let us find the distance between the vectors $$ \mathbf{X}^{(2)} = \begin{bmatrix} 0.13208 & 0.10307 & 0.14871 \end{bmatrix} $$ and $$ \mathbf{Y}^{(1)} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594 \end{bmatrix} $$ We have, $$ \mathbf{Z} = \mathbf{X}^{(2)} - \mathbf{Y}^{(1)} = \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} $$ Therefore, $$ \begin{aligned} (d_{21})^2 &= \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix}\begin{bmatrix} 0.324521 \\ -0.105784 \\ 0.063114 \end{bmatrix} \\ &= 3.4722 \,. \end{aligned} $$

Samples with one observations each

Merge the two samples into one and just use the sample covariance instead of the pooled covariance. The result will not make much dense and a simpler distance measure is preferable.

There are several answers on similar questions on StackExchange (see for example Pairwise Mahalanobis distance in R).

Let us assume your matrices are "data matrices", i.e, they have the dimensions $n \times p$ where $n$ is the number of samples and $p$ is the number of variables.

Ideally you should center the data matrices by subtracting the column means from the columns before forming your covariance matrix and inverting it.

Then the cross-covariance matrix can be computed using $$ \begin{aligned} \mathbf{S} & = E\left[(\mathbf{X} - E[\mathbf{X}]) \otimes (\mathbf{Y} - E[\mathbf{Y}])\right] \\ & = E\left[\mathbf{X}\otimes\mathbf{Y} - E[\mathbf{X}]\otimes\mathbf{Y} - \mathbf{X}\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\right] \\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] \end{aligned} $$ In terms of matrix components $$ S_{ij} = E[X_i\,Y_j] - E[X_i]\,E[Y_j] $$ In your example, if $\mathbf{X}^{(1)} = [X_1,X_2,X_3]^{(1)}$ = [10,5,0] and $\mathbf{Y}^{(1)} = [Y_1,Y_2,Y_3]^{(1)}$ = [7,0,9], then $$ (\mathbf{X} \otimes \mathbf{Y})^{(1)} = \mathbf{X}^T \cdot \mathbf{Y} = \begin{bmatrix} 70 & 0 & 90 \\ 35 & 0 & 45 \\ 0 & 0 & 0\end{bmatrix} $$ You can compute this matrix for the other $N=n_x \times n_y$ combinations of samples from $\mathbf{X}$ and $\mathbf{Y}$.

If we assume that the expected value can be estimated by the sample mean, then we can use $$ E[\mathbf{X}] = \frac{1}{n_x}\sum_{i=1}^{n_x} \mathbf{X}^{(i)} ~,~~ E[\mathbf{Y}] = \frac{1}{n_y}\sum_{i=1}^{n_y} \mathbf{Y}^{(i)} $$ and $$ E[\mathbf{X} \otimes \mathbf{Y}] = \frac{1}{N-1}\sum_{i=1}^N (\mathbf{X} \otimes \mathbf{Y})^{(i)} $$ Plug these in to compute $\mathbf{S}$ and then compute its inverse. $\mathbf{S}$ clearly has to be a square matrix.

Now to find the distance (squared) between $\mathbf{X}^{(i)}$ and $\mathbf{Y}^{(j)}$ we use the relation $$ d^2_{ij} = (\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})^T\cdot\mathbf{S}^{-1} \cdot(\mathbf{X}^{(i)} - \mathbf{Y}^{(j)}) $$ First compute $\mathbf{Z} := \mathbf{X}^{(i)} - \mathbf{Y}^{(j)}$. Both vectors have to be the same length for this operation to work. Define $\mathbf{T} := \mathbf{S}^{-1}$. Then $$ d^2_{ij} = \mathbf{Z}^T \cdot \mathbf{T} \cdot \mathbf{Z} $$ In terms of indices, with $p$ as the length of the vector (number of variables), $$ d^2_{ij} = \sum_{i=1}^p \sum_{j=1}^p Z_i T_{ij} Z_j $$ That should give you the distances you seek.

There are several answers on similar questions on StackExchange (see for example Pairwise Mahalanobis distance in R).

Centering the data

Let us assume your matrices are "data matrices", i.e, they have the dimensions $n \times N$ where $n$ is the number of observations in a sample and $N$ is the number of variables.

Ideally you should center the data matrices by subtracting the column means from the columns before forming your covariance matrix and inverting it. Let us assume that, in your example, two samples are $$ \mathbf{X} = \begin{bmatrix} 0.678177 & 0.989365 & 0.558944\\ 0.652491 & 0.799086 & 0.514262\\ 0.230560 & 0.299587 & 0.023458 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} 0.015513 & 0.932740 & 0.351038 \\ 0.400391 & 0.515025 & 0.179851 \end{bmatrix} $$ Let us center the data by subtracting the sample means which are $$ E[\mathbf{X}] = \begin{bmatrix} 0.52041 & 0.69601 & 0.36555 \end{bmatrix} $$ and $$ E[\mathbf{Y}] = \begin{bmatrix} 0.20795 & 0.72388 & 0.26544\end{bmatrix} $$ The centered data are $$ \mathbf{X} = \begin{bmatrix} 0.15777 & 0.29335 & 0.19339\\ 0.13208 & 0.10307 & 0.14871\\ -0.28985 & -0.39643 & -0.34210 \end{bmatrix} $$ and $$ \mathbf{Y} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594\\ 0.192439 & -0.208857 & -0.085594 \end{bmatrix} $$

Covariance

Then the cross-covariance matrix can be computed using $$ \begin{aligned} \mathbf{S}(\mathbf{X},\mathbf{Y}) & = E\left[(\mathbf{X} - E[\mathbf{X}]) \otimes (\mathbf{Y} - E[\mathbf{Y}])\right] \\ & = E\left[\mathbf{X}\otimes\mathbf{Y} - E[\mathbf{X}]\otimes\mathbf{Y} - \mathbf{X}\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\right] \\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] + E[\mathbf{X}]\otimes E[\mathbf{Y}]\\ & = E[\mathbf{X}\otimes\mathbf{Y}] - E[\mathbf{X}]\otimes E[\mathbf{Y}] \end{aligned} $$ In terms of matrix components $$ S_{ij} = E[X_i\,Y_j] - E[X_i]\,E[Y_j] $$ If we assume that the expected value can be estimated by the sample mean, then we can use $$ E[\mathbf{X}] = \frac{1}{n_x}\sum_{i=1}^{n_x} \mathbf{X}^{(i)} ~,~~ E[\mathbf{Y}] = \frac{1}{n_y}\sum_{i=1}^{n_y} \mathbf{Y}^{(i)} $$ and, with $N$ being the number of variables (columns of the two matrices), $$ E[\mathbf{X} \otimes \mathbf{Y}] = \frac{1}{N-1}\sum_{i=1}^N (\mathbf{X} \otimes \mathbf{Y})^{(i)} $$ If the matrices have been centered, $E[X_i] = E[Y_j] = 0$ and we have $$ S_{ij} = E[X_i\,Y_j] \quad \implies \mathbf{S}(\mathbf{X},\mathbf{Y}) = E[\mathbf{X} \otimes \mathbf{Y}]\,. $$

Pooled covariance

The correlations among the data in each sample are estimated using the sample covariance and an estimate of the population covariance is computed using the pooled covariance.

The pooled covariance is given by $$ \mathbf{S} = \frac{n_x}{n_x+n_y} \mathbf{S}(\mathbf{X},\mathbf{X}) + \frac{n_y}{n_x+n_y} \mathbf{S}(\mathbf{Y},\mathbf{Y}) $$

Covariance example

In the example, if $$ \mathbf{X} = \begin{bmatrix} \mathbf{X}^{(1)} \\ \mathbf{X}^{(2)} \\ \mathbf{X}^{(3)} \end{bmatrix} = \begin{bmatrix} [X_1,X_2,X_3]^{(1)} \\ [X_1,X_2,X_3]^{(2)} \\ [X_1,X_2,X_3]^{(3)} \end{bmatrix} \quad \text{and} \quad \mathbf{Y} = \begin{bmatrix} \mathbf{Y}^{(1)} \\ \mathbf{Y}^{(2)} \end{bmatrix} = \begin{bmatrix} [Y_1,Y_2,Y_3]^{(1)} \\ [Y_1,Y_2,Y_3]^{(2)} \end{bmatrix} $$ Then, $$ (\mathbf{X} \otimes \mathbf{X})^{(1)} = (\mathbf{X}^{(1)})^T \cdot \mathbf{X}^{(1)} = \begin{bmatrix} 0.024891 & 0.046281 & 0.030511\\ 0.046281 & 0.086056 & 0.056731\\ 0.030511 & 0.056731 & 0.037399\end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(2)} = (\mathbf{X}^{(2)})^T \cdot \mathbf{X}^{(2)} = \begin{bmatrix} 0.017446 & 0.013614 & 0.019642\\ 0.013614 & 0.010624 & 0.015328 \\ 0.019642 & 0.015328 & 0.022114 \end{bmatrix} $$ $$ (\mathbf{X} \otimes \mathbf{X})^{(3)} = (\mathbf{X}^{(3)})^T \cdot \mathbf{X}^{(3)} = \begin{bmatrix} 0.084013 & 0.114904 & 0.099157\\ 0.114904 & 0.157153 & 0.135616 \\ 0.099157 & 0.135616 & 0.117030 \end{bmatrix} $$ Adding these and dividing by $N - 1 = 2$ gives $$ \mathbf{S}(\mathbf{X},\mathbf{X}) = \begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} $$ Note that $$ \mathbf{X}^T \cdot \mathbf{X} = (\mathbf{X} \otimes \mathbf{X})^{(1)} + (\mathbf{X} \otimes \mathbf{X})^{(2)} + (\mathbf{X} \otimes \mathbf{X})^{(3)} \,. $$ Similarly, $$ \mathbf{S}(\mathbf{Y},\mathbf{Y}) = \begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$

Pooled covariance example

Then the pooled covariance is $$ \mathbf{S} = \frac{3}{5}\begin{bmatrix} 0.063174 & 0.087400 & 0.074654 \\ 0.087400 & 0.126916 & 0.103837 \\ 0.074654 & 0.103837 & 0.088272 \end{bmatrix} + \frac{2}{5}\begin{bmatrix} 0.074066 &-0.080385 &-0.032943 \\ -0.080385 & 0.087243 & 0.035754 \\ -0.032943 & 0.035754 & 0.014653 \end{bmatrix} $$ Therefore $$ \mathbf{S} = \begin{bmatrix} 0.067531 & 0.020286 & 0.031615 \\ 0.020286 & 0.111047 & 0.076604 \\ 0.031615 & 0.076604 & 0.058824 \end{bmatrix} $$

Inverse of the pooled covariance matrix

Now compute the inverse of $\mathbf{S}$ which clearly has to be a square matrix. This is the most computationally intensive part of the calculation and requires special attention for high dimensional data. $$ \mathbf{S}^{-1} = \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix} $$

Mahalanobis distance

Now to find the distance (squared) between $\mathbf{X}^{(i)}$ and $\mathbf{Y}^{(j)}$ we use the relation $$ d^2_{ij} = (\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})\cdot\mathbf{S}^{-1} \cdot(\mathbf{X}^{(i)} - \mathbf{Y}^{(j)})^T $$ First compute $\mathbf{Z} := \mathbf{X}^{(i)} - \mathbf{Y}^{(j)}$. Both vectors have to be the same length for this operation to work. Define $\mathbf{T} := \mathbf{S}^{-1}$. Then $$ d^2_{ij} = \mathbf{Z} \cdot \mathbf{T} \cdot \mathbf{Z}^T $$ In terms of indices, with $N$ as the length of the vector (number of variables), $$ d^2_{ij} = \sum_{i=1}^N \sum_{j=1}^N Z_i T_{ij} Z_j $$ That should give you the distances you seek.

Mahalanobis distance example

Let us find the distance between the vectors $$ \mathbf{X}^{(2)} = \begin{bmatrix} 0.13208 & 0.10307 & 0.14871 \end{bmatrix} $$ and $$ \mathbf{Y}^{(1)} = \begin{bmatrix} -0.192439 & 0.208857 & 0.085594 \end{bmatrix} $$ We have, $$ \mathbf{Z} = \mathbf{X}^{(2)} - \mathbf{Y}^{(1)} = \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} $$ Therefore, $$ \begin{aligned} (d_{21})^2 &= \begin{bmatrix} 0.324521 &-0.105784 & 0.063114 \end{bmatrix} \begin{bmatrix} 84.030 & 155.460 & -247.610 \\ 155.460 & 376.188 & -573.445 \\ -247.610 & -573.445 & 896.850 \end{bmatrix}\begin{bmatrix} 0.324521 \\ -0.105784 \\ 0.063114 \end{bmatrix} \\ &= 3.4722 \,. \end{aligned} $$

Samples with one observations each

Merge the two samples into one and just use the sample covariance instead of the pooled covariance. The result will not make much dense and a simpler distance measure is preferable.