2

所以我有很多（超过 60000 个）数据点 f(x_k) = k，这里k=0,1,2,...,N。

函数是单调递增的，视觉上看起来很平滑。我希望能够找到合适的F(x)，这样对于每个x_k都会发生k <= F(x_k) < k+1。

我应该如何解决这个问题？

数据示例

x       0     1     3     5     8    10    14    16    20    23    27    29    35    37    41
f(x)    0     1     2     3     4     5     6     7     8     9    10    11    12    13    14

阴谋

4

1 回答 1

2

（这看起来有点像查找表。也许是某种图像处理应用程序？我在过去的生活中做过一些需要去舍入的工具。）

这是一次性的问题，还是您会经常这样做，所以您需要速度？

我会把它扔进SLM。由于我没有数据，因此我无法自行测试或给您任何结果，但只要您使用足够数量的结，确保您想要的质量肯定没有问题。您需要在右侧添加额外的结，因为它似乎接近垂直渐近线，因此是奇点。样条曲线通常不喜欢奇点，因为它们本质上仍然是多项式。

更好的是，交换 x 和 y 轴进行拟合，从而拟合 x = f(y)。左端点不是渐近线，所以不再有奇点。现在您需要做的就是将结果限制为单调递增，并向下凹（因此到处都是负二阶导数。）您将需要更少的结来进行逆拟合，但使用足够多的结以使拟合质量足以满足您的需求目标。

要使用逆向拟合，只需反向插值，这是 SLMEVAL 能够做到的。我将看看它对您提供的少量测试数据的影响（仅使用默认的结数）：

x = [0 1 3 5 8 10 14 16 20 23 27 29 35 37 41];
y = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14];
slm = slmengine(y,x,'plot','on','increasing','on');

在此处输入图像描述

所以拟合似乎合理，但我注意到您的数据似乎有点颠簸。可能确实很难获得一个平稳但完全符合您要求的解决方案。

让我们看看它做得如何：

[x;y;slmeval(x,slm,-1)]'
ans =
         0         0    0.0190
    1.0000    1.0000    0.9656
    3.0000    2.0000    2.0522
    5.0000    3.0000    2.9239
    8.0000    4.0000    4.1096
   10.0000    5.0000    4.8419
   14.0000    6.0000    6.1963
   16.0000    7.0000    6.8331
   20.0000    8.0000    8.0638
   23.0000    9.0000    8.9699
   27.0000   10.0000   10.1459
   29.0000   11.0000   10.7088
   35.0000   12.0000   12.2942
   37.0000   13.0000   12.8285
   41.0000   14.0000       NaN

它完全错过了最后一点，拒绝推断。但其余的并不遥远。但是，它们确实不符合您的要求，因为事实并非如此

k <= F(x_k) < k+1

当然，我没有在规范中构建具有这样要求的样条。如果我试图一般地解决这个问题，我可能会编写代码来直接估计曲线上的值，而无需样条中介。然后我可以轻松地强制执行您的约束，找到满足您的误差线要求和单调性的最平滑的点集，这些点也尽可能接近原始数据。当然，这将涉及一个包含 6 万个未知数的大型系统求解。我不知道 lsqlin 将如何处理这么大的问题，但如果时间是一个问题，还有其他解决方案可能会这样做。

同样，以您的测试数据作为小规模示例：

x = [0 1 3 5 8 10 14 16 20 23 27 29 35 37 41]';
n = numel(x);
k = (0:(n-1))';

% The "unrounding" bound constraints
LB = k;
UB = k+1;

% The best fit possible
Afit = speye(n,n);

% And as smooth as possible
ind = 1:(n-2);
% could do this with diff of course
dx1 = x(ind+1) - x(ind);
dx2 = x(ind+2) - x(ind + 1);

% central second finite difference, for unequal spacing
den = dx1.*dx2.*(dx1 + dx2)/2;
Areg = spdiags([dx2./den,-(dx1 + dx2)./den,dx1./den],[0 1 2],n-2,n);
rhs = [k;zeros(n-2,1)];

% monotonicity constraints...
Amono = spdiags(repmat([1 -1],14,1),[0 1],n-1,n);
bmono = zeros(n-1,1);

% choose a value for r, that allows you to control the smoothness
% larger values of r will make the curve smoother, but the bounds
% will always be enforced. I played with it, and r = 5 seemed a
% reasonable compromise here.
r = 5;
yhat = lsqlin([Afit;r*Areg],rhs,Amono,bmono,[],[],LB,UB);

lsqlin 有点不高兴，因为此时它不处理这种形式的稀疏问题。因此它会发出警告，表明它正在将问题转换为完整问题。

Warning: Large-scale algorithm can handle bound constraints only;
    using medium-scale algorithm instead. 
> In lsqlin at 270 
Warning: This problem formulation not yet available for sparse matrices.
Converting to full to solve. 
> In lsqlin at 320 
Optimization terminated.

当然，对于 60k 未知数的问题，这种转换是完全不可接受的。不要在 60k 数据点上尝试！！！！！！！！！！！！！！！！您的计算机将进入深度冻结状态。

它是怎么做的？

disp([x,k,yhat,k+1])
         0         0    0.4356    1.0000
    1.0000    1.0000    1.0000    2.0000
    3.0000    2.0000    2.0504    3.0000
    5.0000    3.0000    3.0000    4.0000
    8.0000    4.0000    4.2026    5.0000
   10.0000    5.0000    5.0000    6.0000
   14.0000    6.0000    6.2739    7.0000
   16.0000    7.0000    7.0000    8.0000
   20.0000    8.0000    8.0916    9.0000
   23.0000    9.0000    9.0000   10.0000
   27.0000   10.0000   10.2497   11.0000
   29.0000   11.0000   11.0000   12.0000
   35.0000   12.0000   12.2994   13.0000
   37.0000   13.0000   13.0000   14.0000
   41.0000   14.0000   14.0594   15.0000

它工作得很好，尽管对于像你这样的大问题来说，这将是一个淫秽的比例。也许还有另一个优化器（可能在 TOMLAB 或其他包中）可以处理大规模稀疏线性问题，受线性和有界约束。您可能还希望将第一个点强制为零，但这很简单。

最后一个选项，如果说 1000 个点是可行的，则使用上述方案一次重新创建 1010 个批次的曲线。lsqlin 应该能够毫无问题地处理这种规模的问题。在末端留下一些重叠，每个重叠区域5个点就足够了。然后平均重叠区域的结果。

于 2013-05-07T15:18:25.357 回答