add support for noncommutative number types (#37)

* add support for noncommutative number types

* delete identitymap.jl

* add comparison of CompositeMap objects

* add (noncommutative) quaternion tests

* improve logic of issymmetric & friends

* add 138 inference tests

Author: dkarrasch

src/LinearMaps.jl

@@ -105,7 +105,7 @@ include("transpose.jl") # transposing linear maps
 include("linearcombination.jl") # defining linear combinations of linear maps
 include("composition.jl") # composition of linear maps
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
-include("identitymap.jl") # the identity map, to be able to make linear combinations of LinearMap objects and I
+include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
 include("functionmap.jl") # using a function as linear map
 
 """

src/composition.jl

@@ -18,36 +18,56 @@ 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
-function LinearAlgebra.issymmetric(A::CompositeMap)
-    N = length(A.maps)
-    if isodd(N)
-        issymmetric(A.maps[div(N+1, 2)]) || return false
-    end
-    for n = 1:div(N, 2)
-        A.maps[n] == transpose(A.maps[N-n+1]) || return false
+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
     end
-    return true
 end
-function LinearAlgebra.ishermitian(A::CompositeMap)
-    N = length(A.maps)
-    if isodd(N)
-        ishermitian(A.maps[div(N+1, 2)]) || return false
-    end
-    for n = 1:div(N, 2)
-        A.maps[n] == adjoint(A.maps[N-n+1]) || return false
-    end
-    return true
+
+# scalar multiplication and division
+function Base.:(*)(α::Number, A::LinearMap)
+    T = promote_type(eltype(α), eltype(A))
+    return CompositeMap{T}(tuple(A, UniformScalingMap(α, size(A, 1))))
 end
-function LinearAlgebra.isposdef(A::CompositeMap)
-    N = length(A.maps)
-    if isodd(N)
-        isposdef(A.maps[div(N+1, 2)]) || return false
+function Base.:(*)(α::Number, A::CompositeMap)
+    T = promote_type(eltype(α), eltype(A))
+    Alast = last(A.maps)
+    if Alast isa UniformScalingMap
+        return CompositeMap{T}(tuple(Base.front(A.maps)..., UniformScalingMap(α * Alast.λ, size(Alast, 1))))
+    else
+        return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1))))
     end
-    for n = 1:div(N, 2)
-        A.maps[n] == adjoint(A.maps[N-n+1]) || return false
+end
+function Base.:(*)(A::LinearMap, α::Number)
+    T = promote_type(eltype(α), eltype(A))
+    return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A))
+end
+function Base.:(*)(A::CompositeMap, α::Number)
+    T = promote_type(eltype(α), eltype(A))
+    Afirst = first(A.maps)
+    if Afirst isa UniformScalingMap
+        return CompositeMap{T}(tuple(UniformScalingMap(Afirst.λ * α, size(Afirst, 1)), Base.tail(A.maps)...))
+    else
+        return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...))
     end
-    return true
 end
+Base.:(\)(α::Number, A::LinearMap) = inv(α) * A
+Base.:(/)(A::LinearMap, α::Number) = A * inv(α)
+Base.:(-)(A::LinearMap) = -1 * A
 
 # composition of linear maps
 function Base.:(*)(A1::CompositeMap, A2::CompositeMap)
@@ -70,11 +90,16 @@ function Base.:(*)(A1::LinearMap, A2::LinearMap)
     T = promote_type(eltype(A1),eltype(A2))
     return CompositeMap{T}(tuple(A2, A1))
 end
+Base.:(*)(A1::LinearMap, A2::UniformScaling{T}) where {T} = A1 * A2[1,1]
+Base.:(*)(A1::UniformScaling{T}, A2::LinearMap) where {T} = A1[1,1] * A2
 
 # special transposition behavior
 LinearAlgebra.transpose(A::CompositeMap{T}) where {T} = CompositeMap{T}(map(transpose, reverse(A.maps)))
 LinearAlgebra.adjoint(A::CompositeMap{T}) where {T}   = CompositeMap{T}(map(adjoint, reverse(A.maps)))
 
+# comparison of LinearCombination objects
+Base.:(==)(A::CompositeMap, B::CompositeMap) = (eltype(A) == eltype(B) && A.maps == B.maps)
+
 # multiplication with vectors
 function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
     # no size checking, will be done by individual maps

src/identitymap.jl

@@ -1,136 +1,86 @@
-struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}, Ts<:Tuple} <: LinearMap{T}
+struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     maps::As
-    coeffs::Ts
-    function LinearCombination{T, As, Ts}(maps::As, coeffs::Ts) where {T, As, Ts}
+    function LinearCombination{T, As}(maps::As) where {T, As}
         N = length(maps)
-        N == length(coeffs) || error("number of coefficients doesn't match number of linear maps")
         sz = size(maps[1])
         for n = 1:N
             size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination"))
-            promote_type(T, eltype(maps[n]), typeof(coeffs[n])) == T || throw(InexactError())
+            promote_type(T, eltype(maps[n])) == T || throw(InexactError())
         end
-        new{T, As, Ts}(maps, coeffs)
+        new{T, As}(maps)
     end
 end
 
-LinearCombination{T}(maps::As, coeffs::Ts) where {T, As, Ts} = LinearCombination{T, As, Ts}(maps, coeffs)
+LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps)
 
 # basic methods
 Base.size(A::LinearCombination) = size(A.maps[1])
 LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) # sufficient but not necessary
-LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) && all(isreal, A.coeffs) # sufficient but not necessary
-LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) && all(isposdef, A.coeffs) # sufficient but not necessary
+LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) # sufficient but not necessary
+LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficient but not necessary
 
 # adding linear maps
 function Base.:(+)(A1::LinearCombination, A2::LinearCombination)
     size(A1) == size(A2) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1.maps..., A2.maps...), tuple(A1.coeffs..., A2.coeffs...))
+    return LinearCombination{T}(tuple(A1.maps..., A2.maps...))
 end
 function Base.:(+)(A1::LinearMap, A2::LinearCombination)
     size(A1) == size(A2) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1, A2.maps...), tuple(one(T), A2.coeffs...))
+    return LinearCombination{T}(tuple(A1, A2.maps...))
 end
 Base.:(+)(A1::LinearCombination, A2::LinearMap) = +(A2, A1)
 function Base.:(+)(A1::LinearMap, A2::LinearMap)
     size(A1)==size(A2) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1, A2), tuple(one(T), one(T)))
+    return LinearCombination{T}(tuple(A1, A2))
 end
-function Base.:(-)(A1::LinearCombination, A2::LinearCombination)
-    size(A1) == size(A2) || throw(DimensionMismatch("-"))
-    T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1.maps..., A2.maps...), tuple(A1.coeffs..., map(-, A2.coeffs)...))
-end
-function Base.:(-)(A1::LinearMap, A2::LinearCombination)
-    size(A1) == size(A2) || throw(DimensionMismatch("-"))
-    T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1, A2.maps...), tuple(one(T), map(-, A2.coeffs)...))
-end
-function Base.:(-)(A1::LinearCombination, A2::LinearMap)
-    size(A1) == size(A2) || throw(DimensionMismatch("-"))
-    T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1.maps..., A2), tuple(A1.coeffs..., -one(T)))
-end
-function Base.:(-)(A1::LinearMap, A2::LinearMap)
-    size(A1) == size(A2) || throw(DimensionMismatch("-"))
-    T = promote_type(eltype(A1), eltype(A2))
-    return LinearCombination{T}(tuple(A1, A2), tuple(one(T), -one(T)))
-end
-
-# scalar multiplication
-Base.:(-)(A::LinearMap) = LinearCombination{eltype(A)}(tuple(A), tuple(-one(eltype(A))))
-Base.:(-)(A::LinearCombination) = LinearCombination{eltype(A)}(A.maps, map(-, A.coeffs))
-
-function Base.:(*)(α::Number, A::LinearMap)
-    T = promote_type(eltype(α), eltype(A))
-    return LinearCombination{T}(tuple(A), tuple(α))
-end
-Base.:(*)(A::LinearMap, α::Number) = *(α, A)
-function Base.:(*)(α::Number, A::LinearCombination)
-    T = promote_type(eltype(α), eltype(A))
-    return LinearCombination{T}(A.maps, map(x->α*x, A.coeffs))
-end
-Base.:(*)(A::LinearCombination, α::Number) = *(α, A)
-
-function Base.:(\)(α::Number, A::LinearMap)
-    T = promote_type(eltype(α), eltype(A))
-    return LinearCombination{T}(tuple(A), tuple(1/α))
-end
-Base.:(/)(A::LinearMap, α::Number) = \(α, A)
-function Base.:(\)(α::Number, A::LinearCombination)
-    T = promote_type(eltype(α), eltype(A))
-    return LinearCombination{T}(A.maps, map(x->α\x, A.coeffs))
-end
-Base.:(/)(A::LinearCombination, α::Number) = \(α, A)
+Base.:(-)(A1::LinearMap, A2::LinearMap) = +(A1, -A2)
 
-# comparison of LinearCombination objects
-Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A)==eltype(B) && A.maps==B.maps && A.coeffs==B.coeffs)
+# comparison of LinearCombination objects, sufficient but not necessary
+Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps)
 
 # special transposition behavior
-LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps), A.coeffs)
-LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(A)}(map(adjoint, A.maps), map(conj, A.coeffs))
+LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps))
+LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(A)}(map(adjoint, A.maps))
 
 # multiplication with vectors
 function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
     # no size checking, will be done by individual maps
     A_mul_B!(y, A.maps[1], x)
-    A.coeffs[1] == 1 || lmul!(A.coeffs[1], y)
     l = length(A.maps)
     if l>1
         z = similar(y)
         for n=2:l
             A_mul_B!(z, A.maps[n], x)
-            axpy!(A.coeffs[n], z, y)
+            y .+= z
         end
     end
     return y
 end
 function At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
     # no size checking, will be done by individual maps
     At_mul_B!(y, A.maps[1], x)
-    A.coeffs[1] == 1 || lmul!(A.coeffs[1], y)
     l = length(A.maps)
     if l>1
         z = similar(y)
         for n = 2:l
             At_mul_B!(z, A.maps[n], x)
-            axpy!(A.coeffs[n], z, y)
+            y .+= z
         end
     end
     return y
 end
 function Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
     # no size checking, will be done by individual maps
     Ac_mul_B!(y, A.maps[1], x)
-    A.coeffs[1] == 1 || lmul!(conj(A.coeffs[1]), y)
     l = length(A.maps)
     if l>1
         z = similar(y)
         for n=2:l
             Ac_mul_B!(z, A.maps[n], x)
-            axpy!(conj(A.coeffs[n]), z, y)
+            y .+= z
         end
     end
     return y

src/linearcombination.jl

@@ -0,0 +1,38 @@
+struct UniformScalingMap{T} <: LinearMap{T} # T will be determined from maps to which this is added
+    λ::T
+    M::Int
+end
+UniformScalingMap(λ::T, M::Int, N::Int) where {T} =
+    (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square"))
+UniformScalingMap(λ::T, sz::Tuple{Int, Int}) where {T} =
+    (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square"))
+
+# properties
+Base.size(A::UniformScalingMap) = (A.M, A.M)
+Base.isreal(A::UniformScalingMap) = isreal(A.λ)
+LinearAlgebra.issymmetric(::UniformScalingMap) = true
+LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A)
+LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ)
+
+# special transposition behavior
+LinearAlgebra.transpose(A::UniformScalingMap) = A
+LinearAlgebra.adjoint(A::UniformScalingMap)   = UniformScalingMap(conj(A.λ), size(A))
+
+# multiplication with vector
+# call of LinearAlgebra.generic_mul! since order of arguments in mul! in stdlib/LinearAlgebra/src/generic.jl
+# TODO: either leave it as is or use mul! (and lower bound on version) once fixed in LinearAlgebra
+A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
+    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, A.λ, x) : throw(DimensionMismatch("A_mul_B!")))
+Base.:(*)(A::UniformScalingMap, x::AbstractVector) = A.λ * x
+
+At_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
+    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, A.λ, x) : throw(DimensionMismatch("At_mul_B!")))
+
+Ac_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) =
+    (length(x) == length(y) == A.M ? LinearAlgebra.generic_mul!(y, conj(A.λ), x) : throw(DimensionMismatch("Ac_mul_B!")))
+
+# combine LinearMap and UniformScaling objects in linear combinations
+Base.:(+)(A1::LinearMap, A2::UniformScaling{T}) where {T} = A1 + UniformScalingMap{T}(A2[1,1], size(A1, 1))
+Base.:(+)(A1::UniformScaling{T}, A2::LinearMap) where {T} = UniformScalingMap{T}(A1[1,1], size(A2, 1)) + A2
+Base.:(-)(A1::LinearMap, A2::UniformScaling{T}) where {T} = A1 - UniformScalingMap{T}(A2[1,1], size(A1, 1))
+Base.:(-)(A1::UniformScaling{T}, A2::LinearMap) where {T} = UniformScalingMap{T}(A1[1,1], size(A2, 1)) - A2

src/uniformscalingmap.jl

@@ -0,0 +1,2 @@
+julia 0.7
+Quaternions