sympy process out of memory calculating series of asin
thomas...@gmail.com
I am using sympy to analyze and run some code involving power series. In general, the code is always running on a partial sum of a power series, and this is a normal sympy symbolic expression. After doing certain calculations, the data is no longer a series, and I need to run a Taylor expansion in order to get a series. For example, if I calculate the quotient of two partial sums, I need Taylor in order to get another partial sum. In sympy, it is very convenient to call series() to do this.
Now it is clear that this can use a lot of memory, because every additional derivative makes the expression more complex and longer. It appears as though series() does not use parallel processing, so the result is 100% usage of one CPU and increasing use of memory. On my laptop (running Windows 11 and Visual Studio), I have 16GB RAM and have sometimes run into memory issues.
To solve this problem, I have allocated an Ubuntu Linux virtual machine on the Internet (Microsoft Azure). The machine has 2 virtual CPUs and 218GB RAM. With this, I have successfully calculated partial sums up to order 24 or 31. This sometimes runs at 100% CPU allocation for more than 24 hours and completes successfully.
Up to now, I have used series() for algebraic expressions of partial sums (quotient, square root, etc.) but recently started applying this to a trigonometric expression. This is because I need to calculate the solution of a cubic equation and the Cardano formula has run into unpleasant situations involving complex numbers, and I have been given a recommendation to try the trigonometric formula instead. Up to now, it looks very good, since the primary solution is real and does not require complex components. By the way, solving a cubic equation worked will with both the standard formula and via solve(). For the cubic equation, I have started with the formula because I can see if each individual step success before proceeding to the next one.
Specifically, my code is
α = asin(α)
α = series(α, m, n=maxK+1) # calculate series up to maxK
On my laptop, the series statement runs for some time and hits the memory limit of my PC, so I tried running this on the VM, which has more memory and does not have the overhead of Windows or Visual Studio.
This first attempt ran for somewhere between 24 and 30 hours and the VM was stopped. This is presumably because the VM is configured to use a “spot price”, meaning that the administrators can stop it at any time they need additional capacity for something else.
The second attempt ran for somewhere between 30 and 48 hours and my output file ended in the statement killed. Running dmesg gave me this information:
[144763.491084] oom-kill:constraint=CONSTRAINT_NONE,nodemask=(null),cpuset=/,mems_allowed=0,global_oom,task_memcg=/user.slice/user-1000.slice/session-1.scope,task=python3.9,pid=1699,uid=1000
[144763.491101] Out of memory: Killed process 1699 (python3.9) total-vm:226036448kB, anon-rss:223754920kB, file-rss:2740kB, shmem-rss:0kB, UID:1000 pgtables:438192kB oom_score_adj:0
[144772.062513] oom_reaper: reaped process 1699 (python3.9), now anon-rss:0kB, file-rss:0kB, shmem-rss:0kB
In this case, my partial sum is of order 11, and it would be very valuable for me to achieve this solution, but I don’t need to go past 11.
Is this to be expected? Since the derivative of arcsin is 1/sqrt(1-x**2), I expect the calculations to be no worse than what I have done in the past.
Is there any reason why series(asin(α)) would be expected to fail?
Should I try anything else?
Tom
(Dr. Thomas S. Ligon)
Frohnloher Str. 6a
81475 Muenchen
Germany
Tel. +49(89)74575075
Jatin Bhardwaj
Hi Tom,
Instead of series(asin(α)), you can use SymPy’s ring_series.rs_asin, which is optimized for power series manipulations and significantly reduces memory usage.
```
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, x = ring('x', domain='QQ')
asin_series = rs_asin(x, x, 11)
print(asin_series)
```
This method is much more efficient and avoids the heavy symbolic differentiation that series() performs.
You can also check out the documentation for the rs_series in `ring series` module.
Jatin Bhardwaj
thomas...@gmail.com
Hi Jatin,
thanks, this looks very promising. I haven’t used polynomials much yet, since simple symbolic expressions do everything I need, and it looks like I have a deficit there. Here is my current code:
from sympy.polys.domains import QQ
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, m = ring('HillSeries2.m', domain='QQ')
α = rs_asin(α, HillSeries2.m, 3) # should be 11
Exception in the last line:
'Add' object has no attribute 'ring'
It looks like my expression needs to be a polynomial
I added a line
α = α.as_poly(domain='QQ')
but that didn’t help.
In the sample code
R, m = ring('HillSeries2.m', domain='QQ')
is not documented very well. It looks like it is defining a ring R, but not using it, but my IDE doesn’t give me a warning for that.
I am using documentation
https://docs.sympy.org/latest/modules/polys/ringseries.html#rs-series
Von: sy...@googlegroups.com <sy...@googlegroups.com> Im Auftrag von Jatin Bhardwaj
Gesendet: Montag, 24. Februar 2025 22:08
An: sympy <sy...@googlegroups.com>
Betreff: [sympy] Re: sympy process out of memory calculating series of asin
Hi Tom,
Instead of series(asin(α)), you can use SymPy’s ring_series.rs_asin, which is optimized for power series manipulations and significantly reduces memory usage.
```
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, x = ring('x', domain='QQ')
asin_series = rs_asin(x, 11)
print(asin_series)
```
This method is much more efficient and avoids the heavy symbolic differentiation that series() performs.
You can also check out the documentation for the rs_series in `ring series` module.
Jatin Bhardwaj
On Monday, 24 February 2025 at 13:52:57 UTC+5:30 thomas...@gmail.com wrote:
--
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.
To view this discussion visit https://groups.google.com/d/msgid/sympy/86f0738d-851f-48e3-9793-f1709a9f1c24n%40googlegroups.com.
Jatin Bhardwaj
Hi Tom,
I noticed a small issue in your implementation. I’m not sure if this is a standalone piece of code or part of a larger codebase, but I’d be happy to help clarify.
The ring function returns a tuple containing the ring itself and all variables defined within it. In your code:
``` R, m = ring('HillSeries2.m', domain=QQ) ```
You’re defining the variable 'HillSeries2.m' and storing it in m. However, using dots in variable names can sometimes cause confusion. It’s generally a good practice to use simple, consistent names for better readability.
1. The first argument of rs_asin should be a function that is either a polynomial or a series, and its variable must be part of the ring. However, α is not defined before this line.
2. The second argument should be the ring variable, but HillSeries2.m is not a separate variable. Instead, you should use m, since that’s where the variable was stored when defining the ring.
Take a look at this code. I hope this helps,
```from sympy.polys.ring_series import rs_asin
R, α = ring('α', domain=QQ) # Define the ring with variable α
s = rs_asin(α, α, 8) # Compute the series expansion of asin(α) up to order 8
thomas...@gmail.com
Hi Jatin,
thanks again for your valuable help. Here is a long explanation for something that should be very simple but might be causing confusion.
I have a total of about 10,000 lines of code dealing with power series in powers of m. m is defined as
m = symbols('m')
Now I have multiple modules and all of them use m, so I want to make sure that they all use the same thing. I could have repeated the definition m = symbols('m') in each module but was not sure that would be using the same thing, so m is defined in HillSeries2.py and when I use it in another module, I refer to it as HillSeries2.m.
Today, I found documentation of the ring function in
https://docs.sympy.org/latest/modules/polys/domainsref.html#sympy.polys.rings.ring
This documentation, like a lot of SymPy documentation, defines the INPUT parameters, but provides no complete documentation of the OUTPUT parameters. For ring, we have:
Construct a polynomial ring returning (ring, x_1, ..., x_n).
But what is the format and meaning of the elements of the output parameter?
It looks like the symbols are defined in the input and repeated in the output, which seems to be a bit redundant, but there is probably a good reason for that.
You write “However, α is not defined before this line.” Well, of course it is defined, I just didn’t copy everything into this discussion. In fact, it is part of the trigonometric formula for solving a cubic equation. If you read that in Wikipedia, you will find the original cubic equation with coefficients a, b, and c, and the reduced equation with coefficients p and q, then Cardano’s formula for the solutions, using p and q. After that, the trigonometric solution is presented, using p, q, cos, and acos. However, I am not using Wikipedia, but a different paper, and for that I need asin instead of acos.
α is a function of p and q, which are functions of a, b, and c. In my code, all are symbolic expressions of power series in m with rational coefficients. I would be happy to copy the code in here, but that looks to me more like a distraction than help. Let me know if you want to see that part.
α is a symbolic expression of a partial sum of a power series in m.
I print the current value of α in Latex format, but I call it ‘arg’, because it is not the final α, but the argument of asin.
The code is
print(latex(Eq(symbols('arg'), α), order='rev-lex'))
and the Latex output is
arg = -1 + \frac{867 m^{2}}{128} - \frac{12393 m^{3}}{512} - \frac{649831 m^{4}}{32768} - \frac{124581731 m^{5}}{589824} + \frac{14200641281 m^{6}}{113246208} - \frac{25578912841 m^{7}}{50331648} + \frac{241589534597669 m^{8}}{173946175488} - \frac{2019731412177637 m^{9}}{521838526464} + \frac{2581283623490644391 m^{10}}{601157982486528} - \frac{806384032365914837207 m^{11}}{36069478949191680}
Now I have activated the line that makes the symbolic expression α into a polynomial. With that, I am getting exception
'Poly' object has no attribute 'ring'
instead of
'Add' object has no attribute 'ring'
My IDE shows me that α is now of type Poly and has a value of
Poly(-806384032365914837207/36069478949191680*m**11 + 2581283623490644391/601157982486528*m**10 - 2019731412177637/521838526464*m**9 + 241589534597669/173946175488*m**8 - 25578912841/50331648*m**7 + 14200641281/113246208*m**6 - 124581731/589824*m**5 - 649831/32768*m**4 - 12393/512*m**3 + 867/128*m**2 - 1, m, domain='QQ')
The code is now
from sympy.polys.domains import QQ
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, HillSeries2.m = ring('HillSeries2.m', domain='QQ')
# α = Poly(α)
α = α.as_poly(domain='QQ')
α = rs_asin(α, HillSeries2.m, 3) # should be 11
Tom
To view this discussion visit https://groups.google.com/d/msgid/sympy/07031884-de5c-4f77-9ba1-f31987a6bcc6n%40googlegroups.com.
thomas...@gmail.com
Hi Jatin,
I would like to give you a quick update. I got help from the GitHub Copilot.
I decided to try using m = symbols('m') in both modules, thus eliminating HillSeries2.m (replacing it by m).
My code now looks like this:
from sympy.polys.domains import QQ
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, m = ring('m', domain='QQ')
# α = Poly(α)
αP = α.as_poly(domain='QQ')
# α = R(α)
# αP = Poly(α, m)
α0 = αP.coeff_monomial(1) # get constant term
αR = αP - α0 # remove it
αR = R(αR) # define with ring
αR = rs_asin(αR, m, 3) # convert to series, 3 should be 11
α = α0 + αR # restore constant term
print(latex(Eq(symbols('α'), α), order='rev-lex'))
The line αR = R(αR) # define with ring produces an exception:
expected Rational object, got Poly
I’ll search more for that tomorrow.
thomas...@gmail.com
Hi Jatin,
finally, I can give you a status update: success! Thanks again for the help. I also got help from Copilot, and the most important suggestion was to run α = R(α). After that, solving each error led to the next one, and I was sometimes running around in circles. Finally, the code is working well and extremely fast. After series() ran for over 24 hours and failed due to memory exhaustion on a VM with 218GB RAM, rs_asin() completed almost immediately on a laptop with 16GB RAM. Congratulations to everyone who helped make that happen.
My code now looks like this:
αP = Poly(α, m) # convert to polynomial
α0 = αP.coeff_monomial(1) # get constant term
αR = αP - α0 # remove it
from sympy.polys.domains import QQ
from sympy.polys.rings import ring
from sympy.polys.ring_series import rs_asin
R, m = ring('m', domain='QQ') # makes m into PolyElement
# α0 = α0.as_poly(domain='QQ')
αR = αR.as_poly(domain='QQ')
coeffs = [Rational(c) for c in αR.all_coeffs()]
αR = R.from_list(coeffs)
αR = R(αR)
αR = rs_asin(αR, m, 11+1) # convert to series
αR = αR.as_expr()
α = α0 + αR # restore constant term
# α = sympy.core.expr.Expr(α)
print(latex(Eq(symbols('α'), α), order='rev-lex'))
m = symbols('m') # revert to original m
Do you have any recommendations for simplifying or streamlining the code? I have a few remarks:
1) I use m a lot and now define it as m = symbols('m') in every module. However, R, m = ring('m', domain='QQ') redefines m as a PolyElement, which is necessary later in the code. Finaly, I revert to the original definition of m after running this part of code and before continuing with the next parts. I find this quite annoying and hope there is a better way to do it, without redefining m.
2) rs_asin() requires a polynomial without a constant term, so I had to remove it before running rs_asin and then restore it afterwards. It would be nice if this could be done automatically as part of rs_asin(). Would that be a good suggestion for GSoC?
Jatin Bhardwaj
Hi Tom,
First of all, I apologize for not being able to reply to your previous message yesterday—I was a bit busy.
Congratulations on solving the issue! There are a few parts of the code that can be simplified or removed:
• αR = αR.as_poly(domain='QQ') – Unnecessary, since αR is already a Poly object, so there’s no need to re-convert it.
• αR = R(αR) – Redundant, as R.from_list() already creates a ring element.
Regarding your question about redefining m, unfortunately, it must be a ring element to work with ring_series, so there’s no way around that.
Now, coming to your second point—I’ll try to implement this before GSoC. However, if I’m unable to do so, I will surely include it in my GSoC project. I’m planning to work on Power Series domains and possibly fix some issues with ring_series.
Recently, one of my PRs (gh-27669) was merged, which added support for arctan in ring_series, allowing it to handle constant terms properly. I’m working on extending this functionality to other trigonometric functions, and I’ll give you a follow-up once it’s done.
Best regards,
Jatin
thomas...@gmail.com
Hi Jatin,
thanks again for your help. I now have all aspects of using rs_asin working, but there is still an important issue left to solve. I’ll explain everything from the beginning here, but in a more condensed version.
I have an expression, let’s call it P for polynomial. It is a cubic equation in X, and each coefficient is a power series in m, specifically a partial sum up to m**11. I want to solve the cubic equation via the trigonometric method as described in
Fuchs, D. and S. Tabachnikov, Mathematical Omnibus. Thirty Lectures on Classic Mathematics, draft of a monograph.
which is available at
In order to do that, I need to calculate the coefficients of the reduced equation, then I need to calculate an arcsine. After that, I need to calculate a sine.
Since I am working with power series, I tried series(asin(P)), but that failed due to performance and memory issues. Now I am using rs_asin and that runs much, much faster.
Unfortunately, rs_asin does not support calculation of a polynomial with a constant term. Now, it is easy enough to find the constant term and separate it from the rest, giving P=P0+PR, where P0 is the constant term and PR is the rest. I can calculate
rs_asin(PR), which runs very fast and then looks plausible. However, I CANNOT recombine the constant term with the result because there is no sum form for arcsine. The sine of a sum is
sin(α + β) = sin α cos β + cos α sin β
but there is no formula for the arcsine of a sum. As a result, I cannot combine P0 with asin(PR) to get asin(P).
Do you have an idea?
To view this discussion visit https://groups.google.com/d/msgid/sympy/6382da15-4cc4-442c-b880-7b4f01693378n%40googlegroups.com.
Jatin Bhardwaj
Hi Tom,
You’re very welcome! I’m glad to hear that you have rs_asin working well now.
Regarding your issue with recombining the constant term P_0 with rs_asin(P_R), you’re absolutely right that there isn’t a direct sum formula for arcsin like there is for sine. However, a workaround can be derived from the Taylor series expansion of arcsin around P_0.
For small P_R, we can approximate:
arcsin(P_0 + P_R) = arcsin(P_0) + rs_asin( P_R /sqrt( 1 - P_0^2))
Unfortunately, this isn’t directly implemented in rs_asin yet. I’ll try to raise a PR for this soon, but for now, you can use this rough implementation of this approach:
- The EX domain in SymPy (Expression Domain) allows handling of symbolic expressions like asin and sqrt as coefficients within polynomial rings. If you’re unfamiliar with it, I recommend checking out the Expression Domain in SymPy’s documentation.
- This code makes use of multiple functions from ring_series.py in SymPy. I highly recommend checking the documentation for ring_series to better understand the internal workings.
Let me know if this approach works for you!
Best regards,
Jatin
thomas...@gmail.com
Hi Jatinj,
thanks for the feedback. This looks like good progress.
However, I have another solution: I simply use the well-known series expansion of arcsine, so I neither need asin nor rs_asin. It looks like it is working, but I am not yet finished testing.
To view this discussion visit https://groups.google.com/d/msgid/sympy/d64402e1-9333-4246-a844-13ea45407cf9n%40googlegroups.com.