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

Question

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

pythonoptimizationscipyminimizescipy-optimize

3

我正试图使用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, :]/camera_projection[2, :]

    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],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)
    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)

- Susie

优化器达到了最大函数评估的限制。如果您通过 minimize(...., options = {'maxiter' : 10**9, 'maxfev': 10**9}) 增加此限制，结果会怎样呢？ - joni

需要使用“nelder-mead”方法吗？你可以用其他的方法吗？另外，“tol”参数似乎太低了，因为它仅略高于浮点数的精度限制（约为1e-16）。尝试将其增加到类似于1e-9的值，并检查是否有助于优化收敛。 - fdermishin

'Powell' optimizer performs better and produces results that are somewhat similar to x_initial, but still fails find the correct solution. It seems that there is a problem in computation of the cost function error(x) - fdermishin

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- fdermishin · Accepted Answer

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

更好的方法是使用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数组，因为在project中创建了float32数组。它应该被替换为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