ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sat, 24 Feb 2024 12:41:12 +0100StandardAlgorithm returns rays that are not extremalhttps://ask.sagemath.org/question/76187/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.Fri, 23 Feb 2024 12:06:38 +0100https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/Comment by tun27 for <p>Hello,</p>
<p>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:</p>
<pre><code>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):
</code></pre>
<p>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$.</p>
<p>The results $res$ only have integer entries, but also have several properties that indicate wrong results:</p>
<ul>
<li>for some test matrices $A$,
$checkRaysExtremal(res)$ identifies
most of the rays in $res.R$ as being
non-extremal</li>
<li><p>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$</p></li>
<li><p>some items in $res.R$ are identical,
i.e., some rays are generated
multiple times</p></li>
</ul>
<p>I have tried several things to fix the problem, among others:</p>
<ul>
<li>define A_sage over the ring $QQ$</li>
<li>explicitly convert $A$ to a
$np.matrix$ with $dtype=int64$</li>
</ul>
<p>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.</p>
<p>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.</p>
https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76216#post-id-76216$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.Sat, 24 Feb 2024 12:41:12 +0100https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76216#post-id-76216Comment by Max Alekseyev for <p>Hello,</p>
<p>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:</p>
<pre><code>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):
</code></pre>
<p>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$.</p>
<p>The results $res$ only have integer entries, but also have several properties that indicate wrong results:</p>
<ul>
<li>for some test matrices $A$,
$checkRaysExtremal(res)$ identifies
most of the rays in $res.R$ as being
non-extremal</li>
<li><p>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$</p></li>
<li><p>some items in $res.R$ are identical,
i.e., some rays are generated
multiple times</p></li>
</ul>
<p>I have tried several things to fix the problem, among others:</p>
<ul>
<li>define A_sage over the ring $QQ$</li>
<li>explicitly convert $A$ to a
$np.matrix$ with $dtype=int64$</li>
</ul>
<p>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.</p>
<p>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.</p>
https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76212#post-id-76212I can suggest using Normaliz backend for polyhedral calculations.Sat, 24 Feb 2024 05:22:32 +0100https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76212#post-id-76212Comment by FrédéricC for <p>Hello,</p>
<p>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:</p>
<pre><code>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):
</code></pre>
<p>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$.</p>
<p>The results $res$ only have integer entries, but also have several properties that indicate wrong results:</p>
<ul>
<li>for some test matrices $A$,
$checkRaysExtremal(res)$ identifies
most of the rays in $res.R$ as being
non-extremal</li>
<li><p>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$</p></li>
<li><p>some items in $res.R$ are identical,
i.e., some rays are generated
multiple times</p></li>
</ul>
<p>I have tried several things to fix the problem, among others:</p>
<ul>
<li>define A_sage over the ring $QQ$</li>
<li>explicitly convert $A$ to a
$np.matrix$ with $dtype=int64$</li>
</ul>
<p>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.</p>
<p>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.</p>
https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76205#post-id-76205Polyhedral geometry and floating point numbers are not really friends of each other.Fri, 23 Feb 2024 19:50:11 +0100https://ask.sagemath.org/question/76187/standardalgorithm-returns-rays-that-are-not-extremal/?comment=76205#post-id-76205