# Calculate the closest point to many hyperbolic paraboloids

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,1))
v  = 0.5 * np.ones((p0.shape,1))

# Init index list of elements to update

niter -= 1

# Calculate delta and update u,v
delta = inner1d(delta,delta)**0.5

# Update the index list only with items which fail the tolerance test

# 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?

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 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.