Move from `SparseArrays` to a package with multidimensional sparse array support?

[akirakyle]

There's currently a sparsematrix.jl file here and a sparse.jl file in QuantumInterface with functions extending SparseArrays with:

• some gemm! and gemv! functions
• _permutedims which works on arrays with ndims > 2 as SparseArrays only implements permutedims for ndims = 2
• ptrace

I wasn't able to find much evidence that SparseArrays has plans to support anything other than vectors and matrices. The docs even list some external packages that handle multidimensional arrays. Several of these seem to offer interfaces that could help simplify the code we have here. For example, anything offering an Einstein summation notation interface would make ptrace a one-liner.

Finch.jl is perhaps the most interesting one on the list for our use cases and potential to offer significant speedups. Also it might be possible to use its lazy feature as part of the Lazy{Ket,Product,Sum,Tensor}. However, it's still lacking some basic features we would want here like reshape.

A better library for multidimensional sparse arrays would be especially helpful for having cleaner code for the all the superoperator support I'm working on as it would hopefully avoid issues like this (relevant issue: JuliaSparse/SparseArrays.jl#24).

[willow-ahrens]

@akirakyle Thanks for your interest in Finch! Let me know if you have any questions about using Finch, I'd like to help!

[akirakyle]

@willow-ahrens Thank you! Hopefully I'll have some time to experiment with Finch here and I'll let you know how it goes

[akirakyle]

Also, Finch might be useful for the generalized stabilizer tableu being developed in QuantumCliffords. Specifically, I was looking at Ted Yoder's paper in more detail today, and the following sentence stood out as describing the sort of optimization that, if I understand correctly, the Galley optimizer in Finch can perform:

Multiplying [the superoperator matrices that are applied to χ] with one another before application to χ might very well be faster, if χ is dense, than multiplying χ by each in turn. Collecting all such superoperators at the end of the circuit, then determining the optimal multiplication order, makes such a simplification possible.

@Krastanov @Fe-r-oz