30
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ Good question! I want it the other way round, i.e. delete points from the list for a path optimization task... $\endgroup$ Commented Mar 9, 2013 at 20:38
  • 1
    $\begingroup$ @YvesKlett you may use Nearest[data,x,n] and use the second element if the first one has been "deleted" $\endgroup$ Commented Mar 9, 2013 at 21:46
  • 1
    $\begingroup$ @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). $\endgroup$ Commented Mar 10, 2013 at 7:09

3 Answers 3

Reset to default
19
$\begingroup$

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]) /@ #][[1]], Timing[Nearest[points] /@ #][[1]]} &@ RandomReal[1, {10, 3}], {1030}]

enter image description here

Improve this answer
$\endgroup$
3
  • 7
    $\begingroup$ 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... $\endgroup$ Commented Mar 9, 2013 at 20:26
  • 2
    $\begingroup$ @DanielLichtblau In any case, I claim for myself being the bigger half $\endgroup$ Commented Mar 9, 2013 at 20:49
  • 1
    $\begingroup$ I thought you told me it usually wasn't worth admitting stupidity :P $\endgroup$ Commented Mar 9, 2013 at 21:21
8
$\begingroup$

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[[5]], 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] newnf = addpoint[nf, {5, -2}]; 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):

DensityPlot

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.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ 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}] *). $\endgroup$ Commented Mar 9, 2013 at 23:22
  • $\begingroup$ @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 $\endgroup$ Commented Mar 10, 2013 at 2:06
  • 2
    $\begingroup$ 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). $\endgroup$ Commented Mar 10, 2013 at 4:55
  • $\begingroup$ @Oleksandr That is correct, according to Daniel. $\endgroup$ Commented Mar 10, 2013 at 5:55
  • $\begingroup$ @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. $\endgroup$ Commented Mar 10, 2013 at 13:30
8
$\begingroup$

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][[1]], {}], {j, 1, nfuncs[[1]]}]; candidates = Join[Take[lonepts, {2, 1 + lonepts[[1]]}], bucketmins] /. {} :> Sequence[]; vals = Map[(# - pt).(# - pt) &, candidates]; candidates[[Ordering[vals]]][[1]] ] SetAttributes[addPoint, HoldRest]; addPoint[pt_, lonepts_, nfuncs_, ptlists_] := Catch[Module[ {allpts, nf, max1 = Length[lonepts] - 1, max2 = Length[ptlists]}, If[lonepts[[1]] < max1, lonepts[[1]] = lonepts[[1]] + 1; lonepts[[lonepts[[1]] + 1]] = pt; Throw[{}]]; lonepts[[1]] = 0; j = 1; While[j <= max2 && nfuncs[[j + 1]] =!= 0, j++]; If[j > max2, Throw[$Failed]]; If[j > nfuncs[[1]], nfuncs[[1]] = 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[ Do[addPoint[points[[j]], lonePoints, nearestFunctions, ptLists], {j, len}]] (* Out[450]= {6.370000, Null} *)

Finding nearest neighbors is somewhat slower but still not unreasonable.

Timing[ nbrs = Map[{#, nearest[#, lonePoints, nearestFunctions]} &, newpts];] (* Out[453]= {3.660000, Null} *) Max[Apply[EuclideanDistance, nbrs, {1}]] (* Out[454]= 0.0245587226349 *)

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

Timing[nf = Nearest[points];]

(* Out[455]= {0.180000, Null} *)

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

(* Out[458]= {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.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.