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?

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

3 Answers 3


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

Mathematica graphics

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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Sep 15, 2016 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 reciprocal distance matrix

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)]]
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Oct 6, 2016 at 11:09

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.