I am repeatedly solving (related) systems of symbolic linear equations. Solve has worked reliably, but I expected LinearSolve to be faster. I ran into this case, which Solve resolves in ~20 seconds, while LinearSolve seems to be stumped. It just keeps running...
Comment added later: C.E. says LinearSolve gives him an answer in 35 seconds. I let it run for several hours on two different machines and ... nothing. Interesting observation, though: The final equation is complicated, but not complicating -- it is a long involved expression defining a variable that appears in no other equation. Without that variable and equation, LinearSolve finds an answer in 8 or 9 seconds (hallelujah!) and it is easy to back substitute to find this last variable. Odd that something like that should so confound LinearSolve.
The system of equations for Solve: (Sorry, this is the smallest example I've found) has the 14 variables:
variables = {x[2], x[3], x[4], x[5], x[6], z[2], z[3], z[4], z[5], w[3],
w[4], w[5], w[6], y[3]};
and the 14 equations
equations = {1/r[2] - (2 E^(-((2 M r[2])/v^2)) r[2] x[2])/v^2 ==
0,
-s[3]/r[3] + s[3]/r[4] + S[4]/r[3] -
S[4]/r[4] + (-E^(-((2 r[3] s[3])/v^2)) +
E^(-((2 r[3] S[4])/v^2))) x[3] -
(-E^(-((2 r[4] s[3])/v^2)) + E^(-((2 r[4] S[4])/v^2))) x[4] == 0,
-s[4]/r[4] + s[4]/r[5] + S[5]/r[4] -
S[5]/r[5] + (-E^(-((2 r[4] s[4])/v^2)) +
E^(-((2 r[4] S[5])/v^2))) x[4] -
(-E^(-((2 r[5] s[4])/v^2)) + E^(-((2 r[5] S[5])/v^2))) x[5] == 0,
-s[5]/r[5] + s[5]/r[6] + S[6]/r[5] -
S[6]/r[6] + (-E^(-((2 r[5] s[5])/v^2)) +
E^(-((2 r[5] S[6])/v^2))) x[5] -
(-E^(-((2 r[6] s[5])/v^2)) + E^(-((2 r[6] S[6])/v^2))) x[6] == 0,
1/r[6] - (2 r[6] x[6])/v^2 == 0,
1/r[2] - 1/r[6] - (2 E^(-((2 r[2] s[2])/v^2)) r[2] x[2])/v^2 +
(2 E^(-((2 r[6] s[2])/v^2)) r[6] x[6])/v^2 - (h z[2])/r[2] +
(4 E^(-((2 r[2] s[2])/v^2)) Cee[2] r[2]^2 z[2])/v^4 +
(h z[2])/r[6] - (4 E^(-((2 r[6] s[2])/v^2)) Cee[6] r[6]^2 z[2])/v^4 == 0,
1/r[3] - 1/r[4] - (2 E^(-((2 r[3] s[3])/v^2)) r[3] x[3])/v^2 +
(2 E^(-((2 r[4] s[3])/v^2)) r[4] x[4])/v^2 - (h z[3])/r[3] +
(4 E^(-((2 r[3] s[3])/v^2)) Cee[3] r[3]^2 z[3])/v^4 +
(h z[3])/r[4] - (4 E^(-((2 r[4] s[3])/v^2)) Cee[4] r[4]^2 z[3])/v^4 == 0,
1/r[4] - 1/r[5] - (2 E^(-((2 r[4] s[4])/v^2)) r[4] x[4])/v^2 +
(2 E^(-((2 r[5] s[4])/v^2)) r[5] x[5])/v^2 - (h z[4])/r[4] +
(4 E^(-((2 r[4] s[4])/v^2)) Cee[4] r[4]^2 z[4])/v^4 +
(h z[4])/r[5] - (4 E^(-((2 r[5] s[4])/v^2)) Cee[5] r[5]^2 z[4])/v^4 == 0,
1/r[5] - 1/r[6] - (2 E^(-((2 r[5] s[5])/v^2)) r[5] x[5])/v^2 +
(2 E^(-((2 r[6] s[5])/v^2)) r[6] x[6])/v^2 - (h z[5])/r[5] +
(4 E^(-((2 r[5] s[5])/v^2)) Cee[5] r[5]^2 z[5])/v^4 +
(h z[5])/r[6] - (4 E^(-((2 r[6] s[5])/v^2)) Cee[6] r[6]^2 z[5])/v^4 == 0,
-(1/r[2]) + 1/r[3] + (h w[3])/r[2] -
(4 E^(-((2 r[2] S[3])/v^2)) Cee[2] r[2]^2 w[3])/v^4 -
(h w[3])/r[3] + (4 E^(-((2 r[3] S[3])/v^2)) Cee[3] r[3]^2 w[3])/v^4 +
(2 E^(-((2 r[2] S[3])/v^2)) r[2] x[2])/v^2 -
(2 E^(-((2 r[3] S[3])/v^2)) r[3] x[3])/v^2 == 0,
-(1/r[3]) + 1/r[4] + (h w[4])/r[3] -
(4 E^(-((2 r[3] S[4])/v^2)) Cee[3] r[3]^2 w[4])/v^4 -
(h w[4])/r[4] + (4 E^(-((2 r[4] S[4])/v^2)) Cee[4] r[4]^2 w[4])/v^4 +
(2 E^(-((2 r[3] S[4])/v^2)) r[3] x[3])/v^2 -
(2 E^(-((2 r[4] S[4])/v^2)) r[4] x[4])/v^2 == 0,
-(1/r[4]) + 1/r[5] + (h w[5])/r[4] -
(4 E^(-((2 r[4] S[5])/v^2)) Cee[4] r[4]^2 w[5])/v^4 -
(h w[5])/r[5] + (4 E^(-((2 r[5] S[5])/v^2)) Cee[5] r[5]^2 w[5])/v^4 +
(2 E^(-((2 r[4] S[5])/v^2)) r[4] x[4])/v^2 -
(2 E^(-((2 r[5] S[5])/v^2)) r[5] x[5])/v^2 == 0,
-(1/r[5]) + 1/r[6] + (h w[6])/r[5] -
(4 E^(-((2 r[5] S[6])/v^2)) Cee[5] r[5]^2 w[6])/v^4 -
(h w[6])/r[6] + (4 E^(-((2 r[6] S[6])/v^2)) Cee[6] r[6]^2 w[6])/v^4 +
(2 E^(-((2 r[5] S[6])/v^2)) r[5] x[5])/v^2 -
(2 E^(-((2 r[6] S[6])/v^2)) r[6] x[6])/v^2 == 0,
y[3] == 1/2 (-2 (-E^(-((2 r[2] s[2])/v^2)) +
E^(-((2 r[2] S[3])/v^2))) x[2] -
((-h v^2 + 2 r[2] (-t + p r[2]) +
h r[2] (s[2] + S[3])) (-w[3] + z[2]))/r[2]^2 -
2 Cee[2] (-((2 E^(-((2 r[2] S[3])/v^2)) r[2] w[3])/v^2) +
(2 E^(-((2 r[2] s[2])/v^2)) r[2] z[2])/v^2) -
((s[2] - S[3]) (-2 r[2] + h r[2] (w[3] + z[2])))/r[2]^2) +
1/2 (2 (-E^(-((2 r[3] s[3])/v^2)) +
E^(-((2 r[3] S[3])/v^2))) x[3] +
((-h v^2 + 2 r[3] (-t + p r[3]) +
h r[3] (s[3] + S[3])) (-w[3] + z[3]))/r[3]^2 +
2 Cee[3] (-((2 E^(-((2 r[3] S[3])/v^2)) r[3] w[3])/v^2) +
(2 E^(-((2 r[3] s[3])/v^2)) r[3] z[3])/v^2) +
((s[3] - S[3]) (-2 r[3] + h r[3] (w[3] + z[3])))/r[3]^2) +
1/2 (2 (E^(-((2 r[4] s[3])/v^2)) - E^(-((2 r[4] s[4])/v^2))) x[4] -
((-h v^2 + 2 r[4] (-t + p r[4]) +
h r[4] (s[3] + s[4])) (z[3] - z[4]))/r[4]^2 +
2 Cee[4] (-((2 E^(-((2 r[4] s[3])/v^2)) r[4] z[3])/v^2) +
(2 E^(-((2 r[4] s[4])/v^2)) r[4] z[4])/v^2) -
((s[3] - s[4]) (-2 r[4] + h r[4] (z[3] + z[4])))/r[4]^2) +
1/2 (2 (-E^(-((2 r[6] s[2])/v^2)) +
E^(-((2 r[6] s[5])/v^2))) x[6] +
((-h v^2 + 2 r[6] (-t + p r[6]) +
h r[6] (s[2] + s[5])) (z[2] - z[5]))/r[6]^2 +
2 Cee[6] ((2 E^(-((2 r[6] s[2])/v^2)) r[6] z[2])/v^2 -
(2 E^(-((2 r[6] s[5])/v^2)) r[6] z[5])/v^2) +
((s[2] - s[5]) (-2 r[6] + h r[6] (z[2] + z[5])))/r[6]^2) +
1/2 (2 (E^(-((2 r[5] s[4])/v^2)) -
E^(-((2 r[5] s[5])/v^2))) x[5] -
((-h v^2 + 2 r[5] (-t + p r[5]) +
h r[5] (s[4] + s[5])) (z[4] - z[5]))/r[5]^2 +
2 Cee[5] (-((2 E^(-((2 r[5] s[4])/v^2)) r[5] z[4])/v^2) +
(2 E^(-((2 r[5] s[5])/v^2)) r[5] z[5])/v^2) -
((s[4] - s[5]) (-2 r[5] + h r[5] (z[4] + z[5])))/r[5]^2)};
It takes 20 seconds or so, but
sol = Solve[equations, variables];
produces the correct answer.
In an effort to speed things along, ToMatrix and GetRHS converts this to a linear system matrix.variables == rhs
ToMatrix[equations_, variables_] :=
Map[Function[poly, Map[Coefficient[poly, #] &, variables]],
Map[(#[[1]] - #[[2]]) &, equations]]
GetRHS[equations_, variables_] :=
Transpose[{Simplify[Map[(#[[2]] - #[[1]]) &, equations] +
ToMatrix[equations, variables].variables]}]
So we can solve the system with LinearSolve
matrix = ToMatrix[equations, variables];
rhs = GetRHS[equations, variables];
lsol = LinearSolve[matrix, rhs];
Sadly, it's still running. The only reason I can think of is LinearSolve is spending enormous time on ZeroTest. But don't the two use the same test?
Solvetakes 21 seconds andLinearSolvetakes 35 seconds. SoLinearSolveis slower, but not that much slower. If I runLinearSolvea second time it returns in 6.3 seconds, due to caching I suppose.Solvetakes 19 seconds if I run it a second time, so it does not benefit nearly as much from caching. I'm using Mathematica 11.1. $\endgroup$CoefficientArrays[]is an easier way to get the matrix and vector you need:{rhs, matrix} = Normal[CoefficientArrays[equations, variables]]; rhs = -rhs;$\endgroup$Methodoption toRowReduce. $\endgroup$