I am writing a program which has a tensor contraction with a 4d symmetrical tensor at it’s core. Needless to say, since it will be used inside a non linear problem solver in the future, I need this tensor contraction to evaluate as fast as possible.
My first attempt was to use the TensorOperations.jl package for which the implementation would look something like this:
using TensorOperations
using BenchmarkTools
n = 10;
@time A = ones(n, n, n, n);
@time B = rand(n, n);
@time u = rand(n);
@time v = rand(n);
@btime @tensor B[a, b] := A[a, b, c, d] * u[c] * v[d];
This runs pretty fast all together:
Output
0.000040 seconds (4 allocations: 78.281 KiB)
0.000010 seconds (1 allocation: 896 bytes)
0.000009 seconds (1 allocation: 144 bytes)
0.000012 seconds (1 allocation: 144 bytes)
21.500 μs (45 allocations: 90.20 KiB)
My second attempt was to use the Tensorial.jl package which allows me to define a static symmtrical tensor to save space and duplications which looks something like this:
using Tensorial
using BenchmarkTools
n = 10;
@time A = Tensor{Tuple{@Symmetry{n,n,n,n}}}(ones(n, n, n, n))
@time u = Tensor{Tuple{n}}(rand(n))
@time v = Tensor{Tuple{n}}(rand(n))
@btime A ⋅ u ⋅ v;
Here the contraction itself is very fast but assigning the Tensor object is very very slow:
Output
131.126936 seconds (3.13 M allocations: 6.134 GiB, 0.31% gc time, 100.00% compilation time)
0.398449 seconds (273.82 k allocations: 18.272 MiB, 2.93% gc time, 99.91% compilation time)
0.000031 seconds (12 allocations: 720 bytes)
369.903 ns (2 allocations: 2.20 KiB)
Here is my original Matlab code from when I originally wrote the project in matlab before rewtriting in Julia:
n = 100;
tic; A = ones(n, n, n, n); toc
tic; u = rand(n); toc
tic; v = rand(n); toc
tic; tensorprod(tensorprod(A, u, 4, 1), v, 3, 1); toc
This seems to already run slower than the TensorOperations.jl Package so at least using Julia will be an improvement either way but it is not the massive advantage I hoped for. (Note for larger n the Matlab code is only slightly slower than the TensorOperations.jl package.)
Output
Elapsed time is 0.000372 seconds.
Elapsed time is 0.000252 seconds.
Elapsed time is 0.000069 seconds.
Elapsed time is 0.000782 seconds.
Now, I would really like to keep the very fast contraction time of the Tensorial.jl package and I really only need to assign this big tensor once and then save it so that would make this the natural choice. Sadly I have to do these calculations for n>100 and I tried to construct this big symmetrical Tensor object but even after I left it running for an hour it did not do anything.
Is there a method I am missing to construct this symmetrical Tensor object faster or am I missing another obvious way to speed up this contraction? The Tensor.jl package does not seem to allow 4d tensors with a dimension bigger than 3 and I can’t seem to do tensor operations on static arrays from the StaticArrays.jl package.
Thanks in advance,
Laura
