Allow vectors for LinearMap storage (#171)

Author: dkarrasch

Project.toml

@@ -1,10 +1,11 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.5.2"
+version = "3.6.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 julia = "1.6"

docs/Project.toml

@@ -6,6 +6,6 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 
 [compat]
-BenchmarkTools = "0.4, 0.5"
-Documenter = "0.25, 0.26"
+BenchmarkTools = "1"
+Documenter = "0.25, 0.26, 0.27"
 Literate = "2"

docs/src/history.md

@@ -1,5 +1,15 @@
 # Version history
 
+## What's new in v3.6
+
+* Support for Julia versions below v1.6 has been dropped.
+* `Block[Diagonal]Map`, `CompositeMap`, `KroneckerMap` and `LinearCombination` type objects
+  can now be backed by a `Vector` of `LinearMap`-type elements. This can be beneficial in
+  cases where these higher-order `LinearMap`s are constructed from many maps where a tuple
+  backend may get inefficient or impose hard work for the compiler at construction.
+  The default behavior, however, does not change, and construction of vector-based
+  `LinearMap`s requires usage of the unexported constructors ("expert usage").
+
 ## What's new in v3.5
 
 * `WrappedMap`, `ScaledMap`, `LinearCombination`, `AdjointMap`, `TransposeMap` and
@@ -24,7 +34,7 @@
 ## What's new in v3.3
 
 * `AbstractVector`s can now be wrapped by a `LinearMap` just like `AbstractMatrix``
-  typed objects. Upon wrapping, there are not implicitly reshaped to matrices. This
+  typed objects. Upon wrapping, they are not implicitly reshaped to matrices. This
   feature might be helpful, for instance, in the lazy representation of rank-1
   operators `kron(LinearMap(u), v') == ⊗(u, v') == u ⊗ v'` for vectors `u` and `v`.
   The action on vectors,`(u⊗v')*x`, is implemented optimally via `u*(v'x)`.

src/LinearMaps.jl

@@ -8,6 +8,8 @@ using LinearAlgebra
 import LinearAlgebra: mul!
 using SparseArrays
 
+import Statistics: mean
+
 using Base: require_one_based_indexing
 
 abstract type LinearMap{T} end
@@ -16,8 +18,21 @@ const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMat{T}}
 const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}}
 const RealOrComplex = Union{Real, Complex}
 
+const LinearMapTuple = Tuple{Vararg{LinearMap}}
+const LinearMapVector = AbstractVector{<:LinearMap}
+const LinearMapTupleOrVector = Union{LinearMapTuple,LinearMapVector}
+
 Base.eltype(::LinearMap{T}) where {T} = T
 
+# conversion to LinearMap
+Base.convert(::Type{LinearMap}, A::LinearMap) = A
+Base.convert(::Type{LinearMap}, A::AbstractVecOrMat) = LinearMap(A)
+
+convert_to_lmaps() = ()
+convert_to_lmaps(A) = (convert(LinearMap, A),)
+@inline convert_to_lmaps(A, B, Cs...) =
+    (convert(LinearMap, A), convert(LinearMap, B), convert_to_lmaps(Cs...)...)
+
 abstract type MulStyle end
 
 struct FiveArg <: MulStyle end
@@ -61,13 +76,20 @@ function check_dim_mul(C, A, B)
     return nothing
 end
 
-# conversion of AbstractVecOrMat to LinearMap
-convert_to_lmaps_(A::AbstractVecOrMat) = 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...)...)
+_front(As::Tuple) = Base.front(As)
+_front(As::AbstractVector) = @inbounds @views As[1:end-1]
+_tail(As::Tuple) = Base.tail(As)
+_tail(As::AbstractVector) = @inbounds @views As[2:end]
+
+_combine(A::LinearMap, B::LinearMap) = tuple(A, B)
+_combine(A::LinearMap, Bs::LinearMapTuple) = tuple(A, Bs...)
+_combine(As::LinearMapTuple, B::LinearMap) = tuple(As..., B)
+_combine(As::LinearMapTuple, Bs::LinearMapTuple) = tuple(As..., Bs...)
+_combine(A::LinearMap, Bs::LinearMapVector) = Base.vect(A, Bs...)
+_combine(As::LinearMapVector, B::LinearMap) = Base.vect(As..., B)
+_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...)
 
 # The (internal) multiplication logic is as follows:
 #  - `*(A, x)` calls `mul!(y, A, x)` for appropriately-sized y
@@ -233,8 +255,6 @@ function _unsafe_mul!(y::AbstractMatrix, A::LinearMap, x::AbstractMatrix, α, β
     return _generic_mapmat_mul!(y, A, x, α, β)
 end
 
-const LinearMapTuple = Tuple{Vararg{LinearMap}}
-
 include("left.jl") # left multiplication by a transpose or adjoint vector
 include("transpose.jl") # transposing linear maps
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties

src/blockmap.jl

@@ -1,33 +1,32 @@
 struct BlockMap{T,
-                As<:LinearMapTuple,
-                Rs<:Tuple{Vararg{Int}},
-                Rranges<:Tuple{Vararg{UnitRange{Int}}},
-                Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
+                As<:LinearMapTupleOrVector,
+                Rs<:Tuple{Vararg{Int}}} <: LinearMap{T}
     maps::As
     rows::Rs
-    rowranges::Rranges
-    colranges::Cranges
+    rowranges::Vector{UnitRange{Int}}
+    colranges::Vector{UnitRange{Int}}
     function BlockMap{T,As,Rs}(maps::As, rows::Rs) where
-                {T, As<:LinearMapTuple, Rs<:Tuple{Vararg{Int}}}
+                {T, As<:LinearMapTupleOrVector, Rs<:Tuple{Vararg{Int}}}
         for TA in Base.Generator(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in BlockMap constructor")
         end
         rowranges, colranges = rowcolranges(maps, rows)
-        Rranges, Cranges = typeof(rowranges), typeof(colranges)
-        return new{T, As, Rs, Rranges, Cranges}(maps, rows, rowranges, colranges)
+        return new{T, As, Rs}(maps, rows, rowranges, colranges)
     end
 end
 
-BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTuple, Rs} =
+BlockMap{T}(maps::As, rows::Rs) where {T, As<:LinearMapTupleOrVector, Rs} =
     BlockMap{T, As, Rs}(maps, rows)
+BlockMap(maps::As, rows::Rs) where {As<:LinearMapTupleOrVector, Rs} =
+    BlockMap{promote_type(map(eltype, maps)...), As, Rs}(maps, rows)
 
 MulStyle(A::BlockMap) = MulStyle(A.maps...)
 
-function _getranges(maps, dim, inds=ntuple(identity, Val(length(maps))))
-    sizes = map(i -> size(maps[i], dim)::Int, inds)
-    ends = cumsum(sizes)
-    starts = (1, (1 .+ Base.front(ends))...)
+function _getranges(maps, dim, inds=1:length(maps))
+    ends = map(i -> size(maps[i], dim)::Int, inds)
+    cumsum!(ends, ends)
+    starts = vcat(1, 1 .+ @views ends[1:end-1])
     return UnitRange.(starts, ends)
 end
 
@@ -40,14 +39,15 @@ block linear map obtained from `hvcat(rows, maps...)`.
 """
 function rowcolranges(maps, rows)
     # find indices of the row-wise first maps
-    firstmapinds = cumsum((1, Base.front(rows)...))
+    firstmapinds = vcat(1, Base.front(rows)...)
+    cumsum!(firstmapinds, firstmapinds)
     # compute rowranges from first dimension of the row-wise first maps
     rowranges = _getranges(maps, 1, firstmapinds)
 
     # compute ranges from second dimension as if all in one row
     temp = _getranges(maps, 2)
     # introduce "line breaks"
-    colranges = ntuple(Val(length(maps))) do i
+    colranges = map(1:length(maps)) do i
         # for each map find the index of the respective row-wise first map
         # something-trick just to assure the compiler that the index is an Int
         @inbounds firstmapind = firstmapinds[something(findlast(<=(i), firstmapinds), 1)]
@@ -277,7 +277,7 @@ end
 ############
 
 Base.:(==)(A::BlockMap, B::BlockMap) =
-    (eltype(A) == eltype(B) && A.maps == B.maps && A.rows == B.rows)
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps) && all(A.rows .== B.rows))
 
 ############
 # multiplication helper functions
@@ -350,9 +350,9 @@ function _transblockmul!(y, A::BlockMap, x, α, β, transform)
         # subsequent block rows of A (block columns of A'),
         # add results to corresponding parts of y
         # TODO: think about multithreading
-        for (row, xi) in zip(Base.tail(rows), Base.tail(xinds))
-            xrow = selectdim(x, 1, xi)
-            for _ in 1:row
+        @inbounds for i in 2:length(rows)
+            xrow = selectdim(x, 1, xinds[i])
+            for _ in 1:rows[i]
                 mapind +=1
                 yrow = selectdim(y, 1, yinds[mapind])
                 _unsafe_mul!(yrow, transform(maps[mapind]), xrow, α, true)
@@ -396,24 +396,22 @@ end
 ############
 # BlockDiagonalMap
 ############
-struct BlockDiagonalMap{T,
-                        As<:LinearMapTuple,
-                        Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T}
+struct BlockDiagonalMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T}
     maps::As
-    rowranges::Ranges
-    colranges::Ranges
-    function BlockDiagonalMap{T, As}(maps::As) where {T, As<:LinearMapTuple}
+    rowranges::Vector{UnitRange{Int}}
+    colranges::Vector{UnitRange{Int}}
+    function BlockDiagonalMap{T, As}(maps::As) where {T, As<:LinearMapTupleOrVector}
         for TA in Base.Generator(eltype, maps)
             promote_type(T, TA) == T ||
                 error("eltype $TA cannot be promoted to $T in BlockDiagonalMap constructor")
         end
         rowranges = _getranges(maps, 1)
         colranges = _getranges(maps, 2)
-        return new{T, As, typeof(rowranges)}(maps, rowranges, colranges)
+        return new{T, As}(maps, rowranges, colranges)
     end
 end
 
-BlockDiagonalMap{T}(maps::As) where {T, As<:LinearMapTuple} =
+BlockDiagonalMap{T}(maps::As) where {T, As<:LinearMapTupleOrVector} =
     BlockDiagonalMap{T,As}(maps)
 BlockDiagonalMap(maps::LinearMap...) =
     BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps)
@@ -478,7 +476,7 @@ LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} =
     BlockDiagonalMap{T}(map(transpose, A.maps))
 
 Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) =
-    (eltype(A) == eltype(B) && A.maps == B.maps)
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps))
 
 for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
     @eval begin
@@ -496,7 +494,7 @@ end
 function _blockscaling!(y, A::BlockDiagonalMap, x, α, β)
     maps, yinds, xinds = A.maps, A.rowranges, A.colranges
     # TODO: think about multi-threading here
-    @inbounds for i in eachindex(yinds, maps, xinds)
+    @inbounds for i in 1:length(maps)
         _unsafe_mul!(selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β)
     end
     return y

src/composition.jl

@@ -1,4 +1,4 @@
-struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
+struct CompositeMap{T, As<:LinearMapTupleOrVector} <: LinearMap{T}
     maps::As # stored in order of application to vector
     function CompositeMap{T, As}(maps::As) where {T, As}
         N = length(maps)
@@ -12,7 +12,12 @@ struct CompositeMap{T, As<:LinearMapTuple} <: LinearMap{T}
         new{T, As}(maps)
     end
 end
-CompositeMap{T}(maps::As) where {T, As<:LinearMapTuple} = CompositeMap{T, As}(maps)
+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))
+Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector{<:LinearMap{T}}) where {T} =
+    CompositeMap{T}(reverse(maps))
 
 # basic methods
 Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
@@ -26,8 +31,18 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
         LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
         $_f(maps::Tuple{}) = true
         $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
-        $_f(maps::Tuple{Vararg{LinearMap}}) =
+        $_f(maps::LinearMapTuple) =
             maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
+        function $_f(maps::LinearMapVector)
+            n = length(maps)
+            if n == 0
+                return true
+            elseif n == 1
+                return ($f(maps[1]))::Bool
+            else
+                return ((maps[end] == $g(maps[1]))::Bool && $_f(@views maps[2:end-1]))
+            end
+        end
         # since the introduction of ScaledMap, the following cases cannot occur
         # function $_f(maps::Tuple{Vararg{LinearMap}}) # length(maps) >= 2
             # if maps[1] isa UniformScalingMap{<:RealOrComplex}
@@ -45,23 +60,34 @@ end
 LinearAlgebra.isposdef(A::CompositeMap) = _isposdef(A.maps)
 _isposdef(maps::Tuple{}) = true # empty product is equivalent to "I" which is pos. def.
 _isposdef(maps::Tuple{<:LinearMap}) = isposdef(maps[1])
-function _isposdef(maps::Tuple{Vararg{LinearMap}})
-    (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) && 
-    isposdef(maps[1]) && _isposdef(Base.front(Base.tail(maps)))
+function _isposdef(maps::LinearMapTuple)
+    (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) &&
+        isposdef(maps[1]) && _isposdef(Base.front(Base.tail(maps)))
+end
+function _isposdef(maps::LinearMapVector)
+    n = length(maps)
+    if n == 0
+        return true
+    elseif n == 1
+        return isposdef(maps[1])
+    else
+        return (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) &&
+            isposdef(maps[1]) && _isposdef(maps[2:end-1])
+    end
 end
 
 # scalar multiplication and division (non-commutative case)
 function Base.:(*)(α::Number, A::LinearMap)
     T = promote_type(typeof(α), eltype(A))
-    return CompositeMap{T}(tuple(A, UniformScalingMap(α, size(A, 1))))
+    return CompositeMap{T}(_combine(A, UniformScalingMap(α, size(A, 1))))
 end
 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)..., α * Alast))
+        return CompositeMap{T}(_combine(_front(A.maps), α * Alast))
     else
-        return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
+        return CompositeMap{T}(_combine(A.maps, UniformScalingMap(α, size(A, 1))))
     end
 end
 # needed for disambiguation
@@ -71,15 +97,15 @@ function Base.:(*)(α::RealOrComplex, A::CompositeMap{<:RealOrComplex})
 end
 function Base.:(*)(A::LinearMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
-    return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
+    return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A))
 end
 function Base.:(*)(A::CompositeMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
     Afirst = first(A.maps)
     if Afirst isa UniformScalingMap
-        return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...))
+        return CompositeMap{T}(_combine(Afirst * α, _tail(A.maps)))
     else
-        return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
+        return CompositeMap{T}(_combine(UniformScalingMap(α, size(A, 2)), A.maps))
     end
 end
 # needed for disambiguation
@@ -110,19 +136,19 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3);
 """
 function Base.:(*)(A₁::LinearMap, A₂::LinearMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂, A₁))
+    return CompositeMap{T}(_combine(A₂, A₁))
 end
 function Base.:(*)(A₁::LinearMap, A₂::CompositeMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂.maps..., A₁))
+    return CompositeMap{T}(_combine(A₂.maps, A₁))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::LinearMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂, A₁.maps...))
+    return CompositeMap{T}(_combine(A₂, A₁.maps))
 end
 function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
     T = promote_type(eltype(A₁), eltype(A₂))
-    return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
+    return CompositeMap{T}(_combine(A₂.maps, A₁.maps))
 end
 # needed for disambiguation
 Base.:(*)(A₁::ScaledMap, A₂::CompositeMap) = A₁.λ * (A₁.lmap * A₂)
@@ -135,7 +161,8 @@ LinearAlgebra.adjoint(A::CompositeMap{T}) where {T} =
     CompositeMap{T}(map(adjoint, reverse(A.maps)))
 
 # comparison of CompositeMap objects
-Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
+Base.:(==)(A::CompositeMap, B::CompositeMap) =
+    (eltype(A) == eltype(B) && all(A.maps .== B.maps))
 
 # multiplication with vectors/matrices
 _unsafe_mul!(y::AbstractVecOrMat, A::CompositeMap, x::AbstractVector) =