diff --git a/README.md b/README.md index fe509123..b265de9e 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,13 @@ A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently. +## What's new in v2.6 +* New feature: "lazy" Kronecker product, Kronecker sums, and powers thereof + for `LinearMap`s. `AbstractMatrix` objects are promoted to `LinearMap`s if + one of the first 8 Kronecker factors is a `LinearMap` object. +* Compatibility with the generic multiply-and-add interface (a.k.a. 5-arg + `mul!`) introduced in julia v1.3 + ## What's new in v2.5.0 * New feature: concatenation of `LinearMap`s objects with `UniformScaling`s, consistent with (h-, v-, and hc-)concatenation of matrices. Note, matrices diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index bd994b33..a8adf542 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -1,6 +1,7 @@ module LinearMaps export LinearMap +export ⊗, kronsum, ⊕ using LinearAlgebra using SparseArrays @@ -14,6 +15,8 @@ end abstract type LinearMap{T} end +const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}} + Base.eltype(::LinearMap{T}) where {T} = T Base.isreal(A::LinearMap) = eltype(A) <: Real @@ -26,6 +29,26 @@ Base.ndims(::LinearMap) = 2 Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions")) Base.length(A::LinearMap) = size(A)[1] * size(A)[2] +# check dimension consistency for y = A*x and Y = A*X +function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector) + m, n = size(A) + (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!")) + return nothing + end + function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) + m, n = size(A) + (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!")) + return nothing + end + +# conversion of AbstractMatrix to LinearMap +convert_to_lmaps_(A::AbstractMatrix) = 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...)...) + Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x) function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false) length(y) == size(A, 1) || throw(DimensionMismatch("mul!")) @@ -68,50 +91,6 @@ A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x) Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x) -# Matrix: create matrix representation of LinearMap -function Base.Matrix(A::LinearMap) - M, N = size(A) - T = eltype(A) - mat = Matrix{T}(undef, (M, N)) - v = fill(zero(T), N) - @inbounds for i = 1:N - v[i] = one(T) - mul!(view(mat, :, i), A, v) - v[i] = zero(T) - end - return mat -end - -Base.Array(A::LinearMap) = Matrix(A) -Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A) -Base.convert(::Type{Array}, A::LinearMap) = Matrix(A) -Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A) - -# sparse: create sparse matrix representation of LinearMap -function SparseArrays.sparse(A::LinearMap{T}) where {T} - M, N = size(A) - rowind = Int[] - nzval = T[] - colptr = Vector{Int}(undef, N+1) - v = fill(zero(T), N) - Av = Vector{T}(undef, M) - - for i = 1:N - v[i] = one(T) - mul!(Av, A, v) - js = findall(!iszero, Av) - colptr[i] = length(nzval) + 1 - if length(js) > 0 - append!(rowind, js) - append!(nzval, Av[js]) - end - v[i] = zero(T) - end - colptr[N+1] = length(nzval) + 1 - - return SparseMatrixCSC(M, N, colptr, rowind, nzval) -end - include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties 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 @@ -119,13 +98,17 @@ include("linearcombination.jl") # defining linear combinations of linear maps 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("conversion.jl") # conversion of linear maps to matrices """ LinearMap(A; kwargs...) + LinearMap(J, M::Int) LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...) 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`, or +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 on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably, @@ -141,20 +124,21 @@ The keyword arguments and their default values for functions `f` are For existing linear maps or matrices `A`, the default values will be taken by calling `issymmetric`, `ishermitian` and `isposdef` on the existing object `A`. -For functions `f`, there is one more keyword arguments +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 of `f` in the method table. """ -LinearMap(A::Union{AbstractMatrix, LinearMap}; kwargs...) = WrappedMap(A; kwargs...) +LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...) +LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M) LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...) 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::Union{AbstractMatrix, LinearMap}; kwargs...) where {T} = WrappedMap{T}(A; kwargs...) +LinearMap{T}(A::MapOrMatrix; kwargs...) where {T} = WrappedMap{T}(A; kwargs...) LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...) end # module diff --git a/src/composition.jl b/src/composition.jl index 0fefc57c..c8a9ed2e 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -1,9 +1,12 @@ +# helper function +check_dim_mul(A, B) = size(A, 2) == size(B, 1) + struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} maps::As # stored in order of application to vector function CompositeMap{T, As}(maps::As) where {T, As} N = length(maps) for n in 2:N - size(maps[n], 2) == size(maps[n-1], 1) || throw(DimensionMismatch("CompositeMap")) + check_dim_mul(maps[n], maps[n-1]) || throw(DimensionMismatch("CompositeMap")) end for n in 1:N promote_type(T, eltype(maps[n])) == T || throw(InexactError()) @@ -47,7 +50,7 @@ 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)..., UniformScalingMap(α * Alast.λ, size(Alast, 1)))) + return CompositeMap{T}(tuple(Base.front(A.maps)..., α * Alast)) else return CompositeMap{T}(tuple(A.maps..., UniformScalingMap(α, size(A, 1)))) end @@ -60,7 +63,7 @@ function Base.:(*)(A::CompositeMap, α::Number) T = promote_type(typeof(α), eltype(A)) Afirst = first(A.maps) if Afirst isa UniformScalingMap - return CompositeMap{T}(tuple(UniformScalingMap(Afirst.λ * α, size(Afirst, 1)), Base.tail(A.maps)...)) + return CompositeMap{T}(tuple(Afirst * α, Base.tail(A.maps)...)) else return CompositeMap{T}(tuple(UniformScalingMap(α, size(A, 2)), A.maps...)) end @@ -87,22 +90,22 @@ julia> LinearMap(ones(Int, 3, 3)) * CS * I * rand(3, 3); ``` """ function Base.:(*)(A₁::LinearMap, A₂::LinearMap) - size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*")) + check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*")) T = promote_type(eltype(A₁), eltype(A₂)) return CompositeMap{T}(tuple(A₂, A₁)) end function Base.:(*)(A₁::LinearMap, A₂::CompositeMap) - size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*")) + check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*")) T = promote_type(eltype(A₁), eltype(A₂)) return CompositeMap{T}(tuple(A₂.maps..., A₁)) end function Base.:(*)(A₁::CompositeMap, A₂::LinearMap) - size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*")) + check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*")) T = promote_type(eltype(A₁), eltype(A₂)) return CompositeMap{T}(tuple(A₂, A₁.maps...)) end function Base.:(*)(A₁::CompositeMap, A₂::CompositeMap) - size(A₁, 2) == size(A₂, 1) || throw(DimensionMismatch("*")) + check_dim_mul(A₁, A₂) || throw(DimensionMismatch("*")) T = promote_type(eltype(A₁), eltype(A₂)) return CompositeMap{T}(tuple(A₂.maps..., A₁.maps...)) end diff --git a/src/conversion.jl b/src/conversion.jl new file mode 100644 index 00000000..f3e9a45e --- /dev/null +++ b/src/conversion.jl @@ -0,0 +1,76 @@ +# Matrix: create matrix representation of LinearMap +function Base.Matrix(A::LinearMap) + M, N = size(A) + T = eltype(A) + mat = Matrix{T}(undef, (M, N)) + v = fill(zero(T), N) + @inbounds for i = 1:N + v[i] = one(T) + mul!(view(mat, :, i), A, v) + v[i] = zero(T) + end + return mat +end +Base.Array(A::LinearMap) = Matrix(A) +Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A) +Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A) +Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A) +Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A) + +# special cases +Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap) +Base.convert(::Type{Matrix}, A::WrappedMap) = convert(Matrix, A.lmap) +function Base.convert(::Type{Matrix}, ΣA::LinearCombination{<:Any,<:Tuple{Vararg{MatrixMap}}}) + if length(ΣA.maps) <= 10 + return (+).(map(A->getfield(A, :lmap), ΣA.maps)...) + else + S = zero(first(ΣA.maps).lmap) + for A in ΣA.maps + S .+= A.lmap + end + return S + end +end +function Base.convert(::Type{Matrix}, AB::CompositeMap{<:Any,<:Tuple{MatrixMap,MatrixMap}}) + B, A = AB.maps + return A.lmap*B.lmap +end +function Base.convert(::Type{Matrix}, λA::CompositeMap{<:Any,<:Tuple{MatrixMap,UniformScalingMap}}) + A, J = λA.maps + return J.λ*A.lmap +end +function Base.convert(::Type{Matrix}, Aλ::CompositeMap{<:Any,<:Tuple{UniformScalingMap,MatrixMap}}) + J, A = Aλ.maps + return A.lmap*J.λ +end + +Base.Matrix(A::BlockMap) = hvcat(A.rows, convert.(AbstractMatrix, A.maps)...) + +# sparse: create sparse matrix representation of LinearMap +function SparseArrays.sparse(A::LinearMap{T}) where {T} + M, N = size(A) + rowind = Int[] + nzval = T[] + colptr = Vector{Int}(undef, N+1) + v = fill(zero(T), N) + Av = Vector{T}(undef, M) + + for i = 1:N + v[i] = one(T) + mul!(Av, A, v) + js = findall(!iszero, Av) + colptr[i] = length(nzval) + 1 + if length(js) > 0 + append!(rowind, js) + append!(nzval, Av[js]) + end + v[i] = zero(T) + end + colptr[N+1] = length(nzval) + 1 + + return SparseMatrixCSC(M, N, colptr, rowind, nzval) +end + +Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A) +Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap) +Base.convert(::Type{SparseMatrixCSC}, A::BlockMap) = hvcat(A.rows, convert.(SparseMatrixCSC, A.maps)...) diff --git a/src/kronecker.jl b/src/kronecker.jl new file mode 100644 index 00000000..17f9511e --- /dev/null +++ b/src/kronecker.jl @@ -0,0 +1,281 @@ +struct KroneckerMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} + maps::As + function KroneckerMap{T, As}(maps::As) where {T, As} + for A in maps + promote_type(T, eltype(A)) == T || throw(InexactError()) + end + return new{T,As}(maps) + end +end + +KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap{T, As}(maps) + +""" + kron(A::LinearMap, B::LinearMap) + kron(A, B, Cs...) + +Construct a `KroneckerMap <: LinearMap` object, a (lazy) representation of the +Kronecker product of two `LinearMap`s. One of the two factors can be an `AbstractMatrix`, +which is then promoted to a `LinearMap` automatically. + +To avoid fallback to the generic [`Base.kron`](@ref) in the multi-map case, +there must be a `LinearMap` object among the first 8 arguments in usage like +`kron(A, B, Cs...)`. + +For convenience, one can also use `A ⊗ B` or `⊗(A, B, Cs...)` (typed as `\\otimes+TAB`) to construct the +`KroneckerMap`, even when all arguments are of `AbstractMatrix` type. + +If `A`, `B`, `C` and `D` are linear maps of such size that one can form the matrix +products `A*C` and `B*D`, then the mixed-product property `(A⊗B)*(C⊗D) = (A*C)⊗(B*D)` +holds. Upon vector multiplication, this rule is checked for applicability. + +# Examples +```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps) +julia> J = LinearMap(I, 2) # 2×2 identity map +LinearMaps.UniformScalingMap{Bool}(true, 2) + +julia> E = spdiagm(-1 => trues(1)); D = E + E' - 2I; + +julia> Δ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator + +julia> Matrix(Δ) +4×4 Array{Int64,2}: + -4 1 1 0 + 1 -4 0 1 + 1 0 -4 1 + 0 1 1 -4 +``` +""" +Base.kron(A::LinearMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B)) +Base.kron(A::LinearMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A, B.maps...)) +Base.kron(A::KroneckerMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B)) +Base.kron(A::KroneckerMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B.maps...)) +Base.kron(A::LinearMap, B::LinearMap, Cs::LinearMap...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(tuple(A, B, Cs...)) +Base.kron(A::AbstractMatrix, B::LinearMap) = kron(LinearMap(A), B) +Base.kron(A::LinearMap, B::AbstractMatrix) = kron(A, LinearMap(B)) +# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product +for k in 3:8 # is 8 sufficient? + Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1)) + # yields (:A1, :A2, :A3, ..., :A(k-1)) + L = :($(Symbol(:A,k))::LinearMap) + # yields :Ak::LinearMap + mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) + # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) + + @eval Base.kron($(Is...), $L, As::MapOrMatrix...) = + kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) +end + +struct KronPower{p} + function KronPower(p::Integer) + p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return new{p}() + end +end + +""" + ⊗(k::Integer) + +Construct a lazy representation of the `k`-th Kronecker power +`A^⊗(k) = A ⊗ A ⊗ ... ⊗ A`, where `A` can be an `AbstractMatrix` or a `LinearMap`. +""" +⊗(k::Integer) = KronPower(k) + +⊗(A, B, Cs...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(convert_to_lmaps(A, B, Cs...)) + +Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {p} = + ⊗(ntuple(n -> convert_to_lmaps_(A), Val(p))...) + +Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) + +LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = issymmetric(A) +LinearAlgebra.ishermitian(A::KroneckerMap) = all(ishermitian, A.maps) + +LinearAlgebra.adjoint(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(transpose, A.maps)) + +Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) + +################# +# multiplication helper functions +################# + +@inline function _kronmul!(y, B, X, At, T) + na, ma = size(At) + mb, nb = size(B) + v = zeros(T, ma) + Ty = eltype(y) + temp1 = Array{Ty}(undef, na) + temp2 = Array{Ty}(undef, nb) + @views @inbounds for i in 1:ma + v[i] = one(T) + A_mul_B!(temp1, At, v) + A_mul_B!(temp2, X, temp1) + A_mul_B!(y[((i-1)*mb+1):i*mb], B, temp2) + v[i] = zero(T) + end + return y +end +@inline function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T) + na, ma = size(At) + mb, nb = size(B) + if nb*ma < mb*na + mul!(reshape(y, (mb, ma)), B, convert(Matrix, X*At)) + else + mul!(reshape(y, (mb, ma)), convert(Matrix, B*X), At isa MatrixMap ? At.lmap : At.λ) + end + return y +end + +################# +# multiplication with vectors +################# + +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} + require_one_based_indexing(y) + @boundscheck check_dim_mul(y, L, x) + A, B = L.maps + X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false) + _kronmul!(y, B, X, transpose(A), T) + return y +end +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} + require_one_based_indexing(y) + @boundscheck check_dim_mul(y, L, x) + A = first(L.maps) + B = kron(Base.tail(L.maps)...) + X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false) + _kronmul!(y, B, X, transpose(A), T) + return y +end +# mixed-product rule, prefer the right if possible +# (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ) +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) + B, A = L.maps + if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps)) + A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x) + else + A_mul_B!(y, LinearMap(A)*B, x) + end +end +# mixed-product rule, prefer the right if possible +# (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} + As = map(AB -> AB.maps[1], L.maps) + Bs = map(AB -> AB.maps[2], L.maps) + As1, As2 = Base.front(As), Base.tail(As) + Bs1, Bs2 = Base.front(Bs), Base.tail(Bs) + apply = all(A -> check_dim_mul(A...), zip(As1, As2)) && all(A -> check_dim_mul(A...), zip(Bs1, Bs2)) + if apply + A_mul_B!(y, kron(prod(As), prod(Bs)), x) + else + A_mul_B!(y, CompositeMap{T}(map(LinearMap, L.maps)), x) + end +end + +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = + A_mul_B!(y, transpose(A), x) + +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = + A_mul_B!(y, adjoint(A), x) + +############### +# KroneckerSumMap +############### +struct KroneckerSumMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T} + maps::As + function KroneckerSumMap{T, As}(maps::As) where {T, As} + for A in maps + size(A, 1) == size(A, 2) || throw(ArgumentError("operators need to be square in Kronecker sums")) + promote_type(T, eltype(A)) == T || throw(InexactError()) + end + return new{T,As}(maps) + end +end + +KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerSumMap{T, As}(maps) + +""" + kronsum(A, B) + kronsum(A, B, Cs...) + +Construct a `KroneckerSumMap <: LinearMap` object, a (lazy) representation of the +Kronecker sum `A⊕B = A ⊗ Ib + Ia ⊗ B` of two square linear maps of type +`LinearMap` or `AbstractMatrix`. Here, `Ia` and `Ib` are identity operators of +the size of `A` and `B`, respectively. Arguments of type `AbstractMatrix` are +automatically promoted to `LinearMap`s. + +For convenience, one can also use `A ⊕ B` or `⊕(A, B, Cs...)` (typed as +`\\oplus+TAB`) to construct the `KroneckerSumMap`. + +# Examples +```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps) +julia> J = LinearMap(I, 2) # 2×2 identity map +LinearMaps.UniformScalingMap{Bool}(true, 2) + +julia> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I); + +julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum + +julia> Δ₂ = kronsum(D, D); + +julia> Δ₃ = D^⊕(2); + +julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃) +true +``` +""" +kronsum(A::MapOrMatrix, B::MapOrMatrix) = + KroneckerSumMap{promote_type(eltype(A), eltype(B))}(convert_to_lmaps(A, B)) +kronsum(A::MapOrMatrix, B::MapOrMatrix, C::MapOrMatrix, Ds::MapOrMatrix...) = + kronsum(A, kronsum(B, C, Ds...)) + +struct KronSumPower{p} + function KronSumPower(p::Integer) + p > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k")) + return new{p}() + end +end + +""" + ⊕(k::Integer) + +Construct a lazy representation of the `k`-th Kronecker sum power `A^⊕(k) = A ⊕ A ⊕ ... ⊕ A`, +where `A` can be a square `AbstractMatrix` or a `LinearMap`. +""" +⊕(k::Integer) = KronSumPower(k) + +⊕(a, b, c...) = kronsum(a, b, c...) + +Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} = kronsum(ntuple(n -> convert_to_lmaps_(A), Val(p))...) + +Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i)) +Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2)) + +LinearAlgebra.issymmetric(A::KroneckerSumMap) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::KroneckerSumMap{<:Real}) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::KroneckerSumMap) = all(ishermitian, A.maps) + +LinearAlgebra.adjoint(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(transpose, A.maps)) + +Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) + +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) + @boundscheck check_dim_mul(y, L, x) + A, B = L.maps + ma, na = size(A) + mb, nb = size(B) + X = reshape(x, (nb, na)) + Y = reshape(y, (nb, na)) + mul!(Y, X, convert(AbstractMatrix, transpose(A))) + mul!(Y, B, X, true, true) + return y +end + +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = + A_mul_B!(y, transpose(A), x) + +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = + A_mul_B!(y, adjoint(A), x) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 7c2de56a..25461858 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -8,7 +8,7 @@ struct UniformScalingMap{T} <: LinearMap{T} end UniformScalingMap(λ::Number, M::Int, N::Int) = (M == N ? UniformScalingMap(λ, M) : error("UniformScalingMap needs to be square")) -UniformScalingMap(λ::Number, sz::Tuple{Int, Int}) = +UniformScalingMap(λ::T, sz::Dims{2}) where {T} = (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) # properties @@ -18,10 +18,17 @@ LinearAlgebra.issymmetric(::UniformScalingMap) = true LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A) LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ) +# comparison of UniformScalingMap objects +Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (A.λ == B.λ && A.M == B.M) + # special transposition behavior LinearAlgebra.transpose(A::UniformScalingMap) = A LinearAlgebra.adjoint(A::UniformScalingMap) = UniformScalingMap(conj(A.λ), size(A)) +# multiplication with scalar +Base.:(*)(α::Number, J::UniformScalingMap) = UniformScalingMap(α * J.λ, size(J)) +Base.:(*)(J::UniformScalingMap, α::Number) = UniformScalingMap(J.λ * α, size(J)) + # multiplication with vector Base.:(*)(A::UniformScalingMap, x::AbstractVector) = length(x) == A.M ? A.λ * x : throw(DimensionMismatch("A_mul_B!")) @@ -30,7 +37,7 @@ if VERSION < v"1.3.0-alpha.115" function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) (length(x) == length(y) == A.M || throw(DimensionMismatch("A_mul_B!"))) if iszero(A.λ) - return fill!(y, 0) + return fill!(y, zero(eltype(y))) elseif isone(A.λ) return copyto!(y, x) else @@ -40,25 +47,31 @@ function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) end end else # 5-arg mul! exists and order of arguments is corrected -function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) - (length(x) == length(y) == A.M || throw(DimensionMismatch("A_mul_B!"))) - λ = A.λ - if iszero(λ) - return fill!(y, 0) - elseif isone(λ) - return copyto!(y, x) - else - return y .= λ .* x - end -end + +A_mul_B!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector) = mul!(y, J, x) + end # VERSION -function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false) +@inline function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false) @boundscheck (length(x) == length(y) == J.M || throw(DimensionMismatch("mul!"))) + _scaling!(y, J, x, α, β) + return y +end + +@inline function LinearAlgebra.mul!(Y::AbstractMatrix, J::UniformScalingMap, X::AbstractMatrix, α::Number=true, β::Number=false) + @boundscheck size(X) == size(Y) || throw(DimensionMismatch("mul!")) + @boundscheck size(X,1) == J.M || throw(DimensionMismatch("mul!")) + _scaling!(Y, J, X, α, β) + return Y +end + +function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false) λ = J.λ - @inbounds if isone(α) + if isone(α) if iszero(β) - A_mul_B!(y, J, x) + iszero(λ) && return fill!(y, zero(eltype(y))) + isone(λ) && return copyto!(y, x) + y .= λ .* x return y elseif isone(β) iszero(λ) && return y diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 5f919be5..8ebabc9e 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -1,16 +1,16 @@ -struct WrappedMap{T, A<:Union{AbstractMatrix, LinearMap}} <: LinearMap{T} +struct WrappedMap{T, A<:MapOrMatrix} <: LinearMap{T} lmap::A _issymmetric::Bool _ishermitian::Bool _isposdef::Bool end -function WrappedMap(lmap::Union{AbstractMatrix{T}, LinearMap{T}}; +function WrappedMap(lmap::MapOrMatrix{T}; issymmetric::Bool = issymmetric(lmap), ishermitian::Bool = ishermitian(lmap), isposdef::Bool = isposdef(lmap)) where {T} WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef) end -function WrappedMap{T}(lmap::Union{AbstractMatrix, LinearMap}; +function WrappedMap{T}(lmap::MapOrMatrix; issymmetric::Bool = issymmetric(lmap), ishermitian::Bool = ishermitian(lmap), isposdef::Bool = isposdef(lmap)) where {T} @@ -33,6 +33,14 @@ if VERSION ≥ v"1.3.0-alpha.115" LinearAlgebra.mul!(y::AbstractVector, A::WrappedMap, x::AbstractVector, α::Number=true, β::Number=false) = mul!(y, A.lmap, x, α, β) +LinearAlgebra.mul!(Y::AbstractMatrix, A::MatrixMap, X::AbstractMatrix, α::Number=true, β::Number=false) = + mul!(Y, A.lmap, X, α, β) + +else + +LinearAlgebra.mul!(Y::AbstractMatrix, A::MatrixMap, X::AbstractMatrix) = + mul!(Y, A.lmap, X) + end # VERSION # properties diff --git a/test/conversion.jl b/test/conversion.jl new file mode 100644 index 00000000..f4fb92aa --- /dev/null +++ b/test/conversion.jl @@ -0,0 +1,41 @@ +using Test, LinearMaps, LinearAlgebra, SparseArrays + +@testset "conversion" begin + A = 2 * rand(ComplexF64, (20, 10)) .- 1 + v = rand(ComplexF64, 10) + w = rand(ComplexF64, 20) + V = rand(ComplexF64, 10, 3) + W = rand(ComplexF64, 20, 3) + α = rand() + β = rand() + M = @inferred LinearMap(A) + N = @inferred LinearMap(M) + + @test Matrix(M) == A + @test Array(M) == A + @test convert(AbstractArray, M) == A + @test convert(AbstractMatrix, M) === A + @test convert(Matrix, M) === A + @test convert(Array, M) === A + @test convert(AbstractMatrix, M) == A + @test convert(Matrix, M) == A + @test convert(Array, M) == A + @test Matrix(M') == A' + @test Matrix(transpose(M)) == copy(transpose(A)) + @test convert(AbstractMatrix, M') isa Adjoint + @test convert(Matrix, M*3I) == A*3 + + # sparse matrix generation/conversion + @test sparse(M) == sparse(Array(M)) + @test convert(SparseMatrixCSC, M) == sparse(Array(M)) + + B = copy(A) + B[rand(1:length(A), 30)] .= 0 + MS = LinearMap(B) + @test sparse(MS) == sparse(Array(MS)) == sparse(B) + + J = LinearMap(I, 10) + @test J isa LinearMap{Bool} + @test sparse(J) == Matrix{Bool}(I, 10, 10) + @test nnz(sparse(J)) == 10 +end diff --git a/test/kronecker.jl b/test/kronecker.jl new file mode 100644 index 00000000..564c5350 --- /dev/null +++ b/test/kronecker.jl @@ -0,0 +1,68 @@ +using Test, LinearMaps, LinearAlgebra + +@testset "kronecker products" begin + @testset "Kronecker product" begin + A = rand(ComplexF64, 3, 3) + B = rand(ComplexF64, 2, 2) + K = kron(A, B) + LA = LinearMap(A) + LB = LinearMap(B) + LK = @inferred kron(LA, LB) + @test @inferred size(LK) == size(K) + for i in (1, 2) + @test @inferred size(LK, i) == size(K, i) + end + @test LK isa LinearMaps.KroneckerMap{ComplexF64} + for transform in (identity, transpose, adjoint) + @test Matrix(transform(LK)) ≈ transform(Matrix(LK)) ≈ transform(kron(A, B)) + @test Matrix(kron(transform(LA), transform(LB))) ≈ transform(kron(A, B)) + @test Matrix(transform(LinearMap(LK))) ≈ transform(Matrix(LK)) ≈ transform(kron(A, B)) + end + @test kron(LA, kron(LA, B)) == kron(LA, LA, LB) + @test kron(kron(LA, LB), kron(LA, LB)) == kron(LA, LB, LA, LB) + @test kron(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) ≈ Matrix(@inferred LA^⊗(3)) ≈ Matrix(@inferred A^⊗(3)) + K = @inferred kron(A, A, A, LA) + @test K isa LinearMaps.KroneckerMap + @test Matrix(K) ≈ kron(A, A, A, A) + @test Matrix(@inferred K*K) ≈ kron(A, A, A, A)*kron(A, A, A, A) + K4 = @inferred kron(A, B, B, LB) + # check that matrices don't get Kronecker-multiplied, but that all is lazy + @test K4.maps[1].lmap === A + @test @inferred kron(LA, LB)' == @inferred kron(LA', LB') + @test (@inferred kron(LA, B)) == (@inferred kron(LA, LB)) == (@inferred kron(A, LB)) + @test @inferred ishermitian(kron(LA'LA, LB'LB)) + A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B) + @test @inferred issymmetric(kron(LA'LA, LB'LB)) + @test @inferred ishermitian(kron(LA'LA, LB'LB)) + # use mixed-product rule + K = kron(LA, LB) * kron(LA, LB) * kron(LA, LB) + @test Matrix(K) ≈ kron(A, B)^3 + # example that doesn't use mixed-product rule + A = rand(3, 2); B = rand(2, 3) + K = @inferred kron(A, LinearMap(B)) + @test Matrix(@inferred K*K) ≈ kron(A, B)*kron(A, B) + A = rand(3, 2); B = rand(4, 3) + @test Matrix(kron(LinearMap(A), B, [A A])*kron(LinearMap(A), B, A')) ≈ kron(A, B, [A A])*kron(A, B, A') + end + @testset "Kronecker sum" begin + for elty in (Float64, ComplexF64) + A = rand(elty, 3, 3) + B = rand(elty, 2, 2) + LA = LinearMap(A) + LB = LinearMap(B) + KS = @inferred kronsum(LA, B) + KSmat = kron(A, Matrix(I, 2, 2)) + kron(Matrix(I, 3, 3), B) + @test Matrix(KS) ≈ Matrix(kron(A, LinearMap(I, 2)) + kron(LinearMap(I, 3), B)) + @test size(KS) == size(kron(A, Matrix(I, 2, 2))) + for transform in (identity, transpose, adjoint) + @test Matrix(transform(KS)) ≈ transform(Matrix(KS)) ≈ transform(KSmat) + @test Matrix(kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) + @test Matrix(transform(LinearMap(kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) + end + @inferred kronsum(A, A, LB) + @test Matrix(@inferred LA^⊕(3)) == Matrix(@inferred A^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) + @test @inferred(kronsum(LA, LA, LB)) == @inferred(kronsum(LA, kronsum(LA, LB))) == @inferred(kronsum(A, A, B)) + @test Matrix(@inferred kronsum(A, B, A, B, A, B)) ≈ Matrix(@inferred kronsum(LA, LB, LA, LB, LA, LB)) + end + end +end diff --git a/test/linearmaps.jl b/test/linearmaps.jl index b1ffe25e..7570da02 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -19,21 +19,6 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools @test ndims(M) == 2 @test_throws ErrorException size(M, 3) @test length(M) == length(A) - # matrix generation/conversion - @test Matrix(M) == A - @test Array(M) == A - @test convert(Matrix, M) == A - @test convert(Array, M) == A - @test Matrix(M') == A' - @test Matrix(transpose(M)) == copy(transpose(A)) - # sparse matrix generation/conversion - @test sparse(M) == sparse(Array(M)) - @test convert(SparseMatrixCSC, M) == sparse(Array(M)) - - B = copy(A) - B[rand(1:length(A), 30)] .= 0 - MS = LinearMap(B) - @test sparse(MS) == sparse(Array(MS)) == sparse(B) end Av = A * v diff --git a/test/runtests.jl b/test/runtests.jl index 8ed8ba4b..6851f7b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,3 +17,7 @@ include("uniformscalingmap.jl") include("numbertypes.jl") include("blockmap.jl") + +include("kronecker.jl") + +include("conversion.jl") diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index 140592c9..ab38ce88 100644 --- a/test/uniformscalingmap.jl +++ b/test/uniformscalingmap.jl @@ -9,7 +9,7 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools LC = @inferred M + N v = rand(ComplexF64, 10) w = similar(v) - Id = @inferred LinearMaps.UniformScalingMap(1, 10) + Id = @inferred LinearMap(I, 10) @test_throws ErrorException LinearMaps.UniformScalingMap(1, 10, 20) @test_throws ErrorException LinearMaps.UniformScalingMap(1, (10, 20)) @test size(Id) == (10, 10) @@ -22,10 +22,11 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools @test (3 * I + 2 * M') * v == 2 * A'v + 3v @test (2 * M' - 3 * I) * v == 2 * A'v - 3v @test (3 * I - 2 * M') * v == -2 * A'v + 3v + @test (3 * I - 2 * M') * v == -2 * A'v + 3v @test transpose(LinearMap(2 * M' + 3 * I)) * v ≈ transpose(2 * A' + 3 * I) * v - @test LinearMap(2 * M' + 3 * I)' * v ≈ (2 * A' + 3 * I)' * v + @test LinearMap(2 * M' + 0I)' * v ≈ (2 * A')' * v for λ in (0, 1, rand()), α in (0, 1, rand()), β in (0, 1, rand()) - Λ = @inferred LinearMaps.UniformScalingMap(λ, 10) + Λ = @inferred LinearMap(λ*I, 10) x = rand(10) y = rand(10) b = @benchmarkable mul!($y, $Λ, $x, $α, $β) @@ -40,4 +41,6 @@ using Test, LinearMaps, LinearAlgebra, BenchmarkTools J = @inferred LinearMap(LinearMaps.UniformScalingMap(λ, 10)) @test transform(J) * x == transform(λ) * x end + X = rand(10, 10); Y = similar(X) + @test mul!(Y, Id, X) == X end