Towards a new SymPy: part 1 - Outline

This is the first part of a series of posts on big changes for SymPy with particular focus on speed. What I am going to describe in these posts are:

  • A description of the of the issues around speed with SymPy currently.

  • Work that I (and others) have done to speed up SymPy after the past few years.

  • Work in particular that I have done over the last year to make SymPy faster.

  • What I think are the next steps to make SymPy faster.

The other posts in this series can be found at Towards a new SymPy.

Over the last year in particular I have been working as a part of a CZI funded project that has three strands. One of those three strands is for me to work on speeding up SymPy. Now that we come to the end of that year, I want to describe what has been done and spell out a vision for the future.

I will be writing this in a series of blog posts. This first post will outline the structure of the foundations of a computer algebra system (CAS) like SymPy, describe some problems SymPy currently has and what can be done to address them. Then subsequent posts will focus in more detail on particular components and the work that has been done and what should be done in the future.

I am writing this with the intention that it should be accessible to someone who is not a SymPy developer although other SymPy developers are the intended audience for many of the points that I will make. Many of the things that I will describe here are not well understood even by many SymPy developers though and a major goal of this series of posts is to help to try to change that.

What is a CAS?

There are many different kinds of computer algebra system (CAS) and they have different intended uses that are driven by their different kinds of authors, contributors and users. Some like SageMath are directed more at pure mathematicians, while others like SymPy are directed more at scientists and engineers.

Most scientific computing is not done with computer algebra systems. Instead numerical libraries (e.g. NumPy, SciPy, etc in the context of Python) are the usual backbone of computing in science and engineering. What those libraries do is to provide fast implementations of numerical computing based on machine precision floating point numbers. When I say machine precision floating point I mean something like a 64-bit binary floating point number like the float type in Python:

>>> from math import pi
>>> pi
3.141592653589793
>>> type(pi)
<class 'float'>

Many things can be done with machine precision floating point numbers, but there are situations where they are not good enough. Programmers learn fairly early on that floating point calculations are not exact e.g.:

>>> (1 / 49) * 49
0.9999999999999999

We can see here the limits of the precision of 64-bit floating point because we should have a result of 1 but in fact we have something that is quite close to 1 but still not exactly 1. Since the answer is quite close most of the time this is usually not a huge problem even if it does cause some confusion for beginners.

However, there are many situations where this is not acceptable:

  • Small errors can accumulate and become large errors.

  • To solve the error accumulation problem some calculations should be done with many more digits of precision than machine precision floating point numbers can allow.

  • Some calculations cannot be done reliably at all unless they are done exactly.

  • Sometimes we want to calculate not in terms of explicit numbers but in terms of symbols to obtain e.g. “analytic” rather than “numerical” results.

Traditionally these other situations are handled by scientists and engineers turning to a CAS and so in the context of scientific computing the usual distinction is between working “numerically” with floating point numbers or otherwise using a CAS to work “symbolically” or “analytically”. There is a lot more to a CAS than this but in simple terms this is the role that SymPy plays within the Python scientific computing ecosystem. Of course there are many people also using SymPy for applications that are more like pure mathematics but most SymPy users and developers are from a scientific background and that has a big influence on the way that SymPy is developed.

The “core” of a CAS

Different people will have different ideas about what the “core” of a CAS like SymPy is. For example users might think of functions like solve or integrate as being the core functions but these are really high level functions and not really the foundation of a CAS. Most SymPy developers would probably think of the sympy.core module as being the core of SymPy but this actually misses a major component of SymPy which is the sympy.polys module. In my mind these three components will be the core of any CAS like SymPy:

  • A robust numerical evaluation subsystem (sympy.core.evalf).

  • A symbolic manipulation subsystem (sympy.core.basic).

  • A computational algebra subsystem (sympy.polys).

I will explain each of these in turn because it is easy to miss the importance of each of them and I think that many SymPy developers do not understand the role each of these components has and how they relate to one another. Of course there are many other components of SymPy that are important but these three are the foundation of almost all of the higher-level things that SymPy does. Importantly when we want to think about the speed of SymPy we need to understand that almost everything depends on the speed of these three components. Higher level algorithms are absolutely important but at least in SymPy those are all built on top of these three pieces as the foundation for everything else.

A robust numerical evaluation subsystem

Above I use the word “numerical” to mean calculations that are done with machine precision floating point numbers. This is the usual intended meaning of the word but this use of the word is tricky because in the context of a CAS we can also do numerical calculations with arbitrary precision approximate or exact numbers and those are also “numerical” calculations. I will define several levels of numerical calculation ranging from machine precision floating point to exact and symbolic calculations. As we go down this list the calculations become slower but more accurate:

  • Machine precision floating point: This is the usual 64-bit binary floating point numbers that are used in most numerical libraries. They are fast and have a wide range of support in hardware and software but are limited in both range and precision. Most scientific Python libraries like NumPy, SciPy etc work (almost) exclusively with machine precision floating point numbers.

  • Arbitrary fixed precision floating point: This means using floating point numbers that can have as many digits of precision as desired. These calculations are much slower than machine precision floating point calculations but by increasing the precision we can make the calculations as accurate as we want. The mpmath library that is used by SymPy is an example of this. With mpmath you can choose to use say 100 digits and then all calculations will be done with 100 digit numbers.

  • “Robust” floating point calculations: These are calculations that are done with arbitrary precision floating point numbers but are done in such a way that the result is guaranteed (perhaps formally or perhaps heuristically) to be correct to within a specified tolerance. Making this work requires using adaptive precision so if the end result should be accurate to 100 digits then intermediate calculations might need to use many more than 100 digits to control the growth of errors. SymPy’s evalf is an example of heuristic robust numerics and the Arb library that I will discuss later is an example of a library that does formal robust numerics.

  • Exact numerical calculations: these are numerical calculations that are done using e.g. exact numbers like rational numbers or more complicated types of exact numbers.

  • Symbolic calculations: these are calculations that are done with symbols rather than just numbers. The result of a symbolic calculation is a symbolic expression that can be manipulated further. This is the kind of calculation that probably most people think of when they think of a CAS.

Let us just quickly demonstrate each of these:

  • Machine precision floating point:

    >>> from math import sqrt
>>> sqrt(2)  # Use machine precision floating point
1.4142135623730951

  • Arbitrary precision floating point:

    >>> from mpmath import mp
>>> mp.dps = 50 # Use 50 digits for the calculation
>>> mp.sqrt(2)
mpf('1.4142135623730950488016887242096980785696718753769')

  • Robust floating point calculations:

    >>> from sympy import sqrt
>>> sqrt(2).evalf(50) # Compute a result that is correct to 50 digits
1.4142135623730950488016887242096980785696718753769

  • Exact numerical calculations:

    >>> from sympy import Rational
>>> Rational(1, 2) + Rational(1, 3) # Compute exactly
5/6

  • Symbolic calculations:

    >>> from sympy import Symbol, sqrt
>>> x = Symbol('x')
>>> sqrt(x)**2   # Compute with symbols
x

The distinction between these different kinds of calculations can be a bit fuzzy but the first point to note is that the vast majority of scientific computing is done with machine precision floating point numbers as mentioned at the top of the list. Everything below that is what a CAS like SymPy is typically used for. It is also possible to do symbolic calculations with e.g. machine precision floating point numbers so perhaps including “symbolic” in this list does not make sense but I think that from a user’s perspective “symbolic” is in some way the next level after exact numerical calculations.

Many SymPy users do not understand the distinctions between these different kinds of numeric calculations but also many SymPy developers do not understand them either. For example I am not sure how many SymPy developers would automatically realise mpmath and evalf are at different levels in this scheme despite the fact that evalf works entirely by using mpmath. As an example to demonstrate the differences we will compute \(\sin(1)^2 + \cos(1)^2 - 1\) using SymPy’s evalf, mpmath and also SymEngine which is a C++ library that recreates the “core” of SymPy in C++:

>>> import sympy as sym
>>> e = sym.cos(1)**2 + sym.sin(1)**2 - 1
>>> e.evalf()
-0.e-124

>>> from mpmath import mp
>>> mp.dps = 100
>>> mp.cos(1)**2 + mp.sin(1)**2 - 1
mpf('0.0')

>>> import symengine as se
>>> (se.cos(1)**2 + se.sin(1)**2 - 1).evalf()
1.11022302462516e-16

The correct answer here is zero because \(\sin^2(x) + \cos^2(x) = 1\) for all \(x\). SymPy’s evalf method returns the strange looking result -0.e-124. This is evalf’s way of saying that the result is smaller than \(10^{-124}\) and is possibly zero but not proven to be zero. The robust numerics in evalf should ensure that we know that the result is very close to zero but can still never prove that it is or is not exactly zero. We do not get an exact result of zero because evalf is careful not to claim that the result is exactly zero if it cannot prove that. Usually evalf will prove that a nonzero expression is nonzero but it can never prove that a zero expression is exactly zero. This is a fundamental limitation of robust (rather than exact) numerics: we can never prove that something is exactly zero. In order to compute this result evalf will have used mpmath but with higher and higher precision until it reached the maximum allowed precision (124 digits) and gave up.

The mpmath result is exactly zero which is exactly correct but it is only exactly correct by “luck”. Here mpmath is using 100 decimal digits but it still does not prove that the result is zero. It does however return something that is indistinguishable from an exact zero unlike evalf which made it clear that it could not prove that the result is exactly zero. This is a fundamental limitation of arbitrary fixed precision numerics as compared to robust numerics: we can use as many digits as we like and the result will likely be accurate but without very careful analysis we generally do not know anything about how accurate the result is. It is always possible that had we used more digits we would have found that the result was not zero after all but mpmath does not indicate any level of uncertainty about the result it returns.

By contrast SymEngine’s evalf has computed the result here using machine precision floating point numbers and so it gives a fast result but the result is not zero and there is no way to know if that is due to a rounding error or not. We can ask SymEngine’s evalf to use more digits but it will still only use arbitrary fixed precision (like mpmath) and not robust numerics (like SymPy’s evalf):

>>> (se.cos(1)**2 + se.sin(1)**2 - 1).evalf(100)
-3.9443045261050590270586428264e-31

There is nothing inherently wrong with SymEngine’s evalf and there are good reasons to use all of the different levels of numeric calculation that I have described above. Many SymPy users would actually be happier if SymPy’s evalf was faster even if less accurate and would prefer the behaviour of SymEngine’s evalf. However these levels of numeric calculation are not interchangeable and SymPy’s other two core systems (the expression system and the computational algebra subsystem) absolutely do depend on SymPy’s evalf giving robust numeric evaluation and not just arbitrary fixed precision numeric evaluation: many things would break if SymPy’s evalf was changed to behave like SymEngine’s evalf.

I will propose later a plan for how to improve SymPy’s numeric evaluation subsystem both in terms of speed and accuracy. For now I just want to note that the original author of both mpmath and SymPy’s evalf is Fredrik Johansson who subsequently went on to create Arb which is a library for doing formal robust numerics. SymPy’s evalf should change to using Arb-like formal robust numerics and should ultimately provide the option to use the Arb library as the basis for this subsystem.

The symbolic expression system

SymPy’s symbolic expression system is the second component of what I call the “core” of SymPy and is located in the sympy.core package. This is what most SymPy users and contributors are used to working with. This system defines expressions in a symbolic tree representation e.g.:

>>> import sympy as sym
>>> x = sym.Symbol("x")
>>> e = x**2 + 2*x + 1
>>> e
x**2 + 2*x + 1
>>> sym.srepr(e)
"Add(Pow(Symbol('x'), Integer(2)), Mul(Integer(2), Symbol('x')), Integer(1))"
>>> sym.print_tree(e, assumptions=False)
Add: x**2 + 2*x + 1
+-One: 1
+-Pow: x**2
| +-Symbol: x
| +-Integer: 2
+-Mul: 2*x
  +-Integer: 2
  +-Symbol: x

Most work that goes on in SymPy is done either on the internals of this expression system or on the other code that operates with these expressions. In many ways this system is nice but there are also many problems with it. Essentially it is designed with an emphasis on what would be nice for end users who are doing simple things and as a result is not well suited as a foundation for building more complex algorithms. What we do not have is any alternative that can be used instead of this system for the internals of SymPy.

Many of the problems with SymPy and in particular its performance stem from overuse of this expression system. Unfortunately the prominent exposure of the symbolic subsystem in the SymPy API makes it very difficult to change and so realistically the best path forward is to reduce the usage of this system at least in the internals of SymPy. That is difficult though because we usually do not have a clear alternative to use instead and most SymPy developers do not understand how to use SymPy apart from by using this system.

When I say that the symbolic subsystem is overused I should be clear about what the alternatives would be:

  • Use a symbolic subsystem that has a very different design.

  • Use the computational algebra subsystem instead.

Much of the work that I have done recently in SymPy has been to try to expand the computational algebra subsystem, to make it faster and to make more use of it for heavier algorithms in the internals of things like solve, integrate etc. Many things in SymPy (e.g. matrices) can almost immediately be made a lot faster simply by having them use the computational algebra subsystem instead of the symbolic subsystem. I will talk more about this later.

Of course being a CAS that is primarly intended for symbolics there are many things in SymPy that do need to use a symbolic subsystem but the current design of core symbolics in SymPy is not suitable for most of those things. The problems manifest both in terms of:

  • speed: many things are much slower than they could be.

  • behaviour: many things being more difficult to do.

  • features: many users want to do things that the design cannot really support.

  • bugs: the system is hard to use robustly and this leads to bugs.

  • maintainability: making changes to any part of the expression system (including just fixing obvious bugs) is difficult because the effects of any change are far reaching and unpredictable.

What about SymEngine?

In terms of speed one approach that has been mooted for making the symbolic subsystem faster is to rewrite the “core” in e.g. C++ and this is essentially what SymEngine is. So one possibility to make SymPy faster would be to use SymEngine instead of SymPy’s symbolic subsystem. However SymEngine is not a drop-in replacement for SymPy’s symbolic subsystem and it could never be made to be so. There are two basic problems with attempting to use SymEngine for the core of SymPy:

  • On the one hand SymEngine is not sufficiently similar to SymPy’s existing symbolic subsystem so simply switching to use it inside SymPy would break all sorts of things.

  • On the other hand SymEngine is too much like SymPy’s symbolic subsystem so using it would not solve many of the problems that SymPy has and in fact would make it much harder to solve those problems.

The part of SymPy that SymEngine aims to emulate is literally made up hundreds of thousands of lines of code and has poorly defined semantics in many cases. It is basically impossible to make any drop-in replacement for this system that would not differ in ways that would break things. At the same time some of the SymEngine interface is deliberately different from SymPy (e.g. evalf above) so that changing SymPy’s behaviour to be like SymEngine would break compatibility for users and downstream projects.

It is also not possible to adapt the behaviour of SymEngine to be more like what SymPy needs because it is not extensible from Python. Previously Julia used both SymPy and SymEngine for symbolics but subsequently they decided to move away from both and create SymbolicUtils.jl and Symbolics.jl to be the foundation for symbolics in Julia. The reasons given at the time that neither SymPy nor SymEngine were suitable were that:

  • SymPy worked for what they needed but was too slow.

  • SymEngine was fast enough but too inflexible and not extensible from Julia.

These same two considerations actually apply to SymPy itself: SymPy’s symbolic subsystem is too slow for SymPy itself and SymEngine is too inflexible for SymPy to use it as a replacement.

There is a third problem which the Julia people seemed to overlook which is that both SymPy’s symbolic subsystem and SymEngine’s reimplementation of it are based on a design that is not suitable for the kinds of things that SymPy needs either for its internals or also for what many users and downstream libraries would like. There are many aspects to this design problem but the most fundamental one is the problem of automatic evaluation. What I mean by automatic evaluation is basically this:

>>> import sympy as sym
>>> sym.cos(sym.pi/4)
sqrt(2)/2

Quite simply you asked to create the expression \(\cos(\pi/4)\) and SymPy instead gave you the expression \(\sqrt{2}/2\). Maybe that is what you wanted but maybe it is not. The problem is that you do not have any real control over this. There is evaluate=False but it does not work in general and it cannot be made to work in general.

In terms of speed the problem with automatic evaluation is that it makes it impossible to control the performance of higher-level algorithms because every time any expression is created a huge amount of computational work is done in order to try to “evaluate” the expression. Every now and again someone will ask a question like “what is the computational complexity of <some SymPy function>” but if this function uses the symbolic subsystem then it is entirely impossible to answer this question: unbounded computation can occur just during the creation of an expression.

We can try to reduce the cost of the computational work during automatic evaluation but actually we really need to be able to reduce it to zero which is what it would be if there were no automatic evaluation. Creating an expression needs to be so cheap that we think of it as free compared to doing any actual computation (many SymPy contributors do think of it as free which leads them to write very inefficient code).

Many other problems with the design of SymPy’s symbolic subsystem could in principle be fixed without most users really noticing the internal changes. Almost all code that uses the symbolic subsystem relies on automatic evaluation though and so simply making a change to disable evaluation would break almost everything.

The reason this makes it difficult for SymPy to use SymEngine is that SymEngine is also based on automatic evaluation and in such a way that is impossible to control from the outside. In the case of SymEngine evaluate=False is not an option and I believe there are even plans for SymEngine’s internal data structures to change so that it would be impossible to implement something like evaluate=False. SymEngine broadly does “less” automatic evaluation than SymPy which is one reason why it is faster but there is still some. Even just swapping the order of the terms here is unacceptable if there is no way to disable it:

>>> import symengine as se
>>> x = se.symbols('x')
>>> se.exp(x) + x
x + exp(x)

A new symbolic subsystem

What SymPy needs is a symbolic subsystem that is not based on automatic evaluation. Of course many users would still want automatic evaluation and there should be a way to provide that. It is easy though to build a system that has automatic evaluation on top of a system that does not whereas it is impossible to build a system that does not have automatic evaluation on top of a system that does. The current design of Basic and Expr in SymPy and all of their hundreds of subclasses builds automatic evaluation into the core data structures of the symbolic subsystem. This is a fundamental design flaw that causes all kinds of problems for speed, behaviour and extensibility. Using SymEngine (in its current form) can make this subsystem faster but it would then make all of the other problems unsolvable.

The solution to this problem is to build a new symbolic subsystem but it really needs to be built from the ground up: there is no viable way to get there through incremental changes to the existing system. This can be done for the internals of SymPy to get a lot of the benefit in terms of speed and behaviour without breaking compatibility for users and downstream projects (which could still use the existing system).

At some point though it would be better to switch the user facing symbolic system to an implementation that would be based internally on that new subsystem. We can try to minimise the noticeable impact of this change but it would definitely be a breaking change. There has been some discussion of a SymPy 2.0 that would make some backwards incompatible changes and this is the number one thing that should be done. In my opinion we are not ready for SymPy 2.0 until we are ready to switch out the internals of the symbolic subsystem. I don’t see any other potential change that is important enough to warrant any major compatibility break.

Of course if SymPy were redesigned to have a hypothetical new symbolic subsystem with a different design then either SymEngine or something similar (let’s call it SymEngineX) could be designed to provide a faster implementation of that subsystem. The new system though should be designed so that it is possible for a faster implementation to be created without all of the problems that would apply to SymPy using SymEngine as it is now. In particular:

  • The parts that would be provided by SymEngineX must be a small part of the code that makes up the whole system (rather than a rewrite of everything).

  • Those parts would need to have well defined semantics such that it is actually possible to make multiple compatible implementations (at least the SymEngineX implementation would need to be compatible with the pure Python SymPy implementation).

  • The behaviour of the overall system still needs to be controllable from Python so none of the evaluation etc rules of SymEngineX could be hard-coded (in the way that SymPy and SymEngine both hard-code everything right now).

If a new symbolic subsystem was designed in this way then it would become much easier to make an alternate core that would play the role that SymEngine was intended to play in relation to SymPy. It could be SymEngine itself that still provides that core but the parts that SymPy would need from it would look very different to what SymEngine provides right now.

I want to be clear here that I do not want to criticise the excellent work done by many contributors to both SymPy and SymEngine working on these symbolic systems. There is a lot of great code in both projects that make up these systems but there are also problems that are just not fixable without a complete rebuild of the internals starting from the foundations. The existing code and contributions would not be wasted but much of that code would need to be adapted over time to a new framework.

I will talk about my proposals for the symbolic subsystem more concretely in a separate blog post.

The computational algebra subsystem

SymPy’s computational algebra subsystem is the third component of what I call the “core” of SymPy and is located in the sympy.polys package. This subsystem is underused and underappreciated but is absolutely critical to much SymPy’s ability to do anything nontrivial. In a CAS that was aimed more at pure mathematics this subsystem would be the most prominent feature but in SymPy it is mostly hidden away and not well understood. Few SymPy users or developers know how to use the lower levels of this system and few are aware of the features that it provides or why you might want to use them instead of the “higher level” symbolic subsystem.

To give a clear example of how this subsystem is used I will show how to convert an expression from the symbolic subsystem into the computational algebra subsystem using the high-level Poly representation:

>>> import sympy as sym
>>> x = sym.Symbol("x")
>>> e = (x + 1)**2
>>> p = sym.Poly(e, x)
>>> p
Poly(x**2 + 2*x + 1, x, domain='ZZ')
>>> p.as_expr()
x**2 + 2*x + 1

The Poly representation is the highest level representation in the computational algebra subsystem but its internals are still very different from the symbolic subsystem. We can see this by looking at the internal representation:

>>> p.rep.rep
[1, 2, 1]

This is a representation of the polynomial as a list of coefficients. Within the computational algebra subsystem this is referred to as the dense univariate polynomial (DUP) representation. In this representation we cannot represent the unexpanded power (x + 1)**2 but instead we have to expand it out to x**2 + 2*x + 1 because we can only store the coefficients of the expanded form. This is a very simple example but this computational algebra subsystem can represent much more complex expressions. I have described this subsystem in more detail in this page of the SymPy documentation:

https://docs.sympy.org/latest/modules/polys/domainsintro.html

Many important SymPy functions like factor, cancel etc work by converting from the symbolic subsystem to the computational algebra subsystem, doing some computation and then converting back to the symbolic subsystem. These conversions back and forth between the two systems are quite expensive and often just converting back to the symbolic subsystem can be more expensive than the actual calculation that is done in the computational algebra system. This is because just creating expressions is slow in the symbolic system (mainly because of automatic evaluation). It is much more efficient to perform an entire calculation in the computational algebra subsystem without converting back and forth but any code that does this needs to be written by someone who has some understanding of how to use the computational algebra subsystem. At least with the conversions someone can write code that seems like it only uses the symbolic subsystem without needing to know anything about the computational algebra subsystem. The price for that convenience is that it makes various things in SymPy slower than they need to be.

So many things could be made faster if SymPy just used the computational algebra subsystem more. The problem is that a certain mathematical background is needed just to understand how to use it and most of the people who use and contribute to SymPy do not have that background. This is also not helped by the fact that the system is not well documented: there is a lot of documentation but mostly it has not been written so that someone who does not have a background in computational algebra could understand even what the basic functions do (I wrote the doc page linked above to try to address this partially).

Speaking of documentation I think it is quite telling about both SymPy users and developers that while here I am describing the sympy.polys module as one of the three pillars in the foundation of all of SymPy the docs do not even bother to mention it at top level (as of SymPy 1.12):

https://docs.sympy.org/latest/reference/index.html#reference

The main front-page of the API docs consider that “polynomials” should just be listed under “topics” unlike ntheory (number theory) which is a whole section to itself. The tutorial does not mention the polys module directly, there is nothing in the “guides” section about it etc. To give some context here is a very rough count of how many lines of code are in each of SymPy’s top-level modules (this includes tests, comments, docstrings etc):

polys            100651
physics           80805
core              58583
printing          48638
solvers           43966
functions         42605
matrices          32942
utilities         32248
combinatorics     26144
integrals         25588
parsing           25213
stats             22356
tensor            19663
simplify          18296
geometry          15534
assumptions       11524
ntheory           11464
series            11258
plotting          11034
sets              10400
concrete           7093
logic              6929
vector             6772
codegen            6489
categories         4730
testing            4610
holonomic          4299
crypto             3958
calculus           3665
diffgeom           3046
liealgebras        2189
external           2182
algebras           2086
discrete           1813
interactive        1416
strategies         1406
multipledispatch   1239

The other things listed as “topics” alongside polys are geometry, holonomic, liealgebras, categories, crypto, diffgeom, plotting, and stats. Some of these are barely used by anyone while others like plotting are widely used and should also not be listed here. The “polynomials” module is the largest module in SymPy making up over 10% of the all of the code and is used by everything else but is sorely undervalued by the SymPy community.

The basic problem here is that SymPy is a project that is mostly used and developed by people who are not pure mathematicians and the polys module does not look like something that would be relevant to them. We need to understand though that a CAS needs to have a computational algebra subsystem in order to do the things that scientists and engineers want it to do as well. Higher level functions that like solve, integrate etc that are used by scientists and engineers need to be built on top of a computational algebra subsystem in order to be fast and reliable.

Maintenance of the computational algebra subsystem

When I started working on SymPy I did not have a background in computational algebra and I also did not understand how to use this subsystem (the polys module). Over time I have learned more about it and I have also learned how to use it and I have made various improvements to it. To reach this point I had to read a lot of textbooks and papers and look at what other (more algebra-oriented) CASs do.

To begin with what I found was that there was this large almost unmaintained subsystem that was critical to SymPy but yet did not have any active maintainer working on it. Kalevi Suominen was the only SymPy developer who seemed to really understand this part of the codebase and was critical to ensuring that at least pull requests could be reviewed but was also not actively working on major improvements or any redesign. Kalevi’s guidance made it possible for me to make improvements to this subsystem while in the process of learning about it.

The code in the polys module is well designed. Of course there can always be improvements but the fundamentals of the design are good. Looking back to before my own involvement in SymPy it looks like it was originally developed by Mateusz Paprocki and others. It is clear from the code that it was at some point well maintained and loved by people who understood it well and cared about it and had some vision for how it should work and what the plan was for the future. Then it looks though as if those people at some point just moved on leaving the code with various things in a sort of work-in-progress state or in the middle of a redesign that was never finished. This is probably the hardest part about trying to work on it because there are lots of pieces there but no one remembers what the plan was supposed to be for the future.

More recently though there has been renewed interest in this subsystem and there are several maintainers who can understand the code and make improvements. The most significant improvements recently are

  • The addition of the DomainMatrix class which has been mostly work done by myself (I will talk about this more later).

  • The significant redesign and expansion of the sympy.polys.numberfields module by Steve Keiffer who has also added some significant algorithms and improved others in parts of the code that had been untouched for years.

There is another contribution that I want to draw attention to that is very relevant which is this PR by another maintainer Sangyub Lee:

The significance of this PR is that it transfers the calculation of a particular operation (here the Matrix.eigenvects method) from the symbolic subsystem to the computational algebra subsystem. This is a very good example of exactly what needs to be done more in SymPy and how things that many SymPy users could relate to (eigenvectors in this example) can be computed much faster and better simply by using the computational algebra subsystem more.

The main point here though is just that the polys module is not (almost-)unmaintained any more. There are multiple maintainers who can understand it and make improvements to it and we are now in a better position to improve its capabilities and expand its use within SymPy.

Improvements to the computational algebra subsystem

Above I casually suggested that the numeric subsystem should be completely rewritten to use Arb-like semantics. I then went on to say that the symbolic subsystem should be rewritten (which is a major undertaking for SymPy!). I would not say the same thing about the computational algebra subsystem:

  • The fundamentals of the design are good.

  • Much of the code is written exactly how it should be written and is well micro-optimized.

  • The more widely used parts are mature and well tested.

Of course there can be improvements like adding more features, improving algorithms and so on, but this system can be improved incrementally. There is however a major limitation here. While I say that much of the code is written exactly how it should be written, it is all Python code and the lower level parts of this subsystem are precisely the sort of thing that should be written in a lower level language like C. It is clear from the code that this was always part of the development plan but then it just didn’t happen. Happily we do not need to rewrite everything in C ourselves because there is now already a C library that does exactly what we need: Flint. I have been working on making it possible for SymPy to leverage Flint which is a library for doing fast computations with polynomials, matrices and other things in the same way that SymPy’s computational algebra subsystem does.

As for the other topics I mentioned above, I will talk about this more in a subsequent post.

Summary

Above I have described what I consider to be the subsystems that make up the core of SymPy and what the current state of each of them is, as well as what I think needs to be done to improve them. In subsequent posts I will talk about each of these in more detail to explain what has been done and what should be done going forward. Briefly though:

  • The numeric subsystem should be rewritten to use Arb-like semantics and should be able to use Arb directly at least as an optional backend.

  • The symbolic subsystem should be rebuilt from the ground up based on a non-evaluating representation of expressions.

  • The computational algebra subsystem should be made faster both by improving its algorithms and also by leveraging Flint.

  • The computational algebra subsystem should be used more by the rest of SymPy.