Scipy 优化无法找到正确的结果
Scipy optimize unable to find the correct results
我正在尝试使用 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)
您正在求解最小二乘问题,但尝试使用最小化标量函数的求解器对其进行优化。虽然它可以解决问题,但效率很低。它可能需要更多的迭代或根本无法收敛。
更好的方法是使用least_squares
instead of 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)
当前 returns 个 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
我正在尝试使用 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)
您正在求解最小二乘问题,但尝试使用最小化标量函数的求解器对其进行优化。虽然它可以解决问题,但效率很低。它可能需要更多的迭代或根本无法收敛。
更好的方法是使用least_squares
instead of 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)
当前 returns 个 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