Store identity matrices in KroneckerSums

[elisno]

Should the KroneckerSum have extra fields for the identity matrices?

Something like the following:

struct KroneckerSum{T<:Any, TA<:AbstractMatrix, TB<:AbstractMatrix} <: AbstractKroneckerSum{T}
    A::TA
    B::TB
    IA::TA
    IB::TB
    function KroneckerSum(A::AbstractMatrix{T},
                            B::AbstractMatrix{V}) where {T, V}
        (issquare(A) && issquare(B)) || throw(DimensionMismatch(
                                "KroneckerSum only applies to square matrices"))
        return new{promote_type(T, V), typeof(A), typeof(B)}(A, B, oneunit(A), oneunit(B))
    end
end

IA and IB would be used in Base.collect(K::KroneckerSum).

Then it would make sense to "reuse" a KroneckerSum, when repeated collections are required.

It would also complement in-place versions of collect with kron! (see #61).

[MichielStock]

I don't see the value of generating and storing two dense matrices? If we collect as you say, we better generate them on the fly once? Especially for in-place collecting, It does not make much sense to me to use the dense identity matrices. For other processes, an UniformScaling would be more appropriate, no?

[elisno]

I don't see the value of generating and storing two dense matrices? ... ... , an UniformScaling would be more appropriate, no?

I may have jumped the gun with reusing the types TA and TB, as well as using oneunit in the example to initialize the identity matrices, whoops!

I agree that storing the identity matrices as dense arrays is highly inefficient. Although I don't show it in my example, my idea was that if A or B are dense, the corresponding identity matrices would be Diagonal. If A and B are sparse, then the identity matrices should be stored as such as well. Allowing this flexibility would speed up sums with smaller matrices (dense)

I'll admit that I'm not familiar enough with UniformScaling. Were you thinking of using

(I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))

to generate the identity matrices?

Collecting dense products with smaller matrices would be faster than by making them sparse, so there should be some flexibility in the storage of these matrices.

Should the user be responsible for selecting the storage type, or should we set some condition for choosing between types?

---

If we collect as you say, we better generate them on the fly once?

That is the key idea! I like this approach.

[MichielStock]

I agree that storing the identity matrices as dense arrays is highly inefficient. Although I don't show it in my example, my idea was that if A or B are dense, the corresponding identity matrices would be Diagonal. If A and B are sparse, then the identity matrices should be stored as such as well. Allowing this flexibility would speed up sums with smaller matrices (dense)

I'll admit that I'm not familiar enough with UniformScaling. Were you thinking of using

(I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))

to generate the identity matrices?

The uniform scaling is an efficient way of using the identity matrix, as A + I and I * A will be computed efficiently. This might not be needed for our purposes, since the underlying code can work without them being generated. For example A ⊗ I can be collected using for-loops or more abstract functions to form the block-diagonal matrix.

Collecting dense products with smaller matrices would be faster than by making them sparse, so there should be some flexibility in the storage of these matrices.

Should the user be responsible for selecting the storage type, or should we set some condition for choosing between types?

I struggled with this question in other parts of the code as well. Some storage types will be more efficient depending on the storage type, and I found that the theoretical time complexities do not always match reality as it also depends on hardware stuff. I lean toward letting the responsibility to the user, where they have to determine which form they use. A good documentation would be key here.

A fun thing to implement might be a macro @flops that computes the number of element-wise computations for a given expression using Kronecker. This could guide the uses in choosing the right function.

If we collect as you say, we better generate them on the fly once?

That is the key idea! I like this approach.

But this is kind of already done in most of the code.

It might help if you give me some examples of what kind of things you want to compute (how many matrices and what sizes). My background is not physics. The main interest (in my research) I have is in forming systems with a limited number of large matrices. What do you compute with Kronecker sums? It might be, for example, be interesting to store many terms of pure Kronecker products, these would allow one to lazily add Kronecker sums and mixtures. Many linear algebra operations (det, eigen,...) cannot be computed efficiently anymore though.