equation solving - Long Execution Time of Reduce

I have a set of linear homogeneous equations with 8580 variables. I want to estimate the time Mathematica takes to solve the system, by solving smaller sets of equations. For example, I can take three equations, apply Reduce and see the AbsoluteTiming. Then, I repeat the same thing with four, five, ten etc.. equations and I can fit the execution time of Reduce with a polynomial or exponential function.

However, the system is very big, only three equations store about 2MB in a text-file and indeed Reduce takes too much time when applied to only one equation (which instead should be very easy, since it should just resolve on one unknown ). Instead, Solve is pretty fast but I have some experiences that Solve does not find all the solutions, so I would like to rely on the algorithms used by Reduce. For example, you can find one equation here (I am forced to refer to external links, as the output of the system is a messy and wouldn't fit here) that you can save in a file and import with Get. Using Solve just takes 0.46 seconds, while Reduce takes more than 1minute (then, I aborted the command).

How can I speed up the execution time of Reduce in a clever way, when dealing with a lot of variables? I guess that a system of homogenous equations in 8580 unknowns is challenging to solve. Any suggestion?

EDIT

I have constructed an explicit example that anyone can run on his/her laptop.

f1[a_, b_] := Det[{\[Lambda][a], \[Lambda][b]}]
generate\[Lambda][] := Module[{},
Clear[\[Lambda]];
\[Lambda][a_] :=  \[Lambda][a] = 1/RandomInteger[{1, 4}] RandomInteger[{-30, 30}, 2];
Table[\[Lambda][a], {a, 1, 5}];]

func = Sum[Subscript[A, i, j, k, m, n, p] f[i, j] f[k, m] f[n, p], {i, 1, 5}, {j, 1, 5}, {k, 1, 5}, {m, 1, 5}, {n, 1, 5}, {p, 1, 5}];

sys = {};

For[i = 1, i <= 250, i++,
generate\[Lambda][];
AppendTo[sys, (func /. f -> f1)==0]]

The function func generate the equations of 8000 variables, which are stored in sys. The variables $A_{i,j,k,m,n,p}$ are the unknowns in terms of which I want to solve the system. For the system I generate, I get a 601486920 bytes system from ByteCount[sys]. What is the most efficient way to get a solution of this system?

Answer

Just too long for a comment.

This matrix encodes 261 homogeneous equations in 8580 variables.

A = RandomReal[{-1, 1}, {261, 8580}];

Its size is

UnitConvert[Quantity[N@ByteCount[A], "Byte"], "Megabytes"]

Quantity[17.9152, "Megabytes"]

Determining a basis for its nullspace:

nullspace = NullSpace[A]; // AbsoluteTiming // First

1.72552

Takes 1.7 seconds. For all equations at once.

Towards OP's edit that brought up the code for generating such a system of equations.

I am sorry to say this, but this is really one of the worst written pieces of code that I have ever seen. It took me more than 37 GB to run it. Actually I stopped it after 10 minutes or so because I had no hope that it will terminate in the near future.

There are several things that are done in a very wrong way. Most notably:

• 
  

  Using AppendTo for building the list sys where Table would have been entirely sufficient: Each time you append, the whole list has to be copied -- and because there is so large data in it, it takes forever.

  
  
  


• 
  

  Creating aboundant, large symbolic expressions like in func and to use ReplaceAll on it hundreds of times.

  
  


• 
  

  Recomputing numbers over and over again (the results of f[i,j] can be recycled!).

  
  

On the other hand, the coefficient matrix A of your homogeneous system can be computed in machine precision as follows within 10 ms (milliseconds):

First we need a CompiledFunction for the numbercrunching:

cf = With[{Part = Compile`GetElement},
   Compile[{{λ, _Real, 2}},
    Block[{f, mm},
     mm = Length[λ];
     f = Table[λ[[i, 1]] λ[[j, 2]] - λ[[i, 2]] λ[[j, 1]], {i, 1, mm}, {j, 1, mm}];
     Flatten@Table[
       f[[i, j]] f[[k, m]] f[[n, p]], {i, 1, mm}, {j, 1, mm}, {k, 1, mm}, {m, 1, mm}, {n, 1, mm}, {p, 1, mm}
       ]
     ],
    CompilationTarget -> "C",

    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];

Next, we generate all random λ at once and feed them to cf:

A = cf[
     Divide[
      N[RandomInteger[{-30, 30}, {250, 5, 2}]],

      N[RandomInteger[{1, 4}, {250, 5}]]
      ]
     ]; // AbsoluteTiming // First

0.009569

Now we can compute the nullspace:

nullspace = NullSpace[A]; // AbsoluteTiming // First

4.10173

Takes only about 4 seconds.

So, if it is about efficiency, one has also to think about efficient storage formats for the equations. And machine precision matrices are the canonical way of storing linear equations.