Add ScaledMap (#94)

Co-authored-by: Daniel Karrasch <[email protected]>

Author: JeffFessler

Project.toml

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

src/LinearMaps.jl

@@ -16,6 +16,7 @@ end
 abstract type LinearMap{T} end
 
 const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}}
+const RealOrComplex = Union{Real,Complex}
 
 Base.eltype(::LinearMap{T}) where {T} = T
 
@@ -117,6 +118,7 @@ include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby al
 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("scaledmap.jl") # multiply by a (real or complex) scalar
 include("composition.jl") # composition of linear maps
 include("functionmap.jl") # using a function as linear map
 include("blockmap.jl") # block linear maps
@@ -132,10 +134,10 @@ Construct a linear map object, either from an existing `LinearMap` or `AbstractM
 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
+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 by numbers of type `T`.
+whether the action of `f` on a vector will be similar to, e.g., multiplying by numbers of type `T`.
 If not specified, the devault value `T=Float64` will be assumed. Optionally, a corresponding
 function `fc` can be specified that implements the transpose/adjoint of `f`.
 
@@ -150,7 +152,7 @@ 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
+    The default value is guessed by looking at the number of arguments of the first occurrence
     of `f` in the method table.
 """
 LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...)

src/composition.jl

@@ -21,27 +21,28 @@ Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
 Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessary
 
 # the following rules are sufficient but not necessary
-const RealOrComplex = Union{Real,Complex}
 for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
                     (:ishermitian, :_ishermitian, :adjoint),
                     (:isposdef, :_isposdef, :adjoint))
     @eval begin
         LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
         $_f(maps::Tuple{}) = true
         $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
-        function $_f(maps::Tuple{Vararg{<:LinearMap}}) # length(maps) >= 2
-            if maps[1] isa UniformScalingMap{<:RealOrComplex}
-                return $f(maps[1]) && $_f(Base.tail(maps))
-            elseif maps[end] isa UniformScalingMap{<:RealOrComplex}
-                return $f(maps[end]) && $_f(Base.front(maps))
-            else
-                return maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
-            end
-        end
+        $_f(maps::Tuple{Vararg{<:LinearMap}}) = maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
+        # 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}
+            #     return $f(maps[1]) && $_f(Base.tail(maps))
+            # elseif maps[end] isa UniformScalingMap{<:RealOrComplex}
+            #     return $f(maps[end]) && $_f(Base.front(maps))
+            # else
+                # return maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
+            # end
+        # end
     end
 end
 
-# scalar multiplication and division
+# 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))))
@@ -55,6 +56,11 @@ function Base.:(*)(α::Number, A::CompositeMap)
         return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
     end
 end
+# needed for disambiguation
+function Base.:(*)(α::RealOrComplex, A::CompositeMap)
+    T = Base.promote_op(*, typeof(α), eltype(A))
+    return ScaledMap{T}(α, A)
+end
 function Base.:(*)(A::LinearMap, α::Number)
     T = promote_type(typeof(α), eltype(A))
     return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
@@ -68,13 +74,18 @@ function Base.:(*)(A::CompositeMap, α::Number)
         return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
     end
 end
+# needed for disambiguation
+function Base.:(*)(A::CompositeMap, α::RealOrComplex)
+    T = Base.promote_op(*, typeof(α), eltype(A))
+    return ScaledMap{T}(α, A)
+end
+
 Base.:(\)(α::Number, A::LinearMap) = inv(α) * A
 Base.:(/)(A::LinearMap, α::Number) = A * inv(α)
-Base.:(-)(A::LinearMap) = -1 * A
 
 # composition of linear maps
 """
-    A::LinearMap * B::LinearMap
+    *(A::LinearMap, B::LinearMap)
 
 Construct a `CompositeMap <: LinearMap`, a (lazy) representation of the product
 of the two operators. Products of `LinearMap`/`CompositeMap` objects and
@@ -109,8 +120,6 @@ function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap)
     T = promote_type(eltype(A₁), eltype(A₂))
     return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...))
 end
-Base.:(*)(A₁::LinearMap, A₂::UniformScaling) = A₁ * A₂.λ
-Base.:(*)(A₁::UniformScaling, A₂::LinearMap) = A₁.λ * A₂
 
 # special transposition behavior
 LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(transpose, reverse(A.maps)))

src/conversion.jl

@@ -45,6 +45,10 @@ Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
 
 # special cases
 
+# ScaledMap
+Base.Matrix(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) = convert(Matrix, A.λ*A.lmap.lmap)
+SparseArrays.sparse(A::ScaledMap{<:Any,<:Any,<:MatrixMap}) = convert(SparseMatrixCSC, A.λ*A.lmap.lmap)
+
 # UniformScalingMap
 Base.Matrix(J::UniformScalingMap) = Matrix(J.λ*I, size(J))
 Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, size(J, 1)))
@@ -77,14 +81,22 @@ for (TT, T) in ((Type{Matrix}, Matrix), (Type{SparseMatrixCSC}, SparseMatrixCSC)
         B, A = AB.maps
         return convert($T, A.lmap*B.lmap)
     end
-    @eval function Base.convert(::$TT, λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
-        A, J = λA.maps
-        return convert($T, J.λ*A.lmap)
-    end
-    @eval function Base.convert(::$TT, Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
-        J, A = Aλ.maps
-        return convert($T, A.lmap*J.λ)
-    end
+end
+function Base.Matrix(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
+    A, J = λA.maps
+    return convert(Matrix, J.λ*A.lmap)
+end
+function SparseArrays.sparse(λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}})
+    A, J = λA.maps
+    return convert(SparseMatrixCSC, J.λ*A.lmap)
+end
+function Base.Matrix(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
+    J, A = Aλ.maps
+    return convert(Matrix, A.lmap*J.λ)
+end
+function SparseArrays.sparse(Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}})
+    J, A = Aλ.maps
+    return convert(SparseMatrixCSC, A.lmap*J.λ)
 end
 
 # BlockMap & BlockDiagonalMap

src/linearcombination.jl

@@ -23,7 +23,7 @@ LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficien
 
 # adding linear maps
 """
-    A::LinearMap + B::LinearMap
+    +(A::LinearMap, B::LinearMap)
 
 Construct a `LinearCombination <: LinearMap`, a (lazy) representation of the sum
 of the two operators. Sums of `LinearMap`/`LinearCombination` objects and

src/scaledmap.jl

@@ -0,0 +1,68 @@
+struct ScaledMap{T, S<:RealOrComplex, A<:LinearMap} <: LinearMap{T}
+    λ::S
+    lmap::A
+    function ScaledMap{T,S,A}(λ::S, lmap::A) where {T, S <: RealOrComplex, A <: LinearMap}
+        Base.promote_op(*, S, eltype(lmap)) == T || throw(InexactError())
+        new{T,S,A}(λ, lmap)
+    end
+end
+
+# constructor
+ScaledMap{T}(λ::S, lmap::A) where {T,S<:RealOrComplex,A<:LinearMap} =
+    ScaledMap{Base.promote_op(*, S, eltype(lmap)),S,A}(λ, lmap)
+
+# show
+function Base.show(io::IO, A::ScaledMap{T}) where {T}
+    println(io, "LinearMaps.ScaledMap{$T}, scale = $(A.λ)")
+    show(io, A.lmap)
+end
+
+# basic methods
+Base.size(A::ScaledMap) = size(A.lmap)
+Base.isreal(A::ScaledMap) = isreal(A.λ) && isreal(A.lmap)
+LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap)
+LinearAlgebra.ishermitian(A::ScaledMap) = ishermitian(A.lmap)
+LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap)
+
+Base.transpose(A::ScaledMap) = A.λ * transpose(A.lmap)
+Base.adjoint(A::ScaledMap) = conj(A.λ) * adjoint(A.lmap)
+
+# comparison (sufficient, not necessary)
+Base.:(==)(A::ScaledMap, B::ScaledMap) =
+    (eltype(A) == eltype(B) && A.lmap == B.lmap) && A.λ == B.λ
+
+# scalar multiplication and division
+function Base.:(*)(α::RealOrComplex, A::LinearMap)
+    T = Base.promote_op(*, typeof(α), eltype(A))
+    return ScaledMap{T}(α, A)
+end
+function Base.:(*)(A::LinearMap, α::RealOrComplex)
+    T = Base.promote_op(*, typeof(α), eltype(A))
+    return ScaledMap{T}(α, A)
+end
+
+Base.:(*)(α::Number, A::ScaledMap) = (α * A.λ) * A.lmap
+Base.:(*)(A::ScaledMap, α::Number) = A.lmap * (A.λ * α)
+# needed for disambiguation
+Base.:(*)(α::RealOrComplex, A::ScaledMap) = (α * A.λ) * A.lmap
+Base.:(*)(A::ScaledMap, α::RealOrComplex) = (A.λ * α) * A.lmap
+Base.:(-)(A::LinearMap) = -1 * A
+
+# composition (not essential, but might save multiple scaling operations)
+Base.:(*)(A::ScaledMap, B::ScaledMap) = (A.λ * B.λ) * (A.lmap * B.lmap)
+Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
+Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ
+
+# multiplication with vectors
+function A_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector)
+    # no size checking, will be done by map
+    mul!(y, A.lmap, x, A.λ, false)
+end
+
+function LinearAlgebra.mul!(y::AbstractVector, A::ScaledMap, x::AbstractVector, α::Number, β::Number)
+    # no size checking, will be done by map
+    mul!(y, A.lmap, x, A.λ * α, β)
+end
+
+At_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x)
+Ac_mul_B!(y::AbstractVector, A::ScaledMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x)

src/uniformscalingmap.jl

@@ -28,8 +28,13 @@ LinearAlgebra.transpose(A::UniformScalingMap) = A
 LinearAlgebra.adjoint(A::UniformScalingMap)   = UniformScalingMap(conj(A.λ), size(A))
 
 # multiplication with scalar
+Base.:(*)(A::UniformScaling, B::LinearMap) = A.λ * B
+Base.:(*)(A::LinearMap, B::UniformScaling) = A * B.λ
 Base.:(*)(α::Number, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J))
 Base.:(*)(J::UniformScalingMap, α::Number) = UniformScalingMap(J.λ * α, size(J))
+# needed for disambiguation
+Base.:(*)(α::RealOrComplex, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J))
+Base.:(*)(J::UniformScalingMap, α::RealOrComplex) = UniformScalingMap(J.λ * α, size(J))
 
 # multiplication with vector
 Base.:(*)(A::UniformScalingMap, x::AbstractVector) =

test/composition.jl

@@ -1,11 +1,14 @@
 using Test, LinearMaps, LinearAlgebra
 
 @testset "composition" begin
-    F = @inferred LinearMap(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
-    FC = @inferred LinearMap{ComplexF64}(cumsum, y -> reverse(cumsum(reverse(x))), 10; ismutating=false)
+    F = @inferred LinearMap(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
+    FC = @inferred LinearMap{ComplexF64}(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
     FCM = LinearMaps.CompositeMap{ComplexF64}((FC,))
+    L = LowerTriangular(ones(10,10))
     A = 2 * rand(ComplexF64, (10, 10)) .- 1
     B = rand(size(A)...)
+    H = LinearMap(Hermitian(A'A))
+    S = LinearMap(Symmetric(real(A)'real(A)))
     M = @inferred 1 * LinearMap(A)
     N = @inferred LinearMap(B)
     v = rand(ComplexF64, 10)
@@ -14,7 +17,11 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred (F * A) * v == @inferred F * (A * v)
     @test @inferred (A * F) * v == @inferred A * (F * v)
     @test @inferred A * (F * F) * v == @inferred A * (F * (F * v))
-    @test @inferred (F * F) * (F * F) * v == @inferred F * (F * (F * (F * v)))
+    F2 = F*F
+    FC2 = FC*FC
+    F4 = FC2 * F2
+    @test length(F4.maps) == 4
+    @test @inferred F4 * v == @inferred F * (F * (F * (F * v)))
     @test @inferred Matrix(M * transpose(M)) ≈ A * transpose(A)
     @test @inferred !isposdef(M * transpose(M))
     @test @inferred isposdef(M * M')
@@ -23,9 +30,12 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred !issymmetric(M' * M)
     @test @inferred ishermitian(M' * M)
     @test @inferred issymmetric(F'F)
+    @test @inferred issymmetric(F'*S*F)
     @test @inferred ishermitian(F'F)
+    @test @inferred ishermitian(F'*H*F)
     @test @inferred !issymmetric(FC'FC)
     @test @inferred ishermitian(FC'FC)
+    @test @inferred ishermitian(FC'*H*FC)
     @test @inferred isposdef(transpose(F) * F * 3)
     @test @inferred isposdef(transpose(F) * 3 * F)
     @test @inferred !isposdef(-5*transpose(F) * F)

test/conversion.jl

@@ -1,4 +1,4 @@
-using Test, LinearMaps, LinearAlgebra, SparseArrays
+using Test, LinearMaps, LinearAlgebra, SparseArrays, Quaternions
 
 @testset "conversion" begin
     A = 2 * rand(ComplexF64, (20, 10)) .- 1
@@ -31,6 +31,26 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays
     JM = convert(AbstractMatrix, J)
     @test JM == Diagonal(fill(α, 10))
 
+    # ScaledMap
+    q = rand(ComplexF64)
+    A = rand(ComplexF64, 3, 3)
+    @test convert(Matrix, q*LinearMap(A)) ≈ q*A
+    qAs = convert(SparseMatrixCSC, q*LinearMap(A))
+    @test qAs ≈ q*A
+    @test qAs isa SparseMatrixCSC
+
+    # CompositeMap of MatrixMap and UniformScalingMap
+    q = Quaternion(rand(4)...)
+    A = Quaternion.(rand(3,3), rand(3,3), rand(3,3), rand(3,3))
+    for mat in (Matrix, SparseMatrixCSC)
+        Aqmat = convert(mat, LinearMap(A)*q)
+        @test Aqmat ≈ A*q
+        @test Aqmat isa mat
+        qAmat = convert(mat, q*LinearMap(A))
+        @test qAmat ≈ q*A
+        @test qAmat isa mat
+    end
+
     # sparse matrix generation/conversion
     @test sparse(M) == sparse(Array(M))
     @test convert(SparseMatrixCSC, M) == sparse(Array(M))

test/functionmap.jl

@@ -1,4 +1,4 @@
-using Test, LinearMaps, LinearAlgebra
+using Test, LinearMaps, LinearAlgebra, BenchmarkTools
 
 @testset "function maps" begin
     N = 100
@@ -58,7 +58,7 @@ using Test, LinearMaps, LinearAlgebra
     @test_throws ErrorException adjoint(CS!) * v
     @test_throws ErrorException mul!(similar(v), CS!', v)
     @test_throws ErrorException mul!(similar(v), transpose(CS!), v)
-    CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y)), 10; ismutating=true)
+    CS! = LinearMap{ComplexF64}(cumsum!, (y, x) -> (copyto!(y, x); reverse!(y); cumsum!(y, y); reverse!(y)), 10; ismutating=true)
     @inferred adjoint(CS!)
     @test @inferred LinearMaps.ismutating(CS!)
     @test @inferred CS! * v == cv
@@ -67,6 +67,11 @@ using Test, LinearMaps, LinearAlgebra
     @test @inferred CS' * v == reverse!(cumsum(reverse(v)))
     @test @inferred mul!(similar(v), transpose(CS), v) == reverse!(cumsum(reverse(v)))
     @test @inferred mul!(similar(v), adjoint(CS), v) == reverse!(cumsum(reverse(v)))
+    u = similar(v)
+    b = @benchmarkable mul!($u, $(3*CS!), $v)
+    @test run(b, samples=3).allocs == 0
+    b = @benchmarkable mul!($u, $(3*CS!'), $v)
+    @test run(b, samples=3).allocs == 0
 
     # Test fallback methods:
     L = @inferred LinearMap(x -> x, x -> x, 10)