したがって、かなりの数 (60000 以上) のデータ ポイント f(x_k) = kがあり、ここではk=0,1,2,...,Nです。
関数は単調に増加し、視覚的にはかなり滑らかに見えます。すべてのx_kに対して、k <= F(x_k) < k+1となるようなフィッティングF(x)を見つけることができれば幸いです。
この問題にどのようにアプローチすればよいですか？
データ例
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
(これはルックアップ テーブルに少し似ています。おそらく何らかの画像処理アプリケーションでしょうか?前世で、丸めを解除する必要のあるツールをいくつか使用しました。)
これは 1 回限りの問題ですか、それとも頻繁に行うのでスピードが必要ですか?
私はそれをSLMに投げ込みます。データがないので、自分でテストしたり結果を出すことはできませんが、十分な数のノットを使用する限り、希望する品質の確実な適合に問題はありません。垂直漸近線に近づくように見えるため、右側に追加のノットが必要になるため、特異点になります。スプラインは基本的にまだ多項式であるため、一般に特異点を好まない傾向があります。
さらに良いことに、x 軸と y 軸を交換して近似を行い、x = f(y) を近似します。左の終点は漸近線ではないため、特異点はなくなります。ここで必要なことは、結果が単調に増加し、下にくぼむように制約することです (したがって、どこでも負の 2 次導関数になります)。目標。
逆フィットを使用するには、単純に逆方向に補間します。これは、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
もちろん、私は仕様でそのような要件を持つスプラインを作成しませんでした. この問題を一般的に解決しようとすると、スプラインを介さずに、曲線上の値を直接推定するコードを書くことができます。次に、エラーバーの要件と単調性を満たし、可能な限り元のデータに近い点の最も滑らかなセットを見つけて、制約を簡単に適用できます。もちろん、これには 60k の未知数を伴う大規模なシステム ソルブが必要になります。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 つのポイントがあれば十分です。次に、オーバーラップ領域の結果を平均します。