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/CAHVvXxQV_DycJk%2BqzqkQKcP78wpEXsd_%3DO4D2M5Ovb7XOfTz9g%40mail.gmail.com.
