Select seems unusually slow

Question

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.

Answer 1

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] *)

Answer 2

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}}}  *)

Answer 3

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