Define non-vectorized multiplication for Kronecker

[elisno]

Due to the reshape in mul! https://github.com/MichielStock/Kronecker.jl/blob/b4101121ace4ad8834fed81504a2082a23dbd6ce/src/base.jl#L190-L193

, calling kronecker(M,N) * vec(V) is unnecessarily slow when V has some special structure.

In the example below, V::Diagonal.

using LinearAlgebra, Kronecker, Test
N = 256
A = rand(N,N)
B = rand(N,N)
V = Diagonal(rand(N))
K = kronecker(transpose(B),A)

A*V*B # Warm-up
reshape(K * vec(V),(N,N)) # Warm-up

@time m1 = A*V*B # 0.001687 seconds (9 allocations: 1.000 MiB)
@time m2 = reshape(K * vec(V),(N,N)) # 0.148023 seconds (84 allocations: 5.005 MiB)
@test isapprox(m1,m2)

If we define something like

function mul!(x::AbstractMatrix K::AbstractKroneckerProduct, V::AbstractMatrix)
    M, N = getmatrices(K)
    ...
    x .= N * V * transpose(M)
    return x
end

then vec(kronecker(M,N) * V) would be equivalent to kronecker(M,N) * vec(V).

[MichielStock]

This is definitely a good issue. The only thing that I don't like about this is dat vec(kronecker(M,N) * V) is mathematically incorrect, as it would only hold for N^2 x M matrices to be conformable.

Maybe a better solution is to define a function vecmult(X::Abstractmatrix, K::AbstractKroneckerProduct) which computes kronecker(M,N) * vec(V) without reshaping V.

(This could actually be a job for metaprogramming, but might be overkill).

[elisno]

vec(kronecker(M,N) * V) is mathematically incorrect, as it would only hold for N^2 x M matrices to be conformable.

That's a good point. I wasn't thinking generically enough.

For my use case of this, I don't need to use vec. I need to store the original matrix multiplications to solve a matrix equation containing many commutators of the form:

Eagerly evaluating these smaller kronecker products/sums and storing them in an AbstractKroneckerProduct is relatively fast and it simplifies the equations.

I'm not sure if any performance is gained by writing a similar scheme for nested (higher order) kronecker products reshape( (transpose((I⊗t')) ⊗ (I⊗t))vec(V) , size(V) ).

[MichielStock]

higher order - vector matrix products are evaluated correctly by recursion, but I think the performance can drastically be improved (for example using (TensorOperations.jl)[https://github.com/Jutho/TensorOperations.jl]. This is definitely a working point!

[elisno]

which computes kronecker(M,N) * vec(V) without reshaping V

Right now it's unclear to me, but vec(A::AbstractArray) returns a ReshapedArray{eltype(A),1,typeof(A),...} in most cases. There, the parent array A is still stored.

Perhaps defining the Kronecker mul! for v::ReshapedArray{T<:Number,1} which returns

x .= vec(N * v.parent * transpose(M))

Maybe a better solution is to define a function vecmult(X::Abstractmatrix, K::AbstractKroneckerProduct)

It makes sense that the "vec-trick" without vec should use a different function name.

[MichielStock]

This is a splendid idea and is implemented in f935bc2719a945160d69a103a6f76549aaca52d1!

[elisno]

28 closes this issue.