Merge pull request #980 from JuliaReach/schillic/962

#962 #769 - Invertible linear map

docs/src/lib/representations.md

@@ -439,6 +439,7 @@ tosimplehrep(::HPoly{N}) where {N<:Real}
 tohrep(::HPoly{N}) where {N<:Real}
 isempty(::HPoly{N}) where {N<:Real}
 cartesian_product(::HPoly{N}, ::HPoly{N}) where {N<:Real}
+linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}}
 tovrep(::HPoly{N}) where {N<:Real}
 polyhedron(::HPoly{N}) where {N<:Real}
 remove_redundant_constraints
@@ -450,9 +451,6 @@ Inherited from [`LazySet`](@ref):
 * [`radius`](@ref radius(::LazySet, ::Real))
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 
-Inherited from [`AbstractPolytope`](@ref):
-* [`linear_map`](@ref linear_map(::AbstractMatrix{N}, ::AbstractPolytope{N}) where {N<:Real})
-
 #### Polytopes in constraint representation
 
 The following methods are specific for `HPolytope`.

docs/src/lib/utils.md

@@ -7,32 +7,30 @@ end
 
 # Utility functions
 
-```@docs
-sign_cadlag
-ispermutation
-@neutral
-@absorbing
-@declare_array_version
-```
-
 ## Helpers for internal use only
 
 ### Functions and Macros
 
 ```@docs
-@neutral_absorbing
-@array_neutral
-@array_absorbing
-cross_product(::AbstractMatrix{N}) where {N<:Real}
-get_radius!
 an_element_helper
-σ_helper
 binary_search_constraints
+cross_product(::AbstractMatrix{N}) where {N<:Real}
+get_radius!
+isinvertible
+ispermutation
 nonzero_indices
 samedir
+sign_cadlag
 _random_zero_sum_vector
 remove_duplicates_sorted!
 reseed
+σ_helper
+@neutral
+@absorbing
+@neutral_absorbing
+@declare_array_version
+@array_neutral
+@array_absorbing
 ```
 
 ### Types

docs/src/man/set_operations.md

@@ -82,8 +82,8 @@ The table entries have the following meaning.
 | `EmptySet`                   | x   | i | x | x          | x | x       | x         |            | x    | x      | x        |
 | `HalfSpace`                  | x   | x | x | x          | x | x       | x         |            |      |        | i        |
 | `HPolygon`/`HPolygonOpt`     | i   | i | x | i          | i | i       | i         | i          |      |        | i        |
-| `HPolyhedron`                | x   | x | x | i          | x | x       | x         |            |      |        | i        |
-| `HPolytope`                  | x   | x | x | i          | x | x       | i         | i          |      |        | i        |
+| `HPolyhedron`                | x   | x | x | i          | x | x       | x         | x          |      |        | i        |
+| `HPolytope`                  | x   | x | x | i          | x | x       | i         | x          |      |        | i        |
 | `Hyperplane`                 | x   | x | x | x          | x | x       | x         |            |      |        | i        |
 | `Hyperrectangle`             | i   | i | i | i          | i | i       | i         | i          | i    | i      | i        |
 | `Interval`                   | x   | i | x | x          | x | i       | i         | i          | i    | i      | i        |

src/HPolyhedron.jl

@@ -2,8 +2,7 @@ using MathProgBase, GLPKMathProgInterface
 
 import Base: isempty,
              rand,
-             convert,
-             copy
+             convert
 
 export HPolyhedron,
        dim, σ, ∈,
@@ -16,6 +15,7 @@ export HPolyhedron,
        vertices_list,
        singleton_list,
        isempty,
+       linear_map,
        remove_redundant_constraints,
        remove_redundant_constraints!,
        constrained_dimensions
@@ -437,8 +437,10 @@ A new polyhedron obtained by removing the redundant constraints in `P`.
 See `remove_redundant_constraints!`. 
 """
 function remove_redundant_constraints(P::PT;
-                                      backend=GLPKSolverLP()) where {N, PT<:HPoly{N}}
-    return remove_redundant_constraints!(copy(P), backend=backend)
+                                      backend=GLPKSolverLP()
+                                     ) where {N, PT<:HPoly{N}}
+    return remove_redundant_constraints!(PT(copy(constraints_list(P))),
+                                         backend=backend)
 end
 
 """
@@ -499,6 +501,49 @@ function remove_redundant_constraints!(P::PT;
     return P
 end
 
+"""
+    linear_map(M::AbstractMatrix{N}, P::PT; [cond_tol=DEFAULT_COND_TOL]::Number)
+        where {N<:Real, PT<:HPoly{N}}
+
+Concrete linear map of a polyhedron in constraint representation.
+
+### Input
+
+- `M`        -- matrix
+- `P`        -- polyhedron in constraint representation
+- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
+                the matrix is invertible)
+
+### Output
+
+A polyhedron of the same type as the input (`PT`).
+
+### Algorithm
+
+If the matrix ``M`` is invertible (which we check with a sufficient condition),
+then ``y = M x`` implies ``x = \\text{inv}(M) y`` and we transform the
+constraint system ``A x ≤ b`` to ``A \\text{inv}(M) y ≤ b``.
+"""
+function linear_map(M::AbstractMatrix{N},
+                    P::PT;
+                    cond_tol::Number=DEFAULT_COND_TOL
+                   ) where {N<:Real, PT<:HPoly{N}}
+    if !isinvertible(M; cond_tol=cond_tol)
+        if P isa HPolyhedron
+            error("linear maps for polyhedra need to be invertible")
+        end
+        # use the implementation for general polytopes
+        return invoke(linear_map, Tuple{typeof(M), AbstractPolytope{N}}, M, P)
+    end
+    # matrix is invertible
+    invM = inv(M)
+    constraints = Vector{LinearConstraint{N}}(undef, length(constraints_list(P)))
+    for c in constraints_list(P)
+        push!(constraints, LinearConstraint(vec(c.a' * invM), c.b))
+    end
+    return PT(constraints)
+end
+
 # ========================================================
 # External methods that require Polyhedra.jl to be loaded
 # ========================================================

src/LinearMap.jl

@@ -201,13 +201,15 @@ function ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real}
 end
 
 """
-    isbounded(lm::LinearMap)::Bool
+    isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool
 
 Determine whether a linear map is bounded.
 
 ### Input
 
-- `lm` -- linear map
+- `lm`       -- linear map
+- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
+                the matrix is invertible)
 
 ### Output
 
@@ -216,12 +218,18 @@ Determine whether a linear map is bounded.
 ### Algorithm
 
 We first check if the matrix is zero or the wrapped set is bounded.
+If not, we perform a sufficient check whether the matrix is invertible.
+If the matrix is invertible, then the map being bounded is equivalent to the
+wrapped set being bounded, and hence the map is unbounded.
 Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
 """
-function isbounded(lm::LinearMap)::Bool
+function isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool
     if iszero(lm.M) || isbounded(lm.X)
         return true
     end
+    if isinvertible(lm.M; cond_tol=cond_tol)
+        return false
+    end
     return isbounded_unit_dimensions(lm)
 end
 

src/compat.jl

@@ -8,7 +8,7 @@ using Compat: copyto!, axes, argmax, @warn
 import Compat.String
 using Compat.LinearAlgebra
 import Compat.LinearAlgebra: norm, checksquare, LAPACKException,
-                             SingularException, ×
+                             SingularException, ×, cond
 import Compat.InteractiveUtils.subtypes
 export _At_mul_B
 

src/helper_functions.jl

@@ -1,3 +1,6 @@
+# default tolerance for matrix condition number (see 'isinvertible')
+const DEFAULT_COND_TOL = 1e6
+
 """
     sign_cadlag(x::N)::N where {N<:Real}
 
@@ -86,6 +89,32 @@ function ispermutation(u::AbstractVector{T}, v::AbstractVector{T})::Bool where T
     return true
 end
 
+"""
+    isinvertible(M::AbstractMatrix; [cond_tol]::Number=DEFAULT_COND_TOL)
+
+A sufficient check of a matrix being invertible (or nonsingular).
+
+### Input
+
+- `M`        -- matrix
+- `cond_tol` -- (optional, default: `DEFAULT_COND_TOL`) tolerance of matrix
+                condition
+
+### Output
+
+If the result is `true`, `M` is invertible.
+If the result is `false`, this function could not conclude.
+
+### Algorithm
+
+We check whether the
+[matrix condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices)
+`cond(M)` is below some prescribed tolerance.
+"""
+function isinvertible(M::AbstractMatrix; cond_tol::Number=DEFAULT_COND_TOL)
+    return cond(M) < cond_tol
+end
+
 """
     remove_duplicates_sorted!(v::AbstractVector)
 

test/unit_LinearMap.jl

@@ -53,9 +53,10 @@ for N in [Float64, Rational{Int}, Float32]
     @test lm1_copy.X == lm1.X
 
     # boundedness
-    @test isbounded(lm)
-    @test isbounded(zeros(N, 2, 2) * HalfSpace(N[1, 1], N(1)))
-    @test !isbounded(ones(N, 2, 2) * HalfSpace(N[1, 1], N(1)))
+    @test isbounded(ones(N, 2, 2) * Singleton(N[1, 2]))  # bounded set
+    @test isbounded(zeros(N, 2, 2) * HalfSpace(N[1, 1], N(1)))  # zero map
+    @test !isbounded(N[2 3; 1 2] * HalfSpace(N[1, 1], N(1)))  # invertible matrix
+    @test !isbounded(N[2 3; 0 0] * HalfSpace(N[1, 1], N(1)))  # singular matrix (expensive check)
 
     # isempty
     @test !isempty(lm)

test/unit_Polyhedron.jl

@@ -54,6 +54,9 @@ for N in [Float64, Rational{Int}, Float32]
     @test constrained_dimensions(
         HPolyhedron{N}([LinearConstraint(N[1, 0], N(1))])) == [1]
 
+    # concrete linear map with invertible matrix
+    linear_map(N[2 3; 1 2], p)
+
     if test_suite_polyhedra
         # conversion to and from Polyhedra's VRep data structure
         cl = constraints_list(HPolyhedron(polyhedron(p)))
@@ -69,6 +72,9 @@ for N in [Float64, Rational{Int}, Float32]
         @test !isempty(P)
         addconstraint!(P, LinearConstraint(N[-1, 0], N(-1)))  # x >= 1
         @test isempty(P)
+
+        # concrete linear map with noninvertible matrix throws an error
+        @test_throws ErrorException linear_map(N[2 3; 0 0], P)
     end
 end
 

test/unit_Polytope.jl

@@ -106,6 +106,12 @@ for N in [Float64, Rational{Int}, Float32]
     @test BallInf(N[0, 0], N(1)) ⊆ P
     @test !(BallInf(N[0, 0], N(1.01)) ⊆ P)
 
+    # concrete linear map
+    linear_map(N[2 3; 1 2], P)  # invertible matrix
+    if test_suite_polyhedra
+        linear_map(N[2 3; 0 0], P)  # noninvertible matrix
+    end
+
     # -----
     # V-rep
     # -----

test/unit_util.jl

@@ -14,4 +14,8 @@ for _dummy_ in 1:1 # avoid global variable warnings
         push!(vectors, copy(v))
     end
     @test vectors == [[1, 2, 3, 4], [1, 2, 3, 5], [1, 2, 4, 5], [1, 3, 4, 5], [2, 3, 4, 5]]
+
+    # invertible matrix
+    @test LazySets.isinvertible([2 3; 1 2])
+    @test !LazySets.isinvertible([2 3; 0 0])
 end