# Fast pairwise difference

I'm trying to construct a matrix $$A_{ij}=X_i-Y_j$$ where $$X$$ and $$Y$$ are vectors with thousands of real numbers. The fastest and closest thing I found is DistanceMatrix[X, Y], 0.4 seconds for 20000x10000 matrix, but the result only valid if $$X_i>Y_j$$, otherwise, the sign should be flipped. There is a way to flip sign relatively fast with LowerTriangularize (roughly takes another 0.5 seconds), but it works only if $$X$$ and $$Y$$ are of the same size. Maybe there is a way to calculate $$A_{ij}$$ as fast as DistanceMatrix does its thing? Or is there a way to quickly flip sign for matrix of arbitrary height and width? Any cycles work very slow.

Running on an M1 Max (essentially an 8 core machine). Not only does Compile compile faster than FunctionCompile; it also allows parallelization via OpenMP:

f = FunctionCompile[Function[
{
Typed[list1, "PackedArray"::["Integer64", 1]],
Typed[list2, "PackedArray"::["Integer64", 1]]
},
Table[i - list2, {i, list1}]],
UseEmbeddedLibrary -> True
]; // AbsoluteTiming


{1.6023, Null}

cf = Compile[{{a, _Integer}, {b, _Integer, 1}},
Table[a - CompileGetElement[b, j], {j, 1, Length[b]}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
]; // AbsoluteTiming


{0.306628, Null}

m = 20000;
n = 10000;
a = RandomInteger[{1, 10}, m];
b = RandomInteger[{1, 10}, n];

A = f[a, b]; // RepeatedTiming
B = cf[a, b]; // RepeatedTiming

A == B


{0.225791, Null}

{0.0697286, Null}

True

Alas, you need to have a C compiler installed on your system.

• I am stupid, I upvoted and forgot to ask: do you have any good recommendations for C-compiler for mac? With the apple silicon I mean. Not sure if you would like to include it as part of the reply here or maybe leave a comment in the chat.
– bmf
Commented Feb 24, 2023 at 12:33
• I typically use Apple Clang shipped with XCode. I have never had any serious problems with it. Good alternatives are, of course, the (non-Apple) clang and gcc. You can install both of them with homebrew. All of them produce good code with similar performance. I use Apple Clang because (i) it is nicely integrated in the system and because (ii) it is generally a good idea to use the compiler produced by the hardware manufacturer unless you have good reasons not to do so. Commented Feb 24, 2023 at 14:38
• Thanks for this. I will make sure to give it a try :-)
– bmf
Commented Feb 24, 2023 at 14:40
• Holly, that's fast Commented Feb 25, 2023 at 8:40
• Unfortunately ParallelTable is still not supported in FunctionCompile... Commented Jul 6, 2023 at 5:17

Two more ways are suggested below. The output is the same as in @Syed's answer.

alist = RandomInteger[{1, 10}, 10];
blist = RandomInteger[{1, 10}, 4];


Then

😎= # - blist &;
😱 = Map[Subtract[#, blist] &];


and we do

😎/@ alist // MatrixForm
😱@alist // MatrixForm


Edit 1: taking the comment by @ Joshua Schrier into consideration, this is even faster

foo = Compile[{{a, _Integer, 1}, {b, _Integer, 1}},
Map[# - b &, a]];
foo[alist, blist]


Edit 2: taking the comment by @Ben Izd into consideration, for versions after 12 we can use

foo2 = FunctionCompile[
Function[{Typed[list1, "PackedArray"::["Integer64", 1]],
Typed[list2, "PackedArray"::["Integer64", 1]]},
Table[i - list2, {i, list1}]]];
foo2[alist, blist] // MatrixForm


Edit 3:

Comparing the RepeatedTimings of the various approaches so far.

$Version  With the lists alist = RandomInteger[{1, 10}, 20000]; blist = RandomInteger[{1, 10}, 10000];  we compare the following syed[l1_List, l2_List] := Outer[Subtract, l1, l2](*taking the comment by @Roman into account*) 😎 = # - blist &; 😱 = Map[Subtract[#, blist] &]; foo = Compile[{{a, _Integer, 1}, {b, _Integer, 1}}, Map[# - b &, a]]; foo2 = FunctionCompile[ Function[{Typed[list1, "PackedArray"::["Integer64", 1]], Typed[list2, "PackedArray"::["Integer64", 1]]}, Table[i - list2, {i, list1}]]];  and now we measure syed[alist, blist]; // RepeatedTiming 😎 /@ alist; // RepeatedTiming 😱@alist; // RepeatedTiming foo[alist, blist]; // RepeatedTiming foo2[alist, blist]; // RepeatedTiming  • Thanks a lot, that works just as fast as DistanceMatrix. Commented Feb 24, 2023 at 1:13 • @VsevolodA. glad I was able to help :-) – bmf Commented Feb 24, 2023 at 1:15 • You can shave off an extra 2-3x time by:  foo = Compile[{{a, _Integer, 1}, {b, _Integer, 1}}, Map[# - b &, a]] foo[alist, blist]  (change to _Real if you want reals, etc.) Commented Feb 24, 2023 at 1:51 • If you're using version 12.0+, you can halve the timing of Compile by using FunctionCompile: FunctionCompile[ Function[{Typed[list1, "PackedArray"::["Integer64", 1]], Typed[list2, "PackedArray"::["Integer64", 1]]}, Table[i - list2, {i, list1}] ]] Commented Feb 24, 2023 at 5:29 • I would rather use cf =Compile[{{a, _Integer, 1}, {b, _Integer}}, Table[CompileGetElement[a, i] - b, {i, 1, Length[a]}], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> "Speed" ]... Commented Feb 24, 2023 at 7:10 As Henrik has demonstrated in his excellent answer, the good old C compiler interface is still super fast -- the fastest if considering the compiling time. But the new compiler (FunctionCompile) still stands some chance. We see the high level way FunctionCompile[..., Table[...], ...] performs not so well. Intuitively, using ParallelTable instead should have accelerated it, except unfortunately that's not supported yet. So we'll have to go a bit lower level and endure some cumbersome code: fparal = FunctionCompile[Function[ { Typed[a, "PackedArray"["Integer64", 1]] , Typed[b, "PackedArray"["Integer64", 1]] } , Module[{ca, cb, c, cc, la = Length@a, lb = Length@b} , c = ArrayNewPackedArray[TypeSpecifier["Integer64"], LiteralType[2], {la, lb}] ; ca = ArrayGetData[a] ; cb = ArrayGetData[b] ; cc = ArrayGetData[c] ; ParallelParallelDo[ (* <- might use OpenMPParallelDo instead *) Do[ ToRawPointer[cc, i lb + j, FromRawPointer[ca, i] - FromRawPointer[cb, j]] , {j, 0, lb - 1}] , {i, 0, la - 1}] ; c] ]]  20000*10000 is a bit too much for my RAM, so I'll test on smaller example (here cf and f are directly taken from Henrik's answer): $Version
(* Out[]= 13.3.0 for Microsoft Windows (64-bit) (June 1, 2023) *)

m = 2000;
n = 1000;
a = RandomInteger[{1, 10}, m];
b = RandomInteger[{1, 10}, n];

cf[a, b] == f[a, b] == fparal[a, b]
(* Out[]= True *)

RepeatedTiming[     cf[a, b]; , 5]
RepeatedTiming[      f[a, b]; , 5]
RepeatedTiming[ fparal[a, b]; , 5]
(* Out[]=
{0.00354544, Null}
{0.00818736, Null}
{0.00342079, Null}
*)


So fparal performs as well as cf on my machine with an Intel i7 CPU. It would be interesting to see how fparal performs on Apple Silicon though.

• Jepp, I can confirm that fparal has pretty much the same timing as cf (the Compile version) on my Apple M1 Max. God job! Commented Jul 6, 2023 at 5:58
• I'd like to mention that the task at hand is very memory bound which is why not much speed-up can be seen (typical already few CPU cores fully saturate the memory bandwith). Would be cool to see a comparison also for some more compute-intense task... Commented Jul 6, 2023 at 6:01
• @HenrikSchumacher Thanks for your valuable comments! I would surely like to see how the new compiler competes in various compute-intense tasks. Unfortunately I don't have access to cutting-edge powerful computers for that kind of test. Frankly, on my daily work I still mainly use Compile because of the incompatible compiling time, also because the still-lack-of-docs status of FunctionCompile (like I have to search the source code to guess the usage of NewPackedArray). Commented Jul 6, 2023 at 6:16
• All fine. I myself also always use Compile and will continue to do so until FunctionCompile is more stable and better documented. =) "like I have to search the source code to guess the usage of NewPackedArray" Very good. Thank you for you effort! Commented Jul 6, 2023 at 6:28

Using Outer:

alist = RandomInteger[{1, 10}, 10];
blist = RandomInteger[{1, 10}, 4];

{alist, blist}


{{2, 5, 1, 8, 1, 1, 9, 7, 1, 5}, {2, 9, 6, 2}}

(res = Outer[(#1 - #2 &), alist, blist]) // MatrixForm


$$\left( \begin{array}{cccc} 0 & -7 & -4 & 0 \\ 3 & -4 & -1 & 3 \\ -1 & -8 & -5 & -1 \\ 6 & -1 & 2 & 6 \\ -1 & -8 & -5 & -1 \\ -1 & -8 & -5 & -1 \\ 7 & 0 & 3 & 7 \\ 5 & -2 & 1 & 5 \\ -1 & -8 & -5 & -1 \\ 3 & -4 & -1 & 3 \\ \end{array} \right)$$

• The issue is it will about 10 times slower than DistanceMatrix. Commented Feb 24, 2023 at 1:08
• Outer[Subtract, alist, blist] is easier on the eye. Commented Feb 24, 2023 at 6:38
• A request could be put in for SignedDistanceMatrix.
– Syed
Commented Feb 24, 2023 at 8:27
• @Syed signed distance only exists for 1d elements, while DistanceMatrix works for any vectors. Commented Feb 26, 2023 at 0:33

It's worth noting that Outer[Plus, v1, -v2] is much faster than Outer[Subtract, v1, v2], it's not the fastest, but very simple.

m = 20000;
n = 10000;
a = RandomInteger[{1, 10}, m];
b = RandomInteger[{1, 10}, n];

r1 = Outer[Plus, a, -b]; // AbsoluteTiming
r2 = Outer[Subtract, a, b]; // AbsoluteTiming
r1===r2


{0.614791, Null}
{26.8156, Null}
True

Henrik's cf on my computer

cf[a, b]; // RepeatedTiming


{0.218258, Null}

• (+1) This timing difference is so weird that it should be filed as bug report... Commented Jul 6, 2023 at 15:21