Calculating the reciprocal distance matrix without inflicting `ComplexInfinity`

Given a list of coordinates r (r[[i]]!=r[[j]]), I’d like to know the reciprocals of distances of all pairs in the list, and for the convenience subsequent operations, the trace of the resulting matrix should be all zero. I feel that this should be a frequent need, but I can’t do it optimally.

My code:

R = Outer[Norm, r, r, 1];
rR = Quiet[1/R] /. {ComplexInfinity -> 0.}

But this is not such a good idea as ReplaceAll is significantly slower than the other calculations in this code. Is it a good idea to use For or Table and loop over all indices, or is there a better way to do this?




If you know all of the r values are different, why not just set the diagonal to 0 afterward? Do UpperTriangularize[#, 1] + LowerTriangularize[#, -1] &@Quiet[1/R]. See here.
– march
Sep 6 at 5:41



Or what’s your definition of distances of all pairs?
– goodbye_M.SE
Sep 6 at 14:07



I think Outer[Norm, r, r, 1] should be Outer[EuclideanDistance, r, r, 1]?
– xzczd
Sep 6 at 14:41


3 Answers


If the coordinates are machine reals and speed is an issue, I would Compile a function:

reciprocalDist = Compile[{{r, _Real, 2}},
R = Outer[Subtract, r, r, 1]; (* Outer[Norm,..] is not supported in Compile *)
Map[If[# == 0, 0, 1/#] &@Norm[#] &, R, {2}]

SeedRandom[0]; (* for reproducibility *)
reciprocalDist[RandomReal[{-1, 1}, {4, 3}]] // MatrixForm

Theoretically one could cut the speed in half by calculating the upper triangular part. One could preallocate a zero matrix and fill the distances two at a time with a Do loop, for instance.



Why DV? I can’t see what’s wrong, other than the stated limitations. Since no comment, obviously a cowardly troll not sincerely interested in improving the site.
– Michael E2
Sep 15 at 23:46

Since DistanceMatrix[] is built-in, it seems natural to use it for this problem. Using Michael’s example, let me present two approaches:

BlockRandom[SeedRandom[0]; (*for reproducibility*)
pts = RandomReal[{-1, 1}, {4, 3}]];

(* method 1 *)
With[{id = IdentityMatrix[Length[pts]]}, 1/(DistanceMatrix[pts] + id) – id]

(* method 2 *)
DistanceMatrix[pts, DistanceFunction -> (With[{d = EuclideanDistance[##]},
If[d == 0, 0, 1/d]] &)]

Both should return

reci[R_] :=
With[{l = Length@R}, {m = SparseArray[Band[{1, 1}] -> 1, {l, l}]}, (1 – m)/(R + m)]

If you’re before v10.4:

reci[R_] :=
With[{l = Length@R},
With[{m = SparseArray[Band[{1, 1}] -> 1, {l, l}]}, (1 – m)/(R + m)]]



To the downvoter, I am interested in what was missing from my answer. Did I give too little detail? Or, was it considered incorrect in some way? Either way, would you please elaborate? I’m not trying to complain here, I’m just curious about what I could have done instead.
– xzczd
Oct 6 at 11:09