16
$\begingroup$

I haven't been doing nearly as much Mathematica as I really should for the past few months, and in order to get my head back in the game I've been using it to run roughshod over easy Project Euler problems. In this case, one involved finding Pythagorean triples which totaled to 1000. The number of possibilities is small enough that there's no reason not to brute-force:

(triples = IntegerPartitions[1000, {3}]) // Length

(* 83333 *)

Remembering that IntegerPartitions lists the parts from largest to smallest, Select should work, and it does, but way more slowly than I expected!

Select[triples, 
  Apply[Function[{c, b, a}, a^2 + b^2 == c^2]]] // AbsoluteTiming
(* {0.257232, {{425, 375, 200}}} *)

In the grand scheme of things, a quarter-second isn't a long time, but it's still taking microseconds per list element, which seems like a lot. Switching to Cases is a little faster:

Cases[triples, {c_, a_, b_} /; c^2 == a^2 + b^2] // AbsoluteTiming
(* {0.144264, {{425, 375, 200}}} *)

Still, a seventh of a second is a long time with a modern computer. Heck, writing my own Do-based solution is faster:

Reap[
   Do[
    If[MatchQ[triple, {c_, a_, b_} /; c^2 == a^2 + b^2],
     Sow[triple]],
    {triple, triples}]][[-1, 1]] // AbsoluteTiming
(* {0.176861, {{425, 375, 200}}} *)

After a few minutes tweaking, I did come up with something a couple orders of magnitude faster, but it's certainly not the most perspicuous code I've ever written:

Pick[
 triples,
 {1, -1, -1}.Transpose[triples^2], 0]
(* {0.006377, {{425, 375, 200}}} *)

I assume this a matter of packing and unpacking arrays, but the whole thing is way off from my intuitions about Mathematica performance; I'd expect Select to be a bit slower than the optimized Pick-based solution, but nothing like this.

(Even transposing before squaring makes a substantial difference in performance.)

What's going on?

FWIW, this is with Mathematica 11.1.0 on a slightly elderly MacBook Pro.

$\endgroup$
  • 3
    $\begingroup$ With my BoolEval package you could do {a, b, c} = Transpose[triples]; BoolPick[triples, a^2 == b^2 + c^2]. It returns {{425, 375, 200}}. It takes about 0.0014 s on a 2014 MacBook Pro, as measured by RepeatedTiming in version 11.1. In spirit, it is very similar to your last solution. The role of the package is to provide a simple, human-readable notation. $\endgroup$ – Szabolcs Mar 23 '17 at 19:03
  • $\begingroup$ For comparison, your last solution takes 0.0054 seconds on my machine, which (surprisingly) seems slower than BoolEval. However, I just realized that {a, b, c} = Transpose[triples]; alone takes 0.013 seconds. What BoolEval transforms the computation to is 1 - Unitize[a^2 - b^2 - c^2], as can be seen with Clear[a, b, c]; BoolEval[a^2 == b^2 + c^2]. $\endgroup$ – Szabolcs Mar 23 '17 at 19:08
  • 2
    $\begingroup$ Sorry for so many comments ... another thing I just realized was that IntegerPartitions does not return a packed array. Manually packing it speeds up these types of vectorized calculations quite significantly. Now {a, b, c} = Transpose[triples]; is down to 0.0004 s, the BoolPick stays the same, and your Pick is down to 0.0014 s. The Select is up to 0.30 s from 0.25 s (due to unpacking). $\endgroup$ – Szabolcs Mar 23 '17 at 19:13
  • $\begingroup$ Another observation: I compared 10.4.1 with 11.0.1 and 11.1.0. Generally, things are either the same of faster in newer versions. Sometimes they are significantly faster. One exception is {a,b,c} = Transpose[triples]; with a 2x slowdown in 11.1.0 (compared to the other two versions) from 0.007 to about 0.014. $\endgroup$ – Szabolcs Mar 23 '17 at 19:27
  • 1
    $\begingroup$ On this transpose test version 9.0.1 only takes 0.0026 seconds. That's 5x faster than 11.1. If you're curious, the Select is about the same. (It had to be rewritten slightly for compatibility.) Also note that Select can be compiled, which makes it run in 0.012 s here. $\endgroup$ – Szabolcs Mar 23 '17 at 19:44
11
$\begingroup$

I do not know why Select is so slow, but I have observed similar performance problems in the past. I usually try to make use of one two methods to speed up the computation significantly.

But before going into those methods, let us observe that triples is not a packed array and pack it:

triples = Developer`ToPackedArray[triples];

This packing took about 0.0045 s on my machine, and will improve the performance of the below solutions.

Compile

Select is compilable. Let's compile it!

cf = Compile[{{triples, _Integer, 2}},
   Select[triples, #[[1]]^2 == #[[2]]^2 + #[[3]]^2 &]
   ];

cf[triples] // RepeatedTiming
(* {0.013, {{425, 375, 200}}} *)

Vectorization

Vectorization is doing arithmetic on whole arrays instead of individual elements. This way the whole array operation can be implemented in efficient, low-level code and take advantage of SIMD instructions and multiprocessing.

Select cannot use such techniques because it needs to invoke the Mathematica evaluator once for each element of triples.

What you are doing with your Pick[...] solution is very effective vectorization. Packing triples speeds it up further.

Pick[triples, {1, -1, -1}.Transpose[triples^2], 0] // RepeatedTiming
(* {0.0014, {{425, 375, 200}}} *)

For reference, the original ran in 0.0054 s on my machine.

Your solution can be improved slightly by getting rid of Transpose:

Pick[triples, (triples^2).{1, -1, -1}, 0] // RepeatedTiming
(* {0.0012, {{425, 375, 200}}} *)

As you note, the problem with such code is that it is difficult to read and difficult to write (without mistakes). My BoolEval package attempts to tackle this problem by providing a convenient notation.

We can solve this problem using BoolEval as follows:

{a, b, c} = Transpose[triples]; // RepeatedTiming
(* {0.00038, Null} *)

BoolPick[triples, a^2 == b^2 + c^2] // RepeatedTiming
(* {0.00059, {{425, 375, 200}}} *)

Behind the scenes, BoolEval just converts a^2 == b^2 + c^2 to fast array arithmetic. You can see this by running it on symbolic variables:

Clear[a, b, c];
BoolEval[a^2 == b^2 + c^2]
(* 1 - Unitize[a^2 - b^2 - c^2] *)
$\endgroup$
11
$\begingroup$
Solve[{a^2 + b^2 == c^2, a + b + c == 1000, a > 0, b > 0, c > 0, 
   a <= b}, {a, b, c}, Integers] // AbsoluteTiming

(*  {0.09877, {{a -> 200, b -> 375, c -> 425}}}  *)
$\endgroup$
  • 2
    $\begingroup$ I love that using Solve is faster. $\endgroup$ – Pillsy Mar 23 '17 at 19:33
0
$\begingroup$

If one ignores the time for compression also this works very fast at the selection process...

triples = Compress[triples];
Pick[triples, (triples^2).{1, -1, -1}, 0] // RepeatedTiming
$\endgroup$

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.