平滑曲线同时保持其下面的区域不变的算法

Question

考虑由点定义的离散曲线 (x1,y1), (x2,y2), (x3,y3), ... ,(xn,yn)

定义一个常量SUM = y1+y2+y3+...+yn.假设我们改变某些k个y点(增加或减少)的值,使得这些变化点的总和小于或等于常数SUM.

在给定以下两个条件的情况下,调整其他y点的最佳方式是什么:

1. y点的总和(y1'+ y2'+ ... + yn')应保持不变,即SUM.
2. 曲线应保留尽可能多的原始形状.

一个简单的解决方案是定义一些delta如下:

delta = (ym1' + ym2' + ym3' + ... + ymk') - (ym1 + ym2 + ym3 + ... + ymk')

并平均分配delta其余的点数.下面ym1'是修改后修改点ym1的值,是修改前修改点的值,将delta作为修改的总差异.

然而,这不能确保完全平滑的曲线,因为变化点附近的区域看起来很粗糙.这个问题是否存在更好的解决方案/算法？

Answer 1

我使用了以下方法，尽管它有点OTT。

考虑将 d[i] 添加到 y[i]，以获得平滑值 s[i]。我们力求最大限度地减少

S = Sum{ 1<=i<N-1 | sqr( s[i+1]-2*s[i]+s[i-1] } + f*Sum{ 0<=i<N | sqr( d[i])}

第一项是曲线的（近似）二阶导数的平方和，第二项惩罚偏离原始曲线。f 是一个（正）常数。一点代数将其改写为

S = sqr( ||A*d - b||)

其中矩阵 A 具有良好的结构，并且实际上 A'*A 是五对角的，这意味着可以有效地求解正规方程组（即 d = Inv(A'*A)*A'*b）。请注意，d 是直接计算的，无需对其进行初始化。

给定该问题的解 d，我们可以计算同一问题的解 d^，但具有约束 One'*d = 0（其中 One 是所有 1 的向量），如下所示

d^ = d - (One'*d/Q) * e
e = Inv(A'*A)*One
Q = One'*e

f 使用什么值？一个简单的方法是在各种 fs 的样本曲线上尝试此过程，并选择一个看起来不错的值。另一种方法是选择平滑度的估计值，例如二阶导数的均方根，然后选择应达到的值，然后搜索给出该值的 f。作为一般规则，f 越大，平滑曲线的平滑度就越低。

这一切的一些动机。目的是找到一条“接近”给定曲线的“平滑”曲线。为此，我们需要平滑度度量（S 中的第一项）和接近度度量（第二项）。为什么采用这些度量？嗯，每个度量都很容易计算，并且每个变量都是二次的（d[]） ；这将意味着该问题成为线性最小二乘的一个实例，对此有有效的算法可用。此外，每个总和中的每一项都取决于变量的附近值，这反过来意味着“逆协方差”（A' *A) 将具有带状结构，因此可以有效地解决最小二乘问题。为什么要引入 f？好吧，如果我们没有 f（或将其设置为 0），我们可以通过设置 d[i] = 来最小化 S -y[i], 得到一条完全平滑的曲线 s[] = 0, 与 y 曲线无关. 另一方面，如果 f 很大，那么为了最小化 s 我们应该专注于第二项，并设置d[i] = 0，我们的“平滑”曲线就是原始曲线。因此，可以合理地假设，当我们改变 f 时，相应的解将在非常平滑但远离 y（小 f）和接近 y 之间变化。 y 但有点粗糙（大 f）。

人们常说，我在这里提倡使用正规方程，但它是解决最小二乘问题的糟糕方法，这通常是正确的。然而，对于“好的”带状系统（就像这里的系统一样），使用正规方程所带来的稳定性损失并不是那么大，而速度的增益却很大。我已经使用这种方法在合理的时间内平滑了具有数千个点的曲线。

要了解 A 是什么，请考虑我们有 4 个点的情况。那么我们的 S 表达式可以归结为：

sqr( s[2] - 2*s[1] + s[0]) + sqr( s[3] - 2*s[2] + s[1]) + f*(d[0]*d[0] + .. + d[3]*d[3]).

例如，如果我们将 s[i] = y[i] + d[i] 代入其中，我们会得到：

s[2] - 2*s[1] + s[0] = d[2]-2*d[1]+d[0] + y[2]-2*y[1]+y[0]

所以我们看到，对于 sqr( ||A*db||) 我们应该采取

A = ( 1 -2  1  0)
    ( 0  1 -2  1)
    ( f  0  0  0)
    ( 0  f  0  0)
    ( 0  0  f  0)
    ( 0  0  0  f)
and
b = ( -(y[2]-2*y[1]+y[0]))
    ( -(y[3]-2*y[2]+y[1]))
    ( 0 )
    ( 0 )
    ( 0 )
    ( 0 )

但在实现中，您可能不想形成 A 和 b，因为它们仅用于形成正规方程项 A'*A 和 A'*b。直接积累这些会更简单。