# Efficiently appending points to a NearestFunction

I have a program which works with a list of points (in 3D, but it could be in 1D as well, it won't make a difference). Within a loop, it keeps adding new points to the list. To calculate the coordinates of a new point, it needs to use a NearestFunction generated from the list of existing points.

When a new point is added to the list, the NearestFunction needs to be rebuilt as well. This is rather time consuming if I do it using the list of existing points.

Is it possible to more efficiently append a new point to an already constructed NearestFunction? (Perhaps by making use of the knowledge of the internal structure of a NearestFunction)

Here's an example to illustrate the question better:

points = RandomReal[1, {1000, 3}]; (* we have lots of points ... *)
nf = Nearest[points];  (* ... and the corresponding NearestFunction *)
AppendTo[points, RandomReal[1, 3]]; (* we add an extra point, ... *)
nf = Nearest[points]; (* so we need to rebuild the NearestFunction,
but doing it from scratch is inefficient *)


Rebuilding the NearestFunction for a large number of points is not very efficient. Is there a way to add a single point more efficiently, making use of the existing NearestFunction?

• Good question! I want it the other way round, i.e. delete points from the list for a path optimization task... Mar 9 '13 at 20:38
• @YvesKlett you may use Nearest[data,x,n] and use the second element if the first one has been "deleted" Mar 9 '13 at 21:46
• @belisarius true. For a long list that is whittled down to zero I' ll have to repeatedly return more candidates which takes longer (but I need to benchmark to verify, been a while). Mar 10 '13 at 7:09

I know this isn't exactly what you want, but just a stupid idea:

ClearAll[newf];
points = RandomReal[1, {1000000}];(*we have lots of points...*)
nf = Nearest[points];(*... and the corresponding NearestFunction*)
newf[oldf_, newpoints_List] := (Nearest[Union[oldf[#], Nearest[newpoints][#]], #] &);

newf[nf, {3, 4, 5}][1.98]


Edit

Here is a version that reconstructs the NearestFunction[] only after some iterations:

ClearAll[newf];
points = RandomReal[1, {10, 3}];(*we have lots of points...*)
nf = Nearest[points];(*... and the corresponding NearestFunction*)
newf[oldf_, allpoints_, newpoints_List, max_] :=
If[StringCount[ToString@oldf, "Nearest"] > max,
Nearest[allpoints],
(Nearest[Union[oldf[#], Nearest[newpoints][#]], #] &)]

newP = RandomReal[1, {10, 3}];
max = 5;
Do[
Print@(nf = newf[nf,
points = Union[points, newP],
newP = RandomReal[1, {10, 3}],
max])[{1, 1, 1}],
{30}]


Edit

Here you can see the efficiency limits for the method. We start with 2 10^5 triplets and add a point at a time. We use the above method and the standard Nearest Function recalculation to evaluate the distances to ten random points. We see that after adding 10^3 points, the suggested method efficiency gain is almost lost, but can be regained with only one recalc.

ClearAll[newf];
$HistoryLength = 0;$RecursionLimit = 2500;

newf[oldf_, allpoints_, newpoints_List, max_] :=
If[StringCount[ToString@oldf, "Nearest"] > 2 max,
Nearest[allpoints],
(Nearest[Union[oldf[#], Nearest[newpoints][#]], #] &)]

points = RandomReal[1, {200000, 3}];(*we have lots of points...*)
nf = Nearest[points];(*... and the corresponding NearestFunction*)

newP = RandomReal[1, {1, 3}];
max = 1000;(*Evaluate NearestFun each 1000 points)*)
ListLinePlot@
Transpose@
Table[{Timing[(nf = newf[nf, points = Union[points, newP],
newP = RandomReal[1, {1, 3}], max]) /@ #][],
Timing[Nearest[points] /@ #][]} &@
RandomReal[1, {10, 3}], {1030}] • This is the method that I would advocate, so that's two of us. I'm not sure if that makes it half or twice as stupid... Mar 9 '13 at 20:26
• @DanielLichtblau In any case, I claim for myself being the bigger half Mar 9 '13 at 20:49
• I thought you told me it usually wasn't worth admitting stupidity :P Mar 9 '13 at 21:21

NearestFunction seems to have a structure that can be guessed at:

pts = {{1, 0}, {0, 2}, {3, 1}, {-1, 4}};
nf = Nearest[pts];

List @@ nf

{1, (* unknown *)
{4, 2}, (* number of points, dimension *)
3, (* unknown *)
{{1., 0., 3., -1.}, {0., 2., 1., 4.}}, (* transpose of points as Real *)
{{1, 0}, {0, 2}, {3, 1}, {-1, 4}}, (* original points *)
None, (* unknown *)
Automatic, (* distance function *)
Hold[Nearest[{{1, 0}, {0, 2}, {3, 1}, {-1, 4}}]]} (* original code *)


I don't know what all these things stand for, so what follows might not work in all cases. In any case, what I present below will suggest that adding a point and redoing Nearest is probably not a bad way.

This will change what obviously needs changing:

addpoint[nf_, newpoint_] :=
MapAt[# + 1 &, (* update number of points *)
MapAt[
MapIndexed[Append[#, N@newpoint[[First[#2]]]] &, #] &, (* add coordinates to transposed list *)
ReplacePart[nf, {{5}, {-1, 1, 1}} -> Append[nf[], newpoint]], (* add new point to list of points in both places it appears *)
4],
{2, 1}]


It works on points.

pts = {{1, 0}, {0, 2}, {3, 1}, {-1, 4}};
nf = Nearest[pts]

DensityPlot[newnf[{x, y}], {x, -0, 6}, {y, -3, 5}]

mannf = Nearest[Append[pts, {5, -2}]];
DensityPlot[mannf[{x, y}], {x, -0, 6}, {y, -3, 5}]


Both give the same output (with the same artifacts): The main disappointment is that it is no faster. Timings on 1000, 10000 points of the following were about the same.

points = RandomReal[1, {1000, 3}];(*we have lots of points...*)
nf = Nearest[points];(*... and the corresponding NearestFunction*)

(AppendTo[points, {1, 2, 3}];(*we add an extra point,...*)
anf = Nearest[points];) // Timing

(bnf = addpoint[nf, {1, 2, 3}]) // Timing


Comparing the outputs:

anf == bnf
(* -> True *)


The first (recomputing Nearest) was quicker in more trials than mine. Since the first is definitely safe, I'd say stick with it. Perhaps someone may come up with a more efficient way of updating the parts of NearestFunction.

• Unfortunately, your function is not working. The last element of NearestFunction is a fallback that rebuilds the function. It is being invoked each time. Consider ReplacePart[addpoint[nf,point],-1->Hold[banana]][{0,0,0}]: (* banana[{0,0,0}] *). Mar 9 '13 at 23:22
• @Xerxes Thanks for testing it. I'm not sure I understand. I found a error in the last example (from pasting out of context), but the function addpoint seems to work for me. I get the same NearestFunction whether I use AppendTo or addpoint Mar 10 '13 at 2:06
• It's not too surprising that this isn't faster. The internal structure of a NearestFunction is most likely a k-d tree or similar. The FullForm is probably only used to (re)build the internal representation when it is missing (e.g. when the NearestFunction is loaded from a package). Mar 10 '13 at 4:55
• @Oleksandr That is correct, according to Daniel. Mar 10 '13 at 5:55
• @OleksandrR. Actually from the beginning I was doubtful of any way to do what's being asked -- belisarius' is a clever way around it. But I still thought it worth sharing a negative result, even of a naive approach. Mar 10 '13 at 13:30

Here is a variation that is a bit on the slow side but will at least scale fairly well. The idea is to use buckets of NearestFunction objects, as well as a list of "lone" points, and take closest neighbors from amongst neighbors obtained from these separately. Adding new points gets amortized in the sense that usually we add to the "lone" list, then empty that at a threshold and create a new NearestFunction. Which size to create depends on what we already have; if all lists to date are full, we throw them all together 9effectively doubling the size of the largest) and create a new function from those. We then empty the lower slots. When we next need to create a NearestFunction we will fill in those lower slots.

minLength = 32;

nearest[pt_, lonepts_, nfuncs_] := Module[
{bucketmins, candidates, vals},
bucketmins =
Table[If[nfuncs[[j + 1]] =!= 0, nfuncs[[j + 1]][pt][], {}], {j,
1, nfuncs[]}];
candidates =
Join[Take[lonepts, {2, 1 + lonepts[]}], bucketmins] /. {} :>
Sequence[];
vals = Map[(# - pt).(# - pt) &, candidates];
candidates[[Ordering[vals]]][]
]

addPoint[pt_, lonepts_, nfuncs_, ptlists_] := Catch[Module[
{allpts, nf, max1 = Length[lonepts] - 1, max2 = Length[ptlists]},
If[lonepts[] < max1,
lonepts[] = lonepts[] + 1;
lonepts[[lonepts[] + 1]] = pt;
Throw[{}]];
lonepts[] = 0;
j = 1;
While[j <= max2 && nfuncs[[j + 1]] =!= 0, j++];
If[j > max2, Throw[\$Failed]];
If[j > nfuncs[], nfuncs[] = j];
allpts =
Join[{pt}, Rest[lonepts],
Flatten[Table[ptLists[[k]], {k, j - 1}], 1]];
Do[ptlists[[k]] = {}; nfuncs[[k + 1]] = 0, {k, j - 1}];
ptlists[[j]] = allpts;
nfuncs[[j + 1]] = Nearest[allpts];
]]


Example from @belisarius:

lonePoints = ConstantArray[0, minLength + 1];
nearestFunctions = ConstantArray[0, minLength + 1];
ptLists = ConstantArray[{}, minLength];

Timing[
len}]]

(* Out= {6.370000, Null} *)


Finding nearest neighbors is somewhat slower but still not unreasonable.

Timing[
nbrs = Map[{#, nearest[#, lonePoints, nearestFunctions]} &, newpts];]

(* Out= {3.660000, Null} *)

Max[Apply[EuclideanDistance, nbrs, {1}]]

(* Out= 0.0245587226349 *)


This is no match for using Nearest directly, if we don't need to ever add points.

Timing[nf = Nearest[points];]

(* Out= {0.180000, Null} *)

Timing[nbrs2 = Map[{#, nf[#][]} &, newpts];]

(* Out= {0.130000, Null} *)

But it will behave well if we have to find neighbors and then add the new points to the set. Possibly this all could be sped up with judicious use of Compile.