Add linear_map between MatrixZonotope and SPZ

docs/src/lib/lazy_operations/LinearMap.md

@@ -6,7 +6,7 @@ CurrentModule = LazySets
 
 ```@docs
 LinearMap
-*(::Union{AbstractMatrix, UniformScaling, AbstractVector, Real}, ::LazySet)
+*(::Union{AbstractMatrix, UniformScaling, AbstractVector, Real, AbstractMatrixZonotope}, ::LazySet)
 dim(::LinearMap)
 ρ(::AbstractVector, ::LinearMap)
 σ(::AbstractVector, ::LinearMap)

docs/src/lib/sets/SparsePolynomialZonotope.md

@@ -17,6 +17,8 @@ genmat_indep(::SparsePolynomialZonotope)
 indexvector(::SparsePolynomialZonotope)
 polynomial_order(::SparsePolynomialZonotope)
 rand(::Type{SparsePolynomialZonotope})
+linear_map(::MatrixZonotope, ::SparsePolynomialZonotope)
+linear_map(::MatrixZonotopeProduct, ::SparsePolynomialZonotope)
 ```
 ```@meta
 CurrentModule = LazySets

src/Approximations/overapproximate_zonotope.jl

@@ -1251,3 +1251,68 @@ function load_overapproximate_ICP()
         end
     end
 end  # load_overapproximate_ICP()
+
+"""
+	overapproximate(lm::LinearMap{N,S,NM,MAT},
+                         ::Type{Zonotope}) where {N,S<:AbstractZonotope{N},NM,
+                                                  MAT<:MatrixZonotope{NM}}
+
+Overapproximate the linear map of a zonotope through a matrix zonotope,
+following of Proposition 4 of [AlthoffGCKH11](@citet).
+
+### Input
+
+- `lm` -- a linear map of a zonotope through a matrix zonotope
+
+### Output
+
+A zonotope overapproximating the linear map.
+
+"""
+function overapproximate(lm::LinearMap{N,S,NM,MAT},
+                         ::Type{<:Zonotope}) where {N,S<:AbstractZonotope{N},NM,
+                                                  MAT<:MatrixZonotope{NM}}
+    MZ = matrix(lm)
+    Z = set(lm)
+    T = promote_type(N, NM)
+
+    n = dim(Z)
+    w = ngens(MZ)
+    h = ngens(Z)
+
+    c = mapreduce(x -> x*center(Z), +, generators(MZ), init= center(MZ) * center(Z))
+
+    # generator 
+    G = Matrix{T}(undef, n, h * (w + 1))
+    G[:, 1:h] = center(MZ) * genmat(Z)
+    @inbounds for (i, A) in enumerate(generators(MZ))
+        G[:, h * i + 1 : h * (i + 1)] = A * genmat(Z)
+    end
+
+    return Zonotope(c, G)
+end
+
+"""
+    overapproximate(lm::LinearMap{N,S,NM,MAT},
+                         ::Type{<:Zonotope}) where {N,S<:AbstractZonotope{N},NM,
+                                                    MAT<:MatrixZonotopeProduct{NM}}
+
+Overapproximate the linear map of a zonotope through a product of matrix zonotopes,
+by recursively applying the overapproximation rule from the inside out.
+
+### Input
+
+- `lm` -- a linear map of a zonotope through a `MatrixZonotopeProduct`
+- `Zonotope` -- target type
+
+### Output
+
+An overapproximation of the linear map as a zonotope.
+"""
+function overapproximate(lm::LinearMap{N, S, NM, MAT},
+                         T::Type{<:Zonotope}) where {N, S<:AbstractZonotope{N}, NM, MAT<:MatrixZonotopeProduct{NM}}
+    MZs = factors(matrix(lm))
+    P = set(lm)
+    reduced = foldr((A, acc) -> overapproximate(A * acc, T), MZs; init = P)
+    return reduced
+end

src/LazyOperations/LinearMap.jl

@@ -4,14 +4,14 @@ export LinearMap,
        Projection
 
 """
-    LinearMap{N, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}}
-        <: AbstractAffineMap{N, S}
+    LinearMap{N,S<:LazySet{N},NM,
+              MAT<:Union{AbstractMatrix{NM},AbstractMatrixZonotope{NM}}} <: AbstractAffineMap{N,S}
 
 Type that represents a linear transformation ``M⋅X`` of a set ``X``.
 
 ### Fields
 
-- `M` -- matrix/linear map
+- `M` -- linear map; can be a concrete matrix (`AbstractMatrix`) or a set-valued matrix (`AbstractMatrixZonotope`)
 - `X` -- set
 
 ### Notes
@@ -96,13 +96,13 @@ EmptySet{Int64}(3)
 ```
 """
 struct LinearMap{N,S<:LazySet{N},NM,
-                 MAT<:AbstractMatrix{NM}} <: AbstractAffineMap{N,S}
+                 MAT<:Union{AbstractMatrix{NM},AbstractMatrixZonotope{NM}}} <: AbstractAffineMap{N,S}
     M::MAT
     X::S
 
     # default constructor with dimension check
     function LinearMap(M::MAT, X::S) where {N,S<:LazySet{N},NM,
-                                            MAT<:AbstractMatrix{NM}}
+                                            MAT<:Union{AbstractMatrix{NM},AbstractMatrixZonotope{NM}}}
         @assert dim(X) == size(M, 2) "a linear map of size $(size(M)) cannot " *
                                      "be applied to a set of dimension $(dim(X))"
         return new{N,S,NM,MAT}(M, X)
@@ -115,22 +115,22 @@ isconvextype(::Type{<:LinearMap{N,S}}) where {N,S} = isconvextype(S)
 
 """
 ```
-    *(M::Union{AbstractMatrix, UniformScaling, AbstractVector, Real},
+    *(M::Union{AbstractMatrix, UniformScaling, AbstractVector, Real, AbstractMatrixZonotope},
       X::LazySet)
 ```
 
 Alias to create a `LinearMap` object.
 
 ### Input
 
-- `M` -- linear map
+- `M` -- matrix or matrix zonotope
 - `X` -- set
 
 ### Output
 
 A lazy linear map, i.e., a `LinearMap` instance.
 """
-function *(M::Union{AbstractMatrix,UniformScaling,AbstractVector,Real},
+function *(M::Union{AbstractMatrix,UniformScaling,AbstractVector,Real,AbstractMatrixZonotope},
            X::LazySet)
     return LinearMap(M, X)
 end

src/Sets/SparsePolynomialZonotope/SparsePolynomialZonotopeModule.jl

@@ -4,15 +4,17 @@ using Reexport, Requires
 
 using ..LazySets: AbstractSparsePolynomialZonotope, AbstractReductionMethod,
                   genmat, GIR05, order, _remove_redundant_generators_polyzono,
-                  @validate
+                  MatrixZonotope, MatrixZonotopeProduct, ngens, generators,
+                  factors, @validate
 import IntervalArithmetic as IA
+using LinearAlgebra: I
 using Random: AbstractRNG, GLOBAL_RNG
 using ReachabilityBase.Arrays: remove_zero_columns
 using ReachabilityBase.Distribution: reseed!
 using ReachabilityBase.Require: require
 
 @reexport import ..API: center, isoperationtype, rand, scale, translate,
-                        translate!, exact_sum
+                        translate!, exact_sum, linear_map
 @reexport import ..LazySets: expmat, genmat_dep, genmat_indep, indexvector,
                              ngens_dep, ngens_indep, nparams, polynomial_order,
                              reduce_order, remove_redundant_generators
@@ -31,6 +33,7 @@ include("rand.jl")
 include("scale.jl")
 include("translate.jl")
 include("exact_sum.jl")
+include("linear_map.jl")
 
 include("expmat.jl")
 include("genmat_dep.jl")

src/Sets/SparsePolynomialZonotope/linear_map.jl

@@ -0,0 +1,82 @@
+"""
+	linear_map(MZ::MatrixZonotope, P::SparsePolynomialZonotope)
+
+Compute the linear map of a sparse polynomial zonotope through a matrix zonotope,
+following Proposition 1 of [HuangLBS2025](@citet).
+
+### Input
+
+- `lm` -- a linear map of a sparse polynomial zonotope through a matrix zonotope
+
+### Output
+
+A sparse polynomial zonotope representing the linear map ``MZ ⋅ P```.
+
+"""
+function linear_map(MZ::MatrixZonotope, P::SparsePolynomialZonotope)
+    if ngens_indep(P) > 0
+        error("an exact expression for the linear map is only available for " *
+              "`SparsePolynomialZonotope`s with no independent generators. " *
+              "Try using `overapproximate` instead")
+    end
+
+    @assert size(MZ, 2) == dim(P) "a linear map of size $(size(M)) cannot " *
+                                  "be applied to a set of dimension $(dim(X))"
+
+    T = promote_type(eltype(MZ), eltype(P))
+
+    n = dim(P)
+    w = ngens(MZ)
+    h = ngens_dep(P)
+
+    c = center(MZ) * center(P)
+
+    # compute matrix of dependent generators
+    G = Matrix{T}(undef, n, h + w + h * w)
+    G[:, 1:h] = center(MZ) * genmat_dep(P)
+    @inbounds for (i, A) in enumerate(generators(MZ))
+        G[:, h + i] = A * center(P)
+        G[:, (h + w + (i - 1) * h + 1):(h + w + i * h)] = A * genmat_dep(P)
+    end
+
+    Gi = Matrix{T}(undef, n, 0)
+
+    # compute exponent
+    Imat = Matrix{Int}(I, w, w)
+    Ê₁, Ê₂, idx = merge_id(indexvector(P), indexvector(MZ), expmat(P), Imat)
+    pₖ = size(Ê₁, 1)
+    E = Matrix{Int}(undef, pₖ, h + w + h * w)
+    E[:, 1:h] = Ê₁
+    E[:, (h + 1):(h + w)] = Ê₂
+    ones_h = ones(Int, 1, h)
+    @inbounds for l in 1:w
+        cstart = (h + w) + (l - 1) * h + 1
+        cend = (h + w) + l * h
+        col = @view Ê₂[:, l]
+        E[:, cstart:cend] = col * ones_h .+ Ê₁
+    end
+
+    return SparsePolynomialZonotope(c, G, Gi, E, idx)
+end
+
+"""
+	linear_map(MZP::MatrixZonotopeProduct, P::SparsePolynomialZonotope)
+
+Compute the linear map of a sparse polynomial zonotope through a matrix zonotope product
+by recursively applying the overapproximation rule from the inside out.
+
+### Input
+
+- `MZ` -- a matrix zonotope product
+- `P` -- a sparse polynomial zonotope 
+
+### Output
+
+A sparse polynomial zonotope representing the linear map ``MZ ⋅ P```.
+
+"""
+function linear_map(MZP::MatrixZonotopeProduct, P::SparsePolynomialZonotope)
+    MZs = factors(MZP)
+    reduced = foldr((A, acc) -> linear_map(A, acc), MZs; init=P)
+    return reduced
+end

test/Approximations/overapproximate.jl

@@ -154,6 +154,18 @@ for N in @tN([Float64, Float32, Rational{Int}])
     P = overapproximate(S, VPolytope)
     @test ispermutation([N[-1.5, -2.5], N[2.5, -0.5], N[3.5, 0.5], N[-0.5, 2.5], N[-1.5, 1.5]],
                         vertices_list(P))
+
+    # overapproximate the lm of a matrix zonotope and a zonotope
+    MZ = MatrixZonotope(N[1 1; -1 1], [N[1 0; 1 2]])
+    Z = Zonotope(N[2, 0], N[-1 2; -1 0])
+    res = overapproximate(MZ * Z, Zonotope)
+    @test center(res) == N[4, 0]
+    @test genmat(res) == hcat(N[-2 2; 0 -2], N[-1 2; -3 2])
+
+    MZ2= MatrixZonotope(N[1 0; 0 1], [N[2 0; 1 -1]])
+    MZP = MZ2 * MZ
+    res2 = overapproximate(MZP * Z, Zonotope)
+    @test res2 == overapproximate(MZ2 * res, Zonotope)
 end
 
 # tests that do not work with Rational{Int}

test/Sets/SparsePolynomialZonotope.jl

@@ -206,6 +206,30 @@ for N in @tN([Float64, Float32, Rational{Int}])
                              0 0 0 0 0 1 2 1 2
                              1 0 0 1 1 0 0 1 1]
     end
+
+    # linear map with matrix zonotope
+    MZ = MatrixZonotope(N[1 1; -1 1], [N[1 0; 1 2]])
+    P = SparsePolynomialZonotope(N[1, -1], N[1 1; 0 -1], Matrix{N}(undef, 2, 0), [2 1; 0 1; 1 0], [1, 2, 3])
+    res = linear_map(MZ, P)
+    @test center(res) == [0, -2]
+    @test genmat_dep(res) == hcat(N[1 0; -1 -2], N[1, -1], N[1 1; 1 -1])
+    @test genmat_indep(res) == Matrix{N}(undef, 2, 0)
+    @test expmat(res) == hcat([2 1; 0 1; 1 0], [1; 0; 0], [3 2; 0 1; 1 0])
+
+    # case: 0 gens in MZ
+    MZ = MatrixZonotope(N[1 1; -1 1], Vector{Matrix{N}}())
+    res = linear_map(MZ, P)
+    @test res == linear_map(N[1 1; -1 1], P)
+
+    #case: error 
+    P_err = SparsePolynomialZonotope(N[1, -1], N[1 1; 0 -1], hcat(N[0, 1]), [2 1; 0 1; 1 0], [1, 2, 3])
+    @test_throws ErrorException linear_map(MZ, P_err)
+
+    #MZP linear map
+    MZ2 = MatrixZonotope(N[1.1 0.9; -1.1 1.1], [N[1.1 -0.1; 0.9 2.1]])
+    res = linear_map(MZ * MZ2, P)  
+    P_in = linear_map(MZ2, P)
+    @test res == linear_map(MZ, P_in)
 end
 
 for N in @tN([Float64, Float32])