Levi-Civita symbol/tensor

[Blake Johnson]

I keep writing code that needs the Levi-Civita symbol, i.e.:

\epsilon_{ijk} = { 1 if (i,j,k) is an even permutation, -1 if (i,j,k) is an odd permutation, otherwise 0 }

It is used frequently in physics, so I keep expecting it to be in Julia's stdlib, perhaps under a different name (in Mathematica, it is called Signature[]). But, if it is there, I haven't found it yet. Is it there and I just can't find it?

[Jiahao Chen]

> But, if it is there, I haven't found it yet

Me neither. I think Levi-Civita would be useful to have in Base if you'd like to write an implementation.

[Stefan Karpinski]

I agree that it's useful but what should it be called?

[Blake Johnson]

Some possibilities: levicivita, parity, sign, signature, or sgn.

"sign" is already an exported method, but sign(::Permutation) is not defined, yet.

[Jiahao Chen]

I like "sign(::Permutation)".

Thanks,

Jiahao Chen
Staff Research Scientist
MIT Computer Science and Artificial Intelligence Laboratory

[Dawid Crivelli]

I have a short code for checking the sign of a permutation for personal use. It's O(L^2) where L is the length of the permutation vector, but it's perfectly fine for L in the few tens range.

function sign{T<:Integer}(perm::AbstractVector{T})
    L = length(perm)
    crosses = 0
    for i = 1:L
        for j = i+1 : L
            crosses += perm[j] < perm[i]
        end
    end
    return iseven(crosses) ? 1 : -1    
end

[Stefan Karpinski]

Thanks for the code – would you be willing to contribute it under the MIT license?

[Dawid Crivelli]

Thanks for the code – would you be willing to contribute it under the MIT license?

Of course, I'm glad it can be of use

[Hans W Borchers]

The parity should be zero if the vector is not a true permutation, i.e.

    julia> sign([3, 4, 5, 4, 1])
    -1

would at least be misleading. Could result in an ERROR, but I still
suggest to return 0.

[Hans W Borchers]

Please also note that there is an `parity` function in the Permutations package that does exactly that:

    julia> using Permutations

    julia> p = Permutations([3,4,5,2,1])
    (1,3,5)(2,4)

    julia> parity(p)
    1

It returns 1 for odd permutations and 0 for even ones, and an error in case the vector is not a permutation of [1..n]. It computes the parity/sign by adding the length of the cycles.

[Pablo Zubieta]

Hi all!

I believe that

sign{T<:Integer}(perm::AbstractVector{T})

is already defined and returns something like

[sign(n) for n in perm]

So, I think it would be better to call it by other name, or with a different signature.

By looking at the code at the code used in Permutations, I came up with two solutions that should run in O(n) time. They are optimized versions of
Permutations.parity and can be found here:

https://gist.github.com/pabloferz/01675f1bf4c8be359767#file-levicivita-jl

The second option, is almost the same as the first one, but instead of first checking if the Vector is a permutation at the begining, it does the check in-place. There's probably not much difference between the two, I have not analyzed this deeper. Anyway, I believe both should work fine.

Cheers!

[Pablo Zubieta]

Hi again,

There were some bugs in my implementations. I updated the gist with the corrected versions and added a simpler looking function (but of O(n²) running time).

I did some tests and found (with my slow processor) that for permutations of length <= 5 the quadratic implementation (levicivita_simple) performs as fast as the (levicivita_inplace_check). For lengths from 5 to 15, levicivita_inplace_check is the fastest, followed by levicivita_simple. For lengths from 15 to 25 levicivita_simple and levicivita perform the same (but slower than levicivita_inplace_check). For more than 25 elements levicivita_inplace_check is always the fastest, 2x faster than levicivita and n times faster than levicivita_simple.

For people wanting the 3D Levi-Civita tensor, levicivita_simple and levicivita_inplace_check should be the same. For people wanting the parity of a permutation for long permutations levicivita_inplace_check should work the best.

Greetings!

[Blake Johnson]

Thanks for posting these, Pablo. For my most frequent use case I care about n = 3, but I suppose the O(n) algorithms would be more appropriate in Base.

You are also correct that sign(::AbstractVector) currently does an element-wise sign(). I didn't realize before writing my post that combinatorics.jl defines Permutations and Combinations types, but not the singular equivalents. So, that still leaves the naming issue unresolved.

Pablo, would you mind opening an issue or pull request to continue the discussion?

--Blake

[Pablo Zubieta]

Done!

https://github.com/JuliaLang/julia/issues/10172

[Andrew McLean]

For small n I would have expected a lookup-table (LUT) approach to be the best approach. Obviously there is a memory-cpu trade-off, but for n=3 there would only need to be 6 table entries! My suggestion would be to employ memoisation and build up the LUT on the fly. Also, I would suggest only caching the results for inputs that are actually permutations. I think a slower code path for inputs that do not represent permutations would be acceptable.

[lapeyre....@gmail.com]

Maybe its too late, but here is a fast implementation giving the signature of the output of randperm(n)

https://github.com/jlapeyre/PermPlain.jl/blob/master/src/PermPlain.jl#L271

[lapeyre....@gmail.com]

I should note that a couple of routines in PermPlain are not MIT license compatible.  But all the code to compute
the signature of a permutation most definitely is MIT compatible. I wrote it, although the algorithm is well known.

There is also this code which calls PermPlain:
https://github.com/jlapeyre/PermutationsA.jl