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:

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]];
 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 = ",
    If[test > threshold, NotCollinear = False;, NotCollinear = True;];
  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.

  • $\begingroup$ I wonder if collinearity can't be worked into the distance metric so that while Nearest may select collinear points when you ask for the six closest points, it will be unlikely unless all other points are far away. This would automatically sort out the compromise that you talk of as I understand it. At the same time you would not be checking the distance between an unnecessary amount of points because Nearest would use KDTree or whatever method to check distances intelligently. $\endgroup$
    – C. E.
    Commented Aug 27, 2017 at 21:51
  • $\begingroup$ That actually sounds like a very good idea! Although I see two problems: you would have to check with all pairs of points found until that point, the second one, you would have to somehow exclude the original grid point (the square grid in 2D has always 3 collinear points, but it doesn't matter - we have to exclude the central point from this). Plus, I don't know how the condition of non-collinearity is added to Nearest. I know that a criterion can be added, but this seems very complicated. $\endgroup$
    – user16320
    Commented Aug 27, 2017 at 21:55
  • 2
    $\begingroup$ (1) Det is not reliable in this setting. You want to know if the matrix has one or more "small" singular values. If you use say 10 neighbors and have at least five singular values that are not small, then you should be ble to recover the "good" points by zeroing the small singular values and reversing the SVD obtained by SingularValuesDecomposition. $\endgroup$ Commented Aug 27, 2017 at 22:09
  • $\begingroup$ Sorry, aren't all 2D points coplanar by definition? I guess your $C$ matrix is supposed to be your $M$? Is the 5th column supposed to be 1/2 the second column implying $\det(C)= 0$ for all points chosen? $\endgroup$ Commented Aug 27, 2017 at 22:09
  • 1
    $\begingroup$ I don't see how there could be no acceptable solution. EDIT: Oh crap, of course, it's meant to be $1/2 \delta y_i^2$. Just a clarification: rows of the matrix $M$ are coefficients in the Taylor series $\delta f_i = M_{11} \frac{\partial f}{\partial x} + M_{12} \frac{\partial f}{\partial y} + \cdots$, so it's clear that the last three elements should be of the second order in $\delta x$, $\delta y$. $\endgroup$
    – user16320
    Commented Aug 27, 2017 at 22:12

1 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.

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

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

(* 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.

  • $\begingroup$ Hello and thank you very much for the answer! I ran some of your examples and your algorithm unfortunately tends to pick points that are not very suitable. In fact, it seems that it always picks the farthermost out of the $numbers$ set of points. In other words, if I choose $numbers = 30$, it picks points that are 2 layers far from the grid point. If I set $numbers = 100$, those points are so they're no longer related to the grid point in any reasonable manner (not for the FEM purposes). $\endgroup$
    – user16320
    Commented Aug 28, 2017 at 10:44
  • $\begingroup$ It is fixable by fiddling with the norm in the lingandquad function - we can for example divide first two with normsq and the rest with normsq^3/2. This will favor the closest points. However, it won't work on square grid. I changed it like this: grid = Partition[Flatten@Table[{i, j}, {i, 0, 1, 0.1}, {j, 0, 1, 0.1}], 2]; now if I do it for n = 5, numbers = any, the result is five points among which 3 are collinear. How is that so? Aren't 3 collinear points supposed to give 3 linearly dependent rows? $\endgroup$
    – user16320
    Commented Aug 28, 2017 at 10:44
  • $\begingroup$ The use of different norms in different columns might be throwing things off. $\endgroup$ Commented Aug 28, 2017 at 15:02
  • $\begingroup$ I edited my post to include my clumsy code that however does exactly what I wanted. $\endgroup$
    – user16320
    Commented Aug 28, 2017 at 19:43
  • $\begingroup$ Your approach might be as good as can be done. It is tailored (so to speak) to the problem at hand. I'm just not sure about the scaling issues possibly throwing things, in terms of either making a suitable set of five neighbors seem unsuitable, or vice versa. $\endgroup$ Commented Aug 28, 2017 at 21:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.