变坐标系多项式拟合

1 背景
在一些场景里，会对多组二维数据使用一个模型去表达，这时多项式模型成为了常见的选择对象，常规多项式模型的参数求解已经很成熟。然而，有些场景的数据，在当前坐标系下并不能使用多项式模型进行表达，而在另一个坐标系下却能很好的使用多项式进行表达，因此本文的目的是，在已知数据的情况下，找到一个最佳坐标系，以及对应多项式模型系数，进而对模型数据进行更准确表达。

2 理论
2.1 约定
这里约定原始采集\(m\)组二维数据\((x_1,y_1),(x_2,y_2),…,(x_m,y_m)\)。现在需要寻找一个变换矩阵\((R,\mathbf{t})\)使得数据在变换后能尽量贴近多项式模型。这里约定如下变量
(1)
$$
\begin{cases}
\mathbf{p}_i &=\begin{bmatrix}
x_i
\\
y_i
\end{bmatrix}
\\
\\
\mathbf{p}_i’ &=\begin{bmatrix}
x_i’
\\
y_i’
\end{bmatrix}
\\
\\
R&=\begin{bmatrix}
\cos(\theta) &-\sin(\theta)
\\
\sin(\theta)&\cos(\theta)
\end{bmatrix}
\\
\\
\mathbf{t} &=\begin{bmatrix}
t_x
\\
t_y
\end{bmatrix}
\\
\\
\mathbf{p}_i’&=R\mathbf{p}_i+\mathbf{t}
\end{cases}
$$
上面公式里\(\mathbf{p}’\)表示在新坐标系下的点位，这里约定贴近的\(n\)次多项式模型如下
(2)
$$
y’=\sum_{k=0}^na_k(x’)^k
$$
因此，本文的最终目的是找到一个最佳参数\(\theta,t_x,t_y,a_0,a_1,…,a_n\)。其中\(\theta\)反映的新坐标系旋转关系，\(t_x,t_y\)反映的是偏移关系。

由于是数据找模型的问题，因此需要约定什么样的模型可以获得最佳参数。在这里，约定满足如下条件时，获得的参数为最佳参数
(3)
$$
\mathbf{min}:e=\dfrac{1}{2}\sum_{i=1}^m f_i^2
$$
其中\(f_i\)为如下表达式
(4)
$$
\begin{cases}
f_i&=\sum_{k=0}^na_k(x_i’)^k-y_i’
\\
\\
\mathbf{f} &=\begin{bmatrix}
f_1
\\
f_2
\\
…
\\
f_m
\end{bmatrix}
\end{cases}
$$
也就说，这里约定，变换后的数据在新坐标系下，通过构建的多项式，在\(y\)方向上的误差平方和最小即为我们认为的最佳。

2.2 推导
2.2.1 变量降维
从(1)可以发现，\(\mathbf{t}\)变量只影响数据在新坐标系偏移点位，而不影响曲线形态，因此本质上来说\(\mathbf{t}\)设置任何值，都能得到(3)的相同最小化结果，区别只是模型所在新坐标系位置有所不同，而曲线形态完全一致。因此，这里直接将\(\mathbf{t}\)确定出来。考虑到\(t_x\)影响的是所有\(a\)参数,而最优解里\(a_0-t_y\)为固定值，因此这里结合编程时变量范围与精度最终取如下值
(5)
$$
\begin{cases}
t_x &=-\dfrac{1}{n}\sum_{i=1}^ms_i
\\
\\
t_y&=0
\end{cases}
$$
这里\(s_i\)为如下值
(6)
$$
\begin{bmatrix}
s_i
\\
y_i’
\end{bmatrix}=R\begin{bmatrix}
x_i
\\
y_i
\end{bmatrix}
$$

2.2.2 线性方程组求解
通过2.2.1还剩下\(\theta,a_0,a_1,…,a_n\)参数待确定，这里可以发现，当\(\theta\)确定后，(3)最终等价于求如下线下方程组的最小二乘解
(7)
$$
H\mathbf{a} =\mathbf{y}’
$$
其中
(8)
$$
\begin{cases}
\mathbf{a} &=\begin{bmatrix}
a_0
\\
a_1
\\
a_2
\\
…
\\
a_n
\end{bmatrix}
\\
\\
\mathbf{y}’ &=\begin{bmatrix}
y_1′
\\
y_2′
\\
…
\\
y_m’
\end{bmatrix}
\\
\\
\mathbf{H} &=\begin{bmatrix}
1 & x_1′ & x_1’^2 & … & x_1’^n
\\
1 & x_2′ & x_2’^2 & … & x_2’^n
\\
…
\\
1 & x_m’ & x_m’^2 & … & x_m’^n
\end{bmatrix}
\end{cases}
$$

2.2.3 雅可比矩阵
因此整个问题的求解简化为一个变量\(\theta\)的求解。由于(3)是一个非线性问题，因此这里需要得到\(f_i\)对\(\theta\)的导数，进而获得雅可比矩阵。

这里通过(1)(4)(5)(6)可以得到
(9)
$$
\begin{cases}
\dfrac{\partial s_i}{\partial \theta} &=-\sin(\theta)x_i-\cos(\theta)y_i
\\
\\
\dfrac{\partial t_x}{\partial \theta} &=-\dfrac{1}{n}\sum_{i=1}^m\dfrac{\partial s_i}{\partial \theta}=\sin(\theta)u_x+\cos(\theta)u_y
\\
\\
\dfrac{\partial x_i’}{\partial \theta} &=\dfrac{\partial s_i}{\partial \theta}+\dfrac{\partial t_x}{\partial \theta}=\sin(\theta)(u_x-x_i)+\cos(\theta)(u_y-y_i)
\\
\\
\dfrac{\partial y_i’}{\partial \theta} &=\cos(\theta)x_i-\sin(\theta)y_i
\\
\\
\dfrac{\partial f_i}{\partial \theta}&=\dfrac{\partial a_0}{\partial \theta}+\sum_{k=1}^n(x_i’^k\dfrac{\partial a_k}{\partial \theta}+ka_kx_i’^{k-1}\dfrac{\partial x_i}{\partial \theta})-\dfrac{\partial y_i}{\partial \theta}
\end{cases}
$$
上面公式里\(u_x,u_y\)为如下值
(10)
$$
\begin{cases}
u_x&=\dfrac{1}{m}\sum_{i=1}^mx_i
\\
\\
u_y&=\dfrac{1}{m}\sum_{i=1}^my_i
\end{cases}
$$
因此只需要计算得到\(\dfrac{\partial f_i}{\partial \theta}\)就可以得到雅可比矩阵。要计算这个雅克比矩阵，还需要对\(\dfrac{\partial a_i}{\partial \theta}\)参数进行计算。

由于\(\mathbf{a}\)是根据(3)(7)计算得到，因此可以有如下过程
(11)
$$
\begin{cases}
H^TH\mathbf{a} &=H^T\mathbf{y}
\\
&\Downarrow
\\
\dfrac{\partial (H^TH\mathbf{a})}{\partial \theta}&=\dfrac{\partial (H^T\mathbf{y})}{\partial \theta}
\\
&\Downarrow
\\
H^TH\dfrac{\partial \mathbf{a}}{\partial\theta}&=(\dfrac{\partial H}{\partial\theta})^T\mathbf{y}+H^T\dfrac{\partial \mathbf{y}}{\partial\theta}-\dfrac{\partial (H^TH)}{\partial\theta}\mathbf{a}
\\
\\
&=(\dfrac{\partial H}{\partial\theta})^T(\mathbf{y}-H\mathbf{a})+H^T(\dfrac{\partial \mathbf{y}}{\partial\theta}-\dfrac{\partial H}{\partial\theta}\mathbf{a})
\\
\\
&=H^T(\dfrac{\partial \mathbf{y}}{\partial\theta}-\dfrac{\partial H}{\partial\theta}\mathbf{a})-(\dfrac{\partial H}{\partial\theta})^T\mathbf{f}
\\
&\Downarrow
\\
\dfrac{\partial \mathbf{a}}{\partial\theta} &=(H^TH)^{-1}(H^T(\dfrac{\partial \mathbf{y}}{\partial\theta}-\dfrac{\partial H}{\partial\theta}\mathbf{a})-(\dfrac{\partial H}{\partial\theta})^T\mathbf{f})
\end{cases}
$$
上面公式里\(H\)偏导如下
(12)
$$
(\dfrac{\partial H}{\partial\theta})_{ij}=jx_i’^{j-1}\dfrac{\partial x_i’}{\partial\theta}
$$

因此，将(10),(11)带入(9)可以得到最终关于\(\theta\)的雅可比矩阵。

2.3 进阶
① 从2.2可以发现，整个问题就是对一个\(\theta\)的变量求解，由于角度在\(2\pi\)一个循环，因此可以将\(\theta\)求解范围限制在\([0,2\pi]\)内。

② 然而这里本质是拟合曲线形态，\(\theta\)本质反映的是某个坐标系，而镜像坐标系下曲线形态不变，多项式模型x方向平移y方向取反，这时它们对(3)目标一致，因此可以进一步将\(\theta\)范围限制在\([0,\pi]\)内。

③ 由于有雅可比，为了快速求解这种非线性问题，一般会结合凸优化算法进行迭代。这类算法有个特点，当解在迭代值附近时，一般会快速收敛到图凸局部，相反，如果解不在附近，则当前值附近大概率没有最佳解，利用这个特性，且上面两点范围限制，可以在迭代过程中加入求解区域历史记录，进而加速迭代求解。

3 应用
MathSword已经将此功能内置到函数PolyRFit,根据上面理论可知，存在无数个最佳解，这里只返回其中一组最佳参数。