Fix order of multiplication in mixed-product rule (#231)

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.11.2"
+version = "3.11.3"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/LinearMaps.jl

@@ -99,6 +99,9 @@ _combine(As::LinearMapVector, Bs::LinearMapTuple) = Base.vect(As..., Bs...)
 _combine(As::LinearMapTuple, Bs::LinearMapVector) = Base.vect(As..., Bs...)
 _combine(As::LinearMapVector, Bs::LinearMapVector) = Base.vect(As..., Bs...)
 
+_reverse!(As::LinearMapTuple) = reverse(As)
+_reverse!(As::LinearMapVector) = reverse!(As)
+
 # The (internal) multiplication logic is as follows:
 #  - `*(A, x)` calls `mul!(y, A, x)` for appropriately-sized y
 #  - `mul!` checks consistency of the sizes, and calls `_unsafe_mul!`,

src/composition.jl

@@ -15,9 +15,9 @@ end
 CompositeMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} = CompositeMap{T, As}(maps)
 
 Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTupleOrVector) =
-    CompositeMap{promote_type(map(eltype, maps)...)}(reverse(maps))
+    CompositeMap{promote_type(map(eltype, maps)...)}(_reverse!(maps))
 Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector{<:LinearMap{T}}) where {T} =
-    CompositeMap{T}(reverse(maps))
+    CompositeMap{T}(reverse!(maps))
 
 MulStyle(A::CompositeMap) = MulStyle(A.maps...) === TwoArg() ? TwoArg() : ThreeArg()
 
@@ -158,9 +158,9 @@ Base.:(*)(A₁::CompositeMap, A₂::ScaledMap) = (A₁ * A₂.lmap) * A₂.λ
 
 # special transposition behavior
 LinearAlgebra.transpose(A::CompositeMap{T}) where {T} =
-    CompositeMap{T}(map(transpose, reverse(A.maps)))
+    CompositeMap{T}(map(transpose, _reverse!(A.maps)))
 LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} =
-    CompositeMap{T}(map(adjoint, reverse(A.maps)))
+    CompositeMap{T}(map(adjoint, _reverse!(A.maps)))
 
 # comparison of CompositeMap objects
 Base.:(==)(A::CompositeMap, B::CompositeMap) =
@@ -169,7 +169,7 @@ Base.:(==)(A::CompositeMap, B::CompositeMap) =
 # multiplication with vectors/matrices
 function Base.:(*)(A::CompositeMap, x::AbstractVector)
     MulStyle(A) === TwoArg() ?
-        foldr(*, reverse(A.maps), init=x) :
+        foldr(*, _reverse!(A.maps), init=x) :
         invoke(*, Tuple{LinearMap, AbstractVector}, A, x)
 end
 

src/kronecker.jl

@@ -272,7 +272,7 @@ function _unsafe_mul!(y,
     Bs1, Bs2 = _front(Bs), _tail(Bs)
     apply = all(_iscompatible, zip(As1, As2)) && all(_iscompatible, zip(Bs1, Bs2))
     if apply
-        _unsafe_mul!(y, kron(prod(As), prod(Bs)), x)
+        _unsafe_mul!(y, kron(prod(_reverse!(As)), prod(_reverse!(Bs))), x)
     else
         _unsafe_mul!(y, CompositeMap{T}(map(LinearMap, L.maps)), x)
     end

test/kronecker.jl

@@ -86,6 +86,10 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
         Kv = LinearMaps.CompositeMap{ComplexF64}(fill(LA ⊗ LB, 3))
         @test kron(A, B)^3 * ones(6) ≈ Kv * ones(6)
         @test Matrix(K) ≈ kron(A, B)^3
+        A = [0 1; 0 0]
+        B = [0 0; 1 0]
+        J = LinearMap(I, 1)
+        @test Matrix(kron(J, A*B)) == Matrix(kron(J, A) * kron(J, B))
         # example that doesn't use mixed-product rule
         A = rand(3, 2); B = rand(2, 3)
         K = @inferred kron(A, LinearMap(B))