# Solving large linear systems of equations efficiently?

I need to solve linear systems of equations of approximate size $(n!)\times(n!)$ as efficiently as possible for as high parameter n as possible. All the entries will be rational, or equivalently integer. To test the routine provided by Mathematica I do the following:

SetAttributes[timeAvg, HoldFirst]
timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]
times = Table[
size = j!;
mat = Table[RandomInteger[{-1000, 1000}], {i, 1, size}, {j, 1, size }];
b = Table[RandomInteger[{-1000, 1000}], {i, 1, size}];
timeAvg[LinearSolve[mat, b]]
, {j, 0, 6}];


Where timeAvg is the mathematica.stackexchange site standard time averaging function. The result is

ListLogPlot[Table[{i-1, times[[i]]}, {i, 1, times//Length}], Filling->Bottom, Joined->True, PlotRange->{{0, 6}, {10^-6, 1.5}}, BaseStyle->{FontWeight->"Bold", FontSize->16}]


I show the plot up to the size 6! because that is the size which still evaluates fine. If I go to 7! Mathematica calculates indefinitely and never produces a result. My question is, how can I push the boundary and get results even at 7!,8! or even 9!? Is that feasible at all? I realize that these sizes quickly become astronomical, but maybe there is a trick or some powerful technique beyond the standard routine? Thanks for any suggestion!

EDIT:

Dr. belisarius suggested to look at numerical solutions. If I try 7! with numerical entries, it evaluates in a few seconds. However, 8! again becomes way too slow. Any hint for 8! and higher?

• Are you after exact or numerical solutions? – Dr. belisarius Feb 25 '16 at 14:40
• @Dr.belisarius I am interested in both! – Kagaratsch Feb 25 '16 at 14:41
• Can you try changing the setting of LinearSolve[] to Method -> "OneStepRowReduction" or Method -> "DivisionFreeRowReduction"? – J. M. is in limbo Feb 25 '16 at 14:56
• How much memory do you have? An (8!)^2 matrix is like 12Gb just to store the matrix at machine precision. – george2079 Feb 25 '16 at 15:08
• LinearSolve with default method appears to use 2x+ the matrix size. 8! quickly maxed out my 32Gb windows machine. If you really needed to do this I'd imagine there are in-place solvers that wouldn't need 2x memory. – george2079 Feb 26 '16 at 16:21

(1) I don't see much hope for this. The algorithm in use is, I believe, quite good for dense integer matrices. It will scale at around O(n^4) and realistically you won't be able to improve on that. So given the time it takes at 6!, expect a factor of 7^4 increase to get to that next level.

(2) The numeric case will be better behaved in that the scale is O(n^3) and also the size of result scales as O(n) (for the integer case expect O(n^2) due to the way numerators and denominators will tend to grow).

The experiments below should give an idea of what to expect for the integer input case.

size = 1000;
mat = Table[
RandomInteger[{-1000, 1000}], {i, 1, size}, {j, 1, size}];
b = Table[RandomInteger[{-1000, 1000}], {i, 1, size}];
Timing[ByteCount[LinearSolve[mat, b]]]

(* Out[503]= {6.234739, 3696200} *)

size = 2000;
mat = Table[
RandomInteger[{-1000, 1000}], {i, 1, size}, {j, 1, size}];
b = Table[RandomInteger[{-1000, 1000}], {i, 1, size}];
Timing[ByteCount[LinearSolve[mat, b]]]

(* Out[507]= {52.200251, 14560392} *)


size = 4000; mat = Table[ RandomInteger[{-1000, 1000}], {i, 1, size}, {j, 1, size}]; b = Table[RandomInteger[{-1000, 1000}], {i, 1, size}]; Timing[ByteCount[LinearSolve[mat, b]]]

(* Out[511]= {413.295154, 61888776} *)

And we have not even quite hit 7! in that last example. At 8! expect to have problems with memory as well as speed. Actually this is true for the numeric case as well. Really at that point you can only hope to use sparse methods on sparse matrices.

• Thank you for the information. I wonder if there is any place I can look up the order of computational complexity increase for any mathematica function? I don't think it is mentioned in the function library. – Kagaratsch Feb 25 '16 at 23:40
• No such of which I'm aware. Fixed precision linear algebra is fairly standard O(n^3) despite there being methods that are faster in theory. For integer matrices we now use methods based on prime fields followed by Hensel lifting, and complexity there is covered in the relevant literature. – Daniel Lichtblau Feb 26 '16 at 0:37