StandardAlgorithm returns rays that are not extremal
Hello,
I am currently trying to use sage's implementation of the double description algorithm to find all extremal rays of a given pointed cone $C={x: Ax \geq 0}$. My code looks as follows:
from sage.geometry.polyhedron.double_description import StandardAlgorithm, Problem
import numpy as np
def check_rays_extremal(res):
'''Check if all rays in res are actually extremal rays.'''
for i in range(len(res.R)):
if not res.is_extremal(res.R[i]):
print(f"ray{i} is not extremal")
return
def check_rays_projected_to_zero(res, A_bar):
'''Check if all extremal rays are actually projected onto 0 using A_bar.'''
for i in range(len(res.R)):
ray = res.R[i]
npar = np.array(ray)
mul = np.matmul(A_bar, npar)
max_entry = max(mul)
min_entry = min(mul)
if max_entry > 10**(-10) or min_entry < -10**(-10):
print(f"ray {it}: result is not zero: max value {max_entry}, min value {min_entry}")
return
def check_rays_geq_zero(res, A):
'''Check if all rays are in the pointed cone they are supposed to describe.'''
for i in range(len(res.R)):
ray = res.R[i]
npar = np.array(ray)
mul = np.matmul(A, npar)
min_entry = min(mul)
if min(mul) < -10**(-10):
print(f"ray {it}: result is not zero: min value {min(mul)}")
return
A = np.concatenate((A_bar, -A_bar, np.identity(A_bar.shape[1])))
A_sage = matrix(RR, A)
res = StandardAlgorithm(A_sage).run()
check_rays_extremal(res)
check_rays_projected_to_zero(res, A_bar)
check_rays_geq_zero(res, A):
Here, $\bar{A}$ is a numpy.matrix object with full rank and only integer entries (even though the datatype used is $A.dtype = float64$). To be precise, $A$ is an inequality description of a cone of shape ${x \geq 0: \bar{A}x = 0}$. The functions defined above are used to check the reasonability of $res.R$.
The results $res$ only have integer entries, but also have several properties that indicate wrong results:
- for some test matrices $A$, $checkRaysExtremal(res)$ identifies most of the rays in $res.R$ as being non-extremal
some entries of $res.R$ also neither satisfy $\bar{A}x = 0$ nor $Ax\geq 0$, hence they seem to not even be in the cone that is described by $A$
some items in $res.R$ are identical, i.e., some rays are generated multiple times
I have tried several things to fix the problem, among others:
- define A_sage over the ring $QQ$
- explicitly convert $A$ to a $np.matrix$ with $dtype=int64$
Nevertheless, I still can't fix the ascribed issues. Has one of you faced the same problem and managed to fix it? I am using SageMath 9.3 notebook on a Windows device.
If need be I can attach some more detailed information on inputs and outputs, but the matrices I am using usually have more than 200 rows and 100 columns.
Polyhedral geometry and floating point numbers are not really friends of each other.
I can suggest using Normaliz backend for polyhedral calculations.
$A$ only has integer entries between 0 and 10. Initializing it as an integer matrix via $A=np.matrix(..., dtype=int)$ did not solve the issue. I am going to try the Normaliz backend. Still, it seems off to me that sage returns wrong results for a well-conditioned, integer matrix. I am happy to provide further details into the input matrices if anyone wants to investigate the issue further.