feat(matrices): Add fraction free solvers for DomainMatrix #25346
Conversation
|
✅ Hi, I am the SymPy bot (v170). I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like:
This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13. Click here to see the pull request description that was parsed.
Update The release notes on the wiki have been updated. |
|
Benchmark results from GitHub Actions Lower numbers are good, higher numbers are bad. A ratio less than 1 Significantly changed benchmark results (PR vs master) Significantly changed benchmark results (master vs previous release) before after ratio
[8059df73] [46558492]
<sympy-1.12^0>
- 84.6±0.9ms 55.3±0.3ms 0.65 integrate.TimeIntegrationRisch02.time_doit(10)
- 84.9±1ms 53.8±0.2ms 0.63 integrate.TimeIntegrationRisch02.time_doit_risch(10)
Full benchmark results can be found as artifacts in GitHub Actions |
|
I think that it's mostly ready, but the only problems are that it is polyfilled for sparse matrix, |
sympy/polys/matrices/ddm.py
Outdated
| def rref_den(a): | ||
| """Reduced-row echelon form of a and list of pivots""" | ||
| b = a.copy() |
There was a problem hiding this comment.
| def rref_den(a): | |
| """Reduced-row echelon form of a and list of pivots""" | |
| b = a.copy() | |
| def rref_den(a): | |
| """Return ``(ra*d, d, p)`` where ``ra`` is the reduced-row echelon form of $a$ and ``p`` is a list of pivots""" | |
| b = a.copy() |
There was a problem hiding this comment.
Or something like that. The denominator being returned should be mentioned.
I would expect to expose this for Matrix by having some logic around exactly what domain or format should be used. Sparse matrices are mostly only useful in the case of simple domains where the elementary ring arithmetic is relatively cheap. Here that means In [455]: %time ok = M.to_DM().to_field().rref()
CPU times: user 1.65 s, sys: 0 ns, total: 1.65 s
Wall time: 1.65 s
In [456]: M = randMatrix(100)
In [457]: %time ok = M.to_DM().to_field().rref()
CPU times: user 1.71 s, sys: 0 ns, total: 1.71 s
Wall time: 1.7 sWhen the matrix is sparse there can be benefits to using the sparse representation but the benefits are usually not huge: In [458]: M = randMatrix(100, percent=1)
In [459]: %time ok = M.to_DM().to_field().rref()
CPU times: user 2.83 ms, sys: 0 ns, total: 2.83 ms
Wall time: 2.85 ms
In [460]: %time ok = M.to_DM().to_field().to_dense().rref()
CPU times: user 9.11 ms, sys: 49 µs, total: 9.16 ms
Wall time: 9.15 msWhen you have a more complicated domain like In [468]: M = eye(10)*(x+1)**10
In [469]: %time ok = M.to_DM().det()
CPU times: user 209 ms, sys: 0 ns, total: 209 ms
Wall time: 207 ms
In [470]: %time ok = M.to_DM().to_dense().det()
CPU times: user 203 ms, sys: 0 ns, total: 203 ms
Wall time: 201 ms
In [471]: %time ok = M.to_DM().to_field().rref()
CPU times: user 24.8 ms, sys: 0 ns, total: 24.8 ms
Wall time: 23.2 ms
In [472]: %time ok = M.to_DM().to_field().to_dense().rref()
CPU times: user 22.7 ms, sys: 0 ns, total: 22.7 ms
Wall time: 22.1 msIn this example the nonzero elements are so expensive and the zero elements are so cheap that the benefit of using sparse routines is gone. Within SymPy there are a mixture of cases for these matrix operations but usually when things are really slow it is because of expression blowup rather than anything that sparse matrix techniques would help with. The example in gh-25009 is a good one.: In [486]: p, q, r, s, t, u, v, w, x, y = symbols('p, q, r, s, t, u, v, w, x, y', positive=True)
...:
...: X = Matrix([
...: [1/y + 1/w + 1/t + 1/s + 1/p, 0, -1/w, 0, 0, -1/t, 1],
...: [0, 1/u + 1/q, -1/q, 0, 0, 0, 0],
...: [-1/w, -1/q, 1/w + 1/q + 1/x, -1/x, 0, 0, -1],
...: [0, 0, -1/x, 1/r + 1/x, -1/r, 0, 0],
...: [0, 0, 0, -1/r, 1/v + 1/r, -1/v, 0],
...: [-1/t, 0, 0, 0, -1/v, 1/v + 1/t, 0],
...: [1, 0, -1, 0, 0, 0, 0]
...: ])
In [487]: print(repr(X))
Matrix([
[1/y + 1/w + 1/t + 1/s + 1/p, 0, -1/w, 0, 0, -1/t, 1],
[ 0, 1/u + 1/q, -1/q, 0, 0, 0, 0],
[ -1/w, -1/q, 1/x + 1/w + 1/q, -1/x, 0, 0, -1],
[ 0, 0, -1/x, 1/x + 1/r, -1/r, 0, 0],
[ 0, 0, 0, -1/r, 1/v + 1/r, -1/v, 0],
[ -1/t, 0, 0, 0, -1/v, 1/v + 1/t, 0],
[ 1, 0, -1, 0, 0, 0, 0]])This matrix is 7x7 but computing the inverse is extremely slow because it requires very large expressions. We are almost at the point of being to make this happen easily for users: In [497]: %time denom, Xnum = X.to_DM().clear_denoms(convert=True)
CPU times: user 34.8 ms, sys: 0 ns, total: 34.8 ms
Wall time: 33.8 ms
In [498]: %time Xinv, denom2 = Xnum.inv_den()
CPU times: user 601 ms, sys: 0 ns, total: 601 ms
Wall time: 600 ms
In [499]: %time Xinvf = Xinv.to_field()
CPU times: user 2.25 ms, sys: 0 ns, total: 2.25 ms
Wall time: 2.27 ms
In [500]: %time Xinv = Xinvf * (Xinvf.domain.convert(denom.element) / Xinvf.domain.convert(denom2))
CPU times: user 216 ms, sys: 0 ns, total: 216 ms
Wall time: 215 ms
In [502]: %time (Xinv * X.to_DM()).to_Matrix()
CPU times: user 817 ms, sys: 0 ns, total: 817 ms
Wall time: 815 ms
Out[502]:
⎡1 0 0 0 0 0 0⎤
⎢ ⎥
⎢0 1 0 0 0 0 0⎥
⎢ ⎥
⎢0 0 1 0 0 0 0⎥
⎢ ⎥
⎢0 0 0 1 0 0 0⎥
⎢ ⎥
⎢0 0 0 0 1 0 0⎥
⎢ ⎥
⎢0 0 0 0 0 1 0⎥
⎢ ⎥
⎣0 0 0 0 0 0 1⎦I think that it is this kind of case: small matrices but complicated expressions that is mostly significant for SymPy users. There are certainly applications for sparse matrices but it is mostly limited to either extremely sparse matrices or to simple domains. Division free approaches are mostly significant where they allow preventing growth of the denominator expressions. |
🟠Hi, I am the SymPy bot (v170). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it. This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75. The following commits add new files:
If these files were added/deleted on purpose, you can ignore this message. |
|
The pyodide tests are failing, specifically job 3: That has happened repeatedly CC @brocksam @eagleoflqj |
|
This is ready for review now. |
|
I've added release notes. |
This seems to have passed in the latest run so maybe it is an intermittent problem. The same job failed more than once before. |
| >>> from sympy.polys.matrices.ddm import ddm_iinv, ddm_imatmul | ||
| >>> from sympy import QQ | ||
| >>> a = [[QQ(1), QQ(2)], [QQ(3), QQ(4)]] | ||
| >>> ainv = [[None, None], [None, None]] |
There was a problem hiding this comment.
| >>> ainv = [[None, None], [None, None]] | |
| >>> ainv = [] |
There was a problem hiding this comment.
It is not really intended that that would be able to work...
There was a problem hiding this comment.
In retrospect it is probably a mistake that these operations work in place. You can see in this example that it does not give any performance benefit but it does complicate the interface.
I would rather not change that right now though.
There was a problem hiding this comment.
I don't understand how it complicates the interface. You need somewhere to put the result if you aren't going to return it. I am only suggesting that you change the docstring to show that you don't have to worry about what shape or what values must be in the result (mutable list) that you send: just send an empty list and it will be filled.
There was a problem hiding this comment.
just send an empty list and it will be filled.
I don’t think that there can be a guarantee that that may work in the future though
There was a problem hiding this comment.
Even if we do want to work with list of lists then still it would be better just to have ainv = ddm_inv(a, K) that returns a new list.
I generally advice to do type hinting for these sort if possible, than having to be worrying about if you wrote something that consistently returns or not.
I don’t prefer that kind of pattern, because it doesn’t make it clear the if the input list references the output list or the output is new.
There was a problem hiding this comment.
I'm not going to change this right now but a future PR could add out of place versions for any current in place functions. Possibly the in place versions could be removed or otherwise the docstring could just be changed to call the out of place version of the function.
There was a problem hiding this comment.
Not sure if anyone will read this now that the PR has merged (and I might not even see the answer) but if you pass a list of lists as is done in the docstring, then it is pointless to fill it with anything. Neither the content nor the shape of the contents have any impact on what is returned -- and what is returned is a list of lists. All I was suggesting that be changed was the docstring to show a list instead of a list filled with lists filled with values (none of which matter and should not even be thought of when calling the routine that only returns a list of lists that replaces the contents of what was passed).
There was a problem hiding this comment.
Usually, the things to put on the docstrings should be the best practices, or the things that are invariant to the code
There was a problem hiding this comment.
But what is the best practice when the routine does not look at the contents or shape of what you are passing to it? In one docstring it is filled like this ainv = [[None]*3 for i in range(3)] and in another like ainv = [[Q(1)]*3 for i in range(3)]. None of this matters if the routine modifies like this ainv[:] = [row[:n] for row in smthng]. We are making the user reading the docstring think about the contents of ainv when it doesn't matter.
|
I think I have responded to all comments. I agree with all except #25346 (comment) Thanks for the reviews! |
|
Thanks, I merge this |
|
Thanks all for the reviews! |
References to other Issues or PRs
Partial implementation of the fraction-free methods in gh-21834.
Brief description of what is fixed or changed
Adds various methods like
solve_den,rref_den,inv_den,rref_dentoDomainMatrix. These methods provide fraction-free or division-free algorithms for matrix inverse and linear solves over rings rather than fields.This will speed up various operations including
linsolve,Matrix.solveandMatrix.invif used by those methods in future. Example:The speed different is significant for large matrices or more symbols. Here is the current Matrix method timing:
Here is the same with the PR using DomainMatrix:
Obviously this example scales very badly by any technique but DomainMatrix is orders of magnitude faster.
Other comments
There are a few things to finish and a few tests to add and release notes before this is ready for merge.
Release Notes
DomainMatrix.rref_den: A fraction-free method for computing the reduced row echelon form (RREF) with denominator was added.DomainMatrix.adj_det: A division-free method for computing the adjugate and determinant of a square matrix was added.DomainMatrix.adjugate: A division-free method for computing the adjugate of a square matrix was added.DomainMatrix.inv_den: Fraction-free or division-free method for computing the inverse of a matrix with denominator was added.DomainMatrix.solve_den: Fraction-free or division-free method for solving a matrix equationA*x = bwas addedDomainMatrix.adj_poly_det: A method for computing the polynomialpsuch thatp(A) = adjugate(A)was added.DomainMatrix.eval_poly: A method for evaluating a polynomial function of a matrixp(A)was added.DomainMatrix.eval_poly_mul: A method for effciently evaluating a polynomial matrix productp(A)*bwas added.