Nearest non-collinear/non-coplanar points

Question

This is a long post. Most of it is just an introduction so the question itself is understandable. The question is in bold. Just skip to it if you want to hear the point. Everything below the question is motivation.

I'm interested in the following: I have a grid of points, let's say for the simplicity:

Needs["NDSolve`FEM`"]
dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/4)^2, {{x, 0, 1}, {y, 0, 1}}];
grid = ToElementMesh[dom, "MeshOrder" -> 1, MaxCellMeasure -> 0.0005]["Coordinates"]

Now grid stores positions of points in the form {{x1,y1},...}.

Usually, I need to find exactly 5 closest neighbours to some point in 2D or 9 neighbours in 3D (you can easily see the pattern: together with the point itself, it's triangle number $T_n$ points in nD, but I'm particularly interested in 2D and 3D, as those are obviously somehow most relevant).

If you run the function

Nearest[grid, grid[[i]], 6]

Through every grid point 1 <= i <= Length@grid, you will find out, that some of the points are degenerate. Now I will explain what I mean by "degenerate".

Of course, this function selects the i-th grid point itself, so we throw it out:

Drop[Nearest[grid, grid[[i]], 6], 1]

(drops the first element, which is the gridpoint grid[[i]] itself)

In physics, it's often very very important, that the set of those 5 points do not contain any 3 that lie on the same line (collinear), or, equivalently, in 3D, no 5 points lie on the same plane (coplanar). If I think hard about this, I should be able to work out some condition for higher dimensions, but I'm not interested in more than 3D.

I can work out two criteria about collinearity. The first one is very intuitive. Let $p_1, p_2, p_3$ be the points of interest. Let

$$v_1 = p_1 - p_2$$ $$v_2 = p_1 - p_3$$

Now the criterion is:

$$\left| \frac{v_1 \cdot v_2}{||v_1|| ||v_2||} \right| \; \; \; \text{"is close to one"}$$

of course, we should provide some numerical threshold to define "being close to one" (equivalent to angle being somehow close to 0 or $\pi$), but that's usually the case with numerics - strictly speaking, float-type points relatively rarely end up being precisely on the same line together.

The second criterion, much more useful to me, is as follows: let $p_1, ..., p_5$ be those five points besides the original $p_0$.

Let's define:

$$\begin{aligned} \delta x_i &\equiv p_{i1} - p_{01} \\ \delta y_i &\equiv p_{i2} - p_{02} \end{aligned} $$

The following matrix is vey important:

$$ M = \begin{pmatrix} \delta x_1 & \delta y_1 & \frac{1}{2} \delta x_1^2 & \delta x_1 \delta y_1 & \frac{1}{2} \delta y_1^2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \delta x_5 & \delta y_5 & \frac{1}{2} \delta x_5^2 & \delta x_5 \delta y_5 & \frac{1}{2} \delta y_5^2 \end{pmatrix}$$

The set of 5 points $p_1, ..., p_5$ is chosen badly if

$$\det M = \text{"very small"}$$

Again, "very small" is necessary as we're dealing with numerics, not precise numbers.

My idea is:

Run Nearest[] for the grid point. Check if the matrix $M$ is not badly conditioned. If yes, add another point (so basically, run Nearest[..., 7]). Now this is kinda tricky: I would have to select all possible combinations of five points (keep in mind that the original grid point is fixed and does not contribute) and compare all matrices $M$ for them and see which one has "normal" determinant - not too small. If this won't work out, then we add one more point (Nearest[..., 8]) and again, select all 5-tuples with the original grid point fixed, compare matrices, etc... This process would have to stop eventually; if you can't find good set of 5 points and the next you're adding is on the other side of the domain, there is some problem with the domain and user should be noticed about that. Also, there would need to be some kind of compromise between how nice the matrix $M$ should be and how far are those better neighbours: every matrix' determinant is bigger and bigger as we choose farther and farther points, but this is actually not suitable for the applications. The idea is only to slightly fiddle with the nearest neighbourhood, not to venture across the domain for possible "neighbours".

Now this seems like a very complicated process that could take ages to pull off on a grid (grids can be fairly massive - from 1000 to millions of points). How can this be implemented effectively? Isn't this already implemented? Hasn't someone already done that? Thank you for any ideas thrown my way.

You can ignore the next section, it's just for an incentive.

The matrix $M$ is the basis for determining the first and second-order partial derivatives. Namely, it's inverse gives straightforward gradient and Hessian matrix. Let the points be $p_0$ for the fixed grid point and $p_1, ..., p_5$ its "correctly chosen" neighbours. Moreover, let $f_0$ and $f_1, ..., f_5$ be the functional values of a function $f$ at those points. Then:

$$ \begin{aligned} \left. \frac{\partial f}{\partial x} \right|_{p_0} &\approx M^{-1}_{1j} \delta f_j \\ \left. \frac{\partial f}{\partial y} \right|_{p_0} &\approx M^{-1}_{2j} \delta f_j \end{aligned} \; \; \; \; \; \begin{aligned} \left. \frac{\partial^2 f}{\partial x^2} \right|_{p_0} &\approx M^{-1}_{3j} \delta f_j \\ \left. \frac{\partial^2 f}{\partial x \partial y} \right|_{p_0} &\approx M^{-1}_{4j} \delta f_j \\ \left. \frac{\partial^2 f}{\partial y^2} \right|_{p_0} &\approx M^{-1}_{5j} \delta f_j \end{aligned}$$ where $\delta f_i = f_i - f_0$

Now it's obvious why I'm interested in $M$ that has a well-defined inversion.

Just for the sake of completeness: in 3D this would look similar. The matrix $M$ would have 9 rows, and 9 columns, with additional $\delta z$, $\delta z^2$, $\delta x \delta z$ and $\delta y \delta z$ entries. The criterion is the same: degenerate matrix = poorly chosen 9 neighbours.

As the 2D and 3D cases share some similarities, I think the solution in 2D is applicable in 3D as well.

1D case is very easy, the only bad thing that could happen is that two points would lie very close to each other (yes, derivative is small/small, but in numerics, let that not be too small).

The last note: problems usually emerge at the boundaries, as we don't have enough points at the other side, it's hard to approximate second-order derivative. The solution could be to add one more point from the domain, perpendicular to the boundary line, so the second derivative in the perpendicular direction could be well approximated. The question then is, which point should we throw out from the rest of them?

** EDIT 28.8.2017 **

This is my try on the problem (after I hadn't gotten any better idea and response from the community that would solve it in a way I require it to be solved):

n = 50;
nf = {};
nearest = Drop[Nearest[grid, grid[[n]], 30], 1];
i = 1;
threshold = 0.9;
IsCollinear = False;
While[(Length@nf < 5) && (i <= Length@nearest),
 adept = nearest[[i]];
 (*Print[adept];
 Pause[0.5];*)
 If[Length@nf >= 2,
  NotCollinear = True;
  j = 1;
  While[(j <= Length@nf - 1) && NotCollinear,
   k = j + 1;
   While[(k <= Length@nf) && NotCollinear,
    v1 = adept - nf[[j]];
    v2 = adept - nf[[k]];
    test = (v1.v2)^2/(v1.v1 v2.v2);
    (*Print["j= ",j,", k = ",k,", v1 = ",v1,", v2 = ",v2,", test = ",
    test];*)
    If[test > threshold, NotCollinear = False;, NotCollinear = True;];
    k++;
    ];
   j++;
   ];
  If[NotCollinear, nf = Append[nf, adept]; i++, i++];
  , nf = Append[nf, adept]; i++;]
 ]
deltas = Table[nf[[i]] - grid[[n]], {i, 1, Length@nf}];

It does exactly what I want it to: for example, points 1-11 are on the boundary, so it does the proper job of selecting neighbours that can properly approximate derivatives.

(for $n = 5$):

    SetUp[x_, y_] := {x, y, 1/2 x^2, x y, 1/2 y^2};
    Inverse[SetUp @@@ deltas]
    Inverse[SetUp @@@ Table[nearest[[i]] - grid[[n]], {i, 1, Length@nf}];]
    {{0., 0., 20., 0., -5.}, {-5., 5., 0., 0., 0.}, {0., 0., -200., 0., 
  100.}, {100., 1.42109*10^-14, 100., -100., 0.}, {100., 100., 0., 0., 0.}}

Inverse::sing: Matrix {{0.,-0.1,0.,0.,0.005},{0.,0.1,0.,0.,0.005},{0.1,0.,<<21>>,0.,0.},{0.1,-0.1,0.005,-0.01,0.005},{0.1,0.1,0.005,0.01,0.005}} is singular.

So the first list nf is "properly chosen" points, with an inverse matrix, the second one (just neighbours) is not available due to its matrix being singular.

The code above is terrible and it looks like what I'd be writing in c++. I think this is a perfect chance to learn how to code in Mathematica elegantly and how to make code brief (and so it doesn't look like a standard language code). Would someone care to rewrite this little code in Mathematica-friendly way? Thanks.

Answer 1

First I will observe that the columns in the matrix are in a sense of two different scales, and the code below normalizes by dividing linear terms by distance to center point, and quadratic points by square of said distance. I'm not positive this is the right thing to do for purposes of best assessing whether five neighbors are usable for their intended puspose, but I believe it is.

I will advocate a method that is not necessarily optimal but should at least not be pessimal (and in general seems to behave). The idea is to take 30 or so neighbors, form a matrix of the deltas normalized as above, compute its LU decomposition, and use the rows corresponding to the first five elements of the permutation vector used in the decomposition (in effect it records row repositioning).

Here is the example from the post.

Needs["NDSolve`FEM`"]
dom = ImplicitRegion[(x - 1/2)^2 + (y - 1/2)^2 >= (1/4)^2, {{x, 0, 
     1}, {y, 0, 1}}];
grid = ToElementMesh[dom, "MeshOrder" -> 1, MaxCellMeasure -> 0.0005][
   "Coordinates"];

These next two utility routines give differences between a given point and a set of points, and form the five component vectors from a pair of two such differences.

getDeltas[{pt_, ptlist_}] := Map[# - pt &, ptlist]
linAndQuad[{x_, y_}] := 
 With[{normsq = x^2 + y^2}, {x/Sqrt[normsq], y/Sqrt[normsq], 
   x^2/normsq, x*y/normsq, y^2/normsq}]

For each point in our set, find its 30 nearest neighbors. Then form the 20x5 matrix of (normalized) first and second order differences.

nf = Nearest[grid];
numnbrs = 30;
nbrs = Map[{First[#], Rest[#]} &, nf[grid, numnbrs + 1]];
deltas = Map[getDeltas, nbrs];
delMats = Map[linAndQuad, deltas, {2}];

Now we look at a specific such matrix, the one arising from neighbors of the 44th point in the grid.

n = 44;
{lumat, perm, cond} = 
  LUDecomposition[Join[delMats[[n]], IdentityMatrix[numnbrs], 2]];

Find the five "best" points, as assessed by the partial pivoting strategy of LUDecomposition. They may not be "best" for this particular purpose but typically they will at least be reasonable, and better than blindly using the first five. We'll look at the singular values for the differences matrix formd by these five.

dmat = delMats[[n, perm[[1 ;; 5]]]];
SingularValueList[dmat]

(* Out[237]= {2.01935, 1.77039, 1.29386, 0.709237, 0.227315} *)

Compare to the matrix formed from the five nearest neighbors' differences.

SingularValueList[delMats[[n, 1 ;; 5]]]

(* Out[238]= {2.16909, 1.58235, 1.26157, 0.714453, 0.0807256} *)

This seems to be fairly typical. When the five nearest neighbors give a modestly well-behaved set of differences, smallest singular value for the "better" points is around 1.5-3 times bigger (in this example, and using 30 neighbors). When the five closest give a bad matrix then the quality difference is more marked.