`kron` between Diagonal and non-numeric matrix gives output with `undef`s

[lkapelevich]

on 1.7-dev:

julia> using JuMP, LinearAlgebra

julia> X = ones(AffExpr, 2, 2);

julia> kron(I(2), X)
4×4 Matrix{AffExpr}:
    1       1    #undef  #undef
    1       1    #undef  #undef
 #undef  #undef     1       1
 #undef  #undef     1       1

Instead of dispatching to https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/diagonal.jl#L540 which initializes the output with zeros, we go into https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/dense.jl#L414

[nikhilnarang2001]

Instead of using
kron(I(2), X);

you should use:
A = [1 0;0 1];
kron (A,X)

when you use A = I(2)
A = 1 undef
undef 1
this why it is showing undef when you use kron with I(n).

[JATippit]

Instead of using

kron(I(2), X);

you should use:

A = [1 0;0 1];

kron (A,X)

when you use A = I(2)

A = 1 undef

    undef      1

this why it is showing undef when you use kron with I(n).

I don't receive undef when replicating the original issue - I get 0s as expected.

[lkapelevich]

I don't receive undef when replicating the original issue - I get 0s as expected.

Really? Here is an example without JuMP. On both 1.6.0 and 1.7.0 (updated version in OP):

julia> using LinearAlgebra; X = [ones(1,1) for i in 1:2, j in 1:2];

julia> kron(I(2), X)
4×4 Matrix{Matrix{Float64}}:
    [1.0]     [1.0]  #undef    #undef
    [1.0]     [1.0]  #undef    #undef
 #undef    #undef       [1.0]     [1.0]
 #undef    #undef       [1.0]     [1.0]

[JATippit]

Really? Here is an example without JuMP. On both 1.6.0 and 1.7.0 (updated version in OP):

julia> using LinearAlgebra; X = [ones(1,1) for i in 1:2, j in 1:2];



julia> kron(I(2), X)

4×4 Matrix{Matrix{Float64}}:

    [1.0]     [1.0]  #undef    #undef

    [1.0]     [1.0]  #undef    #undef

 #undef    #undef       [1.0]     [1.0]

 #undef    #undef       [1.0]     [1.0]

I realized I tested with v1.5.4. Retried with v1.6.0 and get #undef.

The issue is because the diagonal function is not being dispatched since
Matrix{AffExpr} <: AbstractMatrix{<:Number}
returns false.

Not using JuMP, the way you've written it returns an Array of Arrays, which is defaulting to the same function. Try X = ones(2,2)

Removing the <:Number from S in diagonal.jl corrects the issue. Neither T nor S have this stipulation in dense.jl. Is there any reason we couldn't remove it here?

The difference seems to come from the fact X=ones(AffExpr,2,2) returns Array{GenericAffExpr{Float64,VariableRef},2} in Julia v1.5.4 but 2x2 Matrix{AffExpr} in v1.6.0.

[nikhilnarang2001]

The difference seems to come from the fact X=ones(AffExpr,2,2) returns Array{GenericAffExpr{Float64,VariableRef},2} in Julia v1.5.4 but 2x2 Matrix{AffExpr} in v1.6.0.

Actually, in v1.6.0 the return type is the same as v1.5.4 just its alias name has been used.

julia> @show typeof(ones(AffExpr,2,2))
typeof(ones(AffExpr, 2, 2)) = Matrix{AffExpr}
Matrix{AffExpr} (alias for Array{GenericAffExpr{Float64, VariableRef}, 2})

I think the output varries according to the data type we are using, I tried some different methods and found that it works when we use this,

julia> @show typeof(ones(2,2))
typeof(ones(2, 2)) = Matrix{Float64}
Matrix{Float64} (alias for Array{Float64, 2})

and shows this with kron,

julia> kron(I(2),ones(2,2))
4×4 Matrix{Float64}:
 1.0  1.0  0.0  0.0
 1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0
 0.0  0.0  1.0  1.0

but when we use this,

julia> X = [ones(1,1) for i in 1:2, j in 1:2];

julia> @show typeof(X)
typeof(X) = Matrix{Matrix{Float64}}
Matrix{Matrix{Float64}} (alias for Array{Array{Float64, 2}, 2})

julia> kron(I(2),X)
4×4 Matrix{Matrix{Float64}}:
    [1.0]     [1.0]  #undef    #undef
    [1.0]     [1.0]  #undef    #undef
 #undef    #undef       [1.0]     [1.0]
 #undef    #undef       [1.0]     [1.0]

This doesn't work. The problem is the datatype we use, it shows different for different type of data.

[vtjnash]

IIUC, this is a general bug in some of the diagonal kron methods. I think this would fix it, but needs someone to write a test:

diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl
index e83c5274ac..4a1dedc076 100644
--- a/stdlib/LinearAlgebra/src/diagonal.jl
+++ b/stdlib/LinearAlgebra/src/diagonal.jl
@@ -472,14 +472,13 @@ rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
     invoke(\, Tuple{Union{QR,QRCompactWY,QRPivoted}, AbstractVecOrMat}, A, B)
 
 
-@inline function kron!(C::AbstractMatrix{T}, A::Diagonal, B::Diagonal) where T
-    fill!(C, zero(T))
+@inline function kron!(C::AbstractMatrix, A::Diagonal, B::Diagonal)
     valA = A.diag; nA = length(valA)
     valB = B.diag; nB = length(valB)
     nC = checksquare(C)
     @boundscheck nC == nA*nB ||
         throw(DimensionMismatch("expect C to be a $(nA*nB)x$(nA*nB) matrix, got size $(nC)x$(nC)"))
-
+    isempty(A) || isempty(B) || fill!(C, zero(A[1,1] * B[1,1]))
     @inbounds for i = 1:nA, j = 1:nB
         idx = (i-1)*nB+j
         C[idx, idx] = valA[i] * valB[j]
@@ -497,10 +496,13 @@ end
 
 @inline function kron!(C::AbstractMatrix, A::Diagonal, B::AbstractMatrix)
     Base.require_one_based_indexing(B)
-    (mA, nA) = size(A); (mB, nB) = size(B); (mC, nC) = size(C);
+    (mA, nA) = size(A)
+    (mB, nB) = size(B)
+    (mC, nC) = size(C)
     @boundscheck (mC, nC) == (mA * mB, nA * nB) ||
         throw(DimensionMismatch("expect C to be a $(mA * mB)x$(nA * nB) matrix, got size $(mC)x$(nC)"))
     m = 1
+    isempty(A) || isempty(B) || fill!(C, zero(A[1,1] * B[1,1]))
     @inbounds for j = 1:nA
         A_jj = A[j,j]
         for k = 1:nB
@@ -517,10 +519,13 @@ end
 
 @inline function kron!(C::AbstractMatrix, A::AbstractMatrix, B::Diagonal)
     require_one_based_indexing(A)
-    (mA, nA) = size(A); (mB, nB) = size(B); (mC, nC) = size(C);
+    (mA, nA) = size(A)
+    (mB, nB) = size(B)
+    (mC, nC) = size(C)
     @boundscheck (mC, nC) == (mA * mB, nA * nB) ||
         throw(DimensionMismatch("expect C to be a $(mA * mB)x$(nA * nB) matrix, got size $(mC)x$(nC)"))
     m = 1
+    isempty(A) || isempty(B) || fill!(C, zero(A[1,1] * B[1,1]))
     @inbounds for j = 1:nA
         for l = 1:mB
             Bll = B[l,l]
@@ -535,18 +535,6 @@ end
     return C
 end
 
-function kron(A::Diagonal{T}, B::AbstractMatrix{S}) where {T<:Number, S<:Number}
-    (mA, nA) = size(A); (mB, nB) = size(B)
-    R = zeros(Base.promote_op(*, T, S), mA * mB, nA * nB)
-    return @inbounds kron!(R, A, B)
-end
-
-function kron(A::AbstractMatrix{T}, B::Diagonal{S}) where {T<:Number, S<:Number}
-    (mA, nA) = size(A); (mB, nB) = size(B)
-    R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
-    return @inbounds kron!(R, A, B)
-end
-
 conj(D::Diagonal) = Diagonal(conj(D.diag))
 transpose(D::Diagonal{<:Number}) = D
 transpose(D::Diagonal) = Diagonal(transpose.(D.diag))

edit: this code added in JuliaLang/julia#31069, which removed the calls to fill!(zero)