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

2

I think Outer[Norm, r, r, 1] should be Outer[EuclideanDistance, r, r, 1]?

– xzczd

Sep 6 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}},

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.

1

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

1

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