使用scipy.optimize.minimize求解旋转角度结果错误问题排查

我有两组x、y坐标，其中一组围绕原点(0,0)旋转后得到另一组，需要确定对应的旋转角度。但运行脚本后得到了约-45°的错误结果，而该示例的正确角度应为5°。

我的代码 #

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# function to calculate the error between both datasets
def error(rotation_angle, data1, data2):

    theta = np.radians(rotation_angle)
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)],[np.sin(theta), np.cos(theta)]])
    
    # rotate data1
    rotated_data1 = np.matmul(data1, rotation_matrix)
    error = np.sum((rotated_data1 - data2) ** 2)

    return error

# calculated sample data for a rotation of 5°
data1 = np.array([[0, 0], [0, 1], [0, 2], [0, 3]])
data2 = np.array([[0, 0], [0.08715574, 0.9961947],[0.17431149, 1.9923894],[0.26146723, 2.98858409]])

# this value is low as expected (6.741376595818406e-17)
print(error(5, data1, data2))

initial_guess = 3

# minimize error to obtain the angle
result = minimize(error, initial_guess, args=(data1, data2), method='L-BFGS-B')

print(result)

fitted_rotation_angle = result.x[0]
print("rotation angle:", fitted_rotation_angle)

运行结果 #

message: Optimization terminated successfully.
       success: True
        status: 0
           fun: 13.101511152940214
             x: [-4.500e+01]
           nit: 27
          nfev: 54
 final_simplex: (array([[-4.500e+01],
                       [-4.500e+01]]), array([ 1.310e+01,  1.310e+01]))
rotation angle: -45.0000000000001

补充信息 #

误差函数的图像看起来是正常的：

angle = np.arange(4, 6., 0.01)
plt.plot(angle,[error(a, data1, data2) for a in angle])
plt.xlabel("rotation angle [°]")
plt.ylabel("error")
plt.show()

请问我哪里出错了？

核心问题：旋转矩阵的乘法顺序错误 #

你的旋转矩阵乘法方向搞反了。当前代码里用rotated_data1 = np.matmul(data1, rotation_matrix)，这是行向量右乘旋转矩阵，对应的是逆时针旋转的逆操作（即顺时针旋转），而你实际需要的是将data1逆时针旋转5°得到data2，或者说用旋转矩阵正确地将data1映射到data2。

另外，误差函数在全局范围内存在多个局部最小值，L-BFGS-B作为局部优化算法，从初始值3°出发，可能会跳到距离更近的局部最小值（-45°），而非全局最小值5°。

解决步骤： #

1. 修正旋转矩阵的乘法顺序：
   把行向量转置为列向量，左乘旋转矩阵后再转置回来，这是标准的列向量左乘旋转矩阵操作，对应逆时针旋转θ角：

   rotated_data1 = np.matmul(rotation_matrix, data1.T).T
   

2. 约束优化范围（可选但推荐）：
   为避免优化算法跳到其他局部最小值，给角度设置合理范围（比如-180到180°），在minimize中添加bounds参数：

   result = minimize(error, initial_guess, args=(data1, data2), method='L-BFGS-B', bounds=[(-180, 180)])
   

修正后的完整代码 #

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

def error(rotation_angle, data1, data2):
    theta = np.radians(rotation_angle)
    rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)],
                               [np.sin(theta), np.cos(theta)]])
    # 修正旋转矩阵乘法顺序
    rotated_data1 = np.matmul(rotation_matrix, data1.T).T
    error = np.sum((rotated_data1 - data2) ** 2)
    return error

data1 = np.array([[0, 0], [0, 1], [0, 2], [0, 3]])
data2 = np.array([[0, 0], [0.08715574, 0.9961947],[0.17431149, 1.9923894],[0.26146723, 2.98858409]])

print(error(5, data1, data2))

initial_guess = 3
# 添加角度范围约束
result = minimize(error, initial_guess, args=(data1, data2), method='L-BFGS-B', bounds=[(-180, 180)])

print(result)
fitted_rotation_angle = result.x[0]
print("rotation angle:", fitted_rotation_angle)

验证结果 #

运行修正后的代码，会得到接近5°的正确结果，误差也会降到极小值。

内容的提问来源于stack exchange，提问作者Florian