MichielStock / Kronecker.jl

A general-purpose toolbox for efficient Kronecker-based algebra.
MIT License
86 stars 14 forks source link

Resolve ambiguity in multiplying Diagonal kronecker products #113

Closed jishnub closed 2 years ago

jishnub commented 2 years ago

Fixes the following multiplications:

julia> A = randn(4, 4);

julia> B = Array{Float64, 2}([1 2 3;
                4 5 6;
                7 -2 9]);

julia> K = A ⊗ B;

julia> Kd = Diagonal(A) ⊗ Diagonal(B);

julia> Kd2 = Diagonal(B) ⊗ Diagonal(A);

julia> Kd * Kd2 ≈ Diagonal(Kd) * Diagonal(Kd2)
true

julia> Kd2 * Kd ≈ Diagonal(Kd2) * Diagonal(Kd)
true

julia> K * Kd ≈ collect(K) * Diagonal(Kd)
true

julia> Kd * K ≈ Diagonal(Kd) * collect(K)
true
codecov-commenter commented 2 years ago

Codecov Report

Merging #113 (aaa0b29) into master (70fc882) will increase coverage by 0.04%. The diff coverage is 96.87%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #113      +/-   ##
==========================================
+ Coverage   90.43%   90.48%   +0.04%     
==========================================
  Files          10       10              
  Lines         805      809       +4     
==========================================
+ Hits          728      732       +4     
  Misses         77       77              
Impacted Files Coverage Δ
src/base.jl 92.88% <96.55%> (+0.28%) :arrow_up:
src/vectrick.jl 94.53% <100.00%> (-0.29%) :arrow_down:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 70fc882...aaa0b29. Read the comment docs.

MichielStock commented 2 years ago

I've noticed that Diagonal(K) returns are regular Diagonal matrix rather than a KronProdDiagonal. Would it make sence to also extend this function a la

LinearAlgebra.Diagonal(K::AbstractKronecker) = mapreduce(Diagonal, ⊗, getmatrices(K))
jishnub commented 2 years ago

Diagonal refers to the constructor, therefore it should return an object of the same type. Perhaps we may add a function diagonal to Kronecker that will return the appropriate type? Ideally there could have been such a function in LinearAlgebra that we may overload here, but as of now there's only the constructor that exists.

MichielStock commented 2 years ago

For many cases, the difference between a diagonal or a Kronecker type might be minimal. Though it might be useful to add two Kronecker produces?

Agreed that we should not mess with constructors! LinearAlgebra.diagm might fit the bill though.

jishnub commented 2 years ago

diagm constructs a matrix given the vector that corresponds to the diagonal. Perhaps we're looking for something like diagm ∘ diag, but at that point it's probably easier to define our own function

MichielStock commented 2 years ago

You are right, my bad! I thought it also created a diagonal matrix from a square matrix. If we agree that we want the behaviour as discussed, I am in favor of our own function

diagonal(A::AbstractMatrix) = Diagonal(A)
diagonal(K::AbstractKroneckerProduct) = mapreduce(Diagonal, ⊗, getmatrices(K))