SymmetricTensors.jl

SymmetricTensors.jl provides the SymmetricTensors{T, N} type used to store fully symmetric tensors in more efficient way, without most of redundant repetitions. It uses blocks of Array{T, N} stored in Union{Array{Float,N}, Nothing} structure. Repeated blocks are replaced by Void. The module introduces SymmetricTensors{T, N} type and some basic operations on this type. As of 01/01/2017 @kdomino is the lead maintainer of this package.

Within Julia, just use run

pkg> add SymmetricTensors

to install the files. Julia 1.0 or later is required.

julia> data = ones(4,4);


julia> SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0]; nothing [1.0 1.0; 1.0 1.0]], 2, 2, 4, true)

From Array{T, N} to SymmetricTensors{T, N}

julia> SymmetricTensors(data::Array{T, N}, bls::Int = 2)

where bls is the size of a block. It is a parameter affecting the compuational speed of cumulants. The block size must fulfill bls ∈ {1, 2,..., dats} where dats = size(data, 1) = ... = size(data, N) otherwise error is risen.

julia> data = ones(4,4);


julia> convert(SymmetricTensor, data, 2)
SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[1.0 1.0; 1.0 1.0] [1.0 1.0; 1.0 1.0]; nothing [1.0 1.0; 1.0 1.0]], 2, 2, 4, true)

From SymmetricTensors{T, N} to Array{T, N}

julia> Array(st::SymmetricTensors{T, N})

Wrong block size:

julia> SymmetricTensor(ones(4,4), 5)
ERROR: DimensionMismatch("bad block size 5 > 4")

• frame::ArrayNArrays{T,N}: stores data, where ArrayNArrays{T,N} = Array{Union{Array{T, N}, Nothing}, N}
• bls::Int: the size of ordinary block (the same in each direction),
• bln::Int: maximal number of blocks in each direction,
• dats::Int: the size of data stored (the same in each direction),
• sqr::Bool: if all blocks are squares (N-squares).

Suppose we have N = 2 and dats = 6 and bls = 3 hence data are symmetric matrix of size 6 x 6. Data are stored in the form:

|B11   B12 | 
|null  B22 | 

here bls = 2 and size of B11, B12, and B22 are 3 x 3. Bear in mind, that B11 and B22 his to be symmetric. As B12 (the last block) is square, sqr = True.

Suppose now N = 2 and dats = 5 and bls = 3 hence data are symmetric matrix of size 5 x 5. Data are stored in similar form:

|B11   B12 | 
|null  B22 | 

here bls = 2 and size of B11 is 3 x 3, but size of B12 is 2 x 3, and size B22 is 2 x 2 . Again B11 and B22 his to be symmetric. As B12 (the last block) is not square, sqr = False.

For N = 3 we have analogical pyramid representation, and for N > 3 hyper-pyramid representation.

Elementwise addition: +, - is supported between many SymmetricTensors{T, N} while elementwise substraction: - between two SymmetricTensors{T, N}. Addition substraction multiplication and division +, -, *, / is supported between SymmetricTensors{T, N} and a number.

julia> x = SymmetricTensor(ones(4,4));

julia> y = SymmetricTensor(2*ones(4,4));

julia> x+y
SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[3.0 3.0; 3.0 3.0] [3.0 3.0; 3.0 3.0]; #undef [3.0 3.0; 3.0 3.0]], 2, 2, 4, true)

julia> x*10
SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[10.0 10.0; 10.0 10.0] [10.0 10.0; 10.0 10.0]; #undef [10.0 10.0; 10.0 10.0]], 2, 2, 4, true)

The function diag returns a Vector{T}, of all super-diagonal elements of a SymmetricTensor.

julia> data = ones(5,5,5,5);

julia> st = SymmetricTensor(data);

julia> diag(st)
5-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.0

To generate random Symmetric Tensor with random elements of typer T form a uniform distribution on [0,1) use rand(SymmetricTensor{T, N}, n::Int, b::Int = 2). Here n denotes size of each mode and b denotes block size. Eg. for N = 4 we would have n x n x n x n tensor.

julia> using Random

julia> Random.seed!(42)

julia> x = rand(SymmetricTensor{Float64, 3}, 2)
SymmetricTensor{Float64, 3}(Union{Nothing, Array{Float64, 3}}[[0.5331830160438613 0.4540291355871424; 0.4540291355871424 0.017686826714964354]

[0.4540291355871424 0.017686826714964354; 0.017686826714964354 0.17293302893695128]], 2, 1, 2, true)

julia> Array(x)
2×2×2 Array{Float64, 3}:
[:, :, 1] =
 0.533183  0.454029
 0.454029  0.0176868

[:, :, 2] =
 0.454029   0.0176868
 0.0176868  0.172933

julia> Random.seed!(42)

julia> x = rand(SymmetricTensor{Float64, 2}, 3)
SymmetricTensor{Float64, 2}(Union{Nothing, Matrix{Float64}}[[0.5331830160438613 0.4540291355871424; 0.4540291355871424 0.017686826714964354] [0.17293302893695128; 0.9589258763297348]; nothing [0.9735659798036858]], 2, 2, 3, false)

julia> Array(x)
3×3 Matrix{Float64}:
 0.533183  0.454029   0.172933
 0.454029  0.0176868  0.958926
 0.172933  0.958926   0.973566

julia> using Random

julia> Random.seed!(42)

julia> st = rand(SymmetricTensor{Float64, 2}, 2)
SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[0.533183 0.454029; 0.454029 0.0176868]], 2, 1, 2, true)

julia> st[1,2]
0.4540291355871424

julia> st[2,1]
0.4540291355871424


julia> pyramidindices(st)
3-element Vector{Tuple{Int64, Int64}}:
 (1, 1)
 (1, 2)
 (2, 2)

Function pyramidindices(st::SymmetricTensor) returns the indices of the unique element of the give symmetric tensor

setindex!(st::SymmetricTensor, x::Float, mulind::Int...) changes all symmetric tensor's elements indexed by mulind to x.

julia> st[1,2] = 10.

julia> convert(Array, st)
2×2 Array{Float64,2}:
  0.533183  10.0      
 10.0        0.0176868

julia> unfold(data::Array{T,N}, mode::Int)

unfolds data in a given mode

julia> a = reshape(collect(1.:8.), (2,2,2))
2×2×2 Array{Float64,3}:
[:, :, 1] =
 1.0  3.0
 2.0  4.0

[:, :, 2] =
 5.0  7.0
 6.0  8.0

julia> unfold(a, 1)
2×4 Array{Float64,2}:
 1.0  3.0  5.0  7.0
 2.0  4.0  6.0  8.0

julia> unfold(a, 2)
2×4 Array{Float64,2}:
 1.0  2.0  5.0  6.0
 3.0  4.0  7.0  8.0

julia> unfold(a, 3)
2×4 Array{Float64,2}:
 1.0  2.0  3.0  4.0
 5.0  6.0  7.0  8.0

The block usage is motivated by the paper M. D. Schatz, T. M. Low, R. A. van de Geijn, and T. G. Kolda, "Exploiting symmetry in tensors for high performance: Multiplication with symmetric tensors", SIAM Journal on Scientific Computing, 36 (2014), pp. C453–C479 https://doi.org/10.1137/130907215. There only the meaningful part of the symmetric tensor is stored in blocks to decrease the memory and computational overhead.

For application of this representation to compute cumulants, see: K. Domino, P. Gawron, Ł. Pawela "Efficient Computation of Higher-Order Cumulant Tensors", SIAM Journal on Scientific Computing, 40 (2018) https://doi.org/10.1137/17M1149365. The selection of the optimal block size is not straight forward, however in most cases concerning cumulants one can use 2 or 3.

This project was partially financed by the National Science Centre, Poland – project number 2014/15/B/ST6/05204.