Trying to speed up finding solutions of quadratic form related to the weight of a matrix

Question

Say I have an integer matrix of size $5 \times 5$ - I am trying to find the solutions to the equation $n^4 wtM = n w^2 + x1^2 + x2^2 +x3^2 +x4^2 + x5^2$ subject to the constraints $x1+x2+x3+x4+x5 = 0$ and $w>0$ where, in this case $n=5$ and $wtM$ is the weight of the matrix i.e. $n^4 wtM$ is just some integer - typically with values in the range one thousand to five thousand.

Now I am using the Mathematica function PowersRepresentations to solve $n^4 wtM - n w^2$ for w ranging from 1 until $n^4 wtM - n w^2$ is still positive.

This gives me a list of lists $w,x1,x2,x3,x4,x5$ for each w. I then combine all these solutions together which I call $reps$ in the following code.

Now PowersRepresentations only gives positive values. I need to allow the $x1,x2,x3,x4,x5$ to be positive or negative subject to the constraint $x1+x2+x3+x4+x5 = 0$

So I take each list $w,x1,x2,x3,x4,x5$ and multiply it by vectors of permutations of plus and minus 1 (always having the first element in the vector 1 so as to keep w positive). I don't bother with vectors of all 1's or all -1's as this would automatically violate my summation constraint.

I then select the cases using a compiled function $cf1$ that makes sure I still obey the constraint $x1+x2+x3+x4+x5 = 0$. I also have some other constraints that I apply (via some other compiled functions) such as $x1\le x2 \le x3 \le x4$ and $x1 \ge 0 \space x2 \ge 0 \space x3 \ge 0 \space x4 \ge 0$.

Now the code I have written tries to make use of parallelisation but for anything above $7 \times 7$ takes a very long time to run.

I have perused many articles on here but I am still unsure of the various methods I have employed. The code works but I am concerned I am missing some tricks that would allow much faster execution. I have used Mathematica quite a bit in the past but primarily in using NDSolve and the like but not parallelising and compiling to C.

Any recommendations and pointing out speed increases (and memory usage reduction) would be much appreciated. As I said the code takes a long time and uses a lot of memory for $n \ge 7$. Notice I am not even saving the solutions but just counting them.

The example used below was for $n = 7$ and $wtM = 54/49$ giving $n^4 wtM = 2646$

One that runs in about 12 seconds is $n = 6$ and $wtM = 4/3$ giving $n^4 wtM = 1728$

n = 7;
wtM = 54/49;


progressreporting = True

ClearAll[results1, results2, results3, results4, data1, data2, data3, 
  data4];

ordering1 = "";
Do[ordering1 = 
  StringJoin[ordering1, "#[[" <> ToString[i] <> "]]+"], {i, 2, 
  n + 1}]; ordering1 = 
 ToExpression[StringJoin[StringDrop[ordering1, -1], "==0 &"]];

cf1 = With[{ordering = ordering1}, 
   Compile[{{triples, _Integer, 2}}, 
    DeleteDuplicates[Select[triples, ordering]], 
    RuntimeAttributes -> {Listable}, Parallelization -> True, 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"]];

ordering2 = "";
Do[ordering2 = 
   StringJoin[ordering2, "#[[" <> ToString[i] <> "]]<= "], {i, 2, 
   n}];
ordering2 = 
  ToExpression[StringJoin[StringDrop[ordering2, -3], " &"]];
If[n == 2, ordering2 = ordering1];

cf2 = With[{ordering = ordering2}, 
   Compile[{{triples, _Integer, 2}}, Select[triples, ordering], 
    RuntimeAttributes -> {Listable}, Parallelization -> True, 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"]];

ordering3 = "";
Do[ordering3 = 
  StringJoin[ordering3, "#[[" <> ToString[i] <> "]]>=0 &&"], {i, 2, 
  n}]; ordering3 = 
 ToExpression[StringJoin[StringDrop[ordering3, -2], " &"]];

cf3 = With[{ordering = ordering3}, 
   Compile[{{triples, _Integer, 2}}, Select[triples, ordering], 
    RuntimeAttributes -> {Listable}, Parallelization -> True, 
    CompilationTarget -> "C", RuntimeOptions -> "Speed"]];

DistributeDefinitions[cf1, cf2, cf3];
AbsoluteTiming[
 reps = Join @@ 
   Table[
    Map[Prepend[#, i] &, 
     PowersRepresentations[n^4 wtM - n i^2, n, 2]], {i, 1, 
     IntegerPart[Sqrt[(n^4 wtM)/n]]}];
 vc = Map[Insert[#, 1, 1] &, Take[Tuples[{-1, 1}, n], {2, 2^n - 1}]];
 
 ParallelEvaluate[foo1 = {}]; ParallelEvaluate[foo2 = {}]; 
 ParallelEvaluate[foo3 = {}]; ParallelEvaluate[foo4 = {}];
 
 sow[x_] := (foo1 = {foo1, x};); sow2[x_] := (foo2 = {foo2, x};); 
 sow3[x_] := (foo3 = {foo3, x};); sow4[x_] := (foo4 = {foo4, x};);
 
 ParallelDo[
  
  temp = Extract[reps[[j]], 1]; temp2 = Drop[reps[[j]], 1];
  test0 = Permutations[temp2];
  test1 = Map[Prepend[#, temp] &, test0];
  lst = {}; 
  lst = Table[
    Map[Times[#, test1[[i]]] &, vc], {i, 1, Length[test1]}];
  test2 = Join @@ lst;
  
  test3 = cf1[test2];
  test32 = cf2[test3];
  test33 = cf3[test3];
  test34 = cf3[test32];
  
  sow[Length@test3]; sow2[Length@test32]; sow3[Length@test33]; 
  sow4[Length@test34];
  
  , {j, 1, Length[reps]}, 
  ProgressReporting -> progressreporting,
  Method -> "ItemsPerEvaluation" -> 1000];
 
 results1 = Join @@ ParallelEvaluate[Flatten@foo1];
 results2 = Join @@ ParallelEvaluate[Flatten@foo2];
 results3 = Join @@ ParallelEvaluate[Flatten@foo3];
 results4 = Join @@ ParallelEvaluate[Flatten@foo4];
 
Print["\!\(\*SuperscriptBox[\(n\), \(4\)]\)wtM = ", n^4 wtM];
  
  data1 = Total[results1];
  Print[data1, " solutions in Z with no ordering"];
  If[n > 2, data2 = Total[results2];
    Print[data2, " solutions in Z with ordering"]];
  data3 = Total[results3];
  Print[data3, 
    " solutions in \!\(\*SubscriptBox[\(N\), \(0\)]\) with no \
ordering"];
  If[n > 2, data4 = Total[results4];
    Print[data4, 
     " solutions in \!\(\*SubscriptBox[\(N\), \(0\)]\) with \
ordering"]];] // First

n^4wtM = 2646

342611822 solutions in Z with no ordering

584282 solutions in Z with ordering

169915 solutions in Subscript[N, 0] with no ordering

556 solutions in Subscript[N, 0] with ordering

1688.18

Answer 1

It's not an answer but a long comment.

Without understanding your code, these are some improvements just by using built-in functions:

• ordering1
  You'd used string concatenation to build your expression, which could be improved in both readability and performance?

(* use *)
Plus @@ Take[#, {2, n+1}] == 0 &
(* or *)
Total[Take[#, {2, n+1}]] == 0 &

(* instead of *)

ordering1 = "";
Do[ordering1 = 
  StringJoin[ordering1, "#[[" <> ToString[i] <> "]]+"], {i, 2, 
  n + 1}]; ordering1 = 
 ToExpression[StringJoin[StringDrop[ordering1, -1], "==0 &"]];

(* example: #[[2]]+#[[3]]+#[[4]]+#[[5]]+#[[6]]==0 & *)

• ordering2
  like ordering1, just use the power of built-ins

OrderedQ[Take[#, {2, n}]]&

(* instead of *)

ordering2 = "";
Do[ordering2 = 
   StringJoin[ordering2, "#[[" <> ToString[i] <> "]]<= "], {i, 2, 
   n}];
ordering2 = 
  ToExpression[StringJoin[StringDrop[ordering2, -3], " &"]];

(* example: #[[2]]<= #[[3]]<= #[[4]]<= #[[5]]<= #[[6]] & *)

• ordering3
AllTrue needs version 10

AllTrue[NonNegative[Take[#, {2, n}]], TrueQ] &;

(* instead of *)

ordering3 = "";
Do[ordering3 = 
  StringJoin[ordering3, "#[[" <> ToString[i] <> "]]>=0 &&"], {i, 2, 
  n}]; ordering3 = 
 ToExpression[StringJoin[StringDrop[ordering3, -2], " &"]];

(* example: #1[[2]] >= 0 && #1[[3]] >= 0 && #1[[4]] >= 0 && #1[[5]] >= 0 && #1[[6]] >= 0 & *)

• lst

Outer[Times, test1, vc]

(* instead of *)

Table[Map[Times[#, test1[[i]]] &, vc], {i, 1, Length[test1]}]

• Also, you can use PadLeft instead of Map[Prepend.../Map[Insert[...,1,1] ...