Scipy优化无法找到正确的结果

问题描述

我正在尝试使用 scipy.optimize.minimize 来拟合多元函数的参数，但是，无论我向优化器提供多少无噪声数据点，优化器都无法收敛到正确（或接近）回答。

我想知道我使用优化器的方式是否有错误，但我一直在挠头寻找错误。我将不胜感激任何建议或猜测，谢谢！

import numpy as np
from scipy.optimize import minimize
import math

def get_transform(ai,aj,ak,x,y,z):

    i,j,k = 0,1,2
    si,sj,sk = math.sin(ai),math.sin(aj),math.sin(ak)
    ci,cj,ck = math.cos(ai),math.cos(aj),math.cos(ak)
    cc,cs = ci*ck,ci*sk
    sc,ss = si*ck,si*sk
    M = np.identity(4)
    M[i,i] = cj*ck
    M[i,j] = sj*sc-cs
    M[i,k] = sj*cc+ss
    M[j,i] = cj*sk
    M[j,j] = sj*ss+cc
    M[j,k] = sj*cs-sc
    M[k,i] = -sj
    M[k,j] = cj*si
    M[k,k] = cj*ci
    M[0,3] = x
    M[1,3] = y
    M[2,3] = z
    
    return M

def camera_intrinsic(fx,ppx,fy,ppy):
    K = np.zeros((3,3),dtype='float64')
    K[0,0],K[0,2] = fx,ppx
    K[1,1],K[1,2] = fy,ppy

    K[2,2] = 1

    return K

def apply_transform(p,matrix):
    rotation = matrix[0:3,0:3]
  
    T = np.array([matrix[0][3],matrix[1][3],matrix[2][3]])
    transformed = (np.dot(rotation,p.T).T)+T
    return transformed

def project(points_3D,internal_calibration):
    points_3D = points_3D.T
    projections_2d = np.zeros((2,points_3D.shape[1]),dtype='float32')
    camera_projection = (internal_calibration).dot(points_3D)
    projections_2d[0,:] = camera_projection[0,:]/camera_projection[2,:]
    projections_2d[1,:] = camera_projection[1,:]

    return projections_2d.T

    

def error(x):
    global points,pixels
    transform = get_transform(x[0],x[1],x[2],x[3],x[4],x[5])
    points_transfered = apply_transform(points,transform)
    internal_calibration = camera_intrinsic(x[6],x[7],x[8],x[9])
    projected = project(points_transfered,internal_calibration)
    # print(((projected-pixels)**2).mean())
    return ((projected-pixels)**2).mean()


def generate(points,x):

    transform = get_transform(x[0],internal_calibration)
    return projected


points = np.random.rand(100,3)
x_initial = np.random.rand(10)
pixels = generate(points,x_initial)
x_guess = np.random.rand(10)
results = minimize(error,x_guess,method='nelder-mead',tol = 1e-15)
x = results.x
print(x_initial)
print(x)

解决方法

您正在解决最小二乘问题，但尝试使用最小化标量函数的求解器对其进行优化。虽然它可以解决问题，但它的效率非常低。它可能需要更多的迭代，或者根本无法收敛。

更好的方法是使用 least_squares 而不是 minimize。

为了使其正常工作，您应该通过返回一维 numpy 数组而不是标量来修改 error 函数：

def error(x):
    ...
    return (projected-pixels).flatten()

然后调用 least_squares：

results = least_squares(error,x_guess)
x = results.x
print(x_initial)
print(x)
print('error:',np.linalg.norm(error(x)))

此外，error(x) 当前返回 float32 数组，因为 float32 数组是在 project 中创建的。应该用float64代替，否则最小化无法收敛，因为使用32位精度时大部分梯度都为零。

def project(points_3D,internal_calibration):
    ...
    projections_2d = np.zeros((2,points_3D.shape[1]),dtype='float64')

通过这些修改，求解器在大多数情况下会收敛到解，但有时可能无法收敛。这是因为您随机生成问题，因此在某些情况下，问题可能已退化或没有物理意义。此类案件应自行调查。

使用稳健损失（例如'arctan'）代替线性损失也有帮助：

results = least_squares(error,x_guess,loss='arctan')

结果：

[0.68589904 0.68782115 0.83299068 0.02360941 0.19367124 0.54715374
 0.37609235 0.62190714 0.98824796 0.88385802]
[0.68589904 0.68782115 0.83299068 0.02360941 0.19367124 0.54715374
 0.37609235 0.62190714 0.98824796 0.88385802]
error: 1.2269443642313758e-12