Improve speed with a double for loop

I am trying to calculate the Euclidean distance of pairs of points $$(x_i,y_i)$$ with $$x$$ and $$y$$ given as separate lists.

sqdiff = {};
For[i = 1, i < Length[x] - 1, i++,
{
For[j = i + 1, j < Length[x], j++,
{
AppendTo[sqdiff, (x[[i]] - x[[j]])^2 + (y[[i]] - y[[j]])^2];
}
]
}
]


Even for 800 points, this loops takes orders of magnitudes longer than an identical loop written in MATLAB. Are there any improvements I can make? Thanks.

• "Are there any improvements I can make? ". Yes. Don't program Mathematica as if it were Matlab. Use DistanceMatrix, vectorize it, etc. – ciao Mar 3 '20 at 21:00
• Even in MATLAB this were bad programming style. Typically, one is taught to use vectorization there as well. Due to the use of For/for and the to the successive AppendTo (cat in MATLAB). – Henrik Schumacher Mar 3 '20 at 21:05

Even in MATLAB, this would not be good programming style because successive concatenation is awfully slow. (And for is slow, too.) Better use Table:

n = 800;
x = RandomReal[{-1, 1}, n];
y = RandomReal[{-1, 1}, n];

Table[
With[{xi = x[[i]], yi = y[[i]]},
Table[
(xi - x[[j]])^2 + (yi - y[[j]])^2
, {j, 1, n}]
],
{i, 1, n}
]


Better:

Table[
Sqrt[(x[[i]] - x)^2 + (y[[i]] - y)^2],
{i, 1, n}
];


Even better: Use DistanceMatrix:

DistanceMatrix[Transpose[{x, y}]];

• 100% it is bad style, but it still was done in under 5 seconds in matlab with 500k points, so I was shocked. Thanks for the help! – Gregory Mar 3 '20 at 21:07
• With 500k? Impossible. That would make a $500k \times 500k$ matrix or an $(500k)^2$ vector. I doubt that your MATLAB code is equivalent to this on. – Henrik Schumacher Mar 3 '20 at 21:08
• 1000 random points -> ~ 500,000 distinct pairs iterated. Took a few seconds. – Gregory Mar 3 '20 at 21:10
• Ah, I see. That makes sense. The major point here is AppendTo, because it involves a copy operation of the full array. So your loop is actually of complexity $O(n^3)$. – Henrik Schumacher Mar 3 '20 at 21:10
(Subtract @@@ Subsets[Most@x, {2}])^2 + (Subtract @@@Subsets[Most@y, {2}])^2


will produce precisely the output of your OP with appropriate performance. Since your output is a small subset of all-points distances, it should also outperform things like DistanceMatrix.

A comparison of OP and this up to 300 length for timing:

• Can you please explain some of the syntax in greater detail if you have time. – Gregory Mar 3 '20 at 22:52