#1500 - Option to apply linear_map in H-rep without invertibility check

[schillic]

Closes #1500.

CC @ueliwechsler

[mforets]

The logic of the PR looks good to me. On the other hand, after successive extensions to the concrete linear map interface, i think it has become rather obscure... no?

I suggest to either:

(i) add more examples that illustrate the different use cases, or
(ii) make slight changes of the keyword names.

Thoughts for point (ii):

Problem:

To compute linear_map(A, X) where X is a polyhedral set. There are two ways:

• vrep
  • slower
  • can always be performed if X is polyhedral
• hrep
  • faster
  • A needs to be invertible
    • needs to check invertibility
    • needs to compute inverse

Keywords:

• use_vrep: if the argument is true, use vrep, otherwise try to use hrep or fallback to vrep if one of the preconditions for using the hrep is not fullfilled
• check_invertibility: argument to control whether or not the invertibility check is performed. this is useful if you know that the matrix is invertible (but have not computed its inverse).
• inverse: argument to pass the matrix inverse, if it is known. in this case, the linear map matrix M can be ommited (nothing)

Additional arguments:

• use_inv: use inv or not
• cond_tol: tolerance for invertibility check

Defaults:

• use_vrep is false by default, since we want to try the hrep first (it is relatively fast to check if hrep applies or not)
• check_invertibility is true by default
• inverse is nothing by default (i dont know the inverse a priori)

Additional arguments:

• use_inv:
• cond_tol: use the constant DEFAULT_COND_TOL

Pseudocode:

function linear_map(M::Union{AbstractMatrix{N}, Nothing},
                    P::AbstractPolyhedron{N};
                    use_vrep::Bool=false,
                    inverse::Union{Nothing, AbstractMatrix{N}}=nothing,
                    check_invertibility::Bool=(inverse == nothing),
                    cond_tol::Number=DEFAULT_COND_TOL,
                    use_inv::Bool=(inverse != nothing || !issparse(M))
                   ) where {N<:Real}
    @assert # check that dimensions match, should work if M is Nothing as well, by checking the dimensions of inverse

    # if use_vrep is activated, the rest of the arguments can be safely ignored
    use_vrep && return _linear_map_vrep(M, P)
    
    is_invertible = true
    if check_invertibility    # check if we can invert
        is_invertible = isinvertible(M; cond_tol=cond_tol)
    end

    if is_invertible
       # if we can invert => use hrep
        return _linear_map_hrep(M, P, use_inv; inverse=inverse)
    else
        # otherwise, use vrep
        return _linear_map_vrep(M, P)
    end
end

Comments on the pseudocode:

• there is a possible source of confusion if the user passes linear_map(A, M, use_vrep=false) since this does not necessarily use hrep: it depends on A being invertible. this should be documented.
• the code gives more control on the algorithm than the one in master
• if the inverse is known, one can pass nothing in the field M
• if inverse is passed explicitly, one doesn't need to check for invertibility hence check_invertibility is false by default in that setting

[schillic]

I also thought about this but did not want to break the current API. But since this seems okay, let me add another proposal. Since passing the inverse is a special use case, I think we can add a new method for this purpose.

function linear_map(M::Union{AbstractMatrix{N}, Nothing},
                    P::AbstractPolyhedron{N};
                    inverse::AbstractMatrix) where {N<:Real}
    return _linear_map_hrep(M, P, true; inverse=inverse)
end

function linear_map(M::AbstractMatrix{N}, P::AbstractPolyhedron{N};
                    algorithm::Symbol=(issparse(M) ? :division : :inverse),
                    check_invertibility::Bool=true,
                    cond_tol::Number=DEFAULT_COND_TOL) where {N<:Real}
    @assert dim(P) == size(M, 2)
    @assert algorithm in [:vrep, :inverse, :division]

    # check invertibility
    if algorithm != :vrep && check_invertibility && !isinvertible(M; cond_tol=cond_tol)
        algorithm = :vrep
    end

    if algorithm == :vrep
        return _linear_map_vrep(M, P)
    else
        use_inv = algorithm == :inverse
        return _linear_map_hrep(M, P, use_inv)
    end
end

[schillic]

We agreed to use the version in #1713 (comment) with strings instead of symbols.

[schillic]

Unfortunately the idea with the two methods does not work. For some reason, Julia does not find the method. It only works when passing nothing for M.

julia> using LazySets

julia> M = rand(2, 2)
2×2 Array{Float64,2}:
 0.0509302  0.107692
 0.0785203  0.835182

julia> linear_map(nothing, Singleton([1., 1.]); inverse=inv(M))
HPolygon{Float64}(HalfSpace{Float64,VN} where VN<:AbstractArray{Float64,1}[HalfSpace{Float64,Array{Float64,1}}([-2.3040029056097, 1.4944331500210524], 1.0), HalfSpace{Float64,Array{Float64,1}}([-24.50655542281636, 3.1599899249110077], -1.0), HalfSpace{Float64,Array{Float64,1}}([2.3040029056097, -1.4944331500210524], -1.0), HalfSpace{Float64,Array{Float64,1}}([24.50655542281636, -3.1599899249110077], 1.0)])

julia> linear_map(M, Singleton([1., 1.]); inverse=inv(M))
ERROR: MethodError: no method matching linear_map(::Array{Float64,2}, ::Singleton{Float64,Array{Float64,1}}; inverse=[24.50655542281636 -3.1599899249110077; -2.3040029056097 1.4944331500210524])

EDIT: When I replace Union{AbstractMatrix{N}, Nothing} by AbstractMatrix{N}, then Julia complains that it overrides a method. So apparently the keyword argument, even if required, does not matter for the dispatch.

Should we instead assign a new name to the "inverse" function?
For instance linear_map_inverse(M, P)

[mforets]

ok then i tend to prefer having a single linear_map function with several use cases. passing the inverse is just one use case, as there are others. we give the possibility of passing nothing if the inverse is known and you don't want to pass the matrix M itself (or don't want to compute it). the pseudocode i gave in comment #1713 (comment) covers these cases.

[schillic]

Finally the build passes.