2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – C. E. Mar 26 '17 at 22:34
  • $\begingroup$ BTW: CoefficientArrays[] is an easier way to get the matrix and vector you need: {rhs, matrix} = Normal[CoefficientArrays[equations, variables]]; rhs = -rhs; $\endgroup$ – J. M. will be back soon Mar 27 '17 at 0:50
  • $\begingroup$ 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! $\endgroup$ – JV3 Mar 27 '17 at 4:19
  • $\begingroup$ 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 $\endgroup$ – Yair M Mar 27 '17 at 5:54
  • $\begingroup$ It might be a difference in default setting of Method option to RowReduce. $\endgroup$ – Daniel Lichtblau Mar 27 '17 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.