Add EmbeddedMap (#174)

Co-authored-by: Jeff Fessler <fessler@umich.edu>

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.6.2"
+version = "3.7.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/history.md

@@ -1,5 +1,22 @@
 # Version history
 
+## What's new in v3.7
+
+* A new map type called `EmbeddedMap` is introduced. It is a wrapper of a "small" `LinearMap`
+  (or a suitably converted `AbstractVecOrMat`) embedded into a "larger" zero map. Hence,
+  the "small" map acts only on a subset of the coordinates and maps to another subset of
+  the coordinates of the "large" map. The large map `L` can therefore be imagined as the
+  composition of a sampling/projection map `P`, of the small map `A`, and of an embedding
+  map `E`: `L = E ⋅ A ⋅ P`. It is implemented, however, by acting on a view of the vector
+  `x` and storing the result into a view of the result vector `y`. Such maps can be
+  constructed by the new methods:
+  * `LinearMap(A::MapOrVecOrMat, dims::Dims{2}, index::NTuple{2, AbstractVector{Int}})`,
+    where `dims` is the dimensions of the "large" map and index is a tuple of the `x`- and
+    `y`-indices that interact with `A`, respectively;
+  * `LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})`, where the keyword
+    argument `offset` determines the dimension of a virtual upper-left zero block, to which
+    `A` gets (virtually) diagonally appended.
+
 ## What's new in v3.6
 
 * Support for Julia versions below v1.6 has been dropped.

src/LinearMaps.jl

@@ -266,6 +266,7 @@ 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("embeddedmap.jl") # embedded linear maps
 include("conversion.jl") # conversion of linear maps to matrices
 include("show.jl") # show methods for LinearMap objects
 
@@ -274,12 +275,19 @@ include("show.jl") # show methods for LinearMap objects
     LinearMap(A::AbstractVecOrMat; kwargs...)::WrappedMap
     LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
+    LinearMap(A::MapOrVecOrMat, dims::Dims{2}, index::NTuple{2, AbstractVector{Int}})::EmbeddedMap
+    LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})::EmbeddedMap
 
-Construct a linear map object, either from an existing `LinearMap` or `AbstractVecOrMat` `A`,
-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
+Construct a linear map object, either
+
+1. from an existing `LinearMap` or `AbstractVecOrMat` `A`, with the purpose of redefining
+  its properties via the keyword arguments `kwargs`;
+2. a `UniformScaling` object `J` with specified (square) dimension `M`;
+3. from a function or callable object `f`;
+4. from an existing `LinearMap` or `AbstractVecOrMat` `A`, embedded in a larger zero map.
+
+In the case of item 3, 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, also the `eltype` `T` of the corresponding matrix representation needs to be
 specified, i.e., whether the action of `f` on a vector will be similar to, e.g., multiplying
@@ -305,6 +313,13 @@ For the function-based constructor, there is one more keyword argument:
     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
     occurrence of `f` in the method table.
+
+For the `EmbeddedMap` constructors, `dims` specifies the total dimensions of the map. The
+`index` argument specifies two collections of indices `inds1` and `inds2`, such that for
+the big zero map `L` (thought of as a matrix), one has `L[inds1,inds2] == A`. In other
+words, `inds1` specifies the output indices, `inds2` specifies the input indices.
+Alternatively, `A` may be shifted by `offset`, such that (thinking in terms of matrices
+again) `L[offset[1] .+ axes(A, 1), offset[2] .+ axes(A, 2)] == A`.
 """
 LinearMap(A::MapOrVecOrMat; kwargs...) = WrappedMap(A; kwargs...)
 LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M)
@@ -313,7 +328,12 @@ 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::MapOrMatrix; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
+LinearMap(A::MapOrVecOrMat, dims::Dims{2}, index::NTuple{2, AbstractVector{Int}}) =
+    EmbeddedMap(convert(LinearMap, A), dims, index[1], index[2])
+LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2}) =
+    EmbeddedMap(convert(LinearMap, A), dims; offset=offset)
+
+LinearMap{T}(A::MapOrVecOrMat; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
 LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)
 
 end # module

src/embeddedmap.jl

@@ -0,0 +1,83 @@
+struct EmbeddedMap{T, As <: LinearMap, Rs <: AbstractVector{Int},
+                    Cs <: AbstractVector{Int}} <: LinearMap{T}
+    lmap::As
+    dims::Dims{2}
+    rows::Rs # typically i1:i2 with 1 <= i1 <= i2 <= size(map,1)
+    cols::Cs # typically j1:j2 with 1 <= j1 <= j2 <= size(map,2)
+
+    function EmbeddedMap{T}(map::As, dims::Dims{2}, rows::Rs, cols::Cs) where {T,
+                    As <: LinearMap, Rs <: AbstractVector{Int}, Cs <: AbstractVector{Int}}
+        check_index(rows, size(map, 1), dims[1])
+        check_index(cols, size(map, 2), dims[2])
+        return new{T,As,Rs,Cs}(map, dims, rows, cols)
+    end
+end
+
+EmbeddedMap(map::LinearMap{T}, dims::Dims{2}; offset::Dims{2}) where {T} =
+    EmbeddedMap{T}(map, dims, offset[1] .+ (1:size(map, 1)), offset[2] .+ (1:size(map, 2)))
+EmbeddedMap(map::LinearMap, dims::Dims{2}, rows::AbstractVector{Int}, cols::AbstractVector{Int}) =
+    EmbeddedMap{eltype(map)}(map, dims, rows, cols)
+
+@static if VERSION >= v"1.8-"
+    Base.reverse(A::LinearMap; dims=:) = _reverse(A, dims)
+    function _reverse(A, dims::Integer)
+        if dims == 1
+            return EmbeddedMap(A, size(A), reverse(axes(A, 1)), axes(A, 2))
+        elseif dims == 2
+            return EmbeddedMap(A, size(A), axes(A, 1), reverse(axes(A, 2)))
+        else
+            throw(ArgumentError("invalid dims argument to reverse, should be 1 or 2, got $dims"))
+        end
+    end
+    _reverse(A, ::Colon) = EmbeddedMap(A, size(A), map(reverse, axes(A))...)
+    _reverse(A, dims::NTuple{1,Integer}) = _reverse(A, first(dims))
+    function _reverse(A, dims::NTuple{M,Integer}) where {M}
+        dimrev = ntuple(k -> k in dims, 2)
+        if 2 < M || M != sum(dimrev)
+            throw(ArgumentError("invalid dimensions $dims in reverse!"))
+        end
+        ax = ntuple(k -> dimrev[k] ? reverse(axes(A, k)) : axes(A, k), 2)
+        return EmbeddedMap(A, size(A), ax...)
+    end
+end
+
+function check_index(index::AbstractVector{Int}, dimA::Int, dimB::Int)
+    length(index) != dimA && throw(ArgumentError("invalid length of index vector"))
+    minimum(index) <= 0 && throw(ArgumentError("minimal index is below 1"))
+    maximum(index) > dimB && throw(ArgumentError(
+        "maximal index $(maximum(index)) exceeds dimension $dimB"
+        ))
+    # _isvalidstep(index) || throw(ArgumentError("non-monotone index set"))
+    nothing
+end
+
+# _isvalidstep(index::AbstractRange) = step(index) > 0
+# _isvalidstep(index::AbstractVector) = all(diff(index) .> 0)
+
+Base.size(A::EmbeddedMap) = A.dims
+
+# sufficient but not necessary conditions
+LinearAlgebra.issymmetric(A::EmbeddedMap) =
+    issymmetric(A.lmap) && (A.dims[1] == A.dims[2]) && (A.rows == A.cols)
+LinearAlgebra.ishermitian(A::EmbeddedMap) =
+    ishermitian(A.lmap) && (A.dims[1] == A.dims[2]) && (A.rows == A.cols)
+
+Base.:(==)(A::EmbeddedMap, B::EmbeddedMap) =
+    (eltype(A) == eltype(B)) && (A.lmap == B.lmap) &&
+    (A.dims == B.dims) && (A.rows == B.rows) && (A.cols == B.cols)
+
+LinearAlgebra.adjoint(A::EmbeddedMap) = EmbeddedMap(adjoint(A.lmap), reverse(A.dims), A.cols, A.rows)
+LinearAlgebra.transpose(A::EmbeddedMap) = EmbeddedMap(transpose(A.lmap), reverse(A.dims), A.cols, A.rows)
+
+for (In, Out) in ((AbstractVector, AbstractVecOrMat), (AbstractMatrix, AbstractMatrix))
+    @eval function _unsafe_mul!(y::$Out, A::EmbeddedMap, x::$In)
+        fill!(y, zero(eltype(y)))
+        _unsafe_mul!(selectdim(y, 1, A.rows), A.lmap, selectdim(x, 1, A.cols))
+        return y
+    end
+    @eval function _unsafe_mul!(y::$Out, A::EmbeddedMap, x::$In, alpha::Number, beta::Number)
+        LinearAlgebra._rmul_or_fill!(y, beta)
+        _unsafe_mul!(selectdim(y, 1, A.rows), A.lmap, selectdim(x, 1, A.cols), alpha, !iszero(beta))
+        return y
+    end
+end

test/embeddedmap.jl

@@ -0,0 +1,57 @@
+using Test, LinearMaps, LinearAlgebra, SparseArrays
+
+@testset "embeddedmap" begin
+    m = 6; n = 5
+    M = 10(1:m) .+ (1:n)'; L = LinearMap(M)
+    offset = (3,4)
+
+    BM = [zeros(offset...) zeros(offset[1], size(M,2));
+        zeros(size(M,1), offset[2]) M]
+    BL = @inferred LinearMap(L, size(BM); offset=offset)
+    s1, s2 = size(BM)
+    @test (@inferred Matrix(BL)) == BM
+    @test (@inferred Matrix(BL')) == BM'
+    @test (@inferred Matrix(transpose(BL))) == transpose(BM)
+
+    @test_throws UndefKeywordError LinearMap(M, (10, 10))
+    @test_throws ArgumentError LinearMap(M, (m, n), (0:m, 1:n))
+    @test_throws ArgumentError LinearMap(M, (m, n), (0:m-1, 1:n))
+    @test_throws ArgumentError LinearMap(M, (m, n), (1:m, 1:n+1))
+    @test_throws ArgumentError LinearMap(M, (m, n), (1:m, 2:n+1))
+    @test_throws ArgumentError LinearMap(M, (m, n), offset=(3,3))
+    # @test_throws ArgumentError LinearMap(M, (m, n), (m:-1:1, 1:n))
+    # @test_throws ArgumentError LinearMap(M, (m, n), (collect(m:-1:1), 1:n))
+    @test size(@inferred LinearMap(M, (2m, 2n), (1:2:2m, 1:2:2n))) == (2m, 2n)
+    @test @inferred !ishermitian(BL)
+    @test @inferred !issymmetric(BL)
+    @test @inferred LinearMap(L, size(BM), (offset[1] .+ (1:m), offset[2] .+ (1:n))) == BL
+    Wc = @inferred LinearMap([2 im; -im 0]; ishermitian=true)
+    Bc = @inferred LinearMap(Wc, (4,4); offset=(2,2))
+    @test (@inferred ishermitian(Bc))
+
+    x = randn(s2); X = rand(s2, 3)
+    y = BM * x; Y = zeros(s1, 3)
+
+    @test @inferred BL * x ≈ BM * x
+    @test @inferred BL' * y ≈ BM' * y
+
+    for α in (true, false, rand()),
+            β in (true, false, rand()),
+            t in (identity, adjoint, transpose)
+
+        @test t(BL) * x ≈ mul!(copy(y), t(BL), x) ≈ t(BM) * x
+        @test Matrix(t(BL) * X) ≈ mul!(copy(Y), t(BL), X) ≈ t(BM) * X
+        y = randn(s1); Y = randn(s1, 3)
+        @test (@inferred mul!(copy(y), t(BL), x, α, β)) ≈ mul!(copy(y), t(BM), x, α, β)
+        @test (@inferred mul!(copy(Y), t(BL), X, α, β)) ≈ mul!(copy(Y), t(BM), X, α, β)
+    end
+
+    if VERSION >= v"1.8"
+        M = rand(3,4)
+        L = LinearMap(M)
+        @test Matrix(reverse(L)) == reverse(M)
+        for dims in (1, 2, (1,), (2,), (1, 2), (2, 1), :)
+            @test Matrix(reverse(L, dims=dims)) == reverse(M, dims=dims)
+        end
+    end
+end

test/runtests.jl

@@ -33,3 +33,5 @@ include("left.jl")
 include("fillmap.jl")
 
 include("nontradaxes.jl")
+
+include("embeddedmap.jl")