diff --git a/docs/src/lib/operations.md b/docs/src/lib/operations.md
index 6933a73b18..94fc9cf77e 100644
--- a/docs/src/lib/operations.md
+++ b/docs/src/lib/operations.md
@@ -101,10 +101,16 @@ monotone_chain!
 Intersection
 ∩(::LazySet, ::LazySet)
 dim(::Intersection)
+ρ(::AbstractVector{N}, ::Intersection{N}) where {N<:Real}
+ρ(::AbstractVector{N}, ::Intersection{N, <:LazySet{N}, <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}}) where {N<:Real}
+ρ(::AbstractVector{N}, ::Intersection{N, <:LazySet{N}, <:AbstractPolytope{N}}) where {N<:Real}
+isempty(::Intersection)
 σ(::AbstractVector{Real}, ::Intersection{Real})
 ∈(::AbstractVector{Real}, ::Intersection{Real})
-isempty(::Intersection)
-ρ(::AbstractVector{N}, ::Intersection{N, <:LazySet, S}) where {N<:AbstractFloat, S<:Union{HalfSpace, Hyperplane, Line}}
+isempty_known(::Intersection)
+set_isempty!(::Intersection, ::Bool)
+_line_search
+_projection
 ```
 
 Inherited from [`LazySet`](@ref):
@@ -113,6 +119,12 @@ Inherited from [`LazySet`](@ref):
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 * [`an_element`](@ref an_element(::LazySet{Real}))
 
+#### Intersection cache
+
+```@docs
+IntersectionCache
+```
+
 ### ``n``-ary Intersection
 
 ```@docs
diff --git a/src/Intersection.jl b/src/Intersection.jl
index 3397f9821f..4b3cf2c192 100644
--- a/src/Intersection.jl
+++ b/src/Intersection.jl
@@ -1,9 +1,46 @@
 import Base: isempty, ∈, ∩
 
 export Intersection,
+       isempty_known,
+       set_isempty!,
        IntersectionArray,
        array
 
+"""
+    IntersectionCache
+
+Container for information cached by a lazy `Intersection` object.
+
+### Fields
+
+- `isempty` -- is the intersection empty? There are three possible states,
+               encoded as `Int8` values -1, 0, 1:
+
+    * ``-1`` - it is currently unknown whether the intersection is empty or not
+    *  ``0`` - intersection is not empty
+    *  ``1`` - intersection is empty
+"""
+mutable struct IntersectionCache
+    isempty::Int8
+
+    # default constructor
+    IntersectionCache() = new(Int8(-1))
+end
+
+function isempty_known(c::IntersectionCache)::Bool
+    return c.isempty != Int8(-1)
+end
+
+function isempty(c::IntersectionCache)::Bool
+    @assert isempty_known(c) "'isempty_known' only works if 'isempty' " *
+        "returns 'true'"
+    return c.isempty == Int8(1)
+end
+
+function set_isempty!(c::IntersectionCache, isempty::Bool)
+    c.isempty = isempty ? Int8(1) : Int8(0)
+end
+
 """
     Intersection{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
 
@@ -11,13 +48,15 @@ Type that represents the intersection of two convex sets.
 
 ### Fields
 
-- `X` -- convex set
-- `Y` -- convex set
+- `X`     -- convex set
+- `Y`     -- convex set
+- `cache` -- internal cache for avoiding recomputation; see
+             [`IntersectionCache`](@ref)
 
 ### Examples
 
-Create an expression, ``Z``, that lazily represents the intersection of two squares
-``X`` and ``Y``:
+Create an expression, ``Z``, which lazily represents the intersection of two
+squares ``X`` and ``Y``:
 
 ```jldoctest lazy_intersection
 julia> X, Y = BallInf([0,0.], 0.5), BallInf([1,0.], 0.65);
@@ -38,8 +77,8 @@ julia> isempty(Z)
 false
 ```
 
-Do not confuse `Intersection` with the concrete operation, that is computed with
-the lowercase `intersection`:
+Do not confuse `Intersection` with the concrete operation, which is computed
+with the lowercase `intersection` function:
 
 ```jldoctest lazy_intersection
 julia> W = intersection(X, Y)
@@ -49,13 +88,14 @@ Hyperrectangle{Float64}([0.425, 0.0], [0.075, 0.5])
 struct Intersection{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
     X::S1
     Y::S2
+    cache::IntersectionCache
 
     # default constructor with dimension check
     function Intersection{N, S1, S2}(X::S1, Y::S2) where
             {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
         @assert dim(X) == dim(Y) "sets in an intersection must have the same " *
             "dimension"
-        return new{N, S1, S2}(X, Y)
+        return new{N, S1, S2}(X, Y, IntersectionCache())
     end
 end
 
@@ -74,6 +114,42 @@ Alias for `Intersection`.
 ∩(X::LazySet, Y::LazySet) = Intersection(X, Y)
 
 
+# --- cache propagation functions ---
+
+
+"""
+    isempty_known(cap::Intersection)
+
+Ask whether the status of emptiness is known.
+
+### Input
+
+- `cap` -- intersection of two convex sets
+
+### Output
+
+`true` iff the emptiness status is known.
+In this case, `isempty(cap)` can be used to obtain the status.
+"""
+function isempty_known(cap::Intersection)
+    return isempty_known(cap.cache)
+end
+
+"""
+    set_isempty!(cap::Intersection, isempty::Bool)
+
+Set the status of emptiness in the cache.
+
+### Input
+
+- `cap`     -- intersection of two convex sets
+- `isempty` -- new status of emptiness
+"""
+function set_isempty!(cap::Intersection, isempty::Bool)
+    return set_isempty!(cap.cache, isempty)
+end
+
+
 # --- LazySet interface functions ---
 
 
@@ -139,14 +215,13 @@ end
 
 function ρ_helper(d::AbstractVector{N},
                   cap::Intersection{N},
-                  algorithm::String,
-                  check_intersection::Bool;
+                  algorithm::String;
                   kwargs...) where N<:Real
     X = cap.X # compact set
     H = cap.Y # halfspace or hyperplane or line
 
     # if the intersection is empty => stop
-    if check_intersection && is_intersection_empty(X, H, false)
+    if is_intersection_empty(X, H, false)
         error("the intersection is empty")
     end
 
@@ -170,7 +245,6 @@ end
                         <:LazySet{N},
                         <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
       [algorithm]::String="line_search",
-      [check_intersection]::Bool=true,
       [kwargs...]) where N<:Real
 
 Return the support function of the intersection of a compact set and a
@@ -194,10 +268,6 @@ half-space/hyperplane/line in a given direction.
                          generated by the given direction `d` and the
                          hyperplane's normal vector `n`
 
-- `check_intersection` -- (optional, default: `true`) if `true`, check if the
-                          intersection is empty before actually calculating the
-                          support function
-
 ### Output
 
 The scalar value of the support function of the set `cap` in the given
@@ -207,9 +277,6 @@ direction.
 
 It is assumed that the set `cap.X` is compact.
 
-The `check_intersection` flag can be useful if it is known in advance that the
-intersection is non-empty.
-
 Any additional number of arguments to the algorithm backend can be passed as
 keyword arguments.
 
@@ -237,9 +304,8 @@ function ρ(d::AbstractVector{N},
                              <:LazySet{N},
                              <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
            algorithm::String="line_search",
-           check_intersection::Bool=true,
            kwargs...) where N<:Real
-    return ρ_helper(d, cap, algorithm, check_intersection, kwargs...)
+    return ρ_helper(d, cap, algorithm; kwargs...)
 end
 
 # symmetric method
@@ -248,9 +314,8 @@ function ρ(d::AbstractVector{N},
                              <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                              <:LazySet{N}};
            algorithm::String="line_search",
-           check_intersection::Bool=true,
            kwargs...) where N<:Real
-    return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection, kwargs...)
+    return ρ_helper(d, cap.Y ∩ cap.X, algorithm; kwargs...)
 end
 
 """
@@ -320,9 +385,8 @@ function ρ(d::AbstractVector{N},
                              <:AbstractPolytope{N},
                              <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}}};
            algorithm::String="line_search",
-           check_intersection::Bool=true,
            kwargs...) where N<:Real
-    return ρ_helper(d, cap, algorithm, check_intersection; kwargs...)
+    return ρ_helper(d, cap, algorithm; kwargs...)
 end
 
 # symmetric method
@@ -331,9 +395,8 @@ function ρ(d::AbstractVector{N},
                              <:Union{HalfSpace{N}, Hyperplane{N}, Line{N}},
                              <:AbstractPolytope{N}};
            algorithm::String="line_search",
-           check_intersection::Bool=true,
            kwargs...) where N<:Real
-    return ρ_helper(d, cap.Y ∩ cap.X, algorithm, check_intersection; kwargs...)
+    return ρ_helper(d, cap.Y ∩ cap.X, algorithm; kwargs...)
 end
 
 """
@@ -370,9 +433,22 @@ Return if the intersection is empty or not.
 ### Output
 
 `true` iff the intersection is empty.
+
+### Notes
+
+The result will be cached, so a second query will be fast.
 """
 function isempty(cap::Intersection)::Bool
-    return is_intersection_empty(cap.X, cap.Y)
+    if isempty_known(cap)
+        # use cached result
+        return isempty(cap.cache)
+    end
+    # compute result
+    empty_intersection = is_intersection_empty(cap.X, cap.Y)
+    # update cache
+    set_isempty!(cap, empty_intersection)
+
+    return empty_intersection
 end
 
 
@@ -380,6 +456,7 @@ end
 # intersection of an array of sets
 # ================================
 
+
 """
     IntersectionArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
 
diff --git a/test/unit_Intersection.jl b/test/unit_Intersection.jl
index ff3e45de19..826fe777c1 100644
--- a/test/unit_Intersection.jl
+++ b/test/unit_Intersection.jl
@@ -18,6 +18,9 @@ for N in [Float64, Rational{Int}, Float32]
     @test ∈(ones(N, 2), I) && !∈(N[5, 5], I)
 
     # emptiness of intersection
+    @test !isempty_known(I)
+    @test !isempty(I)
+    @test isempty_known(I)
     @test !isempty(I)
 
     # =================