# Solve resolves symbolic linear system in 20 seconds, LinearSolve just grinds on and on

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];


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?

• For me, Solve takes 21 seconds and LinearSolve takes 35 seconds. So LinearSolve is slower, but not that much slower. If I run LinearSolve a second time it returns in 6.3 seconds, due to caching I suppose. Solve takes 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. Commented Mar 26, 2017 at 22:34
• BTW: CoefficientArrays[] is an easier way to get the matrix and vector you need: {rhs, matrix} = Normal[CoefficientArrays[equations, variables]]; rhs = -rhs; Commented Mar 27, 2017 at 0:50
• C.E. Thanks so much for the help! I'm running Mathematica 10.0 and in that version LinearSolve never did come back with an answer, even after several hours. Surprisingly, removing the last variable and equation did the trick. Without those, it solved in 8 seconds. And it's easy to compute that last variable from the answer. J.M. Thanks for the tip. I'll use it! Commented Mar 27, 2017 at 4:19
• Just a thought as I'm no expert. If you think the ZeroTest is the source of the difference use the ZeroTest option to use your own version of ZeroTest, and enforce the two approaches to use the same one, thus checking your hypothesis Commented Mar 27, 2017 at 5:54
• It might be a difference in default setting of Method option to RowReduce. Commented Mar 27, 2017 at 16:39