On Fri, 25 Aug 2023 at 18:38, Ralf Schlatterbeck <[email protected]> wrote: > > Hello Oscar, > First, thanks for the very quick answer! > > On Fri, Aug 25, 2023 at 03:19:31PM +0100, Oscar Benjamin wrote: > > > > > > Trying to solve this with the 'LU' solver produces a really huge > > > solution (but the solution is quite fast): > ... > > > The size of the result created by 'LU' prevents further computations to > > > happen in reasonable time. > > > > I have been working on improving this and will write about this in > > more detail soon. For now if you use the current SymPy master branch > > then it is possible to get a result like the GJ one but a lot faster > > if you use DomainMatrix instead of Matrix. With the current master > > branch Matrix has a .to_DM() method for converting a Matrix to a > > DomainMatrix: > > > > In [7]: %time sol_num, den = A.to_DM().solve_den(b.to_DM()) > > CPU times: user 20.6 s, sys: 236 ms, total: 20.8 s > > Wall time: 22.2 s > > So currently that takes 20 seconds on this (not very powerful) laptop. > ... > > In [8]: %time sol1 = (sol_num.to_field() / den).to_Matrix() > > CPU times: user 369 ms, sys: 2.67 ms, total: 372 ms > > Wall time: 372 ms > ... > > Nice. > On my machine this takes 1.1 to 1.3 seconds including the sol1 > computation.
Actually I was timing with an out of date branch. With actual current master it takes 2.4 seconds total for sol1 (on the same laptop). > So this is really what I was looking for! Great! > Has the DM implementation been changed recently? Because lcapy already > has code in the last incarnation to use DM but it still takes ages for > me with sympy 1.12, the code for older sympy uses: > > from sympy.polys.domainmatrix import DomainMatrix > dM = DomainMatrix.from_list_sympy( > *M.shape, rows=M.tolist()).to_field() > return dM.inv().to_Matrix() > > and then multiplies the inverted matrix with b. The implementation has changed. For one thing the solve_den method is new. Also there is now inv_den which is currently much faster than inv in most cases. It is all in a bit of a state of flux at the moment but if interested in further reading then see e.g.: https://github.com/sympy/sympy/pull/25336 https://github.com/sympy/sympy/pull/25346 https://github.com/sympy/sympy/pull/25367 https://github.com/sympy/sympy/pull/25395 https://github.com/sympy/sympy/pull/25443 https://github.com/sympy/sympy/pull/25458 https://github.com/sympy/sympy/pull/25495 This one is not yet finished: https://github.com/sympy/sympy/pull/25452 See also: https://github.com/sympy/sympy/issues/25403 https://github.com/sympy/sympy/issues/25410 Also I wrote the first of several blog posts here: https://oscarbenjamin.github.io/blog/czi/post1.html A future blog post will explain what is going on with DomainMatrix etc. The situation with all of the matrix types is a bit complicated now and I am not sure even how many SymPy developers understand how it all fits together any more. Ultimately the intention is that the Matrix methods will just call through to DomainMatrix automatically so that hopefully downstream code does not need to use DomainMatrix directly and everything just becomes faster, simpler, less bug-prone, etc. > > The usual method for measuring expression size in SymPy is with > > "operation count". There is nothing particularly special about this as > > a measure but it is at least easy to compute because a method is ready > > made for computing it: > > > > In [19]: print(sol1.applyfunc(lambda e: e.count_ops())) > > Matrix([[0], [2472], [2030], [1823], [1795], [1683], [1672], [1611], > > [1610], [1783], [1664], [1602], [2059]]) > > OK, I've timed that measure on the LU solution and it takes > significantly longer than my simple recursive node counter. > > For x [2]: > My node counter: > 2: size: 884884 time: 0.22558188438415527 > count_ops: > 2: size: 766394 time: 6.558408737182617 > > I had hoped for something that could tell me quickly how long that > result would be *without* the risk that I wait for hours. An occasion I > waited for hours is the solutions produced by the 'CH' and 'LDL' solver > methods (whatever they do) which seem to produce so large solutions that > my node counter doesn't finish in hours :-) Yes, count_ops is slow. It tries to do too much and then ends up being slow for the more important common case of just wanting a quick measure. It is also possible to make something *much* faster than your recursive version when there are many repeating subexpressions as there will be in these cases (more on this later but a simple fix is just to use a cache like lru_cache). -- 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/CAHVvXxQ7eUYOuR0nDxPCfoCh9TLjL8MgLCD6p%3DVv%3D0z_oCH1UQ%40mail.gmail.com.
