SciPy minimiza para encontrar a função inversa?

Uma ideia que descobri e que ajudou muito foi converter isso de um problema de minimização em um problema de descoberta de raízes .

Os localizadores de raízes não suportam restrições, então você precisa alterar a função objetivo para forçar x a ser normalizado L1 antes de converter em coordenadas cartesianas.

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

Depois de fazer isso, você pode encontrar o ponto onde 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

(Também alterei o palpite inicial do centro da face do octaedro para np.arcsin(tsp) / (np.pi / 2). Descobri que isso reduziu o tempo que levou para convergir.)

Descobri que isso reduziu o erro da sua solução atual em um fator de cerca de 10^7, enquanto usava menos chamadas para ak() do que sua solução anterior.

Código de teste completo

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

Nota: O parâmetro 𝛂 tem muito mais dígitos do que o necessário. Quando o Python faz cálculos com esse número, ele ignora todos os dígitos que não se encaixam em um float de precisão dupla, o que significa que o número é convertido para 3,2278062371438843 internamente.

Outras abordagens tentadas

Eu também tentei criar um Jacobiano analítico para essa função usando SymPy, mas tive muita dificuldade para fazer a diferenciação automática funcionar. Se você conseguisse encontrar isso, conseguiria diferenciar a função de forma mais rápida e precisa.

Também tentei uma solução baseada em scipy.interpolate.RbfInterpolator. Ela obteve precisão e velocidade semelhantes ao método minimizar.

@Reinderian apontou que eu poderia ter usado code-review, mas quando escrevi a pergunta, não tinha nada para revisar. Estou cauteloso com cross-posting, então vou deixar esta pergunta aqui.

A resposta de @nick-odell foi particularmente útil - e adotei seus pensamentos de todo o coração.

A resposta a seguir é uma atualização da resposta original (que deixei abaixo, especialmente porque foi comentada por Nick)

Isso segue principalmente o código de Nick, mas adicionou o jacobian que ele sugeriu, então pensei que ele (ou uma pessoa interessada nessa área) gostaria de ver onde cheguei. @Nick, estou preocupado que você tenha revelado a localização de muitos submarinos - embora um definitivamente pareça estar preso em um pequeno lago em Nunavut, Canadá.

Os erros de ida e volta do MDC (longitude/latitude) estão bem abaixo de 1 µm, o que é muito mais impressionante do que eu esperava.

A sugestão de Nick para initial_guess parece perder muita viabilidade com N=1e+5 (mesma semente) - e eu substituí por , np.sign(tsp) * 1./3.que teve menos chamadas e parece mais estável com jac.

Uma configuração de tolerância de 1e-12 está retornando um erro máximo de 0,75 µm para N=1e+5, com três cálculos exigindo 32 chamadas para convergência.
Modificar a tolerância para 1e-15 é bom para 4 nm e requer até 41 chamadas para convergência.

Coloquei as rotinas em uma classe (estática) para arrumação e também porque estava com dificuldades para fazer a jacfunção ser chamada corretamente. Além disso, a análise sympy jacobian não é rápida, então é bom gerá-la apenas uma vez.

O Charles Karney's GeographicLibfoi importado para obter métricas de distância confiáveis ​​- este é o sistema geodésico subjacente usado por (quase) todas as bibliotecas GIS. Ele é usado apenas para medir a precisão.

Créditos e elogios ao Nick e também ao Anders.

Estou geralmente interessado em projeções autálicas (ou quase autálicas) entre o octaedro unitário e a esfera unitária - especialmente em uma grade triangular, especialmente aquelas que são fáceis de ver, e onde não há muito desvio do triângulo equilátero (nem preciso dizer que estou ciente de que os cantos do vértice da projeção da esfera estão em 90º, enquanto no octaedro estão em 60º). Estou igualmente interessado em melhorias (velocidade, precisão, eficiência) no código.

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}')

Resposta anterior

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.

E a saída...

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