SciPy 最小化以找到反函数？

我发现一个很有帮助的想法是将其从最小化问题转换为求根问题。

根查找器不支持约束，因此您需要更改目标函数以强制 x 在转换为笛卡尔坐标之前进行 L1 归一化。

def fn_root(op, tx):         # octa_point, target_sphere_point
    norm = np.linalg.norm(op, ord=1)
    return ak(op / norm) - tx

一旦你这样做了，你就可以找到 的点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

（我还将八面体表面中心的初始猜测改为np.arcsin(tsp) / (np.pi / 2)。我发现这减少了收敛所需的时间。）

我发现这将当前解决方案的误差减少了约 10^7 倍，同时对 ak() 的调用次数比之前的解决方案要少。

完整测试代码

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()

注意：参数 𝛂 的位数比实际需要的多得多。当 Python 对这个数字进行数学运算时，它会忽略双精度浮点数中所有不适合的数字，这意味着该数字在内部被转换为 3.2278062371438843。

尝试过其他方法

我还尝试使用 SymPy 为该函数提出解析雅可比矩阵，但在自动微分方面遇到了太多麻烦。如果你能找到这个，你就能更快、更准确地微分该函数。

我也尝试了基于的解决方案scipy.interpolate.RbfInterpolator。它获得了与最小化方法类似的准确性和速度。

@Reinderian 指出我可以使用代码审查，但当我撰写问题时，我没有什么可审查的。我担心交叉发布，所以将这个问题留在这里。

@nick-odell 的回答特别有用- 我全心全意地采纳了他的想法。

以下答案是对原始答案的更新（我将其留在下面，特别是因为 Nick 对此进行了评论）

这主要遵循了 Nick 的代码，但添加了他建议的雅可比矩阵，所以我想他（或对这个领域感兴趣的人）可能想看看我到了哪里。@Nick，我很担心你泄露了许多潜艇的位置 - 虽然其中一艘肯定被困在加拿大努纳武特的一个小湖里。

GCD（经度/纬度）的往返误差远低于 1µm，这比我所希望的要令人印象深刻得多。

Nick 对 initial_guess 的建议似乎在使用 N=1e+5（相同种子）时失去了可行性，而我用它替换了，np.sign(tsp) * 1./3.它调用次数更少，而且看起来更稳定jac。

容差设置为 1e-12 时，N=1e+5 的最大误差为 0.75µm，三次计算需要调用 32 次才能收敛。
将容差修改为 1e-15 可达到 4nm，最多需要调用 41 次才能收敛。

我将例程放入 (静态) 类中，以便整洁，也因为我很难jac正确调用该函数。另外，sympy 雅可比分析速度不快，所以最好只生成一次。

GeographicLib为了获得可靠的距离度量，已导入Charles Karney 的距离度量 - 这是 (几乎) 每个 GIS 库使用的底层测地线系统。它仅用于测量精度。

向 Nick 和 Anders 表示感谢和敬意。

我通常对单位八面体和单位球体之间的等积（或近等积）投影感兴趣 - 尤其是在三角网格上，尤其是那些看起来很舒服并且与等边三角形偏差不大的投影（不用说，我知道球体投影的顶点角为 90º，而在八面体上为 60º）。我同样对代码的改进（速度、准确性、效率）感兴趣。

import sympy as sp
import numpy as np
from geographiclib.geodesic import Geodesic
from scipy.optimize import root


class AK:
    _ALPHA = 3.2278062371438842603765  # 𝛂 - vis. Kaseorg.
    _jac_fn = None

    # def __init__(self):
    #     self.jac_fn = None

    @classmethod
    def xyz_ll(cls, 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)

    @classmethod
    def ll_xyz(cls, 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)

    @classmethod
    def ak(cls, 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 octahedron.
        # output: UVW on a unit sphere.
        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 + cls._ALPHA * w2 * v2) ** 0.25
        y1p = xv * (u2 + w2 + cls._ALPHA * u2 * w2) ** 0.25
        y2p = xw * (u2 + v2 + cls._ALPHA * u2 * v2) ** 0.25
        pv = np.array([y0p, y1p, y2p])
        return pv / np.linalg.norm(pv, keepdims=True)

    @classmethod
    def fn_root(cls, op, tx):         # octa_point, target_sphere_point
        norm = np.linalg.norm(op, ord=1)
        return cls.ak(op / norm) - tx

    @classmethod
    def constraint(cls, op):  # Constraint: |u|+|v|+|w|=1 (surface of the unit octahedron)
        return np.sum(np.abs(op)) - 1

    @classmethod
    def set_jac(cls):
        u, v, w = sp.symbols('u v w')  # Define symbolic variables for inputs
        tan_u = sp.tan(sp.pi * u / 2)
        tan_v = sp.tan(sp.pi * v / 2)
        tan_w = sp.tan(sp.pi * w / 2)

        u2 = tan_u**2
        v2 = tan_v**2
        w2 = tan_w**2

        y0p = tan_u * (v2 + w2 + cls._ALPHA * w2 * v2)**0.25
        y1p = tan_v * (u2 + w2 + cls._ALPHA * u2 * w2)**0.25
        y2p = tan_w * (u2 + v2 + cls._ALPHA * u2 * v2)**0.25

        # Combine outputs into a vector
        y = sp.Matrix([y0p, y1p, y2p])

        # Normalize the vector (divide by its magnitude)
        norm = sp.sqrt(y[0]**2 + y[1]**2 + y[2]**2)
        y_normalized = y / norm

        variables = [u, v, w]
        jacobian = y_normalized.jacobian(variables)
        cls._jac_fn = sp.lambdify(sp.Matrix(variables), jacobian, modules=['numpy'])

    @classmethod
    def inverse_ak_root(cls, tsp):  # Inverse function using numerical optimization
        initial_guess = np.sign(tsp) * 1./3.

        def wrapped_jac(x, _):
            return cls._jac_fn(*x)

        result = root(
            cls.fn_root,
            initial_guess,
            args=(tsp,),
            jac=wrapped_jac,
            method='hybr', tol=1e-13
        )
        assert result.success, result
        result.x /= np.linalg.norm(result.x, ord=1)
        return result


if __name__ == '__main__':
    geod = Geodesic.WGS84
    jk = AK()
    jk.set_jac()
    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 = jk.ll_xyz(et)
        result2 = jk.inverse_ak_root(sph[0])
        octal2 = result2.x
        spherical2 = jk.ak(octal2)
        rx2 = jk.xyz_ll(np.array([spherical2]))
        g = geod.Inverse(rx2[0][0], rx2[0][1], et[0][0], et[0][1])
        print(f'{i:02} inverse:{rx2}; Distance:{1000000.*g['s12']}µm; Calls: {result2.nfev}')

先前的答案

import numpy as np
from scipy.optimize import minimize


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 octahedron.
    # 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 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
    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-20, 'disp': False}
    )
    if result.success:
        return result.x
    else:
        return None


if __name__ == '__main__':
    et = np.array([[48.85837183409, 2.294492063847]])  # in the city of lights.
    sph = ll_xyz(et)
    octal = inverse_ak(sph[0])
    spherical = ak(octal)
    rx = xyz_ll(np.array([spherical]))
    print(f'Original Pt:{et},\nVia inverse:{rx},\n Difference:{rx-et}')  # within a few metres.

并且输出...

Original Pt:[[48.85837183  2.29449206]],
Via inverse:[[48.85837185  2.29449121]],
 Difference:[[ 1.89213694e-08 -8.54396489e-07]]