#750 - Allow unbounded ρ for Hyperplane/HalfSpace

[schillic]

See #750.

This PR brings only the support function.
Intuitively, it is clear what the support vector should look like, but I do not see the general formula.

[mforets]

lgtm

[mforets]

Something that we could add is a note on the algorithm behind σ_helper.

It is a consequence of weak-duality in LPs: if the primal is unbounded then the dual is infeasible. Since there is only 1 constraint, the feasible set of the dual problem is a*y == d, y >= 0 with a and d being the normal vector of the halfspace and d the given direction respectively (the objective function is b*y). It is easy to see that this problem is infeasible whenever a is not parallel to d, or they are parallel with negative combination coefficient (dot product -1, after normalization).

[mforets]

Intuitively, it is clear what the support vector should look like, but I do not see the general formula.

Can't we just let:

function σ(d::AbstractVector{N}, hp::Hyperplane{N})::N where {N<:Real}
    v, unbounded = σ_helper(d, hp, error_unbounded=true)
    return v
end

EDIT: i just saw that this is precisely the function defined in lines 95-98 😄

I think it is good as is, no? For those cases that ρ returns Inf, σ is undefined.

[schillic]

I think it is good as is, no? For those cases that ρ returns Inf, σ is undefined.

We also have an implementation for HPolyhedron. There we use ±Inf in every unbounded dimension. So I thought we would do something in this line.
Since the support vector is not unique, to me every solution that is consistent with ρ should be fine.

[schillic]

Something that we could add is a note on the algorithm behind σ_helper.

Done.

[mforets]

Alright, i had forgotten that some components of the support vector of a HPolyhedron may return Inf, if it is unbounded. I don't know if this answer has some practical value. That is to say that seems natural to me to make error_unbounded=true by default.

If σ is unbounded, then a and d are not parallel. Does it make sense to turn the non-zero entries of d[i] to be +-Inf, like:

julia> d
3-element Array{Float64,1}:
 -1.0
  0.0
  2.0

julia> dInf = map(x -> x != 0 ? sign(x) * Inf : 0.0, d)
3-element Array{Float64,1}:
 -Inf  
    0.0
  Inf