Specialized form of Outer[f, list, list] for symmetric f

Question

Consider a function f that is symmetric with regard to its two arguments i.e. f[x,y] == f[y,x]. Using Outer one can for instance build a correspondence matrix that computes some distance measure f between each two points in a list using

Outer[f, list, list]

This however computes every entry f[i, j] literally twice (since f[i, j] == f[j, i]) which is unnecessary. Is there a way to beat Outer (which is quite fast already) using this knowledge?

I tried

Table[Plus[list[[i]], list[[j]]], {i, 1, Length[list]}, {j, i, Length[list]}]

which only computes each entry once but is orders of magnitude slower even for simple examples like f = Plus; list = Range[1000].

It makes no sense to have an iterator {I, list}, since list is a List and not an integer. Perhaps you mean Length[list]. If you want to use Table, then try: Table[f[I,j], {I, Length[list},{I,j}] to compute the unique values, and then fill in the rest with a simple Table call. — David G. Stork, Jan 08 '18 at 06:53
@DavidG.Stork Table[expr, {i, {i1, i2, i3, ...}}] and thus Table[expr_, {i, someList}] is a valid syntax variant for Table. Its the 5th entry in the syntax overview in the documentation. — Sascha, Jan 08 '18 at 06:56
@DavidG.Stork Never mind, I only tried the Table variant with input like list=Range[10] and for this the Part specification works out. — Sascha, Jan 08 '18 at 07:00
Plus is special in a number of way; and Outer[Plus,...] is a special case. It might not be a good test example? — Michael E2, Jan 08 '18 at 21:11
When working with symmetric arrays, it may be useful to consider SymmetrizedArray, particularly in cases with high symmetry (with deep arrays). You may use SymmetrizedArray as constructor, or you can use standard array methods and then use Symmetrize to produce SymmetrizedArray results. This may also help in case in which storage space is important, because SymmetrizedArray only keeps one copy of each independent element. — jose, Jan 09 '18 at 18:36

score 6 · Answer 1 · answered Jan 08 '18 at 21:48

Outer[Plus,...] is a special case that is very fast. It's not a good example to compare with Table. I'm not sure what assumptions one might make about f, which might lead to particular optimizations, but I'll pick Norm. It's neither Listable nor compilable. The half Table saves some time over Outer:

list = Range[1., 1000.];

m1 = Outer[Norm@*List, list, list]; // AbsoluteTiming
(*  {3.31615, Null}  *)

Table[Norm[{list[[i]], list[[j]]}], {i, 1, Length[list]}, {j, i, 
    Length[list]}]; // AbsoluteTiming
(*  {2.18397, Null}  *)

Here are a couple of ways to fill out the rest of the array inspired by Best way to create symmetric matrices and Generating random symmetric matrix:

(mat = PadLeft[
     Table[
      Norm[{list[[i]], list[[j]]}], {i, 1, Length[list]}, {j, i, 
       Length[list]}],
     Automatic, 0.];
   m2 = mat + Transpose[UpperTriangularize[mat, 1]]); // AbsoluteTiming
(*  {2.32359, Null}  *)

(upper = Flatten[
     Table[
      Norm[{list[[i]], list[[j]]}], {i, 1, Length[list]}, {j, i + 1, Length[list]}]];
   diag = Table[Norm[{x, x}], {x, list}];
   m3 = Statistics`Library`VectorToSymmetricMatrix[upper, diag, Length@list]
); // AbsoluteTiming
(*  {2.22123, Null}  *)

m1 == m2 == m3
(*  True  *)

You can also use Sqrt@Outer[Plus, #, #] &[list^2] – chyanog Feb 23 '23 at 07:56 — chyanog, Feb 23 '23 at 07:56

score 5 · Answer 2 · answered Jan 08 '18 at 08:43

For such a simple thing like Plus, I can beat Outer on my quad core CPU by a factor of 2 with

cf1 = With[{Part = Compile`GetElement},
   Compile[{{a, _Real}, {b, _Real, 1}},
    Table[a + b[[i]], {i, 1, Length[b]}],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];

We can do this for the upper triangular part of the output matrix in a similar fashion with

cf2 = With[{Part = Compile`GetElement},
   Compile[{{a, _Real}, {b, _Real, 1}, {i, _Integer}},
    Table[a + b[[j]], {j, i, Length[b]}],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ];

Here are the timings:

n = 10000;
list = RandomReal[{-1, 1}, {n}];
r0 = Outer[Plus, list, list]; // RepeatedTiming // First
r1 = cf1[list, list]; // RepeatedTiming // First
r2 = cf2[list, list, Range[n]]; // RepeatedTiming // First
r0 == r1

0.513

0.20

0.11

True

Indeed, cf2 is faster than cf1 by almost a factor of 2 (as expected). However, note that cf2 produces a ragged list of lists and not a square array. Getting this back into a square matrix is not that inexpensive. So it may be still be a good idea to compute some entries of the output twice. Note that often (and in particular in this case) it is not the floating point computation that limits the execution speed; it is memory access.

score 0 · Answer 3 · answered Jan 08 '18 at 18:55

The very first lines of the question bring to mind the Orderless attribute.

"...Consider a function f that is symmetric with regard to its two arguments i.e. f[x,y] == f[y,x]"

One possible way to overcome the problem reported below

"...This however computes every entry f[i, j] literally twice (since f[i, j] == f[j, i])"

might be to memoize the function f

ClearAll[f,g]
SetAttributes[f, Orderless]
f[x__] := f[x] = g[x]

where g[] stands for an appropriate expression for f.

Consider the following example:

i = 1;
g[x__] := (Print[i]; i += 1; h[x])
list = {1, 2, 3}
Outer[g, list, list, 1]

Evaluation of that, will print numbers 1 through 9 and then return

{{h[1, 1], h[1, 2], h[1, 3]}, {h[2, 1], h[2, 2], h[2, 3]}, {h[3, 1], h[3, 2], h[3, 3]}}

If, on the other hand, we were to evaluate

i = 1;
Outer[f, list, list, 1]

then, the output we'd receive, would be

{{h[1, 1], h[1, 2], h[1, 3]}, {h[1, 2], h[2, 2], h[2, 3]}, {h[1, 3], h[2, 3], h[3, 3]}}

after printing numbers 1 through 6.

Please note that if you evaluate this last example code once more, nothing will get printed and the result will be the same as it was after the first evaluation.

score 0 · Answer 4 · answered Mar 15 '18 at 14:04

Without any attribute Outer[g, {1, 2, 3}, {1, 2, 3}] yields

 {{h[1, 1], h[1, 2], h[1, 3]}, {h[2, 1], h[2, 2], h[2, 3]}, {h[3, 1], h[3, 2], h[3, 3]}}

Using a pure function with the Orderless attribute:

g = Function[{x, y}, h[x, y], Orderless]

we get

Outer[g, {1, 2, 3}, {1, 2, 3}]

{{h[1, 1], h[1, 2], h[1, 3]}, {h[1, 2], h[2, 2], h[2, 3]}, {h[1, 3], h[2, 3], h[3, 3]}}

Hence, if the code for h is efficient, so is the application of Outer.

Specialized form of Outer[f, list, list] for symmetric f

4 Answers4