Hi, sorry for not replying before I ended up solving my problem using mathematica since I needed to compute this things quite a few times. But if you could perhaps give out some references and some pointers as to how this "gcd should be implemented using *sparse* polynomials" could be implemented, like maybe some paper on the topic. Since I've never contributed to open source and don't really know anything about sympy internals maybe I wouldn't be able to work on it as a as a "GSOC project " but I would love to work on the topic , this issue is the main reason I stoped using sympy lately.
you mentioned that first step to take might be 1. Add code in inv to check for domains that can be handled efficiently by DomainMatrix (e.g. at least ZZ, QQ, ZZ_I, QQ_I) and delegate to the lower-level methods. Any documentation on which kind of domains can be handled efficiently and how such check could be done ? or should it be benchmarked? the _rep works amazingly fast but the expressions are quite complicated to simplify when compared to the Mathematica output, so In my case even though the actual answer is simple I couldn't get it from the _rep inverse, same with the connected components approach El miércoles, 12 de enero de 2022 a las 2:35:35 UTC+1, Oscar escribió: > On Tue, 11 Jan 2022 at 20:48, Gerardo Suarez <[email protected]> wrote: > > > > Sorry, you're right I made a mistake while deriving that equation in > python the actual matrix that I computed using mathematica is > > > > e = Matrix([[-gamma_c*(n_c + 1) - gamma_h*(n_h + 1), 0, 0, 0, 0, > gamma_c*n_c, 0, 0, 0, 0, gamma_h*n_h, 0, 0, 0, 0], [0, -I*epsilon_c - > gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*(n_h + 1), I*g, 0, 0, 0, 0, > 0, 0, 0, 0, gamma_h*n_h, 0, 0, 0], [0, I*g, -I*epsilon_h - gamma_c*(n_c + > 1) - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0, 0, 0, 0, gamma_c*n_c, 0, 0, 0, > 0, 0, 0, 0], [0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 > - gamma_h*(n_h + 1)/2 - I*(epsilon_c + epsilon_h), 0, 0, 0, 0, 0, 0, 0, 0, > 0, 0, 0], [0, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*(n_h > + 1) - I*(epsilon_h - (epsilon_c + epsilon_h)), 0, 0, 0, -I*g, 0, 0, 0, 0, > 0, gamma_h*n_h], [gamma_c*(n_c + 1) - gamma_h*n_h, 0, 0, 0, 0, -gamma_c*n_c > - gamma_h*n_h - gamma_h*(n_h + 1), I*g, 0, 0, -I*g, -gamma_h*n_h, 0, 0, 0, > 0], [0, 0, 0, 0, 0, I*g, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - > gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(-epsilon_c + epsilon_h), 0, 0, 0, > -I*g, 0, 0, 0, 0], [0, 0, gamma_c*(n_c + 1), 0, 0, 0, 0, -I*epsilon_h - > gamma_c*n_c - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0, 0, 0, -I*g, 0, 0, 0], > [0, 0, 0, 0, -I*g, 0, 0, 0, -gamma_c*(n_c + 1) - gamma_h*n_h/2 - > gamma_h*(n_h + 1)/2 - I*(epsilon_c - (epsilon_c + epsilon_h)), 0, 0, 0, 0, > gamma_c*n_c, 0], [0, 0, 0, 0, 0, -I*g, 0, 0, 0, -gamma_c*n_c/2 - > gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(epsilon_c - > epsilon_h), I*g, 0, 0, 0, 0], [-gamma_c*n_c + gamma_h*(n_h + 1), 0, 0, 0, > 0, -gamma_c*n_c, -I*g, 0, 0, I*g, -gamma_c*n_c - gamma_c*(n_c + 1) - > gamma_h*n_h, 0, 0, 0, 0], [0, gamma_h*(n_h + 1), 0, 0, 0, 0, 0, -I*g, 0, 0, > 0, -I*epsilon_c - gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h, 0, 0, > 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + > 1)/2 - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 + I*(epsilon_c + epsilon_h), 0, > 0], [0, 0, 0, 0, 0, 0, 0, 0, gamma_c*(n_c + 1), 0, 0, 0, 0, I*epsilon_h - > gamma_c*n_c - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, I*g], [0, 0, 0, 0, > gamma_h*(n_h + 1), 0, 0, 0, 0, 0, 0, 0, 0, I*g, I*epsilon_c - gamma_c*n_c/2 > - gamma_c*(n_c + 1)/2 - gamma_h*n_h]]) > > > > Using the one I provided in mathematica yields an error right away. > Sorry about it. The _rep approach seems to work nicely and definitely > something I will use a lot from now on. The answer is a little convoluted > and I'm trying to get it simplified just to compare it with mathematica > > To be clear the _rep attribute is private. I showed it because you asked > if there is interest in improving the performance of Matrix.inv and I > wanted to explain the internals and that it is there. If you write code > that uses the _rep attribute then you cannot depend on that code working in > future versions of SymPy. Note that this is a common convention in Python > that names with a leading underscore are private. Nothing prevents you from > using private attributes but you do so on the basis of the "consenting > adults" principle: you were warned. > > The matrix you show now no longer has a column of zeros but it still has > the sort of structure that I mentioned before: > > In [13]: e.print_nonzero() > > [X X X ] > [ XX X ] > [ XX X ] > [ X ] > [ X X X] > [X XX XX ] > [ XX X ] > [ X X X ] > [ X X X ] > [ X XX ] > [X XX XX ] > [ X X X ] > [ X ] > [ X XX] > [ X XX] > > > In [14]: e.connected_components() > Out[14]: [[0, 5, 6, 10, 9], [1, 2, 7, 11], [3], [4, 8, 13, 14], [12]] > > The output of connected components shows subsets of the rows and columns > that are related by entries in the matrix. The different subsets are > effectively independent and represent subproblems that can be solved > independently. Put another way it shows a permutation of the rows and > columns that can bring the matrix into block diagonal form: > > In [16]: p = sum(e.connected_components(), []) > > In [17]: p > Out[17]: [0, 5, 6, 10, 9, 1, 2, 7, 11, 3, 4, 8, 13, 14, 12] > > In [18]: e[p,p].print_nonzero() > [XX X ] > [XXXXX ] > [ XXX ] > [XXXXX ] > [ X XX ] > [ XX X ] > [ XXX ] > [ XXX ] > [ X XX ] > [ X ] > [ XX X ] > [ XXX ] > [ XXX ] > [ X XX ] > [ X] > > What this means is that the inverse can be computed by inverting smaller > matrices and then putting the entries back together again into the larger > matrix. For example you can get the first component with e[cc[0],cc[0]] > invert that and then distribute the elements of the inverse back to the > original rows and columns of the inverse of the full matrix. This is a > significant optimisation that should be implemented (I think it is already > used in some places). > > The existing routines don't find the inverse of the largest 5x5 component > efficiently but it should be doable. The characteristic polynomial is > calculated quickly and from there it's not much work to get the inverse: > > In [90]: from sympy.polys.matrices.domainmatrix import DomainMatrix, > DomainScalar > > In [91]: dM = DomainMatrix.from_Matrix(e[cc[0],cc[0]]) > > In [92]: %time p = dM.charpoly() > CPU times: user 52 ms, sys: 0 ns, total: 52 ms > Wall time: 50.4 ms > > In [93]: detA = DomainScalar(p[5], dM.domain) > > In [94]: %time cofactors = -(p[4]*dM**0 + p[3]*dM**1 + p[2]*dM**2 + > p[1]*dM**3 + p[0]*dM**4) > CPU times: user 264 ms, sys: 4 ms, total: 268 ms > Wall time: 269 ms > > So far so good. It's taken about 0.3 seconds. Now just need to divide the > matrix of cofactors by the determinant and here is where it slows down: > > In [95]: %time Ainv = cofactors / detA > CPU times: user 59.7 s, sys: 0 ns, total: 59.7 s > Wall time: 59.7 s > > The reason this is slow is because of the slow polynomial gcd > implementation that SymPy uses for multivariate polynomials. It is slow > because it uses the dense multivariate polynomial (DMP) representation. > Instead gcd should be implemented using *sparse* polynomials. Hence the > issue I linked earlier: > https://github.com/sympy/sympy/issues/20874 > DMP and sparse representations are explained briefly here: > > https://docs.sympy.org/latest/modules/polys/domainsintro.html#dmp-representation > > Altogether though we should be able to get the inverse here in a few > minutes even with the slow gcd by using the connected components and then > the charpoly-inverse method for small-ish matrices as I showed above but > it's definitely possible to do much better than that. Improving gcd > calculations would be a *major* boost for SymPy. We should list it more > prominently as a GSOC project because it really is a top priority. > > Separately adding flint as an optional dependency would make short work > of many problems like this. It already supports fast implementations for > inverting matrices of sparse polynomials. I'm not sure that it can handle > QQ_I (Gaussian rationals) as coefficients so it's not entirely clear to me > how it could be used in this specific case but it would definitely speed up > other things enormously. > > -- > Oscar > -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2148d866-b724-4084-9db2-92076026f08bn%40googlegroups.com.
