Class representation for series expansions

[Shivam Vats]

Hi

Sympy currently lacks a class based representation for series expansions.
Though the current approach works well, it has its problems, like- all types of
series are lumped under one name, issues with infinite series, etc. Moreover, such
an implementation, if successful, can also be ported to C++, which will be much
faster. Hence, I would like to give series expansion a class representation.

There has been a lot of dicussion here :
https://groups.google.com/forum/#!searchin/sympy/series/sympy/hiRuIHa8ImA/JLOBsMr9yUcJ

This approach builds on ring_series.py and Ondrej's implementation here
https://github.com/certik/sympy/blob/59492520b443ea5f0ef31fc018e9bc700b93b818/a.py

is even faster than the existing one based on ring_series.

I found another approach here https://github.com/sympy/sympy/wiki/UD-Sequences-and-formal-power-series-prototype

which uses sequences to construct power series. There was an attempt to follow
this approach in https://github.com/sympy/sympy/pull/7846 using formula based
power series, but the PR is yet to be merged.

How do the two approaches compare? What are the expectations from such a class
based representation?

Regards

[Ondřej Čertík]

Hi Shivam,

On Sun, Mar 8, 2015 at 5:13 AM, Shivam Vats <shiv...@gmail.com> wrote:
> Hi
>
> Sympy currently lacks a class based representation for series expansions.
> Though the current approach works well, it has its problems, like- all types
> of
> series are lumped under one name, issues with infinite series, etc.

And it is slow.

> Moreover, such
> an implementation, if successful, can also be ported to C++, which will be
> much
> faster.

It should be done in SymPy first, it will speed it up a lot, as well
as allow us to figure out the best design. Then we should translate to
C++ to get additional speed.

> Hence, I would like to give series expansion a class representation.
>
> There has been a lot of dicussion here :
> https://groups.google.com/forum/#!searchin/sympy/series/sympy/hiRuIHa8ImA/JLOBsMr9yUcJ
>
> This approach builds on ring_series.py and Ondrej's implementation here
> https://github.com/certik/sympy/blob/59492520b443ea5f0ef31fc018e9bc700b93b818/a.py
>
> is even faster than the existing one based on ring_series.

It's a little bit faster, but I think that's just because I sort the
dictionaries first:

https://github.com/certik/sympy/blob/59492520b443ea5f0ef31fc018e9bc700b93b818/a.py#L114

This can be done easily in ring_series as well. The two approaches are
equivalent, as far as I can tell.

>
> I found another approach here
> https://github.com/sympy/sympy/wiki/UD-Sequences-and-formal-power-series-prototype
>
> which uses sequences to construct power series. There was an attempt to
> follow
> this approach in https://github.com/sympy/sympy/pull/7846 using formula
> based
> power series, but the PR is yet to be merged.

How fast is that approach?

>
> How do the two approaches compare? What are the expectations from such a
> class
> based representation?

I think we should use ring_series in a new Series class, and implement
any missing features in ring_series as needed. See the other thread
you cited above.

I think this is the most optimal way to represent and implement series
expansion. It's what Piranha uses (except that Piranha does some
clever optimizations in C++, but the overall structure, i.e. a
hashtable is the same). The only other way that I know of is to use
some kind of a dense representation (polys can do dense representation
too, but overall the sparse representation seems about as fast or
faster for most problems). So I would however use sparse
representation, i.e. ring_series. So just doing the above will provide
a very fast implementation, and it seems to me that any other
datastructure/design will only provide minor algorithmic speedup. Once
this is done, we should implement this in C++ and use some of the
tricks from Piranha, that will provide a good constant speedup.

You can try to implement the Kronecker trick with packing the
exponents in Python as well. That provides good speedup, but it's an
additional complexity.

Ondrej

[Shivam Vats]

Hi Ondrej

I tried to compare a naive implementation of Knonecker in Python with what Karatsuba that
Sympy currently uses for multiplication. I have modified  pernici's code for Kronecker.
https://github.com/shivamvats/sympy/commit/99b1fa5fecccb7a46b19d0863b4e8f3315f555da

```
res = 2*x+1
t1 = clock()
for i in range(1000):
    res = R.dup_pack_mul(2*x+1, res)
t2 = clock()
print(t2-t1)

31.3395597935
```

```res = 2*x+1
t1 = clock()
for i in range(1000):
    res = R.dup_mul(2*x+1, res)
t2 = clock()
print(t2-t1)

6.16460609436
```

Karatsuba, and even normal multiplication is faster than Kronecker by a large factor.
The main reason is of course, that we should use a fast library for the integer multiplication.
Moreover, how much padding needs to be done, is also not a trivial question.

So, I too, think we should use sparse representation to optimise Sympy
for no.

Shivam

 

[Fredrik Johansson]

On Thursday, March 12, 2015 at 3:05:28 AM UTC+1, Shivam Vats wrote:

Hi Ondrej

I tried to compare a naive implementation of Knonecker in Python with what Karatsuba that
Sympy currently uses for multiplication. I have modified  pernici's code for Kronecker.
https://github.com/shivamvats/sympy/commit/99b1fa5fecccb7a46b19d0863b4e8f3315f555da

```
res = 2*x+1
t1 = clock()
for i in range(1000):
    res = R.dup_pack_mul(2*x+1, res)
t2 = clock()
print(t2-t1)

31.3395597935
```

```res = 2*x+1
t1 = clock()
for i in range(1000):
    res = R.dup_mul(2*x+1, res)
t2 = clock()
print(t2-t1)

6.16460609436
```

Karatsuba, and even normal multiplication is faster than Kronecker by a large factor.
The main reason is of course, that we should use a fast library for the integer multiplication.
Moreover, how much padding needs to be done, is also not a trivial question.

The integer multiplication could be an issue, but on top of that the packing and unpacking is done with O(n^2) complexity.

Fredrik

[Shivam Vats]

Hi Fredrik

Indeed! Packin/unpacking when repeated 1000 times, takes too much time. And anyway,
we don't need to pack and unpack every time.

```
         res = 2*x+1
    ...: for i in range(1000):
    ...:     res = R.dup_mul(2*x+1, res)

       : t1 = clock()
    ...: res2 = R.dup_mul(res, res)
    ...: t2 = clock()
    ...: print(t2-t1)
    ...:
1.30838489532

       : t1 = clock()
    ...: res2 = R.dup_pack_mul(res, res)
    ...: t2 = clock()
    ...: print(t2-t1)
    ...:
0.579982995987
```
This should be a better comparison between the two methods. So, Kronecker does outperform
Karatsuba.

Thanks a lot for pointing out!

Regards
Shivam

[Ondřej Čertík]

Thanks for coding this up!

My experiments with Piranha show that packing the exponents speed
things up about 7x. It might be that packing is much faster in C++
than in Python.

Ondrej

[Shivam Vats]

Hi

I wrote a FormalPowerSeries class on top of `rings_series` here: FormalPowerSeries
Currently it can handle only univariate power series, but it can easily be
extended to multivariate and laurent series (by multiplying by a suitable
power of x).

I'm planning to have the following class structure :

• FormalSeries
• FormalPowerSeries
• FormaLaurentSeries
• FormalPuiseuxSeries

What do you think?

However, I am not so sure about how to integrate PuiseuxSeries with

ring_series.

Regards

Shivam

[Shivam Vats]

Hi

I've updated my proposal for Series expansion here GSoC Proposal

It is still  a work in progress. I'd be glad if you could review it.

Regards 

Shivam Vats

[Alexey U. Gudchenko]

21.03.2015 15:43, Shivam Vats wrote:
> Hi
> I've updated my proposal for Series expansion here GSoC Proposal

> <https://github.com/sympy/sympy/wiki/GSoC-2015-Application--Shivam-Vats>

> It is still a work in progress. I'd be glad if you could review it.
>
> Regards
> Shivam Vats
>

Hi, Shivam

I have updated the title of this page.

Does it plan to implement the multi variables series (Phase II or III)?
It is more complicated case, but it must be taken into account when
designing classes.


--
Alexey Gudchenko

> --
> You received this message because you are subscribed to the Google
> Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to sympy+un...@googlegroups.com
> <mailto:sympy+un...@googlegroups.com>.
> To post to this group, send email to sy...@googlegroups.com
> <mailto:sy...@googlegroups.com>.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/9c00e59f-e025-4057-9e8e-4504451e5926%40googlegroups.com
> <https://groups.google.com/d/msgid/sympy/9c00e59f-e025-4057-9e8e-4504451e5926%40googlegroups.com?utm_medium=email&utm_source=footer>.
> For more options, visit https://groups.google.com/d/optout.

[Shivam Vats]

Hi Alexey

Thanks a lot!

Does it plan to implement the multi variables series

Yes it does. The series class will be based on `ring_series` which, itself

makes use of monomials and tuples to represent terms. However, the 

currently implemented functions in it (like `rs_log`) are univariate. I will 

implement their multivariate series as well. However I am not yet sure of 

best way to do it. I have Modern Computer Arithmetic" by Richard Brent 

and Paul Zimmermann , which describes argument reduction for elementary

series. 

I realised I wasn't clear about it in the proposal. Thanks for pointing out.

I have created a PR https://github.com/sympy/sympy/pull/9182 for discussing

my proposal. I hope it makes the discussions easier.

Regards

Shivam Vats

[Shivam Vats]

Hi

What is your opinion about having different series classes (which

I have proposed) as opposed to having a single series class?

I have read your  ideas about series representation that you 

wrote here https://github.com/sympy/sympy/wiki/UD-Sequences-and-formal-power-series-prototype

Do you see any issue with using ring_series?

Regards

Shivam