1 前言
很多工程都会遇到对数据进行平滑降噪的问题。常规的降噪除了使用滤波方式外，也可以通过优化的方式进行处理。如果场景特殊，也可以自行根据场景构建非线性模型进行优化降噪。这篇文章介绍一种构建方式，以作参考。
2 算法
文献[1]作者使用优化的方式提出了Whittaker平滑器，对于一维输入序列\(\mathbf{x}\)，假设降噪处理后的序列为\(\mathbf{y}\)，则Whittaker算法需要最小化如下目标函数
(1)
$$
e = (\mathbf{x}-\mathbf{y})^TW_m(\mathbf{x}-\mathbf{y})+(D\mathbf{y})^TW_s(D\mathbf{y})
$$
上面公式中\(W_m\)为逼近原始数据的权重对称矩阵, \(W_s\)为度量平滑程度的权重对称矩阵。\(D\mathbf{y}\)表示对\(\mathbf{y}\)求二阶导数。很明显，这个目标函数有两部分，前面部分表示输出数据与原始数据的逼近程度，后面部分表示输出数据的平滑程度。
因为这里是不带其它信息的离散数据求导(默认每个数据点之间的间隔相等)，文献[1]根据前后两点斜率作为一阶结果的方式得到二阶\(D\)为如下值
(2)
$$
D=\begin{bmatrix}
1 & -2 & 1
\\
& 1 & -2 & 1
\\
&& 1 & -2 & 1
\\
&&& 1 & -2 & 1
\\
&&&& 1 & -2 & 1
\\
&&&&& 1 & -2 & 1
\end{bmatrix}
$$
使用一阶判断方法很容易判断(1)里的问题为凸问题，要最小化这个凸问题，令其一阶导数为0，则化简得到如下结果：
(3)
$$
\mathbf{y} = (W_m+D^TW_sD)^{-1}W_m\mathbf{x}
$$
因此公式(3)就是最终得到的平滑算法。
PS: Whittaker算法默认每个点间隔一致，如果输入序列附加额外的属性，比如附带距离属性，则上面的最小化公式不再适用，建议单独针对附加属性来构建相关约束，进而求解。
3 求解
一般地，我们会这样\(W_m=wI,W_s=\lambda I\)来设置权重矩阵，这里设\(G=I+\dfrac{\lambda}{w} H\)，这时公式(3)变成
(4)
$$
\begin{aligned}
\mathbf{y} &=w(wI+\lambda D^TD)^{-1}\mathbf{x}
\\
&=(I+\dfrac{\lambda }{w}D^TD)^{-1} \mathbf{x}
\\
&=(I+\dfrac{\lambda}{w}H)^{-1} \mathbf{x}
\\
&=G^{-1}\mathbf{x}
\end{aligned}
$$
可以发现，参数\(w\)与\(\lambda\)同时缩放，不失一般性，这里设\(\eta = \dfrac{\lambda}{w}\)，每次计算时只需设置参数\(\eta\)即可，因此可以得到 \(G=I+\eta H\)，根据(2)的结果，计算出\(H\)为如下值
(5)
$$
H=\begin{bmatrix}
1 & -2 & 1
\\
-2 & 5 & -4 & 1
\\
1 & -4 & 6 & -4 & 1
\\
& 1 & -4 & 6 & -4 & 1
\\
&&1& -4 & 5 & -2
\\
&&& 1 & -2 & 1
\end{bmatrix}
$$
这里进一步对\(G\)做如下三角分解(\(L\)为下三角矩阵，\(U\)为单位上三角矩阵)
(6)
$$
G=LU
$$
很明显, 公式中\(G\)为对称带状矩阵，总存在如下形式的\(L\)、\(U\)矩阵，使得\(L\)、\(U\)满足公式(6)
(7)
$$
L=\begin{bmatrix}
l_{00}
\\
l_{10} & l_{11}
\\
l_{20} & l_{21} & l_{22}
\\
&l_{31} & l_{32} & l_{33}
\\
&&l_{42} & l_{43} & l_{44}
\\
&&&l_{53} & l_{54} & l_{55}
\end{bmatrix}
$$
$$
U=\begin{bmatrix}
1 & u_{01} & u_{02}
\\
& 1 & u_{12} & u_{13}
\\
&& 1 & u_{23} & u_{24}
\\
&&& 1 & u_{34} & u_{35}
\\
&&&& 1 & u_{45}
\\
&&&&& 1
\end{bmatrix}
$$
因此对于\(n\)阶\(G\)矩阵，根据公式(6)(7)可以得到如下分解公式
(8)
$$\begin{aligned}
&\begin{cases}
\begin{cases}
l_{0,0} = g_{0,0}
\\
u_{0,1} =\dfrac{g_{0,1}}{g_{0,0}}
\\
u_{0,2} =\dfrac{g_{0,2}}{g_{0,0}}
\\
l_{1,0} =g_{1,0}
\\
l_{1,1} =g_{1,1}-g_{1,0}u_{0,1}
\\
u_{1,2}=\dfrac{g_{1,2}-g_{1,0}u_{0,2}}{l_{1,1}}
\\
u_{1,3}=\dfrac{g_{1,3}}{l_{1,1}}
\end{cases}
\\
\\
\begin{cases}
l_{i,i-2}=g_{i,i-2}
\\
l_{i,i-1}=g_{i,i-1}-l_{i,i-2}u_{i-2,i-1}
\\
l_{i,i}=g_{i,i}-l_{i,i-2}u_{i-2,i}-l_{i,i-1}u_{i-1,i}
\\
u_{i,i+1}=\dfrac{g_{i,i+1}-l_{i,i-1}u_{i-1,i+1}}{l_{i,i}}
\\
u_{i,i+2}=\dfrac{g_{i,i+2}}{l_{i,i}}
\end{cases},\ i = 2,3,4,…,n-3
\\
\\
\begin{cases}
l_{n-2,n-4}=g_{n-2,n-4}
\\
l_{n-2,n-3}=g_{n-2,n-3}-l_{n-2,n-4}u_{n-4,n-3}
\\
l_{n-2,n-2}=g_{n-2,n-2}-l_{n-2,n-4} u_{n-4,n-2} – l_{n-2,n-3} u_{n-3,n-2}
\\
u_{n-2,n-1}=\dfrac{g_{n-2,n-1}-l_{n-2,n-3} u_{n-3,n-1}}{l_{n-2,n-2}}
\\
l_{n-1,n-3}=g_{n-1,n-3}
\\
l_{n-1,n-2}=g_{n-1,n-2}-l_{n-1,n-3}u_{n-3,n-2}
\\
l_{n-1,n-1}=g_{n-1,n-1} -l_{n-1,n-2}u_{n-2,n-1}-l_{n-1,n-3}u_{n-3,n-1}
\end{cases}
\end{cases}
\end{aligned}$$
通过公式(8)最终求得三角矩阵\(L\)、\(U\)，最终先求解如下方程得到\(\mathbf{z}\)
(9)
$$
L\mathbf{z}=\mathbf{x}
$$
然后按如下方式求得\(\mathbf{y}\)
(10)
$$
U\mathbf{y}=\mathbf{z}
$$
4 压缩求解
如果直接构建矩阵方式求解公式(8)比较耗时，仔细观察上面数据\(H\)、\(L\)、\(U\)的分布情况，可以发现，这个问题可以采用压缩存储的方式进行计算。结合公式(8)的下标，这里使用新的\(\mathbf{a},\mathbf{b}\)来存储\(L\)中的元素，\(\mathbf{c},\mathbf{d}\)来存储\(U\)中的元素，对LU分解可采用如下方式进行
(11)
$$\begin{cases}
\begin{cases}
b_0 = 1+\eta
\\
c_0 =\dfrac{-2\eta}{b_0}
\\
d_0 =\dfrac{\eta}{b_0}
\\
a_1 =-2\eta
\\
b_1 =1+(5+2c_0)\eta
\\
c_1=\dfrac{2d_0-4}{b_1}\eta
\\
d_1=\dfrac{\eta}{b_1}
\end{cases}
\\
\\
\begin{cases}
a_i=-(4 + c_{i-2})\eta
\\
b_i=1+(6-d_{i-2})\eta -a_i c_{i-1}
\\
c_i=\dfrac{-4\eta-a_i d_{i-1}}{b_i}
\\
d_i=\dfrac{\eta}{b_i}
\end{cases},\ i = 2,3,4,…,n-3
\\
\\
\begin{cases}
a_{n-2}=-(4+c_{n-4})\eta
\\
b_{n-2}=1+(5-d_{n-4})\eta – a_{n-2} c_{n-3}
\\
c_{n-2}=\dfrac{-2\eta-a_{n-2} d_{n-3}}{b_{n-2}}
\\
a_{n-1}=-(2+c_{n-3})\eta
\\
b_{n-1}=1+(1-d_{n-3})\eta -a_{n-1}c_{n-2}
\end{cases}
\end{cases}$$
接着按如下公式计算\(\mathbf{z}\)
(12)
$$\begin{cases}
z_0 = \dfrac{x_0}{b_0}
\\
z_1 = \dfrac{x_1-a_1z_0}{b_1}
\\
\\
z_i=\dfrac{x_i-a_iz_{i-1}-\eta z_{i-2}}{b_i},\ i=2,3,4,…,n-1
\end{cases}$$
最后按如下公式计算\(\mathbf{y}\)
(13)
$$\begin{cases}
y_{n-1} = z_{n-1}
\\
y_{n-2}=z_{n-2}-c_{n-2}y_{n-1}
\\y_i=z_i-c_iy_{i+1}-d_iy_{i+2},\ i=n-3,n-4,n-5,…,0
\end{cases}$$
最终，Whittaker算法平滑降噪的步骤为
step1：设置输入序列\(\mathbf{x}\)以及权重参数比\(\eta\)。
step2：根据公式(11)计算中间变量\(\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d}\)。
step3：根据公式(12)求得中间变量\(z\)。
step4：根据公式(13)求得最终的平滑降噪序列\(\mathbf{y}\)。
6 应用
对于一维问题，直接使用步骤进行平滑；对于二维问题，先对每一列单独进行平滑，然后在所有列平滑完之后，对新的数据每一行再进行平滑即可。
参考:
[1] Zuliana S U , Perperoglou A . Two dimensional smoothing via an optimised Whittaker smoother[J]. Big Data Analytics, 2017, 2(1):6.
[2]https://eigenvector.com/wp-content/uploads/2020/01/WhittakerSmoother.pdf