# 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 '16 at 5:41
• Or what's your definition of distances of all pairs? – xyz Sep 6 '16 at 14:07
• I think Outer[Norm, r, r, 1] should be Outer[EuclideanDistance, r, r, 1]? – xzczd Sep 6 '16 at 14:41

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

reciprocalDist = Compile[{{r, _Real, 2}},
Module[{R},
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 '16 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 '16 at 11:09