Issues converting HPolyhedron to VPolytope

[chsahit]

What happened?

I am trying to convert between a 2-dimensional HPolyhedron in 3D ambient space and the corresponding VPolytope. If the 2D plane lies on the origin, the conversion fails (I get two points instead of 4). If I offset the plane a little bit the conversion succeeds. I'm not sure if this is simply an issue of numerical precision or not - but I doubt it because then the bug shouldn't be fixed by applying a small offset to the plane. Here's a minimal reproduction of the issue.

import numpy as np
from pydrake.all import HPolyhedron, VPolytope

A = np.array([[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]])

b_linear = np.array([0.54, -0.46, 0.04, 0.04, 0, 0])
H_linear = HPolyhedron(A, b_linear)
V_linear = VPolytope(H_linear).vertices().T
print(f"Vertex representation of square: \n {V_linear}")

b_affine = np.array([0.54, -0.46, 0.04, 0.04, 1e-4, -1e-4])
H_affine = HPolyhedron(A, b_affine)
V_affine = VPolytope(H_affine).vertices().T
print(f"Vertex representation of shifted square: \n {V_affine}")

(It may be possible that #21055 is related)

Version

1.27.0

What operating system are you using?

Ubuntu 22.04

What installation option are you using?

pip install drake

Relevant log output

No response

[chsahit]

Some additional notes: (1) just updated - the issue persists on version 1.28. (2) Given that this is a pip install I believe the underlying solver is Gurobi

[RussTedrake]

Thanks for reporting this. On my system (also ubuntu 22.04), I get the output

Vertex representation of unit square: 
 [[ 0.54  0.04  0.  ]
 [ 0.46  0.04  0.  ]
 [ 0.54 -0.04  0.  ]
 [ 0.46 -0.04  0.  ]]
Vertex representation of shifted unit square: 
 [[ 5.4e-01  4.0e-02  1.0e-04]
 [ 4.6e-01  4.0e-02  1.0e-04]
 [ 5.4e-01 -4.0e-02  1.0e-04]
 [ 4.6e-01 -4.0e-02  1.0e-04]]

That seems correct to me? Can you post the output that you're getting which you think is incorrect?

I've confirmed that we have H_linear.ambient_dimension() == 3, H_linear.IsEmpty() == False, and affine_hull.AffineDimension() == 2, so it's dropping into the "Next, handle the case where the HPolyhedron is not full dimensional." code by @cohnt .

I don't think Gurobi is in play here (this code is mostly based on QHull), unless something in the AffineSubspace calls are triggering it. But if you installed via pip then one thing we know is that it's not using Gurobi. Just in case... did you install the Mosek license?

[chsahit]

Hi Russ! The output I get on my machine is as follows:

Vertex representation of square:
 [[ 4.60000000e-01 -2.12132034e-11  1.73241161e-20]
 [ 4.99207921e-01  4.00000000e-02  1.76776938e-22]]
Vertex representation of shifted square:
 [[ 4.60000000e-01  4.00000000e-02  1.00000000e-04]
 [ 4.60000000e-01 -4.00000000e-02  9.99999998e-05]
 [ 5.40000000e-01  4.00000000e-02  1.00000001e-04]
 [ 5.40000000e-01 -4.00000000e-02  1.00000000e-04]]

Your output does seem correct. Also interestingly enough, my bug doesn't reproduce on deepnote. I did not install the mosek license (and i double checked that moseklm_license_file was unset).

[RussTedrake]

Looking more, AffineSubspace does call a solver in this code path (It should be a linear program). It should use the first available solver from this list : GurobiSolver, MosekSolver, ClpSolver, ClarabelSolver, .... That should be ClpSolver for any pip install that doesn't BYO mosek license.

I've confirmed that if I disable gurobi and mosek (to get clp) on my ubuntu 22.04 install-from-source, then I see the failure.

Vertex representation of square: 
 [[ 4.60000000e-01 -2.12132034e-11  1.73241161e-20]
 [ 4.99207921e-01  4.00000000e-02  1.76776938e-22]]
Vertex representation of shifted square: 
 [[ 4.60000000e-01  4.00000000e-02  1.00000000e-04]
 [ 4.60000000e-01 -4.00000000e-02  9.99999998e-05]
 [ 5.40000000e-01  4.00000000e-02  1.00000001e-04]
 [ 5.40000000e-01 -4.00000000e-02  1.00000000e-04]]

To summarize, on ubuntu 22.04:

• gurobi enabled => correct
• mosek enabled => correct
• both disabled => incorrect

[cohnt]

@RussTedrake @hongkai-dai now that I'm back in town, I'm happy to take the lead on fixing this.

[cohnt]

I've been investigating this issue, and I have 3 observations.

1. On my machine, calling the test case from python reveals the problem. But running it in C++ as a test case for VPolytope follows the expected behavior. See here. Should we be concerned about this?
2. The HPolyhedron->VPolytope constructor has an optional argument in C++ that wasn't exposed in the Python bindings. I've fixed that here.
3. Setting a looser tolerance fixes this problem. See here. Perhaps we should have a looser default tolerance? Should we try to change the algorithm so it has more intuitive behavior when the tolerance is set too tightly?

I can go ahead and PR these changes, but 1. and 3. above are still open questions in my mind.

[RussTedrake]

Is https://github.com/RobotLocomotion/drake/pull/21527 related in any way?

I'm surprised by 1. Yes, I think we should definitely make sure we understand that.

[cohnt]

Don't think that PR (or its associated issue) is connected, and I've verified the behavior of the two test cases (with and without Gurobi/Mosek) is the same since the merge.

[AlexandreAmice]

I don't think 1. is a problem, I can replicate it in C++ fine on my machine.

The bug happens because Mosek and Gurobi pick the affine hull basis to be [1, 0] C = [0, 1] [0, 0]

The result is that after changing coordinates to the affine hull, the last two rows of the polytope read 0 = 0 exactly.

Meanwhile CLP chooses the basis [-1/√2, -1/√2] [-1/√2, 1/√2] [0, 0].

However, since 1/√2 is not rational, the [0,0] get set to -1e-11, 1e-9. Because of this (after rescaling) the last two rows polytope in the affine hull coordinates read [-1, -100] [x] <= [0] [1, 100] [y] [0] Which is a line.