Lazy sums of KroneckerProduct

[vsaase]

Hi,

how would you implement a lazy sum of KroneckerProducts?

size(A) == size(C)
size(B) == size(D)
E = kronecker(A,B)
F = kronecker(C,D)
G = E+F

The problem is that when computing G, I am loosing the memory efficiency of kronecker, and I only want to do efficient Matrix vector multiplication down the line (G*v as E*v + F*v), without keeping track of all summands

Does this look like a reasonable addition to Kronecker.jl, or is this use case maybe too special?

[dkarrasch]

This functionality exists in LinearMaps.jl, which has, by the way, it's own lazy Kronecker product implemented. So you have two options:

# use kronecker from Kronecker.jl and the linear combination from LinearMaps.jl
using Kronecker, LinearMaps
E, F = ... as above
G = LinearMap(E) + F

# use LinearMaps.jl for both the Kronecker product and +
E = kron(LinearMap(A), B)
F = kron(LinearMap(C), D)
G = E + F

I guess performance should be identical, but I'd love to know if not. Since LinearMaps.jl can wrap any AbstractMatrix and Kronecker.jl returns such, I don't think it makes sense to include addition of KroneckerProduct <: AbstractMatrix in this package.

[MichielStock]

The most straightforward way is to use the distributive property E * v + F * v. Though, if you don't want to keep track of this, it should be easy to extend this. I made a working example:

struct SumOfKroneckers{T<:Any, TA<:GeneralizedKroneckerProduct, TB<:AbstractMatrix} <: GeneralizedKroneckerProduct{T}
    KA::TA
    KB::TB
    function SumOfKroneckers(KA::GeneralizedKroneckerProduct{T}, KB::AbstractMatrix{V}) where {T, V}
        return new{promote_type(T, V), typeof(KA), typeof(KB)}(KA, KB)
    end
end

Base.size(K::SumOfKroneckers) = size(K.KA)

Base.getindex(K::SumOfKroneckers, i1::Integer, i2::Integer) = K.KA[i1,i2] + K.KB[i1,i2] 

function Base.:+(KA::GeneralizedKroneckerProduct, KB::AbstractMatrix)
    @assert size(KA) == size(KB) "`A` and `B` have to be conformable"
    return SumOfKroneckers(KA, KB)
end

Base.:*(K::SumOfKroneckers, V::VecOrMat) = K.KA * V .+ K.KB * V


A = rand(10, 10)
B = rand(Bool, 5, 5)
C = randn(10, 10)
D = rand(1:100, 5, 5)

KA = A ⊗ B
KB = C ⊗ D

K = KA + KB

v = randn(50)

This might overlap with LinearMap, but it might be worthwhile adding to Kronecker. If so, I can do it within a couple of days.

[MichielStock]

See #77

[vsaase]

wow, thanks for the quick implementation. I will test it tomorrow