Why is matrix inversion so slow?

[Gerardo Suarez]

Hi, 

So I was working in some problem and I had to switch to mathematica due to the time the inverse of a matrix was taking, in mathematica it was solved in less than a second while the computation in sympy has been running for a day and hasn't finished.

Is there anyway to speed up matrix inversion in sympy?

Given that it is so much faster in mathematica is there any interest in implementing some method (maybe theirs) that yields the inverse of a matrix faster? if so I'd like to do it maybe someone could give me some pointers to get me started 

The matrix  is 15*15 and attached as an image, but any example of how to speed things up will do. Thanks a lot!

[Oscar]

The image is not a useful way to share this. Can you show simple code to make this matrix? 

(i.e. just the code to make the symbols and then something like the repr() of the matrix but make sure it's fully runnable code without any missing bits)

--

Oscar

[Gerardo Suarez]

Sure, sorry  about the image, I was unaware of how easy is to generate python code with the print methods

from sympy import Symbol, Matrix,I
gamma_c = Symbol('gamma_c')
n_c = Symbol('n_c')
gamma_h = Symbol('gamma_h')
n_h = Symbol('n_h')
epsilon_c = Symbol('epsilon_c')
g = Symbol('g')
epsilon_h = Symbol('epsilon_h')
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], [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], [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, 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, 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, 0], [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, 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, 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, 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, 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, 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], [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, 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, 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, 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, 0], [gamma_c*n_c + gamma_h*n_h, 0, 0, 0, 0, gamma_c*n_c + gamma_h*n_h + gamma_h*(n_h + 1), 0, 0, 0, 0, gamma_c*n_c + gamma_c*(n_c + 1) + gamma_h*n_h, 0, 0, 0, 0, 0]])

However, I care more about how to speed up matrix inversion in general, as it was terribly fast in Mathematica and it never finished in sympy. I guess there is a more efficient way than using 

e.inv()

[Oscar Benjamin]

On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <gerardo...@gmail.com> wrote:
>
> Sure, sorry about the image, I was unaware of how easy is to generate python code with the print methods

Thanks for that. I just wanted to test it to see if I could see what
was slow or whether there is already a faster way to do this in SymPy.
I think that the approach that is needed for this sort of thing is
along the lines of what you see here:
https://github.com/sympy/sympy/issues/21834

Basically there already is a faster matrix implementation that is used
internally by each matrix but it's faster routines are not used most
of the time yet. The faster implementation is called DomainMatrix and
every Matrix has an internal _rep attribute that is a DomainMatrix:

In [44]: M = randMatrix(17)

In [45]: %time Mok = M.inv()
CPU times: user 21.9 s, sys: 32 ms, total: 21.9 s
Wall time: 21.9 s

In [46]: %time Mok = M._rep.to_field().inv().to_Matrix()
CPU times: user 20 ms, sys: 0 ns, total: 20 ms
Wall time: 20.2 ms

For this example you see that DomainMatrix is 1000x faster but the
difference grows as the matrices get larger. See here for more about
DomainMatrix and the domain system:
https://docs.sympy.org/latest/modules/polys/domainmatrix.html
https://docs.sympy.org/latest/modules/polys/domainsintro.html

Note that usually the _rep for a Matrix is a DomainMatrix with the
EXRAW domain. For faster results you really need a domain that is not
EX or EXRAW. In your particular case the domain is more complicated so
you need to do something like this to discover the domain
automatically:

In [50]: from sympy.polys.matrices import DomainMatrix

In [51]: dM = DomainMatrix.from_Matrix(e)

In [52]: dM.domain
Out[52]: QQ_I[g,epsilon_c,epsilon_h,gamma_c,gamma_h,n_c,n_h]

Unfortunately in this domain calculations are not faster and that's
due to multivariate gcd being slow because it's based on dense
multivariate polynomials which are very inefficient when you have many
symbols.

The gcd only comes in at the final step so basically what has been
computed is the matrix of cofactors and you're dividing by the
determinant. Each involves large polynomials and it is trying to
cancel out common factors. It could be possible to skip that step and
just return a matrix where every element has a large and complicated
determinant expression in the denominator.

The things that should be done to improve this are:

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.
2. Add the fraction-free methods for DomainMatrix that are referred to
in the issue above.
3. Implement proper sparse polynomial gcd:
https://github.com/sympy/sympy/issues/20874
4. Improve pivoting in the solver methods for inv, det etc. I had some
ideas in this PR (although it is now out of date):
https://github.com/sympy/sympy/pull/20483
5. Delegate the core polynomial manipulation and some matrix
operations to a fast C library like flint:
https://www.flintlib.org/

The other thing that is really needed is good benchmarks to be added
to the benchmarks repo. We need to make sure that we can keep track of
when things are getting faster or slower every time changes are made.
Otherwise hard-won performance gains can disappear because of
seemingly unimportant modifications in the codebase.

There absolutely is interest in working on this. I plan to do some of
these things but probably not for some time. If you want to have a go
at any of this then I'm happy to help.

--
Oscar

> However, I care more about how to speed up matrix inversion in general, as it was terribly fast in Mathematica and it never finished in sympy. I guess there is a more efficient way than using
>
> e.inv()
> El martes, 11 de enero de 2022 a las 3:21:33 UTC+1, Oscar escribió:
>>
>> On Monday, 10 January 2022 at 15:51:49 UTC gerardo...@gmail.com wrote:
>>>
>>> Hi,
>>>
>>> So I was working in some problem and I had to switch to mathematica due to the time the inverse of a matrix was taking, in mathematica it was solved in less than a second while the computation in sympy has been running for a day and hasn't finished.
>>>
>>> Is there anyway to speed up matrix inversion in sympy?
>>>
>>> Given that it is so much faster in mathematica is there any interest in implementing some method (maybe theirs) that yields the inverse of a matrix faster? if so I'd like to do it maybe someone could give me some pointers to get me started
>>>
>>> The matrix is 15*15 and attached as an image, but any example of how to speed things up will do. Thanks a lot!
>>
>>
>> The image is not a useful way to share this. Can you show simple code to make this matrix?
>>
>> (i.e. just the code to make the symbols and then something like the repr() of the matrix but make sure it's fully runnable code without any missing bits)
>>
>> --
>> 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 sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/5fb279aa-81e9-49ca-bbed-0a65e184c507n%40googlegroups.com.

[Hanspeter Schmid]

Hmm, the last column of your martrix is zero, so it cannot be inverted because its determiant (e.det()) is zero.  So how would Mathematica be able to calculate an inverse?

[Oscar Benjamin]

On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <gerardo...@gmail.com> wrote:
>

> Sure, sorry about the image, I was unaware of how easy is to generate python code with the print methods
>
> from sympy import Symbol, Matrix,I
> gamma_c = Symbol('gamma_c')
> n_c = Symbol('n_c')
> gamma_h = Symbol('gamma_h')
> n_h = Symbol('n_h')
> epsilon_c = Symbol('epsilon_c')
> g = Symbol('g')
> epsilon_h = Symbol('epsilon_h')
> 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], [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], [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, 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, 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, 0], [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, 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, 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, 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, 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, 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], [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, 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, 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, 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, 0], [gamma_c*n_c + gamma_h*n_h, 0, 0, 0, 0, gamma_c*n_c + gamma_h*n_h + gamma_h*(n_h + 1), 0, 0, 0, 0, gamma_c*n_c + gamma_c*(n_c + 1) + gamma_h*n_h, 0, 0, 0, 0, 0]])
>
> However, I care more about how to speed up matrix inversion in general, as it was terribly fast in Mathematica and it never finished in sympy. I guess there is a more efficient way than using
>
> e.inv()

The determinant of this matrix is zero so it isn't invertible. I
checked this using DomainMatrix:

In [87]: from sympy.polys.matrices import DomainMatrix

In [88]: dM = DomainMatrix.from_Matrix(e)

In [89]: %time ok = dM.charpoly()
CPU times: user 2min 36s, sys: 100 ms, total: 2min 36s
Wall time: 2min 36s

In [90]: ok[-1] # constant term of charpoly is the determinant
Out[90]: (0 + 0*I)

You can check more quickly just by substituting random values for the symbols:

In [91]: vals = dict(zip(e.free_symbols, randMatrix(7, 1)))

In [92]: e.subs(vals).det()
Out[92]: 0

In [93]: vals
Out[93]: {ε_c: 6, εₕ: 92, g: 44, γ_c: 84, γₕ: 63, n_c: 37, nₕ: 42}

Are you sure that this matrix is the same as the one that you used in
Mathematica?

Did Mathematica actually return an inverse or did it just say that the
matrix is not invertible?

Also I can see now how inv can be made much faster for cases like
this. The matrix has a number of different components:

In [101]: e.connected_components()
Out[101]: [[0, 5, 6, 10, 9, 15], [1, 2, 7, 11], [3], [4, 8, 13, 14], [12]]

The inverse can be computed from the submatrices generated by the
numbered rows and columns. The determinant is the product of their
determinants and the first submatrix has a zero determinant:

In [110]: ccs = e.connected_components()

In [111]: e[ccs[0],ccs[0]].det()
Out[111]: 0

In fact this is zero because its right column is zero which in turn is
because column 15 of the matrix is zero:

In [121]: list(e[:,15])
Out[121]: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Not sure why I didn't notice that straight away!

In [127]: 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 ]
[X X X ]

That's definitely a case that could be handled more efficiently by the
various routines but it makes me wonder if you showed the right
matrix.

--
Oscar

[Peter Stahlecker]

Would Mathematic return a Penrose inverse, if det(M) =0? I do not know Mathematica at all.

--
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 on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTtngh4LHPu2aid9bqycn%2B7cvg0cZ8Ey0_r5ms4CYzQig%40mail.gmail.com.

--

Best regards,

Peter Stahlecker

[Gerardo Suarez]

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

I'll try some of the other suggestions. Thanks a lot Oscar

[Oscar Benjamin]

On Tue, 11 Jan 2022 at 20:48, Gerardo Suarez <gerardo...@gmail.com> 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

[Gerardo Suarez]

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

[Oscar Benjamin]

On Thu, 17 Feb 2022 at 19:04, Gerardo Suarez <gerardo...@gmail.com> wrote:
>
> 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.

I've added a description of what needs to be done in relation to
sparse polynomial gcd here:
https://github.com/sympy/sympy/issues/23131

> 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?

It should just be benchmarked. If you test out the domains yourself
though you will quickly discover what the differences are. The domain
system is discussed here:
https://docs.sympy.org/latest/modules/polys/domainsintro.html

> 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

It shouldn't be too hard to simplify these so maybe there's some other
problem. Can you show the full code that you used and what you
expected as output? If the expressions are quite complicated it might
be better just to focus on the expected output say for the top-left
element of the matrix rather than every element.

--
Oscar

[Oscar Benjamin]

On Sun, 20 Feb 2022 at 13:33, Oscar Benjamin <oscar.j....@gmail.com> wrote:
>
> On Thu, 17 Feb 2022 at 19:04, Gerardo Suarez <gerardo...@gmail.com> wrote:
> >
> > 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.
>
> I've added a description of what needs to be done in relation to
> sparse polynomial gcd here:
> https://github.com/sympy/sympy/issues/23131

I've added this to the list of high priority ideas for a GSOC project:

https://github.com/sympy/sympy/wiki/GSoC-Ideas#polynomial-gcd

--
Oscar