Add overapproximate for matrix zonotope exponential (#4000)

* add conversion from `IntervalMatrix` to `MatrixZonotope`

* add mz exp overapproxination

* Apply suggestions from code review - part 1

Co-authored-by: Christian Schilling <[email protected]>

* apply suggestions from code review -part 2

* fix tests and docstring

* Update test/Approximations/overapproximate_expmap.jl

Co-authored-by: Christian Schilling <[email protected]>

* Update docs/src/lib/conversion.md

Co-authored-by: Christian Schilling <[email protected]>

---------

Co-authored-by: Christian Schilling <[email protected]>

Author: alecarraro

docs/src/lib/conversion.md

@@ -43,3 +43,11 @@ convert(::Type{SimpleSparsePolynomialZonotope}, ::AbstractSparsePolynomialZonoto
 convert(::Type{SparsePolynomialZonotope}, ::AbstractZonotope{N}) where {N}
 convert(::Type{SparsePolynomialZonotope}, ::SimpleSparsePolynomialZonotope{N}) where {N}
 ```
+
+```@meta
+CurrentModule = LazySets.MatrixZonotopeModule
+```
+
+```@docs
+convert(::Type{MatrixZonotope}, ::IntervalMatrices.IntervalMatrix)
+```

src/Approximations/overapproximate_expmap.jl

@@ -156,6 +156,59 @@ function load_intervalmatrices_overapproximation_expmap()
             res = overapproximate(E * Z, Zonotope)
             return res
         end
+
+        """
+            overapproximate(expA::MatrixZonotopeExp{N,T}, ::Type{<:MatrixZonotope},
+                                 k::Int=2) where {N,T<:AbstractMatrixZonotope{N}}
+        
+        Overapproximate the matrix zonotope exponential ``exp(\\mathcal{A})``
+
+        ### Input
+
+        - `expA` -- `MatrixZonotopeExp`
+        - `MatrixZonotope` -- target type
+        - `k` -- (default: `2`) the order of the Taylor expansion
+
+        ### Output 
+        
+        A matrix zonotope overapproximating the matrix zonotope exponential
+
+        ### Algorithm
+        
+        The expansion
+
+        ```math 
+        exp(\\mathcal{A}) ⊆ \\sum_i^k \\frac{\\mathcal{A}^i}{i!} + E_k
+        ```
+
+        is computed by overapproximating the matrix zonotope powers ``A^i`` 
+        for ``i=0, …, k``. 
+        The remainder term ``E_k`` is computed through interval arithmetic 
+        following Proposition 4.1 by [AlthoffKS11](@citet). 
+        """
+        function overapproximate(expA::MatrixZonotopeExp{N,T}, ::Type{<:MatrixZonotope},
+                                 k::Int=2) where {N,T<:AbstractMatrixZonotope{N}}
+            # overapproximate the exponent, which can be a product A*B*…
+            MZP = MatrixZonotopeProduct(expA.M)
+            MZ = overapproximate(MZP, MatrixZonotope)
+
+            # compute the Taylor expansion
+            powers = Vector{typeof(MZ)}(undef, k)
+            powers[1] = MZ
+            @inbounds for i in 2:k
+                term = overapproximate(MZ * powers[i - 1], MatrixZonotope)
+                powers[i] = scale(1 / i, term)
+            end
+            W = reduce(minkowski_sum, powers)
+            W = MatrixZonotope(center(W) + Matrix{N}(I, size(W)), generators(W))
+
+            # overapproximate mat zon by interval matrix and overapproximate remainder
+            IM = overapproximate(MZ, IntervalMatrix)
+            E = IntervalMatrices._exp_remainder(IM, N(1), k)
+            
+            res = minkowski_sum(W, convert(MatrixZonotope, E))
+            return remove_redundant_generators(res)
+        end
     end
 end
 

src/MatrixSets/MatrixZonotopeModule.jl

@@ -1,6 +1,6 @@
 module MatrixZonotopeModule
 
-using Reexport
+using Reexport, Requires
 
 using Random: AbstractRNG, GLOBAL_RNG
 using ReachabilityBase.Distribution: reseed!
@@ -26,10 +26,13 @@ include("minkowski_sum.jl")
 include("reduce_order.jl")
 include("remove_redundant_generators.jl")
 include("reshape.jl")
+include("convert.jl")
 
 include("generators.jl")
 include("indexvector.jl")
 include("ngens.jl")
 include("order.jl")
 
+include("init.jl")
+
 end  # module

src/MatrixSets/convert.jl

@@ -0,0 +1,40 @@
+function load_intervalmatrices_conversion()
+    return quote
+        using .IntervalMatrices: IntervalMatrix, mid, sup, radius
+
+        """
+            convert(::Type{MatrixZonotope}, IM::IntervalMatrix)
+        
+        Convert an interval matrix to a matrix zonotope
+
+        ### Input
+
+        - `MatrixZonotope` -- target type
+        - `IM` -- an interval matrix
+
+        ### Output
+
+        A matrix zonotope with one generator
+
+        ### Example 
+
+        ```jldoctest
+        julia> using LazySets, IntervalMatrices
+
+        julia> IM = IntervalMatrix([interval(-1.1, -0.9) interval(-4.1, -3.9);
+                    interval(3.9, 4.1) interval(-1.1, -0.9)])
+        2×2 IntervalMatrix{Float64, IntervalArithmetic.Interval{Float64}, Matrix{IntervalArithmetic.Interval{Float64}}}:
+         [-1.10001, -0.9]  [-4.1, -3.89999]
+          [3.89999, 4.1]   [-1.10001, -0.9]
+        
+        julia> MZ = convert(MatrixZonotope, IM)
+        MatrixZonotope{Float64, Matrix{Float64}}([-1.0 -4.0; 4.0 -1.0], [[0.10000000000000009 0.10000000000000009; 0.10000000000000009 0.10000000000000009]], [1])
+        ```
+        """
+        function Base.convert(::Type{MatrixZonotope}, IM::IntervalMatrix)
+            c = mid(IM)
+            G = [radius(IM)]
+            return MatrixZonotope(c, G)
+        end
+    end
+end

src/MatrixSets/init.jl

@@ -0,0 +1,3 @@
+function __init__()
+    @require IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9" eval(load_intervalmatrices_conversion())
+end

test/Approximations/overapproximate_expmap.jl

@@ -1,4 +1,5 @@
 using LazySets, Test, SparseArrays
+using LazySets.MatrixZonotopeModule: vectorize
 if !isdefined(@__MODULE__, Symbol("@tN"))
     macro tN(v)
         return v
@@ -51,6 +52,22 @@ for N in @tN([Float32, Float64])
         pts = sample(Pex, 10; sampler=sampler)
         @test all(p ∈ Zex for p in pts)
 
+        # overapproximate MatrixZonotopeExp
+        
+        # test degenerate case
+        c = N[1 -2; 2 -1]
+        A = MatrixZonotope(c, [zeros(N, 2, 2)])
+        expA = MatrixZonotopeExp(A)
+        res = overapproximate(expA, MatrixZonotope, 20) #for large k it converges to exp(c)
+        @test isapprox(center(res), exp(c))
+        @test ngens(res)== 0
+
+        # degenerate case + product 
+        B = MatrixZonotope(N[1 0; 0 1], Matrix{N}[])
+        expAB = MatrixZonotopeExp(A * B)
+        res2 = overapproximate(expA, MatrixZonotope, 20)
+        @test isapprox(center(res), center(res2))
+
         @static if isdefined(@__MODULE__, :ExponentialUtilities) || isdefined(@__MODULE__, :Expokit)
             # matrix
             mzexp = SparseMatrixExp(sparse(M))
@@ -61,3 +78,16 @@ for N in @tN([Float32, Float64])
         end
     end
 end
+
+for N in @tN([Float64])
+    # This test is only run with Float64 since it is expensive
+    # and requires high Taylor order `k` to pass on Float32
+    @static if isdefined(@__MODULE__, :IntervalMatrices)
+        #inclusion 
+        C = MatrixZonotope(N[1 -2; 2 -1], [N[0.1 0.05; 0 0.1]])
+        expC = MatrixZonotopeExp(C)
+        res = reduce_order(overapproximate(expC, MatrixZonotope, 4), 1)
+        res2 = reduce_order(overapproximate(expC, MatrixZonotope, 8), 1)
+        @test vectorize(res2) ⊆ vectorize(res)       
+    end
+end

test/MatrixSets/MatrixZonotope.jl

@@ -1,4 +1,6 @@
 using LazySets, Test
+IA = LazySets.IA
+using LazySets.IA: interval
 using LazySets.MatrixZonotopeModule: vectorize, matrixize
 if !isdefined(@__MODULE__, Symbol("@tN"))
     macro tN(v)
@@ -122,6 +124,22 @@ for N in @tN([Float64, Float32, Rational{Int}])
     @test size(expMZ) == size(A_mz)
     @test_throws AssertionError MatrixZonotopeExp(B_mz)
 
+    # convert IM to MZ
+    @static if isdefined(@__MODULE__, :IntervalMatrices)
+        using IntervalMatrices: IntervalMatrix
+        IM = IntervalMatrix([interval(-N(1.1), -N(0.9)) interval(-N(4.1), -N(3.9));
+                            interval(N(3.9), N(4.1)) interval(-N(1.1), -N(0.9))])
+        MZ = convert(MatrixZonotope, IM)
+        if N==Rational{Int} # isapprox is sensitive to minor rounding using Rational{Int}
+            T = Float64
+            @test isapprox(convert(Matrix{T}, center(MZ)), T[-1.0 -4.0; 4.0 -1.0])
+            @test isapprox(convert(Vector{Matrix{T}}, generators(MZ)), [T[0.1 0.1; 0.1 0.1]])
+        else
+            @test isapprox(center(MZ), N[-1.0 -4.0; 4.0 -1.0])
+            @test isapprox(generators(MZ), [N[0.1 0.1; 0.1 0.1]])
+        end
+    end
+
     # vectorize and matrixize
     c = N[1 0; 0 3]
     gens = [N[1 -1; 0 2]]