0

I want to contract large, n dimensional vectors with the Levi Civita tensor. If I want to use Numpy's einsum function, I have to define the Levi Civita tensor in advance, which quickly blows up my computer's memory at n dimensions. I would just have to get the einsum function to accept a callable so that I can do something like this: 

import operator
from functools import reduce

import numpy as np
from math import factorial

def prod(a):
    """
    Return product of elements of a. Start with int 1 so if only ints are included then an int result is
    returned.
    """
    return reduce(operator.mul, a, 1)


def eijk(args):
    """
    Represent the Levi-Civita symbol. For even permutations of indices it returns 1, for odd
    permutations -1, and for everything else (a repeated index) it returns 0. Thus it represents an
    alternating pseudotensor.

    Parameters
    ----------
    args : tuple with int
        A tuple with indices.

    Returns
    -------
    int
    """
    n = len(args)

    return prod(prod(args[j] - args[i] for j in range(i + 1, n)) / factorial(i) for i in range(n))


a = np.array([1., 0., 0., 8.])
b = np.array([1., 5., 8., 4.])

c = np.einsum("ijkl, j, k -> ij", eijk, a, b)

This will lead to ValueError: einstein sum subscripts string contains too many subscripts for operand 0. It is only working when the Levi Civita tensor is predefined like the following:

# The array `eijk_array` has a shape of (4, 4, 4, 4). That's why I won't type the entire tensor. 
eijk_array = array([[[[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0.,  1.],
                      [ 0.,  0., -1.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  0., -1.],
                      [ 0.,  0.,  0.,  0.],
                      [ 0.,  1.,  0.,  0.]],
                     [[ 0.,  0.,  0.,  0.],
                      [ 0.,  0.,  1.,  0.],
                      [ 0., -1.,  0.,  0.],
                      [ 0.,  0.,  0.,  0.]]],
                      ...

a = np.array([1., 0., 0., 8.])
b = np.array([1., 5., 8., 4.])

c = np.einsum("ijkl, j, k -> ij", eijk_array, a, b)

The handling of four dimensions is still going well. But when working with 10, 20 oder 30 dimensional vectors like in my case it can not be calculated in this form. For a three-dimensional vector I solved it as the following example shows:

a = np.array([1., 0., 0.])
b = np.array([1., 5., 8.])

args = [a, b]

# The dimension of the output. This is the amount of `.. -> ij` `np.einsum` subscriptions 
dim = a.shape[0] - len(args)

# Start Contraction (`eijk` is the predefined function above)
c = np.zeros(np.repeat(a.shape[0], dim))

for k in range(a.shape[0]):
    for i, j in product(range(a.shape[0]), repeat=2):
        c[k] = c[k] + a[i] * b[j] * eijk((i, j, k))

I don't know how to extend this loop to n dimensions, yet. Before I try to solve the problem with much effort, I wanted to ask if there is already a solution for this problem.

Edit

I have tried to trick Numpy by creating a subclass of Numpy.ndarray and modifying the magic method __getitem__ such like:

class LeviCivita(np.ndarray):
    def __new__(...):
        # ....

    def __getitem__(self, item):
        return eijk(item)

Unfortunately this did not work either.

Petermailpan
  • 103
  • 1
  • 10
  • 1
    What dimensions does the arrays you want to contract with usually have (always 1D) and how many at once? Even if it would be possible to pass a python levi-civita generating function to einsum it would be terribly slow (repeated call of python function). Also consider that the levi-civita tensor on 20 or 30 dims is huge and quite sparse. In general this looks to be a problem which can be solved with numba, cython, or any package specialized on this maybe written in (C, Fortran). – max9111 Apr 20 '20 at 17:53
  • 1
    Thank you for your valuable comment @max9111. I have written all functions you see above in Cython. For a better overview of the problem, I have converted the codes back to pure Python codes here. But you just gave me an idea. The contraction with the Levi-Civita tensor is in my case nothing else but calculating the determinants of the matrices, for which there are, as you already said accurately, many better and different solutions and algorithms. The dimensions of my array are (n, 2) or (n, 4). – Petermailpan Apr 20 '20 at 19:25

0 Answers0