10
$\begingroup$

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.

$\endgroup$

5 Answers 5

10
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – bmf
    Feb 24, 2023 at 12:33
  • 3
    $\begingroup$ 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. $\endgroup$ Feb 24, 2023 at 14:38
  • $\begingroup$ Thanks for this. I will make sure to give it a try :-) $\endgroup$
    – bmf
    Feb 24, 2023 at 14:40
  • $\begingroup$ Holly, that's fast $\endgroup$ Feb 25, 2023 at 8:40
  • $\begingroup$ Unfortunately ParallelTable is still not supported in FunctionCompile... $\endgroup$
    – Silvia
    Jul 6, 2023 at 5:17
13
$\begingroup$

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

list

Edit 3:

Comparing the RepeatedTimings of the various approaches so far.

$Version

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

timings

$\endgroup$
8
  • 1
    $\begingroup$ Thanks a lot, that works just as fast as DistanceMatrix. $\endgroup$ Feb 24, 2023 at 1:13
  • $\begingroup$ @VsevolodA. glad I was able to help :-) $\endgroup$
    – bmf
    Feb 24, 2023 at 1:15
  • 4
    $\begingroup$ 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.) $\endgroup$ Feb 24, 2023 at 1:51
  • 2
    $\begingroup$ 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}] ]] $\endgroup$
    – Ben Izd
    Feb 24, 2023 at 5:29
  • 3
    $\begingroup$ I would rather use cf =Compile[{{a, _Integer, 1}, {b, _Integer}}, Table[Compile`GetElement[a, i] - b, {i, 1, Length[a]}], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True, RuntimeOptions -> "Speed" ]... $\endgroup$ Feb 24, 2023 at 7:10
7
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Jepp, I can confirm that fparal has pretty much the same timing as cf (the Compile version) on my Apple M1 Max. God job! $\endgroup$ Jul 6, 2023 at 5:58
  • 1
    $\begingroup$ 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... $\endgroup$ Jul 6, 2023 at 6:01
  • 1
    $\begingroup$ @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). $\endgroup$
    – Silvia
    Jul 6, 2023 at 6:16
  • $\begingroup$ 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! $\endgroup$ Jul 6, 2023 at 6:28
6
$\begingroup$

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

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

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}

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

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.