A constrained optimization problem

Problem statement

Consider a function which is the sum of two Gaussian bumps. For given values of the mean and standard deviation for each bump, that sum is a function of two variables.

F (x, y) := pdf_normal (x, m1, s1) + pdf_normal (y, m2, s2) $

Let's take a look at the contours of that function. We'll choose the mean and standard deviation of the Gaussian bumps. I pulled these values out of a hat.

[m1, s1] : [-2, 1.8] $

[m2, s2] : [1.3, 2.4] $

load (distrib) $

ratprint : false $

contour_plot (F, [x, -5, 5], [y, -5, 5]);

The problem is this: staying on one of the contours of F, try to get as close as possible to a given point (x1, y1). This problem is stated in this post on Stack Overflow: Calculate inverse of normal function that minimizes another function

Motivation of solution

At the point on the contour closest to (x1, y1), a circle centered on (x1, y1) is tangent to the contour, and therefore just touches the contour. As one travels around that circle, the value of F rises and falls, and has a local extremum at the tangent point. When the circle is too small, it doesn't touch the chosen contour. When the circle is too big, it crosses the contour (not tangent). So to find the right circle, and therefore the tangent location, we'll vary the size until there is a local extremum equal to the chosen contour. This is expressed by maximizing or minimizing the value of F on the circle (function fmax_circular) within solving the equation max/min of F = contour (function find_root). fmax_circular makes use of fargmax1, which is an implementation of golden section search (purpose-built, not built-in). find_root is built-in.

Solution

Choose the point to which we want to get close.

[x1, y1] : [1.82, -0.24] $

Choose the contour level on which we want to stay while we try to get close to (x1, y1).

z0 : 0.2 $

Load a script to define the functions needed to solve the problem. (You can find the script here.)

load ("SO_constrained_optimization.mac") $

Find the distance, r0, from (x1, y1) to the nearest point on the contour F(x, y) = z0.

fpprintprec : 4 $

r0 : find_root (lambda ([r], G (r, x1, y1) - z), r, 0.001, 5), z = z0, numer;

⇒

$0.8446$

Find the tangent point on the circle of radius r0 centered on (x1, y1). fargmax_circular returns a list; assume it's just one element. That's not guaranteed to work, since [t0] : ... fails when right-hand side has 2 or more elements. I expect that happens rarely, if ever.

[t0] : fargmax_circular (lambda ([t], F (x1 + r0 * cos (t), y1 + r0 * sin (t))));

⇒

$\left[2.605\right]$

(x0, y0) is the point on z0 contour of F, nearest to (x1, y1). r0 is distance from (x1, y1) to (x0, y0).

[x0, y0] : [x1 + r0 * cos (t0), y1 + r0 * sin (t0)];

⇒

$[1.094,0.1919]$

F(x0, y0) should be equal to z0.

F (x0, y0), numer;

z0;

⇒

$0.2$

⇒

$0.2$

If we draw the contour F(x, y) = z0, we should find it just touches the circle of radius r0 centered on (x1, y1).

set_plot_option ([same_xy, true]) $

set_plot_option ([legend, false]) $

load (implicit_plot) $

eq1 : F (x, y) = z0 $

eq2 : (x - x1)^2 + (y - y1)^2 = r0^2 $

implicit_plot ([eq1, eq2], [x, -5, 5], [y, -5, 5]), numer;