I have some black-box function that takes lists of 2D points, and I would like to find some practical methods of searching for collections of points that maximize it.


I've come across the need to search for sets of points in the plane that satisfy/optimize constraints so many times, but have never adequately found workable strategies in Mathematica. To name a few of these times:

  • When packing geometric entities in various ways
  • Searching for graph embeddings with certain properties
  • Various combinatorial optimization problems

Minimal Example

Kobon triangles (mentioned in comments), is a nice example of this type of problem: searching for arrangements of lines that make the most triangles. Jason B. provided a brute-force solution to this problem by converting the lines to a graph, which is great, but to illustrate my question let's look at it this way: consider the following function f that takes (flattened) lines as arguments

(* the ith line is `InfiniteLine[{{xi1, yi1}, {xi1, yi2}}] *)
f[x11, y11, x12, y12, x21, y21, x22, y22, ...] := ...

and returns the number of non-overlapping triangles. Here are some inputs pts where f is 1:

pts = RandomReal[1, {3, 2, 2}];
lines = RegionIntersection[#, Disk[{0, 0}, 10]] & /@ (InfiniteLine /@ pts);
tri = triangles[lines];
Graphics[{lines, LightBlue, Triangle /@ tri}, Frame -> False, PlotRange -> Full]

enter image description here

Where I get stuck

I first tried using FindMaximum on the simplest case of only 3 lines (i.e. find 1 triangle from 12 coordinates). As you can see below, it fails to find a set of 6 points in the plane, even when supplied with the correct coordinates as starting values (which was very strange to me):

d = Disk[{0, 0}, 20];
s = FindMaximum[{f[x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, 
    y6], {x1, y1} \[Element] d && {x2, y2} \[Element] 
     d && {x3, y3} \[Element] d && {x4, y4} \[Element] 
     d && {x5, y5} \[Element] d && {x6, y6} \[Element] d}, 
    {{x1, 0.7242101930120808`}, {y1, 0.029677803428552307`}, {x2, 0.00816459598186059`}, {y2, 0.47699819634429996`}, {x3, 0.06540238199229464`}, {y3, 0.69214964488526`}, {x4,0.6253484480271967`}, {y4, 0.9545360094125843`}, {x5, 0.12098288193445206`}, {y5, 0.7013188710966709`}, {x6, 0.2806641729642241`}, {y6, 0.926496649677153`}}]

enter image description here

Neither NMaximize nor FindInstance seem up to the task either:

NMaximize[{f[x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, 
    y6], {x1, y1} \[Element] d && {x2, y2} \[Element] 
     d && {x3, y3} \[Element] d && {x4, y4} \[Element] 
     d && {x5, y5} \[Element] d && {x6, y6} \[Element] d}, {x1, y1, 
   x2, y2, x3, y3, x4, y4, x5, y5, x6, y6}, 
  Method -> "DifferentialEvolution"]

enter image description here

FindInstance[f[x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, x6, y6] == 1, {x1, y1, x2,
y2, x3, y3, x4, y4, x5, y5, x6, y6}, Reals, 1]

enter image description here

When the variables are continuous and the function you are looking to maximize is integer-valued, this is no doubt a hard problem. However, are there are definitely better ways out there to solve this than grid/random search?

Code for f in the above example

f[x1_, y1_, x2_, y2_, x3_, y3_, x4_, y4_, x5_, y5_, x6_, y6_] := 
 Module[{l, h, lines},
  l = InfiniteLine /@ {{{x1, y1}, {x2, y2}}, {{x3, y3}, {x4, 
       y4}}, {{x5, y5}, {x6, y6}} };
  lines = RegionIntersection[#, Disk[{0, 0}, 10]] & /@ l;

triangles[lines:{__}]:= Module[
    {lineSegments, vertices, edges, triangles},

    lineSegments = Flatten[
        Map[Function @ Partition[Sort @ #, 2, 1],
                List @@@ Map[RegionIntersection[Part[lines, n], #]&, Delete[lines, n]]},
                Flatten[Select[l, Depth[#] >= 3&], 1]
                {n, Length @ lines}

   vertices = Flatten[lineSegments, 1] // DeleteDuplicates;
   edges = lineSegments /. MapIndexed[#1 -> First@#2 &, vertices];
   triangles = FindCycle[Graph[#1 \[UndirectedEdge] #2 & @@@ edges], {3}, All];
   triangles = triangles[[All, All, 1]];
   vertices[[#,1]] & /@ triangles

Additional Details

  • A valid assumption is that all the points are in some Disk[{0, 0}, n] or another fixed region.

  • Since this is useful for finding geometric configurations subject to constraints - perhaps FindGeometricConjectures would be a better fit to try (when v12 comes out)?

  • 3
    $\begingroup$ This seems familiar. Are you trying to use FindMaximum to find the maximum Kobon triangles? $\endgroup$
    – Jason B.
    Sep 6, 2018 at 20:19
  • 1
    $\begingroup$ No just trying to see how to get Mathematica to move points around in 2d to optimize any function $\endgroup$
    – M.R.
    Sep 19, 2018 at 16:27
  • 2
    $\begingroup$ It is unlikely you are going to have much luck using these methods when the function you are tackling has only integer values. $\endgroup$
    – kirma
    Sep 19, 2018 at 16:51
  • 1
    $\begingroup$ In my case, there is no way to solve for the vertices and anyway I'd like a black box technique for similar problems. At the very least, optimizing an integer valued function over configurations in a region could be done better than random search with BayesianMaximize no? $\endgroup$
    – M.R.
    Sep 25, 2018 at 22:28
  • 1
    $\begingroup$ So the lists can be of arbitrary length, correct? You can't just parameterize a list of 2D points length $L$ with $2L$ parameters? $\endgroup$
    – MikeY
    Jan 19, 2019 at 17:04


Your Answer

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