Minimum Distance to f(x) with Scipy Optimize

 

Shortest Distance Between f(x) and a Point¶

In this notebook, I show a simple demonstration of using Scipy Optimize to find shortest distance from a given point $(𝑥_p,𝑦_p)$ to a curve: 𝑦 = 𝑓(𝑥).

I will need to import a few Python libraries for this.

In [1]:

import numpy as np
import scipy.optimize as optimize
import matplotlib.pyplot as plt

Scipy Optimize has an impressive array of functions to support common optimization problems. For this problem, I will use scipy.minimize to find the minimum distance. minimize requires an objective function and an initial guess for $x_{opt}$. f(x) is the "curve" we want to find the shortest distance to. This can be any continuous function of x.

In [2]:

def f(x):
    return np.cos(x) # + np.sin(x)


def dist(x):
    """
    Returns the Euclidean distance from (xp, yp) to (x, f(x))
    """  
    y = f(x)
    dist = ((x-xp)**2 + (y-yp)**2)**.5
    
    return dist

I define a reference point (xp, yp) from which our distance is calculated and an initial guess to get the optimizer going. Note that this problem has no constraints or bounds.

In [3]:

xp, yp = .5, .5  # reference point
x0 = 1  # initial guess

Now I call optimize.minimize and the magic happens. The default algorithm for optmize is BFGS (Broyden-Fletcher-Goldfarb-Shanno) when no constraints or bounds are applied. minimize can solve with many other methods but BFGS works good for this simple problem.

In [7]:

result = optimize.minimize(dist, x0)                        
if result.success:
    xopt = result.x[0]
    print('Solution: {:.4f}'.format(xopt))
    print('dist(x) @ optimum: {:.3f}'.format(dist(xopt)))
    print('Iterations: {}'.format(result.nit))
    
else:
    raise ValueError(result.message)

Solution: 0.6754
dist(x) @ optimum: 0.331
Iterations: 4

In 4 iterations the BFGS method has found the optimal x = 0.675 with a minimum distance of 0.331. Here is a graphical view of the solution.

In [8]:

x = np.arange(xp - 1, xp + 1, .1)
plt.plot(x, f(x), '-g', label='f(x)' )
plt.plot(xp, yp, '+r', label='ref point')
plt.plot(xopt, f(xopt), 'ob', label='optimum')
plt.legend()
plt.axis('scaled')
plt.tight_layout()

ax = plt.gca()
circle1 = plt.Circle((xp, yp), dist(xopt), color='k', ls=':', fill=False, )
ax.add_patch(circle1);

I added a circle to highlight the validity of the solution. If f(x) is differentiable you may be able to solve this problem with calculus. In this case, the derivative of the distance squared is $2x-2x_p-2sin(x)cos(x)+2sin(x)y_p$ Substitute $x_p$ and $y_p$, set equal to zero and we can find x = 0.6754
Scipy optimize can do a whole lot more. Check out some examples here.

In [ ]: