diff --git a/src/liftedrep.jl b/src/liftedrep.jl index b1808e5d..8d1c1027 100644 --- a/src/liftedrep.jl +++ b/src/liftedrep.jl @@ -18,7 +18,17 @@ end FullDim(rep::LiftedHRepresentation) = size(rep.A, 2) - 1 similar_type(::Type{LiftedHRepresentation{S, MT}}, ::FullDim, ::Type{T}) where {S, T, MT} = LiftedHRepresentation{T, similar_type(MT, T)} -hvectortype(p::Type{LiftedHRepresentation{T, MT}}) where {T, MT} = vectortype(MT) +hvectortype(::Type{LiftedHRepresentation{T, MT}}) where {T, MT} = vectortype(MT) +fulltype(::Type{LiftedHRepresentation{T,MT}}) where {T,MT} = LiftedHRepresentation{T,MT} +dualtype(::Type{LiftedHRepresentation{T,MT}}) where {T,MT} = LiftedVRepresentation{T,MT} +function dualfullspace(::LiftedHRepresentation, d::FullDim, ::Type{T}) where {T} + N = fulldim(d) + return LiftedVRepresentation{T}( + [SparseArrays.spzeros(T, N) one(T) * I], + BitSet(1:N), + ) +end + LiftedHRepresentation{T}(A::AbstractMatrix{T}, linset::BitSet=BitSet()) where {T} = LiftedHRepresentation{T, typeof(A)}(A, linset) LiftedHRepresentation{T}(A::AbstractMatrix, linset::BitSet=BitSet()) where {T} = LiftedHRepresentation{T}(AbstractMatrix{T}(A), linset) diff --git a/src/projection.jl b/src/projection.jl index 3557b741..2a66b50c 100644 --- a/src/projection.jl +++ b/src/projection.jl @@ -20,6 +20,56 @@ struct FourierMotzkin <: EliminationAlgorithm end BlockElimination Computation of the projection by computing the H-representation and applying the block elimination algorithm to it. + +The idea is as follows. +Consider a H-representation given by ``Ax + By \\le c``. +For any ``d``, the projection on the `x` variables has +belongs to the halfspace ``d^\\top x \\le \\delta`` where ``\\delta`` +is +```math +\\begin{aligned} +\\max \\; & d^\\top x \\ +\\text{s.t.} & A x + By \\le c \\ +\\end{aligned} +``` +By duality, this is equal to +```math +\\begin{aligned} +\\min \\; & c^\\top \\nu \\ +\\text{s.t.} & A^\\top \\nu = d \\ +& B^\\top \\nu = 0 \\ +& \\nu \\ge 0 +\\end{aligned} +``` +Let `R` be the matrix for which the columns are the extreme rays of the +polyhedron with defined by +``B^\\top \\nu = 0, \\nu \\ge 0``. +The program can be rewritten as +```math +\\begin{aligned} +\\min \\; & c^\\top \\nu \\ +\\text{s.t.} & A^\\top \\nu = d \\ +& \\nu = R\\lambda \\ +& \\lambda \\ge 0 +\\end{aligned} +``` +which is equivalent to +```math +\\begin{aligned} +\\min \\; & c^\\top R \\lambda \\ +\\text{s.t.} & A^\\top R \\lambda = d \\ +& \\lambda \\ge 0 +\\end{aligned} +``` +Dualizing, we obtain +```math +\\begin{aligned} +\\max \\; & d^\\top x \\ +\\text{s.t.} & R^\\top A x \\le R^\\top c \\ +\\end{aligned} +``` +where we see that ``R^\\top A x \\le R^\\top c`` +is the H-representation of the polyhedron. """ struct BlockElimination <: EliminationAlgorithm end @@ -49,8 +99,19 @@ project(p::Polyhedron, pset, algo) = eliminate(p, setdiff(1:fulldim(p), pset), a supportselimination(p::Polyhedron, ::FourierMotzkin) = false eliminate(p::Polyhedron, delset, ::FourierMotzkin) = error("Fourier-Motzkin elimination not implemented for $(typeof(p))") -supportselimination(p::Polyhedron, ::BlockElimination) = false -eliminate(p::Polyhedron, delset, ::BlockElimination) = error("Block elimination not implemented for $(typeof(p))") +supportselimination(p::Polyhedron, ::BlockElimination) = true + +struct ProjectionMatrix <: AbstractMatrix{Bool} + indices::Vector{Int} + n::Int +end +Base.size() +Base.adjoint(P::ProjectionMatrix) = transpose(P) + +function eliminate(p::Polyhedron, delset, ::BlockElimination) + h_eliminate = ProjectionMatrix(sort(collect(delset)), fulldim(p)) \ hrep(p) + return h_eliminate +end eliminate(p::Polyhedron, algo::EliminationAlgorithm) = eliminate(p, BitSet(fulldim(p)), algo)