In this question I asked for a way to compute the closest projected point to a hyperbolic paraboloid using python.

Using the iterative approximation answer, I'm able to use the code below to calculate the closest point to multiple paraboloids.

import numpy as np

# This function calculate the closest projection on a hyperbolic paraboloid
# As Answered by @Jaime https://stackoverflow.com/questions/18858448/speeding-up-a-closest-point-on-a-hyperbolic-paraboloid-algorithm

def solve_uv(p0, p1, p2, p3, p, tol=1e-6, niter=100):
    a = p1 - p0
    b = p3 - p0
    c = p2 - p3 - p1 + p0
    p_ = p - p0
    u = 0.5
    v = 0.5
    error = None
    while niter and (error is None or error > tol):
        niter -= 1

        u_ = np.dot(p_ - v*b, a + v*c) / np.dot(a + v*c, a + v*c)
        v_ = np.dot(p_ - u*a, b + u*c) / np.dot(b + u*c, b + u*c)
        error = np.linalg.norm([u - u_, v - v_])
        u, v = u_, v_
    return np.array([u, v])

# Generate random hyperbolic paraboloids
COUNT = 1000
p0 = np.random.random_sample((COUNT,3))
p1 = np.random.random_sample((COUNT,3))
p2 = np.random.random_sample((COUNT,3))
p3 = np.random.random_sample((COUNT,3))
p = np.random.random_sample(3)

iterated = np.empty((COUNT,2))
for i in xrange(COUNT):
    iterated[i,:] = solve_uv(p0[i], p1[i], p2[i], p3[i], p)

This works fine, but obviously doesn't scale well and becomes very slow when dealing with large data sets.

To remedy this I made some changes to vectorize the function. It will stop iterating until all the given items pass the tolerance test or reach the maximum iteration count. I also use inner1d from numpy.core.umath_tests instead of np.dot and np.linalg.norm because in my tests inner1d performs better with large data sets.

import numpy as np
from numpy.core.umath_tests import inner1d

def solver_vectorized(p0, p1, p2, p3, p, tol=1e-6, niter=100):
    a = p1 - p0
    b = p3 - p0
    c = p2 - p3 - p1 + p0
    p_ = p - p0
    u  = 0.5 * np.ones((p0.shape[0],1))
    v  = 0.5 * np.ones((p0.shape[0],1))

    # Init index list of elements to update
    mask = np.arange(p0.shape[0])

    while niter and mask.shape[0]:
        niter -= 1

        avc = a[mask] + v[mask]*c[mask]
        buc = b[mask] + u[mask]*c[mask]
        u_ = (inner1d(p_[mask] - v[mask]*b[mask], avc) / inner1d(avc,avc))[:,np.newaxis]
        v_ = (inner1d(p_[mask] - u[mask]*a[mask], buc) / inner1d(buc,buc))[:,np.newaxis]

        # Calculate delta and update u,v
        delta = np.hstack([u[mask]-u_, v[mask]-v_])
        delta = inner1d(delta,delta)**0.5
        u[mask], v[mask] = u_, v_

        # Update the index list only with items which fail the tolerance test
        mask = mask[np.where(delta > tol)]

    # Return uv's
    return np.hstack([u, v])

vectorized = solver_vectorized(p0, p1, p2, p3, p)
print np.allclose(iterated,vectorized) #Returns: True


On very large data sets (over 100,000 items) the vectorized version will solve in ~1.5 seconds. But I am sure there are some inefficiencies that could be improved. For example the square root could probably be dropped for a small performance increase. Or by using masks differently. How would you change the function to improve its performance?


1 Answer 1


How would you change the function to improve its performance?

I agree that using masks differently could prove fruitful. Here is my suggestion. Currently the tolerance is a constant. Suppose we substitute it with cur_tol which starts out as e.g. 100 * tol, and change the loop exit criterion to (niter and mask.shape[0] AND cur_tol <= tol)? Then perform cur_tol /= 2 in each iteration. The idea is to have mask reference fewer elements during early updates, and for each iteration to simultaneously get closer to the answer AND include a potentially larger number of input values that will drive us closer to that answer. So we are lazy about evaluating a smaller number of points when we clearly are far from the answer, and cast a wider net as we approach closer to it.


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.