Ускорить Least_square scipy

Скорее всего, во время scipy.optimize.least_squares() большая часть времени уходит на вычисление остатков, которые в вашем случае принимают форму кода в minReproj().

Однако код, который вы предоставили в minReproj(), похоже, имеет неоптимальное управление памятью, которое можно значительно улучшить с помощью предварительного выделения. Это приведет к значительному увеличению скорости.

Например, vstack() и hstack() страдают от необходимости копировать свои аргументы в память их конечного результата. Учтите, первые несколько строк minReproj(), которые я перегруппировал в gen_affine_OP(). Их можно переписать как gen_affine() с гораздо лучшими таймингами (я также переписал gen_euler_zxz(), чтобы не выделять новую память):

import numpy as np
from math import sin, cos


def gen_euler_zxz(result, psi, theta, sigma):
    c1 = cos(psi)
    s1 = sin(psi)
    c2 = cos(theta)
    s2 = sin(theta)
    c3 = cos(sigma)
    s3 = sin(sigma)
    result[0,0] = (c1 * c3) - (s1 * c2 * s3)
    result[0,1] = (-c1 * s3) - (s1 * c2 * c3)
    result[0,2] = (s1 * s2)
    result[1,0] = (s1 * c3) + (c1 * c2 * s3)
    result[1,1] = (-s1 * s3) + (c1 * c2 * c3)
    result[1,2] = (-c1 * s2)
    result[2,0] = (s2 * s3)
    result[2,1] = (s2 * c3)
    result[2,2] = c2
    return result


def gen_affine(dof):
    result = np.zeros((4, 4))
    gen_euler_zxz(result[:3, :3], dof[0], dof[1], dof[2])
    result[:3, 3] = dof[3:]
    result[3, 3] = 1
    return result


def gen_affine_OP(dof):
    Rmat = gen_euler_zxz(np.empty((3, 3)), dof[0], dof[1], dof[2])
    translationArray = np.array([[dof[3]], [dof[4]], [dof[5]]])
    temp = np.hstack((Rmat, translationArray))
    return np.vstack((temp, [0, 0, 0, 1]))


dof = 1, 2, 3, 4, 5, 6
%timeit gen_affine_OP(dof)
# 100000 loops, best of 3: 16.6 µs per loop
%timeit gen_affine(dof)
# 100000 loops, best of 3: 5.02 µs per loop

Точно так же можно избежать вызова residual = np.vstack((errorA,errorB)), определив больший остаток и предоставив его представление для errorA и errorB.

Другой пример — при создании Ja, Jb, Wa, Wb:

def func_OP(numPoints):
    for i in range(numPoints):
        Ja = np.ones((3))
        Jb = np.ones((3))
        Wa = np.ones((4))
        Wb = np.ones((4))


def func(n):
    Ja = np.empty(3)
    Jb = np.empty(3)
    Wa = np.empty(3)
    Wb = np.empty(3)
    for i in range(n):
        Ja[:] = 1
        Jb[:] = 1
        Wa[:] = 1
        Wb[:] = 1


%timeit func_OP(1000)
# 100 loops, best of 3: 9.31 ms per loop
%timeit func(1000)
# 100 loops, best of 3: 2.2 ms per loop

Кроме того, .flatten() делает копию, которая вам не нужна, а .ravel() просто возвращает представление, но этого достаточно для ваших нужд и происходит намного быстрее:

a = np.ones((300, 300))
%timeit a.ravel()
# 1000000 loops, best of 3: 215 ns per loop
%timeit a.flatten()
# 10000 loops, best of 3: 113 µs per loop

Последние замечания касаются ускорения вашего основного цикла. Я не разрабатываю для этого простую векторизованную перезапись, но вы можете ускорить процесс с помощью Numba (при условии, что он компилируется в необъектном режиме).

Для этого вам также необходимо JIT-декорировать gen_euler_zxz() и заменить вызовы np.matmul() на np.dot() (поскольку np.matmul() не поддерживается Numba.

В конце концов, исправленный minReproj() будет выглядеть так:

from math import sin, cos
import numpy as np
import numba as nb


matmul = np.dot


@nb.jit
def gen_euler_zxz(result, psi, theta, sigma):
    c1 = cos(psi)
    s1 = sin(psi)
    c2 = cos(theta)
    s2 = sin(theta)
    c3 = cos(sigma)
    s3 = sin(sigma)
    result[0, 0] = (c1 * c3) - (s1 * c2 * s3)
    result[0, 1] = (-c1 * s3) - (s1 * c2 * c3)
    result[0, 2] = (s1 * s2)
    result[1, 0] = (s1 * c3) + (c1 * c2 * s3)
    result[1, 1] = (-s1 * s3) + (c1 * c2 * c3)
    result[1, 2] = (-c1 * s2)
    result[2, 0] = (s2 * s3)
    result[2, 1] = (s2 * c3)
    result[2, 2] = c2
    return result


@nb.jit
def minReproj(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix):
    perspectiveProj = np.zeros((4, 4))
    gen_euler_zxz(perspectiveProj[:3, :3], dof[0], dof[1], dof[2])
    perspectiveProj[:3, 3] = dof[3:]
    perspectiveProj[3, 3] = 1

    numPoints = d2dPoints1.shape[0]
    residual = np.empty((2 * numPoints, 3))

    forwardProj = matmul(w2cMatrix, perspectiveProj)
    backwardProj = matmul(w2cMatrix, np.linalg.inv(perspectiveProj))

    Ja = np.empty(3)
    Jb = np.empty(3)
    Wa = np.empty(4)
    Wb = np.empty(4)
    for i in range(numPoints):
        j = i + numPoints

        Ja[2] = Jb[2] = 1.0
        Wa[3] = Wb[3] = 1.0
        Ja[0:2] = d2dPoints1[i, :]
        Jb[0:2] = d2dPoints2[i, :]
        Wa[0:3] = d3dPoints1[i, :]
        Wb[0:3] = d3dPoints2[i, :]

        JaPred = matmul(forwardProj, Wb)
        JaPred /= JaPred[-1]
        residual[i, :] = Ja - JaPred

        JbPred = matmul(backwardProj, Wa)
        JbPred /= JbPred[-1]
        residual[j, :] = Jb - JbPred
    return residual.ravel()

Тестирование на некоторых игрушечных данных показывает значительное ускорение, которое становится драматичным при правильном использовании Numba:

n = 10000
dof = 1, 2, 3, 4, 5, 6
d2dPoints1 = np.random.random((n, 2))
d2dPoints2 = np.random.random((n, 2))
d3dPoints1 = np.random.random((n, 3))
d3dPoints2 = np.random.random((n, 3))
w2cMatrix = np.random.random((3, 4))
x = minReproj_OP(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
y = minReproj(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
z = minReproj_nb(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
print(np.allclose(x, y))
# True
print(np.allclose(x, z))
# True


%timeit minReproj_OP(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
# 1 loop, best of 3: 277 ms per loop
%timeit minReproj(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
# 10 loops, best of 3: 177 ms per loop
%timeit minReproj_nb(dof, d2dPoints1, d2dPoints2, d3dPoints1, d3dPoints2, w2cMatrix)
# 100 loops, best of 3: 5.69 ms per loop