3
$\begingroup$

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.

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

This is going to be an IO-bound operation, so there is little practical advantage in looking for alternative formulations. There are two things you can do to speed that up:

  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.
  2. compile that loop. It doesn't matter much what language you use, Cython, Numba, Tensorflow (using a GPU may matter, of course, under the right conditions), C, Julia, whatever. Just not plain-old interpreted, weakly-typed Python.

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.