Add FillMaps (#113)

Author: dkarrasch

docs/src/index.md

@@ -16,23 +16,24 @@ Let
 
     A = LinearMap(rand(10, 10))
     B = LinearMap(cumsum, reverse∘cumsum∘reverse, 10)
-    
+
 be a matrix- and function-based linear map, respectively. Then the following code just works,
 indistinguishably from the case when `A` and `B` are both `AbstractMatrix`-typed objects.
 
-```
-3.0A + 2B
-A*B'
-[A B; B A]
-kron(A, B)
-```
+    3.0A + 2B
+    A + I
+    A*B'
+    [A B; B A]
+    kron(A, B)
 
 The `LinearMap` type and corresponding methods combine well with the following packages:
+
 * [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl): iterative eigensolver
   `eigs` and SVD `svds`;
 * [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl): iterative
   solvers, eigensolvers, and SVD;
-* [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl): Krylov-based algorithms for linear problems, singular value and eigenvalue problems
+* [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl): Krylov-based algorithms for linear
+  problems, singular value and eigenvalue problems
 * [TSVD.jl](https://github.com/andreasnoack/TSVD.jl): truncated SVD `tsvd`.
 
 ```julia
@@ -95,34 +96,34 @@ operator in the special case of a square matrix).
 
 The LinearMaps package provides the following functionality:
 
-1.  A `LinearMap` type that shares with the `AbstractMatrix` type that it
-    responds to the functions `size`, `eltype`, `isreal`, `issymmetric`,
-    `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication
-    with a vector using both `*` or the in-place version `mul!`. Linear algebra
-    functions that use duck-typing for their arguments can handle `LinearMap`
-    objects similar to `AbstractMatrix` objects, provided that they can be
-    written using the above methods. Unlike `AbstractMatrix` types, `LinearMap`
-    objects cannot be indexed, neither using `getindex` or `setindex!`.
-
-2.  A single function `LinearMap` that acts as a general purpose
-    constructor (though it is only an abstract type) and allows to construct
-    linear map objects from functions, or to wrap objects of type
-    `AbstractMatrix` or `LinearMap`. The latter functionality is useful to
-    (re)define the properties (`isreal`, `issymmetric`, `ishermitian`,
-    `isposdef`) of the existing matrix or linear map.
-
-3.  A framework for combining objects of type `LinearMap` and of type
-    `AbstractMatrix` using linear combinations, transposition, composition,
-    concatenation and Kronecker product/sums,
-    where the linear map resulting from these operations is never explicitly
-    evaluated but only its matrix-vector product is defined (i.e. lazy
-    evaluation). The matrix-vector product is written to minimize memory
-    allocation by using a minimal number of temporary vectors. There is full
-    support for the in-place version `mul!`, which should be preferred for
-    higher efficiency in critical algorithms. In addition, it tries to recognize
-    the properties of combinations of linear maps. In particular, compositions
-    such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A`
-    and `B` and positive definite `C` are recognized as being positive definite
-    and hermitian. In case a certain property of the resulting `LinearMap`
-    object is not correctly inferred, the `LinearMap` method can be called to
-    redefine the properties.
+1. A `LinearMap` type that shares with the `AbstractMatrix` type that it
+   responds to the functions `size`, `eltype`, `isreal`, `issymmetric`,
+   `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication
+   with a vector using both `*` or the in-place version `mul!`. Linear algebra
+   functions that use duck-typing for their arguments can handle `LinearMap`
+   objects similar to `AbstractMatrix` objects, provided that they can be
+   written using the above methods. Unlike `AbstractMatrix` types, `LinearMap`
+   objects cannot be indexed, neither using `getindex` or `setindex!`.
+
+2. A single function `LinearMap` that acts as a general purpose
+   constructor (though it is only an abstract type) and allows to construct
+   linear map objects from functions, or to wrap objects of type
+   `AbstractMatrix` or `LinearMap`. The latter functionality is useful to
+   (re)define the properties (`isreal`, `issymmetric`, `ishermitian`,
+   `isposdef`) of the existing matrix or linear map.
+
+3. A framework for combining objects of type `LinearMap` and of type
+   `AbstractMatrix` using linear combinations, transposition, composition,
+   concatenation and Kronecker product/sums,
+   where the linear map resulting from these operations is never explicitly
+   evaluated but only its matrix-vector product is defined (i.e. lazy
+   evaluation). The matrix-vector product is written to minimize memory
+   allocation by using a minimal number of temporary vectors. There is full
+   support for the in-place version `mul!`, which should be preferred for
+   higher efficiency in critical algorithms. In addition, it tries to recognize
+   the properties of combinations of linear maps. In particular, compositions
+   such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A`
+   and `B` and positive definite `C` are recognized as being positive definite
+   and hermitian. In case a certain property of the resulting `LinearMap`
+   object is not correctly inferred, the `LinearMap` method can be called to
+   redefine the properties.

docs/src/related.md

@@ -3,30 +3,25 @@
 The following open-source packages provide similar or even extended functionality as
 `LinearMaps.jl`.
 
-*   [`Spot`: A linear-operator toolbox for Matlab](https://github.com/mpf/spot),
-    which seems to have heavily inspired the Julia package
-    [`LinearOperators.jl`](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl)
-    and the Python package [`PyLops`](https://github.com/equinor/pylops)
-
-*   [`fastmat`: fast linear transforms in Python](https://pypi.org/project/fastmat/)
-
-*   [`FunctionOperators.jl`](https://github.com/hakkelt/FunctionOperators.jl)
-    and [`LinearMapsAA.jl`](https://github.com/JeffFessler/LinearMapsAA.jl)
-    also support mappings between `Array`s, inspired by the `fatrix` object type in the
-    [Matlab version of the Michigan Image Reconstruction Toolbox (MIRT)](https://github.com/JeffFessler/mirt).
+* [`Spot`: A linear-operator toolbox for Matlab](https://github.com/mpf/spot),
+  which seems to have heavily inspired the Julia package
+  [`LinearOperators.jl`](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl)
+  and the Python package [`PyLops`](https://github.com/equinor/pylops)
+* [`fastmat`: fast linear transforms in Python](https://pypi.org/project/fastmat/)
+* [`FunctionOperators.jl`](https://github.com/hakkelt/FunctionOperators.jl)
+  and [`LinearMapsAA.jl`](https://github.com/JeffFessler/LinearMapsAA.jl)
+  also support mappings between `Array`s, inspired by the `fatrix` object type in the
+  [Matlab version of the Michigan Image Reconstruction Toolbox (MIRT)](https://github.com/JeffFessler/mirt).
 
 As for lazy array manipulation (like addition, composition, Kronecker products and concatenation),
 there exist further related packages in the Julia ecosystem:
 
-*   [`LazyArrays.jl`](https://github.com/JuliaArrays/LazyArrays.jl)
-
-*   [`BlockArrays.jl`](https://github.com/JuliaArrays/BlockArrays.jl)
-
-*   [`BlockDiagonals.jl`](https://github.com/invenia/BlockDiagonals.jl)
-
-*   [`Kronecker.jl`](https://github.com/MichielStock/Kronecker.jl)
+* [`LazyArrays.jl`](https://github.com/JuliaArrays/LazyArrays.jl)
+* [`BlockArrays.jl`](https://github.com/JuliaArrays/BlockArrays.jl)
+* [`BlockDiagonals.jl`](https://github.com/invenia/BlockDiagonals.jl)
+* [`Kronecker.jl`](https://github.com/MichielStock/Kronecker.jl)
+* [`FillArrays.jl`](https://github.com/JuliaArrays/FillArrays.jl)
 
-*   [`FillArrays.jl`](https://github.com/JuliaArrays/FillArrays.jl)
 Since these packages provide types that are subtypes of Julia `Base`'s `AbstractMatrix` type,
 objects of those types can be wrapped by a `LinearMap` and freely mixed with, for instance,
 function-based linear maps. The same applies to custom matrix types as provided, for instance,

docs/src/types.md

@@ -84,6 +84,14 @@ Base.cat
 SparseArrays.blockdiag
 ```
 
+### `FillMap`
+
+Type for lazily representing constantly filled matrices.
+
+```@docs
+LinearMaps.FillMap
+```
+
 ## Methods
 
 ### Multiplication methods
@@ -103,28 +111,28 @@ as in the usual matrix case: `transpose(A) * x` and `mul!(y, A', x)`, for instan
 
 ### Conversion methods
 
-*   `Array`, `Matrix` and associated `convert` methods
+* `Array`, `Matrix` and associated `convert` methods
 
-    Create a dense matrix representation of the `LinearMap` object, by
-    multiplying it with the successive basis vectors. This is mostly for testing
-    purposes or if you want to have the explicit matrix representation of a
-    linear map for which you only have a function definition (e.g. to be able to
-    use its `transpose` or `adjoint`). This way, one may conveniently make `A`
-    act on the columns of a matrix `X`, instead of interpreting `A * X` as a
-    composed linear map: `Matrix(A * X)`. For generic code, that is supposed to
-    handle both `A::AbstractMatrix` and `A::LinearMap`, it is recommended to use
-    `convert(Matrix, A*X)`.
+  Create a dense matrix representation of the `LinearMap` object, by
+  multiplying it with the successive basis vectors. This is mostly for testing
+  purposes or if you want to have the explicit matrix representation of a
+  linear map for which you only have a function definition (e.g. to be able to
+  use its `transpose` or `adjoint`). This way, one may conveniently make `A`
+  act on the columns of a matrix `X`, instead of interpreting `A * X` as a
+  composed linear map: `Matrix(A * X)`. For generic code, that is supposed to
+  handle both `A::AbstractMatrix` and `A::LinearMap`, it is recommended to use
+  `convert(Matrix, A*X)`.
 
-*   `convert(Abstract[Matrix/Array], A::LinearMap)`
+* `convert(Abstract[Matrix/Array], A::LinearMap)`
 
-    Create an `AbstractMatrix` representation of the `LinearMap`. This falls
-    back to `Matrix(A)`, but avoids explicit construction in case the `LinearMap`
-    object is matrix-based.
+  Create an `AbstractMatrix` representation of the `LinearMap`. This falls
+  back to `Matrix(A)`, but avoids explicit construction in case the `LinearMap`
+  object is matrix-based.
 
-*   `SparseArrays.sparse(A::LinearMap)` and `convert(SparseMatrixCSC, A::LinearMap)`
+* `SparseArrays.sparse(A::LinearMap)` and `convert(SparseMatrixCSC, A::LinearMap)`
 
-    Create a sparse matrix representation of the `LinearMap` object, by
-    multiplying it with the successive basis vectors. This is mostly for testing
-    purposes or if you want to have the explicit sparse matrix representation of
-    a linear map for which you only have a function definition (e.g. to be able
-    to use its `transpose` or `adjoint`).
+  Create a sparse matrix representation of the `LinearMap` object, by
+  multiplying it with the successive basis vectors. This is mostly for testing
+  purposes or if you want to have the explicit sparse matrix representation of
+  a linear map for which you only have a function definition (e.g. to be able
+  to use its `transpose` or `adjoint`).

src/LinearMaps.jl

@@ -2,6 +2,7 @@ module LinearMaps
 
 export LinearMap
 export ⊗, kronsum, ⊕
+export FillMap
 
 using LinearAlgebra
 import LinearAlgebra: mul!
@@ -239,18 +240,22 @@ 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("fillmap.jl") # linear maps representing constantly filled matrices
 include("conversion.jl") # conversion of linear maps to matrices
 include("show.jl") # show methods for LinearMap objects
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap
     LinearMap(A::AbstractMatrix; kwargs...)::WrappedMap
     LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
+    LinearMap(λ::Number, M::Int, N::Int) = FillMap(λ, (M, N))::FillMap
+    LinearMap(λ::Number, dims::Dims{2}) = FillMap(λ, dims)::FillMap
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
 
 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`;
-a `UniformScaling` object `J` with specified (square) dimension `M`; or
+a `UniformScaling` object `J` with specified (square) dimension `M`; from a `Number`
+object to lazily represent filled matrices; 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`).

src/conversion.jl

@@ -171,3 +171,6 @@ function SparseArrays.sparse(L::KroneckerSumMap)
     IB = sparse(Diagonal(ones(Bool, size(B, 1))))
     return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
 end
+
+# FillMap
+Base.Matrix{T}(A::FillMap) where {T} = fill(T(A.λ), size(A))

src/fillmap.jl

@@ -0,0 +1,65 @@
+"""
+    FillMap(λ, (m, n))::FillMap
+    FillMap(λ, m, n)::FillMap
+
+Construct a (lazy) representation of an operator whose matrix representation
+would be an m×n-matrix filled constantly with the value `λ`.
+"""
+struct FillMap{T} <: LinearMap{T}
+    λ::T
+    size::Dims{2}
+    function FillMap(λ::T, dims::Dims{2}) where {T}
+        (dims[1]>=0 && dims[2]>=0) || throw(ArgumentError("dims of FillMap must be non-negative"))
+        return new{T}(λ, dims)
+    end
+end
+
+FillMap(λ, m::Int, n::Int) = FillMap(λ, (m, n))
+
+# properties
+Base.size(A::FillMap) = A.size
+MulStyle(A::FillMap) = FiveArg()
+LinearAlgebra.issymmetric(A::FillMap) = A.size[1] == A.size[2]
+LinearAlgebra.ishermitian(A::FillMap) = isreal(A.λ) && A.size[1] == A.size[2]
+LinearAlgebra.isposdef(A::FillMap) = (size(A, 1) == size(A, 2) == 1 && isposdef(A.λ))
+Base.:(==)(A::FillMap, B::FillMap) = A.λ == B.λ && A.size == B.size
+
+LinearAlgebra.adjoint(A::FillMap) = FillMap(adjoint(A.λ), reverse(A.size))
+LinearAlgebra.transpose(A::FillMap) = FillMap(transpose(A.λ), reverse(A.size))
+
+function Base.:(*)(A::FillMap, x::AbstractVector)
+    T = typeof(oneunit(eltype(A)) * (zero(eltype(x)) + zero(eltype(x))))
+    return fill(iszero(A.λ) ? zero(T) : A.λ*sum(x), A.size[1])
+end
+
+function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector)
+    return fill!(y, iszero(A.λ) ? zero(eltype(y)) : A.λ*sum(x))
+end
+
+function _unsafe_mul!(y::AbstractVecOrMat, A::FillMap, x::AbstractVector, α::Number, β::Number)
+    if iszero(α)
+        !isone(β) && rmul!(y, β)
+    else
+        temp = A.λ * sum(x) * α
+        if iszero(β)
+            y .= temp
+        elseif isone(β)
+            y .+= temp
+        else
+            y .= y .* β .+ temp
+        end
+    end
+    return y
+end
+
+Base.:(+)(A::FillMap, B::FillMap) = A.size == B.size ? FillMap(A.λ + B.λ, A.size) : throw(DimensionMismatch())
+Base.:(-)(A::FillMap) = FillMap(-A.λ, A.size)
+Base.:(*)(λ::Number, A::FillMap) = FillMap(λ * A.λ, size(A))
+Base.:(*)(A::FillMap, λ::Number) = FillMap(A.λ * λ, size(A))
+Base.:(*)(λ::RealOrComplex, A::FillMap) = FillMap(λ * A.λ, size(A))
+Base.:(*)(A::FillMap, λ::RealOrComplex) = FillMap(A.λ * λ, size(A))
+
+function Base.:(*)(A::FillMap, B::FillMap)
+    check_dim_mul(A, B)
+    return FillMap(A.λ*B.λ*size(A, 2), (size(A, 1), size(B, 2)))
+end

src/show.jl

@@ -39,6 +39,9 @@ end
 function _show(io::IO, A::ScaledMap{T}, i) where {T}
     " with scale: $(A.λ) of\n" * map_show(io, A.lmap, i+2)
 end
+function _show(io::IO, A::FillMap{T}, _) where {T}
+    " with fill value: $(A.λ)"
+end
 
 # helper functions
 function _show_typeof(A::LinearMap{T}) where {T}

test/fillmap.jl

@@ -0,0 +1,40 @@
+using LinearMaps, LinearAlgebra, Test
+
+@testset "filled maps" begin
+    M, N = 2, 3
+    μ = rand()
+    for λ in (true, false, 3, μ, μ + 2im)
+        L = FillMap(λ, (M, N))
+        @test L == FillMap(λ, M, N)
+        @test occursin("$M×$N FillMap{$(typeof(λ))} with fill value: $λ", sprint((t, s) -> show(t, "text/plain", s), L))
+        @test LinearMaps.MulStyle(L) === LinearMaps.FiveArg()
+        A = fill(λ, (M, N))
+        x = rand(typeof(λ) <: Real ? Float64 : ComplexF64, 3)
+        X = rand(typeof(λ) <: Real ? Float64 : ComplexF64, 3, 4)
+        w = similar(x, 2)
+        W = similar(X, 2, 4)
+        @test size(L) == (M, N)
+        @test adjoint(L) == FillMap(adjoint(λ), (3,2))
+        @test transpose(L) == FillMap(λ, (3,2))
+        @test Matrix(L) == A
+        @test L * x ≈ A * x
+        @test mul!(w, L, x) ≈ A * x
+        @test mul!(W, L, X) ≈ A * X
+        for α in (true, false, 1, 0, randn()), β in (true, false, 1, 0, randn())
+            @test mul!(copy(w), L, x, α, β) ≈ fill(λ * sum(x) * α, M) + w * β
+            @test mul!(copy(W), L, X, α, β) ≈ λ * reduce(vcat, sum(X, dims=1) for _ in 1:2) * α + W * β
+        end
+    end
+    @test issymmetric(FillMap(μ + 1im, (3, 3)))
+    @test ishermitian(FillMap(μ + 0im, (3, 3)))
+    @test isposdef(FillMap(μ, (1,1))) == isposdef(μ)
+    @test !isposdef(FillMap(μ, (3,3)))
+    α = rand()
+    β = rand()
+    @test FillMap(μ, (M, N)) + FillMap(α, (M, N)) == FillMap(μ + α, (M, N))
+    @test FillMap(μ, (M, N)) - FillMap(α, (M, N)) == FillMap(μ - α, (M, N))
+    @test α*FillMap(μ, (M, N)) == FillMap(α * μ, (M, N))
+    @test FillMap(μ, (M, N))*α == FillMap(μ * α, (M, N))
+    @test FillMap(μ, (M, N))*FillMap(μ, (N, M)) == FillMap(μ^2*N, (M, M))
+    @test Matrix(FillMap(μ, (M, N))*FillMap(μ, (N, M))) ≈ fill(μ, (M, N))*fill(μ, (N, M))
+end

test/numbertypes.jl

@@ -45,6 +45,7 @@ using Test, LinearMaps, LinearAlgebra, Quaternions
     @test M == (L * L * γ) * β == (L * L * α) * β == (L * L * α) * β.λ
     @test length(M.maps) == 3
     @test M.maps[1].λ == γ*β.λ
+    @test γ*FillMap(γ, (3, 4)) == FillMap(γ^2, (3, 4)) == FillMap(γ, (3, 4))*γ
 
     # exercise non-RealOrComplex scalar operations
     @test Array(γ * (L'*L)) ≈ γ * (A'*A) # CompositeMap

test/runtests.jl

@@ -30,3 +30,5 @@ include("kronecker.jl")
 include("conversion.jl")
 
 include("left.jl")
+
+include("fillmap.jl")