Memory leaks in libsingular polynomial evaluation and substitution

[simon-king-jena]

The following two examples leak memory (in sagemath from 8.9 to 9.3.rc5)

sage: R = PolynomialRing(ZZ, 'x', 50)
sage: d = {str(g): g for g in R.gens()}
sage: p = sum(d.values())
sage: import gc
sage: mem0 = get_memory_usage()
sage: for i in range(20):
....:     for _ in range(50):
....:         _ = p.subs(**d)
....:     _ = gc.collect()
....:     mem1 = get_memory_usage()
....:     if mem1 > mem0:
....:         print(i, mem1-mem0)
....:         mem0 = mem1

and

sage: R.<x, y> = ZZ[]
sage: p = (x + y)**100
sage: mem0 = get_memory_usage()
sage: for i in range(20):
....:     _ = p(x + y, y)
....:     _ = gc.collect()
....:     mem1 = get_memory_usage()
....:     if mem1 > mem0:
....:         print(i, mem1-mem0)

This was reported on sage-devel and ask.sagemath.org.

We fix .subs (the first example above) by modifying the appropriate singular call. We identified two possibly related memory leaks in singular

• singular issue #1089
• singular issue #1090

In the mean time, we use the non-satisfactory solution of going via naive substitution in Python. To do so, we moved the generic implementation on the base class MPolynomial.

See also #13447 also related to memory handling in singular.

Upstream: Reported upstream. No feedback yet.

CC: @nbruin @malb

Component: memleak

Keywords: libsingular polynomial memleak

Author: Vincent Delecroix, ​Dima Pasechnik

Branch/Commit: u/vdelecroix/27261 @ edf8847

Reviewer: Dima Pasechnik

Issue created by migration from https://trac.sagemath.org/ticket/27261

[simon-king-jena]

comment:1

For testing, I changed the code so that the id() of any newly created MPolynomial_libsingular instance is stored in some set, and is removed from the set as soon as it becomes garbage collected.

Result: In the above examples, the MPolynomial_libsingular instances are correctly garbage collected.

Conclusion: The memleak occurs in Sage's usage of libsingular. Apparently the underlying data of a polynomial are not correctly freed when it is deallocated.

[simon-king-jena]

comment:2

Or, another possibility: Some temporary libsingular data created during the computation of Q is not freed.

To detect this, I will try to make the example more fine-grained: Test if the leak occurs if the only arithmetic operation involved is _add_, _mul_ or __pow__, respectively. (X+Y)<sup>120*Y</sup>100 mixes these three operations.

[simon-king-jena]

comment:3

Interestingly, creating a copy of P from scratch does not leak:

sage: import gc
sage: R.<X,Y> = ZZ[]
sage: P = (X+Y)^120*Y^100
sage: mem0 = get_memory_usage()
sage: for i in range(500):
....:     Q = (X+Y)^120*Y^100
....:     del Q
....:     _ = gc.collect()
....:     mem1 = get_memory_usage()
....:     if mem1>mem0:
....:         print(i, mem1-mem0)
....:         mem0 = mem1
....:         
(0, 0.1484375)
sage: 

So, it seems that calling the polynomial (i.e., P(X,Y)) is what is leaking!

[simon-king-jena]

comment:4

Note that there previously was a leak in polynomial evaluation, as by comment in sage.libs.singular.polynomial.singular_polynomial_call - see #16958. So, perhaps it didn't get completely fixed?

[simon-king-jena]

comment:5

singular_polynomial_call copies the input data before doing the actual call. This, I think, is a waste of time. I tested that taking the copy is not causing the leak, but while we are at it, it could as well be improved, too.

[simon-king-jena]

comment:6

Also, singular_polynomial_call has an argument poly *(*get_element)(object) that is (at least in the Sage library) always assigned with the function MPolynomial_libsingular_get_element, which does nothing more than return the ._poly of an MPolynomial_libsingular.

So, basically the get_element argument is void. Should it be removed?

[embray]

comment:8

Ticket retargeted after milestone closed (if you don't believe this ticket is appropriate for the Sage 8.8 release please retarget manually)

[embray]

comment:9

As the Sage-8.8 release milestone is pending, we should delete the sage-8.8 milestone for tickets that are not actively being worked on or that still require significant work to move forward. If you feel that this ticket should be included in the next Sage release at the soonest please set its milestone to the next release milestone (sage-8.9).

[videlec]

comment:10

Stumble upon a similar problem using .subs()

sage: R.<X,Y> = ZZ[] 
....: import gc 
....: mem0 = get_memory_usage() 
....: for i in range(1000): 
....:     for _ in range(100): 
....:         _ = (X + Y).subs(X=X, Y=Y) 
....:     _ = gc.collect() 
....:     mem1 = get_memory_usage() 
....:     if mem1>mem0: 
....:         print(i, mem1-mem0) 
....:         mem0 = mem1                                                                                                                                                                      
78 0.1328125
100 0.1328125
121 0.12890625
142 0.12890625
164 0.1328125
172 2.0
185 0.12890625
206 0.12890625
228 0.1328125
...

See also https://ask.sagemath.org/question/29444/high-memory-usage-when-substituting-variables/

[videlec]

comment:11

Since __call__ calls subs() I gues the latter is the culprit.

[videlec]

Commit: ee9a54a

[videlec]

New commits:

ee9a54a

27261: fix memory leak in polynomial substitution

[videlec]

Changed dependencies from #13447 to none

[videlec]

Branch: u/vdelecroix/27261

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ee9a54a to 00b85c8

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

00b85c8

27261: doctest subs leak

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

bd4d41f

27261: doctest subs and `__call__` leak

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 00b85c8 to bd4d41f

[videlec]

Description changed:

--- 
+++ 
@@ -1,138 +1,36 @@
-The following leaks (I tested in Py3, but in Py2 it's the same):
+The following two examples leak memory (in sagemath from 8.9 to 9.3.rc5)

-sage: R.<X,Y> = ZZ[] -sage: P = (X+Y)^120*Y^100 +sage: R = PolynomialRing(ZZ, 'x', 50) +sage: d = {str(g): g for g in R.gens()} +sage: p = sum(d.values()) sage: import gc sage: mem0 = get_memoryusage() -sage: for i in range(1000): -....: Q = P(X,Y) -....: del Q +sage: for i in range(20): +....: for in range(50): +....: = p.subs(**d) ....: = gc.collect() ....: mem1 = get_memory_usage() -....: if mem1>mem0: +....: if mem1 > mem0: ....: print(i, mem1-mem0) ....: mem0 = mem1 -....:
-0 0.25 -98 0.12890625 -128 0.12890625 -157 0.12890625 -187 0.12890625 -216 0.12890625 -246 0.12890625 -275 0.12890625 -305 0.12890625 -329 2.0 -331 0.12890625 -360 0.12890625 -390 0.12890625 -420 0.1328125 -450 0.12890625 -479 0.12890625 -509 0.12890625 -539 0.1328125 -569 0.12890625 -599 0.1328125 -629 0.12890625 -659 0.1328125 -673 2.0 -689 0.12890625 -718 0.12890625 -748 0.12890625 -778 0.1328125 -808 0.12890625 -822 0.1953125 -866 0.12890625 -896 0.12890625 -925 0.12890625 -955 0.12890625 -984 0.12890625 + +and + + +sage: R.<x, y> = ZZ[] +sage: p = (x + y)**100 +sage: mem0 = get_memoryusage() +sage: for i in range(20): +....: = p(x + y, y) +....: _ = gc.collect() +....: mem1 = get_memory_usage() +....: if mem1 > mem0: +....: print(i, mem1-mem0)

-With #13447 in Py2, the leak is reduced, but still present:
+This was reported [on sage-devel](https://groups.google.com/forum/#!topic/sage-devel/iO1SzoW0kcw) and [ask.sagemath.org](https://ask.sagemath.org/question/29444/high-memory-usage-when-substituting-variables/).

-```
-sage: mem0 = get_memory_usage()
-sage: for i in range(1000):
-....:     Q = P(X,Y)
-....:     del Q
-....:     _ = gc.collect()
-....:     mem1 = get_memory_usage()
-....:     if mem1>mem0:
-....:         print(i, mem1-mem0)
-....:         mem0 = mem1
-....:         
-(147, 0.03125)
-(176, 0.12890625)
-(206, 0.12890625)
-(235, 0.12890625)
-(265, 0.12890625)
-(294, 0.12890625)
-(324, 0.12890625)
-(329, 2.0)
-(350, 0.12890625)
-(379, 0.12890625)
-(409, 0.12890625)
-(438, 0.12890625)
-(468, 0.12890625)
-(497, 0.12890625)
-(527, 0.12890625)
-(556, 0.12890625)
-(586, 0.12890625)
-(615, 0.12890625)
-(645, 0.12890625)
-(672, 0.1328125)
-(673, 2.0)
-(705, 0.12890625)
-(734, 0.12890625)
-(764, 0.12890625)
-(793, 0.12890625)
-(823, 0.12890625)
-(852, 0.12890625)
-(882, 0.12890625)
-(911, 0.12890625)
-(941, 0.12890625)
-(970, 0.12890625)
-```
+A way to get rid of this problem is to use naive substitution in Python and avoid using low-level calls to singular library.

-There is a worse leak (even with #13447):
-
-```
-sage: for i in range(1000):
-....:     Q = P(X+Y,Y)
-....:     del Q
-....:     _ = gc.collect()
-....:     mem1 = get_memory_usage()
-....:     if mem1>mem0:
-....:         print(i, mem1-mem0)
-....:         mem0 = mem1
-....:         
-(0, 0.546875)
-(2, 0.12890625)
-(3, 0.2578125)
-(4, 0.12890625)
-(5, 0.2578125)
-(6, 2.2578125)
-(7, 0.12890625)
-(8, 0.2578125)
-(9, 0.2578125)
-(10, 0.2578125)
-(11, 0.12890625)
-(12, 0.2578125)
-(13, 2.265625)
-(14, 0.12890625)
-(15, 0.2578125)
-(16, 0.2578125)
-(17, 0.12890625)
-(18, 0.2578125)
-(19, 0.2578125)
-(20, 0.12890625)
-(21, 2.265625)
-(22, 0.2578125)
-(23, 0.12890625)
-(24, 0.2578125)
-(25, 0.2578125)
-```
-
-This was reported [on sage-devel](https://groups.google.com/forum/#!topic/sage-devel/iO1SzoW0kcw)
+See also #13447 also related to memory handling in singular.

[videlec]

Author: Vincent Delecroix

[dimpase]

comment:16

Does gc.collect() trigger Singular's GC - does it even have such an option? It does its own memory management.

[videlec]

comment:17

Replying to @dimpase:

Does gc.collect() trigger Singular's GC - does it even have such an option? It does its own memory management.

What does this have to do with the problem? The gc.collect() can be removed from the snippets. They are here to guarantee that there is no undeleted temporary Python objects lying around and get consistent tests. You can crash sage with

sage: R = PolynomialRing(ZZ, 'x', 50)
sage: d = {str(g): g for g in R.gens()}
sage: p = sum(d.values())
sage: while True:
....:     _ = p.subs(**d)

[dimpase]

comment:18

The following 1-line patch fixes the leak from comment:10 and the 1st leak reported on the ticket.

--- a/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
+++ b/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
@@ -3620,9 +3620,8 @@ cdef class MPolynomial_libsingular(MPolynomial):
                 res_id = fast_map_common_subexp(from_id, _ring, to_id, _ring)
                 _p = res_id.m[0]

-                from_id.m[0] = NULL
+                p_Delete(&from_id.m[0], _ring)
                 res_id.m[0] = NULL
-
                 id_Delete(&from_id, _ring)
                 id_Delete(&res_id, _ring)

it does not fix the 2nd leak on the ticket, but I guess it might be that one needs to clean more of these from_id parts...

[dimpase]

comment:19

It'd be good to have more eyes on this, subs() code, preferably of people who wrote chunks of it.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from bd4d41f to b4e22a4

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

19881aa	27261: fix memory leak in polynomial substitution
b4e22a4	27261: doctest subs and `__call__` leak

[videlec]

comment:21

thx for your input!

[dimpase]

comment:22

Replying to @videlec:

Replying to @dimpase:

Does gc.collect() trigger Singular's GC - does it even have such an option? It does its own memory management.

What does this have to do with the problem? The gc.collect() can be removed from the snippets. They are here to guarantee that there is no undeleted temporary Python objects lying around and get consistent tests. You can crash sage with

sage: R = PolynomialRing(ZZ, 'x', 50)
sage: d = {str(g): g for g in R.gens()}
sage: p = sum(d.values())
sage: while True:
....:     _ = p.subs(**d)

with the patch I just proposed this can be run seemingly forever, without showing any memory increase after few minutes of CPU time.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

fbc22ce

27261: fix leak tests

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from b4e22a4 to fbc22ce

[videlec]

comment:24

Replying to @dimpase:

Replying to @videlec:

Replying to @dimpase:

Does gc.collect() trigger Singular's GC - does it even have such an option? It does its own memory management.

What does this have to do with the problem? The gc.collect() can be removed from the snippets. They are here to guarantee that there is no undeleted temporary Python objects lying around and get consistent tests. You can crash sage with

sage: R = PolynomialRing(ZZ, 'x', 50)
sage: d = {str(g): g for g in R.gens()}
sage: p = sum(d.values())
sage: while True:
....:     _ = p.subs(**d)

with the patch I just proposed this can be run seemingly forever, without showing any memory increase after few minutes of CPU time.

Very good. This is implemented in commit 19881aa.

[videlec]

Changed author from Vincent Delecroix to Vincent Delecroix, ​Dima Pasechnik

[dimpase]

comment:25

seems to be a similar bug, from_id, to_id stuff in singular_polynomial_call. Here is how to trigger it there, (needs --long)

--- a/src/sage/libs/singular/polynomial.pyx
+++ b/src/sage/libs/singular/polynomial.pyx
@@ -167,12 +167,12 @@ cdef int singular_polynomial_call(poly **ret, poly *p, ring *r, list args, poly
         sage: import gc
         sage: F.<a> = GF(7^2)
         sage: R.<x,y> = F[]
-        sage: p = x+2*y
+        sage: p = (x+2*y)^3
         sage: def leak(N):
         ....:     before = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
         ....:     gc.collect()
         ....:     for i in range(N):
-        ....:         _ = p(a, a)
+        ....:         _ = p(a+x, a)
         ....:     after = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
         ....:     return (after - before) * 1024   # ru_maxrss is in kilobytes

gives

File "src/sage/libs/singular/polynomial.pyx", line 187, in sage.libs.singular.polynomial.singular_polynomial_call
Failed example:
    for i in range(30):  # long time
        n = leak(10000)
        print("Leaked {} bytes".format(n))
        if n == 0:
            zeros += 1
            if zeros >= 6:
                break
        else:
            zeros = 0
Expected:
    Leaked...
    Leaked 0 bytes
    Leaked 0 bytes
    Leaked 0 bytes
    Leaked 0 bytes
    Leaked 0 bytes
Got:
    Leaked 6180864 bytes
    Leaked 3956736 bytes
    Leaked 3923968 bytes
    Leaked 3923968 bytes
...

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

3ffc6dd

27261: move call test in sage/libs/singular

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from fbc22ce to 3ffc6dd

[videlec]

comment:27

True. I moved my test over ZZ in the same file.

[dimpase]

comment:28

unfortunately I don't know how to fix the leak in singular_polynomial_call(). Trying to do something similar leads tp segfaults.

[videlec]

Description changed:

--- 
+++ 
@@ -31,6 +31,6 @@

 This was reported [on sage-devel](https://groups.google.com/forum/#!topic/sage-devel/iO1SzoW0kcw) and [ask.sagemath.org](https://ask.sagemath.org/question/29444/high-memory-usage-when-substituting-variables/).

-A way to get rid of this problem is to use naive substitution in Python and avoid using low-level calls to singular library.
+We fix `.subs` (the first example above) by modifying the appropriate singular call. For `__call__` however we use the non-satisfactory solution of going via naive substitution in Python.

 See also #13447 also related to memory handling in singular.

[videlec]

comment:29

Thanks for trying. I update the description. Should we agree on the current status? Should we make it for sage-9.3?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 3ffc6dd to 1bf94f4

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

1bf94f4

27261: coerce argument to parent if possible

[mwageringel]

comment:31

I believe this indicates that the problem is upstream, in fast_map_common_subexp. I have replaced the function by one that does nothing but return a copy of the original polynomial.

--- a/src/sage/libs/singular/polynomial.pyx
+++ b/src/sage/libs/singular/polynomial.pyx
@@ -173,8 +173,6 @@ cdef int singular_polynomial_call(poly **ret, poly *p, ring *r, list args, poly
         ....:     for i in range(N):
         ....:         _ = p(x, y)
         ....:         _ = p(x + y, y)
-        ....:         _ = p(1, -1)
-        ....:         _ = p(0, 0)
         ....:     gc.collect()
         ....:     after = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
         ....:     return (after - before) * 1024   # ru_maxrss is in kilobytes
@@ -185,7 +183,6 @@ cdef int singular_polynomial_call(poly **ret, poly *p, ring *r, list args, poly
         ....:     gc.collect()
         ....:     before = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss
         ....:     for i in range(N):
-        ....:         _ = p(a, a)
         ....:         _ = p(x, y)
         ....:         _ = p(x + y, y)
         ....:         _ = p(a + x, a)
@@ -228,7 +225,7 @@ cdef int singular_polynomial_call(poly **ret, poly *p, ring *r, list args, poly
     from_id.m[0] = p

     rChangeCurrRing(r)
-    cdef ideal *res_id = fast_map_common_subexp(from_id, r, to_id, r)
+    cdef ideal *res_id = dummy_map_common_subexp(from_id, r, to_id)
     ret[0] = res_id.m[0]

     # Unsure why we have to normalize here. See #16958
@@ -243,6 +240,12 @@ cdef int singular_polynomial_call(poly **ret, poly *p, ring *r, list args, poly

     return 0

+# can fail if the result is expected to be constant
+cdef ideal *dummy_map_common_subexp(ideal *from_id, ring *r, ideal *to_id):
+    cdef ideal *res_id = idInit(1, 1)
+    res_id.m[0] = p_Copy(from_id.m[0], r)
+    return res_id
+
 cdef int singular_polynomial_cmp(poly *p, poly *q, ring *r):
     """
     Compare two Singular elements ``p`` and ``q`` in ``r``.

Got:
    Leaked 0 and 0 bytes
    Leaked 0 and 0 bytes
    Leaked 0 and 0 bytes
    Leaked 0 and 0 bytes
    Leaked 0 and 0 bytes
    Leaked 0 and 0 bytes

[dimpase]

comment:32

How does one convince the upstream? Would writing a C++ program using libsingular suffice? Or a Singular example would be required? It could be that Singular does not use this function the way it's used here. In fact, in Singular it's only called at one place, once, in in maMapIdeal() And the latter is only called once:

./Singular/ipshell.cc:  v->data=maMapIdeal(IDIDEAL(w), src_ring, (ideal)theMap, currRing,nMap);

in Singular's shell, so it's probably not too hard to figure out the Singular command calling it.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

f90202c

27261: be more careful about parent of result

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 1bf94f4 to f90202c

[dimpase]

comment:34

Singular's subst() leaks memory very fast, I don't know why we expect anything better from libsingular. E.g.

ring r=0,(x,y),dp;
poly p=(x+y)^10;
while (1) {
subst(p,x,x+y);
}

takes 30 sec to reach many gigabytes.

[videlec]

Upstream: Reported upstream. No feedback yet.

[videlec]

Description changed:

--- 
+++ 
@@ -31,6 +31,6 @@

 This was reported [on sage-devel](https://groups.google.com/forum/#!topic/sage-devel/iO1SzoW0kcw) and [ask.sagemath.org](https://ask.sagemath.org/question/29444/high-memory-usage-when-substituting-variables/).

-We fix `.subs` (the first example above) by modifying the appropriate singular call. For `__call__` however we use the non-satisfactory solution of going via naive substitution in Python.
+We fix `.subs` (the first example above) by modifying the appropriate singular call. For `__call__`, we identified a memory leak in singular, see [singular issue #1089](https://github.com/Singular/Singular/issues/1089). In the mean time, we use the non-satisfactory solution of going via naive substitution in Python.

 See also #13447 also related to memory handling in singular.

[videlec]

comment:35

I opened an issue on the singular side https://github.com/Singular/Singular/issues/1089.

[videlec]

comment:36

Sage allows some weird things

sage: R.<x,y> = ZZ[]
sage: S.<q> = ZZ[]
sage: (x+y).subs(x=q)  # expected
Traceback (most recent call last):
...
TypeError: unsupported operand parent(s) for +:
  'Univariate Polynomial Ring in q over Integer Ring' and
  'Multivariate Polynomial Ring in x, y over Integer Ring'
sage: x.subs(x=q)  # why does it work!?
q