Refactor Hilbert_basis() and replace a slow test #40387

orlitzky · 2025-07-08T10:18:39Z

Initially my goal was to replace one slow Hilbert_basis() test, but I did some refactoring along the way -- none of which really affects the performance of the Hilbert basis calculation. The commit message lists these changes.

To replace the test, I made up some examples and fed them to normaliz. If the sage answer agrees with the normaliz answer, they must both be right, right?

mantepse · 2025-07-08T13:41:32Z

I don't understand this. On my computer, in a fresh sage, current develop branch, I have

sage: cone = Cone([[1,2,3,4], [0,1,0,7], [3,1,0,2], [0,0,1,0]]).dual()
sage: timeit("cone.Hilbert_basis()")
5 loops, best of 3: 45.6 ns per loop

Are you saying that this is too slow? Or does this depend on certain packages which may or may not be installed?

src/sage/geometry/cone.py

orlitzky · 2025-07-08T14:04:10Z

I don't understand this. On my computer, in a fresh sage, current develop branch, I have
sage: cone = Cone([[1,2,3,4], [0,1,0,7], [3,1,0,2], [0,0,1,0]]).dual()
sage: timeit("cone.Hilbert_basis()")
5 loops, best of 3: 45.6 ns per loop
Are you saying that this is too slow? Or does this depend on certain packages which may or may not be installed?

The result is cached so running it in a loop and picking the fastest one makes it look very fast indeed. Runs 2,3,... are instantaneous.

When I run sage -t --long, this test raises a warning because it takes about 40s measured in CPU time. Your CPU may fare better (until we can normalize the CPU time, #33022) but in any case this test is relatively slow and does not (AFAIK, which is not very far) test anything that other faster examples cannot test.

mantepse · 2025-07-08T14:35:04Z

I don't understand this. On my computer, in a fresh sage, current develop branch, I have
sage: cone = Cone([[1,2,3,4], [0,1,0,7], [3,1,0,2], [0,0,1,0]]).dual()
sage: timeit("cone.Hilbert_basis()")
5 loops, best of 3: 45.6 ns per loop
Are you saying that this is too slow? Or does this depend on certain packages which may or may not be installed?
The result is cached so running it in a loop and picking the fastest one makes it look very fast indeed. Runs 2,3,... are instantaneous.

Oh, you are absolutely right. However, ...

When I run sage -t --long, this test raises a warning because it takes about 40s measured in CPU time. Your CPU may fare better (until we can normalize the CPU time, #33022) but in any case this test is relatively slow and does not (AFAIK, which is not very far) test anything that other faster examples cannot test.

... I find this hard to believe:

sage: cone = Cone([[1,2,3,4], [0,1,0,7], [3,1,0,2], [0,0,1,0]]).dual()
sage: timeit("cone.Hilbert_basis()", number=1, repeat=1)
1 loop, best of 1: 1.49 s per loop

This would mean that my CPU would be 30 times faster. cat /proc/cpuinfo says Intel(R) Core(TM) Ultra 5 125H. Is this such a good computer?

src/sage/geometry/cone.py

mantepse · 2025-07-08T16:51:21Z

src/sage/geometry/cone.py

+        else:
+            # Avoid the PointCollection overhead if nothing was
+            # added to the irreducible list beyond self.rays().
+            return self.rays()
+


Suggested change

else:

# Avoid the PointCollection overhead if nothing was

# added to the irreducible list beyond self.rays().

return self.rays()

# Avoid the PointCollection overhead if nothing was

# added to the irreducible list beyond self.rays().

return self.rays()

Good idea, thanks. I also made a similar change at the beginning of the function. Defining L = () unconditionally does not take any time and lets us eliminate one branch of the if/else.

mantepse · 2025-07-08T17:01:39Z

One tiny change that seems to make a slight difference is

diff --git a/src/sage/geometry/cone.py b/src/sage/geometry/cone.py
index bcaa578c339..d8aed5e222e 100644
--- a/src/sage/geometry/cone.py
+++ b/src/sage/geometry/cone.py
@@ -1685,7 +1685,7 @@ class ConvexRationalPolyhedralCone(IntegralRayCollection, Container, ConvexSet_c
         need_strict = region.endswith("interior")
         M = self.dual_lattice()
         for c in self._PPL_cone().minimized_constraints():
-            pr = M(c.coefficients()) * point
+            pr = M(*(c.coefficients())) * point
             if c.is_equality():
                 if pr != 0:
                     return False

It may be possible to speed up this line further, or perhaps ToricLattice_generic.__call__ by bypassing some checks, but I don't know enough about this.

orlitzky · 2025-07-09T12:46:37Z

One tiny change that seems to make a slight difference is...

I tried this but the results were inconsistent, it gets slower as the lattice gets bigger. I let myself get carried away and spent a few hours trying to optimize this method last night. I was only able to obtain a very small improvement by reorganizing the conditionals inside of the loop (see the latest commit). The speedup is consistent though.

src/sage/geometry/cone.py

orlitzky · 2025-07-09T13:04:06Z

... I find this hard to believe:
...
This would mean that my CPU would be 30 times faster. cat /proc/cpuinfo says Intel(R) Core(TM) Ultra 5 125H. Is this such a good computer?

A factor of 30x is not outrageous. Keep in mind that timeit() is measuring wall time and not CPU time. Wall time is inherently unreliable -- you can slow it down by watching a movie on your PC at the same time the test is running, speed it up by running the code in parallel, etc. The "slow test" warnings use CPU time to avoid some of that inconsistency, but CPU time can still vary without a normalizing factor (that no one has gotten around to yet).

CPU time will be affected by the options used to compile sage and its dependencies, the vector features (SSE, AVX, etc.) of the CPU, hardware mitigations for spectre and meltdown, and some other internal details of the processor. I have an old Core 2 Duo Thinkpad with all of the vulnerability mitigations turned on, and this test takes about 28s on it. The computer where it takes 40s is actually brand new, but it has 64 processors with the trade-off being that each of them is individually pretty slow.

So long as it leads to useful refactorings and performance improvements I'm not too worried about accidentally speeding up a test that might not technically be considered slow once we have the normalizing factor.

Some refactorings to the implementation of Hilbert_basis(): * Construct a cone L from our linear_subspace() so that "y in L" works as intended (currently we try to coerce y into a vector space in a try/except block). This is not any faster, but it makes the code easier to read. * Remove the irreducibles from "gens" as we construct it. * Negate a condition in a loop to avoid bailing out with "continue" as part of the normal control flow. * Use a boolean indicator to check if the list of irreducibles was modified, rather than recomputing its length.

We have one Hilbert_basis() test that is raising "slow test!" warnings at around 40s. Here we replace it with three tests, each of which runs relatively quickly. The trio completes in about 15s. Since I know very little about Hilbert bases, I have checked the results using Normaliz. For example, $ cat cone.in amb_space 4 cone 4 1 0 1 0 -1 0 1 0 0 1 1 0 0 -1 1 0 $ normaliz --HilbertBasis cone $ cat cone.out ... (the resulting basis is written to cone.out).

Tighten the whitespace around list/generator comprehensions, and simplify the control flow in two instances by eliminating an if/else branch. (Only one of these was suggested by the reviewer, but the other is in a similar spirit.) Thanks to Martin Rubey for the suggestions.

mantepse · 2025-07-09T13:25:44Z

One tiny change that seems to make a slight difference is...

I tried this but the results were inconsistent, it gets slower as the lattice gets bigger.

Ah, you are right! Calling a function with many arguments is more expensive.

mantepse · 2025-07-09T13:30:40Z

Interestingly, the new version is about 10% slower on my computer on the example above :-(

Edit: Hm, not sure. develop is the same speed, so it's probably noise.

orlitzky · 2025-07-09T13:55:42Z

If you use %timeit -c (built-in to ipython) instead of sage's own timeit(), it will use the CPU time and show you the standard deviation as well as the mean. It's a bit more reliable.

mantepse · 2025-07-09T14:54:26Z

The following seems to make a big difference, but is, apparently, not correct - there are failing tests in fan_morphism which I do not understand.

diff --git a/src/sage/geometry/cone.py b/src/sage/geometry/cone.py
index c34bc9242a9..377da5e6216 100644
--- a/src/sage/geometry/cone.py
+++ b/src/sage/geometry/cone.py
@@ -1684,8 +1684,9 @@ class ConvexRationalPolyhedralCone(IntegralRayCollection, Container, ConvexSet_c
             return False
         need_strict = region.endswith("interior")
         M = self.dual_lattice()
+        E = M.element_class
         for c in self._PPL_cone().minimized_constraints():
-            pr = M(c.coefficients()) * point
+            pr = E(M, [ZZ(e) for e in c.coefficients()]) * point
             if pr < 0:
                 return False
             elif pr > 0:

mantepse · 2025-07-09T14:55:05Z

src/sage/geometry/cone.py

                return False
+            elif pr > 0:


Suggested change

elif pr > 0:

if pr > 0:

I have dropped this commit as well for now. It turns out that the new version can be slower in some pathological cases (like the empty cone in a big space) where most of the constraints are equality and the equality constraints are listed first. It might be possible to work around, but I have other things I should be doing instead of trying to shave nanoseconds off of this method :)

orlitzky · 2025-07-09T19:59:49Z

The following seems to make a big difference, but is, apparently, not correct - there are failing tests in fan_morphism which I do not understand.

Quotient lattices have a different element_class even though they may have the same number of coordinates. The M() constructor does some work to dispatch correctly.

When constructing the subcone that represents the given cone's linear_subspace(), we don't need to check that the generators are valid or minimal -- fewer generators might work, but no subset will work.

When the cone K has no rays, the test "x in K" can be done quickly by checking if x is zero. This can be a significant improvement if the lattice is large, and risks wasting only as much time as it takes to compare an integer to zero (i.e. nothing compared to how long the rest of the containment test is going to take).

orlitzky · 2025-07-09T23:11:32Z

Two more small improvements:

The Cone() in Hilbert_basis() can use check=False because we know that the generators are good.
I noticed that _contains() does not have a special case for the trivial cone (the set {0}). This can be checked instantaneously and can save a lot of time when the lattice is large.

mantepse · 2025-07-10T09:50:14Z

Looks good, the failures seem to be unrelated. Thank you for your patience!

orlitzky · 2025-07-10T14:14:12Z

Thank you for the careful review!

sagemathgh-40387: Refactor Hilbert_basis() and replace a slow test Initially my goal was to replace one slow `Hilbert_basis()` test, but I did some refactoring along the way -- none of which really affects the performance of the Hilbert basis calculation. The commit message lists these changes. To replace the test, I made up some examples and fed them to normaliz. If the sage answer agrees with the normaliz answer, they must both be right, right? URL: sagemath#40387 Reported by: Michael Orlitzky Reviewer(s): Martin Rubey, Michael Orlitzky

github-actions bot added the s: needs review label Jul 8, 2025

orlitzky requested review from vbraun and user202729 July 8, 2025 10:19

orlitzky force-pushed the simpler-hilbert-basis branch from da05285 to 3542b29 Compare July 8, 2025 11:50

mantepse reviewed Jul 8, 2025

View reviewed changes