Select seems unusually slow

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

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