Calculating the reciprocal distance matrix without inflicting `ComplexInfinity`

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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