6
$\begingroup$

I’m looking for code that is well-suited to solving a fairly simple minimization problem:

I have a reference point $\mathbf p$ in 3D space, and I want to minimize $\|\mathbf x - \mathbf p\|^2$ subject to a constraint $F(\mathbf x)=0$.

The functions $F$ are very smooth, but typically have no useful convexity properties, and in some cases (not shown below) they might also be very complex and expensive to evaluate. The simple aspects (I think) are that there are only three independent variables, and the objective function is just a quadratic form. The problem comes from engineering geometry or CAD — I’m trying to find the minimum distance from the point $\mathbf p$ to the surface $\{\mathbf x : F(\mathbf x) = 0\}$.

Some examples of typical surfaces are:
Ellipsoid: $F(x,y,z) = 4x^2 + 9y^2 + z^2 - 1$.
Clebsch: $F(x,y,z) =64x^3 + 48x^2z - 192y^2x + 48y^2z - 31z^3 - 54z^2 - 24z$
ThreeHoles: $F(x,y,z) = x(x^2 - 3y^2) - z(z^2 - 1)$
SchwartzP: $F(x,y,z) = \cos(x) + \cos(y) + \cos(z)$
Scherk: $F(x,y,z) = \sin(z) - \tfrac12 \sinh(x) \sinh(y)$
Costa: $F(x,y,z) = (x^2 + y^2 - 1)z - (x^2 - y^2)$

And here are a few pictures: enter image description here

I know that I can attack this problem by solving the equation $(\mathbf x - \mathbf p) \times \nabla F(\mathbf x) = 0$, but I thought using a minimization algorithm might work better.

My question: There are dozens of optimization algorithms, and I’m hoping that someone can recommend one (or more) that is well-suited to my problem, perhaps by taking advantage of its simplicity.

Ideally, I’d like to get working code, as opposed to just a mathematical description. I don’t really care what programming language is used. I don’t need great accuracy (5 or 6 good digits would be OK), but speed and reliability are important.

Improve this question
$\endgroup$
16
  • 1
    $\begingroup$ Is it guaranteed that $\{x:F(x)=0\}$ is a surface, that it has co-dimension $1$? Otherwise one can encode any system $\phi_(x)=0$ via sum-of-squares $F(x)=\sum\phi_i(x)^2$. If it is a surface, then at the optimum point the line to $r$ is orthogonal to the tangent plane, $F'(x)=\lambda (x-r)$, giving generically a full system, but will all critical points of the distance function as solutions. $\endgroup$
    – Lutz Lehmann
    Commented Sep 8, 2022 at 15:05
  • 1
    $\begingroup$ The approach with Grad(F) in general would have multiple solutions, so to find a true distance to the surface one would need to find all those solutions and compare them. But geometrically solving a problem like this would be convenient putting a sphere centered at the test point, and increasing the radius of the sphere gradually until it touches the surface. This approach can be probably used as the basis of an iterative solution algorithm here. $\endgroup$
    – Maxim Umansky
    Commented Sep 9, 2022 at 3:09
  • 1
    $\begingroup$ Is the function F(x) given with an analytical expression? Do you want to solve many points for one surface? $\endgroup$
    – Amit Hochman
    Commented Sep 10, 2022 at 14:05
  • 1
    $\begingroup$ @AmitHochman. Not necessarily. But we can assume that we have a software function that returns $F(x,y,z)$ at any input point $(x,y,z)$. And yes, I typically do want to calculate the distance from many points. $\endgroup$
    – bubba
    Commented Sep 10, 2022 at 22:23
  • 1
    $\begingroup$ @hardmath. Thanks. But the problem you mentioned is much simpler than mine, and can be solved algebraically. I’m interested in numerical methods that will work on a broader range of problems. $\endgroup$
    – bubba
    Commented Sep 10, 2022 at 22:27

1 Answer 1

Reset to default
6
+50
$\begingroup$

You could try with gradient projection, here is a quick implementation in python:

import numpy as np

def project(p, f, gradf, tol=1e-2):
    u = p
    for i in range(100):
        g = gradf(u)
        v = u - f(u) * g / np.linalg.norm(g)**2
        if np.linalg.norm(u - v) < tol:
            return v
        else:
            u = v
    return v

For example, using it to project a point $p = (1, 1, 1)^T$ onto Ellipsoid from your first example:

def f(x):
    return 4 * x[0]**2 + 9 * x[1]**2 + x[2]**2 - 1

def gradf(x):
    return np.array([8 * x[0], 18 * x[1], 2 * x[2]])

p = np.array([1, 1, 1])
q = project(p, f, gradf)
q, np.linalg.norm(p - q)

we get a closest point $q$ on the surface and distance $\left\lVert p - q \right\rVert_{2}$:

(array([0.31380786, 0.0296923 , 0.77342578]), 1.2098316328577217)

Similarly, for Clebsch surface from your second example:

def f(x):
    return (64 * x[0]**3 + 48 * x[0]**2 * x[2] - 192 * x[1]**2 * x[0]
            + 48 * x[1]**2 * x[2] - 31 * x[2]**3 - 54 * x[2]**2 - 24 * x[2])

def gradf(x):
    return np.array([192 * x[0]**2 - 192 * x[1]**2 + 96 * x[0] * x[2],
                     -384 * x[0] * x[1] + 96 * x[1] * x[2],
                     -24 + 48 * x[0]**2 + 48 * x[1]**2 - 108 * x[2] - 93 * x[2]**2])

p = np.array([1, 1, 1])
q = project(p, f, gradf)
q, np.linalg.norm(p - q)

gives:

(array([1.1304619 , 0.62158043, 0.83049637]), 0.4346874290585481)

And here is a figure of surface plots with points and projections: examples

Improve this answer
$\endgroup$

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.