Lazy Intersections

In this section we illustrate the use of lazy intersection in LazySets. We will use the ellipsoid set type.

• Lazy Intersections

An ellipsoid $E$ can be created by giving its center $c$ and its shape matrix $Q$, which should be positive definite, i.e. its eigenvalues must be positive. Mathematically, it is the set

Let's make two rotated ellipsoids and plot them in the same pair of axes.

using Plots, LazySets, LazySets.Approximations

E₁ = Ellipsoid(zeros(2), [1 0; 0 2.])
E₂ = Ellipsoid(ones(2), [2 0; 0 1.])

pell = plot(E₁, 1e-3, aspectratio=1, alpha=.5)
pell = plot!(pell, E₂, 1e-3, alpha=.5)

Now let's take the lazy intersection of the ellipses:

Z = E₁ ∩ E₂
typeof(Z)

Intersection{Float64,Ellipsoid{Float64},Ellipsoid{Float64}}

On the other hand, the concrete intersection of sets, called intersection in LazySets, is not yet available for ellipsoids:

hasmethod(intersection, Tuple{typeof(E₁), typeof(E₂)})

false

So how can we work with the intersection of the ellipsoids?

One way is to overapproximate them by polygons (or polytopes in higher dims) and then take their intersection, because this function is defined, whose return type is again a HPolytope:

hasmethod(intersection, Tuple{HPolytope{Float64}, HPolytope{Float64}})

true

import LazySets.Approximations.overapproximate

# the parameter epsilon controls the accuracy of the iterative refinement,
# with respect to the Hausdorff distance
H₁(ε) = overapproximate(E₁, HPolygon, ε)
H₂(ε) = overapproximate(E₂, HPolygon, ε)

# using the concrete hpolytope-hpolytope intersection here
Hint(ε) = intersection(convert.(HPolytope, [H₁(ε), H₂(ε)])...);

Hint (generic function with 1 method)

pell = plot(E₁, 1e-3, aspectratio=1, alpha=.5)
pell = plot!(pell, E₂, 1e-3, alpha=.5)
pεsmaller = plot!(pell, convert(HPolygon, Hint(0.5)), alpha=.4)

pell = plot(E₁, 1e-3, aspectratio=1, alpha=.5)
pell = plot!(pell, E₂, 1e-3, alpha=.5)
pεbigger = plot!(pell, convert(HPolygon, Hint(0.05)), alpha=.4)

plot(pεsmaller, pεbigger, layout=(1, 2))

Note how dividing the $\varepsilon$ threshold by 10 makes the polygonal overapproximation of the intersection tighter.

Yet another approach is to directly query the directions of the lazy intersection E₁ ∩ E₂. We can overapproximate using template directions, such as a box, an octagon, or other.

The idea behind the template overapproximation method is to use the property that the support function of the intersection of two convex sets is upper bounded by the min of the support function of each set. We can see in the following experiments that the resulting set is quite tight.

import LazySets.Approximations.overapproximate
using LazySets.Approximations, Polyhedra

# overapproximate the lazy intersection using a box
Xbox = overapproximate(E₁ ∩ E₂, BoxDirections(2))

# overapproximate the lazy intersection using octagonal directions
Xoct = overapproximate(E₁ ∩ E₂, OctDirections(2))

pell = plot(E₁, 1e-3, aspectratio=1, alpha=.5)
pell = plot!(pell, E₂, 1e-3, alpha=.5)
pbox = plot!(pell, Xbox, alpha=.4)

pell = plot(E₁, 1e-3, aspectratio=1, alpha=.5)
pell = plot!(pell, E₂, 1e-3, alpha=.5)
poct = plot!(pell, Xoct, alpha=.4)

plot(pbox, poct, layout=(1, 2))

Using support function evaluations over a set of fixed directions is in general more efficient than iterative refinement, but the drawback is that one does not have control on the overapproximation error. Moreover, iterative refinement is currently only available in two dimensions, but overapproximation with template directions can be used in any dimension.

Let's time it!

using BenchmarkTools

@btime overapproximate($E₁ ∩ $E₂, BoxDirections(2))
@btime overapproximate($E₁ ∩ $E₂, OctDirections(2));

HPolytope{Float64}(HalfSpace{Float64,VN} where VN<:AbstractArray{Float64,1}[HalfSpace{Float64,SparseVector{Float64,Int64}}(  [1]  =  1.0
  [2]  =  1.0, 1.73205), HalfSpace{Float64,SparseVector{Float64,Int64}}(  [1]  =  1.0
  [2]  =  -1.0, 1.73205), HalfSpace{Float64,SparseVector{Float64,Int64}}(  [1]  =  -1.0
  [2]  =  1.0, 1.73205), HalfSpace{Float64,SparseVector{Float64,Int64}}(  [1]  =  -1.0
  [2]  =  -1.0, -0.267949), HalfSpace{Float64,SingleEntryVector{Float64}}([1.0, 0.0], 1.0), HalfSpace{Float64,SingleEntryVector{Float64}}([0.0, 1.0], 1.41421), HalfSpace{Float64,SingleEntryVector{Float64}}([0.0, -1.0], 0.0), HalfSpace{Float64,SingleEntryVector{Float64}}([-1.0, 0.0], 0.414214)])

We can work with higher dimensional ellipsoids as well:

using LinearAlgebra

# a random ellipsoid in n-dimensions
function rand_ellipsoid(n)
    A = rand(n,n)
    Q = (A+transpose(A))/2 + n * I
    Ellipsoid(rand(n), Q)
end;

rand_ellipsoid (generic function with 1 method)

for n in [2, 5, 50, 100]
    println("\nn = $n\n")
    global E₁, E₂ = rand_ellipsoid(n), rand_ellipsoid(n)

    # overapproximate the lazy intersection using an n-dimensional box
    @btime overapproximate($E₁ ∩ $E₂, BoxDirections($n))

    # overapproximate the lazy intersection using octagonal directions in R^n
    @btime overapproximate($E₁ ∩ $E₂, OctDirections($n))
end;

n = 2

  1.754 μs (61 allocations: 1.44 KiB)
  5.226 μs (134 allocations: 4.31 KiB)

n = 5

  4.152 μs (145 allocations: 3.36 KiB)
  36.561 μs (802 allocations: 33.33 KiB)

n = 50

  56.066 μs (1405 allocations: 32.27 KiB)
  7.935 ms (79918 allocations: 11.28 MiB)

n = 100

  147.663 μs (2805 allocations: 64.23 KiB)
  99.805 ms (319820 allocations: 72.15 MiB)