StandardAlgorithm returns rays that are not extremal

asked 2024-02-23 12:06:38 +0200

tun27 gravatar image

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.

edit retag flag offensive close merge delete

Comments

Polyhedral geometry and floating point numbers are not really friends of each other.

FrédéricC gravatar imageFrédéricC ( 2024-02-23 19:50:11 +0200 )edit

I can suggest using Normaliz backend for polyhedral calculations.

Max Alekseyev gravatar imageMax Alekseyev ( 2024-02-24 05:22:32 +0200 )edit

$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.

tun27 gravatar imagetun27 ( 2024-02-24 12:41:12 +0200 )edit