Improved order reduction methods (#2220)

* improved order reduction methods * fix doctest * add missing method * Update src/Approximations/overapproximate.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Approximations/overapproximate.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Interfaces/AbstractZonotope.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Interfaces/AbstractZonotope.jl Co-authored-by: Christian Schilling <[email protected]> * Update src/Sets/Zonotope.jl Co-authored-by: Christian Schilling <[email protected]> * review comments * fix test * fix tests Co-authored-by: Christian Schilling <[email protected]>
JuliaReach · Jul 10, 2022 · 3382153 · 3382153
1 parent ce69f5c
commit 3382153
Show file tree

Hide file tree

Showing 7 changed files with 185 additions and 51 deletions.
diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md
@@ -356,6 +356,8 @@ vertices_list(::AbstractZonotope)
 order(::AbstractZonotope)
 togrep(::AbstractZonotope)
 remove_redundant_generators(::AbstractZonotope)
+reduce_order(::AbstractZonotope, ::Number, ::AbstractReductionMethod)
+LazySets.AbstractReductionMethod
 ```
 
 ### Implementations

diff --git a/docs/src/lib/sets/Zonotope.md b/docs/src/lib/sets/Zonotope.md
@@ -14,7 +14,6 @@ scale(::Real, ::Zonotope)
 scale!(::Real, Z::Zonotope)
 ngens(::Zonotope)
 togrep(::Zonotope)
-reduce_order(::Zonotope, ::Union{Integer, Rational})
 low(::Zonotope, ::Int)
 high(::Zonotope, ::Int)
 remove_zero_generators(::Zonotope)
@@ -52,3 +51,11 @@ Inherited from [`AbstractZonotope`](@ref):
 * [`constraints_list`](@ref constraints_list(::AbstractZonotope{N}; ::Bool=true) where {N<:AbstractFloat})
 * [`vertices_list`](@ref vertices_list(::AbstractZonotope))
 * [`order`](@ref order(::AbstractZonotope))
+* [`reduce_order`](@ref reduce_order(::AbstractZonotope, ::Number, ::AbstractReductionMethod))
+
+## Order reduction methods
+
+```@docs
+LazySets.COMB03
+LazySets.GIR05
+```
diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl
@@ -913,7 +913,7 @@ Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([1.0, -2.1], [2.5, 6.5
 julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom)
 [-1.5, 3.5] × [-8.60001, 4.40001]
 
-julia> H = convert(Hyperrectangle, Y) # this IntevalBox is the same as X
+julia> H = convert(Hyperrectangle, Y) # this IntervalBox is the same as X
 Hyperrectangle{Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}([1.0, -2.1000000000000005], [2.5, 6.500000000000001])
 ```
 However, the zonotope returns better results if we want to approximate the `TM`,
@@ -1400,53 +1400,24 @@ end
 """
     overapproximate(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N}
 
-Reduce the order of a zonotope by overapproximating with a zonotope with less
-generators.
+Reduce the order of a zonotope by overapproximating with a zonotope with fewer generators.
 
 ### Input
 
-- `Z` -- zonotope
+- `Z`        -- zonotope
 - `Zonotope` -- desired type for dispatch
-- `r` -- desired order
+- `r`        -- desired order
 
 ### Output
 
 A new zonotope with less generators, if possible.
 
 ### Algorithm
 
-This function implements the algorithm described in A. Girard's
-*Reachability of Uncertain Linear Systems Using Zonotopes*, HSCC. Vol. 5. 2005.
-
-If the desired order is smaller than one, the zonotope is *not* reduced.
+This function falls back to `reduce_order` with the default algorithm.
 """
 function overapproximate(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N}
-    c, G = Z.center, Z.generators
-    d, p = dim(Z), ngens(Z)
-
-    if r * d >= p || r < 1
-        # do not reduce
-        return Z
-    end
-
-    h = zeros(N, p)
-    for i in 1:p
-        h[i] = norm(G[:, i], 1) - norm(G[:, i], Inf)
-    end
-    ind = sortperm(h)
-
-    m = p - floor(Int, d * (r - 1)) # subset of ngens that are reduced
-    rg = G[:, ind[1:m]] # reduced generators
-
-    # interval hull computation of reduced generators
-    Gbox = Diagonal(sum(abs, rg, dims=2)[:])
-    if m < p
-        Gnotred = G[:, ind[m+1:end]]
-        Gred = [Gnotred Gbox]
-    else
-        Gred = Gbox
-    end
-    return Zonotope(c, Gred)
+    reduce_order(Z, r)
 end
 
 """

diff --git a/src/Initialization/init_StaticArrays.jl b/src/Initialization/init_StaticArrays.jl
@@ -4,3 +4,4 @@ eval(load_split_static())
 eval(load_genmat_2D_static())
 eval(load_genmat_ballinf_static())
 eval(load_genmat_hyperrectangle_static())
+eval(load_reduce_order_static())
diff --git a/src/Interfaces/AbstractZonotope.jl b/src/Interfaces/AbstractZonotope.jl
@@ -8,7 +8,8 @@ export AbstractZonotope,
        translate,
        translate!,
        togrep,
-       split!
+       split!,
+       reduce_order
 
 """
     AbstractZonotope{N} <: AbstractCentrallySymmetricPolytope{N}
@@ -731,3 +732,47 @@ By default this function returns the input zonotope. Subtypes of
 function remove_redundant_generators(Z::AbstractZonotope)
     return Z  # fallback implementation
 end
+
+"""
+    AbstractReductionMethod
+
+Abstract supertype for zonotope order reduction methods.
+"""
+abstract type AbstractReductionMethod end
+
+"""
+    reduce_order(Z::AbstractZonotope, r::Number, method::AbstractReductionMethod=GIR05())
+
+Reduce the order of a zonotope by overapproximating with a zonotope with fewer generators.
+
+### Input
+
+- `Z`      -- zonotope
+- `r`      -- desired order
+- `method` -- (optional, default: `GIR05()`) the reduction method used
+
+### Output
+
+A new zonotope with fewer generators, if possible.
+
+### Algorithm
+
+The available algorithms are:
+
+```jldoctest; setup = :(using LazySets: subtypes, AbstractReductionMethod)
+julia> subtypes(AbstractReductionMethod)
+2-element Vector{Any}:
+ LazySets.COMB03
+ LazySets.GIR05
+```
+
+See the documentation of each algorithm for references. These methods split the given zonotope `Z` into two zonotopes,
+`K` and `L`, where `K` contains the most "representative" generators and `L` contains the generators that are reduced, `Lred`,
+using a box overapproximation. We follow the notation from [YS18]. See also [KSA17].
+
+- [KSA17] Kopetzki, A. K., Schürmann, B., & Althoff, M. (2017). *Methods for order reduction of zonotopes.* In 2017 IEEE 56th Annual CDC (pp. 5626-5633). IEEE.
+- [YS18] Yang, X., & Scott, J. K. (2018). *A comparison of zonotope order reduction techniques.* Automatica, 95, 378-384.
+"""
+function reduce_order(Z::AbstractZonotope, r::Number, method::AbstractReductionMethod=GIR05())
+    reduce_order(convert(Zonotope, Z), r, method)
+end
diff --git a/src/Sets/Zonotope.jl b/src/Sets/Zonotope.jl
@@ -329,30 +329,135 @@ function scale!(α::Real, Z::Zonotope)
     return Z
 end
 
+# ==================================
+# Zonotope order reduction methods
+# ==================================
+
 """
-    reduce_order(Z::Zonotope, r::Union{Integer, Rational})
+    GIR05 <: AbstractReductionMethod
 
-Reduce the order of a zonotope by overapproximating with a zonotope with fewer
-generators.
+Zonotope order reduction method from [GIR05].
 
-### Input
+- [G05] A. Girard. *Reachability of Uncertain Linear Systems Using Zonotopes*, HSCC. Vol. 5. 2005.
+"""
+struct GIR05 <: AbstractReductionMethod end
 
-- `Z` -- zonotope
-- `r` -- desired order
+"""
+    COMB03 <: AbstractReductionMethod
 
-### Output
+Zonotope order reduction method from [COMB03].
 
-A new zonotope with fewer generators, if possible.
+- [COMB03] C. Combastel. *A state bounding observer based on zonotopes.* In Proc. of the European Control Conference, p. 2589–2594, 2003.
+"""
+struct COMB03 <: AbstractReductionMethod end
 
-### Algorithm
+function reduce_order(Z::Zonotope, r::Number, method::AbstractReductionMethod=GIR05())
+    r >= 1 || throw(ArgumentError("the target order should be at least 1, but it is $r"))
+    c = Z.center
+    G = Z.generators
+    n, p = size(G)
 
-See `overapproximate(Z::Zonotope, ::Type{<:Zonotope}, r::Union{Integer, Rational})`
-for details.
-"""
-function reduce_order(Z::Zonotope, r::Union{Integer, Rational})
-    return overapproximate(Z, Zonotope, r)
+    # r is bigger than the order of Z => don't reduce
+    (r * n >= p) && return Z
+
+    if isone(r)
+        # if r = 1 => m = 0 and the generators need not be sorted
+        Lred = _interval_hull(G, 1:p) 
+        return Zonotope(c, Lred)
+    end
+
+    # sort generators
+    indices = Vector{Int}(undef, p)
+    _weighted_gens!(indices, G, method)
+
+    # the first m generators have greatest weight
+    m = floor(Int, n * (r - 1))
+
+    # compute interval hull of L
+    Lred = _interval_hull(G, view(indices, (m+1):p))
+
+    # concatenate non-reduced and reduced generators
+    Gred = _hcat_KLred(G, view(indices, 1:m), Lred)
+
+    return Zonotope(c, Gred)
+end
+
+# Return the indices of the generators in G (= columns) sorted according to decreasing 2-norm.
+# The generator index with highest score goes first.
+function _weighted_gens!(indices, G::AbstractMatrix{N}, ::COMB03) where {N}
+    p = size(G, 2)
+    weights = Vector{N}(undef, p)
+    @inbounds for j in 1:p
+        v = view(G, :, j)
+        weights[j] = norm(v, 2)
+    end
+    sortperm!(indices, weights, rev=true, initialized=false)
+    return indices
 end
 
+# Return the indices of the generators in G (= columns) sorted according to ||⋅||₁ - ||⋅||∞ difference.
+# The generator index with highest score goes first.
+function _weighted_gens!(indices, G::AbstractMatrix{N}, ::GIR05) where {N}
+    n, p = size(G)
+    weights = Vector{N}(undef, p)
+    @inbounds for j in 1:p
+        aux_norm_1 = zero(N)
+        aux_norm_inf = zero(N)
+        for i in 1:n
+            abs_Gij = abs(G[i, j])
+            aux_norm_1 += abs_Gij
+            if aux_norm_inf < abs_Gij
+                aux_norm_inf = abs_Gij
+            end
+        end
+        weights[j] = aux_norm_1 - aux_norm_inf
+    end
+    sortperm!(indices, weights, rev=true, initialized=false)
+    return indices
+end
+
+# compute interval hull of the generators of G (= columns) corresponding to `indices`
+function _interval_hull(G::AbstractMatrix{N}, indices) where {N}
+    n, p = size(G)
+    Lred = zeros(N, n, n)
+    @inbounds for i in 1:n
+        for j in indices
+            Lred[i, i] += abs(G[i, j])
+        end
+    end
+    return Lred
+end
+
+# given an n x p matrix G and a vector of m integer indices with m <= p,
+# concatenate the columns of G given by `indices` with the matrix Lred
+function _hcat_KLred(G::AbstractMatrix, indices, Lred::AbstractMatrix)
+    K = view(G, :, indices)
+    return hcat(K, Lred)
+end
+
+function load_reduce_order_static()
+return quote
+
+# implementation for static arrays
+function _interval_hull(G::SMatrix{n, p, N, L}, indices) where {n, p, N, L}
+    Lred = zeros(MMatrix{n, n, N})
+    @inbounds for i in 1:n
+        for j in indices
+            Lred[i, i] += abs(G[i, j])
+        end
+    end
+    return SMatrix{n, n}(Lred)
+end
+
+# implementation for static arrays
+function _hcat_KLred(G::SMatrix{n, p, N, L1}, indices, Lred::SMatrix{n, n, N, L2}) where {n, p, N, L1, L2}
+    m = length(indices)
+    K = SMatrix{n, m}(view(G, :, indices))
+    return hcat(K, Lred)
+end
+
+end end # quote / load_reduce_order_static
+
 # ============================
 # Zonotope splitting methods
 # ============================

diff --git a/test/Sets/Zonotope.jl b/test/Sets/Zonotope.jl
@@ -122,6 +122,7 @@ for N in [Float64, Rational{Int}, Float32]
     @test ngens(Z) == 4
     @test genmat(Z) == Z.generators
     @test ispermutation(collect(generators(Z)), [genmat(Z)[:, j] for j in 1:ngens(Z)])
+
     # test order reduction
     Zred1 = reduce_order(Z, 1)
     @test ngens(Zred1) == 2
@@ -137,6 +138,8 @@ for N in [Float64, Rational{Int}, Float32]
     @test ngens(Znogen) == 0
     @test genmat(Znogen) == Matrix{N}(undef, 2, 0)
     @test collect(generators(Znogen)) == Vector{N}()
+    Zs = Zonotope(SVector{2}(Z.center), SMatrix{2, 6}(Z.generators))
+    @test reduce_order(Zs, 2) isa Zonotope{N, SVector{2, N}, SMatrix{2, 4, N, 8}}
 
     # conversion from zonotopic sets
     Z = Zonotope(N[0, 0], hcat(N[1, 1]))