Matrix cannot be solved in reasonable time

63 views
Skip to first unread message

Ralf Schlatterbeck

unread,
Aug 25, 2023, 9:08:35 AM8/25/23
to sy...@googlegroups.com
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.

What I'm trying to achieve here is to compute the transfer function "H"
of the circuit, this is V(9)/V(1) (see circuit for points). It is x [8]
in the solution above (if it could be computed with sympy). I can do
this in sympy as below (without computing the voltages using the matrix,
instead the solution involves viewing the circuit as a cascade of
voltage dividers) in less than a second. The result is of reasonable
size (less than 1k).

I can successfully solve the matrix above in Maxima in 8 seconds
yielding a reasonably-sized result that can be further processed (with
ratsimp in my case). With ratsimp (x [1]) I indeed get 1 (indexes are
1-based in Maxima). And x [9] yields the correct transfer function.

Now my question: Should I expect such a huge matrix solution?
Is there another solver and/or another way to produce a reasonably-sized
result in reasonable time (several seconds)? Or can I force repeated
cancellation during matrix solving to prevent it from getting so big?

The state of the art in all intro texts on circuit simulation seems to
be to solve a matrix for the currents or voltages, should we use a
different method here (like the voltage divider solution below which is
applicable only to that particular ladder topology)?

Fast code to compute transfer function H for circuit below using voltage
dividers for the ladder components:

from sympy import solve, symbols
R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3, U1 = symbols(
'R1, R2, C1, C2, C3, C4, C5, C6, C7, L1, L2, L3, U1', positive=True,
real=True)
s = symbols('s')

R_V = R1 + 1 / (s * C1)

R_U7 = 1/(1/(1/(s*C7)+s*L3)+1/(1/(s*C6)+R2))
R_U5 = 1/(1/(1/(s*C4)+s*L2)+1/(1/(s*C5)+R_U7))
R_U3 = 1/(1/(1/(s*C2)+s*L1)+1/(1/(s*C3)+R_U5))

U3 = (R_U3 / (R_V + R_U3)) * U1
U5 = (R_U5 / (1 / (s * C3) + R_U5)) * U3
U7 = (R_U7 / (1 / (s * C5) + R_U7)) * U5
U9 = (R2 / (1 / (s * C6) + R2)) * U7

H = U9 / U1

[This should be identical to x [8] in the matrix solution above which so
far I've not been able to compute with sympy in a useable form]


Circuit ASCII-Art:

1 +----+ 2 C1 3 C3 5 C5 7 C6 9
o--| R1 |--+--| |---+---| |---+---| |---+---| |---+-----o
+----+ | | | |
| C2 | C4 | C7 |
--- --- --- ---
--- --- --- | | R2
|4 |6 |8 | |
> > > ---
> L1 > L2 > L3 |
> > > |
| | | |
o-------------------+---------+---------+---------+-----o

Thanks
Ralf
--
Dr. Ralf Schlatterbeck Tel: +43/2243/26465-16
Open Source Consulting www: www.runtux.com
Reichergasse 131, A-3411 Weidling email: off...@runtux.com

Oscar Benjamin

unread,
Aug 25, 2023, 10:19:51 AM8/25/23
to sy...@googlegroups.com
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

Ralf Schlatterbeck

unread,
Aug 25, 2023, 1:38:04 PM8/25/23
to sy...@googlegroups.com
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.
So this is really what I was looking for!

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 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 :-)

I'm looking forward to your further improvements!

Thanks again!

Oscar Benjamin

unread,
Aug 26, 2023, 5:26:21 AM8/26/23
to sy...@googlegroups.com
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

Ralf Schlatterbeck

unread,
Sep 4, 2023, 8:03:38 AM9/4/23
to sy...@googlegroups.com
On Sat, Aug 26, 2023 at 10:26:02AM +0100, Oscar Benjamin wrote:
>
> It is all in a bit of a state of flux at the moment but if interested
> in further reading then see e.g.:

[...]
>
> 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.

I've read throught this over the weekend.

The blog-post is the most accessible to me, so I'd appreciate when you
continue this series.
For me it was escpecially nice to see that what you write in the blog
post makes the difference if solving a matrix (in my case 13x13 as in a
previous message in this thread) is possible or not.

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

Yes, that would be very useable.

Thanks!
Reply all
Reply to author
Forward
0 new messages