Solve slower in versions 9 and 10

Question

With some systems of equations, Solve is much slower in versions 9 and 10 than in earlier versions, apparently because it is doing more simplification of the results.

With the following example linear system:

eqs = {0 == g1 - 2 g1 r[1, 1] - I oa r[1, 3] + I oa r[3, 1] + g2 r[3, 3],
   0 == -I da r[1, 2] + I db r[1, 2] - 2 g1 r[1, 2] - I ob r[1, 3] + I oa r[3, 2],
   0 == -I oa r[1, 1] - I ob r[1, 2] - I da r[1, 3] - 2 g1 r[1, 3] - g2 r[1, 3] + I oa r[3, 3],
   0 == I da r[2, 1] - I db r[2, 1] - 2 g1 r[2, 1] - I oa r[2, 3] + I ob r[3, 1],
   0 == g1 - 2 g1 r[2, 2] - I ob r[2, 3] + I ob r[3, 2] + g2 r[3, 3],
   0 == -I oa r[2, 1] - I ob r[2, 2] - I db r[2, 3] - 2 g1 r[2, 3] - g2 r[2, 3] + I ob r[3, 3],
   0 == I oa r[1, 1] + I ob r[2, 1] + I da r[3, 1] - 2 g1 r[3, 1] - g2 r[3, 1] - I oa r[3, 3],
   0 == I oa r[1, 2] + I ob r[2, 2] + I db r[3, 2] - 2 g1 r[3, 2] - g2 r[3, 2] - I ob r[3, 3],
   0 == I oa r[1, 3] + I ob r[2, 3] - I oa r[3, 1] - I ob r[3, 2] - 2 g1 r[3, 3] - 2 g2 r[3, 3]
   };
vars = {r[1, 1], r[1, 2], r[1, 3], r[2, 1], r[2, 2], r[2, 3], r[3, 1], r[3, 2], r[3, 3]};

{$VersionNumber, Round[First@Timing[sol = Solve[eqs, vars];], .01], LeafCount@sol}

I get the results:

\begin{array}{ccc} \text{Version} & \text{Timing} & \text{LeafCount} \\ 6. & 0.19 & 219384 \\ 7. & 0.08 & 227942 \\ 8. & 1.82 & 86317 \\ 9. & 63.32 & 29452 \\ 10.2 & 30.84 & 82043 \\ \end{array}

Version 10 was run on a 2.3GHz MacBook Pro, earlier versions were run on an older 2.4GHz PC. (Note that if the Solve command is executed twice, the second time is faster, presumably due to some caching of results. The timings are all for the first evaluation.)

So there is an inverse correlation between Timing and LeafCount, which makes sense if extra time is being taken for some sort of simplification. But for my purposes the ~3x reduction in LeafCount between versions 7 and 10 is not worth the ~400x slowdown.

Does anyone happen to know if there's an undocumented Method setting to get Solve to use a method from an earlier version? Or maybe some other workaround?

Update 1:

Solve with Method -> "Legacy" (suggested by Alexey Popkov) and System`Private`OldSolve (suggested by Guess who) both give me a timing of 18.7 s and a LeafCount of 44087 on 10.2. So it seems that they are both doing the same thing. It's supposed to cause Solve to use the version 7 algorithm, but it's clearly not using the exact algorithm from version 7. Using Method -> "Legacy" in different versions I get

\begin{array}{ccc} \text{Version} & \text{Timing} & \text{LeafCount} \\ 8. & 1.78 & 48622 \\ 9. & 21.92 & 27358 \\ 10.2 & 18.77 & 44087 \\ \end{array}

So for this system of equations it's an improvement in both timing and leaf count in all versions, but still nowhere close to version 7's speed.

Update 2:

As pointed out in Michael E2's answer, we can do much better with LinearSolve, and we can transfer this benefit to Solve by setting the options for RowReduce, as suggested by Daniel Lichtblau.

After setting

SetOptions[RowReduce, Method -> "CofactorExpansion"]

Solve gives the results

\begin{array}{ccc} \text{Version} & \text{Timing} & \text{LeafCount} \\ 6. & 0.42 & 64004 \\ 7. & 0.42 & 64004 \\ 8. & 0.5 & 67879 \\ 9. & 0.5 & 67879 \\ 10.2 & 0.27 & 67879 \\ \end{array}

Now it's reasonably fast, and stable across versions. It doesn't beat the default Solve in version 7 for speed, or version 9 for leaf count, but it seems like a good overall compromise.

Answer 1

This is best I can do so far. The system is linear so LinearSolve is a natural thing to try.

arrays = CoefficientArrays[eqs, vars]
(*  {SparseArray[< 2 >, {9}], SparseArray[< 35 >, {9, 9}]}  *)

solv1 = Thread[
    vars -> LinearSolve[arrays[[2]], -arrays[[1]], 
      Method -> "CofactorExpansion"]]; // AbsoluteTiming
(*  {0.28347, Null}  *)

I cannot figure out how to get Solve to make this call or its equivalent.

---

With Daniel and Simon's help, here is how to do it with Solve:

With[{opts = Options[RowReduce]},
 Internal`WithLocalSettings[
  SetOptions[RowReduce, Method -> "CofactorExpansion"],
  Solve[eqs, vars]; // AbsoluteTiming,
  SetOptions[RowReduce, opts]
  ]]
(*  {0.282098, Null}  *)