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?


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[{},
\[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++,
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?

  • 1
    $\begingroup$ "homogeneous equations with 8580 variables" That's something you are supposed to pack into a matrix and to run Nullspace on it. (Hopefully the coefficients of your system are floating point numbers.) $\endgroup$ Jul 23 '19 at 17:32
  • $\begingroup$ @HenrikSchumacher Nullaspace does not seems to be effecient even for one equation, you can check it with the equation posted in privatebin.net/… $\endgroup$
    – apt45
    Jul 23 '19 at 21:38
  • $\begingroup$ Do you require exact solutions of your homogeneous equations or does machine precision suffice for you? $\endgroup$ Jul 23 '19 at 21:46
  • $\begingroup$ @HenrikSchumacher I really would like exact solutions. Then, if that's too computationally expensive, I can be happy with machine precision. Notice that If I export 261 equation in a text file, then I store more the 180MB of memory, you can download an example of such equations here dropbox.com/s/u8at3zu3h5074o8/system.dat?dl=0 $\endgroup$
    – apt45
    Jul 23 '19 at 21:49

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


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}];
       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[
      N[RandomInteger[{-30, 30}, {250, 5, 2}]],
      N[RandomInteger[{1, 4}, {250, 5}]]
     ]; // AbsoluteTiming // First


Now we can compute the nullspace:

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


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.

  • $\begingroup$ Thanks.. however this does not work so easily for my case.. try to download the system of equations I have uploaded here dropbox.com/s/u8at3zu3h5074o8/system.dat?dl=0 it takes too much time. $\endgroup$
    – apt45
    Jul 23 '19 at 21:58
  • $\begingroup$ The point I wanted to make: You stored your system in a really, really inefficient way. I do not even know how to retrieve that data from file! And honestly, I have not time to find it out. The data is not stored as a valid Mathematica expression! $\endgroup$ Jul 23 '19 at 22:03
  • $\begingroup$ Thank you for your answer. I have used the most naive command Export to export the system of equations. Maybe it is not the most efficient way, but I thought it was the most obvious one. $\endgroup$
    – apt45
    Jul 23 '19 at 22:06
  • 1
    $\begingroup$ When Mathematica translates the code into C, it adds certain bound checks whenever there is a read access to an array done with Part. The reasons is that the operation system might quit the copiled library and the Mathematica kernel as its parent process whenever it tries to read out of scope of the arrays (if it is a good operating system). So this is a robustness (and security) measure. But if one knows a priorily that the indexing does not go out of bounds, these checks just waste time. $\endgroup$ Jul 25 '19 at 17:04
  • 1
    $\begingroup$ One can deactivate the checks by using Compile`GetElement along with the option RuntimeOptions -> "Speed". Be careful though: This can lead to runtime errors that are incredibly hard to debug, in particular if you have complex programs with many calls to compiled libraries. See here for example. $\endgroup$ Jul 25 '19 at 17:07

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.