0
$\begingroup$

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


$\endgroup$
0
$\begingroup$

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] ...
$\endgroup$
10
  • $\begingroup$ Thank you for the reply. Let me take a while to digest these. Much appreciated. This is the kind of thing I haven't done in Mathematica before and the learning curve is steep. - FYI I am using version 13.0 $\endgroup$
    – 1729taxi
    Jan 4 at 11:34
  • $\begingroup$ I found that your replacement for ordering1 slows my code down by a factor of 2. The replacement for ordering2 gives the following error ("Compilation of \ LessEqual@@Take[Compile`FunctionVariable$184772,{2,n}] is not \ supported for the function argument LessEqual.") And the replacement for ordering3 slows me down by at least a factor of 3. I am puzzled as to the slowdowns. $\endgroup$
    – 1729taxi
    Jan 4 at 12:04
  • $\begingroup$ Hi again. Your suggestions all improve readability (though I can't compile with Take as a function argument for LessEqual) - but the performance is poor. I lose abut a factor of two for each change including the Outer[Times,vc,test1,1,1]. For an n=6 I have been testing upon that runs in 25 seconds takes several minutes with your suggestions. I appreciate your input though - always good to see a different way of doing things. $\endgroup$
    – 1729taxi
    Jan 4 at 14:05
  • $\begingroup$ Thanks again but I think the problem from a performance standpoint is that Compile does not like these. I get slowdown factors of 2 or 3 for all of them. Also your edited orddering3 gives a Compile error to boot. I might try these without Compiling the functions but I think I'll get a big slowdown in ordering1 because of not compiling Select. $\endgroup$
    – 1729taxi
    Jan 4 at 14:32
  • $\begingroup$ @1729taxi Yes, you're right, the code is simplified but it seems we should cross the performance word. $\endgroup$
    – Ben Izd
    Jan 4 at 14:35

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.