diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md index a77f594ae4..c8295505ce 100644 --- a/docs/src/lib/representations.md +++ b/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`. diff --git a/docs/src/lib/utils.md b/docs/src/lib/utils.md index 9c23be1553..cb655153ae 100644 --- a/docs/src/lib/utils.md +++ b/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 diff --git a/docs/src/man/set_operations.md b/docs/src/man/set_operations.md index ff9ca48e2c..498ad3c604 100644 --- a/docs/src/man/set_operations.md +++ b/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 | diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 45c75b4135..846e322321 100644 --- a/src/HPolyhedron.jl +++ b/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 # ======================================================== diff --git a/src/LinearMap.jl b/src/LinearMap.jl index 468dae5233..0b4d156371 100644 --- a/src/LinearMap.jl +++ b/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 diff --git a/src/compat.jl b/src/compat.jl index 150787b8fc..f1ee86fdd7 100644 --- a/src/compat.jl +++ b/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 diff --git a/src/helper_functions.jl b/src/helper_functions.jl index 908dc140ac..eb0ee89e02 100644 --- a/src/helper_functions.jl +++ 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) diff --git a/test/unit_LinearMap.jl b/test/unit_LinearMap.jl index 78572e70cf..481cce021d 100644 --- a/test/unit_LinearMap.jl +++ b/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) diff --git a/test/unit_Polyhedron.jl b/test/unit_Polyhedron.jl index efe5f84e0f..e48e84dbec 100644 --- a/test/unit_Polyhedron.jl +++ b/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 diff --git a/test/unit_Polytope.jl b/test/unit_Polytope.jl index e8845b94cf..1afe9c310c 100644 --- a/test/unit_Polytope.jl +++ b/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 # ----- diff --git a/test/unit_util.jl b/test/unit_util.jl index 63356d1db8..1a06f150d4 100644 --- a/test/unit_util.jl +++ b/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