Largest Empty Rotated Square: lots of trigonometry

Work in progress...

Here, I delve into the problem of finding the largest empty rotated square. Given \(n\) data points, find the largest square rotated by an angle \(\theta\), such that the square does not contain any of the points.  In [1], I discussed two special cases: 

• an axis-aligned square (angle is 0°),
• a diamond shape, which is a square rotated by 45°.

Here we focus on arbitrary angles.

1. Exogenous Angle

Here, we assume we know the rotation angle \({\color{darkblue}\theta}\) in advance.

We need to check whether a point is outside the rotated square. That is a bit difficult. So, instead, we rotate all data points by \(-{\color{darkblue}\theta}\).

We use the same data set as in [1]:

Instead of rotating our square, we rotate the data. That will result in something like:

Here, the angle is 25°.

Now, we can simply compare the new points against an axis-aligned square. This rotation by \(-{\color{darkblue}\theta}\) is implemented as follows:  \[\begin{align} & {\color{darkblue}p}'_{i,x} := {\color{darkblue}p}_{i,x} \cos({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkblue}\theta}) \\ & {\color{darkblue}p}'_{i,y} := -{\color{darkblue}p}_{i,x} \sin({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkblue}\theta})  \end{align}\] These formulas can be easily derived from the standard rotation matrix [2] and using \(\sin(-x)=-\sin(x)\) and \(\cos(-x)=\cos(x)\). Now we can just try to find a large standard square, as documented in [1]. It is noted that this is just a data-manipulation step; it is not a model equation.

One complication is that we need to ensure the square is within our outer box. This is now much more complicated. We need to construct all corner points (that is easy), and then rotate them back to our original coordinate system. Let's denote by \({\color{darkred}q}_{j,c}'\) the coordinates of the four corner points in the transformed space (so \(j=1,\dots,4\)). Then rotate them back to \({\color{darkred}q}_{j,c}\) and make sure that all corner points are inside \([0,{\color{darkblue}{\mathit{maxSize}}}] \times [0,{\color{darkblue}{\mathit{maxSize}}}]\). A complete model can look like:

Exogenous Angle Model

\[\begin{align}     \max\> & {\color{darkred}s} \\&{\color{darkblue}p}'_{i,x} \le {\color{darkred}x}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\&{\color{darkblue}p}'_{i,x} \ge {\color{darkred}x}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,2} && \forall i \\&{\color{darkblue}p}'_{i,y} \le {\color{darkred}y}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,3} && \forall i \\&{\color{darkblue}p}'_{i,y} \ge {\color{darkred}y}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i \\&{\color{darkred}q}_{j,x}' = {\color{darkred}x}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,x} && \forall j \\&{\color{darkred}q}_{j,y}' = {\color{darkred}y}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,y} && \forall j  \\&{\color{darkred}q}_{j,x} ={\color{darkred}q}_{j,x}' \cos({\color{darkblue}\theta}) - {\color{darkred}q}_{j,y}' \sin({\color{darkblue}\theta}) && \forall j \\&{\color{darkred}q}_{j,y} ={\color{darkred}q}_{j,x}' \sin({\color{darkblue}\theta}) + {\color{darkred}q}_{j,y}' \cos({\color{darkblue}\theta}) && \forall j \\&{\color{darkred}q}_{j,c} \in [0, {\color{darkblue}{\mathit{maxSize}}}] \\&{\color{darkred}s} \in [0,{\color{darkblue}{\mathit{maxSize}}}] \\&{\color{darkred}\delta}_{i,k} \in \{0,1\} \\&{\color{darkblue}p}'_{i,x} := {\color{darkblue}p}_{i,x} \cos({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkblue}\theta}) \\&{\color{darkblue}p}'_{i,y} := -{\color{darkblue}p}_{i,x} \sin({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkblue}\theta}) \\&{\color{darkblue}\Delta} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \\  1 & 0 \\ 1 & 1\end{pmatrix} \end{align}\]

  

Notes:

• \(\sin({\color{darkblue}\theta})\) and \(\cos({\color{darkblue}\theta})\) are constants. Unbelievably, this is actually a linear MIP model!
• \({\color{darkred}x}',{\color{darkred}y}'\) denotes the lower-left corner of the axis-aligned square (like in [1]). 
• \({\color{darkred}q}_{j,c}', j=1,\dots,4\) are the four corner points (in transformed space). We map them back to \({\color{darkred}q}_{j,c}\), which form a rotated square.
• Below, we plot the solution both in the transformed space and in the original coordinate system.
• This is a very easy MIP:
  
     Proven optimal solution
     MIP Solution:            3.085938    (657 iterations, 0 nodes)

The solution in transformed coordinates looks like:

When we transform this back to our original coordinate system, we see:

This solution is slightly better (i.e., we have a larger optimal square) than the simple axis-aligned square in [1].

Obviously, we already have two good test cases for this model from [1]: 0° and 45°. 

Next thing to try: make \(\theta\) endogenous.

References

1. Yet Another Math Programming Consultant: Largest empty shapes
2. Rotation matrix - Wikipedia

Largest empty shapes

Finding the largest empty shapes

This post is about finding empty regions (square, rectangle, or circle) in a (large) collection of given points. These are interesting little optimization problems.

Data

I generated \(n=100\) data points \((x_i,y_i), i=1,\dots,n\) drawn from a uniform distribution \[\begin{align}&x_i \sim U(0,10) \\ & y_i \sim U(0,10)\end{align} \] 

It is always a good idea to have a careful look at the data. I.e., print and plot it (for very large datasets, print or plot a subset), and collect some elementary statistics (min, max, average, etc.). In practical modeling, it is surprising how many problems arise from simple data errors.  

Convex hull models

Convex hull as an optimization problem

In the previous posts, the construction of a convex hull played a significant role: it was an easy way to reduce the size of the data sets, often by a large amount. Somehow, it escaped me that we can try to formulate this as what turns out to be a conceptually rather simple optimization problem. This has probably little or no practical value. But it remains an interesting exercise.

Minimum enclosing circle/ellipse 2

In [1] where I discussed how to find the minimum enclosing circle and minimum enclosing ellipse around a set of points. This is a follow-up post where I extend this to sets of circles and ellipses. 

1. MINIMUM ENCLOSING CIRCLE

Here our data is a set of \(n\) circles (or disks) of different size. We want to find the smallest circle that contains all these circles. 

An example data set with random circles can look like:

----     30 PARAMETER circles  coordinates of center and radius

                   x           y           r

circle1        4.294      21.082       1.663
circle2       13.759       7.528       4.014
circle3        7.305       5.601       1.961
circle4        8.746      21.407       6.235
circle5        1.678      12.505       2.591
circle6       24.953      14.468       2.715
circle7       24.778      19.056       4.564
circle8        3.267      15.993       5.336
circle9        3.988       6.252       4.769
circle10      16.723      10.884       3.783
circle11       8.993       8.786       3.480
circle12       3.287       3.753       1.706
circle13      14.728      20.772       2.885
circle14       5.770      16.643       1.279
circle15      19.396       7.591       3.031

Minimum enclosing ellipse

Minimum encompassing circle and ellipse

This is again about finding the smallest geometric shape containing all our data points. Here I focus on easy, convex cases: a circle and an ellipse. After the last post, where I was struggling mightily with a non-convex version of this model, this should be a breeze. 

I'll discuss SDP (semidefinite programming), SOCP (second order cone programming), rotated SOCP and traditional NLP approaches using gradient based solvers. So there is a lot of ground to cover.

First let's do circles.

Revisiting a crazy global NLP problem

Minimum encompassing triangle

This looks like a simple problem. Given \(n\) 2d points, find the smallest encompassing triangle. I follow the formulation from [1]. This post is self-contained, so you don't have to go back to [1]. The main contribution here is how one could work around the performance issues of the original formulation. 

Summary

The problem can be formulated as a nonconvex quadratic problem. However, this leads to models that are very difficult for global NLP solvers: global optimality can not be proven in reasonable time, even for very small data sets. One possible way to attack the problem is to use a different measure for size of the triangle: circumference instead of the more obvious area of a triangle. Although there are explicit algorithms for this problem [3,4], here we are trying to use an optimization model.

Experience with NLP solvers on a simple economic growth model

Growth Models

In this post, we formulate and solve a relatively simple economic growth model where we want to find an optimal savings rate. Consumers can either use their income for direct consumption or they can save. Savings lead to investment and to more income down the road. The resulting NLP (nonlinear programming) model is not too complicated, and it is interesting to see how we can implement this model using different tools.

I often get the question: why not use Python or R for my (economic) models? This note tries to give an answer using a rather smallish model. Although the model is small and simple, it actually illustrates some of the problems we can face when using say Python or R code. Some of the code I present here is not very obvious, or self-explanatory. One reason is that some algorithm developers look at NLP problems as something like: \[\begin{align}\min\>&f(x) \\ &g(x)=0 \\& h(x)\ge 0\end{align}\] Unfortunately, such a representation is close to the solver, but far away from practical models. It is just the wrong abstraction level. The model implementation impedance becomes even larger if we need to provide analytic gradients. Time we could have spent on developing better models is instead wasted on error-prone low-level detailed programming. This problem is even more costly, when we think about maintenance such as bug-fixing and adding new features to the model. We'll see how this looks like for this model.  

The underlying economic model is often called a Ramsey growth model, after Frank Plumpton Ramsey [1].

Crack the passcode

In this puzzle [1,2], we need to determine what the 3-digit passcode is, using a few hints. Each digit is an integer between 0 and 9. The hints are:

Let's see if we can shoehorn this into a MIP model.

Sorting: minimize number of swaps

In this post, I want to delve further into sorting. In a question on or.stackexchange.com [1], the subject was minimizing the number of swaps in sorting algorithms. A swap, i.e., an interchange of two items, is a basic operation in sorting. We usually don't pay much attention to this. First, we assume swaps are cheap. If they are not, we can sort not the real data (which can be complex, and thus more complicated and expensive to swap), but rather pointers to the data or keys. But what if we really have to reorder bulky physical things? Then the number of swaps is a more meaningful concept. Can we minimize the work to reorder say \(n\) boxes (or containers to make it more dramatic) by minimizing the number of swaps (assuming interchanging the position of two boxes is the way reordering works).

Clock problem

From [1]:

The hour, minute and second hands of this clock are all the same length and move smoothly in a circle. The dial contains hour and minute markers, but the numbers are missing. Therefore, it’s impossible to tell which one of the 12 hour markers belongs to the 12. The two hands on the left are positioned exactly on hour markers, and the hand on the right is positioned between a minute and an hour marker. What time does the clock show?

It is possible to solve this without really any math, but, of course, here I try to model this as a mathematical programming model.