Zero Volume for 14D polytope

[alhaas]

Dear polytope devs,

I hope this is the right place for this. I have created an example similar to the one in the tutorial for a convex polytope in H-representation, however it gives me a volume of 0.0 and also, there is no warning in case the dimensions do not match. I think I am missing something here, as the same problem gives a nonzero volume with the Volesti package in R.

My Python version is 3.7.6.

import numpy as np
import polytope as pt

A = np.array([[0.01281, -0.00435, -0.00303,  0.00307, -0.00161, -0.00153,
         0.00728, -0.00249,  0.02289,  0.00303,  0.00186, -0.01623,
        -0.00289, -0.00246],
       [ 0.01246, -0.00457, -0.0033 ,  0.00266, -0.00236, -0.00184,
         0.0068 , -0.00278,  0.0233 ,  0.00258,  0.00129,  0.11378,
        -0.00316, -0.00274],
       [ 0.03503, -0.0095 , -0.013  , -0.01412,  0.15911, -0.01669,
        -0.00333, -0.01422,  0.06059, -0.02439, -0.07807,  0.1101 ,
        -0.01338, -0.01452],
       [-0.03503,  0.0095 ,  0.013  ,  0.01412, -0.15911,  0.01669,
         0.00333,  0.01422, -0.06059,  0.02439,  0.07807, -0.1101 ,
         0.01338,  0.01452],
       [-0.05537,  0.01528,  0.01336, -0.00266,  0.19761,  0.01101,
        -0.0184 ,  0.01247, -0.09187,  0.00186,  0.02847, -0.12141,
         0.01314,  0.01252]])

b = np.array([62.2,   62.2, 1416. ,  978. , 1059.])

p = pt.Polytope(A, b)       
print(p.volume)

[johnyf]

Thank you for reporting this error. I confirm that zero volume is computed using commit f981f43. The reason appears to be a missing else branch at (within the computation of the bounding box, which is in turn used to compute the polytope's volume):

polytope/polytope/polytope.py

Lines 1304 to 1306 in f981f43

	if sol['status'] == 0:
	x = sol['x']
	l[i] = x[i]

The consequence of this missing branch is that if the optimization does not succeed, then the value of 0 is used for all bounds, because this is how the bounds are initialized within that function:

polytope/polytope/polytope.py

Line 1296 in f981f43

l = np.zeros([n, 1])

Running the lines of code from that area in ipython, when using scipy I get:

>>> sol
{'status': 3,
 'x': array([-6.37083e+07,  4.27238e+07,  4.38813e+06, -1.10456e+08,
        -1.77619e+08, -3.76229e+07, -1.46560e+08, -1.02404e+07,
        -1.70243e+08, -1.55426e+08, -3.99459e+08, -6.16730e+07,
         3.36853e+05, -7.07112e+08]),
 'fun': -707111604.4008977}

and when using cvxopt, I get:

>>> sol
{'status': 3, 'x': None, 'fun': None}

Quoting from the documentation of the function scipy.optimize.linprog:

3 : Problem appears to be unbounded.

and the documentation of the function cvxopt.solvers.lp says:

With the GLPK option, the solver does not return certificates of primal or dual infeasibility: if the status is 'primal infeasible' or 'dual infeasible', all entries of the output dictionary are None.

which is relevant to the lines:

polytope/polytope/solvers.py

Line 107 in f981f43

sol = cvx.solvers.lp(

polytope/polytope/solvers.py

Lines 115 to 116 in f981f43

	elif sol['status'] == 'dual infeasible':
	result['status'] = 3

So:

• the bounding box computation needs to be corrected to take into account solver failure
• the bounding box computation could be updated to handle unbounded polytopes, though that may better be discussed as a separate issue
• if the polytope described above (Zero Volume for 14D polytope #69 (comment)) is indeed unbounded, then its volume would be unbounded. So I presume that the Volesti package would be expected to return an infinite volume. However, I am not familiar with the Volesti package, and so I do not know whether it supports unbounded polytopes. The manual of Volesti does not mention unbounded polytopes. Its volume computation is either exact for zonotopes and simplices in V-representation (which are bounded sets), or approximate for other cases. Moreover, in the literature it is not unusual for "polytope" to be be defined as a bounded set, i.e., one that can be described in V-representation. A polytope in H-representation can be unbounded.

I had been aware that unbounded polytopes are currently unsupported by the package polytope, and intended to consider that topic in a separate issue at some point in the future.

In conclusion, the simplest approach for now seems to be to raise an exception from within the function polytope.polytope.bounding_box if no solution is found (until unbounded polytopes become supported).

Regarding the polytope given above (#69 (comment)), if indeed it is unbounded, then its volume is infinite.

[alhaas]

Thank you very much for the exhaustive answer. Indeed, the problem was unboundend, sorry for the confusion, I tried to prepare a minimal working example which ended up being too minimal as the bounds stayed in the pulp variable declaration.

I now created a file including bounds which results in a nonzero volume, however a different value is computed each time the script is executed. Is this behavior problem specific and/or is there a way to change the solving algorithm or it's parameters in order to better "suit" this problem?
example.txt

[necozay]

Volume algorithm implemented in polytope package is a randomized algorithm that tries to give an estimate of the volume. It might be unreliable especially in higher dimensions. You can try to increase the number of samples on the following line (this number should increase exponentially with the dimension to get similar accuracy):

polytope/polytope/polytope.py

Line 1467 in f981f43

N = 10000

to improve accuracy. This should reduce the variation but there will still be some variation due to randomness. Increasing the samples will slow down the computation. We currently don't have an alternative deterministic volume computation algorithm implemented.

[johnyf]

Thank you @alhaas for the example. Thank you @necozay for the comment.

I fixed the bug due to the missing else branch mentioned above (#69 (comment)) in commit 7f59645, which is now on branch master.

I added a parameter nsamples to the function polytope.polytope.volume in commit 2bd0a5b, which is now on branch dev. This can be used to set the number N mentioned above (#69 (comment)), which controls the behavior of the randomized algorithm. Also, I changed the function polytope.polytope.volume in commit bfc6b0b to not use the precomputed value of volume (this is already done within the property volume of the classes Polytope and Region).

So now the function volume can be called multiple times, and also after modifying a polytope (not recommended to modify a Polytope or Region, though possible because these classes are mutable), to recompute the volume (without the need to set _volume to None to trigger recomputation of the volume). There are a number of other changes too on branch dev.

I confirm that the example code given above (#69 (comment)) returns varying volume values in each run. Nonetheless, some of the values computed by different runs are: 8.071350273957078e+59, 7.371284178766923e+59, 9.059678878931413e+59. This is quite large an exponent. Even though the value varies, the variation is still a relatively small percentage of the computed value.

Using commit 5fb8a70 and changing the example's (#69 (comment)) last lines to:

p = pt.Polytope(A, b)
pt.volume(p, nsamples=10**6)  # the default value is 10**4
print(p.volume)

returns the value 8.524334217903648e+59 with variability < 2 %. So the changes on branch dev appear to address the issue for user code.

[alhaas]

Thanks again for the feedback, from my side the issue is solved. If possible, maybe it would be an enhancement for the future to also give a seed to the volume function, this way results should become reproducible.

[johnyf]

Thank you for the suggestion @alhaas. I implemented this approach by adding in commit 99e2090 a parameter named seed to the function polytope.polytope.volume.

In more detail, the randomness within the function polytope.polytope.volume originates from the function numpy.random.rand:

polytope/polytope/polytope.py

Line 1501 in 7f59645

+ np.random.rand(n, N)

The documentation of the function numpy.random.rand reads:

This is a convenience function for users porting code from Matlab, and wraps random_sample.

The function numpy.random.rand does not include a parameter for defining a seed. The documentation of the function numpy.random.random_sample reads:

New code should use the random method of a default_rng() instance instead;

The function numpy.random.random_sample does not include a parameter for defining a seed. The documentation of the function numpy.random.default_rng does include a parameter seed that can be an integer, among other types. Also relevant is the documentation of the class numpy.random.SeedSequence.

So the line of code linked-to above has been replaced by the line:

polytope/polytope/polytope.py

Line 1562 in 99e2090

+ np.random.default_rng(seed).random((n, N))

where seed is provided as an argument to the function polytope.polytope.volume.

Aside

The CI tests for Python 2.7 for commit 99e2090 failed, because the function numpy.random.default_rng does not exist, as it was added to numpy in a recent version.

The reason is that the last numpy version compatible with Python 2.7 was released about 1.5 years ago, on Dec 29, 2019 numpy == 1.16.6. Larger versions of numpy do not support Python 2.7 since Jul 26, 2019 (numpy == 1.17.0). I do not think there is any reason for polytope to remain compatible with Python 2.7. For that reason I opened #70.

[johnyf]

As of commit 9ee3f96, support for Python 2.7 has been removed, and also support for Python 3.5 and 3.6 has been removed. For motivation, please read the commit messages on branch dev.

The latest setuptools and numpy are now required on branch dev, the supported Python versions are 3.7, 3.8, 3.9 (, and soon 3.10), and the CI tests pass.

[johnyf]

As of commit a23ddcf, the changes discussed above have been merged to the mainline branch, which addresses this issue.

Support for unbounded polytopes is an orthogonal issue, and will be considered in a separate issue that I plan to open in time.