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
numpy.core.umath_tests instead of
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 mask = np.arange(p0.shape) while niter and mask.shape: 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?