Kronecker product does not preserve sparsity

[johnbcoughlin]

Passing a NamedDimsArray into the kron function does not preserve sparsity:

julia> using NamedDims, SparseArrays
julia> perm = sortperm(rand(20))
julia> p = NamedDimsArray{(:x, :y)}(permute(sparse(I(20)), 1:20, perm))
julia> K = kron(I(10), p)
julia> typeof(K)
Array{Bool,2}

Note that this is in contrast to passing an unwrapped sparse array:

julia> typeof(kron(I(10), parent(p)))
SparseMatrixCSC{Bool,Int64}

[mcabbott]

I think this package doesn't overload kron at all, but it does prevent dispatch to the sparse-sparse method, and so you get the default method:

julia> @which kron(p, p)
kron(a::AbstractMatrix{T}, b::AbstractMatrix{S}) where {T, S} in LinearAlgebra at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/dense.jl:416

julia> @which kron(p.data, p.data)
kron(A::SparseArrays.AbstractSparseMatrixCSC{T1, S1}, B::SparseArrays.AbstractSparseMatrixCSC{T2, S2}) where {T1, S1, T2, S2} in SparseArrays at /Users/me/.julia/dev/julia/usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:1385

The usual thing here would be to overload kron to unwrap named arguments, call it on those, and re-wrap. The question is what names the result should have, since each dimension combines two dimensions of the arguments. Should kron and vec behave the same here?

julia> nda = NamedDimsArray(rand(Int8,3,4), (:x, :y))
3×4 NamedDimsArray(::Matrix{Int8}, (:x, :y)):
     → y
↓ x  -69  107   52  -58
      24  125  -34  -80
      60   84    2   94

julia> xs = nda[y=1]; ys = nda[x=1];

julia> kron(xs, ys)
12-element Vector{Int8}:
 -103
   41
   -4
  -94
 ...
 
 julia> vec(nda)
12-element Vector{Int8}:
 -69
  24
  60
 107
...

[johnbcoughlin]

Yes, definitely for now I think the result should be an unnamed but sparse array.

In the process of using this library I was doing a lot of kron, reshape, etc. operations and found myself wishing that there was some way to represent the order of nesting of "combined" dimensions like you get after a kron operation. Something like this:

shape = NamedDimsArray{(:y, :y)}(I(ny))
block = NamedDimsArray{(:x, :x)}(spdiagm(-1 => -ones(nx-1), 1 => ones(nx-1)) ./ 2dx
kron(shape, block) # would be of type NamedDimsArray{((:y, :x), (:y, :x))}

Often the eventual destination of an N-dimensional named array will be a 2-dimensional MatOrVec, since you want to do some linear algebra on it eventually. In such cases it would be nice to keep track of whether you're working with an ":x-major" or a ":y-major" nested dimension.

[mcabbott]

In other packages... I messed around a bit with NamedPlus's split / join. And there is this:

julia> using Kronecker, NamedDims

julia> ab = NamedDimsArray(rand(Int8,3,4), (:a, :b));

julia> cd = NamedDimsArray(rand(Int8,3,4), (:c, :d));

julia> kronecker(ab, cd)
9×16 NamedDimsArray{(:aᵡc, :bᵡd), Int8, 2, Kronecker.KroneckerProduct{Int8, Matrix{Int8}, Matrix{Int8}}}:
  -72   123  -93     9   -80  -66   -50  -86     8    61   85   -33    40    49   -87    91
   68   -15    3    -1   -24   58  -114   38    92    55  -11    89   -52    19   -55   -67
...

julia> kronecker(cd, ab)
9×16 NamedDimsArray{(:cᵡa, :dᵡb), Int8, 2, Kronecker.KroneckerProduct{Int8, Matrix{Int8}, Matrix{Int8}}}:
  -72   -80     8    40   123  -66    61    49  -93   -50   85   -87     9  -86   -33    91
    0  -128    96    72    96   80    60   101   96   -48   92    61    32  112  -108  -105
...