[ANN] PermutationSymmetricTensors.jl: Multidimensional arrays that are symmetric under index permutation

Package Announcements
1

I am happy to announce my first package: PermutationSymmetricTensors.jl.

Feel free to try it out and let me know what you think. Suggestions, issues, feature and/or pull requests are very welcome. Note that, although I have written many unit tests, the package has not really been battle-tested yet. Below is a summary of the documentation:

PermutationSymmetricTensors provides an efficient framework for the use of multidimensional arrays that are symmetric under any permutation of their indices, implemented in pure Julia. Such symmetric tensors are implemented in the SymmetricTensor{T, N, dim} type, where T is the element type, dim is the number of indices required to index into the tensor, and N is maximal index for each dimension. For example, to index a SymmetricTensor{ComplexF64, 20, 6} , you need 6 indices between 1 and 20.

This package exports basic constructors of SymmetricTensors, and a few convenience functions for working with them. The main advantage of using a SymmetricTensor is that it requires much less memory to store than the full Array would.

Usage

A SymmetricTensor is a mutable subtype of AbstractArray, and thus supports most operations that we are used to (see the documentation for notable exceptions).

Example:

julia> using PermutationSymmetricTensors

julia> a = rand(SymmetricTensor{Int8, 3, 3})
3×3×3 SymmetricTensor{Int8, 3, 3}:
[:, :, 1] =
 -115  -31  117
  -31  110   95
  117   95  -57

[:, :, 2] =
 -31  110    95
 110  -30    33
  95   33  -106

[:, :, 3] =
 117    95   -57
  95    33  -106
 -57  -106    87

julia> a[1,2,3]
95

julia> a[3,1,2]
95

julia> a[3,2,2] = 6
6

julia> a
3×3×3 SymmetricTensor{Int8, 3, 3}:
[:, :, 1] =
 -115  -31  117
  -31  110   95
  117   95  -57

[:, :, 2] =
 -31  110    95
 110  -30     6
  95    6  -106

[:, :, 3] =
 117    95   -57
  95     6  -106
 -57  -106    87

julia> a[:, 1, 1] .= 0
3-element view(::SymmetricTensor{Int8, 3, 3}, :, 1, 1) with eltype Int8:
 0
 0
 0

julia> a
3×3×3 SymmetricTensor{Int8, 3, 3}:
[:, :, 1] =
 0    0    0
 0  110   95
 0   95  -57

[:, :, 2] =
   0  110    95
 110  -30     6
  95    6  -106

[:, :, 3] =
   0    95   -57
  95     6  -106
 -57  -106    87

As, you can see, mutating one element also mutates all elements that are equivalent under index permutation. SymmetricTensors leverage symmetry to minimize memory usage. It is easy to create SymmetricTensors that would have more elements than typemax(Int64) if stored naively:

julia> @time d = rand(SymmetricTensor{Float64, 14, 17});
  0.793732 seconds (23 allocations: 913.699 MiB, 1.09% gc time)

julia> println("This tensor requires ", round(sizeof(d)/2^30, digits=2), "GB memory")
This tensor requires 0.89GB memory

julia> println("a full array of this shape would require ", length(d)*8/2^30, "GB memory.")
a full array of this shape would require 2.27e11GB memory.

julia> println("it would have ", 14^17, " elements.")
it would have -6402141418087907328 elements.

julia> println("oops, I meant ", big(14)^17)
oops, I meant 30491346729331195904

In the second to last line, the computation 14^17 overflowed, and therefore returned the wrong result.

Please read the documentation for a more complete account of the implemented functionality.

See also

There are two packages with comparable functionality, SymmetricTensors.jl and Tensors.jl.

Tensors.jl provides immutable, stack-allocated 1-, 2-, and 4-dimensional symmetric tensors. This package is preferable if the tensors are small, that is, when they have fewer than roughly 100 elements. It is also more full-featured, implementing many different Linear Algebra operations on the tensors instead of just the basic functionality.

SymmetricTensors.jl provides a SymmetricTensor type just like the one exported in this package. Its implementation is based on a blocked memory pattern, sacrificing performance for cache locality. Additionally, SymmetricTensors.jl is less general since it only allows elements with supertype AbstractFloat, while PermutationSymmetricTensors.jl allows arbitrary element types.

A quick benchmark:

import SymmetricTensors
using BenchmarkTools

Tensor creation:

julia> T = Float64;
julia> N = 30;
julia> dim = 5;

# SymmetricTensors.jl
julia> a  = @btime rand(SymmetricTensors.SymmetricTensor{$T, $dim}, $N);
  2.601 s (21697571 allocations: 1.00 GiB)

julia> sizeof(a.frame)
6075000

# This Package
julia> b  = @btime rand(PermutationSymmetricTensors.SymmetricTensor{$T, $N, $dim});
  642.400 μs (10 allocations: 2.12 MiB)

julia> sizeof(b)
2227288

getindex and setindex!:

# SymmetricTensors.jl
julia> @btime $a[1,5,2,5,21];
  2.189 μs (24 allocations: 1.38 KiB)

julia> @btime $a[1,5,2,5,21] = 6.0;
  50.800 μs (767 allocations: 51.06 KiB)

# This Package
julia> @btime $b[1,5,2,5,21];
  4.200 ns (0 allocations: 0 bytes)

julia> @btime $b[1,5,2,5,21] = 6.0;
  4.500 ns (0 allocations: 0 bytes)

I decided to publish a new package rather than file PRs to SymmetricTensors.jl since this is a completely different implementation with a different memory layout (and because I wanted to learn how to publish a package). There might well be situations where SymmetricTensors.jl is more performant or convenient, but since I don’t fully understand their code, I find that hard to judge.

Please try it out and let me know what you think! Since I started using Julia relatively recently, I might be making some weird mistakes. Any feedback on that front is also welcome.

4 Likes