From 77a70a653a2b534b84b8c0935a5c6108f73e6748 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 25 Aug 2019 19:25:34 +0200 Subject: [PATCH 01/33] add Kronecker products for LinearMaps --- src/LinearMaps.jl | 1 + src/kronecker.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+) create mode 100644 src/kronecker.jl diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index bd994b33..f6236e43 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -119,6 +119,7 @@ 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 """ LinearMap(A; kwargs...) diff --git a/src/kronecker.jl b/src/kronecker.jl new file mode 100644 index 00000000..5ffa1991 --- /dev/null +++ b/src/kronecker.jl @@ -0,0 +1,46 @@ +struct KroneckerMap{T, As<:Tuple{LinearMap,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{LinearMap,LinearMap}} = KroneckerMap{T, As}(maps) +Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) +Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) +Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), B) + +Base.size(A::KroneckerMap) = size(A.maps[1]) .* size(A.maps[2]) + +LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = all(issymmetric, A.maps) +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) + +function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap, x::AbstractVector) + A, B = L.maps + ma, na = size(A) + mb, nb = size(B) + (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("kronecker product")) + X = LinearMap(reshape(x, (nb, na))) + M = B * X * transpose(A) + T = eltype(L) + v = zeros(T, ma) + @views @inbounds for i in 1:ma + v[i] = one(T) + mul!(y[((i-1)*mb+1):i*mb], M, v) + v[i] = zero(T) + end + return y +end + +LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) + +LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) From 94b000665ef6821a5d7a71098ee17be26bad6b27 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 28 Aug 2019 15:44:13 +0200 Subject: [PATCH 02/33] announce feature in README --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index fe509123..6b56b92b 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,11 @@ 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 of `LinearMap`s. `AbstractMatrix`es are + promoted to `LinearMap`s if one of the first 8 Kronecker factors is a + `LinearMap` object. + ## 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 From 33b6a0a25db637b47aa91870548aee9f8c67741d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 28 Aug 2019 15:47:48 +0200 Subject: [PATCH 03/33] add promotion of matrices, tests --- src/kronecker.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/kronecker.jl b/src/kronecker.jl index 5ffa1991..da4a3a7c 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -12,6 +12,13 @@ KroneckerMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerM Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), B) +for k in 3:8 # is 8 sufficient? + Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1)) + L = :($(Symbol(:A,k))::LinearMap) + mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) + + @eval Base.kron($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = kron($(mapargs...), $(Symbol(:A,k)), As...) +end Base.size(A::KroneckerMap) = size(A.maps[1]) .* size(A.maps[2]) From e2ea192c23410776c937b27d622de9cdfc462610 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 28 Aug 2019 17:04:40 +0200 Subject: [PATCH 04/33] add comments, use mixed-product trick --- src/kronecker.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index da4a3a7c..8c312a52 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -12,12 +12,17 @@ KroneckerMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerM Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), 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 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::Union{LinearMap,AbstractMatrix}...) = kron($(mapargs...), $(Symbol(:A,k)), As...) + @eval Base.kron($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = + kron($(mapargs...), $(Symbol(:A,k)), As...) end Base.size(A::KroneckerMap) = size(A.maps[1]) .* size(A.maps[2]) @@ -29,6 +34,14 @@ 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)) +function Base.:(*)(A::KroneckerMap, B::KroneckerMap) + if (size(A.maps[1], 2) == size(B.maps[1], 1) && size(A.maps[2], 2) == size(B.maps[2], 1)) + return kron(A.maps[1]*B.maps[1], A.maps[2]*B.maps[2]) + else + return CompositeMap{T}(tuple(B, A)) + end +end + Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap, x::AbstractVector) From d7e027f9b53c432a3afddd410f4833c8298c8eeb Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 2 Sep 2019 13:25:42 +0200 Subject: [PATCH 05/33] add size(A, i) for kronecker --- src/kronecker.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 8c312a52..a08025cf 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -25,7 +25,8 @@ for k in 3:8 # is 8 sufficient? kron($(mapargs...), $(Symbol(:A,k)), As...) end -Base.size(A::KroneckerMap) = size(A.maps[1]) .* size(A.maps[2]) +Base.size(A::KroneckerMap, i) = prod(size.(A.maps, i)) +Base.size(A::KroneckerMap) = (size(A, 1), size(A, 2)) LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = all(issymmetric, A.maps) From cca7f879805cc3faf2d1729febc0a675b64cbe32 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 30 Aug 2019 12:38:09 +0200 Subject: [PATCH 06/33] export UniformScalingMaps, add scalar mul, tests --- src/LinearMaps.jl | 2 +- src/uniformscalingmap.jl | 29 ++++++++++++++++++++++++++++- 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index f6236e43..db3257a5 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -1,6 +1,6 @@ module LinearMaps -export LinearMap +export LinearMap, UniformScalingMap using LinearAlgebra using SparseArrays diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 7c2de56a..520314d2 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -1,3 +1,22 @@ +""" + UniformScalingMap(λ::Number, n::Int) + +Construct an `n×n` linear map that acts on vectors by uniformly scaling with `λ`. +This should only be used in the construction of linear maps in which the size of +the uniform scaling cannot be inferred automatically, like in the Kronecker product. +In all other cases, it is recommended to use `LinearAlgebra`'s `I::UniformScaling`. + +## Examples +```jldoctest; setup = :(using LinearMaps) +julia> J = UniformScalingMap(2.0, 10); + +julia> size(J) +(10, 10) + +julia> eltype(J) +Float64 +``` +""" struct UniformScalingMap{T} <: LinearMap{T} λ::T M::Int @@ -8,20 +27,28 @@ 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 +Base.size(A::UniformScalingMap, n) = (n==1 || n==2 ? A.M : error("LinearMap objects have only 2 dimensions")) 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.λ) +# comparison of UniformScalingMap objects, sufficient but not necessary +Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (eltype(A) == eltype(B) && 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!")) From 2e6795549e571620b6fd58e831a27d2c9fc07928 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 30 Aug 2019 12:43:56 +0200 Subject: [PATCH 07/33] simplify scalar multiplication of composites --- src/composition.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index 0fefc57c..c32553dc 100644 --- a/src/composition.jl +++ b/src/composition.jl @@ -47,7 +47,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 +60,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 From 8fd6ac0b2bed16b263c248ffc63f9b6afb9b9206 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 2 Sep 2019 13:45:06 +0200 Subject: [PATCH 08/33] adopt to new test suite --- test/kronecker.jl | 32 ++++++++++++++++++++++++++++++++ test/runtests.jl | 2 ++ 2 files changed, 34 insertions(+) create mode 100644 test/kronecker.jl diff --git a/test/kronecker.jl b/test/kronecker.jl new file mode 100644 index 00000000..1b824e66 --- /dev/null +++ b/test/kronecker.jl @@ -0,0 +1,32 @@ +using Test, LinearMaps, LinearAlgebra + +@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(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) + K = @inferred kron(A, A, A, LA) + @test K isa LinearMaps.KroneckerMap + @test Matrix(K) ≈ kron(A, A, A, A) + K4 = kron(A, B, B, LB) + @test K4.maps[1].maps[1].maps[1].lmap === A + @test @inferred kron(LA, LB)' == @inferred kron(LA', LB') + @test kron(LA, B) == kron(LA, LB) == kron(A, LB) + @test ishermitian(kron(LA'LA, LB'LB)) + A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B) + @test issymmetric(kron(LA'LA, LB'LB)) + @test ishermitian(kron(LA'LA, LB'LB)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 8ed8ba4b..64fc208d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,3 +17,5 @@ include("uniformscalingmap.jl") include("numbertypes.jl") include("blockmap.jl") + +include("kronecker.jl") From 5a4f073cbc13144787715df47dd8ce5d0f9deb9a Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 2 Sep 2019 13:45:33 +0200 Subject: [PATCH 09/33] prepare for require_one_based_indexing --- src/kronecker.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/kronecker.jl b/src/kronecker.jl index a08025cf..56be2a08 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -46,6 +46,7 @@ end Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap, x::AbstractVector) + # require_one_based_indexing(y, x) A, B = L.maps ma, na = size(A) mb, nb = size(B) From 2301a318c8c938cd5735740463b55fa1065d80f2 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 3 Sep 2019 17:48:33 +0200 Subject: [PATCH 10/33] revert export of UnifformScalingMap, add LinearMap constructor instead --- src/LinearMaps.jl | 3 ++- src/uniformscalingmap.jl | 19 ------------------- test/linearmaps.jl | 5 +++++ 3 files changed, 7 insertions(+), 20 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index db3257a5..d635e117 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -1,6 +1,6 @@ module LinearMaps -export LinearMap, UniformScalingMap +export LinearMap using LinearAlgebra using SparseArrays @@ -150,6 +150,7 @@ For functions `f`, there is one more keyword arguments of `f` in the method table. """ LinearMap(A::Union{AbstractMatrix, LinearMap}; 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...) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 520314d2..2f2f981e 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -1,22 +1,3 @@ -""" - UniformScalingMap(λ::Number, n::Int) - -Construct an `n×n` linear map that acts on vectors by uniformly scaling with `λ`. -This should only be used in the construction of linear maps in which the size of -the uniform scaling cannot be inferred automatically, like in the Kronecker product. -In all other cases, it is recommended to use `LinearAlgebra`'s `I::UniformScaling`. - -## Examples -```jldoctest; setup = :(using LinearMaps) -julia> J = UniformScalingMap(2.0, 10); - -julia> size(J) -(10, 10) - -julia> eltype(J) -Float64 -``` -""" struct UniformScalingMap{T} <: LinearMap{T} λ::T M::Int diff --git a/test/linearmaps.jl b/test/linearmaps.jl index b1ffe25e..47f8cc48 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -34,6 +34,11 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools 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 Av = A * v From d9e3fc9b3c398f04816396424781491c78bd96a3 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 3 Sep 2019 17:56:20 +0200 Subject: [PATCH 11/33] update docstring --- src/LinearMaps.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index d635e117..ad145dfd 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -123,10 +123,12 @@ include("kronecker.jl") # Kronecker product of linear maps """ 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, @@ -142,7 +144,7 @@ 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`). From 5c71ae929001c36e2837c41853fcad61b8ba4830 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 6 Sep 2019 09:51:09 +0200 Subject: [PATCH 12/33] add docstring to kron --- src/kronecker.jl | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/kronecker.jl b/src/kronecker.jl index 56be2a08..7c1cfeee 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -9,6 +9,38 @@ struct KroneckerMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T} end KroneckerMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerMap{T, As}(maps) + +""" + kron(A::LinearMap, B::LinearMap) + +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` +object, which then is promoted to `LinearMap` automatically. To avoid fallback to +the generic [`Base.kron`](@ref), there must be a `LinearMap` object among the +first 8 arguments in usage like `kron(A, B, Cs...)`. + +Note: 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)` +holds. Upon multiplication of Kronecker products, this rule is checked for +applicability, which leads to type-instability with a union of two types. + +# 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{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), B) From e00632b9474db8b031a7b6e5200e600cbf3934dc Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 7 Sep 2019 14:29:23 +0200 Subject: [PATCH 13/33] add KroneckerSumMap, yet untested --- src/kronecker.jl | 92 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 1 deletion(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 7c1cfeee..54caff67 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -78,7 +78,7 @@ end Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap, x::AbstractVector) - # require_one_based_indexing(y, x) + # require_one_based_indexing(y) A, B = L.maps ma, na = size(A) mb, nb = size(B) @@ -98,3 +98,93 @@ end LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) + +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::LinearMap, B::LinearMap) + +Construct a `KroneckerSumMap <: LinearMap` object, a (lazy) representation of the +Kronecker sum `A⊕B = kron(A, Ib) + kron(Ia, B)` of two square `LinearMap`s. Here, +`Ia` and `Ib` are identity One of the two factors can be an `AbstractMatrix` +object, which then is promoted to `LinearMap` automatically. To avoid fallback to +the generic [`Base.kron`](@ref), there must be a `LinearMap` object among the +first 8 arguments in usage like `kron(A, B, Cs...)`. + +Note: 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)` +holds. Upon multiplication of Kronecker products, this rule is checked for +applicability, which leads to type-instability with a union of two types. + +# 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.kronsum(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) +Base.kronsum(A::LinearMap, Bs::LinearMap...) = kronsum(A, kronsum(Bs...)) +Base.kronsum(A::LinearMap, B::AbstractArray) = kronsum(A, LinearMap(B)) +Base.kronsum(A::AbstractArray, B::LinearMap) = kronsum(LinearMap(A), 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 + mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) + # yields (:LinearMap(A1), :LinearMap(A2), ..., LinearMap(A(k-1))) + + @eval Base.kronsum($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = + kron($(mapargs...), $(Symbol(:A,k)), As...) +end + +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) + +function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) + A, B = L.maps + ma, na = size(A) + mb, nb = size(B) + (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("kronecker product")) + X = reshape(x, (nb, na)) + Y = reshape(y, (nb, na)) + mul!(Y, X, convert(AbstractMatrix{eltype(A)}, transpose(A))) + mul!(Y, B, X, true, true) + return y +end + +LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) + +LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) From 8252f6810396b679ad2d17d317645ef980a6edce Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 7 Sep 2019 14:53:01 +0200 Subject: [PATCH 14/33] add KroneckerPowerMap, yet untested --- src/kronecker.jl | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/kronecker.jl b/src/kronecker.jl index 54caff67..59735cb9 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -188,3 +188,46 @@ end LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) + +struct KroneckerPowerMap{T, TA<:LinearMap, N} <: LinearMap{T} + lmap::TA + k::Int + function KroneckerPowerMap(A::TA, k::Int) where {TA<:LinearMap} + k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return new{eltype(A),TA,k}(A, k) + end +end + +""" + kron(A::LinearMap, k::Int) +Construct a `KroneckerPowerMap<:LinearMap` object, a lazy representation of the +`k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` for `k≥2`. +""" +kron(A::LinearMap, k::Int) = KroneckerPowerMap(A, k) + +Base.size(A::KroneckerPowerMap) = size(A.lmap).^A.k + +LinearAlgebra.adjoint(A::KroneckerPowerMap) = KroneckerPowerMap(adjoint(A.lmap), A.k) +LinearAlgebra.transpose(A::KroneckerPowerMap) = KroneckerPowerMap(transpose(A.lmap), A.k) + +function Base.:*(A::KroneckerPowerMap{TA, <:LinearMap, N}, B::KroneckerPowerMap{TB, <:LinearMap, N}) where {TA, TB, N} + if size(A, 2) == size(B, 1) + return KroneckerPowerMap(A.lmap * B.lmap, N) + else + return CompositeMap{promote_type(eltype(A), eltype(B))}((B, A)) + end +end + +getmaps(A::KroneckerPowerMap{T, <:LinearMap, N}) where {T, N} = (A.lmap, KroneckerPowerMap(A.lmap, A.k-1)) +getmaps(A::KroneckerPowerMap{T, <:LinearMap, 2}) where {T} = (A.lmap, A.lmap) + +function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerPowerMap, x::AbstractVector) + A, B = getmaps(L) + K = kron(A, B) + A_mul_B!(y, K, x) + return y +end + +LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) + +LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) From 0fcbde42a759774ae0484aa3c49052300d5bba1d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 9 Sep 2019 11:15:41 +0200 Subject: [PATCH 15/33] KroneckerMap with many maps, optimize mul, more tests --- src/kronecker.jl | 161 +++++++++++++++++++++++----------------------- test/kronecker.jl | 74 +++++++++++++-------- 2 files changed, 127 insertions(+), 108 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 59735cb9..144daf8d 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -1,4 +1,4 @@ -struct KroneckerMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T} +struct KroneckerMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T} maps::As function KroneckerMap{T, As}(maps::As) where {T, As} for A in maps @@ -8,7 +8,7 @@ struct KroneckerMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T} end end -KroneckerMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerMap{T, As}(maps) +KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap{T, As}(maps) """ kron(A::LinearMap, B::LinearMap) @@ -41,7 +41,9 @@ julia> Matrix(Δ) 0 1 1 -4 ``` """ -Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) +# Base.kron(A::LinearMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B)) +# Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) +Base.kron(As::LinearMap...) = KroneckerMap{promote_type(map(eltype, As)...)}(tuple(As...)) Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), B) # promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product @@ -51,14 +53,20 @@ for k in 3:8 # is 8 sufficient? L = :($(Symbol(:A,k))::LinearMap) # yields :Ak mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) - # yields (:LinearMap(A1), :LinearMap(A2), ..., LinearMap(A(k-1))) + # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) @eval Base.kron($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = - kron($(mapargs...), $(Symbol(:A,k)), As...) + kron($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) end -Base.size(A::KroneckerMap, i) = prod(size.(A.maps, i)) -Base.size(A::KroneckerMap) = (size(A, 1), size(A, 2)) +promote_to_lmaps_(A::AbstractMatrix) = LinearMap(A) +promote_to_lmaps_(A::LinearMap) = A +promote_to_lmaps() = () +promote_to_lmaps(A) = (promote_to_lmaps_(A),) +@inline promote_to_lmaps(A, B, Cs...) = + (promote_to_lmaps_(A), promote_to_lmaps_(B), promote_to_lmaps(Cs...)...) + +Base.size(A::KroneckerMap) = (prod(size.(A.maps, 1)), prod(size.(A.maps, 2))) LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = all(issymmetric, A.maps) @@ -68,8 +76,8 @@ LinearAlgebra.adjoint(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(adjoin LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(transpose, A.maps)) function Base.:(*)(A::KroneckerMap, B::KroneckerMap) - if (size(A.maps[1], 2) == size(B.maps[1], 1) && size(A.maps[2], 2) == size(B.maps[2], 1)) - return kron(A.maps[1]*B.maps[1], A.maps[2]*B.maps[2]) + if length(A.maps) == length(B.maps) && all(M -> size(M[1], 2) == size(M[2], 1), zip(A.maps, B.maps)) + return kron(map(prod, zip(A.maps, B.maps))...) else return CompositeMap{T}(tuple(B, A)) end @@ -77,28 +85,58 @@ end Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap, x::AbstractVector) +function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} # require_one_based_indexing(y) + (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A, B = L.maps - ma, na = size(A) + 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 +function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} + # require_one_based_indexing(y) + (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) + 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 + +function _kronmul!(y, B, X, At, T) + na, ma = size(At) mb, nb = size(B) - (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("kronecker product")) - X = LinearMap(reshape(x, (nb, na))) - M = B * X * transpose(A) - T = eltype(L) 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) - mul!(y[((i-1)*mb+1):i*mb], M, v) + 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 +# function _kronmul!(y, B::MatrixMap, X, At::MatrixMap, T) +# na, ma = size(At) +# mb, nb = size(B) +# if (nb + ma) * na < (ma + mb) * nb +# mul!(reshape(y, (size(B, 1), size(At, 2))), B.lmap, X*At.lmap) +# else +# mul!(reshape(y, (size(B, 1), size(At, 2))), B.lmap*X, At.lmap) +# end +# return y +# end LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) LinearMaps.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} @@ -114,51 +152,44 @@ KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = Kroneck """ kronsum(A::LinearMap, B::LinearMap) + kronsum(A, B, Cs...) Construct a `KroneckerSumMap <: LinearMap` object, a (lazy) representation of the Kronecker sum `A⊕B = kron(A, Ib) + kron(Ia, B)` of two square `LinearMap`s. Here, -`Ia` and `Ib` are identity One of the two factors can be an `AbstractMatrix` -object, which then is promoted to `LinearMap` automatically. To avoid fallback to -the generic [`Base.kron`](@ref), there must be a `LinearMap` object among the -first 8 arguments in usage like `kron(A, B, Cs...)`. - -Note: 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)` -holds. Upon multiplication of Kronecker products, this rule is checked for -applicability, which leads to type-instability with a union of two types. +`Ia` and `Ib` are identity operators of the size of `A` and `B`, respectively. +Arguments of type `AbstractMatrix` are promoted to `LinearMap`s as long as there +is a `LinearMap` object among the first 8 arguments. # 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> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I); -julia> Δ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator +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 +julia> Δ₂ = LinearMaps.kronsum(D, D); # same operator as Kronecker sum + +julia> Matrix(Δ₁) == Matrix(Δ₂) +true ``` """ -Base.kronsum(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) -Base.kronsum(A::LinearMap, Bs::LinearMap...) = kronsum(A, kronsum(Bs...)) -Base.kronsum(A::LinearMap, B::AbstractArray) = kronsum(A, LinearMap(B)) -Base.kronsum(A::AbstractArray, B::LinearMap) = kronsum(LinearMap(A), B) -# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product -for k in 3:8 # is 8 sufficient? +kronsum(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerSumMap{promote_type(TA,TB)}((A, B)) +kronsum(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) = kronsum(A, kronsum(B, C, Ds...)) +# kronsum(A::LinearMap, B::AbstractArray) = kronsum(A, LinearMap(B)) +# kronsum(A::AbstractArray, B::LinearMap) = kronsum(LinearMap(A), B) +# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker sum +for k in 1: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 mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) - # yields (:LinearMap(A1), :LinearMap(A2), ..., LinearMap(A(k-1))) + # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) - @eval Base.kronsum($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = - kron($(mapargs...), $(Symbol(:A,k)), As...) + @eval kronsum($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = + kronsum($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) end Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i)) @@ -180,7 +211,7 @@ function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractV (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("kronecker product")) X = reshape(x, (nb, na)) Y = reshape(y, (nb, na)) - mul!(Y, X, convert(AbstractMatrix{eltype(A)}, transpose(A))) + mul!(Y, X, convert(AbstractMatrix, transpose(A))) mul!(Y, B, X, true, true) return y end @@ -189,45 +220,13 @@ LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) -struct KroneckerPowerMap{T, TA<:LinearMap, N} <: LinearMap{T} - lmap::TA - k::Int - function KroneckerPowerMap(A::TA, k::Int) where {TA<:LinearMap} - k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) - return new{eltype(A),TA,k}(A, k) - end -end - """ kron(A::LinearMap, k::Int) -Construct a `KroneckerPowerMap<:LinearMap` object, a lazy representation of the -`k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` for `k≥2`. -""" -kron(A::LinearMap, k::Int) = KroneckerPowerMap(A, k) - -Base.size(A::KroneckerPowerMap) = size(A.lmap).^A.k - -LinearAlgebra.adjoint(A::KroneckerPowerMap) = KroneckerPowerMap(adjoint(A.lmap), A.k) -LinearAlgebra.transpose(A::KroneckerPowerMap) = KroneckerPowerMap(transpose(A.lmap), A.k) -function Base.:*(A::KroneckerPowerMap{TA, <:LinearMap, N}, B::KroneckerPowerMap{TB, <:LinearMap, N}) where {TA, TB, N} - if size(A, 2) == size(B, 1) - return KroneckerPowerMap(A.lmap * B.lmap, N) - else - return CompositeMap{promote_type(eltype(A), eltype(B))}((B, A)) - end -end - -getmaps(A::KroneckerPowerMap{T, <:LinearMap, N}) where {T, N} = (A.lmap, KroneckerPowerMap(A.lmap, A.k-1)) -getmaps(A::KroneckerPowerMap{T, <:LinearMap, 2}) where {T} = (A.lmap, A.lmap) - -function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerPowerMap, x::AbstractVector) - A, B = getmaps(L) - K = kron(A, B) - A_mul_B!(y, K, x) - return y +Construct a lazy representation of the `k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` +for `k≥2`. This function is currently not type-stable! +""" +function Base.kron(A::LinearMap, k::Int) + k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return kron(Base.fill_to_length((), A, Val(k))...) end - -LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) - -LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) diff --git a/test/kronecker.jl b/test/kronecker.jl index 1b824e66..216f56bd 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -1,32 +1,52 @@ using Test, LinearMaps, LinearAlgebra -@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) +@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(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) ≈ Matrix(kron(LA, 3)) + K = @inferred kron(A, A, A, LA) + @test K isa LinearMaps.KroneckerMap + @test Matrix(K) ≈ kron(A, A, A, A) + @test K*K isa LinearMaps.KroneckerMap{ComplexF64,<:Tuple{Vararg{LinearMaps.CompositeMap}}} + 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 kron(LA, B) == kron(LA, LB) == kron(A, LB) + @test ishermitian(kron(LA'LA, LB'LB)) + A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B) + @test issymmetric(kron(LA'LA, LB'LB)) + @test ishermitian(kron(LA'LA, LB'LB)) 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)) + @testset "Kronecker sum" begin + A = rand(ComplexF64, 3, 3) + B = rand(ComplexF64, 2, 2) + LA = LinearMap(A) + LB = LinearMap(B) + KS = @inferred LinearMaps.kronsum(LA, B) + KSmat = kron(A, Matrix(I, 2, 2)) + kron(Matrix(I, 3, 3), B) + @test_skip 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_skip Matrix(transform(KS)) ≈ transform(Matrix(KS)) ≈ transform(KSmat) + @test_skip Matrix(LinearMaps.kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) + @test_skip Matrix(transform(LinearMap(LinearMaps.kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) + end + @inferred LinearMaps.kronsum(A, A, LB) end - @test kron(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) - K = @inferred kron(A, A, A, LA) - @test K isa LinearMaps.KroneckerMap - @test Matrix(K) ≈ kron(A, A, A, A) - K4 = kron(A, B, B, LB) - @test K4.maps[1].maps[1].maps[1].lmap === A - @test @inferred kron(LA, LB)' == @inferred kron(LA', LB') - @test kron(LA, B) == kron(LA, LB) == kron(A, LB) - @test ishermitian(kron(LA'LA, LB'LB)) - A = rand(3, 3); B = rand(2, 2); LA = LinearMap(A); LB = LinearMap(B) - @test issymmetric(kron(LA'LA, LB'LB)) - @test ishermitian(kron(LA'LA, LB'LB)) end From 0a8d31b67eabee672a20363766b4a842fdacc699 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 11 Sep 2019 13:28:51 +0200 Subject: [PATCH 16/33] add dim-checker helper, resolve type instability, tests --- src/composition.jl | 13 ++++++++----- src/kronecker.jl | 44 +++++++++++++++++++++++++++++++------------- test/kronecker.jl | 17 ++++++++++++----- 3 files changed, 51 insertions(+), 23 deletions(-) diff --git a/src/composition.jl b/src/composition.jl index c32553dc..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()) @@ -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/kronecker.jl b/src/kronecker.jl index 144daf8d..ea42179c 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -41,11 +41,13 @@ julia> Matrix(Δ) 0 1 1 -4 ``` """ -# Base.kron(A::LinearMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B)) -# Base.kron(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerMap{promote_type(TA,TB)}((A, B)) -Base.kron(As::LinearMap...) = KroneckerMap{promote_type(map(eltype, As)...)}(tuple(As...)) -Base.kron(A::LinearMap, B::AbstractArray) = kron(A, LinearMap(B)) -Base.kron(A::AbstractArray, B::LinearMap) = kron(LinearMap(A), B) +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)) @@ -75,14 +77,6 @@ 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)) -function Base.:(*)(A::KroneckerMap, B::KroneckerMap) - if length(A.maps) == length(B.maps) && all(M -> size(M[1], 2) == size(M[2], 1), zip(A.maps, B.maps)) - return kron(map(prod, zip(A.maps, B.maps))...) - else - return CompositeMap{T}(tuple(B, A)) - end -end - Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} @@ -102,6 +96,30 @@ function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractV _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ᵣ) +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(prod, zip(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ᵣ) +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 function _kronmul!(y, B, X, At, T) na, ma = size(At) diff --git a/test/kronecker.jl b/test/kronecker.jl index 216f56bd..4d6cbef1 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -22,16 +22,23 @@ using Test, LinearMaps, LinearAlgebra K = @inferred kron(A, A, A, LA) @test K isa LinearMaps.KroneckerMap @test Matrix(K) ≈ kron(A, A, A, A) - @test K*K isa LinearMaps.KroneckerMap{ComplexF64,<:Tuple{Vararg{LinearMaps.CompositeMap}}} + @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 kron(LA, B) == kron(LA, LB) == kron(A, LB) - @test ishermitian(kron(LA'LA, LB'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 issymmetric(kron(LA'LA, LB'LB)) - @test ishermitian(kron(LA'LA, LB'LB)) + @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) end @testset "Kronecker sum" begin A = rand(ComplexF64, 3, 3) From ffba20d92b16ac839d9df9d67e63bf3f07bec322 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 26 Sep 2019 22:48:11 +0200 Subject: [PATCH 17/33] add conversion rules --- src/LinearMaps.jl | 45 +---------------------------- src/conversion.jl | 73 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+), 44 deletions(-) create mode 100644 src/conversion.jl diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index ad145dfd..5eaedb49 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -68,50 +68,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 @@ -120,6 +76,7 @@ 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...) diff --git a/src/conversion.jl b/src/conversion.jl new file mode 100644 index 00000000..4a1478f4 --- /dev/null +++ b/src/conversion.jl @@ -0,0 +1,73 @@ +# 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::MatrixMap) = convert(AbstractMatrix, A.lmap) +Base.convert(::Type{Matrix}, A::MatrixMap) = 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 + +# 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) From 3bed1dc52aecdfab50aac86eba939262a7bcc756 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 26 Sep 2019 22:49:10 +0200 Subject: [PATCH 18/33] optimize MatrixMap and scaling multiplication, enable KroneckerSum tests --- src/kronecker.jl | 22 +++++++++++----------- src/uniformscalingmap.jl | 14 +++++++++++++- src/wrappedmap.jl | 3 +++ test/kronecker.jl | 8 ++++---- 4 files changed, 31 insertions(+), 16 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index ea42179c..80f7b8e2 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -83,7 +83,7 @@ function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,Lin # require_one_based_indexing(y) (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A, B = L.maps - X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false) + X = reshape(x, (size(B, 2), size(A, 2))) _kronmul!(y, B, X, transpose(A), T) return y end @@ -137,16 +137,16 @@ function _kronmul!(y, B, X, At, T) end return y end -# function _kronmul!(y, B::MatrixMap, X, At::MatrixMap, T) -# na, ma = size(At) -# mb, nb = size(B) -# if (nb + ma) * na < (ma + mb) * nb -# mul!(reshape(y, (size(B, 1), size(At, 2))), B.lmap, X*At.lmap) -# else -# mul!(reshape(y, (size(B, 1), size(At, 2))), B.lmap*X, At.lmap) -# end -# return y -# end +function _kronmul!(y, B::FreeMap, X, At::FreeMap, T) + na, ma = size(At) + mb, nb = size(B) + if (nb + ma) * na < (ma + mb) * nb + 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 LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 2f2f981e..79439bf4 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -61,8 +61,20 @@ function A_mul_B!(y::AbstractVector, A::UniformScalingMap, x::AbstractVector) end 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 iszero(β) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 5f919be5..3015ec93 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -35,6 +35,9 @@ LinearAlgebra.mul!(y::AbstractVector, A::WrappedMap, x::AbstractVector, α::Numb end # VERSION +LinearAlgebra.mul!(Y::AbstractMatrix, A::MatrixMap, X::AbstractMatrix, α::Number=true, β::Number=false) = + mul!(Y, A.lmap, X, α, β) + # properties Base.size(A::WrappedMap) = size(A.lmap) LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric diff --git a/test/kronecker.jl b/test/kronecker.jl index 4d6cbef1..e82c7c61 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -47,12 +47,12 @@ using Test, LinearMaps, LinearAlgebra LB = LinearMap(B) KS = @inferred LinearMaps.kronsum(LA, B) KSmat = kron(A, Matrix(I, 2, 2)) + kron(Matrix(I, 3, 3), B) - @test_skip Matrix(KS) ≈ Matrix(kron(A, LinearMap(I, 2)) + kron(LinearMap(I, 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_skip Matrix(transform(KS)) ≈ transform(Matrix(KS)) ≈ transform(KSmat) - @test_skip Matrix(LinearMaps.kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) - @test_skip Matrix(transform(LinearMap(LinearMaps.kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) + @test Matrix(transform(KS)) ≈ transform(Matrix(KS)) ≈ transform(KSmat) + @test Matrix(LinearMaps.kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) + @test Matrix(transform(LinearMap(LinearMaps.kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) end @inferred LinearMaps.kronsum(A, A, LB) end From c6fe9cf51ad2f058e660750aa1a3022b495bff11 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 27 Sep 2019 09:50:51 +0200 Subject: [PATCH 19/33] fix matrix map handling across different julia versions --- src/kronecker.jl | 2 +- src/wrappedmap.jl | 9 +++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 80f7b8e2..089fa9c8 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -137,7 +137,7 @@ function _kronmul!(y, B, X, At, T) end return y end -function _kronmul!(y, B::FreeMap, X, At::FreeMap, T) +function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T) na, ma = size(At) mb, nb = size(B) if (nb + ma) * na < (ma + mb) * nb diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 3015ec93..ff6b73ab 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -33,11 +33,16 @@ 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, α, β) -end # VERSION - 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 Base.size(A::WrappedMap) = size(A.lmap) LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric From 649bde50960316f891b386dc9a42a93a2b5305e5 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 27 Sep 2019 15:09:56 +0200 Subject: [PATCH 20/33] minor fixes, improve coverage --- README.md | 8 +++++--- src/LinearMaps.jl | 1 + src/kronecker.jl | 15 +++++++++++++-- src/uniformscalingmap.jl | 10 ++++++---- test/conversion.jl | 38 ++++++++++++++++++++++++++++++++++++++ test/linearmaps.jl | 20 -------------------- test/runtests.jl | 2 ++ test/uniformscalingmap.jl | 4 +++- 8 files changed, 68 insertions(+), 30 deletions(-) create mode 100644 test/conversion.jl diff --git a/README.md b/README.md index 6b56b92b..b265de9e 100644 --- a/README.md +++ b/README.md @@ -8,9 +8,11 @@ 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 of `LinearMap`s. `AbstractMatrix`es are - promoted to `LinearMap`s if one of the first 8 Kronecker factors is a - `LinearMap` object. +* 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, diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 5eaedb49..954158e4 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -1,6 +1,7 @@ module LinearMaps export LinearMap +export kronsum using LinearAlgebra using SparseArrays diff --git a/src/kronecker.jl b/src/kronecker.jl index 089fa9c8..fcf3f35a 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -79,7 +79,7 @@ LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(tran Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} # require_one_based_indexing(y) (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A, B = L.maps @@ -87,7 +87,7 @@ function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,Lin _kronmul!(y, B, X, transpose(A), T) return y end -function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} +function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} # require_one_based_indexing(y) (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A = first(L.maps) @@ -248,3 +248,14 @@ function Base.kron(A::LinearMap, k::Int) k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) return kron(Base.fill_to_length((), A, Val(k))...) end + +""" + kronsum(A::LinearMap, k::Int) + +Construct a lazy representation of the `k`-th Kronecker sum power `A ⊕ A ⊕ ... ⊕ A` +for `k≥2`. This function is currently not type-stable! +""" +function kronsum(A::LinearMap, k::Int) + k > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k")) + return kronsum(Base.fill_to_length((), A, Val(k))...) +end diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 79439bf4..6706f9e4 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -38,7 +38,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 @@ -52,7 +52,7 @@ 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) + return fill!(y, zero(eltype(y))) elseif isone(λ) return copyto!(y, x) else @@ -71,14 +71,16 @@ end @boundscheck size(X) == size(Y) || throw(DimensionMismatch("mul!")) @boundscheck size(X,1) == J.M || throw(DimensionMismatch("mul!")) _scaling!(Y, J, X, α, β) - return y + return Y end function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false) λ = J.λ @inbounds 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/test/conversion.jl b/test/conversion.jl new file mode 100644 index 00000000..42a47e82 --- /dev/null +++ b/test/conversion.jl @@ -0,0 +1,38 @@ +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 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/linearmaps.jl b/test/linearmaps.jl index 47f8cc48..7570da02 100644 --- a/test/linearmaps.jl +++ b/test/linearmaps.jl @@ -19,26 +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) - - J = LinearMap(I, 10) - @test J isa LinearMap{Bool} - @test sparse(J) == Matrix{Bool}(I, 10, 10) - @test nnz(sparse(J)) == 10 end Av = A * v diff --git a/test/runtests.jl b/test/runtests.jl index 64fc208d..6851f7b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,3 +19,5 @@ include("numbertypes.jl") include("blockmap.jl") include("kronecker.jl") + +include("conversion.jl") diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index 140592c9..d51090e3 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) @@ -40,4 +40,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 From c9d654788a8d0146d746a7167c263ccb441c2eeb Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 27 Sep 2019 17:28:42 +0200 Subject: [PATCH 21/33] improve coverage, minor fixes --- src/uniformscalingmap.jl | 17 ++++------------- test/kronecker.jl | 31 ++++++++++++++++++------------- test/uniformscalingmap.jl | 5 +++-- 3 files changed, 25 insertions(+), 28 deletions(-) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index 6706f9e4..adda50f7 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -12,7 +12,6 @@ UniformScalingMap(λ::T, sz::Dims{2}) where {T} = (sz[1] == sz[2] ? UniformScalingMap(λ, sz[1]) : error("UniformScalingMap needs to be square")) # properties -Base.size(A::UniformScalingMap, n) = (n==1 || n==2 ? A.M : error("LinearMap objects have only 2 dimensions")) Base.size(A::UniformScalingMap) = (A.M, A.M) Base.isreal(A::UniformScalingMap) = isreal(A.λ) LinearAlgebra.issymmetric(::UniformScalingMap) = true @@ -48,17 +47,9 @@ 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, zero(eltype(y))) - 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 @inline function LinearAlgebra.mul!(y::AbstractVector, J::UniformScalingMap, x::AbstractVector, α::Number=true, β::Number=false) @@ -76,7 +67,7 @@ end function _scaling!(y, J::UniformScalingMap, x, α::Number=true, β::Number=false) λ = J.λ - @inbounds if isone(α) + if isone(α) if iszero(β) iszero(λ) && return fill!(y, zero(eltype(y))) isone(λ) && return copyto!(y, x) diff --git a/test/kronecker.jl b/test/kronecker.jl index e82c7c61..d24fd417 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -39,21 +39,26 @@ using Test, LinearMaps, LinearAlgebra 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 - A = rand(ComplexF64, 3, 3) - B = rand(ComplexF64, 2, 2) - LA = LinearMap(A) - LB = LinearMap(B) - KS = @inferred LinearMaps.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(LinearMaps.kronsum(transform(LA), transform(LB))) ≈ transform(KSmat) - @test Matrix(transform(LinearMap(LinearMaps.kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) + 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(kronsum(LA, 3)) ≈ Matrix(kronsum(LA, A, A)) end - @inferred LinearMaps.kronsum(A, A, LB) end end diff --git a/test/uniformscalingmap.jl b/test/uniformscalingmap.jl index d51090e3..ab38ce88 100644 --- a/test/uniformscalingmap.jl +++ b/test/uniformscalingmap.jl @@ -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, $α, $β) From e701ff52ab5c25ffd41823b2a9b3c83deb59f91e Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 27 Sep 2019 18:34:32 +0200 Subject: [PATCH 22/33] improve coverage further --- src/kronecker.jl | 2 -- test/conversion.jl | 3 +++ test/kronecker.jl | 3 +++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index fcf3f35a..9555a21c 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -195,8 +195,6 @@ true """ kronsum(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerSumMap{promote_type(TA,TB)}((A, B)) kronsum(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) = kronsum(A, kronsum(B, C, Ds...)) -# kronsum(A::LinearMap, B::AbstractArray) = kronsum(A, LinearMap(B)) -# kronsum(A::AbstractArray, B::LinearMap) = kronsum(LinearMap(A), B) # promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker sum for k in 1:8 # is 8 sufficient? Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1)) diff --git a/test/conversion.jl b/test/conversion.jl index 42a47e82..f4fb92aa 100644 --- a/test/conversion.jl +++ b/test/conversion.jl @@ -17,6 +17,9 @@ using Test, LinearMaps, LinearAlgebra, SparseArrays @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 diff --git a/test/kronecker.jl b/test/kronecker.jl index d24fd417..09b199f1 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -18,6 +18,8 @@ using Test, LinearMaps, LinearAlgebra @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(kron(LA, 3)) K = @inferred kron(A, A, A, LA) @test K isa LinearMaps.KroneckerMap @@ -59,6 +61,7 @@ using Test, LinearMaps, LinearAlgebra end @inferred kronsum(A, A, LB) @test Matrix(kronsum(LA, 3)) ≈ Matrix(kronsum(LA, A, A)) + @test kronsum(LA, LA, LB) == kronsum(LA, kronsum(LA, LB)) end end end From 576304afccfc0289c7a63710d81d98e4bae1dd1c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 1 Oct 2019 08:51:58 +0200 Subject: [PATCH 23/33] remove false type-instability comment, uniform scaling comparison --- src/kronecker.jl | 7 +++---- src/uniformscalingmap.jl | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 9555a21c..08a04df6 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -19,10 +19,9 @@ object, which then is promoted to `LinearMap` automatically. To avoid fallback t the generic [`Base.kron`](@ref), there must be a `LinearMap` object among the first 8 arguments in usage like `kron(A, B, Cs...)`. -Note: 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)` -holds. Upon multiplication of Kronecker products, this rule is checked for -applicability, which leads to type-instability with a union of two types. +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) diff --git a/src/uniformscalingmap.jl b/src/uniformscalingmap.jl index adda50f7..25461858 100644 --- a/src/uniformscalingmap.jl +++ b/src/uniformscalingmap.jl @@ -18,8 +18,8 @@ LinearAlgebra.issymmetric(::UniformScalingMap) = true LinearAlgebra.ishermitian(A::UniformScalingMap) = isreal(A) LinearAlgebra.isposdef(A::UniformScalingMap) = isposdef(A.λ) -# comparison of UniformScalingMap objects, sufficient but not necessary -Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (eltype(A) == eltype(B) && A.λ == B.λ && A.M == B.M) +# comparison of UniformScalingMap objects +Base.:(==)(A::UniformScalingMap, B::UniformScalingMap) = (A.λ == B.λ && A.M == B.M) # special transposition behavior LinearAlgebra.transpose(A::UniformScalingMap) = A From 7e1de301ae21baf15fde25e27f6857b20f8c6f24 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 2 Oct 2019 13:23:55 +0200 Subject: [PATCH 24/33] require one-based indexing, move kronecker power, simplify size --- src/kronecker.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 08a04df6..07474672 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -67,7 +67,18 @@ promote_to_lmaps(A) = (promote_to_lmaps_(A),) @inline promote_to_lmaps(A, B, Cs...) = (promote_to_lmaps_(A), promote_to_lmaps_(B), promote_to_lmaps(Cs...)...) -Base.size(A::KroneckerMap) = (prod(size.(A.maps, 1)), prod(size.(A.maps, 2))) +""" + kron(A::LinearMap, k::Int) + +Construct a lazy representation of the `k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` +for `k≥2`. This function is currently not type-stable! +""" +@inline function Base.kron(A::LinearMap, k::Int) + k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return kron(Base.fill_to_length((), A, Val(k))...) +end + +Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps) LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = all(issymmetric, A.maps) @@ -79,7 +90,7 @@ LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(tran Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} - # require_one_based_indexing(y) + require_one_based_indexing(y) (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A, B = L.maps X = reshape(x, (size(B, 2), size(A, 2))) @@ -87,7 +98,7 @@ function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x return y end function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} - # require_one_based_indexing(y) + require_one_based_indexing(y) (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) A = first(L.maps) B = kron(Base.tail(L.maps)...) @@ -100,7 +111,7 @@ end 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(prod, zip(A.maps, B.maps))...), x) + A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x) else A_mul_B!(y, LinearMap(A)*B, x) end @@ -235,17 +246,6 @@ LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) -""" - kron(A::LinearMap, k::Int) - -Construct a lazy representation of the `k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` -for `k≥2`. This function is currently not type-stable! -""" -function Base.kron(A::LinearMap, k::Int) - k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) - return kron(Base.fill_to_length((), A, Val(k))...) -end - """ kronsum(A::LinearMap, k::Int) From 63b558f78244d908da9c9bd40035ee064b287d30 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 4 Oct 2019 11:19:09 +0200 Subject: [PATCH 25/33] clean up Kronecker (sum) powers --- src/LinearMaps.jl | 2 +- src/kronecker.jl | 76 +++++++++++++++++++++++++++++++---------------- test/kronecker.jl | 4 +-- 3 files changed, 53 insertions(+), 29 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 954158e4..4915699f 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -1,7 +1,7 @@ module LinearMaps export LinearMap -export kronsum +export ⊗, kronsum, ⊕ using LinearAlgebra using SparseArrays diff --git a/src/kronecker.jl b/src/kronecker.jl index 07474672..9d1cf478 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -14,10 +14,12 @@ KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap kron(A::LinearMap, B::LinearMap) 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` -object, which then is promoted to `LinearMap` automatically. To avoid fallback to +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), there must be a `LinearMap` object among the -first 8 arguments in usage like `kron(A, B, Cs...)`. +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`. 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)` @@ -67,16 +69,25 @@ promote_to_lmaps(A) = (promote_to_lmaps_(A),) @inline promote_to_lmaps(A, B, Cs...) = (promote_to_lmaps_(A), promote_to_lmaps_(B), promote_to_lmaps(Cs...)...) +struct KronPower{T<:Integer} + p::T + function KronPower(p::Integer) + p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return new{typeof(p)}(p) + end +end + """ - kron(A::LinearMap, k::Int) + ⊗(k::Integer) -Construct a lazy representation of the `k`-th Kronecker power `A ⊗ A ⊗ ... ⊗ A` -for `k≥2`. This function is currently not type-stable! +Construct a lazy representation of the `k`-th Kronecker power `A^⊗(k) = A ⊗ A ⊗ ... ⊗ A`, +where `A` can be an `AbstractMatrix` or a `LinearMap`. """ -@inline function Base.kron(A::LinearMap, k::Int) - k > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) - return kron(Base.fill_to_length((), A, Val(k))...) -end +⊗(k::Integer) = KronPower(k) + +⊗(a, b, c...) = kron(a, b, c...) + +Base.:(^)(a::Union{LinearMap,AbstractMatrix}, p::KronPower) = kron(ntuple(n->promote_to_lmaps_(a), p.p)...) Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) @@ -183,10 +194,12 @@ KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = Kroneck kronsum(A, B, Cs...) Construct a `KroneckerSumMap <: LinearMap` object, a (lazy) representation of the -Kronecker sum `A⊕B = kron(A, Ib) + kron(Ia, B)` of two square `LinearMap`s. Here, +Kronecker sum `A⊕B = A ⊗ Ib + Ia ⊗ B` of two square `LinearMap`s. Here, `Ia` and `Ib` are identity operators of the size of `A` and `B`, respectively. Arguments of type `AbstractMatrix` are promoted to `LinearMap`s as long as there -is a `LinearMap` object among the first 8 arguments. +is a `LinearMap` object among the first 8 arguments. 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) @@ -195,11 +208,13 @@ 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 +julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum -julia> Δ₂ = LinearMaps.kronsum(D, D); # same operator as Kronecker sum +julia> Δ₂ = kronsum(D, D); -julia> Matrix(Δ₁) == Matrix(Δ₂) +julia> Δ₃ = D^⊕(2); + +julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃) true ``` """ @@ -218,6 +233,26 @@ for k in 1:8 # is 8 sufficient? kronsum($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) end +struct KronSumPower{T<:Integer} + p::T + function KronSumPower(p::Integer) + p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) + return new{typeof(p)}(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::Union{LinearMap,AbstractMatrix}, p::KronSumPower) = kronsum(ntuple(n->promote_to_lmaps_(a), p.p)...) + Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i)) Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2)) @@ -245,14 +280,3 @@ end LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) - -""" - kronsum(A::LinearMap, k::Int) - -Construct a lazy representation of the `k`-th Kronecker sum power `A ⊕ A ⊕ ... ⊕ A` -for `k≥2`. This function is currently not type-stable! -""" -function kronsum(A::LinearMap, k::Int) - k > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k")) - return kronsum(Base.fill_to_length((), A, Val(k))...) -end diff --git a/test/kronecker.jl b/test/kronecker.jl index 09b199f1..3723ea1e 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -20,7 +20,7 @@ using Test, LinearMaps, LinearAlgebra 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(kron(LA, 3)) + @test kron(A, A, A) ≈ Matrix(@inferred kron(LA, LA, LA)) ≈ Matrix(LA^⊗(3)) ≈ Matrix(A^⊗(3)) K = @inferred kron(A, A, A, LA) @test K isa LinearMaps.KroneckerMap @test Matrix(K) ≈ kron(A, A, A, A) @@ -60,7 +60,7 @@ using Test, LinearMaps, LinearAlgebra @test Matrix(transform(LinearMap(kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) end @inferred kronsum(A, A, LB) - @test Matrix(kronsum(LA, 3)) ≈ Matrix(kronsum(LA, A, A)) + @test Matrix(LA^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) @test kronsum(LA, LA, LB) == kronsum(LA, kronsum(LA, LB)) end end From 9c9f2cfb3f4746ba5594bd841e6a2b55a90899ea Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 6 Oct 2019 17:46:44 +0200 Subject: [PATCH 26/33] make kron(sum)powers inferrable --- src/kronecker.jl | 16 +++++++--------- test/kronecker.jl | 4 ++-- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 9d1cf478..f290066c 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -69,11 +69,10 @@ promote_to_lmaps(A) = (promote_to_lmaps_(A),) @inline promote_to_lmaps(A, B, Cs...) = (promote_to_lmaps_(A), promote_to_lmaps_(B), promote_to_lmaps(Cs...)...) -struct KronPower{T<:Integer} - p::T +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{typeof(p)}(p) + return new{p}() end end @@ -87,7 +86,7 @@ where `A` can be an `AbstractMatrix` or a `LinearMap`. ⊗(a, b, c...) = kron(a, b, c...) -Base.:(^)(a::Union{LinearMap,AbstractMatrix}, p::KronPower) = kron(ntuple(n->promote_to_lmaps_(a), p.p)...) +Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronPower{p}) where {p} = kron(Base.fill_to_length((), promote_to_lmaps_(A), Val(p))...) Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) @@ -233,11 +232,10 @@ for k in 1:8 # is 8 sufficient? kronsum($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) end -struct KronSumPower{T<:Integer} - p::T +struct KronSumPower{p} function KronSumPower(p::Integer) - p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) - return new{typeof(p)}(p) + p > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k")) + return new{p}() end end @@ -251,7 +249,7 @@ where `A` can be a square `AbstractMatrix` or a `LinearMap`. ⊕(a, b, c...) = kronsum(a, b, c...) -Base.:(^)(a::Union{LinearMap,AbstractMatrix}, p::KronSumPower) = kronsum(ntuple(n->promote_to_lmaps_(a), p.p)...) +Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronSumPower{p}) where {p} = kronsum(Base.fill_to_length((), promote_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)) diff --git a/test/kronecker.jl b/test/kronecker.jl index 3723ea1e..030c1eed 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -20,7 +20,7 @@ using Test, LinearMaps, LinearAlgebra 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(LA^⊗(3)) ≈ Matrix(A^⊗(3)) + @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) @@ -60,7 +60,7 @@ using Test, LinearMaps, LinearAlgebra @test Matrix(transform(LinearMap(kronsum(LA, LB)))) ≈ Matrix(transform(KS)) ≈ transform(KSmat) end @inferred kronsum(A, A, LB) - @test Matrix(LA^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) + @test Matrix(@inferred LA^⊕(3)) == Matrix(@inferred A^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) @test kronsum(LA, LA, LB) == kronsum(LA, kronsum(LA, LB)) end end From 5362ab23d9ff2c08593fdd55ea7ac10e696e80d1 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 16 Oct 2019 11:44:18 +0200 Subject: [PATCH 27/33] remove obsolete specifiers --- src/kronecker.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index f290066c..e7aab5aa 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -168,9 +168,9 @@ function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{Matrix return y end -LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) ############### # KroneckerSumMap @@ -263,7 +263,7 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) +function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) A, B = L.maps ma, na = size(A) mb, nb = size(B) @@ -275,6 +275,6 @@ function LinearMaps.A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractV return y end -LinearMaps.At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) -LinearMaps.Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) From f438e0dc92db23047c99087d5ebedbfeee3b9352 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 16 Oct 2019 12:46:04 +0200 Subject: [PATCH 28/33] rename promote_to_lmaps -> convert_to_lmaps to avoid conflict with BlockMaps --- src/kronecker.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index e7aab5aa..3f42cb0b 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -59,15 +59,15 @@ for k in 3:8 # is 8 sufficient? # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) @eval Base.kron($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = - kron($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) + kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) end -promote_to_lmaps_(A::AbstractMatrix) = LinearMap(A) -promote_to_lmaps_(A::LinearMap) = A -promote_to_lmaps() = () -promote_to_lmaps(A) = (promote_to_lmaps_(A),) -@inline promote_to_lmaps(A, B, Cs...) = - (promote_to_lmaps_(A), promote_to_lmaps_(B), promote_to_lmaps(Cs...)...) +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...)...) struct KronPower{p} function KronPower(p::Integer) @@ -86,7 +86,7 @@ where `A` can be an `AbstractMatrix` or a `LinearMap`. ⊗(a, b, c...) = kron(a, b, c...) -Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronPower{p}) where {p} = kron(Base.fill_to_length((), promote_to_lmaps_(A), Val(p))...) +Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronPower{p}) where {p} = kron(Base.fill_to_length((), convert_to_lmaps_(A), Val(p))...) Base.size(A::KroneckerMap) = map(*, size.(A.maps)...) @@ -229,7 +229,7 @@ for k in 1:8 # is 8 sufficient? # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) @eval kronsum($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = - kronsum($(mapargs...), $(Symbol(:A,k)), promote_to_lmaps(As...)...) + kronsum($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) end struct KronSumPower{p} @@ -249,7 +249,7 @@ where `A` can be a square `AbstractMatrix` or a `LinearMap`. ⊕(a, b, c...) = kronsum(a, b, c...) -Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronSumPower{p}) where {p} = kronsum(Base.fill_to_length((), promote_to_lmaps_(A), Val(p))...) +Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronSumPower{p}) where {p} = kronsum(Base.fill_to_length((), 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)) From d982a17ecff181a7efd8314ed3c6807b66cfcb9f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 21 Oct 2019 14:13:43 +0200 Subject: [PATCH 29/33] move conversion and check_dim to Main, pepper with at-propagate_inbounds --- src/LinearMaps.jl | 20 ++++++++++ src/kronecker.jl | 97 +++++++++++++++++++++++++---------------------- 2 files changed, 71 insertions(+), 46 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 4915699f..e8d273a8 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -27,6 +27,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!")) diff --git a/src/kronecker.jl b/src/kronecker.jl index 3f42cb0b..7e412262 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -62,13 +62,6 @@ for k in 3:8 # is 8 sufficient? kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) end -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...)...) - struct KronPower{p} function KronPower(p::Integer) p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k")) @@ -99,17 +92,52 @@ LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(tran Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +################# +# multiplication helper functions +################# + +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 +function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T) + na, ma = size(At) + mb, nb = size(B) + if (nb + ma) * na < (ma + mb) * nb + 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 +################# + +@inline function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} require_one_based_indexing(y) - (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) + @boundscheck check_dim_mul(y, L, x) A, B = L.maps X = reshape(x, (size(B, 2), size(A, 2))) _kronmul!(y, B, X, transpose(A), T) return y end -function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} +Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T} require_one_based_indexing(y) - (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("A_mul_B!")) + @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) @@ -118,7 +146,7 @@ function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) wher end # mixed-product rule, prefer the right if possible # (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ) -function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector) +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) @@ -128,7 +156,7 @@ function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap, end # mixed-product rule, prefer the right if possible # (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ) -function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T} +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) @@ -141,36 +169,11 @@ function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerM end end -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 -function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T) - na, ma = size(At) - mb, nb = size(B) - if (nb + ma) * na < (ma + mb) * nb - 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 - -At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = + A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) = + A_mul_B!(y, adjoint(A), x) ############### # KroneckerSumMap @@ -263,11 +266,11 @@ LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(ma Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps) -function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) +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) - (length(y) == size(L, 1) && length(x) == size(L, 2)) || throw(DimensionMismatch("kronecker product")) X = reshape(x, (nb, na)) Y = reshape(y, (nb, na)) mul!(Y, X, convert(AbstractMatrix, transpose(A))) @@ -275,6 +278,8 @@ function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector) return y end -At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, transpose(A), x) +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = + A_mul_B!(y, transpose(A), x) -Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = A_mul_B!(y, adjoint(A), x) +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) = + A_mul_B!(y, adjoint(A), x) From ff62abb60626546625141e7824229918d44f86aa Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 23 Oct 2019 09:45:50 +0200 Subject: [PATCH 30/33] make at-inbounds work properly, add some conversion rules --- src/conversion.jl | 9 ++++++--- src/kronecker.jl | 6 +++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/conversion.jl b/src/conversion.jl index 4a1478f4..f3e9a45e 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -11,7 +11,6 @@ function Base.Matrix(A::LinearMap) 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) @@ -19,8 +18,8 @@ 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::MatrixMap) = convert(AbstractMatrix, A.lmap) -Base.convert(::Type{Matrix}, A::MatrixMap) = convert(Matrix, A.lmap) +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)...) @@ -45,6 +44,8 @@ function Base.convert(::Type{Matrix}, Aλ::CompositeMap{<:Any,<:Tuple{UniformSca 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) @@ -71,3 +72,5 @@ function SparseArrays.sparse(A::LinearMap{T}) where {T} 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 index 7e412262..9ba26522 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -96,7 +96,7 @@ Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps # multiplication helper functions ################# -function _kronmul!(y, B, X, At, T) +@inline function _kronmul!(y, B, X, At, T) na, ma = size(At) mb, nb = size(B) v = zeros(T, ma) @@ -112,7 +112,7 @@ function _kronmul!(y, B, X, At, T) end return y end -function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T) +@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) * na < (ma + mb) * nb @@ -127,7 +127,7 @@ end # multiplication with vectors ################# -@inline function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T} +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 From eee5ae4b9139fd5b1614704d2161b6c10f767fe7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 10 Nov 2019 13:02:17 +0100 Subject: [PATCH 31/33] simplify kronsum constructor --- src/kronecker.jl | 20 ++++++-------------- test/kronecker.jl | 3 ++- 2 files changed, 8 insertions(+), 15 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 9ba26522..a3f1eecd 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -220,20 +220,12 @@ julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃) true ``` """ -kronsum(A::LinearMap{TA}, B::LinearMap{TB}) where {TA,TB} = KroneckerSumMap{promote_type(TA,TB)}((A, B)) -kronsum(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) = kronsum(A, kronsum(B, C, Ds...)) -# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker sum -for k in 1: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 - mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) - # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) - - @eval kronsum($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = - kronsum($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) -end +kronsum(A::Union{LinearMap,AbstractMatrix}, B::Union{LinearMap,AbstractMatrix}) = + KroneckerSumMap{promote_type(eltype(A), eltype(B))}(convert_to_lmaps(A, B)) +kronsum(As::Vararg{Union{LinearMap,AbstractMatrix}}) = + kronsum(kronsum(first(As)), kronsum(convert_to_lmaps(Base.tail(As)...)...)) +kronsum(A::Union{LinearMap,AbstractMatrix}) = convert_to_lmaps_(A) +kronsum() = () struct KronSumPower{p} function KronSumPower(p::Integer) diff --git a/test/kronecker.jl b/test/kronecker.jl index 030c1eed..564c5350 100644 --- a/test/kronecker.jl +++ b/test/kronecker.jl @@ -61,7 +61,8 @@ using Test, LinearMaps, LinearAlgebra end @inferred kronsum(A, A, LB) @test Matrix(@inferred LA^⊕(3)) == Matrix(@inferred A^⊕(3)) ≈ Matrix(kronsum(LA, A, A)) - @test kronsum(LA, LA, LB) == kronsum(LA, kronsum(LA, LB)) + @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 From 1b732ba850e0a5f7ee2643db1b0fe3f365bdec16 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 11 Nov 2019 10:19:21 +0100 Subject: [PATCH 32/33] edits as per review, update docstrings --- src/LinearMaps.jl | 6 ++++-- src/kronecker.jl | 54 +++++++++++++++++++++++++---------------------- src/wrappedmap.jl | 6 +++--- 3 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index e8d273a8..a8adf542 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -15,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 @@ -129,14 +131,14 @@ For functions `f`, there is one more keyword argument 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/kronecker.jl b/src/kronecker.jl index a3f1eecd..1a15b676 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -12,14 +12,18 @@ KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap """ 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), 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`. +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)` @@ -54,11 +58,11 @@ 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 + # 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::Union{LinearMap,AbstractMatrix}...) = + @eval Base.kron($(Is...), $L, As::MapOrMatrix...) = kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) end @@ -72,19 +76,20 @@ 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`. +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, c...) = kron(a, b, c...) +⊗(A, B, Cs...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(convert_to_lmaps(A, B, Cs...)) -Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronPower{p}) where {p} = kron(Base.fill_to_length((), convert_to_lmaps_(A), Val(p))...) +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}) = 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)) @@ -192,16 +197,17 @@ end KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerSumMap{T, As}(maps) """ - kronsum(A::LinearMap, B::LinearMap) + 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 `LinearMap`s. Here, -`Ia` and `Ib` are identity operators of the size of `A` and `B`, respectively. -Arguments of type `AbstractMatrix` are promoted to `LinearMap`s as long as there -is a `LinearMap` object among the first 8 arguments. For convenience, one can -also use `A ⊕ B` or `⊕(A, B, Cs...)` (typed as `\\oplus+TAB`) to construct the -`KroneckerSumMap`. +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) @@ -220,12 +226,10 @@ julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃) true ``` """ -kronsum(A::Union{LinearMap,AbstractMatrix}, B::Union{LinearMap,AbstractMatrix}) = +kronsum(A::MapOrMatrix, B::MapOrMatrix) = KroneckerSumMap{promote_type(eltype(A), eltype(B))}(convert_to_lmaps(A, B)) -kronsum(As::Vararg{Union{LinearMap,AbstractMatrix}}) = - kronsum(kronsum(first(As)), kronsum(convert_to_lmaps(Base.tail(As)...)...)) -kronsum(A::Union{LinearMap,AbstractMatrix}) = convert_to_lmaps_(A) -kronsum() = () +kronsum(A::MapOrMatrix, B::MapOrMatrix, C::MapOrMatrix, Ds::MapOrMatrix...) = + kronsum(A, kronsum(B, C, Ds...)) struct KronSumPower{p} function KronSumPower(p::Integer) @@ -244,7 +248,7 @@ where `A` can be a square `AbstractMatrix` or a `LinearMap`. ⊕(a, b, c...) = kronsum(a, b, c...) -Base.:(^)(A::Union{LinearMap,AbstractMatrix}, ::KronSumPower{p}) where {p} = kronsum(Base.fill_to_length((), convert_to_lmaps_(A), Val(p))...) +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)) diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index ff6b73ab..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} From fd10ea2494301b9915ef93d06df1a6cd8455f676 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 11 Nov 2019 12:44:37 +0100 Subject: [PATCH 33/33] use memory requirement heuristic for branching --- src/kronecker.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/kronecker.jl b/src/kronecker.jl index 1a15b676..17f9511e 100644 --- a/src/kronecker.jl +++ b/src/kronecker.jl @@ -120,7 +120,7 @@ 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) * na < (ma + mb) * nb + 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.λ) @@ -136,7 +136,7 @@ Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T, require_one_based_indexing(y) @boundscheck check_dim_mul(y, L, x) A, B = L.maps - X = reshape(x, (size(B, 2), size(A, 2))) + 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