On Fri, 25 Aug 2023 at 14:08, Ralf Schlatterbeck <r...@runtux.com> wrote:
>
> Hi all,
>
> [This is quite long, tl;dr: I cannot solve a matrix in reasonable time
> that can be solved with Maxima in 8s.]
>
> I'm using Lcapy (which in turn uses sympy) for computing a transfer
> function of a simple circuit (see below for an ascii-art drawing).
>
> I was referred here by the author of Lcapy.
>
> As many of you may be aware one can compute the parameters of a circuit
> by setting up a matrix and solving for the voltages (or currents). See
> at the end for an ascii-art of the circuit.
>
> Lcapy comes up with the following matrix:
>
> from sympy import Matrix, symbols
> R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3 = symbols(
> 'R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3', positive=True,
> real=True)
> s = symbols('s')
> b = Matrix([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1])
> A = Matrix([
> [1/R1, -1/R1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
> [-1/R1, C1*s + 1/R1, -C1*s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
> [0, -C1*s, C1*s + C2*s + C3*s, -C2*s, -C3*s, 0, 0, 0, 0, 0, 0, 0, 0],
> [0, 0, -C2*s, C2*s, 0, 0, 0, 0, 0, 1, 0, 0, 0],
> [0, 0, -C3*s, 0, C3*s + C4*s + C5*s, -C4*s, -C5*s, 0, 0, 0, 0, 0, 0],
> [0, 0, 0, 0, -C4*s, C4*s, 0, 0, 0, 0, 1, 0, 0],
> [0, 0, 0, 0, -C5*s, 0, C5*s + C6*s + C7*s, -C7*s, -C6*s, 0, 0, 0, 0],
> [0, 0, 0, 0, 0, 0, -C7*s, C7*s, 0, 0, 0, 1, 0],
> [0, 0, 0, 0, 0, 0, -C6*s, 0, C6*s + 1/R2, 0, 0, 0, 0],
> [0, 0, 0, 1, 0, 0, 0, 0, 0, -L1*s, 0, 0, 0],
> [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -L2*s, 0, 0],
> [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -L3*s, 0],
> [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
>
> x = A.solve(b, method='LU')
>
> Trying to solve this with the 'LU' solver produces a really huge
> solution (but the solution is quite fast):
>
> def exsize (expr):
> n = 0
> if not expr.args:
> return 1
> for arg in expr.args:
> n += exsize (arg)
> return n
>
> >>> print ([exsize (k) for k in x])
> [52298, 26148, 884884, 128156, 730566, 65218, 510990, 35741, 255579,
> 128154, 65216, 35739, 26146]
>
> [Is there a better way to determine the size of an expression?]
>
> Simplifying the first expression (x [0]) should yield 1 (probably using
> scipy.cancel). But I'm too impatient to wait for it, I've killed it when
> it took 1.2GB resident memory after an hour or so.
>
> Other solvers are *very* slow. Trying to solve with default 'GJ' method
> wasn't completed in several hours, so I left it running overnight and
> had terminated in the morning, unfortunately I had not computed the
> sizes as above.
>
> The size of the result created by 'LU' prevents further computations to
> happen in reasonable time.
These problems are both well known (to me at least). Some solvers like
LU produce large unsimplified output that potentially includes divide
by zero within the large expressions. Others like GJ would produce
correct and simplified output but are extremely slow. See also the
discussion here:
https://github.com/sympy/sympy/issues/24679
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.
This computes a matrix numerator and scalar denominator that when
divided represent the solution to your matrix equation. We can perform
the division either with or without cancellation:
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
In [9]: %time sol2 = sol_num.to_Matrix() / den.as_expr()
CPU times: user 104 ms, sys: 1.88 ms, total: 106 ms
Wall time: 105 ms
I won't show the full output because it is complicated but sol1 is in
the canonical form of a ratio of expanded polynomials with gcd
cancelled. Both solutions have a 1 in the leading entry as expected:
In [14]: sol1[0]
Out[14]: 1
In [15]: sol2[0]
Out[15]: 1
The two matrices are not otherwise equal because the expressions in
sol2 are not fully cancelled:
In [16]: sol1 == sol2
Out[16]: False
In [17]: sol1 == sol2.applyfunc(cancel)
Out[17]: True
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]])
In [20]: print(sol2.applyfunc(lambda e: e.count_ops()))
Matrix([[0], [2595], [2119], [1891], [1871], [1751], [1743], [1679],
[1672], [1859], [1735], [1671], [2140]])
In [22]: print(A.solve(b, method='LU').applyfunc(lambda e: e.count_ops()))
Matrix([[45301], [22650], [766394], [110995], [632734], [56484],
[442554], [30946], [221347], [110993], [56482], [30944], [22648]])
So both sol2 and sol1 are a lot simpler than what LU gives (on this
measure) and they both compute much more quickly than GJ.
Installing gmpy2 can help with timings for DomainMatrix in some cases
(although probably does not make much difference in this example). In
future it will be possible to speed this up a lot using python-flint.
Specifically python-flint will provide fast multiplication and gcd for
sparse multivariate polynomials which are the bottlenecks in the two
steps shown above.
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CAHVvXxQqwsfnmKn%2BvA3%2BLZ3AHytGnuC_to4%2BrFnnu6gg6WSqeg%40mail.gmail.com.
