Calculate the closest point to many hyperbolic paraboloids

Question

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
np.set_printoptions(suppress=True)  


# 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

Question

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?

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.