Fast pairwise difference

Question

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.

Answer 1

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 - Compile`GetElement[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.

Answer 2

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

Answer 3

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 = Array`NewPackedArray[TypeSpecifier["Integer64"], LiteralType[2], {la, lb}]
          ; ca = Array`GetData[a]
          ; cb = Array`GetData[b]
          ; cc = Array`GetData[c]
          ; Parallel`ParallelDo[ (* <- might use OpenMP`ParallelDo 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.

Answer 4

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)$$

Answer 5

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}