Optimizing determinants and traces of Kron (#78)

* Optimize determinant and trace calculation for Kron

* Add tests for det and tr of Kron

* Optimize determinant and trace calculation for Kron

* Add tests for det and tr of Kron

* minor

* Kron Tests: compare determinant calculation against Julia's implementation for kronecker of rectangular factors

Author: emerali

src/LazyArrays.jl

@@ -6,7 +6,7 @@ module LazyArrays
 using Base, Base.Broadcast, LinearAlgebra, FillArrays, StaticArrays, ArrayLayouts
 import LinearAlgebra.BLAS
 
-import Base: AbstractArray, AbstractMatrix, AbstractVector, 
+import Base: AbstractArray, AbstractMatrix, AbstractVector,
         ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
              RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
              getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
@@ -36,8 +36,8 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcas
                         combine_eltypes, DefaultArrayStyle, instantiate, materialize,
                         materialize!, eltypes
 
-import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky, 
-                        norm2, norm1, normInf, normMinusInf
+import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky,
+                        norm2, norm1, normInf, normMinusInf, det, tr
 
 import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 
@@ -55,8 +55,8 @@ if VERSION < v"1.2-"
     import Base: has_offset_axes
     require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
 else
-    import Base: require_one_based_indexing    
-end             
+    import Base: require_one_based_indexing
+end
 
 export Mul, Applied, MulArray, MulVector, MulMatrix, InvMatrix, PInvMatrix,
         Hcat, Vcat, Kron, BroadcastArray, BroadcastMatrix, BroadcastVector, cache, Ldiv, Inv, PInv, Diff, Cumsum,

src/lazyoperations.jl

@@ -23,6 +23,42 @@ axes(a::Kron{<:Any,1}) = (OneTo(size(a,1)),)
 axes(a::Kron{<:Any,2}) = (OneTo(size(a,1)), OneTo(size(a,2)))
 axes(a::Kron{<:Any,N}) where {N} = (@_inline_meta; ntuple(M -> OneTo(size(a, M)), Val(N)))
 
+
+function det(K::Kron{<:Any, 2})
+    (size(K, 1) == size(K, 2)) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(K))"))
+
+    d = 1.
+    s = size(K, 1)
+
+    for A in K.args
+        if size(A, 1) == size(A, 2)
+            dA = det(A)
+            if iszero(dA)
+                return dA
+            end
+            d *= dA^(s ÷ size(A, 1))
+        else
+            # The Kronecker Product of rectangular matrices, if it is square, will
+            # have determinant zero. This can be shown by using the fact that
+            # rank(A ⊗ B) = rank(A)rank(B) and showing that this is strictly less
+            # than the number of rows in the resulting Kronecker matrix. Hence,
+            # since A ⊗ B does not have full rank, its determinant must be zero.
+            return zero(d)
+        end
+    end
+    return d
+end
+
+function tr(K::Kron{<:Any, 2})
+    (size(K, 1) == size(K, 2)) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(K))"))
+    if all(A -> (size(A, 1) == size(A, 2)), K.args)  # check if all component matrices are square
+        return prod(tr.(K.args))
+    else
+        return sum(diag(K))
+    end
+end
+
+
 getindex(K::Kron{T,1,<:Tuple{<:AbstractVector}}, k::Int) where T =
     first(K.args)[k]
 

test/runtests.jl

@@ -58,14 +58,35 @@ include("broadcasttests.jl")
 
     A = randn(3,2)
     B = randn(4,6)
-    K = Kron(A,B)
-    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == kron(A,B)
-    K = Kron(A,B')
-    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B')) == kron(A,B')
-    K = Kron(A',B)
-    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B)) == kron(A',B)
-    K = Kron(A',B')
-    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B')) == kron(A',B')
+    K, k = Kron(A,B), kron(A,B)
+    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B)) == k
+    @test det(K) == 0  # kronecker of rectangular factors
+    @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
+    @test tr(K) ≈ tr(k)
+
+    K, k = Kron(A,B'), kron(A,B')
+    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B')) == k
+    @test_throws DimensionMismatch det(K)
+    @test_throws DimensionMismatch tr(K)
+
+    K, k = Kron(A',B), kron(A',B)
+    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B)) == k
+    @test_throws DimensionMismatch det(K)
+    @test_throws DimensionMismatch tr(K)
+
+    K, k = Kron(A',B'), kron(A',B')
+    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A',B')) == k
+    @test det(K) == 0  # kronecker of rectangular factors
+    @test isapprox(det(k), det(K); atol=eps(eltype(K)), rtol=0)
+    @test tr(K) ≈ tr(k)
+
+    A = randn(3,3)
+    B = randn(6,6)
+    C = randn(2,2)
+    K, k = Kron(A,B,C), kron(A,B,C)
+    @test [K[k,j] for k=1:size(K,1), j=1:size(K,2)] == Array(Kron(A,B,C)) == k
+    @test det(K) ≈ det(k)
+    @test tr(K) ≈ tr(k)
 
     A = randn(3,2)
     B = randn(4,6)
@@ -157,7 +178,7 @@ end
     end
 
     @testset "colsupport past size" begin
-        C = cache(Zeros(5,5)); C[5,1]; 
+        C = cache(Zeros(5,5)); C[5,1];
         @test colsupport(C,1) == Base.OneTo(5)
         @test colsupport(C,3) == 1:0
         @test rowsupport(C,1) == Base.OneTo(1)
@@ -258,12 +279,12 @@ end
         v == BroadcastVector(exp, [1,2,3]) == exp.([1,2,3])
 
     Base.IndexStyle(typeof(BroadcastVector(exp, [1,2,3]))) == IndexLinear()
-    
+
     bc = broadcasted(exp,[1 2; 3 4])
     M = BroadcastArray(exp, [1 2; 3 4])
     @test BroadcastArray(bc) == BroadcastMatrix(bc) == BroadcastMatrix{Float64,typeof(exp),typeof(bc.args)}(bc) ==
         M == BroadcastMatrix(BroadcastMatrix(bc)) == BroadcastMatrix(exp,[1 2; 3 4]) == exp.([1 2; 3 4])
-    
+
     @test exp.(v') isa BroadcastMatrix
     @test exp.(transpose(v)) isa BroadcastMatrix
     @test exp.(M') isa BroadcastMatrix