I have a (non-invertible) function ak([u,v,w]) This takes a point on the surface of the unit octahedron (p: such that |u|+|v|+|w| = 1) and returns a point on the surface of the unit sphere. The function isn't perfect but the intention is to keep the distance between points authalic. I was thinking of using SciPy minimize to provide a numerical inverse, but I cannot wrap my head around it. input: spherical pt [x,y,z], output octahedral pts [u,v,w] such that ak([u,v,w])=[x,y,z] My function ak() is defined like this: def ak(p): # Convert point on octahedron onto the sphere. # Credit to Anders Kaseorg: https://math.stackexchange.com/questions/5016695/ # input: oct_pt is a Euclidean point on the surface of a unit octagon. # output: UVW on a unit sphere. a = 3.227806237143884260376580641604959964752197265625 # 𝛂 - vis. Kaseorg. p1 = (np.pi * p) / 2.0 tp1 = np.tan(p1) xu, xv, xw = tp1[0], tp1[1], tp1[2] u2, v2, w2 = xu ** 2, xv ** 2, xw ** 2 y0p = xu * (v2 + w2 + a * w2 * v2) ** 0.25 y1p = xv * (u2 + w2 + a * u2 * w2) ** 0.25 y2p = xw * (u2 + v2 + a * u2 * v2) ** 0.25 pv = np.array([y0p, y1p, y2p]) return pv / np.linalg.norm(pv, keepdims=True)  This function is based on a post that I made on the Math StackExchange. Any hints?
One idea I found that helped significantly was to convert this from a minimization problem to a root-finding problem. The root-finders don't support constraints, so you need to change the objective function to force x to be L1 normalized before converting into Cartesian coordinates. def fn_root(op, tx): # octa_point, target_sphere_point norm = np.linalg.norm(op, ord=1) return ak(op / norm) - tx  Once you do that, you can find the point where ak(op / norm) - tx = 0. # Inverse function using numerical optimization def inverse_ak_root(tsp): initial_guess = np.arcsin(tsp) / (np.pi / 2) result = root(fn_root, initial_guess, args=(tsp,), method='hybr', tol=1e-12) assert result.success, result result.x /= np.linalg.norm(result.x, ord=1) return result  (I also changed the initial guess from the center of the octahedron's face to np.arcsin(tsp) / (np.pi / 2). I found this reduced how long it took to converge.) I found that this reduced the error of your current solution by about a factor of 10^7, while using fewer calls to ak() than your previous solution. Full testing code import numpy as np from scipy.optimize import minimize, root def xyz_ll(spt): x, y, z = spt[:, 0], spt[:, 1], spt[:, 2] lat = np.degrees(np.arctan2(z, np.sqrt(x**2. + y**2.))) lon = np.degrees(np.arctan2(y, x)) return np.stack([lat, lon], axis=1) def ll_xyz(ll): # convert to radians. phi, theta = np.radians(ll[:, 0]), np.radians(ll[:, 1]) x = np.cos(phi) * np.cos(theta) y = np.cos(phi) * np.sin(theta) z = np.sin(phi) # z is 'up' return np.stack([x, y, z], axis=1) def ak(p): # Convert point on octahedron onto the sphere. # Credit to Anders Kaseorg: https://math.stackexchange.com/questions/5016695/ # input: oct_pt is a Euclidean point on the surface of a unit octagon. # output: UVW on a unit sphere. a = 3.227806237143884260376580641604959964752197265625 # 𝛂 - vis. Kaseorg. p1 = (np.pi * p) / 2.0 tp1 = np.tan(p1) xu, xv, xw = tp1[0], tp1[1], tp1[2] u2, v2, w2 = xu ** 2, xv ** 2, xw ** 2 y0p = xu * (v2 + w2 + a * w2 * v2) ** 0.25 y1p = xv * (u2 + w2 + a * u2 * w2) ** 0.25 y2p = xw * (u2 + v2 + a * u2 * v2) ** 0.25 pv = np.array([y0p, y1p, y2p]) return pv / np.linalg.norm(pv, keepdims=True) def fn(op, tx): # octa_point, target_sphere_point return np.sum((ak(op) - tx) ** 2.) def fn_root(op, tx): # octa_point, target_sphere_point norm = np.linalg.norm(op, ord=1) return ak(op / norm) - tx def constraint(op): # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron) return np.sum(np.abs(op)) - 1 # Inverse function using numerical optimization def inverse_ak(tsp): initial_guess = np.sign(tsp) * [0.5, 0.5, 0.5] # the centre of the current side # initial_guess = np.arcsin(tsp) / (np.pi / 2) constraints = {'type': 'eq', 'fun': constraint} result = minimize( # maxiter a bit ott maybe, but works. fn, initial_guess, args=(tsp,), constraints=constraints, bounds=[(-1., 1.), (-1., 1.), (-1., 1.)], method='SLSQP', options={'maxiter': 1000, 'ftol': 1e-15, 'disp': False} ) # assert result.success, result return result # Inverse function using numerical optimization def inverse_ak_root(tsp): initial_guess = np.arcsin(tsp) / (np.pi / 2) result = root(fn_root, initial_guess, args=(tsp,), method='hybr', tol=1e-12) assert result.success, result result.x /= np.linalg.norm(result.x, ord=1) return result if __name__ == '__main__': N = 50 np.random.seed(42) lat = np.random.uniform(-90, 90, N) lon = np.random.uniform(-180, 180, N) all_points = np.column_stack([lat, lon]) for i in range(len(all_points)): et = np.array([all_points[i]]) sph = ll_xyz(et) result = inverse_ak(sph[0]) octal = result.x result2 = inverse_ak_root(sph[0]) octal2 = result2.x spherical = ak(octal) rx = xyz_ll(np.array([spherical])) spherical2 = ak(octal2) rx2 = xyz_ll(np.array([spherical2])) print(f'Old Approach') print(f'Original Pt:{et}') print(f'Via inverse:{rx}') print(f'Difference:{rx-et}') print(f'Calls: {result.nfev}') print(f'New Approach') print(f'Via inverse:{rx2}') print(f'Difference:{rx2-et}') print(f'Ratio: {np.linalg.norm(rx-et) / (np.linalg.norm(rx2-et) + 1e-15):g}') print(f'Calls: {result2.nfev}') print()  Note: The parameter 𝛂 has many more digits than necessary. When Python does math on this number, it ignores all digits that do not fit in a double-precision float, which means the number is converted to 3.2278062371438843 internally. Other approaches tried I also tried coming up with an analytic Jacobian for this function using SymPy, but had too much trouble getting automatic differentiation to work. If you were able find this, you'd be able to differentiate the function more quickly and accurately. I also tried a solution based on scipy.interpolate.RbfInterpolator. It got similar accuracy and speed to the minimize method.