coercion, categories, and slow code

[Martin Albrecht]

Hi there,

the following code is really, really really (REALLY!) slow:

sage: IP = InfinitePolynomialRing(QQ)
sage: x = IP.gen()
sage: x10000 = x[10000]
sage: %time 1/2*x10000
CPU times: user 0.86 s, sys: 0.00 s, total: 0.86 s
Wall time: 0.87 s
1/2*x10000

Alright, so that was slow. Let's try something new:

sage: x10001 = x[10001]
sage: %time 1/2*x10000
CPU times: user 49.56 s, sys: 0.03 s, total: 49.60 s
Wall time: 50.62 s
1/2*x10000

50 seconds!

The second try is a bit faster (I assume some maps were already constructed):

sage: %time 1/2*x10000
CPU times: user 7.74 s, sys: 0.00 s, total: 7.74 s
Wall time: 7.86 s
1/2*x10000

I lost track of coercion/categories, thus any pointers how to fix this would
be greatly appreciated. Patches even more ;)

Cheers,
Martin

PS: This stuff slows down the linear programming stuff quite a bit.

--
name: Martin Albrecht
_pgp: http://pgp.mit.edu:11371/pks/lookup?op=get&search=0x8EF0DC99
_otr: 47F43D1A 5D68C36F 468BAEBA 640E8856 D7951CCF
_www: http://www.informatik.uni-bremen.de/~malb
_jab: martinr...@jabber.ccc.de

[Robert Bradshaw]

On Dec 1, 2009, at 10:02 AM, Martin Albrecht wrote:

> Hi there,
>
> the following code is really, really really (REALLY!) slow:
>
> sage: IP = InfinitePolynomialRing(QQ)
> sage: x = IP.gen()
> sage: x10000 = x[10000]
> sage: %time 1/2*x10000
> CPU times: user 0.86 s, sys: 0.00 s, total: 0.86 s
> Wall time: 0.87 s
> 1/2*x10000

Ouch, that is pretty bad. Looks like something changed in 4.2, it's
doing a huge amount of regular expressions stuff...

sage: %prun 1/2*x10000
160277 function calls in 1.002 CPU seconds
Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno(function)
2 0.719 0.359 1.000 0.500
infinite_polynomial_ring.py:595(__getitem__)
20002 0.133 0.000 0.189 0.000 {built-in method sub}
20002 0.047 0.000 0.059 0.000 re.py:229(_compile)
20002 0.033 0.000 0.281 0.000 re.py:144(sub)
20002 0.028 0.000 0.056 0.000 re.py:271(_subx)
20002 0.019 0.000 0.028 0.000 re.py:251(_compile_repl)
40004 0.016 0.000 0.016 0.000 {method 'get' of 'dict'
objects}
20018 0.006 0.000 0.006 0.000 {isinstance}
1 0.001 0.001 1.002 1.002 <string>:1(<module>)
5 0.000 0.000 0.000 0.000
infinite_polynomial_element.py:1053(__init__)
1 0.000 0.000 0.000 0.000
infinite_polynomial_element.py:1151(_add_)

- Robert

[Simon King]

Hi Martin!

On 1 Dez., 19:02, Martin Albrecht <m...@informatik.uni-bremen.de>
wrote:

> Hi there,
>
> the following code is really, really really (REALLY!) slow:

Well, the default implementation creates a polynomial ring with 10000
variables in the background when you say x[10000], and another ring
with 10001 variables if you then say x[10001].

My original motivation for the InfinitePolynomialRings was the
computation of Symmetric Gröbner Bases. But from a couple of recent
posts I learn that people use it for different applications. And for
these applications, the sparse implementation might be better. So,
perhaps one should change the default?

Anyway, it is not slow at all when doing it sparsely:

sage: IP = InfinitePolynomialRing(QQ,implementation='sparse')

sage: x = IP.gen()
sage: x10000 = x[10000]
sage: %time 1/2*x10000

CPU times: user 0.00 s, sys: 0.02 s, total: 0.02 s
Wall time: 0.04 s
1/2*x10000
sage: timeit('y=1/2*x10000')
625 loops, best of 3: 133 µs per loop

Here, you just have one variable in the background:

sage: x10000.polynomial().parent()
Univariate Polynomial Ring in x10000 over Rational Field

So, who votes for changing the default?

Cheers,
Simon

[Simon King]

Hi Robert!

On 1 Dez., 19:18, Robert Bradshaw <rober...@math.washington.edu>
wrote:
[...]

> Ouch, that is pretty bad. Looks like something changed in 4.2, it's  
> doing a huge amount of regular expressions stuff...

It does a lot of regular expressions in the background. I did not see
a way to do it better. However, the code did not change since 4.2,
AFAIK.
I really think it is because the current default implementation is not
good if you need big variable indices.

Cheers,
Simon

[Martin Albrecht]

On Tuesday 01 December 2009, Simon King wrote:
> Hi Martin!
>
> On 1 Dez., 19:02, Martin Albrecht <m...@informatik.uni-bremen.de>
>
> wrote:
> > Hi there,
> >
> > the following code is really, really really (REALLY!) slow:
>
> Well, the default implementation creates a polynomial ring with 10000
> variables in the background when you say x[10000], and another ring
> with 10001 variables if you then say x[10001].

Sure, but up to 50 seconds for a simple coercion seems way way too much even
in that case.

Martin

[Simon King]

Hi Martin!

On 1 Dez., 20:30, Martin Albrecht <m...@informatik.uni-bremen.de>
wrote:
[...]

> Sure, but up to 50 seconds for a simple coercion seems way way too much even
> in that case.

Agreed. Let's try to find out what happens here.

My first thought was that it is due to the huge polynomial rings. But
the following seems to indicate that it is really a problem with
coercion and categories, after all:

sage: sage: IP = InfinitePolynomialRing(QQ)
sage: sage: x = IP.gen()
sage: sage: x10000 = x[10000]
sage: sage: %time 1/2*x10000
CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s
Wall time: 0.00 s
1/2*x10000
sage: sage: x10001 = x[10001]
sage: %prun 1/2*x10000
11924289 function calls (11804267 primitive calls) in 46.270

CPU seconds

Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno
(function)

1 19.536 19.536 46.270 46.270
infinite_polynomial_element.py:1190(_mul_)
8 10.051 1.256 24.464 3.058 pushout.py:506(pushout)
40004 3.854 0.000 3.854 0.000 {method 'intersection'
of 'set' objects}
520078 2.978 0.000 4.041 0.000 {built-in method sub}
240115 2.410 0.000 3.178 0.000 category.py:74(__init__)
80012/40008 1.147 0.000 7.033 0.000 pushout.py:193
(__mul__)
520078 0.927 0.000 1.158 0.000 re.py:229(_compile)
120028 0.717 0.000 4.295 0.000 pushout.py:144(__init__)
520078 0.624 0.000 5.823 0.000 re.py:144(sub)
520078 0.566 0.000 1.063 0.000 re.py:271(_subx)
520078 0.317 0.000 0.497 0.000 re.py:251(_compile_repl)
1040203 0.314 0.000 0.314 0.000 {method 'get' of 'dict'
objects}
1240409 0.303 0.000 0.303 0.000 {isinstance}
40008 0.273 0.000 0.395 0.000 pushout.py:49(__init__)
160020 0.258 0.000 0.641 0.000 pushout.py:175(__cmp__)
...

So, most of the time is spent for "pushout.py" and "category.py". It
isn't explicitely used in InfinitePolynomialRing/Element, so, I guess
there is indeed some coercion problem.

Nevertheless, the regular expression business isn't good either. I'll
see what I can do -- recent sage-devel/sage-support threads indicated
some improvements.

Best regards,
Simon

[Nick Alexander]

> Nevertheless, the regular expression business isn't good either. I'll
> see what I can do -- recent sage-devel/sage-support threads indicated
> some improvements.

I'm certain you are aware, but there is an art to optimizing regular
expressions. It might be that a tuned regex is necessary, rather than
avoiding a regex altogether.

Nick

[Simon King]

Hi!

I just found that part of the problem is coercion -- or actually
conversion -- for classical polynomial rings:

sage: R1 = PolynomialRing(QQ,'x',10001)
sage: R2 = PolynomialRing(QQ,'x',10002)
sage: x10000 = R1('x10000')
sage: %prun a = R2(x10000)
160026 function calls in 5.073 CPU seconds

Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno
(function)

1 4.849 4.849 5.073 5.073 <string>:1(<module>)
20003 0.114 0.000 0.157 0.000 {built-in method sub}
20003 0.035 0.000 0.045 0.000 re.py:229(_compile)
20003 0.023 0.000 0.225 0.000 re.py:144(sub)
20003 0.022 0.000 0.043 0.000 re.py:271(_subx)
20003 0.013 0.000 0.020 0.000 re.py:251(_compile_repl)
40006 0.013 0.000 0.013 0.000 {method 'get' of 'dict'
objects}
20003 0.004 0.000 0.004 0.000 {isinstance}
1 0.000 0.000 0.000 0.000 {method 'disable' of
'_lsprof.Profiler' objects}

This is too long, I think. And it is done behind the scenes when
multiplying 1/2*x10000 after creating x[10001] in your example.

Of course, this does not explain everything, but still it could be
worth while to improve.

My plan would now be to analyse your example further and then open
various tickets, also for the polynomial conversion.

Cheers,
Simon

[Simon King]

On 1 Dez., 21:57, Nick Alexander <ncalexan...@gmail.com> wrote:
[...]

> I'm certain you are aware, but there is an art to optimizing regular  
> expressions.  It might be that a tuned regex is necessary, rather than  
> avoiding a regex altogether.

Can you suggest some manual that can teach me about optimizing regular
expressions?

Thanks,
Simon

[Simon King]

Hi!

Here is a first ticket, with patch: http://trac.sagemath.org/sage_trac/ticket/7578

It reduces the computation time of the original example of this thread
to the time that is needed to convert the underlying finite
polynomials. Hence, as long as the conversion does not improve, it
seems to me that the arithmetic of the infinite polynomials can not be
further improved.

Or, positively: There seems to be considerable room for improving the
polynomial conversion, and the infinite polynomials would directly
profit from any progress. But this should be a different ticket.

Cheers,
Simon

[William Stein]

WTF? Regular expressions?!?!

> --
> To post to this group, send an email to sage-...@googlegroups.com
> To unsubscribe from this group, send an email to sage-devel-...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/sage-devel
> URL: http://www.sagemath.org



--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org

[Mike Hansen]

On Wed, Dec 2, 2009 at 12:57 PM, William Stein <wst...@gmail.com> wrote:
> WTF?  Regular expressions?!?!

The following messages are probably relevant for the fast conversion
between singular polynomial rings:

On Sat, Oct 18, 2008 at 2:55 AM, Michael Brickenstein
<bricke...@mfo.de> wrote:
> In Singular the same thing is essentially done from the interpreter
> level by the more general command fetch.
> I had a look, what it does internally and came to the conclusion,
> that it just calls
> poly prCopyR(poly p, ring src_r, ring dest_r)
> in your simple case (same coefficient domains).
> So first, you should setup a new ring and
> then map the polynomial via
> prCopyR
>
> Michael

On Mon, Oct 20, 2008 at 8:43 PM, <han...@mathematik.uni-kl.de> wrote:
> if the monomial ordering is really the same,
> you may also use
> poly prCopyR_NoSort(poly p, ring src_r, ring dest_r)
> which avoids the sorting the polynomial after mapping each monomial.
> There are also corresponding routines for ideals
> (ideal idrCopyR(ideal id, ring src_r, ring dest_r),
> ideal idrCopyR_NoSort(ideal id, ring src_r, ring dest_r)
> )
>

--Mike

[William Stein]

On Wed, Dec 2, 2009 at 1:01 AM, Mike Hansen <mha...@gmail.com> wrote:
> On Wed, Dec 2, 2009 at 12:57 PM, William Stein <wst...@gmail.com> wrote:
>> WTF?  Regular expressions?!?!
>
> The following messages are probably relevant for the fast conversion
> between singular polynomial rings:

Thanks. There's also regular expression stuff done directly in
InfinitePolynomialRing, which I think doesn't involve singular
directly (see the __init__ method of InfinitePolynomialRing).

Oh, I made a ticket with a few small bugs I found in
InfinitePolynomialRing when playing around with it just now:
http://trac.sagemath.org/sage_trac/ticket/7580

>
> On Sat, Oct 18, 2008 at 2:55 AM, Michael Brickenstein
> <bricke...@mfo.de> wrote:
>> In Singular the same thing is essentially done from the interpreter
>> level by the more general command fetch.
>> I had a look, what it does internally and came to the conclusion,
>> that it just calls
>> poly prCopyR(poly p, ring src_r, ring dest_r)
>> in your simple case (same coefficient domains).
>> So first, you should setup a new ring and
>> then map the polynomial via
>> prCopyR
>>
>> Michael
>
> On Mon, Oct 20, 2008 at 8:43 PM,  <han...@mathematik.uni-kl.de> wrote:
>> if the monomial ordering is really the same,
>> you may also use
>> poly prCopyR_NoSort(poly p, ring src_r, ring dest_r)
>> which avoids the sorting the polynomial after mapping each monomial.
>> There are also corresponding routines for ideals
>> (ideal idrCopyR(ideal id, ring src_r, ring dest_r),
>> ideal idrCopyR_NoSort(ideal id, ring src_r, ring dest_r)
>> )
>>
>
> --Mike
>

[Simon King]

Hi!

On Dec 2, 6:40 am, William Stein <wst...@gmail.com> wrote:

> On Wed, Dec 2, 2009 at 1:01 AM, Mike Hansen <mhan...@gmail.com> wrote:
> > On Wed, Dec 2, 2009 at 12:57 PM, William Stein <wst...@gmail.com> wrote:
> >> WTF?  Regular expressions?!?!

There are regular expressions in InfinitePolynomialRing, but (at least
after applying my patch) I don't think that this is the problem here.

> > The following messages are probably relevant for the fast conversion
> > between singular polynomial rings:

Thanks! Is there a ticket already for improving conversions of
multivariate polynomial rings?

> Oh, I made a ticket with a few small bugs I found in
> InfinitePolynomialRing when playing around with it just now:
>    http://trac.sagemath.org/sage_trac/ticket/7580

Thanks! I'll look into it.

Cheers,
Simon

[Nathann Cohen]

Hello everybody !!!!

Concerning the use of InfinitePolynomialRing in Sage, it was discussed
in another thread and I since wrote a patch (#7561) to change it. As
mentionned, I need nothing of what this class has been built for, and
now that it is replaced with plain "var", it is a thousand times
faster ( a friend of mine needed to solve a LP with ~2000 variables
and could not even generate the LP because of this slowness in
InfinitePolynomialRing.

Besides the slowness, there was a bug coming for
InfinitePolynomialRing which I have not been able to reproduce yet
( the reason why there was no bug report for it ) when the LP got big
enough :some "string" problem prevented us from using the generators
anymore, and restarting Sage ( ? ) sometimes solved the problem.

If I can produce this bug again, you'll have a bug report for it, but
I you know the class by heart perhaps you will be helped to know that
you need a large number of variables to create this problem.

Thanks !!!!!!

Nathann

[William Stein]

On Wed, Dec 2, 2009 at 4:36 AM, Nathann Cohen <nathan...@gmail.com> wrote:
> Hello everybody !!!!
>
> Concerning the use of InfinitePolynomialRing in Sage, it was discussed
> in another thread and I since wrote a patch (#7561)  to change it. As
> mentionned, I need nothing of what this class has been built for, and
> now that it is replaced with plain "var", it is a thousand times
> faster ( a friend of mine needed to solve a LP with ~2000 variables
> and could not even generate the LP because of this slowness in
> InfinitePolynomialRing.
>
> Besides the slowness, there was a bug coming for
> InfinitePolynomialRing which I have not been able to reproduce yet
> ( the reason why there was no bug report for it ) when the LP got big
> enough :some "string" problem prevented us from using the generators
> anymore, and restarting Sage ( ? ) sometimes solved the problem.
>

I'm not surprised. Looking through the code, its use of strings and
regular expressions is fairly delicate -- I wouldn't use regular
expressions at all to implement the same functionality (and more). But
I'm not rewriting anything (too lazy), so I shouldn't complain.

By the way, I wrote "WTF? Regular expressions?!?!" above, which was
very rude. Please accept my apology for that.

William

[Simon King]

Hi William!

On 2 Dez., 17:26, William Stein <wst...@gmail.com> wrote:
...

> I'm not surprised.  Looking through the code, its use of strings and
> regular expressions is fairly delicate -- I wouldn't use regular
> expressions at all to implement the same functionality (and more). But
> I'm not rewriting anything (too lazy), so I shouldn't complain.

IIRC, I tried various other methods (without strings), but they were
all slower. However, I don't remember any concrete examples.
So, it would help me if you commented on http://trac.sagemath.org/sage_trac/ticket/7580
what I should try instead. I am overhauling the code anyway.

> By the way, I wrote "WTF?  Regular expressions?!?!" above, which was
> very rude.  Please accept my apology for that.

Sure!

Cheers,
Simon

[Simon King]

Hi!

On 2 Dez., 17:47, Simon King <simon.k...@nuigalway.ie> wrote:
[...]

> IIRC, I tried various other methods (without strings), but they were
> all slower. However, I don't remember any concrete examples.

> So, it would help me if you commented onhttp://trac.sagemath.org/sage_trac/ticket/7580

> what I should try instead. I am overhauling the code anyway.

Here is one, that is actually needed internally:

# Create a polynomial ring
# (the typical underlying finite polynomial ring of a densely
implemented InfinitePolynomialRing)
sage: Vars = ['x_'+str(n) for n in range(50)]+['y'+str(n) for n in
range(50)]
sage: R = PolynomialRing(QQ,Vars)
# Create a big random element
sage: p = R.random_element()
sage: p *= R.random_element()
sage: p *= R.random_element()
sage: p *= R.random_element()
sage: p *= R.random_element()
# Some RE to analyse the leading monomial
sage: import re
sage: RE = re.compile('([a-zA-Z0-9]+)_([0-9]+)\^?([0-9]*)')
sage: RE.findall(str(p.lm()))
[('x', '1', '2'),
('x', '13', '2'),
('x', '16', ''),
('x', '23', ''),
('x', '45', '')]
# Essentially the same information is contained in the exponent vector
# of the leading monomial. Say:
sage: zip(Vars,p.lm().exponents()[0])
[('x_0', 0),
('x_1', 2),
('x_2', 0),
('x_3', 0),
('x_4', 0),
('x_5', 0),
('x_6', 0),
...
# Now compare the timings:
sage: timeit('L = RE.findall(str(p.lm()))')
625 loops, best of 3: 23.8 µs per loop
sage: timeit('L = zip(Vars,p.lm().exponents()[0])')
625 loops, best of 3: 40.5 µs per loop

The obvious reason is that only few variables of R actually appear in
p.lm().
*But*: This is the typical situation. And this is why I chose to work
with regular expressions rather than "cleaner" methods.

Cheers,
Simon

[YannLC]

first, you might be interested by this:

L = zip(Vars,p.lm().exponents()[0].sparse_iter())

its faster but still not enough...
then you might look at http://trac.sagemath.org/sage_trac/ticket/7587,
apply it ( review it ;) )
and do

L = [(Vars[i],e) for i,e in enumerate(p.lm().exponents
(as_ETuples=False)[0][:50]) if e]

to get what you want.

Yann