Comment #20 on issue 2015 by [email protected]: Hangs attempting to solve a
system of linear equations
http://code.google.com/p/sympy/issues/detail?id=2015
I did some experimentation and using solve_lin_sys() from #1850 and sparse
rational functions from sparse-polys branch I was able to get the same
result as Sage gives for the simplified case in 1.6 seconds and as Maple
gives for the original case in 3 seconds (on i7-2600K @ 3.4 GHz):
In [1]: from sympy.polys.solving import solve_lin_sys
In [2]: F, d,r,e,f,g,i,j,l,o,m,p,q,c,h,k,n,b =
field("d,r,e,f,g,i,j,l,o,m,p,q,c,h,k,n,b", QQ)
In [3]: G = [b + q/d - c/d, c*(1/d + 1/e + 1/g) - f/g - q/d, f*(1/g + 1/i +
1/j) - c/g - h/i, h*(1/i + 1/l + 1/m) - f/i - k/m, k*(1/m + 1/o + 1/p) -
h/m - n/p, n/p - k/p]
In [4]: %time S = solve_lin_sys(G, [c,f,h,k,n,b], F)
CPU times: user 1.61 s, sys: 0.04 s, total: 1.65 s
Wall time: 1.60 s
In [5]: %paste
T = {b: (e*i*l*q + e*i*m*q + e*i*o*q + e*j*l*q + e*j*m*q + e*j*o*q +
e*l*m*q + e*l*o*q + g*i*l*q + g*i*m*q + g*i*o*q + g*j*l*q + g*j*m*q +
g*j*o*q + g*l*m*q + g*l*o*q + i*j*l*q + i*j*m*q + i*j*o*q + j*l*m*q +
j*l*o*q)/(-d*e*i*l - d*e*i*m - d*e*i*o - d*e*j*l - d*e*j*m - d*e*j*o -
d*e*l*m - d*e*l*o - d*g*i*l - d*g*i*m - d*g*i*o - d*g*j*l - d*g*j*m -
d*g*j*o - d*g*l*m - d*g*l*o - d*i*j*l - d*i*j*m - d*i*j*o - d*j*l*m -
d*j*l*o - e*g*i*l - e*g*i*m - e*g*i*o - e*g*j*l - e*g*j*m - e*g*j*o -
e*g*l*m - e*g*l*o - e*i*j*l - e*i*j*m - e*i*j*o - e*j*l*m - e*j*l*o),
c: (-e*g*i*l*q - e*g*i*m*q - e*g*i*o*q - e*g*j*l*q - e*g*j*m*q - e*g*j*o*q
- e*g*l*m*q - e*g*l*o*q - e*i*j*l*q - e*i*j*m*q - e*i*j*o*q - e*j*l*m*q -
e*j*l*o*q)/(-d*e*i*l - d*e*i*m - d*e*i*o - d*e*j*l - d*e*j*m - d*e*j*o -
d*e*l*m - d*e*l*o - d*g*i*l - d*g*i*m - d*g*i*o - d*g*j*l - d*g*j*m -
d*g*j*o - d*g*l*m - d*g*l*o - d*i*j*l - d*i*j*m - d*i*j*o - d*j*l*m -
d*j*l*o - e*g*i*l - e*g*i*m - e*g*i*o - e*g*j*l - e*g*j*m - e*g*j*o -
e*g*l*m - e*g*l*o - e*i*j*l - e*i*j*m - e*i*j*o - e*j*l*m - e*j*l*o),
f: (-e*i*j*l*q - e*i*j*m*q - e*i*j*o*q - e*j*l*m*q - e*j*l*o*q)/(-d*e*i*l
- d*e*i*m - d*e*i*o - d*e*j*l - d*e*j*m - d*e*j*o - d*e*l*m - d*e*l*o -
d*g*i*l - d*g*i*m - d*g*i*o - d*g*j*l - d*g*j*m - d*g*j*o - d*g*l*m -
d*g*l*o - d*i*j*l - d*i*j*m - d*i*j*o - d*j*l*m - d*j*l*o - e*g*i*l -
e*g*i*m - e*g*i*o - e*g*j*l - e*g*j*m - e*g*j*o - e*g*l*m - e*g*l*o -
e*i*j*l - e*i*j*m - e*i*j*o - e*j*l*m - e*j*l*o),
h: (-e*j*l*m*q - e*j*l*o*q)/(-d*e*i*l - d*e*i*m - d*e*i*o - d*e*j*l -
d*e*j*m - d*e*j*o - d*e*l*m - d*e*l*o - d*g*i*l - d*g*i*m - d*g*i*o -
d*g*j*l - d*g*j*m - d*g*j*o - d*g*l*m - d*g*l*o - d*i*j*l - d*i*j*m -
d*i*j*o - d*j*l*m - d*j*l*o - e*g*i*l - e*g*i*m - e*g*i*o - e*g*j*l -
e*g*j*m - e*g*j*o - e*g*l*m - e*g*l*o - e*i*j*l - e*i*j*m - e*i*j*o -
e*j*l*m - e*j*l*o),
k: e*j*l*o*q/(d*e*i*l + d*e*i*m + d*e*i*o + d*e*j*l + d*e*j*m + d*e*j*o +
d*e*l*m + d*e*l*o + d*g*i*l + d*g*i*m + d*g*i*o + d*g*j*l + d*g*j*m +
d*g*j*o + d*g*l*m + d*g*l*o + d*i*j*l + d*i*j*m + d*i*j*o + d*j*l*m +
d*j*l*o + e*g*i*l + e*g*i*m + e*g*i*o + e*g*j*l + e*g*j*m + e*g*j*o +
e*g*l*m + e*g*l*o + e*i*j*l + e*i*j*m + e*i*j*o + e*j*l*m + e*j*l*o),
n: e*j*l*o*q/(d*e*i*l + d*e*i*m + d*e*i*o + d*e*j*l + d*e*j*m + d*e*j*o +
d*e*l*m + d*e*l*o + d*g*i*l + d*g*i*m + d*g*i*o + d*g*j*l + d*g*j*m +
d*g*j*o + d*g*l*m + d*g*l*o + d*i*j*l + d*i*j*m + d*i*j*o + d*j*l*m +
d*j*l*o + e*g*i*l + e*g*i*m + e*g*i*o + e*g*j*l + e*g*j*m + e*g*j*o +
e*g*l*m + e*g*l*o + e*i*j*l + e*i*j*m + e*i*j*o + e*j*l*m + e*j*l*o)}
## -- End pasted text --
In [6]: S == T
Out[6]: True
In [7]: H = [b + r/d - c/d, c*(1/d + 1/e + 1/g) - f/g - r/d, f*(1/g + 1/i +
1/j) - c/g - h/i, h*(1/i + 1/l + 1/m) - f/i - k/m, k*(1/m + 1/o + 1/p) -
h/m - n/p, n*(1/p + 1/q) - k/p]
In [8]: H
Out[8]:
[(d*b + r - c)/d, (-d*e*f + d*e*c + d*g*c - r*e*g + e*g*c)/d*e*g, (f*g*i +
f*g*j + f*i*j - g*j*h - i*j*c)/g*i*j, (-f*l*m + i*l*h - i*l*k +
i*m*h + l*m*h)/i*l*m, (o*m*k - o*m*n - o*p*h + o*p*k + m*p*k)/o*m*p, (p*n -
q*k + q*n)/p*q]
In [12]: %time S = solve_lin_sys(H, [c,f,h,k,n,b], F)
CPU times: user 3.00 s, sys: 0.00 s, total: 3.00 s
Wall time: 2.93 s
The last result is the same as Maple gives. I translated Maple's output
using macros in vim and compared results in Python (it would be cool to
have a parser for Maple's language).
There is room for improvements because ~90% of time is spent computing
multivariate GCDs (simplifying rational functions). ~90% of that time is
spent dividing multivariate polynomials. This is because HEUGCD is used
which divides polynomials a lot. Better multivariate division and/or GCD
algorithms would help a lot (maybe also smarter simplification of rational
functions which avoids unnecessary GCDs).
On my computer the approach from comment #19 can give a "solution" faster,
but really one can't call this a solution when taking into account sizes of
output expressions:
In [1]: var('d,r,e,f,g,i,j,l,o,m,p,q,c,h,k,n,b')
Out[1]: (d, r, e, f, g, i, j, l, o, m, p, q, c, h, k, n, b)
In [2]: %time ans=solve([b + r/d - c/d, c*(1/d + 1/e + 1/g) - f/g - r/d,
f*(1/g + 1/i + 1/j) - c/g - h/i, h*(1/i + 1/l + 1/m) - f/i - k/m, k*(1/m +
1/o + 1/p) - h/m - n/p, n*(1/p + 1/q) - k/p], b, c, f, h, k, n, manual=True)
CPU times: user 0.60 s, sys: 0.03 s, total: 0.63 s
Wall time: 0.60 s
In [3]: [ expr.count_ops() for expr in ans[0] ]
Out[3]: [68630, 68628, 68595, 22786, 2390, 1195]
Parameter `a' is not used but it doesn't change anything.
Most of this code is available through sparse-polys, but not the solver.
Unfortunately, I had to make a lot of changes to matrices to disable
sympyfication (on demand), so I will publish this to another branch when I
figure out what to do with matrices.
--
You received this message because this project is configured to send all
issue notifications to this address.
You may adjust your notification preferences at:
https://code.google.com/hosting/settings
--
You received this message because you are subscribed to the Google Groups
"sympy-issues" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sympy-issues?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
