[Idea Prompts] Adding More Matrix Decomposition Methods to SymPy for GSOC 2025
Temiloluwa
Hello SymPy Community,
My name is Temiloluwa, and I would like to propose the addition of more matrix decomposition methods (Schur, Polar, and Hermite Decomposition) as my project idea for Google Summer of Code (GSoC) 2025.
I have successfully gotten a PR in for a QR decomposition method for DomainMatrix and I am currently finalizing a PLDU decomposition and a fraction-free QR decomposition (QRD) pull request for DomainMatrix, aimed at improving the GramSchmidt process
From my exploration of the SymPy codebase, I have observed that the few matrix decomposition methods that exist are LU, QR, Cholesky variants and some others to name a few. However, important methods like:
- Schur Decomposition ()
- Polar Decomposition ()
- Hermite Decomposition ()
are yet to be implemented.
Integrating these decompositions would enhance SymPy’s linear algebra capabilities, particularly in eigenvalue computation and numerical stability. I look forward to discussing this further and receiving feedback from the community.
Best regards,
Temiloluwa
Oscar Benjamin
Yes, these kinds of things would be good to have in SymPy. More
important than adding many different factorisations though is having
good implementations of the most important ones. I would say that just
implementing fraction-free and division-based LU algorithms fully
would make a suitable GSOC project. I suspect that you currently
underestimate how much is involved in that even though you have a PR
that implements the algorithm already. The additional things that are
needed to complete the LU implementation are things like:
- Both sparse and dense implementations of LU and FFLU for
DomainMatrix (consider if any very different sparse algorithms are
worth implementing).
- Make use of python-flint for FLINT's FFLU function when possible and
ensure that the outputs are comparable to the SymPy implementation.
- Numerically stable version of LU for floats.
- Matrix methods that call into the DomainMatrix methods making them
accessible to users and speeding up end user calculations.
- Upper and lower triangular solves for LU and fraction-free LU.
- Solving over/under determined systems.
- Make use of LU to compute other things like matrix inverse, RREF, nullspace.
- Lots of timings and benchmarks to compare performance of different
approaches for different domains, shapes, densities etc. Decide when
different algorithms should be used.
- Make it so that solve and Matrix.solve and other functions use
LU-based solve if it is better.
See also this PR which I did not finish but was intended to make all
linear equation solving code go through a single place so that it
could be optimised for all cases:
https://github.com/sympy/sympy/pull/26000
Putting all of that sort of stuff together is more useful than adding
other less frequently used factorisations. A good implementation of LU
is a foundation for many other things so it is worth spending some
time to make it as good as possible. As you have seen other things
like QRD can be almost trivial to implement if they can leverage the
LU decomposition.
In general not all matrix factorisations that can be defined
abstractly are actually suitable for exact or symbolic arithmetic
rather than floating point arithmetic and vice versa. For example I am
not sure that I have seen any implementation of the polar
decomposition for exact numbers. Unless SymPy could provide something
that is more useful than SciPy's polar function I'm not sure it is
worth adding. Note that a high precision implementation should be part
of mpmath rather than SymPy. The only value implementing this in SymPy
would be if it could be exact but I don't know how you would compute
the polar decomposition exactly or whether that is something that only
makes sense for very small matrices.
> 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/008ee22b-e3c3-40b0-8565-18e2f1aef17cn%40googlegroups.com.
Temiloluwa
Hi Oscar,
Thank you for your detailed feedback! I now see that rather than adding multiple matrix decompositions, a well-optimized LU and Fraction-Free LU (FFLU) implementation is far more impactful in sympy. Now that I have gotten a suitable topic that is GSOC worthy and it is worth spending time on during the summer. Here is a draft template of my idea, I await your feedback for refinements.
Title
implementing division based LU & Fraction-Free LU (FFLU) algorithms in SymPy
Idea
LU decomposition is a fundamental operation in linear algebra, used for solving systems of linear equations, computing matrix inverses, and finding determinant values. Currently, SymPy lacks a fully optimized and integrated LU decomposition, especially for fraction-free computations.
This project aims to:
Implement both sparse and dense LU/FFLU in DomainMatrix, ensuring efficient computations for various matrix types.
Make use of python-flint for FLINT's FFLU function when possible and ensure that the outputs are comparable to the SymPy implementation.
Implement numerically stable LU decomposition version for floats
Improve Matrix methods that call into the DomainMatrix methods, making LU/FFLU easily accessible to users and improving overall computational speed.
Develop upper/lower triangular solvers for LU/FFLU.
Solve over/under determined systems
Integrate LU-based solving into core SymPy functions (solve, Matrix.solve, linsolve, _linsolve_aug) for better performance.
Utilize LU for computing matrix inverse, RREF, and null space efficiently.
Add Extensive benchmarking & lots of timing to compare performance of different approaches for different domains, shapes, densities etc and will Decide when different algorithms should be used.
By completing these improvements, LU/FFLU will become the backbone of SymPy’s equation-solving infrastructure, making linear algebra operations significantly faster and more robust.
Status
SymPy currently has basic LU decomposition, but it lacks an optimized fraction-free version and efficient sparse implementations.
Work on QRD using FFLU is in progress (see PR #27423)
(PR #26000) is refactoring SymPy’s linear equation solvers. This project will extend and integrate LU/FFLU into that.
Python-Flint provides FFLU for integer matrices, and this project will ensure SymPy uses FLINT’s version when applicable.
Involved Software
SymPy
python-flint
Difficulty
Advanced
Matrix decomposition algorithms (LU, FFLU, QRD).
Symbolic and exact arithmetic (ZZ, QQ, python-flint).
Sparse & dense matrix operations in SymPy’s DomainMatrix.
Integration with existing linear solvers in linsolve, _linsolve_aug, and Matrix.solve.
Optimizing performance for exact arithmetic (ZZ, QQ) vs floating-point numbers (RR, QQ) is non-trivial.
Prerequisite Knowledge
knowledge of linear algebra, especially LU decomposition and fraction-free algorithms.
Needs understanding of sparse vs. dense matrix operations.
Project Length
175/350hours
To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxTQX8t9LB-cT5RfOv-t1vDZx%2B9PQD55bTht-gQr33XN8w%40mail.gmail.com.
Temiloluwa
If you have any suggestions for improvements, I would be happy to refine the proposal further. Here is the link to my proposal. I have granted your access to the document.
Looking forward to your feedback.
Best regards,
Temiloluwa
To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxTQX8t9LB-cT5RfOv-t1vDZx%2B9PQD55bTht-gQr33XN8w%40mail.gmail.com.