Revision history [back]

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.