add Kronecker products for LinearMaps (#61)

* add Kronecker products for LinearMaps

* announce feature in README

* add promotion of matrices, tests

* add comments, use mixed-product trick

* export UniformScalingMaps, add scalar mul, tests

* simplify scalar multiplication of composites

* adopt to new test suite

* prepare for require_one_based_indexing

* revert export of UnifformScalingMap, add LinearMap constructor instead

* update docstring

* add docstring to kron

* add KroneckerSumMap, yet untested

* add KroneckerPowerMap, yet untested

* KroneckerMap with many maps, optimize mul, more tests

* add dim-checker helper, resolve type instability, tests

* add conversion rules

* optimize MatrixMap and scaling multiplication, enable KroneckerSum tests

* fix matrix map handling across different julia versions

* minor fixes, improve coverage

* improve coverage, minor fixes

* improve coverage further

* remove false type-instability comment, uniform scaling comparison

* require one-based indexing, move kronecker power, simplify size

* clean up Kronecker (sum) powers

* make kron(sum)powers inferrable

* remove obsolete specifiers

* rename promote_to_lmaps -> convert_to_lmaps to avoid conflict with BlockMaps

* move conversion and check_dim to Main, pepper with at-propagate_inbounds

* make at-inbounds work properly, add some conversion rules

* simplify kronsum constructor

* edits as per review, update docstrings

* use memory requirement heuristic for branching

Author: dkarrasch

README.md

@@ -7,6 +7,13 @@
 
 A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently.
 
+## What's new in v2.6
+*   New feature: "lazy" Kronecker product, Kronecker sums, and powers thereof
+    for `LinearMap`s. `AbstractMatrix` objects are promoted to `LinearMap`s if
+    one of the first 8 Kronecker factors is a `LinearMap` object.
+*   Compatibility with the generic multiply-and-add interface (a.k.a. 5-arg
+    `mul!`) introduced in julia v1.3
+
 ## What's new in v2.5.0
 *   New feature: concatenation of `LinearMap`s objects with `UniformScaling`s,
     consistent with (h-, v-, and hc-)concatenation of matrices. Note, matrices

src/LinearMaps.jl

@@ -1,6 +1,7 @@
 module LinearMaps
 
 export LinearMap
+export ⊗, kronsum, ⊕
 
 using LinearAlgebra
 using SparseArrays
@@ -14,6 +15,8 @@ end
 
 abstract type LinearMap{T} end
 
+const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}}
+
 Base.eltype(::LinearMap{T}) where {T} = T
 
 Base.isreal(A::LinearMap) = eltype(A) <: Real
@@ -26,6 +29,26 @@ Base.ndims(::LinearMap) = 2
 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
 Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
 
+# check dimension consistency for y = A*x and Y = A*X
+function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector)
+     m, n = size(A)
+     (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!"))
+     return nothing
+ end
+ function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)
+     m, n = size(A)
+     (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
+     return nothing
+ end
+
+# conversion of AbstractMatrix to LinearMap
+convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A)
+convert_to_lmaps_(A::LinearMap) = A
+convert_to_lmaps() = ()
+convert_to_lmaps(A) = (convert_to_lmaps_(A),)
+@inline convert_to_lmaps(A, B, Cs...) =
+    (convert_to_lmaps_(A), convert_to_lmaps_(B), convert_to_lmaps(Cs...)...)
+
 Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
 function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false)
     length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
@@ -68,64 +91,24 @@ A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector)  = mul!(y, A,
 At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
 Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
 
-# Matrix: create matrix representation of LinearMap
-function Base.Matrix(A::LinearMap)
-    M, N = size(A)
-    T = eltype(A)
-    mat = Matrix{T}(undef, (M, N))
-    v = fill(zero(T), N)
-    @inbounds for i = 1:N
-        v[i] = one(T)
-        mul!(view(mat, :, i), A, v)
-        v[i] = zero(T)
-    end
-    return mat
-end
-
-Base.Array(A::LinearMap) = Matrix(A)
-Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A)
-Base.convert(::Type{Array}, A::LinearMap) = Matrix(A)
-Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
-
-# sparse: create sparse matrix representation of LinearMap
-function SparseArrays.sparse(A::LinearMap{T}) where {T}
-    M, N = size(A)
-    rowind = Int[]
-    nzval = T[]
-    colptr = Vector{Int}(undef, N+1)
-    v = fill(zero(T), N)
-    Av = Vector{T}(undef, M)
-
-    for i = 1:N
-        v[i] = one(T)
-        mul!(Av, A, v)
-        js = findall(!iszero, Av)
-        colptr[i] = length(nzval) + 1
-        if length(js) > 0
-            append!(rowind, js)
-            append!(nzval, Av[js])
-        end
-        v[i] = zero(T)
-    end
-    colptr[N+1] = length(nzval) + 1
-
-    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
-end
-
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
 include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
 include("transpose.jl") # transposing linear maps
 include("linearcombination.jl") # defining linear combinations of linear maps
 include("composition.jl") # composition of linear maps
 include("functionmap.jl") # using a function as linear map
 include("blockmap.jl") # block linear maps
+include("kronecker.jl") # Kronecker product of linear maps
+include("conversion.jl") # conversion of linear maps to matrices
 
 """
     LinearMap(A; kwargs...)
+    LinearMap(J, M::Int)
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)
 
 Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
-with the purpose of redefining its properties via the keyword arguments `kwargs`, or
+with the purpose of redefining its properties via the keyword arguments `kwargs`;
+a `UniformScaling` object `J` with specified (square) dimension `M`; or
 from a function or callable object `f`. In the latter case, one also needs to specify
 the size of the equivalent matrix representation `(M, N)`, i.e. for functions `f` acting
 on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably,
@@ -141,20 +124,21 @@ The keyword arguments and their default values for functions `f` are
 For existing linear maps or matrices `A`, the default values will be taken by calling
 `issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.
 
-For functions `f`, there is one more keyword arguments
+For functions `f`, there is one more keyword argument
 *   ismutating::Bool : flags whether the function acts as a mutating matrix multiplication
     `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
     or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
     The default value is guessed by looking at the number of arguments of the first occurence
     of `f` in the method table.
 """
-LinearMap(A::Union{AbstractMatrix, LinearMap}; kwargs...) = WrappedMap(A; kwargs...)
+LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...)
+LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M)
 LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
 LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)
 LinearMap(f, fc, M::Int; kwargs...) = LinearMap{Float64}(f, fc, M; kwargs...)
 LinearMap(f, fc, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, fc, M, N; kwargs...)
 
-LinearMap{T}(A::Union{AbstractMatrix, LinearMap}; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
+LinearMap{T}(A::MapOrMatrix; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
 LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)
 
 end # module

src/composition.jl

@@ -1,9 +1,12 @@
+# helper function
+check_dim_mul(A, B) = size(A, 2) == size(B, 1)
+
 struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     maps::As # stored in order of application to vector
     function CompositeMap{T, As}(maps::As) where {T, As}
         N = length(maps)
         for n in 2:N
-            size(maps[n], 2) == size(maps[n-1], 1) || throw(DimensionMismatch("CompositeMap"))
+            check_dim_mul(maps[n], maps[n-1]) || throw(DimensionMismatch("CompositeMap"))
         end
         for n in 1:N
             promote_type(T, eltype(maps[n])) == T || throw(InexactError())
@@ -47,7 +50,7 @@ function Base.:(*)(α::Number, A::CompositeMap)
     T = promote_type(typeof(α), eltype(A))
     Alast = last(A.maps)
     if Alast isa UniformScalingMap
-        return CompositeMap{T}(tuple(Base.front(A.maps)..., UniformScalingMap(α * Alast.λ, size(Alast, 1))))
+        return CompositeMap{T}(tuple(Base.front(A.maps)..., α * Alast))
     else
         return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
     end
@@ -60,7 +63,7 @@ function Base.:(*)(A::CompositeMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
     Afirst = first(A.maps)
     if Afirst isa UniformScalingMap
-        return CompositeMap{T}(tuple(UniformScalingMap(Afirst.λ * α, size(Afirst, 1)), Base.tail(A.maps)...))
+        return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...))
     else
         return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
     end
@@ -87,22 +90,22 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
 ```
 """
 function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
-    size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
+    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂, A₁))
 end
 function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
-    size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
+    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂.maps..., A₁))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
-    size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
+    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂, A₁.maps...))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
-    size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*"))
+    check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
 end

src/conversion.jl

@@ -0,0 +1,76 @@
+# Matrix: create matrix representation of LinearMap
+function Base.Matrix(A::LinearMap)
+    M, N = size(A)
+    T = eltype(A)
+    mat = Matrix{T}(undef, (M, N))
+    v = fill(zero(T), N)
+    @inbounds for i = 1:N
+        v[i] = one(T)
+        mul!(view(mat, :, i), A, v)
+        v[i] = zero(T)
+    end
+    return mat
+end
+Base.Array(A::LinearMap) = Matrix(A)
+Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A)
+Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A)
+Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A)
+Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A)
+
+# special cases
+Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
+Base.convert(::Type{Matrix}, A::WrappedMap) = convert(Matrix, A.lmap)
+function Base.convert(::Type{Matrix}, ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}})
+    if length(ΣA.maps) <= 10
+        return (+).(map(A->getfield(A, :lmap), ΣA.maps)...)
+    else
+        S = zero(first(ΣA.maps).lmap)
+        for A in ΣA.maps
+            S .+= A.lmap
+        end
+        return S
+    end
+end
+function Base.convert(::Type{Matrix}, AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}})
+    B, A = AB.maps
+    return A.lmap*B.lmap
+end
+function Base.convert(::Type{Matrix}, λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
+    A, J = λA.maps
+    return J.λ*A.lmap
+end
+function Base.convert(::Type{Matrix}, Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
+    J, A = Aλ.maps
+    return A.lmap*J.λ
+end
+
+Base.Matrix(A::BlockMap) = hvcat(A.rows, convert.(AbstractMatrix, A.maps)...)
+
+# sparse: create sparse matrix representation of LinearMap
+function SparseArrays.sparse(A::LinearMap{T}) where {T}
+    M, N = size(A)
+    rowind = Int[]
+    nzval = T[]
+    colptr = Vector{Int}(undef, N+1)
+    v = fill(zero(T), N)
+    Av = Vector{T}(undef, M)
+
+    for i = 1:N
+        v[i] = one(T)
+        mul!(Av, A, v)
+        js = findall(!iszero, Av)
+        colptr[i] = length(nzval) + 1
+        if length(js) > 0
+            append!(rowind, js)
+            append!(nzval, Av[js])
+        end
+        v[i] = zero(T)
+    end
+    colptr[N+1] = length(nzval) + 1
+
+    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
+end
+
+Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
+Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
+Base.convert(::Type{SparseMatrixCSC}, A::BlockMap) = hvcat(A.rows, convert.(SparseMatrixCSC, A.maps)...)