# Compute distances from a vector to a matrix of vectors

Let $$\vec{a} \in \mathbb{R}^\alpha$$, and let $$H$$ be a rank 3 tensor with dimensions $$M_{[i \in \mathbb{N}]} \times N_{[j \in \mathbb{N}]} \times \alpha_{[k \in \mathbb{N}]}$$ (where the subscript are the indexing of each one of the dimensions).

We can analogously look at the tensor $$H$$ as a matrix $$H'$$ which each element is an $$\mathbb{R}^\alpha$$ vector.

I wish to compute a norm (lets say for sake of example the $$\| \cdot \|_2^2$$ norm, squared here just to avoid taking the root after the computation of differences) of vector $$\vec{a}$$ with each element of matrix $$H'$$

A simple algorithm (using python without numpy as example) for that can be expressed as:

D = [[0 for i in range(M)] for j in range(N)]
for i in range(M):
for j in range(N):
norm = 0
for k in range(alpha):
norm += (a[k] - H[i][j][k])**2
D[i][j] == norm


Mathematically we could write:

$$D_{i,j} = \| H_{i,j} - \vec{a} \|_2^2 = \sum_{k=1}^\alpha \left(H_{{i,j}_k} - a_k \right)^2$$

My question are:

1. There is any mathematical operation that express this more clearly? Any matrix/tensor operation that can express that more concisely?

2. Is there better algorithms, approaches to this computation? Since I will run large tensors with large vectors I will calculate the distance to many input vectores. My idea is to optimize using gpus, probably with tensorflow.

• Is it important that your vectors are arrange over two indices $i,j$? Couldn't you just unroll these indices into one? – Wolfgang Bangerth Feb 26 '20 at 18:22

1. choose the order of indices and for loops to make sure that memory accesses are as contiguous as possible. In your case (assuming Fortran-style memory layout, which should be what Python does, and assuming a takes less space than D), I think you want to reshape H so that k is the first index, and then you want to swap the order of the i and j loops.
As an alternative to doing 2 by hand, look if you can find a library function that does it for you. For instance, a quick Google search returned scipy.spatial.distance.cdist, but I don't know if it is compiled and fast enough for your purposes.
Matlab users would tell you to use bsxfun or a subtraction with singleton expansion reshape(a, [1,1,length(a)]) - H followed by vecnorm, which probably is the fastest thing you can do in Matlab, but is still not optimal because it loops through $$H$$ twice.