# Best current approach for working with sparse tensors?

**URL:** <https://discourse.julialang.org/t/best-current-approach-for-working-with-sparse-tensors/79937>\
**Category:** Specific Domains\
**Tags:** sparse, tensors\
**Created:** [April 24, 2022, 11:26am UTC](https://discourse.julialang.org/t/best-current-approach-for-working-with-sparse-tensors/79937 "2022-04-24T11:26:13Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![xor0110](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xor0110/32/7926_2.png) [@xor0110](https://discourse.julialang.org/u/xor0110)\
**Post date:** [April 24, 2022, 11:26am UTC](https://discourse.julialang.org/t/best-current-approach-for-working-with-sparse-tensors/79937/1 "2022-04-24T11:26:13Z")

</div>

I would like to hear ideas about how to work with this kind of application right now in Julia. Some previous posts seem to be in line with what I’m looking for, but I’m not sure what’s the conclusion, and what’s the state of the art.

> [@Tensor multiplication](https://discourse.julialang.org/t/tensor-multiplication/21450):
>
> Hi all, I am implementing a function to perform a generalization of matrix multiplication to a general N-dimensional array or tensor. This product is denoted as \times\_m to multiply a conformable matrix A with a tensor \mathcal{X} according to dimension n. A working example is given below (note, I already tried several things to make it more performant: @inbounds, efficient looping over array elements etc). function tensormult4!(X, A, Y, T, dim) @assert size(X, dim) == size(A, 2) Tmax…

> [@Sparse high dimensional array support einsum?](https://discourse.julialang.org/t/sparse-high-dimensional-array-support-einsum/47222/6):
>
> Do you mean that you sum over all indices except the first and last? i there looks like it runs from 1 to 10, not 100. (But something else could label which indices don’t get summed, the nth & mth.) You could write this, but there’s little to gain over writing the loops yourself: julia\> inds = rand(1:10, 333, 100); vals = randn(333); julia\> out = zeros(10,10); julia\> using Tullio julia\> @tullio out[inds[r,1], inds[r,100]] += vals[r]

As an example, suppose I’m simply computing the gradient from an image. We can represent the image `img[j,k]` as a vector `img[:]` and each gradient as a sparse matrix `Dx` and `Dy`. It’s useful to keep all of this data in the same array for further linear operations, so we can append these matrices into a `D = [Dx; Dy]`. But now if I want to access my gradient components I need to deal with computing the indices, or reshape my vector.

Ideally, I would like to be able to have my operation as a tensor, and using Einstein summation I can compute `dimg[a,p] = D[a,p,q] * img[q]`. This gets more exciting later if we also have separate color channels. Perhaps we might even preserve the original image shape as well?

In my specific application I don’t have a huge number of dimensions. I could write my code handling indexing by myself if necessary, but I want to be able to rely on Julia for indexing as much as possible.

I’m also not super worried about performance optimization, apart from exploiting the matrix sparsity, and hopefully Julia can optimize some matrix multiplications. In other words, I would like to avoid writing lots of for-loops, but I’m not expecting an amazingly optimized code à la [http://tensor-compiler.org/](http://tensor-compiler.org/).

The issues I’ve run into right now was that I tried to represent this gradient operator as a reshaped  
`SparseMatrixCSC`, but both `Einsum` and `TensorOperations` seemed to have problems with that. And I’m not sure ITensors supports the kind of sparsity I need.

I hope this doesn’t sound like I’m criticizing the current tools because I have some super specific demand that is not fulfilled yet! I’m just thinking other developers might have dealt with similar problems in the past. What’s a good way to deal with this problem right now? For example, is there a different way I could be dealing with sparsity that might help me out? Or is there a handy way I can reshape arrays in my code whenever needed?

---

<div class="post-metadata">

**Author:** ![xor0110](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xor0110/32/7926_2.png) [@xor0110](https://discourse.julialang.org/u/xor0110)\
**Post date:** [April 24, 2022, 11:40am UTC](https://discourse.julialang.org/t/best-current-approach-for-working-with-sparse-tensors/79937/2 "2022-04-24T11:40:21Z")

</div>

Some code for illustrative purposes:

```julia
using SparseArrays
jj,kk=4,4
N=jj*kk
Dx = sparse(1:N, 1:N, ones(N)) + sparse(1:N, [N;1:N-1], -1)
Dy = sparse(1:N, 1:N, ones(N)) + sparse(1:N, [N-jj+1:N;1:N-jj], -1) 
D = reshape([Dx; Dy], N,N,2)
img = rand(N)
dimg = zeros(N,2)

## That works
dimg = hcat(D[:,:,1]*img, D[:,:,2]*img)

## results in `ERROR: MethodError: no method matching Strided.UnsafeStridedView(::SparseMatrixCSC{Float64, Int64})`
using TensorOperations
@tensor begin
    dimg[p,g] = D[p,q,g] * img[q]
end
## Seems to have some performance issue for larger arrays (eg N=128).
using Einsum
@einsum dimg[p,g] = D[p,q,g] * img[q]

```

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 24, 2022, 3:11pm UTC](https://discourse.julialang.org/t/best-current-approach-for-working-with-sparse-tensors/79937/3 "2022-04-24T15:11:09Z")

</div>

I have seen [GitHub - ITensor/ITensors.jl: A Julia library for efficient tensor computations and tensor network calculations](https://github.com/ITensor/ITensors.jl) mentioned around here.